高斯消元是一种利用矩阵的算法，常用来实现多组线性方程组求解，也可求出矩阵的秩或是逆矩阵。本文主要说明高斯消元是如何对线性方程组求解的。
假设有如下三个方程组：
2x+y−z=8
−3x−y+2z=−11
−2x+y−2=−3
利用高斯消元法，首先将三组方程组得到扩充后的增广矩阵
⎣⎢⎡2−3−21−11−12−2∣∣∣8−11−3⎦⎥⎤
①先将行根据第一列 (是第一个未知数) 的绝对值最大行置顶，代码实现：
1
2
3
4
5
6
| for i in range(n):
r = i
for j in range(i,n):
if abs(mat[j][i]) > abs(mat[r][i]):
r = j
mat[r],mat[i] = mat[i],mat[r]
|
②将除第一行外，其他行的第一个未知数的系数消为 0，接着是第二行、第三行直到得到一个倒三角，代码实现如下，例题中即得到
⎣⎢⎡−300−1−502−23∣∣∣−11−13−3⎦⎥⎤
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
| for i in range(n):
r = i
for j in range(i,n):
if abs(mat[j][i]) > abs(mat[r][i]):
r = j
mat[r],mat[i] = mat[i],mat[r]
for j in range(i + 1,n):
tmp = int(gmpy2.lcm(mat[j][i],mat[i][i]))
n1,n2 = tmp // mat[j][i],tmp // mat[i][i]
for k in range(n + 1):
mat[j][k] = mat[j][k] * n1 - mat[i][k] * n2
print()
matrix_print(mat)
|
③之后再循环回代，根据以求得的倒三角矩阵逆序求得每一项的值，即得到矩阵
⎣⎢⎡100010001∣∣∣23−1⎦⎥⎤
该方程组的解为x=2,y=3,z=1
循环回代的代码实现为：
1
2
3
4
5
6
7
| for i in range(n - 1,-1,-1):
tmp = mat[i][i]
for j in range(i + 1,n):
mat[i][n] -= mat[i][j] * mat[j][n]
mat[i][j] = 0
mat[i][n] //= mat[i][i]
mat[i][i] = 1
|
最后贴两个完整版脚本，一个是精度略低的，一个是用公倍数求能有更高的精度，代码实现如下
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
|
def matrix_print(mat):
for i in mat:
print(i)
def select(mat,n):
for i in range(n):
r = i
for j in range(i,n):
if abs(mat[j][i]) > abs(mat[r][i]):
r = j
mat[r],mat[i] = mat[i],mat[r]
for j in range(i + 1,n):
tmp = mat[j][i] / mat[i][i]
for k in range(n + 1):
mat[j][k] -= tmp * mat[i][k]
print()
matrix_print(mat)
return
def Gauss(mat,n):
ans = []
select(mat,n)
for i in range(n - 1,-1,-1):
tmp = 1 / mat[i][i]
for j in range(i,n + 1):
mat[i][j] *= tmp
if i != n - 1:
for j in range(i + 1,n):
mat[i][n] -= mat[i][j] * mat[j][n]
mat[i][j] = 0
print()
matrix_print(mat)
return
mat = [
[2,1,-1,8],
[-3,-1,2,-11],
[-2,1,2,-3]
]
n = len(mat)
matrix_print(mat)
Gauss(mat,n)
|
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
51
52
53
54
55
56
57
| import gmpy2
def matrix_print(mat):
for i in mat:
print(i)
def select(mat,n):
for i in range(n):
r = i
for j in range(i,n):
if abs(mat[j][i]) > abs(mat[r][i]):
r = j
mat[r],mat[i] = mat[i],mat[r]
for j in range(i + 1,n):
tmp = int(gmpy2.lcm(mat[j][i],mat[i][i]))
n1,n2 = tmp // mat[j][i],tmp // mat[i][i]
for k in range(n + 1):
mat[j][k] = mat[j][k] * n1 - mat[i][k] * n2
print()
matrix_print(mat)
return
def Gauss(mat,n):
select(mat,n)
for i in range(n - 1,-1,-1):
tmp = mat[i][i]
for j in range(i + 1,n):
mat[i][n] -= mat[i][j] * mat[j][n]
mat[i][j] = 0
mat[i][n] //= mat[i][i]
mat[i][i] = 1
print()
matrix_print(mat)
return
mat = [
[2,1,-1,8],
[-3,-1,2,-11],
[-2,1,2,-3]
]
n = len(mat)
matrix_print(mat)
Gauss(mat,n)
|
