[算法] 高斯消元详解

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;
}
上一篇:NOIP 模拟 $14\; \text{影魔}$


下一篇:NOIP 模拟 $26\; \rm 降雷皇$