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

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

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

$$\begin{cases} x-2y+3z = 6..........①\\ 4x-5y+6z=12.......②\\ 7x-8y+10z = 21.....③ \end{cases}$$

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

$$\begin{cases} x-2y+3z = 6.............①\\ 0x+3y-6z=-12.......②\\ 0x-6y-11z = -21.....③ \end{cases}$$
2. ②两边同除以3，得到：
$$\begin{cases} x-2y+3z = 6.............①\\ 0x+y-2z=-4..........②\\ 0x-6y-11z = -21.....③ \end{cases}$$
3. ① - ②*(-2)；③-②*6，得到：
$$\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$，再代入②式得到 $y = 2$，再代入到①式得到 $x=1$

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

$$\bbox[yellow,5px,border:2px solid red] { \left[ \begin{array}{ccc|c} 1&-2&3&6\\ 4&-5&6&12\\ 7&-8&10&21 \end{array} \right]}$$

此矩阵为增广矩阵，竖线之前的矩阵为方程组的系数矩阵 $A$，而竖线之后的列矩阵为常数矩阵 $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();

        //把正在处理的未知数的系数变为1
        double 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("input.txt", "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;
}