东南大学研究生-数值分析上机题（2023）Python 3 线性代数方程组数值解法

列主元Gauss消去法

3.1 题目

对于某电路的分析，归结为就求解线性方程组$RI=V$，其中
$$R=\left[\begin{array}{cccccccccc}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{array}\right]$$
$$V^T=(-15,27,-23,0-20,12,-7,7,10)$$
(1) 编写解$n$阶线性方程组$Ax=b$的列主元Gauss消去法的通用程序；
(2) 用所编程序解线性方程组$RI=V$ ，并打印出解向量，保留5位有效数字；
(3) 通过本题程序的编写，你提高了哪些编程能力？

3.2 Python源程序

线性方程组$Ax=b$的列主元Gauss消去法

求解不同的线性方程组$Ax=b$，只需要修改参数$A$与$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位有效数字，则方程$RI=V$的解向量$I$为
$$I=\left[\begin{array}{c}-0.28923\\ 0.34544\\ -0.71281\\ -0.22061\\ -0.43040\\ 0.15431\\ -0.057823\\ 0.20105\\ 0.29023\end{array}\right]$$

3.4 总结感悟

通过本次线性方程组的列主元Gauss消元法编程，熟悉掌握了Python的重要工具Numpy模块的用法，对复杂数据的编程操作的能力进一步提升，对列主元高斯消去法也有了进一步的认识，我们通常可以通过笔算很快地得到一些简单线性方程组的解，但是真正将这一算法通过严谨的程序编写出来很是有很大的难度。列主元 Gauss 消去法可以减小由消去过程中分母相较分子过小（两个相差很大的数相除）引起的误差，列主元Gauss消去法基本上是稳定的。