高斯消元 O ( n 3 ) O(n^3) O(n3)
对于
n
n
n个未知数
n
n
n个方程构成的线性方程组:
{
a
11
x
1
+
a
12
x
2
+
⋯
+
a
1
n
x
n
=
b
1
a
21
x
1
+
a
22
x
2
+
⋯
+
a
2
n
x
n
=
b
2
⋯
a
n
1
x
1
+
a
n
2
x
2
+
⋯
+
a
n
n
x
n
=
b
n
\left\{ \begin{array}{c} a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n=b_1 \\ a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n=b_2 \\ \cdots \\ a_{n1}x_1+a_{n2}x_2+\cdots+a_{nn}x_n=b_n \end{array} \right.
⎩⎪⎪⎨⎪⎪⎧a11x1+a12x2+⋯+a1nxn=b1a21x1+a22x2+⋯+a2nxn=b2⋯an1x1+an2x2+⋯+annxn=bn
其对应的增广矩阵为:
[
a
11
a
12
⋯
a
1
n
b
1
a
21
a
22
⋯
a
2
n
b
2
⋮
⋮
⋱
⋮
⋮
a
n
1
a
n
2
⋯
a
n
n
b
n
]
\left[ \begin{array}{ccccc} a_{11} & a_{12} & \cdots & a_{1n} & b_1 \\ a_{21} & a_{22} & \cdots & a_{2n} & b_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} & b_n \end{array} \right]
⎣⎢⎢⎢⎡a11a21⋮an1a12a22⋮an2⋯⋯⋱⋯a1na2n⋮annb1b2⋮bn⎦⎥⎥⎥⎤
用初等行列变化,我们要把这个方程的增广矩阵消成上三角形(唯一解),或者上阶梯型(无解或无穷解)。
矩阵有3种初等行列变化:
- 把某一行乘一个非0的数 (方程的两边同时乘上一个非0数不改变方程的解);
- 交换某两行 (交换两个方程的位置);
- 把某行的若干倍加到另一行上去 (把一个方程的若干倍加到另一个方程上去);
算法步骤:
从左往右枚举每一列 c c c:
- 找到当前列绝对值最大的元素所在的行;
- 把这一行换到最上面。这里的最上面指的是还未被消成阶梯型的行中的最上面,不是整个矩阵的最上面一行;
- 将该行乘以一个数,使得第一个元素(非阶梯型算起的第一个元素,是矩阵的第 c c c列)变成1;
- 用这个
1
,利用该行把该列下面的元素都消成0,注意一行的元素要一起变化。
上面循环执行完后,顺利的话,我们会得到一个上三角形:
[
1
a
12
′
⋯
a
1
n
′
b
1
′
0
1
⋯
a
2
n
′
b
2
′
⋮
⋮
⋱
⋮
⋮
0
0
⋯
1
b
n
′
]
\left[ \begin{array}{cccccc} 1 & a'_{12} & \cdots & a'_{1n} & b'_1 \\ 0 & 1 & \cdots & a'_{2n} & b'_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 & b'_n \end{array} \right]
⎣⎢⎢⎢⎡10⋮0a12′1⋮0⋯⋯⋱⋯a1n′a2n′⋮1b1′b2′⋮bn′⎦⎥⎥⎥⎤
上三角形对应的唯一解这么求:
- 根据最后一个方程得出 x n = b n ′ x_n=b'_n xn=bn′;
- 将 x n x_n xn代入第 n − 1 n-1 n−1个方程,得出 x n − 1 = b n − 1 ′ − a n − 1 , n ′ x n x_{n-1}=b'_{n-1}-a'_{n-1,n}x_n xn−1=bn−1′−an−1,n′xn;
- 不断向上回代直到第一个方程,得到 x 1 x_1 x1,解毕。
但有时候不能得到上三角形,只能得到一个类似阶梯型。比如我们现在的
n
=
4
n=4
n=4的增广矩阵消到这个情况:
[
1
3
2
4
∣
5
0
1
0
3
∣
4
0
0
0
1
∣
4
0
0
0
4
∣
3
]
\left[ \begin{array}{cccccc} 1 & 3 & 2 & 4 & | & 5 \\ 0 & 1 & 0 & 3 & | & 4 \\ 0 & 0 & 0 & 1 & | & 4 \\ 0 & 0 & 0 & 4 & | & 3 \\ \end{array} \right]
⎣⎢⎢⎡1000310020004314∣∣∣∣5443⎦⎥⎥⎤
此时我们准备消第3列。但是发现,这一列剩下的元素(
(
3
,
3
)
,
(
4
,
3
)
(3,3),(4,3)
(3,3),(4,3))中绝对值最大的元素是0
,也就是说这一列全是0
。那不管我们怎么搞,这一列都整不出来数字1
。所以对角线位置
(
3
,
3
)
(3,3)
(3,3)上不能是1
了。没办法,只能往右边走一列。现在我们走到第4列,这一列的两个元素1,4
,有非零绝对值的,所以可以正常消。消完之后得到:
[
1
3
2
4
∣
5
0
1
0
3
∣
4
0
0
0
1
∣
3
/
4
0
0
0
0
∣
15
/
4
]
\left[ \begin{array}{cccccc} 1 & 3 & 2 & 4 & | & 5 \\ 0 & 1 & 0 & 3 & | & 4 \\ 0 & 0 & 0 & 1 & | & 3/4 \\ 0 & 0 & 0 & 0 & | & 15/4 \\ \end{array} \right]
⎣⎢⎢⎡1000310020004310∣∣∣∣543/415/4⎦⎥⎥⎤
记住,只要中间有一次整列都是0
,那最后底下就会多出一行系数全是0
的方程(比如这里在消第3列的时候出现全是0的情况,所以最后一行对应出现一个系数全为0的方程)。我们要检查底下的所有系数全是0
的方程。如果
- 出现一次诸如 0 = 15 / 4 0=15/4 0=15/4,那方程直接无解。因为这个方程永远无法满足;
- 如果所有方程都是 0 = 0 0=0 0=0,那说明有无效不矛盾方程,方程有无穷解。
代码中,最复杂的是初等行列变化怎么实现,仔细看是不难的:
例题:高斯消元解线性方程组
#include <iostream>
#include <cmath>
using namespace std;
const int N = 110;
const double eps = 1e-6;
int n;
double a[N][N];
int gauss()
{
int c, r; // column, row
for (c = 0, r = 0; c < n; c ++)
{
int t = r; // 从对角线元素开始往下遍历这一列
for (int i = t; i < n; i ++) // 找到该列绝对值最大的元素所在的行号
if (fabs(a[i][c]) > fabs(a[t][c]))
t = i;
if (fabs(a[t][c]) < eps) continue; // 最大绝对值是0,那么这一列剩下的元素全是0,不用管这列。这个地方导致后面不能r ++,也意味着底部将会增加一个系数全0的方程
for (int i = c; i <= n; i ++) swap(a[t][i], a[r][i]); // 把最大绝对值元素所在的行换到未处理行的最上面(即当前要处理的的第r行)
for (int i = n; i >= c; i --) a[r][i] /= a[r][c]; // 把现在的第r行的数字全部除以一个系数,使得左上角a[r][c]变成1
for (int i = r + 1; i < n; i ++) // 把当前列下的所有数都消成0,要对应两行元素一起变化
if (fabs(a[i][c]) > eps) //已经是0的就不用操作了,省点计算
for (int j = n; j >= c; j --)
a[i][j] -= a[r][j] * a[i][c];
r ++;
}
// 上面步骤走完之后,矩阵a[][]扣掉增广的最后一列系数以外,剩下的已经是个上三角阵或者阶梯阵
if (r < n) // 说明有效的方程个数小于n,那要么无穷解,要么无解
{
for (int i = r; i < n; i ++){
if (fabs(a[i][n]) > eps) // a[i][n] = b_i不等于0
return 2; // 无解
return 1; // 都是0 = 0的方程,无穷解
}
}
// 唯一解,从下往上回代,得到方程的解
for (int i = n - 1; i >= 0; i --)
for (int j = i + 1; j < n; j ++)
a[i][n] -= a[i][j] * a[j][n];
return 0;
}
int main(){
cin >> n;
for (int i = 0; i < n; i ++)
for (int j = 0; j <= n; j ++)
cin >> a[i][j];
int t = gauss();
if (t == 0)
for (int i = 0; i < n; i ++) printf("%.2f\n", a[i][n]);
else if (t == 1) puts("Infinite group solutions");
else puts("No solution");
return 0;
}