数值分析——高斯赛德尔（Gauss_Seidel）迭代法（Python实现）

import numpy as np
import math
import sys
#分解矩阵
def DLU(A):
    D=np.zeros(np.shape(A))
    L=np.zeros(np.shape(A))
    U=np.zeros(np.shape(A))
    for i in range(A.shape[0]):
        D[i,i]=A[i,i]
        for j in range(i):
            L[i,j]=-A[i,j]
        for k in list(range(i+1,A.shape[1])):
            U[i,k]=-A[i,k]
    L=np.mat(L)
    D=np.mat(D)
    U=np.mat(U)
    return D,L,U

#迭代
def Jacobi_iterative(A,b,x0,maxN,p):  #x0为初始值，maxN为最大迭代次数，p为允许误差
    D,L,U=DLU(A)
    if len(A)==len(b):
        M=np.linalg.inv(D-L)
        D_inv=np.mat(M)
        B=D_inv * U
        B=np.mat(B)
        f=M*b
        f=np.mat(f)
    else:
        print('维数不一致')
        sys.exit(0)  # 强制退出
    
    a,b=np.linalg.eig(B) #a为特征值集合，b为特征值向量
    c=np.max(np.abs(a)) #返回谱半径
    if c<1:
        print('迭代收敛')
    else:
        print('迭代不收敛')
        sys.exit(0)  # 强制退出
#以上都是判断迭代能否进行的过程，下面正式迭代
    k=0
    while k<maxN:
        x=B*x0+f
        k=k+1
        eps=np.linalg.norm(x-x0,ord=2)
        if eps<p:
            break
        else:
            x0=x
    return k,x

A = np.array([[8,-3,2],[4,11,-1],[5,3,12]])
b = np.array([[20],[33],[36]])
# A = np.mat([[10,3,1],[2,-10,3],[1,3,10]])
# b = np.mat([[14],[-5],[14]])
x0=np.mat([[0],[0],[0]])
maxN=100
p=0.0000000001
print("原系数矩阵a:")
print(A, "\n")
print("原值矩阵b:")
print(b, "\n")
k,x=Jacobi_iterative(A,b,x0,maxN,p)
print("迭代次数")
print(k, "\n")
print("数值解")
print(x, "\n")

输出

系数矩阵a:
[[ 8 -3  2]
 [ 4 11 -1]
 [ 5  3 12]]

原值矩阵b:
[[20]
 [33]
 [36]]

迭代收敛
迭代次数
13

数值解
[[2.95056375]
 [2.04163053]
 [1.26019081]]

精确解为

x=[3,2,1]

Gauss_Seidel迭代法和Jacobi迭代法是很相似的，它们的区别只在于划分D,L,U时，前者取逆的矩阵时D-L，而后者是D。这个细小的区别使得高斯赛德尔迭代法的迭代次数要比雅可比迭代法少的多。以下题为例

A = np.mat([[10,3,1],[2,-10,3],[1,3,10]])
b = np.mat([[14],[-5],[14]])
x0=np.mat([[0],[0],[0]])
maxN=100
p=0.0000000001

Gauss_Seidel迭代法的迭代次数是15，而Jacobi迭代法的迭代次数是22.