高斯列主消元详解及模板

高斯列主消元详解及模板

分类： 高斯消元2013-08-29 19:28 213人阅读 评论(0) 收藏 举报 高斯消元

采用高斯先列主元消元法求解线性方程组AX=b

方法说明（以4阶为例）：

（1）第1步消元——在增广矩阵（A，b）第一列中找到绝对值最大的元素，将其所在行与第一行交换，再对（A，b）

做初等行变换使原方程组转化为如下形式：注：“*”代表非0。

[cpp] view plaincopy

1. *x1+*x2+*x3+*x4=y1;  
2.  0 +*x2+*x3+*x4=y2;  
3.  0 +*x2+*x3+*x4=y3;  
4.  0 +*x2+*x3+*x4=y4;  

（2）第2步消元——在增广矩阵（A，b）中的第二列中（从第二行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A，b）做初等行变换使原方程组转化为：

[cpp] view plaincopy

1. *x1+*x2+*x3+*x4=y1;  
2.  0 +*x2+*x3+*x4=y2;  
3.  0 + 0 +*x3+*x4=y3;  
4.  0 + 0 +*x3+*x4=y4;  

（3）第3步消元——在增广矩阵（A，b）中的第三列中（从第三行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A，b）做初等行变换使原方程组转化为：

[cpp] view plaincopy

1. *x1+*x2+*x3+*x4=y1;  
2.  0 +*x2+*x3+*x4=y2;  
3.  0 + 0 +*x3+*x4=y3;  
4.  0 + 0 + 0 +*x4=y4;  

（4）按x4 ; x3; ; x1 的顺序回代求解出方程组的解

再回代求解即可。

下面是高斯消元模板：

[cpp] view plaincopy

1. #include<iostream>  
2. #include<cstring>  
3. #include<string>  
4. #include<cstdio>  
5. #include<cmath>  
6. #define eps 0.0000001  
7. using namespace std;  
8. double x[102];  
9. double g[102][102];  
10. int flag;  
11.   
12. void gauss(int n,int m)  
13. {  
14.     int row,col,i,j,k;  
15.     for(row=1,col=1;row<n,col<m;row++,col++)  
16.     {  
17.         k=row;  
18.         for(i=row+1;i<=n;i++)  //列主元  
19.             if(fabs(g[i][col])>fabs(g[k][col]))  
20.                 k=i;  
21.         if(k!=row)   //行交换  
22.         {  
23.             for(i=col; i<=m; i++)  
24.                 swap(g[k][i],g[row][i]);  
25.         }  
26.         if(fabs(g[row][col])<eps) //主元是0  
27.         {  
28.             flag=1;  
29.             return;  
30.         }  
31.         for(i=row+1; i<=n; i++)  //主元不是0把下面的行第一个值全部变为0  
32.         {  
33.             if(fabs(g[i][col])<eps)  
34.                 continue;  
35.             double t=g[i][col]/g[row][col];  
36.             g[i][col]=0.0;  
37.             for(j=col+1;j<=m;j++)  
38.                 g[i][j]-=t*g[row][j];  
39.         }  
40.     }  
41.     if(fabs(g[n][n])<eps)  
42.     {  
43.         flag=1;  
44.         return;  
45.     }  
46.     for(i=n;i>=1;i--)  //回代求解  
47.     {  
48.         x[i]=g[i][m];  
49.         for(j=i+1;j<=n;j++)  
50.             x[i]-=x[j]*g[i][j];  
51.         x[i]/=g[i][i];  
52.     }  
53. }  
54.   
55. int main()  
56. {  
57.     int n,m,i,j;  
58.     while(~scanf("%d%d",&n,&m))  
59.     {  
60.         flag=0;  
61.         for(i=1;i<=n;i++)  
62.             for(j=1;j<=m;j++)  
63.               scanf("%lf",&g[i][j]);  
64.                 
65.         gauss(n,m);  
66.         if(flag)  
67.         {  
68.             cout<<"方程组无解或有无数组解"<<endl;  
69.             continue;  
70.         }  
71.         cout<<"方程组的解为:"<<endl;  
72.         for(i=1;i<=n;i++)  
73.             cout<<"x"<<i<<":   "<<x[i]<<endl;  
74.     }  
75.     return 0;  
76. }  
77.   
78. /* 
79. 2 3 
80. 2 3 5 
81. 1 4 6 
82. 2 3 
83. 2 3 5 
84. 0 0 6 
85. 2 3 
86. 2 3 5 
87. 0 4 6 
88. */  

再推荐一首歌：

这一路走来来自杨宗纬