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
8double a[ L ][ L ], b[ L ], l[ L ][ L ], u[ L ][ L ], x[ L ], y[ L ];
9
10
int main() 
{
11
int n, i, j, k, r;
12 scanf( "%d", &n );
13
for ( i = 1; i <= n; ++i ) {
14 for ( j = 1; j <= n; ++j ) {
15 scanf( "%lf", &a[ i ][ j ] );
16
}
17 }
18 for ( i = 1; i <= n; ++i ) {
19 scanf( "%lf", &b[ i ] );
20 }
21 // Dolittle
22 for ( i = 1; i <= n; ++i ) {
23 for ( j = 1; j <= n; ++j ) {
24 l[ i ][ j ] = u[ i ][ j ] = 0.0;
25 }
26 }
27 for ( k = 1; k <= n; ++k ) {
28 for ( j = k; j <= n; ++j ) {
29 u[ k ][ j ] = a[ k ][ j ];
30 for ( r = 1; r < k; ++r ) {
31 u[ k ][ j ] -= l[ k ][ r ] * u[ r ][ j ];
32 }
33 }
34 if ( iszero( u[ k ][ k ] ) ) {
35 // error
36 }
37 for ( i = k + 1; i <= n; ++i ) {
38 l[ i ][ k ] = a[ i ][ k ];
39 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 for ( i = 1; i <= n; ++i ) {
47 y[ i ] = b[ i ];
48 for ( j = 1; j < i; ++j ) {
49 y[ i ] -= l[ i ][ j ] * y[ j ];
50 }
51 }
52 for ( i = n; i > 0; --i ) {
53 x[ i ] = y[ i ];
54 for ( j = i + 1; j <= n; ++j ) {
55 x[ i ] -= u[ i ][ j ] * x[ j ];
56 }
57 x[ i ] /= u[ i ][ i ];
58 }
59 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
5int 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 for ( i = 1; i <= n; ++i ) {
12 for ( j = 1; j <= n; ++j ) {
13 scanf( "%lf", &a[ i ][ j ] );
14 }
15 }
16 for ( i = 1; i <= n; ++i ) {
17 scanf( "%lf", b + i );
18 }
19
20 // A = LU
21 for ( i = 1; i <= n; ++i ) {
22 for ( j = 1; j <= n; ++j ) {
23 L[ i ][ j ] = U[ i ][ j ] = 0;
24 }
25 }
26 for ( k = 1; k <= n; ++k ) {
27 // a[ k ][ k ] 不能为零
28 for ( j = k; j <= n; ++j ) {
29 U[ k ][ j ] = a[ k ][ j ];
30 }
31 L[ k ][ k ] = 1;
32 for ( i = k + 1; i <= n; ++i ) {
33 L[ i ][ k ] = s = a[ i ][ k ] / a[ k ][ k ];
34 for ( j = k; j <= n; ++j ) {
35 a[ i ][ j ] -= a[ k ][ j ] * s;
36 }
37 }
38 }
39
40 // Ly = b
41 for ( i = 1; i <= n; ++i ) {
42 y[ i ] = b[ i ];
43 for ( j = 1; j < i; ++j ) {
44 y[ i ] -= L[ i ][ j ] * y[ j ];
45 }
46 }
47
48 // Ux = y
49 for ( i = n; i >= 1; --i ) {
50 x[ i ] = y[ i ];
51 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 for ( i = 1; i <= n; ++i ) {
59 printf( "%0.2lf\n", x[ i ] );
60 }
61 return 0;
62}