0.前置知识
- 知道如何解三元一次方程组
- 有手，有脑子
1.答案的表示与存储
先解一个方程组：
2x+3y+5z=31
x-4y -z=-6
4x+2y-5z=9
我们把这个方程组写成 机器能读懂 的 表格形式 ：
2 3 5 31
1 -4 -1 -6
4 2 -5 9
第一列代表
x
x
x 的系数，第二列代表
y
y
y 的系数……
注意多出来的第四列是答案的具体数值。
第一行代表式子 1，第二列代表式子 2……
因为老是使用 x , y , z x,y,z x,y,z 表示 n n n 个数不方便，所以我们使用 a 1 , a 2 . . . a n a_1,a_2...a_n a1,a2...an 来表示。
即，上面的表格的意义就是：
第 i ( 1 ≤ i ≤ n ) i(1\le i\le n) i(1≤i≤n) 列代表每个式子中 a i a_i ai 的系数，第 n + 1 n+1 n+1 列表示每个式子的具体数值。（因为存储的时候这个都是普通的数值，所以读者务必辨清）
第 i ( 1 ≤ i ≤ n ) i(1\le i\le n) i(1≤i≤n) 行代表第 i i i 个式子。
接着全文中，使用 a [ i ] [ j ] a[i][j] a[i][j] 表示第 i i i 行第 j j j 列的数值/系数。
2.高斯消元思想
名字这么高级，但是方法很土。具体操作是：加减消元法+代入消元法。
首先，接着上面的系数矩阵：
2 3 5 31
1 -4 -1 -6
4 2 -5 9
然后开始高斯消元。
枚举 i i i ，顺序消除未知数 a i a_i ai 。
然后枚举 j j j ，找到一个绝对值最小不为 0 0 0 的 a [ j ] [ i ] a[j][i] a[j][i] 。
找到那个构成最小 a [ j ] [ i ] a[j][i] a[j][i] 的 j j j 后，将第 j j j 个式子与第 i i i 个式子交换
枚举 j j j ，设 m u l = a [ i ] [ i ] / a [ j ] [ i ] mul=a[i][i]/a[j][i] mul=a[i][i]/a[j][i] ，
将 j j j 式子乘上 m u l mul mul 并且减去 i i i 式子
我们一层层剖析为什么这样做。
首先，枚举
i
i
i ，需要消掉第
i
i
i 个未知数。
然后，选择
j
j
j 使得
a
[
j
]
[
i
]
a[j][i]
a[j][i] 最小，这个其实不是必须的（后面讲为什么）。交换
j
j
j 式子和
i
i
i 式子。
再后来，枚举每一个
j
j
j ，将
j
j
j 式子乘上
a
[
i
]
[
i
]
/
a
[
j
]
[
i
]
a[i][i]/a[j][i]
a[i][i]/a[j][i] ，是因为需要将
j
j
j 式子第
i
i
i 项的系数 调整与
i
i
i 式子第
i
i
i 项的系数一样；
将
j
j
j 式子减去
i
i
i 式子，目的为了消去第
i
i
i 项的系数 。
为什么要选择 绝对值最小而不为
0
0
0 的
a
[
j
]
[
i
]
a[j][i]
a[j][i] 的
j
j
j 呢？因为绝对值最小的
a
[
j
]
[
i
]
a[j][i]
a[j][i] 大致可以构成更多整除的情况，而这样就可以避免掉精度。
不能取
0
0
0 是当然的。
如果还没有理解？那么一层层模拟。
还是以这个方程组为例子：
2 3 5 31
1 -4 -1 -6
4 2 -5 9
首先，消去第一个系数。
1 -4 -1 -6
2 3 5 31
4 2 -5 9
找到最小 a [ 1 ] a[1] a[1] 的式子，调到第一个；
1 -4 -1 -6
1 3/2 5/2 31/2
4 2 -5 9
第二个柿子 × 0.5 \times 0.5 ×0.5 。
1 -4 -1 -6
0 11/2 7/2 43/2
4 2 -5 9
第二个式子减去第一个。
以此类推，最后这个答案的系数矩阵是：（没有对齐，有点丑）
4 0 0 20
0 -4.5 0 -9
0 0 7.611 22.833
然后最后的 n + 1 n+1 n+1 系数除以 第 i i i 行第 i i i 列的数值即可。
就知道你们想看代码。
P3389 【模板】高斯消元法
注释有点少，能看懂就行。
#include<bits/stdc++.h>
#define RI register int
using namespace std;typedef long long LL;const int inf=1073741823;int FL,CH;template<typename T>void in(T&a){for(a=0,FL=1,CH=getchar();!isdigit(CH);CH=getchar())FL=(CH=='-')?-1:1;for(;isdigit(CH);CH=getchar())a=a*10+CH-'0';a*=FL;}template<typename T,typename ...Args>void in(T&a,Args&...args){in(a);in(args...);}
const int maxn=110;const double eps=1e-6;
double a[maxn][maxn];
int n;
int main()
{
cin>>n;
for(RI i=1;i<=n;++i){for(RI j=1;j<=n+1;++j)scanf("%lf",&a[i][j]);}
for(RI rp,i=1;i<=n;++i)
{
rp=i;
for(RI j=i+1;j<=n;++j)if(fabs(a[j][i])>fabs(a[rp][i]))rp=j;
if(fabs(a[rp][i])<=eps)return printf("No Solution"),0;
//寻找绝对值最小的 j
//这个特殊一点，要判无解
if(i!=rp)swap(a[rp],a[i]);
//交换 j 和 i
double div=a[i][i];
for(RI j=1;j<=n;++j)
if(j!=i)
{
div=a[j][i]/a[i][i];
//刚才的 mul
for(RI k=i+1;k<=n+1;++k)
a[j][k]-=a[i][k]*div;
//乘上，再减去
}
}
for(RI i=1;i<=n;++i)
printf("%.2lf\n",a[i][n+1]/a[i][i]);//输出答案
return 0;
}