高斯消元

· · 算法·理论

用于求解线性方程组。

[SDOI2006] 线性方程组

这里不放高斯消元板子是因为其中 No solution 情况分无穷组解和无解两种，没有这题这么详细。

题意：已知 n 元线性一次方程组：

\begin{cases} 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 \\ \cdots \\ a_{n,1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n \end{cases}

其中任意的 a 和 b 都以矩阵方式给定，求解 x_1\dots x_n。

回顾初中知识我们知道了加减消元的过程。具体而言我们会让方程组的某一元的系数一样或变为相反数从而用加减法去掉这一个元。例如这个：

\begin{cases}3x+5y=9\\2x+6y=20\end{cases}

我们让下面的式子 \times\frac{3}{2}，得到：3x+9y=30，这样用上面的式子减去下面的式子就可以得到：-4y=-21，解得 y=\frac{21}{4}，然后再把 y=5.25 带回到之前的式子求出 x=-5.75。

我们尝试模拟这个过程，具体而言：

  1. 枚举 1\to n,我们找到 x_i 最大的系数的那一个方程并把其与当前行交换。最大系数的原因是因为减少误差，详细证明看这里，由于笔者太菜就不证了。
  2. 我们用该方程消去其它方程的 x_i，即乘上 \frac{a_{j,i}}{a_{line,i}}，然后用别的方程减去这个方程。这样最后的结果为仅剩当前方程含有 x_i，而别的 x_i 系数都为 0。
  3. 最后我们就会得到一个只有 a_{i,i} 不为 0 的方程组，x_i=\frac{a_{i,n+1}}{a_{i,i}}。这是方程的唯一节，因为 \frac{a_{i,n+1}}{a_{i,i}} 是固定的。

这个方法看似是很正确的，但是有一个必要前提：每次的元系数都不能为 0。

举个例子：\begin{cases}2x+4y=12\\6x+12y=36\end{cases}。我们尝试用 1 式消 2 式，但是我们发现 2 式消完之后就没有 y 了。这属于无唯一解中的有无穷多解。

如何判断到底是无穷多解还是无解？容易发现若 b_i 为 0 则这行相当于无用的，否则就是无解。因为若不为 0 则会产生 0\times x+0\times y=1 这样的方程，显然是无解。否则如果没有这样的方程即剩余的都是 0=0 那么我们就认为有无穷多组解。

#include <bits/stdc++.h>
using namespace std;

const int _ = 55;
const double eps = 1e-9;

int n;
double a[_][_];

int Gauss() {
  int line = 1;
  for(int i = 1;i <= n;i++) {
    int mx = line;
    for(int j = line + 1;j <= n;j++) {  //找到系数最大的行
      if(fabs(a[j][i]) > fabs(a[mx][i])) mx = j;
    }
    if(fabs(a[mx][i]) < eps) continue; //若误差<1e-9则认为为 0
    for(int j = 1;j <= n + 1;j++) swap(a[line][j],a[mx][j]);
    for(int j = 1;j <= n;j++) { //消元其他的
      if(line == j) continue;
      double tmp = a[j][i] / a[line][i];
      for(int k = i;k <= n + 1;k++) {
        a[j][k] -= a[line][k] * tmp;
      }
    }
    line++;
  }
  for(int i = line;i <= n;i++) {
    if(fabs(a[i][n + 1]) > eps) { //若bi不为0
      return -1;
    }
  }
  return line > n;
}

int main() {
  cin >> n;
  for(int i = 1;i <= n;i++) {
    for(int j = 1;j <= n + 1;j++) {
      cin >> a[i][j];
    }
  }
  int t = Gauss();
  if(t == 1) {
    for(int i = 1;i <= n;i++) {
      cout << 'x' << i << '=' << fixed << setprecision(2) << a[i][n + 1] / a[i][i] << '\n';
    }
  } else {
    cout << t << '\n';
  }
  return 0;
}

P3389 【模板】高斯消元法