介绍
何为高斯消元?就是拿来解多元一次方程(组)的一个算法。
什么是多元一次方程?就是有多个未知数,每一个未知数的次数皆为1的方程。多个多元一次方程构成多元一次方程组。如下是一个三元一次方程组:
{
3
x
1
+
2
x
2
+
x
3
=
3
2
x
1
+
x
2
+
3
x
3
=
6
x
1
+
3
x
2
+
2
x
3
=
9
\begin{cases} 3x_1+2x_2+x_3=3\\ 2x_1+x_2+3x_3=6\\ x_1+3x_2+2x_3=9 \end{cases}
⎩⎪⎨⎪⎧3x1+2x2+x3=32x1+x2+3x3=6x1+3x2+2x3=9
容易通过代入消元法或加减消元法获得唯一解:
{
x
1
=
−
1
x
2
=
2
x
3
=
2
\begin{cases}x_1=-1\\x_2=2\\x_3=2\end{cases}
⎩⎪⎨⎪⎧x1=−1x2=2x3=2
具体题目
题目(只有5个n,m=3的数据,凑合着看看吧)
讲解
还是那个方程:
{
3
x
1
+
2
x
2
+
x
3
=
3
①
2
x
1
+
x
2
+
3
x
3
=
6
②
x
1
+
3
x
2
+
2
x
3
=
9
③
\begin{cases} 3x_1+2x_2+x_3=3\ ①\\ 2x_1+x_2+3x_3=6\ ②\\ x_1+3x_2+2x_3=9\ ③ \end{cases}
⎩⎪⎨⎪⎧3x1+2x2+x3=3 ①2x1+x2+3x3=6 ②x1+3x2+2x3=9 ③
我们可以直接解:
②
∗
3
2
−
*\frac{3}{2}-
∗23−①,③
∗
3
−
*3-
∗3−①得
{
3
x
1
+
2
x
2
+
x
3
=
3
①
−
1
2
x
2
+
7
2
x
3
=
6
④
7
x
2
+
5
x
3
=
24
⑤
\begin{cases} 3x_1+2x_2+x_3=3\ ①\\ -\frac{1}{2}x_2+\frac{7}{2}x_3=6\ ④\\ 7x_2+5x_3=24\ ⑤ \end{cases}
⎩⎪⎨⎪⎧3x1+2x2+x3=3 ①−21x2+27x3=6 ④7x2+5x3=24 ⑤
所以由③
∗
(
−
1
14
)
−
*(-\frac{1}{14})-
∗(−141)−④可以得到
{
3
x
1
+
2
x
2
+
x
3
=
3
①
−
1
2
x
2
+
7
2
x
3
=
6
④
27
7
x
3
=
54
7
⑥
\begin{cases} 3x_1+2x_2+x_3=3\ ①\\ -\frac{1}{2}x_2+\frac{7}{2}x_3=6\ ④\\ \frac{27}{7}x_3=\frac{54}{7}\ ⑥ \end{cases}
⎩⎪⎨⎪⎧3x1+2x2+x3=3 ①−21x2+27x3=6 ④727x3=754 ⑥
于是就可以得到
x
3
=
2
x_3=2
x3=2。
代入④中,得
−
1
2
x
2
+
7
=
6
⇒
x
2
=
2
-\frac{1}{2}x_2+7=6\Rightarrow x_2=2
−21x2+7=6⇒x2=2
代入①中,得
3
x
1
+
4
+
2
=
3
⇒
x
1
=
−
1
3x_1+4+2=3\Rightarrow x_1=-1
3x1+4+2=3⇒x1=−1
于是得解。实际上,我们刚才已经模拟了一个高斯消元的过程。虽然这个过程在数学上看来是很麻烦的,但是计算机并不会觉得麻烦,而且这也是一个计算机中的可行算法。模拟手算中的代入消元法,这就是高斯消元的真谛。
矩阵
先不管矩阵是什么,看着就明白了。
如下,原方程
{
3
x
1
+
2
x
2
+
x
3
=
3
①
2
x
1
+
x
2
+
3
x
3
=
6
②
x
1
+
3
x
2
+
2
x
3
=
9
③
\begin{cases} 3x_1+2x_2+x_3=3\ ①\\ 2x_1+x_2+3x_3=6\ ②\\ x_1+3x_2+2x_3=9\ ③ \end{cases}
⎩⎪⎨⎪⎧3x1+2x2+x3=3 ①2x1+x2+3x3=6 ②x1+3x2+2x3=9 ③可以化为这样一个增广矩阵
[
3
2
1
∣
3
2
1
3
∣
6
1
3
1
∣
9
]
\begin{bmatrix} 3& 2& 1&|&3\\ 2& 1& 3&|&6\\ 1& 3& 1&|&9\\ \end{bmatrix}
⎣⎡321213131∣∣∣369⎦⎤
按照这个道理,可以将
{
3
x
1
+
2
x
2
+
x
3
=
3
①
−
1
2
x
2
+
7
2
x
3
=
6
④
27
7
x
3
=
54
7
⑥
\begin{cases} 3x_1+2x_2+x_3=3\ ①\\ -\frac{1}{2}x_2+\frac{7}{2}x_3=6\ ④\\ \frac{27}{7}x_3=\frac{54}{7}\ ⑥ \end{cases}
⎩⎪⎨⎪⎧3x1+2x2+x3=3 ①−21x2+27x3=6 ④727x3=754 ⑥也转化成一个矩阵,就是
[
3
2
1
∣
3
0
−
1
2
7
2
∣
6
0
0
27
7
∣
54
7
]
\begin{bmatrix} 3& 2& 1&|&3\\ 0& -\frac{1}{2}& \frac{7}{2}&|&6\\ 0& 0& \frac{27}{7}&|&\frac{54}{7}\\ \end{bmatrix}
⎣⎡3002−210127727∣∣∣36754⎦⎤
可以看到,除去最后一行常数,系数构成了一个类似三角形的矩阵,叫做阶梯矩阵(虽然我个人还是习惯叫三角矩阵)
所以我们的目标就是将普通的矩阵化为三角矩阵,这里就是纯粹的模拟了。具体不懂的看注释吧……
#define Abs(x) ( x > 0 ? x : -x )
const double EPS = 1e-6;//这个东西是拿来预防double误差的,一般来说,一个数的绝对值<EPS则这个数为0
int n, m;//n个方程,m个未知数
double a[401][402];//矩阵
double ans[401];//储存答案,即所有x的解
bool fr[401];//自由元,是自由元为0,否则为1
register int now, k;
for(k = 1, now = 1; k <= n && now <= m; k++, now++) {//要遍历每一个方程,每一个未知数
int j = k;
//因为前面的都已经化成了三角矩阵,所以要在还未进行消元的部分,
//即还没有化为三角矩阵的部分进行消元
for(register int i = k + 1; i <= n; i++)
if(Abs(a[j][now]) < Abs(a[i][now])) j = i;
//找出当前这一个未知数中最大的系数,因为这样可以尽可能地避免误差
if(Abs(a[j][now]) < EPS) { k--; continue; }
//如果当前这一个未知数所有的系数都为0,因为最大值为0,则遍历下一个未知数
if(j != k) swap(a[j], a[k]);//不用担心,是不会CE的
//直接把消元的当前这一行放到前面的一行,这样就可以构造一个三角形的矩阵
//处理当前这一行的系数,以构造三角形矩阵
for(register int i = 1; i <= n; i++) {//对每一个方程消元
if(Abs(a[i][now]) < EPS || i == k) continue;//若当前未知数的系数为0,当然就跳过啦
double D = a[i][now] / a[k][now];
//例题是浮点型的,所以用除法消元,如果是整形,就应该用gcd(a[i][now], a[k][now])
for(register int x = now; x <= m + 1; x++) a[i][x] -= D * a[k][x];//消这个未知数
}
}
可以明显看出,高斯消元的复杂度是
O
(
n
3
O(n^3
O(n3+巨大常数
)
)
)的。(啊这真是一个快速的算法)
特殊情况处理
无解
例如有如下方程:
{
x
+
y
=
1
x
+
y
=
2
\begin{cases}x+y=1\\x+y=2\end{cases}
{x+y=1x+y=2
显然冲突,无解。这种情况该如何处理呢?
很显然,如果有这样的一行
[
0
0
0
⋯
0
∣
b
]
(
b
≠
0
)
\begin{bmatrix}0&0&0&\cdots&0&|&b\end{bmatrix}(b\neq0)
[000⋯0∣b](b̸=0),那这个方程显然无解,于是可以用这样的一行来判断一下是否有解:
for(reg int i = k; i <= n; i++)//为什么从k开始?因为前k开始系数都是被消掉的全是0了,这样才满足条件
if(Abs(a[i][m + 1]) > EPS) return !printf("No solution");
//因为有消元了过后,前面全是0,而后面常数项非0的,这种情况当然无解
自由元
自由元是什么?
如下:
{
4
x
1
+
2
x
2
+
x
3
=
3
①
2
x
1
+
x
2
+
x
3
=
6
②
x
3
=
9
③
\begin{cases} 4x_1+2x_2+x_3=3\ ①\\ 2x_1+x_2+x_3=6\ ②\\ x_3=9\ ③ \end{cases}
⎩⎪⎨⎪⎧4x1+2x2+x3=3 ①2x1+x2+x3=6 ②x3=9 ③
可以解出
x
3
=
9
x_3=9
x3=9,但是
x
1
,
x
2
x_1,x_2
x1,x2都没有确定的解,这个就叫做自由元
因为消了
k
k
k个元,所以如果未知数的个数
m
>
k
m>k
m>k,那就可能存在自由元。
if(k <= m) {//存在自由元
for(reg int i = k - 1; i >= 1; i--) {
int Unknown = 0, idx = 0;//统计该方程中未出现的未知数的个数和那个数的下标
for(reg int j = 1; j <= m; j++)//找每一个系数
if(!fr[j] && Abs(a[i][j]) > EPS) Unknown++, idx = j;
//非0,且不是已经有解了的未知数
if(Unknown > 1) continue;//此方程暂时无解
double tmp = a[i][m + 1];
for(reg int j = 1; j <= m; j++)
if(j != idx && Abs(a[i][j]) > EPS) tmp -= a[i][j] * ans[j];//模拟回代不解释
//倒着推回来,进行代入消元
ans[idx] = tmp / a[i][idx];//得到这个未知数的解,还要记住可能存在的系数
fr[idx] = 1;//这个未知数肯定不是自由元啦~
}
}
普通情况
for(reg int i = m; i >= 1; i--) {//m个元,m个含系数的未知式
double tmp = a[i][m + 1];
for(reg int j = i + 1; j <= m; j++)
tmp -= a[i][j] * ans[j];//代入……
ans[i] = tmp / a[i][i];//除以系数,得解
fr[i] = 1;//不是自由元
}//普通的回带不解释
代码总结
#define reg register
#define Abs(a) ( a > 0 ? a : -a )
const double EPS = 1e-6;
int n, m;
double a[401][402];
double ans[401];
bool fr[401];//自由元
inline bool Gauss() {//返回1表示有解,返回0表示无解
register int now, k;
for(k = 1, now = 1; k <= n && now <= m; k++, now++) {//要遍历每一个方程,每一个未知数
int j = k;
//因为前面的都已经化成了三角矩阵,所以要在还未进行消元的部分,
//即还没有化为三角矩阵的部分进行消元
for(register int i = k + 1; i <= n; i++)
if(Abs(a[j][now]) < Abs(a[i][now])) j = i;
//找出当前这一个未知数中最大的系数,因为这样可以尽可能地避免误差
if(Abs(a[j][now]) < EPS) { k--; continue; }
//如果当前这一个未知数所有的系数都为0,因为最大值为0,则遍历下一个未知数
if(j != k) swap(a[j], a[k]);//不用担心,是不会CE的
//直接把消元的当前这一行放到前面的一行,这样就可以构造一个三角形的矩阵
//处理当前这一行的系数,以构造三角形矩阵
for(register int i = 1; i <= n; i++) {//对每一个方程消元
if(Abs(a[i][now]) < EPS || i == k) continue;//若当前未知数的系数为0,当然就跳过啦
double D = a[i][now] / a[k][now];
//例题是浮点型的,所以用除法消元,如果是整形,就应该用gcd(a[i][now], a[k][now])
for(register int x = now; x <= m + 1; x++) a[i][x] -= D * a[k][x];//消这个未知数
}
}
for(reg int i = k; i <= n; i++)//为什么从k开始?因为前k开始系数都是被消掉的全是0了,这样才满足条件
if(Abs(a[i][m + 1]) > EPS) return 0;
//因为有消元了过后,前面全是0,而后面常数项非0的,这种情况当然无解
if(k <= m) {//存在自由元
for(reg int i = k - 1; i >= 1; i--) {
int Unknown = 0, idx = 0;//统计该方程中未出现的未知数的个数和那个数的下标
for(reg int j = 1; j <= m; j++)//找每一个系数
if(!fr[j] && Abs(a[i][j]) > EPS) Unknown++, idx = j;
//非0,且不是已经有解了的未知数
if(Unknown > 1) continue;//此方程暂时无解
double tmp = a[i][m + 1];
for(reg int j = 1; j <= m; j++)
if(j != idx && Abs(a[i][j]) > EPS) tmp -= a[i][j] * ans[j];//模拟回代不解释
//倒着推回来,进行代入消元
ans[idx] = tmp / a[i][idx];//得到这个未知数的解,还要记住可能存在的系数
fr[idx] = 1;//这个未知数肯定不是自由元啦~
}
return 1;
}
for(reg int i = m; i >= 1; i--) {//m个元,m个含系数的未知式
double tmp = a[i][m + 1];
for(reg int j = i + 1; j <= m; j++)
tmp -= a[i][j] * ans[j];//代入……
ans[i] = tmp / a[i][i];//除以系数,得解
fr[i] = 1;//不是自由元
}//普通的回带不解释
return 1;
}
int main() {
scanf("%d%d", &n, &m);
for(reg int i = 1; i <= n; i++)
for(reg int j = 1; j <= m + 1; j++)
scanf("%lf", &a[i][j]);
bool sol = Gauss();
if(!sol) {
printf("No solution");
return 0;
}
for(reg int i = 1; i <= m; i++) {
if(!fr[i]) printf("X[%d] not determined\n", i);
else {
printf("X[%d] = ", i);
if(Abs(ans[i]) < EPS)
ans[i] = 0;
printf("%.4lf\n", ans[i]);
}
}
}
附上:上方代码没有头文件的原因
#include<map>
#include<set>
#include<list>
#include<queue>
#include<deque>
#include<stack>
#include<ctime>
#include<cmath>
#include<vector>
#include<bitset>
#include<cctype>
#include<string>
#include<sstream>
#include<fstream>
#include<cstdlib>
#include<cstring>
#include<climits>
#include<iomanip>
#include<iostream>
#include<algorithm>
#include<tr1/unordered_map>
using namespace std;
using namespace tr1;