Gauss-Seidel method 解释与实现

Gauss-Seidel method

算法定义

定义线性方程组：

$$\{\begin{array}{c}{a}_{11}{x}_1+a_{12}{x}_2+a_{13}{x}_3+a_{14}{x}_4+\cdots +a_{1i}{x}_i=b_1\\ a_{21}{x}_1+a_{22}{x}_2+a_{23}{x}_3+a_{24}{x}_4+\cdots +a_{2i}{x}_i=b_2\\ a_{31}{x}_1+a_{32}{x}_2+a_{33}{x}_3+a_{34}{x}_4+\cdots +a_{3i}{x}_i=b_3\\ a_{41}{x}_1+a_{42}{x}_2+a_{43}{x}_3+a_{44}{x}_4+\cdots +a_{4i}{x}_i=b_4\\ ⋮\\ a_{i1}{x}_1+a_{i2}{x}_2+a_{i3}{x}_3+a_{i4}{x}_4+\cdots +a_{ii}{x}_i=b_i\end{array}$$

经过变换可以得到：

$$\{\begin{array}{c}{x}_1=-\frac{1}{a_{11}}\left(\bcancel{a_{11}{x}_1}+a_{12}{x}_2+a_{13}{x}_3+a_{14}{x}_4+\cdots +a_{1i}{x}_i-b_1\right)\\ \\ x_2=-\frac{1}{a_{22}}\left(a_{21}{x}_1+\bcancel{a_{22}{x}_2}+a_{23}{x}_3+a_{24}{x}_4+\cdots +a_{2i}{x}_i-b_2\right)\\ \\ x_3=-\frac{1}{a_{33}}\left(a_{31}{x}_1+a_{32}{x}_2+\bcancel{a_{33}{x}_3}+a_{34}{x}_4+\cdots +a_{3i}{x}_i-b_3\right)\\ \\ x_4=-\frac{1}{a_{44}}\left(a_{41}{x}_1+a_{42}{x}_2+a_{43}{x}_3+\bcancel{a_{44}{x}_4}+\cdots +a_{4i}{x}_i-b_4\right)\\ ⋮\\ x_i=-\frac{1}{a_{ii}}\left(a_{i1}{x}_1+a_{i2}{x}_2+a_{i3}{x}_3+a_{i4}{x}_4+\cdots +\bcancel{a_{ii}{x}_i}-b_i\right)\end{array}$$

所以迭代的过程可以这样表示：

$$\{\begin{array}{l}{x}_1^{(k+1)}=-\frac{1}{a_{11}}\left(\bcancel{a_{11}{x}_1}+a_{12}{x}_2^{(k)}+a_{13}{x}_3^{(k)}+a_{14}{x}_4^{(k)}+\cdots +a_{1i}{x}_i^{(k)}-b_1\right)\\ \\ x_2^{(k+1)}=-\frac{1}{a_{22}}\left(a_{21}{x}_1^{(k+1)}+\bcancel{a_{22}{x}_2}+a_{23}{x}_3^{(k)}+a_{24}{x}_4^{(k)}+\cdots +a_{2i}{x}_i^{(k)}-b_2\right)\\ \\ x_3^{(k+1)}=-\frac{1}{a_{33}}\left(a_{31}{x}_1^{(k+1)}+a_{32}{x}_2^{(k+1)}+\bcancel{a_{33}{x}_3}+a_{34}{x}_4^{(k)}+\cdots +a_{3i}{x}_i^{(k)}-b_3\right)\\ \\ x_4^{(k+1)}=-\frac{1}{a_{44}}\left(a_{41}{x}_1^{(k+1)}+a_{42}{x}_2^{(k+1)}+a_{43}{x}_3^{(k+1)}+\bcancel{a_{44}{x}_4}+\cdots +a_{4i}{x}_i^{(k)}-b_4\right)\\ ⋮\\ x_i^{(k+1)}=-\frac{1}{a_{ii}}\left(a_{i1}{x}_1^{(k+1)}+a_{i2}{x}_2^{(k+1)}+a_{i3}{x}_3^{(k+1)}+a_{i4}{x}_4^{(k+1)}+\cdots +\bcancel{a_{ii}{x}_i}-b_i\right)\end{array}$$

用矩阵的方式表示：

$$\begin{array}{rl}\mathbf{AX}& =\mathbf{b}\\ \mathbf{LX}+\mathbf{DX}+\mathbf{UX}& =\mathbf{b}\\ {\mathbf{X}}^{(\mathbf{k}+1)}& ={\mathbf{D}}^{-1}\left(\mathbf{b}-\mathbf{L}{\mathbf{X}}^{(\mathbf{k}+1)}-\mathbf{U}{\mathbf{X}}^{(\mathbf{k})}\right)\\ {\mathbf{X}}^{(\mathbf{k}+1)}& ={\left(\mathbf{D}+\mathbf{L}\right)}^{-1}\mathbf{b}-{\left(\mathbf{D}+\mathbf{L}\right)}^{-1}\mathbf{U}{\mathbf{X}}^{(\mathbf{k})}\\ {\mathbf{X}}^{(\mathbf{k}+1)}& =f_G-B_G{\mathbf{X}}^{(\mathbf{k})}\end{array}$$

Python代码实现

采用 Python 的 Numpy 库实现矩阵运算

import numpy as np


# 采用标准式 AX= B // 及使用增广矩阵的形式
def gauss_seidel(a_list: list, init_x_list: list, eps: int) -> (int, list):
    """
    此函数采用增广矩阵的形式表示线性方程 及 AX=B
    :param eps: 限制精度范围
    :param a_list: 输入的增广矩阵
    :param init_x_list: 初始假设的X向量
    :return: 返回收敛数值 在控制台输出
    """
    # ---------------- 获取矩阵
    # 创建矩阵
    a_mat = np.mat(a_list)
    # 系数矩阵
    coe_mat = a_mat[:, 0:-1]
    # 增广列
    b_mat = a_mat[:, -1]
    # 对角矩阵
    d_mat = np.mat(np.diag(np.diag(coe_mat)))
    # 下三角矩阵
    u_mat = np.mat(np.triu(coe_mat, 1))
    # 上三角矩阵
    l_mat = np.mat(np.tril(coe_mat, -1))

    # ---------------- 矩阵计算
    # 计算 (D - L)^-1
    d_sub_l_mat_inv = np.linalg.inv(d_mat + l_mat)
    # 计算 B_G 和 f_G
    b_g = np.matrix.dot(d_sub_l_mat_inv, u_mat)
    f_g = np.matrix.dot(d_sub_l_mat_inv, b_mat)
    # 定义初始值
    init_x_mat = np.mat(init_x_list)
    init_x_mat = init_x_mat.T
    print("The init result")
    for a_cell in init_x_mat.tolist():
        print(a_cell[0], end="  ")

    # ---------------- 开始迭代
    count = 0
    while count < 100000:  # 如果10000 次不收敛则结果发散
        result_x = (f_g - np.matrix.dot(b_g, init_x_mat)).round(eps)
        if np.all(init_x_mat == result_x):
            break
        else:
            print()
            print("The", count + 1, "result:")
            for a_cell in result_x.tolist():
                print(a_cell[0], end="  ")
            init_x_mat = result_x
            count += 1
    print()
    return count, init_x_mat.T.tolist()[0]

# 这里是测试， 
if __name__ == '__main__':
    c = [[-2.045, 1, 0, -100.9],
         [1, -2.045, 1, -0.9],
         [0, 1, -1.0225, -0.45]]
    re_c = gauss_seidel(c, [0, 0, 0], 2)
    print(re_c)

注释

• 为了方便编写程序，这里对原本的算法定义做了一定的更改
  高斯-赛德尔迭代算法介绍（重庆交通大学ppt）
• 为了更好的贴合定义，这里仍然采用了$\mathbf{L},\mathbf{D},\mathbf{U}$的定义方法，并且在写程序时也没有做简化，其实$\mathbf{L}+\mathbf{D}$ 可以合并来看就是下三角矩阵

实验测试

给出线性方程组：

$$\{\begin{array}{l}{t}_3-2.045t_2+100.9=0\\ t_2-2.045t_3+t4+0.9=0\\ t_3-1.0225t_4+0.45=0\end{array}$$

输出的结果：

$$t_2=92.2℃,t_3=87.7℃,t_4=86.2℃$$

The init result
0  0  0  
The 1 result:
49.34  24.57  24.47  
The 2 result:
61.35  42.41  41.91  
The 3 result:
70.08  55.2  54.43  
The 4 result:
76.33  64.38  63.41  
The 5 result:
80.82  70.97  69.85  
The 6 result:
84.04  75.69  74.47  
The 7 result:
86.35  79.08  77.78  
The 8 result:
88.01  81.51  80.16  
The 9 result:
89.2  83.26  81.86  
The 10 result:
90.05  84.51  83.09  
The 11 result:
90.67  85.41  83.97  
The 12 result:
91.11  86.05  84.6  
The 13 result:
91.42  86.51  85.05  
The 14 result:
91.64  86.84  85.37  
The 15 result:
91.8  87.08  85.6  
The 16 result:
91.92  87.25  85.77  
The 17 result:
92.0  87.37  85.89  
The 18 result:
92.06  87.46  85.97  
The 19 result:
92.11  87.52  86.03  
The 20 result:
92.14  87.56  86.08  
The 21 result:
92.16  87.6  86.11  
The 22 result:
92.18  87.62  86.13  
The 23 result:
92.19  87.64  86.15  
The 24 result:
92.2  87.65  86.16  
The 25 result:
92.2  87.66  86.17  
The 26 result:
92.21  87.67  86.18  
(26, [92.21, 87.67, 86.18])

Process finished with exit code 0

需要改进

• 需要对代码进行封装
• 需要列出伪代码式
• 需要扩展到其他语言

结束语

       第一次写 CSDN 的文档，属于是闲来无事记录一下。免得时间久了，便渐渐忘却了。如果有幸可以帮助你解决问题，那么我们也是有缘人了。如果以后有时间的话会多在这个平台上发表一些理解和看法。我现在还是一个计算机小白，许多问题还需要不断去探索，也希望大家能够为我指出错误。