列主元Gauss消去法
3.1 题目
对于某电路的分析,归结为就求解线性方程组
R
I
=
V
\pmb{RI=V}
RI=V,其中
R
=
[
31
−
13
0
0
0
−
10
0
0
0
−
13
35
−
9
0
−
11
0
0
0
0
0
−
9
31
−
10
0
0
0
0
0
0
0
−
10
79
−
30
0
0
0
−
9
0
0
0
−
30
57
−
7
0
−
5
0
0
0
0
0
−
7
47
−
30
0
0
0
0
0
0
0
−
30
41
0
0
0
0
0
0
−
5
0
0
27
−
2
0
0
0
−
9
0
0
0
−
2
29
]
\pmb{R} = \begin{bmatrix} 31 & -13 & 0 & 0 & 0 & -10 & 0 & 0 & 0\\ -13 & 35 & -9 & 0 & -11 & 0 & 0 & 0 & 0 \\ 0 & -9 & 31 & -10 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -10 & 79 & -30 & 0 & 0 & 0 & -9 \\ 0 & 0 & 0 & -30 & 57 & -7 & 0 & -5 & 0 \\ 0 & 0 & 0 & 0 & -7 & 47 & -30 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & -30 & 41 & 0 & 0 \\ 0 & 0 & 0 & 0 & -5 & 0 & 0 & 27 & -2& \\ 0 & 0 & 0 & -9 & 0 & 0 & 0 & -2 & 29 \end{bmatrix}
R=
31−130000000−1335−90000000−931−100000000−1079−30000−90−110−3057−70−50−10000−747−300000000−3041000000−50027−2000−9000−229
V
T
=
(
−
15
,
27
,
−
23
,
0
−
20
,
12
,
−
7
,
7
,
10
)
\pmb{V}^T=(-15,27,-23,0-20,12,-7,7,10)
VT=(−15,27,−23,0−20,12,−7,7,10)
(1) 编写解
n
n
n阶线性方程组
A
x
=
b
\pmb{Ax=b}
Ax=b的列主元Gauss消去法的通用程序;
(2) 用所编程序解线性方程组
R
I
=
V
\pmb{RI=V}
RI=V ,并打印出解向量,保留5位有效数字;
(3) 通过本题程序的编写,你提高了哪些编程能力?
3.2 Python源程序
线性方程组 A x = b \pmb{Ax=b} Ax=b的列主元Gauss消去法
求解不同的线性方程组 A x = b Ax=b Ax=b,只需要修改参数 A \pmb{A} A与 b \pmb{b} b
import numpy as np
A = np.array([[31, -13, 0, 0, 0, -10, 0, 0, 0],
[-13, 35, -9, 0, -11, 0, 0, 0, 0],
[0, -9, 31, -10, 0, 0, 0, 0, 0],
[0, 0, -10, 79, -30, 0, 0, 0, -9],
[0, 0, 0, -30, 57, -7, 0, -5, 0],
[0, 0, 0, 0, -7, 47, -30, 0, 0],
[0, 0, 0, 0, 0, -30, 41, 0, 0],
[0, 0, 0, 0, -5, 0, 0, 27, -2],
[0, 0, 0, -9, 0, 0, 0, -2, 29]], dtype=np.float64)
b = np.array([[-15], [27], [-23], [0], [-20], [12], [-7], [7], [10]], dtype=np.float64)
row_num = A.shape[0] # 系数矩阵A的行数
col_num = A.shape[1] # 系数矩阵A的列数
x = np.zeros(row_num)
# A_表示为增广矩阵
A_ = A.copy()
A_abs = np.absolute(A_)
A_ = np.append(A_, values=b, axis=1)
for i in range(row_num-1):
# 处理主元为0的情况
if A_[i,i] == 0:
m = i + 1
while A_[m, i] == 1:
m += 1
A_[[i, m], :] = A_[[m, i], :]
# 选取列主元
index = A_abs[range(i, row_num), i].argmax()
index += i
A_[[i, index], :] = A_[[index, i], :]
# 消元
for j in range(i, row_num-1):
k = A_[j+1, i]/A_[i, i]
A_[j+1, :] = A_[j+1, :] - k * A_[i, :]
# # 回代求解向量
x[row_num-1] = A_[row_num-1, row_num]/ A_[row_num-1, row_num-1]
for i in range(row_num - 2, -1, -1):
x[i] = (A_[i, row_num] - (A_[i, i+1:row_num].dot(x[i+1:row_num])))/A_[i, i]
print("解向量x为")
for i in range(row_num):
print("{:.5f}".format(x[i]))
3.3 程序运行结果
解向量x为
-0.28923381601575443
0.34543571577911536
-0.7128117310868789
-0.22060851057052858
-0.43040043270402234
0.15430873983831117
-0.05782287328904061
0.2010538948236807
0.29022866187974494
根据题意,保留5位有效数字,则方程
R
I
=
V
\pmb{RI=V}
RI=V的解向量
I
\pmb{I}
I为
I
=
[
−
0.28923
0.34544
−
0.71281
−
0.22061
−
0.43040
0.15431
−
0.057823
0.20105
0.29023
]
\pmb{I} = \begin{bmatrix} -0.28923\\0.34544\\-0.71281\\-0.22061\\-0.43040\\0.15431\\-0.057823\\0.20105\\0.29023 \end{bmatrix}
I=
−0.289230.34544−0.71281−0.22061−0.430400.15431−0.0578230.201050.29023
3.4 总结感悟
通过本次线性方程组的列主元Gauss消元法编程,熟悉掌握了Python的重要工具Numpy模块的用法,对复杂数据的编程操作的能力进一步提升,对列主元高斯消去法也有了进一步的认识,我们通常可以通过笔算很快地得到一些简单线性方程组的解,但是真正将这一算法通过严谨的程序编写出来很是有很大的难度。列主元 Gauss 消去法可以减小由消去过程中分母相较分子过小(两个相差很大的数相除)引起的误差,列主元Gauss消去法基本上是稳定的。