jacobi迭代-python

相关公式：（老师的ppt）

代码：

import numpy as np
#输入A、b
n=3
A=np.array([[10,-2,-1],
            [-2,10,-1],
            [-1,-2,5]],dtype=np.float_)
b=np.array([[3,15,10]]).T

#初始化xk+1，xk
xk=np.zeros((3,1),dtype=np.float_)
xk1=np.zeros((3,1),dtype=np.float_)



#生成D、L、U
D=np.diag(np.diag(A))
#  D=[[10,0,0],
#     [0,10,0],
#     [0,0,5]]

L=-np.tril(A,k=-1)
#  L=[[0,0,0],
#     [-1,0,0],
#     [-1,-1,0]]

U=-np.triu(A,k=1)
#  U=[[0,-1,-2],
#     [0,0,-2],
#     [0,0,0]]


#计算M、g
M=np.linalg.inv(D)@(L+U)
g=np.linalg.inv(D)@b

#迭代
for i in range(10):
    print(xk)
    xk1=M@xk+g
    xk=xk1

#结果
print(xk.T)