[模板] 高斯消元

求解方程

高斯消元法，是一种求解线性方程组的算法.
何为线性方程组?
就是这样一组方程:
$$\{\begin{array}{l}{a}_{1_1}{x}_1+a_{1_2}{x}_2+\cdots +a_{1_n}{x}_n=A_1\\ a_{2_1}{x}_1+a_{2_2}{x}_2+\cdots +a_{2_n}{x}_n=A_2\\ \dots \\ a_{n_1}{x}_1+a_{n_2}{x}_2+\cdots +a_{n_n}{x}_n=A_n\end{array}$$
所有的未知量的指数都是$1$.
这样的方程亦可以用这样的矩阵表示:
$$\left[\begin{array}{ccccc}{a}_{1_1}& a_{1_2}& \dots & a_{1_n}& A_1\\ a_{2_1}& a_{2_2}& \dots & a_{2_n}& A_2\\ \dots \dots & \dots \dots & \dots \dots & \dots \dots & \dots \dots \\ a_{n_1}& a_{n_2}& \dots & a_{n_n}& A_n\end{array}\right]$$
并称这样的矩阵为增广矩阵.
高斯消元法就是在增广矩阵上操作的.

思想

高斯消元法的思路和直接笔算的思路一致，具体(步骤)如下:

1. 从方程集合中选择一个方程.
2. 将这个方程的一个未知量($x_k$)变为$1$.
3. 用这个方程将其他所剩方程的$x_k$均变换为$0$.
4. 将这个方程从方程集合中删去.

变换规则(操作)如下:

• (1)可以将矩阵一行所有元素同时乘一个非零数.(其实也可以乘 $0$，但谁会那样做呢 ? 我会)
• (2)将矩阵一行的元素加到另一行上.
• (3)将矩阵的两行交换.

这样可能比较难以理解，可以结合下例:
$$\left[\begin{array}{cccc}1& 2& -1& -6\\ 2& 1& -3& -9\\ -1& -1& 2& 7\end{array}\right]$$
对应方程组:
$$\{\begin{array}{l}{x}_1+2x_2-x_3=-6\\ 2x_1+x_2-3x_3=-9\\ -x_1-x_2+2x_3=7\end{array}$$
这个方程的解是$\begin{array}{c}{x}_1=1,x_2=-2,x_3=3\end{array}$.(什么?不规范?)
那么观察高斯消元的方法:
我们先通过操作(1)将这个矩阵变为:
$$\left[\begin{array}{cccc}1& 2& -1& -6\\ -1& -0.5& 1.5& 4.5\\ -1& -1& 2& 7\end{array}\right]$$
然后通过操作(2)变为:
$$\left[\begin{array}{cccc}1& 2& -1& -6\\ 0& 1.5& 0.5& -1.5\\ 0& 1& 1& 1\end{array}\right]$$
现在第一行已经没有用了，将其从方程集合中删去:
注意，方程 1 在原增广矩阵中还是存在的.
$$\left[\begin{array}{cccc}0& 1.5& 0.5& -1.5\\ 0& 1& 1& 1\end{array}\right]$$
然后通过操作(1)变为:
$$\left[\begin{array}{cccc}0& 1.5& 0.5& -1.5\\ 0& -1.5& -1.5& -1.5\end{array}\right]$$
然后通过操作(2)和步骤(4)变为:
$$\left[\begin{array}{cccc}0& 0& -1& -3\end{array}\right]$$
在执行操作(1)，变为:
$$\left[\begin{array}{cccc}0& 0& 1& 3\end{array}\right]$$
然后观察原增广矩阵，变成了这样:
$$\left[\begin{array}{cccc}1& 2& -1& -6\\ 0& 1.5& 0.5& -1.5\\ 0& 0& 1& 3\end{array}\right]$$
在经过(1)次操作(1)得:
$$\left[\begin{array}{cccc}1& 2& -1& -6\\ 0& 1& \frac{1}{3}& -1\\ 0& 0& 1& 3\end{array}\right]$$
这个矩阵称为阶梯形矩阵，在通过一波操作之后变成了:
$$\left[\begin{array}{cccc}1& 0& 0& 1\\ 0& 1& 0& -2\\ 0& 0& 1& 3\end{array}\right]$$
对应的就是上面的解.
这个矩阵就是简化阶梯形矩阵.
高斯消元法的本质就是将一般的增广矩阵变换为简化阶梯形矩阵.
别想用 $int$，$double$就好.

特殊情况

当然，这只是“一路顺风”的情况.还有可能“一路逆风”
比如，平常我们解的方程组有可能无解，也有可能有无数组解.
那么在增广矩阵中怎么判断？
首先，如果发现了这样的情况:
$$\left[\begin{array}{cccccc}0& 0& 0& \dots & 0& k\end{array}\right]$$
就肯定不是正常情况了.
很容易想到如下的讨论形式:
$$\{\begin{array}{ll}Infinitenumberofsolutions& k=0\\ Nosolution& others(k̸=0)\end{array}$$
当然，如果坚持做完高斯消元的话，这个矩阵可能会变成这个样子(要用到所有步骤和操作):
$$\left[\begin{array}{cccc}1& 0& 0& k_1\\ 0& 0& 0& 0\\ 0& 0& 0& k_3(k̸=0)\end{array}\right]$$
很显然，这个方程是无解的(别的解再多这个坎就是过不去).
所以我们又将这个搞定了.
顺便说一下，对于$0=0$这样的$INFSolutions$的情况，那些随便取值而不会影响结果的(有无限种取值的)为自由元，只能根据自由元取值的(受自由元约束的)为主元.
举个例子:
$$\{\begin{array}{l}{x}_1=4-x_2\\ x_3=1\end{array}$$
其中$x_2$为自由元，$x_1,x_3$为主元.
或者变换一下方程(1)，$x_1$就会变为自由元，$x_2,x_3$就是主元.

一系列的判断

众所周知我很懒，希望能简化步骤.
这里不想看可以跳过.
所以，给出以下几个判断:

1. 我们可以每次用这个方程的非零系数中最高的那一项进行操作
   主要是因为我们懒，每次反正都要删去系数非零的最高项，顺序可以不一.
2. 当我们发现有一项所有方程都没有时，就可以准备进入$INFSolutions$和 $NoSolution$的判断了.
   很简单，这项取什么都不会影响结果.
   但是万一发生了上面矩阵的情况?所以只能准备，不能直接判断.
3. 当一行所有系数都为$0$时，直接返回$INFSolutions$或$NoSolution$.
   这无需进行讨论…吧.

时间复杂度

选择各项$O(n)$.
选择合适的方程(行)$O(n)$.
处理其他行$O(n)$.
总的时间复杂度$O(n^3)$.
还是很暴力优秀的.

Code

题目传送门
话说这题有点怪…

#include<bits/stdc++.h>
using namespace std;
const int N=105;
int n;
double a[N][N];
bool iszero(double x)
{
    return abs(x)<1e-6;
}
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]);
    }
    for (int i=1;i<=n;i++)
    {
        int p;
        for (p=i;p<=n;p++)
            if (!iszero(a[i][p])) break;
        if (p>n)
        {
            if (!iszero(a[i][n+1])) printf("No Solution");
            else printf("No Solution");//printf("INF Solutions");
            return 0;
        }
        for (int j=n+1;j>=p;j--)
            a[i][j]/=a[i][p];
        for (int j=i+1;j<=n;j++)
            for (int k=n+1;k>=p;k--)
                a[j][k]-=a[j][p]*a[i][k];
    }
    for (int i=n;i>=1;i--)
        for (int j=i-1;j>=1;j--)
        {
            a[j][n+1]-=a[j][i]/a[i][i]*a[i][n+1];
            a[j][i]=0;
        }
    for (int i=1;i<=n;i++) printf("%.2lf\n",a[i][n+1]);
    return 0;
}

感谢奆老关注 qwq ?