高斯消元

高斯消元可能这名字听着挺高大上……但其实……没错他就是挺高大上的……他可以用来解线性方程组(不知道线性方程组是什么的自己去百度吧)!!!至于什么逆矩阵之类的我还没有研究。。。这里我先介绍一下主元高斯消元法

x-2y+3z=6(1)

4x-5y+6z=12(2)

7x-8y+10z=21(3)

这是一个三元一次方程组,先来手动解一下这个方程会对后面的方法有很大帮助!!!

我们看到这个方程,一定很希望让他变成如下的形式

x+0y+0z=a

0x+y+0z=b

0x+0y+z=c

所以我就要变形啦!!!先来(2)-(1)*4,然后再来一发(3)-(1)*7,这样除了(1)中的x都被消掉了。同理,我们可以得到唯一含有y的式子和唯一含有z的式子

x-2y+3z=6

3y-6z=-12(4)

6y-11z=-21(5)

现在我们的target是只让(4)中含有未知数y,消去其他方程中的y,首先把(4)/3

然后就得到了y-2z=-4(6)

最后我们用(1)-(6)*(-2)再来(5)-(6)*6

x-z=-2(7)

y-2z=-4

z=3(8)

然后随便算一算就知道(x=1,y=2,z=3)的方程的解了。这样即使未知数增加,也可以用同样的方法来求

下面我给出代码(其实我在学的时候觉得这个代码是不太好理解的…………

 1 typedef vector<double> vec;
 2 typedef vector<vec> mat;
 3 
 4 vec gauss_jordan(const mat &A,const vec &b){
 5     int n=A.size();
 6     mat B(n,vec(n+1));
 7     for(int i=0;i<n;i++){
 8         for(int j=0;j<n;j++)B[i][j]=A[i][j];
 9     }
10     for(int i=0;i<n;i++)B[i][n]=b[i];
11     //把b和A合并一发,这样比较好处理 
12     for(int i=0;i<n;i++){
13         int pivot=i;
14         for(int j=i;j<n;j++){
15             if(abs(B[j][i])>abs(B[pivot][i]))pivot=j;
16         }
17         swap(B[i],B[pivot]);
18         //选择要消去的未知数系数的绝对值尽量大的方程可以减小误差 
19         for(int j=i+1;j<=n;j++)B[i][j]/=B[i][i];
20         //把正在处理的式子未知数系数变为1 
21         for(int j=0;j<n;j++){
22             if(i!=j){
23                 //第j个式子中消去第i个未知数 
24                 for(int k=i+1;k<=n;k++)B[j][k]-=B[j][i]*B[i][k];
25             }
26         }
27     }
28     vec x(n);
29     for(int i=0;i<n;i++){x[i]=B[i][n];}
30     //最后存放在数组最右边的就是答案 
31     return x;
32 }

其实这个完全可以用二维数组来实现……这个代码其实主要的一件事就是从第一个未知数开始一个一个进行处理,变成我们一开始所希望的形式

x+0y+0z=a

0x+y+0z=b

0x+0y+z=c

posted @ 2016-01-26 11:25  543~  阅读(224)  评论(0编辑  收藏  举报