【数学】高斯消元

1. 算法简介#

高斯消元法（Gauss–Jordan elimination）是求解线性方程组的经典算法。

例如求解下列方程组：

$\{\begin{array}{l}2x+9y-5z=10\\ 4x+20y+z=24\\ x-2y+3z=8\end{array}$

形式化的，高斯消元可用于求解类似于

$\{\begin{array}{l}{a}_{1,1}{x}_1+a_{1,2}{x}_2+\cdots +a_{1,n}{x}_n=b_1\\ a_{2,1}{x}_1+a_{2,2}{x}_2+\cdots +a_{2,n}{x}_n=b_2\\ \dots \\ a_{n,1}{x}_1+a_{n,2}{x}_2+\cdots +a_{n,n}{x}_n=b_n\end{array}$

的方程组。

2. 算法理论#

2.1 高斯消元#

先将线性方程组转化为矩阵：

$\left[\begin{array}{ccccc}{a}_{1,1}& a_{1,2}& \cdots & a_{1,n}& b_1\\ a_{2,1}& a_{2,2}& \cdots & a_{2,n}& b_2\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ a_{n,1}& a_{n,2}& \cdots & a_{n,n}& b_n\end{array}\right]$

显然为一个 $n\times (n+1)$ 的矩阵。

我们希望将它消成上三角的形式，如下：

$\left[\begin{array}{ccccc}{w}_{1,1}& w_{1,2}& \cdots & w_{1,n}& B_1\\ 0& w_{2,2}& \cdots & w_{2,n}& B_2\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & a_{n,n}& B_n\end{array}\right]$

我们可以枚举对角线，依次将对角线下面的所有数消成 $0$。

例如第一次操作利用 $a_{1,1}$ 将 $a_{2,1},a_{3,1}\dots a_{n,1}$ 消成 $0$。

$\left[\begin{array}{ccccc}{a}_{1,1}& a_{1,2}& \cdots & a_{1,n}& b_1\\ a_{2,1}-\frac{a_{2,1}}{a_{1,1}}\cdot a_{1,1}& a_{2,2}-\frac{a_{2,1}}{a_{1,1}}\cdot a_{1,1}& \cdots & a_{2,n}-\frac{a_{n,1}}{a_{1,1}}\cdot a_{1,1}& b_2-\frac{a_{2,1}}{a_{1,1}}\cdot a_{1,1}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ a_{n,1}-\frac{a_{n,1}}{a_{1,1}}\cdot a_{1,1}& a_{n,2}-\frac{a_{n,1}}{a_{1,1}}\cdot a_{1,1}& \cdots & a_{n,n}-\frac{a_{n,1}}{a_{1,1}}\cdot a_{1,1}& b_n-\frac{a_{n,1}}{a_{1,1}}\cdot a_{1,1}\end{array}\right]$

变成：

$\left[\begin{array}{ccccc}{a}_{1,1}& a_{1,2}& \cdots & a_{1,n}& b_1\\ 0& a_{2,2}-\frac{a_{2,1}}{a_{1,1}}\cdot a_{1,1}& \cdots & a_{2,n}-\frac{a_{n,1}}{a_{1,1}}\cdot a_{1,1}& b_2-\frac{a_{2,1}}{a_{1,1}}\cdot a_{1,1}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& a_{n,2}-\frac{a_{n,1}}{a_{1,1}}\cdot a_{1,1}& \cdots & a_{n,n}-\frac{a_{n,1}}{a_{1,1}}\cdot a_{1,1}& b_n-\frac{a_{n,1}}{a_{1,1}}\cdot a_{1,1}\end{array}\right]$

依次类推消掉下三角的所有数，消成上三角的形式：

$\left[\begin{array}{ccccc}{w}_{1,1}& w_{1,2}& \cdots & w_{1,n}& B_1\\ 0& w_{2,2}& \cdots & w_{2,n}& B_2\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & a_{n,n}& B_n\end{array}\right]$

然后回带，依次相减。消成对角线的形式：

$\left[\begin{array}{ccccc}{W}_{1,1}& 0& \cdots & 0& C_1\\ 0& W_{2,2}& \cdots & 0& C_2\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & W_{n,n}& C_n\end{array}\right]$

则答案 $x_i=\frac{B_i}{W_{i,i}}$.

判断解的情况：

1. 若 $\mathrm{\exists }i$ 使得 $W_{i,i}=0$ 且 $C_i\ne 0$，则方程组无解;
2. 若 $\mathrm{\exists }i$ 使得 $W_{i,i}=0$ 且 $C_i=0$，则方程组有无数解;

误差处理：

1. 由于误差传递，要将 $a_{i,i}=0$ 写成 $fabs(a_{i,i})<0$;
2. 为了减小误差，每次将小的数除以大的数再去消元;

2.1 高斯-约旦消元#

对于矩阵

$\left[\begin{array}{ccccc}{a}_{1,1}& a_{1,2}& \cdots & a_{1,n}& b_1\\ a_{2,1}& a_{2,2}& \cdots & a_{2,n}& b_2\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ a_{n,1}& a_{n,2}& \cdots & a_{n,n}& b_n\end{array}\right]$

枚举对角线 $(i,i)$，找到 $a_{i,i}\ne 0$ 的一行，交换至第一行，依次消元;

高斯消元为每一次往下消元，而高斯-约旦消元为每一次将除本行的其他行全部消元。

由于每一次都会将除对角线外的所有数消成 $0$，所以高斯-约旦消元后只剩下了

$\left[\begin{array}{ccccc}{W}_{1,1}& 0& \cdots & 0& C_1\\ 0& W_{2,2}& \cdots & 0& C_2\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & W_{n,n}& C_n\end{array}\right]$

不需要高斯消元的回带过程。

3. 算法实现#

luoguP3389 【模板】高斯消元法#

高斯消元模板

#include<bits/stdc++.h>
#define int long long
#define For(i,l,r) for(int i=l;i<=r;++i)
#define FOR(i,r,l) for(int i=r;i>=l;--i)
#define eqs 1e-7

using namespace std;

const int N = 1e3 + 10;

int n;

double a[N][N];

signed main() {
  ios::sync_with_stdio(0);
  cin.tie(0), cout.tie(0);
  cin >> n;
  For(i,1,n) {
    For(j,1,n+1) {
      cin >> a[i][j];
    }
  }
  for (int i = 1; i <= n - 1; ++i) {
    For(j,i+1,n) if(fabs(a[i][i]) < fabs(a[j][i])) swap(a[i], a[j]);
    if(fabs(a[i][i]) < eqs) continue;
    For(j,i+1,n) {
      double w = a[j][i] / a[i][i];
      For(l,i,n+1) {
        a[j][l] -= w * a[i][l];
      }
    }
  }
  for (int i = n; i >= 2; --i) {
    if(fabs(a[i][i]) < eqs) continue;
    FOR(j,i-1,1) {
      double w = a[j][i] / a[i][i];
      For(l,i,n+1) {
        a[j][l] -= w * a[i][l];
      }
    }
  }
  For(i,1,n) {
    int cnt0 = 0;
    For(j,1,n) {
      cnt0 += (a[i][j] == 0);
    }
    if(cnt0 == n) return puts("No Solution"), 0;
  }
  For(k,1,n) {
    printf("%.2lf\n", a[k][n+1] / a[k][k]);
  }
  return 0;
}

/*
9 3       2         2 
0 3.66667 6.77778   2.77778 
0 0       -1.15152  2.75758 


9 0       2         -13.5526 
0 3.66667 0         19.0088 
0 0       -1.15152  2.75758
 */

luoguP2455 [SDOI2006] 线性方程组#

高斯-约旦模板。

高斯消元的华点在于这样一组样例：

4
0 0 1 1 2
0 0 1 0 1
0 0 0 1 1
0 0 0 0 0

0

直接用高斯消元模板会出现对角线全 $0$ 而无法消元的情况。

#include<bits/stdc++.h>
#define ll long long
#define For(i,l,r) for(int i=l;i<=r;++i)
#define FOR(i,r,l) for(int i=r;i>=l;--i)
#define eps 1e-9

using namespace std;

const int N = 55;

int n;

double a[N][N];

bool st[N];

signed main() {
  ios::sync_with_stdio(0);
  cin.tie(0), cout.tie(0);
  cin >> n;
  For(i,1,n) {
    For(j,1,n+1) {
      cin >> a[i][j];
    }
  }
  For(i,1,n) {
    For(j,1,n) if(fabs(a[i][i]) < fabs(a[j][i]) && !st[j]) swap(a[i], a[j]), swap(st[i], st[j]);
    if(fabs(a[i][i]) < eps) continue;
    st[i] = 1;
    For(j,1,n) {
      if(i == j) continue;
      double w = a[j][i] / a[i][i];
      For(l,1,n+1) {
        a[j][l] -= w * a[i][l];
      }
    }
  }
  For(i,1,n) if(fabs(a[i][i]) < eps && fabs(a[i][n+1]) > eps) return puts("-1"), 0; 
  For(i,1,n) if(fabs(a[i][i]) < eps && fabs(a[i][n+1]) < eps) return puts("0"), 0;
  For(k,1,n) {
    printf("x%d=%.2lf\n", k, a[k][n+1] / a[k][k]);
  }
  return 0;
}