Description
应用LU 分解算法求解n×n的线性方程组Ax=b。
数据保证可以LU分解并且有唯一解。
Input
第1行为一个整数n(0 接下去的n行表示了系数矩阵A,每行有n个整数。
再接下去的n行表示了b,每行只有一个整数。
Output
输出有n行,每行有1个小数(精确到0.01),表示方程组的解。
Sample Input
2
1 3
7 8
4
15
Sample Output
1.00
1.00
我的代码:
Dolittle 算法:
1
#include <stdio.h>
2
3
#define eps 0.0001
4
#define iszero(x) ( (-eps<(x)) && ((x)<eps) )
5
6
#define L 30
7
8
double a[ L ][ L ], b[ L ], l[ L ][ L ], u[ L ][ L ], x[ L ], y[ L ];
9
10![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
int main()
{
11
int n, i, j, k, r;
12
scanf( "%d", &n );
13![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( i = 1; i <= n; ++i )
{
14![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( j = 1; j <= n; ++j )
{
15
scanf( "%lf", &a[ i ][ j ] );
16
}
17
}
18![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( i = 1; i <= n; ++i )
{
19
scanf( "%lf", &b[ i ] );
20
}
21
// Dolittle
22![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( i = 1; i <= n; ++i )
{
23![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( j = 1; j <= n; ++j )
{
24
l[ i ][ j ] = u[ i ][ j ] = 0.0;
25
}
26
}
27![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( k = 1; k <= n; ++k )
{
28![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( j = k; j <= n; ++j )
{
29
u[ k ][ j ] = a[ k ][ j ];
30![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( r = 1; r < k; ++r )
{
31
u[ k ][ j ] -= l[ k ][ r ] * u[ r ][ j ];
32
}
33
}
34![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
if ( iszero( u[ k ][ k ] ) )
{
35
// error
36
}
37![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( i = k + 1; i <= n; ++i )
{
38
l[ i ][ k ] = a[ i ][ k ];
39![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( r = 1; r < k; ++r )
{
40
l[ i ][ k ] -= l[ i ][ r ] * u[ r ][ k ];
41
}
42
l[ i ][ k ] /= u[ k ][ k ];
43
}
44
l[ k ][ k ] = 1.0;
45
}
46![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( i = 1; i <= n; ++i )
{
47
y[ i ] = b[ i ];
48![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( j = 1; j < i; ++j )
{
49
y[ i ] -= l[ i ][ j ] * y[ j ];
50
}
51
}
52![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( i = n; i > 0; --i )
{
53
x[ i ] = y[ i ];
54![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( j = i + 1; j <= n; ++j )
{
55
x[ i ] -= u[ i ][ j ] * x[ j ];
56
}
57
x[ i ] /= u[ i ][ i ];
58
}
59![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( i = 1; i <= n; ++i )
{
60
printf( "%0.2lf\n", x[ i ] );
61
}
62
return 0;
63
}
高斯
1
#include <stdio.h>
2
3
#define N 30
4
5![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
int main()
{
6
double a[ N ][ N ], L[ N ][ N ], U[ N ][ N ], b[ N ], x[ N ], y[ N ];
7
int n, i, j, k;
8
double s;
9
10
scanf( "%d", &n );
11![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( i = 1; i <= n; ++i )
{
12![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( j = 1; j <= n; ++j )
{
13
scanf( "%lf", &a[ i ][ j ] );
14
}
15
}
16![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( i = 1; i <= n; ++i )
{
17
scanf( "%lf", b + i );
18
}
19
20
// A = LU
21![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( i = 1; i <= n; ++i )
{
22![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( j = 1; j <= n; ++j )
{
23
L[ i ][ j ] = U[ i ][ j ] = 0;
24
}
25
}
26![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( k = 1; k <= n; ++k )
{
27
// a[ k ][ k ] 不能为零
28![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( j = k; j <= n; ++j )
{
29
U[ k ][ j ] = a[ k ][ j ];
30
}
31
L[ k ][ k ] = 1;
32![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( i = k + 1; i <= n; ++i )
{
33
L[ i ][ k ] = s = a[ i ][ k ] / a[ k ][ k ];
34![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( j = k; j <= n; ++j )
{
35
a[ i ][ j ] -= a[ k ][ j ] * s;
36
}
37
}
38
}
39
40
// Ly = b
41![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( i = 1; i <= n; ++i )
{
42
y[ i ] = b[ i ];
43![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( j = 1; j < i; ++j )
{
44
y[ i ] -= L[ i ][ j ] * y[ j ];
45
}
46
}
47
48
// Ux = y
49![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( i = n; i >= 1; --i )
{
50
x[ i ] = y[ i ];
51![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( j = n; j > i; --j )
{
52
x[ i ] -= U[ i ][ j ] * x[ j ];
53
}
54
x[ i ] /= U[ i ][ i ];
55
}
56
57
// output
58![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for ( i = 1; i <= n; ++i )
{
59
printf( "%0.2lf\n", x[ i ] );
60
}
61
return 0;
62
}