Python——高斯赛德尔迭代求线性方程组的根

高斯赛德尔迭代

与雅克比迭代对比：Python——雅克比迭代求线性方程组的根

计算代码

矩阵A——方程组的系数。代码中只需改变矩阵系数即可
矩阵B——方程组等号右侧的常数。

#高斯赛德尔迭代
#求解方程组：
#10*x1 - 2*x2 - 1*x3 = 3
#-2*x1 + 10*x2 - 1*x3 = 15
#-1*x1 - 2*x2 + 5*x3 = 10
#AX = B

#只用改变矩阵A、B的系数即可
A = [[10, -2, -1], [-2, 10, -1], [-1, -2, 5]]
B = [[3], [15], [10]]
#初始值
x = [0, 0, 0]
aij = None
n = 0

while True:
    v_max = 0
    n += 1
    print(f'第{n}次迭代')
    for i in range(len(A)):
        a = A[i]
        ax = 0
        for j in range(len(a)):
            if j == i :
                aij = a[j]
            else:
                ax += a[j]*x[j]
        x_new = (B[i][0] - ax) / aij
        v = abs(x[i] - x_new)
        v_max = max(v_max, v)
        x[i] = x_new

    if v_max <= 10**-3:
        print(x)
        print('e=',v_max)
        break
    else:
        print(x)
        print('e=',v_max)

结果：

第1次迭代
[0.3, 1.56, 2.684]
e= 2.684
第2次迭代
[0.8804000000000001, 1.94448, 2.9538719999999996]
e= 0.5804
第3次迭代
[0.9842831999999999, 1.9922438399999998, 2.993754176]
e= 0.10388319999999984
第4次迭代
[0.9978241856000001, 1.9989402547200001, 2.999140939008]
e= 0.013540985600000166
第5次迭代
[0.9997021448448, 1.9998545228697602, 2.999882238116864]
e= 0.0018779592447999427
第6次迭代
[0.9999591283856384, 1.9999800494888142, 2.9999838454726535]
e= 0.0002569835408383625

高斯赛德尔代收敛的充分必要条件

L——矩阵A的下三角矩阵
U——矩阵A的上三角矩阵
D——矩阵A的对角矩阵
例如：

from numpy import *
# numpy.linalg模块包含线性代数的函数。使用这个模块，可以计算逆矩阵、求特征值、解线性方程组以及求解行列式等
from numpy.linalg import *
A = mat([[25, 5, 1], [64, 8, 1], [144, 12, 1]])
B = mat([[106.8], [177.2], [279.2]])
L = tril(A, -1) #下三角矩阵
U = triu(A,1) #上三角矩阵
D = diag(diag(A))  #对家矩阵

G = dot(inv(D+L),-U)
print(eigvals(G))  #特征值

结果：

[0.        0.8668019 4.1531981]

特征值中只要有一个绝对值＞1，用高斯赛德尔迭代求解就不收敛
所以上面的方程组用高斯赛德尔迭代不收敛
再例如：

from numpy import *
# numpy.linalg模块包含线性代数的函数。使用这个模块，可以计算逆矩阵、求特征值、解线性方程组以及求解行列式等
from numpy.linalg import *
A = mat([[2, 1, 1], [1, 2, 1], [1, 1, 2]])
B = mat([[4], [2], [0]])
L = tril(A, -1) #下三角矩阵
U = triu(A,1) #上三角矩阵
D = diag(diag(A))  #对角矩阵

G = dot(inv(D+L),-U)
print(eigvals(G))  #特征值

结果：

[0.    +0.j         0.3125+0.16535946j 0.3125-0.16535946j]

三个特征值的绝对值都<1，所以用高斯赛德尔迭代可以求解这个方程组。

先判断是否收敛再迭代的综合python代码

#高斯赛德尔迭代
#求解方程组：
#10*x1 - 2*x2 - 1*x3 = 3
#-2*x1 + 10*x2 - 1*x3 = 15
#-1*x1 - 2*x2 + 5*x3 = 10
#AX = B
from numpy import *
# numpy.linalg模块包含线性代数的函数。使用这个模块，可以计算逆矩阵、求特征值、解线性方程组以及求解行列式等
from numpy.linalg import *
#只用改变矩阵A、B的系数即可
A = [[10, -2, -1], [-2, 10, -1], [-1, -2, 5]]
B = [[3], [15], [10]]

#判断该矩阵用高斯-赛德尔迭代是否收敛
L = tril(A, -1) #下三角矩阵
U = triu(A,1) #上三角矩阵
D = diag(diag(A))  #对角矩阵

G = dot(inv(D+L),-U)
root = eigvals(G)
flag = True
for r in root:
    if r > 1:
        print('高斯-赛德尔迭代不收敛')
        flag = False
        break
    else:
        pass

#如果迭代收敛，求迭代
if flag:
    # 初始值
    x = [0, 0, 0]
    aij = None
    n = 0
    while True:
        v_max = 0
        n += 1
        print(f'第{n}次迭代')
        for i in range(len(A)):
            a = A[i]
            ax = 0
            for j in range(len(a)):
                if j == i:
                    aij = a[j]
                else:
                    ax += a[j] * x[j]
            x_new = (B[i][0] - ax) / aij
            v = abs(x[i] - x_new)
            v_max = max(v_max, v)
            x[i] = x_new

        if v_max <= 10 ** -3:
            print(x)
            print('e=', v_max)
            break
        else:
            print(x)
            print('e=', v_max)

结果：

第1次迭代
[0.3, 1.56, 2.684]
e= 2.684
第2次迭代
[0.8804000000000001, 1.94448, 2.9538719999999996]
e= 0.5804
第3次迭代
[0.9842831999999999, 1.9922438399999998, 2.993754176]
e= 0.10388319999999984
第4次迭代
[0.9978241856000001, 1.9989402547200001, 2.999140939008]
e= 0.013540985600000166
第5次迭代
[0.9997021448448, 1.9998545228697602, 2.999882238116864]
e= 0.0018779592447999427
第6次迭代
[0.9999591283856384, 1.9999800494888142, 2.9999838454726535]
e= 0.0002569835408383625