python实现雅克比(Jacobi)迭代法_雅克比迭代求解方程组python

import numpy as np
#e为误差
def Jacobi(A,b,x,e,times=100):
    length, width = np.shape(A)
    D = np.mat(np.diag(np.diag(A)))
    L = np.triu(A, 1)
    U = np.tril(A, -1)
    J = -D.I * (L + U)
    H = np.eye(length) - D.I * A
    eig,_ = np.linalg.eig(H)
    spectral_radius = max(abs(eig))
    if spectral_radius < 1:
        print('此方程组收敛,谱半径为',round(spectral_radius,5))
 
        x0=x
        x = J * x + D.I * b
        #x = D.I*(b-(L+U)*x0)
        k = 1
        while k < times:
            if abs(np.max(abs(x-x0), axis=0))>e:
                x0=x
                x = J * x + D.I * b
                #x=D.I*(b-(L+U)*x0)
                k += 1
            else:
                print('当精度为', e, '时,Jacobi在%d次内收敛' % k)
                break
        print('Jacobi迭代解为\n', x)
B = np.mat([[8,-1,1],[2,10,-1],[1,1,-5]])
b = np.mat('1;4;3')
x = np.mat('0;0;0')
e = 0.001
times = 100
Jacobi(B,b,x,e,times)