高斯消元

Posted 2020-06-09

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

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 }

其实这个完全可以用二维数组来实现……

)