高斯消元

高斯消元

板题 [SDOI2006] 线性方程组（别问我为什么不是【模板】高斯消元法，这个太**了）

思路

首先需要引入一个定义增广矩阵。

所以一个 \(n\) 元线性方程组就可以抽象成一个矩阵，\(a\) 为系数，\(b\) 为方程的常数项:

\(\begin{bmatrix}a_{11}&\cdots&a_{1n}&b_{1}\\\vdots&\ddots&\vdots\\a_{n1}&\cdots&a_{nn}&b{n}\end{bmatrix}\)

我们知道对于任意一个一元一次方程，我们都需要将其解成一个类似于 \(ax=b\) 的形式，然后变成 \(x=\frac{b}{a}\)，所以对于一个方程组，其实是一样的也需要解成一堆七七八八的东西，然后将 \(a\) 变成一个单位矩阵，此时右边的 \(b\) 就是对应行的主元的解了。

为了达成这个目标，我们有三种行变换操作，目的是为了在保证主元解不变的情况下进行变化：

1. 交换两行，这个很显然不会变答案
2. 将一行整体缩放一个数，这个很显然可以用等式的性质二来证明（欸，要证明吗，好像本来就满足）
3. 将一行叠加到另一行去，这个可以用等式的性质一来证明（嘿，好像不用证明）

好，此时我们就可以快速口胡一份代码了，我们可以对于每一行都将其下面的全部都变成 \(0\)，因为动了上面的，那么上面的也得行变换，下面的也都得行变换了，注意此时是要用行变换的！！！做完以后就发现最后一行刚好是一个 \(a_{nn}=b_n\) 的东西，所以可以将 \(x_n\) 求出，然后倒着带入求解其它的元即可。

到这里，我们就解决完了所有等式有唯一解的情况了，接下来就要说一些特殊的情况：

1. 等式当前行主元系数为 \(0\)，这个又分两种情况，一个呢是下面的当前主元系数不为 \(0\)，那么此时只要用行变换把它弄上来即可，至于为什么不动上面的，理由和上面一样，二个呢是下面的当前主元系数也为 \(0\)，那么这个主元，就变成了无论是什么，整个式子都成立，所以就成了无唯一解，可此时我们还不能下定结论，因为下面还有更离谱的事情，不过这就告诉我们，我们不能单纯的只维护当前行，应该把当前列也维护了，因为如果当前行的当前主元无唯一解时，我们就需要让当前列往后走，找到一个有唯一解的主元给当前行，那么跑到最后时，当前行和当前列不一定都跑完了所有的方程。
2. 等式接完后有剩余等式，这个时候我们会发现一件事情，就是这些等式与上面的等式是线性相关的，所以解不解都没关系，但是，这里因为已经被上面削成了 \(0\)，所以有可能会出现 \(0\ne 0\) 的情况，所以要判无解，然后再判无唯一解。

好，到这里就全算完了。

code

#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std;

const int MaxN = 110;
const double eps = 1e-4;

double a[MaxN][MaxN], ans[MaxN];
int 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 i, j;
  for (i = 1, j = 1; i <= n && j <= n;) {
    int u = 0;
    for (int k = n; k >= i; k--) {
      if (fabs(a[k][j]) > eps) {
        u = k;
      }
    }
    if (!u) {
      j++;
      continue;
    } 
    for (int k = 1; k <= n + 1; k++) {
      swap(a[i][k], a[u][k]);
    }
    for (int k = i + 1; k <= n; k++) {
      double d = a[k][j] / a[i][j];
      for (int l = j; l <= n + 1; l++) {
        a[k][l] -= a[i][l] * d;
      }
    }
    i++, j++;
  }
  for (int k = i; k <= n; k++) {
    if (fabs(a[k][n + 1]) > eps) {
      cout << "-1" << endl;
      return 0;
    }
  }
  if (j < n || i < n) {
    cout << "0" << endl;
    return 0;
  }
  for (int i = n; i >= 1; i--) {
    double sum = 0;
    for (int j = 1; j <= n; j++) {
      sum += ans[j] * a[i][j];
    }
    if (fabs(a[i][i]) < eps) {
      cout << "0" << endl;
      return 0;
    }
    ans[i] = (a[i][n + 1] - sum) / a[i][i];
  }
  for (int i = 1; i <= n; i++) {
    cout << fixed << setprecision(2) << "x" << i << "=" << ans[i] << endl;
  }
  return 0;
}

未完待续……