coreBugZJ

此 blog 已弃。

应用LU分解算法求解线性方程组——算法作业 1.2,EOJ 1041

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 
10int 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}

posted on 2011-03-23 15:57 coreBugZJ 阅读(494) 评论(0)  编辑 收藏 引用 所属分类: 课内作业


只有注册用户登录后才能发表评论。
【推荐】超50万行VC++源码: 大型组态工控、电力仿真CAD与GIS源码库
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理