列主元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消去法基本上是稳定的。