CSP考了高斯消元，后悔没学，这就来给他学了----苦痛

img

高斯消元

简单来说，高斯消元是一种计算线性方程组的方法

那么说到线性方程组，就不得不提到行列式以及矩阵，而高斯消元正是建立在矩阵或行列式上的一种计算线性方程组的方法，先说矩阵和行列式。

在数学中，矩阵（Matrix）是一个按照长方阵列排列的复数或实数集合，最早来自于方程组的系数及常数所构成的方阵

比如说这样：

\(\begin{vmatrix} 1 & 2 \\ 3 & 4 \\8 & 9 \end{vmatrix}\)

如果行列式中某一行等于两个行列式之和，则这个行列式等于两个行列式之和，这两个行列式分别以这两组数为改行，其余各行相同

\({\begin{vmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{vmatrix}}={\begin{vmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\b_{21}&b_{22}&\cdots&b_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{vmatrix}}+{\begin{vmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\c_{21}&c_{22}&\cdots&c_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{vmatrix}}\)

把一行乘以某数加到另一行，行列式的值不变

高斯消元的重点

若有\(b_{2i}=a_{kj}*\mu\)，则有

\({\begin{vmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{vmatrix}}={\begin{vmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\ b_{21}&b_{22}&\cdots&b_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{vmatrix}}\)

步入正题：

高斯消元是我们求解线性方程组的一种方法，其实它本质上是我们计算行列式的一种方法，略做扩展后，它也可以在矩阵里消元，再利用矩阵来计算线性方程组

计算行列式是有公式的，但复杂度极高

先引入一种特殊的行列式：三角行列式，它和三角矩阵的定义基本一致

\(\begin{vmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ 0 & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a_{nn} \end{vmatrix}\)

这个行列式的值就是主对角线的乘积

实际上高斯消元也就是完成的这个工作，我们通过最后一条性质就可以对整个行列式进行反复的消元，最终达到上三角的形式

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
double z[114][114];//二维数组模拟矩阵 
int n;//n阶矩阵 
 double Guass1()
{
	double ans=1.00;//最终要输出的结果 
	for( LL i=1;i<=n;i++)
	{
		for( LL j=i+1;j<=n;j++)
		{
			double temp=z[j][i]/z[i][i];//计算要计算的比值 
			for( LL u=1;u<=n;u++) z[j][u]-=temp*z[i][u];//把要消去的那一行用当前行消(减) 
		}
	}
	for( LL i=1;i<=n;i++) ans*=z[i][i];//计算主对角线的乘积之和，也就是上三角时行列式的值 
	return  ans;//返回答案 
}

根据这个代码我们可大约得出，这个算法的时间复杂度约为O(\(n^3\))，虽然还是不小，但是相比于指数级还是小了很多

但是这个1.0的Code并不完美，或者说是有许多的漏洞，主要有以下两点：

  1. 若要消去的元素本来就是0时，这个算比值的方法就无法成立了
  2. 精度损失很大

我们先来针对第一个问题进行修改

根据我们前面所讲到的性质Ⅱ“将行列式中的两行（列）对换，行列式变号”，我们可以把对应元素为不为0的那一行和当前行进行调换，然后将行列式取负即可

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
dd a[514][514];
int n;//同1.0 
double Guass1Plus()
{
	double ans=1.00;
	for( LL i=1;i<=n;i++)
	{
		for( LL j=i+1;j<=n;j++)
		{
			if(fabs(z[j][i])>1e-8)//如果当前行的对应元素不为0，把它交换上来 
			{
				ans=-ans;//因为我们交换了两行，所以行列式变号 
				for( LL u=1;u<=n;u++) swap(z[i][u],a[j][u]);//交换这两行 
				break;//当前行已经交换完了 
			}
		}
		if(fabs(z[i][i])<=1e-8) return 0.00;//若还是0，说明这一列都是0，
        //根据前面说到的性质Ⅲ，这个行列式值就为0 
		for (LL j=i+1;j<=n;j++)
		{
			double temp=z[j][i]/z[i][i];
			for(R LL u=1;u<=n;u++) z[j][u]-=temp*z[i][u];//和刚才1.0一样的进行消元 
		}
	}
	for( LL i=1;i<=n;i++) ans*=z[i][i];//计算行列式的值 
	return ans;
}

但是这样之后我们损失精度最严重的地方还是没有解决

下面我们来分析如何解决精度问题

A1:作为一个由整数元素构成的行列式，其值也必定是整数，但是我们这样算必定得出的时实数，而实数在运算的过程当中就不可避免的会出现精度损失

Q2:为什么这个精度损失大？

A2:因为我们没有做任何的处理，直接就让两个元素求比值，这样做的z[i][i]的大小是不确定的，精度损失也就相应的很大

因此我们若想在不改变算法的前提下优化精度问题，就想办法把每列绝对值大的弄到上面来

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
 double Guass1PlusPlus()
{
	double ans=1.00;
	for( LL i=1;i<=n;i++)
	{
		for( LL j=i+1;j<=n;j++)
		{
			if(fabs(z[j][i])>fabs(z[i][i]))//若当前j遍历到的这个元素绝对值较大，交换两行 
			{
				ans=-ans;//交换两行行列式变号 
				for( LL u=1;u<=n;u++) swap(z[a][c],z[b][c]);
			}
			if(fabs(z[i][i])<=1e-8) return 0.00;//若还是0，说明此列都为0，那么此行列式为0 
			for( LL j=i+1;j<=n;j++)
			{
				double temp=z[j][i]/z[i][i];
				for(R LL u=1;u<=n;u++) z[j][i]-=temp*z[i][i];
			} 
		}
	}
	for( LL i=1;i<=n;i++) ans*=z[i][i];
return ans;
}

那么这个解决了两个问题后的这种算法叫做主元高斯消元法

本质上，精度问题没有解决

为了解决掉这个问题，下面将引出辗转相消高斯消元法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
int Guass2()
{
	int ans=1;	
	for(R int i=1;i<=n;i++)
	{
		for(R int j=1;j<=n;j++)
		{
			while(z[j][i]!=0) 
			{
				int temp=z[i][i]/z[j][i];
				for( int u=1;u<=n;u++) z[i][u]-=z[j][u]*temp;
				for( int u=1;u<=n;u++) swap(z[i][u],z[j][u]);
				ans=-ans;
			}
		}
	}
	for( int i=1;i<=n;i++) ans*=z[i][i];
	return ans;
}

模板：

【模板】高斯消元法

题目背景

如果想要更好地测试高斯消元算法模板请在通过此题后尝试通过 SDOI2006 线性方程组 这一题。

题目描述

给定一个线性方程组，对其求解。

\[ \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}\]

输入格式

第一行，一个正整数 \(n\)。

第二至 \(n+1\) 行，每行 \(n+1\) 个整数，为 $ a_1, a_2, ,a_n$ 和 \(b\)，代表一组方程。

输出格式

共 \(n\) 行，每行一个数，第 \(i\) 行为 \(x_i\)（四舍五入保留 \(2\) 位小数）。

如果不存在唯一解或无解，在第一行输出 No Solution.

样例 #1

样例输入 #1

1
2
3
4
3
1 3 4 5
1 4 7 3
9 3 2 2

样例输出 #1

1
2
3
-0.97
5.18
-2.39

提示

本题 special judge 用于处理可能由于浮点数问题输出 -0.00 的情况。若某个 \(x_i\) 的解四舍五入后是 0.00，那么你的程序输出 -0.00 和输出 0.00 都是正确的。

数据范围：$1 n , | a_i | ^4 , |b | ^4 $。保证数据若有解则所有解均满足 \(|x_i|\le 10^3\)，且 \(x_i\pm 10^{-6}\) 和 \(x_i\) 四舍五入后的结果相同（即不会因为较小的精度误差导致四舍五入后的结果不同）。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include<bits/stdc++.h>
using namespace std;
using ll=long long;
#define rep(i,a,n) for(int i=a;i<=n;i++)
#define frep(i,a,n) for(int i=a;i>=n;i--)
#define int long long
#define PII pair<int,int>
#define lowbit(x) (x&(-x))
const int mod=1e9+7;
const double pai=acos(-1.0);
#define ios ios::sync_with_stdio(0); cin.tie(0),cout.tie(0);
#define LF(x)   fixed<<setprecision(x)
const int N=102;
int n;
double z[N][N],ans[N];
bool guass(){
	for(int i=1;i<=n;i++){
		for(int j=i+1;j<=n;j++)
			if(fabs(z[j][i])>fabs(z[i][i]))
				for(int u=1;u<=n+1;u++) 
                swap(z[i][u],z[j][u]);
		if(fabs(z[i][i])<=1e-8) return 0;
		for(int j=i+1;j<=n;j++){
			double k=z[j][i]/z[i][i];
			for(int u=1;u<=n+1;u++) 
            z[j][u]-=z[i][u]*k;
		}
	}
	return 1;
}

signed main()
{
	cin>>n;
	for(int i=1;i<=n;i++)for(int j=1;j<=n+1;j++)scanf("%lf",&z[i][j]);
	if(!guass()){cout<<"No Solution";return 0;}
	for(int i=n;i>=1;i--)
	{
		double x=0.00;
		for(int j=i+1;j<=n;j++){
			x+=ans[j]*z[i][j];
		}
		ans[i]=(z[i][n+1]-x)/z[i][i];
	}
	for(int i=1;i<=n;i++){
		printf("%0.2lf\n",ans[i]);
	}
	return  0;
}

入门到入土 | 高斯消元 - HerikoDeltana - 博客园 (cnblogs.com)

网课学习笔记——矩阵行列式与高斯消元。 - hyl天梦 - 博客园 (cnblogs.com)