Jacobi 迭代法 解方程

Code:

1. #include<stdio.h>   
2. #include<math.h>   
3.   
4. #define N 6   
5. #define M 6   
6.   
7. double max(double a[M]);   
8. double jocabi(double a[][M],double b[N],double x[M],double p);   
9.   
10. /* a 为 系数 矩阵 N 行 M列  
11.  b  为 Ax=b中的 b；M个元素  
12.  x  为初值 维数为 M的 一维数组  
13.  p  为精度  
14. */  
15.   
16. int main()   
17. {   
18.         double x[M],p;   
19.         double a[N][M]={{4,-1,0,-1,0,0},   
20.                         {-1,4,-1,0,-1,0},   
21.                         {0,-1,4,0,0,-1},   
22.                         {-1,0,0,4,-1,0},   
23.                         {0,-1,0,-1,4,-1},   
24.                         {0,0,-1,0,-1,4}};   
25.         double b[N]={0,5,0,6,-2,6};   
26.         int i;   
27.         printf("请输入初值：/n");   
28.         for(i=0;i<M;i++)   
29.         {   
30.           printf("x%d(0)=",i+1);   //初值   
31.   
32.           scanf("%lf",&x[i]);   
33.         }   
34.         printf("请输入求解精度：/n");   
35.         printf("p=");   
36.         scanf("%lf",&p);   
37.         printf("/n");   
38.   
39.         for(i=0;i<M;i++)   
40.         printf("x%d(0)=%lf/t/t",i+1,i,x[i]);  //输出初始值   
41.         printf("/n/n");   
42.         jocabi(a,b,x,p);   
43. }   
44. double max(double a[M])   
45. {   
46.     int i;   
47.     double t;   
48.     for(i=0,t=a[0];i<M;i++)   
49.     if(a[i]>t) t=a[i];// t 始终存储 最大值   
50.     return t;   
51. }   
52. double jocabi(double a[][M],double b[N],double x[M],double p)   
53. {   
54.         int k,i,j;   
55.         double e[M]={0},y[M]={0},temp;  //M为列   
56.         printf("迭代列表如下：/n");   
57.         for(k=0;k<5000; )   
58.         {   
59.               for(i=0;i<M;i++)   
60.                {   
61.                     temp=0;   
62.                     for(j=0;j<M;j++)   
63.                    {   
64.                       if(i==j) continue;// 主对角线的元素 已被 提取来了   
65.                      else  
66.                      {   
67.                        temp+=a[i][j]*x[j];//系数乘以 变量 x[j];   
68.                      }   
69.                    }   
70.                      y[i]=(b[i]-temp)/a[i][i];   
71.                      printf("x%d(%d)=%lf/t/t",i+1,k+1,y[i]);   
72.                      e[i]=fabs(x[i]-y[i]);//e[i] 存储 每次的 x[i]-y[j]差值   
73.                }   
74.             printf("/n");   
75.             if(max(e)<p)//e的最大值 小于精度 p 时 输出结果   
76.             {   
77.                     printf("/n");   
78.                     printf("该方程组的解为：/n");   
79.                     for(j=0;j<N;j++)   
80.                     {   
81.                       printf("x%d=%lf /t /t",j+1,y[j]);//未知数的从 X1 开始，而 数组y[j]从 0开始   
82.                     }   
83.                      printf("/n");   
84.                     break;   
85.             }   
86.                k++;   
87.             for(j=0;j<M;j++)   
88.             x[j]=y[j];   
89.         }   
90.         if(k==5000)//迭代数组次数超过 5000 超过 迭代次数 结束   
91.         printf("此方程迭代失败/n");   
92.         return 0;   
93. }