Gauss-Seidel method
算法定义
定义线性方程组:
{ a 11 x 1 + a 12 x 2 + a 13 x 3 + a 14 x 4 + ⋯ + a 1 i x i = b 1 a 21 x 1 + a 22 x 2 + a 23 x 3 + a 24 x 4 + ⋯ + a 2 i x i = b 2 a 31 x 1 + a 32 x 2 + a 33 x 3 + a 34 x 4 + ⋯ + a 3 i x i = b 3 a 41 x 1 + a 42 x 2 + a 43 x 3 + a 44 x 4 + ⋯ + a 4 i x i = b 4 ⋮ a i 1 x 1 + a i 2 x 2 + a i 3 x 3 + a i 4 x 4 + ⋯ + a i i x i = b i \left\{ \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\\ \vdots\\ 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} \right. ⎩ ⎨ ⎧a11x1+a12x2+a13x3+a14x4+⋯+a1ixi=b1a21x1+a22x2+a23x3+a24x4+⋯+a2ixi=b2a31x1+a32x2+a33x3+a34x4+⋯+a3ixi=b3a41x1+a42x2+a43x3+a44x4+⋯+a4ixi=b4⋮ai1x1+ai2x2+ai3x3+ai4x4+⋯+aiixi=bi
经过变换可以得到:
{ x 1 = − 1 a 11 ( a 11 x 1 + a 12 x 2 + a 13 x 3 + a 14 x 4 + ⋯ + a 1 i x i − b 1 ) x 2 = − 1 a 22 ( a 21 x 1 + a 22 x 2 + a 23 x 3 + a 24 x 4 + ⋯ + a 2 i x i − b 2 ) x 3 = − 1 a 33 ( a 31 x 1 + a 32 x 2 + a 33 x 3 + a 34 x 4 + ⋯ + a 3 i x i − b 3 ) x 4 = − 1 a 44 ( a 41 x 1 + a 42 x 2 + a 43 x 3 + a 44 x 4 + ⋯ + a 4 i x i − b 4 ) ⋮ x i = − 1 a i i ( a i 1 x 1 + a i 2 x 2 + a i 3 x 3 + a i 4 x 4 + ⋯ + a i i x i − b i ) \left\{ \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)\\ \vdots\\ 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} \right. ⎩ ⎨ ⎧x1=−a111(a11x1 +a12x2+a13x3+a14x4+⋯+a1ixi−b1)x2=−a221(a21x1+a22x2 +a23x3+a24x4+⋯+a2ixi−b2)x3=−a331(a31x1+a32x2+a33x3 +a34x4+⋯+a3ixi−b3)x4=−a441(a41x1+a42x2+a43x3+a44x4 +⋯+a4ixi−b4)⋮xi=−aii1(ai1x1+ai2x2+ai3x3+ai4x4+⋯+aiixi −bi)
所以迭代的过程可以这样表示:
{ x 1 ( k + 1 ) = − 1 a 11 ( a 11 x 1 + a 12 x 2 ( k ) + a 13 x 3 ( k ) + a 14 x 4 ( k ) + ⋯ + a 1 i x i ( k ) − b 1 ) x 2 ( k + 1 ) = − 1 a 22 ( a 21 x 1 ( k + 1 ) + a 22 x 2 + a 23 x 3 ( k ) + a 24 x 4 ( k ) + ⋯ + a 2 i x i ( k ) − b 2 ) x 3 ( k + 1 ) = − 1 a 33 ( a 31 x 1 ( k + 1 ) + a 32 x 2 ( k + 1 ) + a 33 x 3 + a 34 x 4 ( k ) + ⋯ + a 3 i x i ( k ) − b 3 ) x 4 ( k + 1 ) = − 1 a 44 ( a 41 x 1 ( k + 1 ) + a 42 x 2 ( k + 1 ) + a 43 x 3 ( k + 1 ) + a 44 x 4 + ⋯ + a 4 i x i ( k ) − b 4 ) ⋮ x i ( k + 1 ) = − 1 a i i ( a i 1 x 1 ( k + 1 ) + a i 2 x 2 ( k + 1 ) + a i 3 x 3 ( k + 1 ) + a i 4 x 4 ( k + 1 ) + ⋯ + a i i x i − b i ) \begin{cases} 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)\\ \vdots\\ 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{cases} ⎩ ⎨ ⎧x1(k+1)=−a111(a11x1 +a12x2(k)+a13x3(k)+a14x4(k)+⋯+a1ixi(k)−b1)x2(k+1)=−a221(a21x1(k+1)+a22x2 +a23x3(k)+a24x4(k)+⋯+a2ixi(k)−b2)x3(k+1)=−a331(a31x1(k+1)+a32x2(k+1)+a33x3 +a34x4(k)+⋯+a3ixi(k)−b3)x4(k+1)=−a441(a41x1(k+1)+a42x2(k+1)+a43x3(k+1)+a44x4 +⋯+a4ixi(k)−b4)⋮xi(k+1)=−aii1(ai1x1(k+1)+ai2x2(k+1)+ai3x3(k+1)+ai4x4(k+1)+⋯+aiixi −bi)
用矩阵的方式表示:
A X = b L X + D X + U X = b X ( k + 1 ) = D − 1 ( b − L X ( k + 1 ) − U X ( k ) ) X ( k + 1 ) = ( D + L ) − 1 b − ( D + L ) − 1 U X ( k ) X ( k + 1 ) = f G − B G X ( k ) \begin{align} \bf AX&= \bf b\\ \bf LX+DX+UX&=\bf b\\ \bf X^{(k+1)}& = \bf D^{-1}\left( b- LX^{(k+1)}-UX^{(k)} \right)\\ \bf X^{(k+1)}& = \bf \left(D+L\right)^{-1}b - \left(D+L\right)^{-1}UX^{(k)}\\ \bf X^{(k+1)}& = f_G - B_G\bf X^{(k)} \end{align} AXLX+DX+UXX(k+1)X(k+1)X(k+1)=b=b=D−1(b−LX(k+1)−UX(k))=(D+L)−1b−(D+L)−1UX(k)=fG−BGX(k)
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) - 为了更好的贴合定义,这里仍然采用了 L , D , U \bf L, D,U L,D,U的定义方法,并且在写程序时也没有做简化,其实 L + D \bf L+ D L+D 可以合并来看就是下三角矩阵
实验测试
给出线性方程组:
{ t 3 − 2.045 t 2 + 100.9 = 0 t 2 − 2.045 t 3 + t 4 + 0.9 = 0 t 3 − 1.0225 t 4 + 0.45 = 0 \begin{cases} 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{cases} ⎩ ⎨ ⎧t3−2.045t2+100.9=0t2−2.045t3+t4+0.9=0t3−1.0225t4+0.45=0
输出的结果:
t 2 = 92.2 ℃ , t 3 = 87.7 ℃ , t 4 = 86.2 ℃ t_2 = 92.2 ℃,\space t_3=87.7℃,\space t_4=86.2℃ t2=92.2℃, t3=87.7℃, t4=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 的文档,属于是闲来无事记录一下。免得时间久了,便渐渐忘却了。如果有幸可以帮助你解决问题,那么我们也是有缘人了。如果以后有时间的话会多在这个平台上发表一些理解和看法。我现在还是一个计算机小白,许多问题还需要不断去探索,也希望大家能够为我指出错误。