【数学】高斯消元
1. 算法简介
高斯消元法（Gauss–Jordan elimination）是求解线性方程组的经典算法。
例如求解下列方程组：
\( \begin{cases}2x+9y-5z=10\\ 4x+20y+z=24\\ x-2y+3z=8 \end{cases} \)
形式化的，高斯消元可用于求解类似于
\( \begin{cases} a_{1,1} x_1+a_{1,2} x_2+\dots+a_{1,n} x_n = b_1\\ a_{2,1} x_1+a_{2,2} x_2+\dots+a_{2,n} x_n = b_2\\ \dots\\ a_{n,1} x_1+a_{n,2} x_2+\dots+a_{n,n} x_n = b_n\\ \end{cases} \)
的方程组。
2. 算法理论
2.1 高斯消元
先将线性方程组转化为矩阵：
\( \begin{bmatrix} a_{1,1}&a_{1,2}&\cdots&a_{1,n}&b_1\\ a_{2,1}&a_{2,2}&\cdots&a_{2,n}&b_2\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ a_{n,1}&a_{n,2}&\cdots&a_{n,n}&b_n\\ \end{bmatrix} \)
显然为一个 \(n\times (n+1)\) 的矩阵。
我们希望将它消成上三角的形式，如下：
\( \begin{bmatrix} w_{1,1}&w_{1,2}&\cdots&w_{1,n}&B_1\\ 0&w_{2,2}&\cdots&w_{2,n}&B_2\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&a_{n,n}&B_n\\ \end{bmatrix} \)
我们可以枚举对角线，依次将对角线下面的所有数消成 \(0\)。
例如第一次操作利用 \(a_{1,1}\) 将 \(a_{2,1},a_{3,1}\dots a_{n,1}\) 消成 \(0\)。
\( \begin{bmatrix} 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}\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ 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{bmatrix} \)
变成：
\( \begin{bmatrix} 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}\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ 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{bmatrix} \)
依次类推消掉下三角的所有数，消成上三角的形式：
\( \begin{bmatrix} w_{1,1}&w_{1,2}&\cdots&w_{1,n}&B_1\\ 0&w_{2,2}&\cdots&w_{2,n}&B_2\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&a_{n,n}&B_n\\ \end{bmatrix} \)
然后回带，依次相减。消成对角线的形式：
\( \begin{bmatrix} W_{1,1}&0&\cdots&0&C_1\\ 0&W_{2,2}&\cdots&0&C_2\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&W_{n,n}&C_n\\ \end{bmatrix} \)
则答案 \(x_i=\frac{B_i}{W_{i,i}}\).
判断解的情况：
- 若 \(\exists i\) 使得 \(W_{i,i}=0\) 且 \(C_i \not= 0\)，则方程组无解;
- 若 \(\exists i\) 使得 \(W_{i,i}=0\) 且 \(C_i = 0\)，则方程组有无数解;
误差处理：
- 由于误差传递，要将 \(a_{i,i}=0\) 写成 \(fabs(a_{i,i})<0\);
- 为了减小误差，每次将小的数除以大的数再去消元;
2.1 高斯-约旦消元
对于矩阵
\( \begin{bmatrix} a_{1,1}&a_{1,2}&\cdots&a_{1,n}&b_1\\ a_{2,1}&a_{2,2}&\cdots&a_{2,n}&b_2\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ a_{n,1}&a_{n,2}&\cdots&a_{n,n}&b_n\\ \end{bmatrix} \)
枚举对角线 \((i,i)\)，找到 \(a_{i,i}\not=0\) 的一行，交换至第一行，依次消元;
高斯消元为每一次往下消元，而高斯-约旦消元为每一次将除本行的其他行全部消元。
由于每一次都会将除对角线外的所有数消成 \(0\)，所以高斯-约旦消元后只剩下了
\( \begin{bmatrix} W_{1,1}&0&\cdots&0&C_1\\ 0&W_{2,2}&\cdots&0&C_2\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&W_{n,n}&C_n\\ \end{bmatrix} \)
不需要高斯消元的回带过程。
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;
}
