「Note」高斯消元（未完）

引入

求解一个线性方程组

\[\begin{cases} a_{1,1}x_1+a_{1,2}x_2+...+a_{1,n}x_n=b_1 \\ a_{2,1}x_1+a_{2,2}x_2+...+a_{2,n}x_n=b_2 \\ \vdots \\ a_{n,1}x_1+a_{n,2}x_2+...+a_{n,n}x_n=b_n \\ \end{cases} \]

可以将方程组进行改写：

\[ \begin{pmatrix} a_{1,1}x_1+a_{1,2}x_2+...+a_{1,n}x_n \\ a_{2,1}x_1+a_{2,2}x_2+...+a_{2,n}x_n \\ \vdots \\ a_{n,1}x_1+a_{n,2}x_2+...+a_{n,n}x_n \\ \end{pmatrix} = \begin{pmatrix} b_{1}\\ b_{2} \\ \vdots \\ b_{n} \end{pmatrix} \]

\[ \begin{pmatrix} a_{1,1} & a_{1,2} & \dots & a_{1,n}\\ a_{2,1} & a_{2,2} & \dots & a_{2,n} \\ \vdots & \vdots& \vdots & \vdots \\ a_{n,1} & a_{n,2} & \dots & a_{n,n} \end{pmatrix} \begin{pmatrix} x_{1}\\ x_{2} \\ \vdots \\ x_{n} \end{pmatrix} = \begin{pmatrix} b_{1}\\ b_{2} \\ \vdots \\ b_{n} \end{pmatrix} \]

\[ 记 A= \begin{pmatrix} a_{1,1} & a_{1,2} & \dots & a_{1,n}\\ a_{2,1} & a_{2,2} & \dots & a_{2,n} \\ \vdots & \vdots& \vdots & \vdots \\ a_{n,1} & a_{n,2} & \dots & a_{n,n} \end{pmatrix} X= \begin{pmatrix} x_{1}\\ x_{2} \\ \vdots \\ x_{n} \end{pmatrix} B= \begin{pmatrix} b_{1}\\ b_{2} \\ \vdots \\ b_{n} \end{pmatrix} \]

\(A\) 表示系数矩阵， \(X\) 表示未知量矩阵， \(B\) 表示常数项矩阵。

所以原方程为 \(AX=B\) 。

相关概念

我们把系数矩阵和常数项矩阵放在一起后的矩阵，称为增广矩阵。

\[(A|B)= \left\{ \begin{array}{cccc|c} a_{1,1} & a_{1,2} & \dots & a_{1,n} & b_1\\ a_{2,1} & a_{2,2} & \dots & a_{2,n} & b_2\\ \vdots & \vdots& \vdots & \vdots & \vdots\\ a_{n,1} & a_{n,2} & \dots & a_{n,n} & b_n \end{array} \right\} \]

思路

对于解方程，最基本的方法就是消元。而高斯消元顾名思义也是要用到这个方法。

首先，要了解到解方程的几个初等变换：

1.互换两个方程的位置
2.用一个非零数乘某个方程的两边
3.将一个方程的两边同乘以某常数加到另一个方程

而我们在解线性方程组的过程中对 \(AX=B\) 进行初等变换实际就是在对 ta 的增广矩阵进行行初等变换。

不妨举个例子：

\[\begin{cases} x_1+2x_2+x_3=12 \\ 2x_1+x_2+...+x_3=11 \\ 2x_1+3x_2+...+3x_3=25 \\ \end{cases} \]

则增广矩阵为：

\[\left\{ \begin{array}{ccc|c} 1 & 2 & 1 & 12 \\ 2 & 1 & 1 & 11 \\ 2 & 3 & 3 & 25 \end{array} \right\} \]

我们可以对其进行以上的变换：

\[将第二，三行减去第一行 \times 2。\\ \left\{ \begin{array}{ccc|c} 2 & 4 & 2 & 24 \\ & -3 & -1 & -13 \\ & -1 & 1 & 1 \end{array} \right\} \\ 将第三行减去第二行 \times \frac{1}{3} 。\\ \left\{ \begin{array}{ccc|c} 2 & 4 & 2 & 24 \\ & -3 & -1 & -13 \\ & & \frac{4}{3} & \frac{16}{3} \end{array} \right\} \\ 此时可以通过第三行解得 x_3=\frac{\frac{16}{3}}{\frac{4}{3}}=4\\ 代入第二行可解得 x_2=\frac{-9}{-3}=3\\ 再将 x_2,x_3 代入第一行，即可解得 x_1=\frac{4}{2}=2 \]

总的来说，我们在第 \(k\) 消元过程中，找定一行为这次消元的主行，则 \(a_{k,k}\) 就为这次消元的主元素。

我们消元的目的就是将第 \(k+1\)~\(n\) 行的 \(x_k\) 的系数消为 \(0\) 。

所以我们将第 \(i\) 行（ \(i \in [k+1,n]\) ），减去第 \(k\) 行乘上 \(\frac{a_{k,i}}{a_{k,k}}\) ，使得两行 \(x_k\) 的系数相同，进而消掉 \(a_{k,i}\) 。

最终到达的效果就是第 \(n\) 行，第 \(1\)~\(n-1\) 列的系数都为 \(0\) 。因此可以直接解得 \(x_n\)。

而 \(n-1\) 行，只有第 \(n-1\) ，\(n\) 列上有系数。由于我们解得 \(x_n\) ，于是可以代入 \(x_n\) 解得 \(x_{n-1}\) 。

以此类推下去，就可以解出整个方程组的解。

解的状态判定

对于方程组，通常有 无解，无数解，唯一解 三种状态。

代码

\\P3389
#include <cmath>
#include <cstdio>
#include <algorithm>

#define Maxn 105
const double Eps = 1e-8;

using namespace std;

int n;
double a[Maxn][Maxn];

int Gauss () {
    int r, c;
    for (r = 1, c = 1; c <= n && r <= n; c ++) {
        int pos = r;
        for (int i = r; i <= n; i ++) {
            if (fabs (a[i][c]) > fabs (a[pos][c])) pos = i;
        }
        if (fabs (a[pos][c]) < Eps) continue;
        for (int i = c; i <= n + 1; i ++) swap (a[pos][i], a[r][i]);
        for (int i = n + 1; i >= c; i --) a[r][i] /= a[r][c];
        for (int i = r + 1; i <= n; i ++)
            for (int j = n + 1; j >= c; j --) a[i][j] -= a[r][j] * a[i][c];
        r ++;
    }
    if (r <= n) {
        for (int i = r + 1; i <= n; i ++) {
            if (fabs (a[i][n + 1]) > Eps) return 2;
        }
        return 1;
    } else {
        for (int i = n; i >= 1; i --)
            for (int j = i + 1; j <= n; j ++) a[i][n + 1] -= a[i][j] * a[j][n + 1];
        return 0;
    }
}

int main () {
    scanf ("%d", &n);

    for (int i = 1; i <= n; i ++) {
        for (int j = 1; j <= n + 1; j ++) {
            scanf ("%lf", &a[i][j]);
        }
    }

    if (Gauss ()) {
        printf ("No Solution");
    } else {
        for (int i = 1; i <= n; i ++) {
            printf ("%.2lf\n", a[i][n + 1]);
        }
    }
    return 0;
}