高斯消元法——求解线性方程组

高斯消元法——求解线性方程组

学习了《挑程》中的高斯消元法，它是求解最基础的线性方程组（未知数个数和方程个数相等，并且有唯一解）的算法。

首先举一个例子：求解如下方程组：

⎧⎩⎨x−2y+3z&#①4x−5y+6z&#②7x−8y+10z&#③ $$begin{cases} x-2y+3z = 6..........①\ 4x-5y+②\ 7x-8y+10z = 21.....③ end{cases}$$

我们手算一下这个方程组，过程如下：
1. ② - 4*①；③ - 7*①，得到如下式子：

⎧⎩⎨x−2y+3z&#①0x+3y−6z=−12.......②0x−6y−11z=−21.....③ $$begin{cases} x-2y+3z = 6.............①\ 0x+3y-6z=-12.......②\ 0x-6y-11z = -21.....③ end{cases}$$
2. ②两边同除以3，得到：
⎧⎩⎨x−2y+3z&#①0x+y−2z=−4..........②0x−6y−11z=−21.....③ $$begin{cases} x-2y+3z = 6.............①\ 0x+y-2z=-4..........②\ 0x-6y-11z = -21.....③ end{cases}$$
3. ① - ②*(-2)；③-②*6，得到：
⎧⎩⎨x+0y−z=−2.............①0x+y−2z=−4...........②0x+0y+z&#③ $$bbox[5px,border:2px solid red]{begin{cases} x+0y-z = -2.............①\ 0x+y-2z=-4...........②\ 0x+0y+z =3..............③ end{cases}}$$
4. 从③式得到 z=3 $z=3$，再代入②式得到 y=2 $y = 2$，再代入到①式得到 x=1 $x=1$

我们这是平时手算的过程，但是一旦学习过线性代数这门课程，我们就不用再写方程式了，直接写矩阵的形式，还是这个例子，我们便可以写成如下矩阵：

⎡⎣⎢⎢147−2−5−8361061221⎤⎦⎥⎥ $$bbox[yellow,5px,border:2px solid red] { left[begin{array}{ccc|c}1&-2&3&6\4&-5&6&12\7&-8&10&21end{array} right]}$$

此矩阵为增广矩阵，竖线之前的矩阵为方程组的系数矩阵 A $A$，而竖线之后的列矩阵为常数矩阵 b$b$。
所以方程组的运算就是对增广矩阵进行相关操作，而在我们的高级语言中，可以通过二维数组来表示此矩阵，当然也可以使用vector来表示。

具体代码如下：

#include<cstdio>
#include<vector>
#include<cmath>
#include<algorithm>
using namespace std;
const double EPS = 1E-8;
typedef vector<double> vec;
typedef vector<vec> mat;vec gauss_jordan(const mat &A, const vec &b)
{int n = A.size();           //n个未知数 mat B(n, vec(n+1));         //初始化增广矩阵 for(int i=0; i<n; i++){ for(int j=0; j<n; j++){ B[i][j] = A[i][j];  //系数矩阵 }}for(int i=0; i<n; i++){B[i][n] = b[i];         //常数矩阵 }   //把正在处理的未知数系数的绝对值最大的式子换到第i行 for(int i=0; i<n; i++){int pivot = i;for(int j=i; j<n; j++){if(abs(B[j][i]) > abs(B[pivot][i])) pivot = j;  }   swap(B[i], B[pivot]);//无解或者有无穷解 if(abs(B[i][i]) < EPS) return vec();//把正在处理的未知数的系数变为1double k = B[i][i];for(int j=i; j<=n; j++) B[i][j] /= k;for(int j=0; j<n; j++){if(i != j){//从第j个式子中消去第i个未知数 for(int k=i+1; k<=n; k++){B[j][k] -= B[j][i] * B[i][k];} }   }
//      for(int j=0; j<n; j++)
//      {
//          for(int k=0; k<n+1; k++)
//          {
//              printf("%f ",B[j][k]);
//          }   
//          printf("n");
//      } 
//      printf("--------------------n");}vec x(n);for(int i=0; i<n; i++){x[i] = B[i][n];}return x;
}int main()
{freopen(&#", "r", stdin);int n;scanf("%d", &n);mat A(n, vec(n));vec b(n);double u;for(int i=0; i<n; i++){for(int j=0; j<n; j++){scanf("%lf", &A[i][j]);}}for(int i=0; i<n; i++){scanf("%lf", &b[i]);}vec c = gauss_jordan(A, b);for(int i=0; i<c.size(); i++){printf("%f ",c[i]);}return 0;
}