数值分析——雅可比（Jacobi）迭代法（Python及MATLAB实现）_雅可比迭代法matlab

一、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):
        D_inv=np.linalg.inv(D)
        D_inv=np.mat(D_inv)
        B=D_inv * (L+U)
        B=np.mat(B)
        f=D_inv*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.00000001
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:
[[ 10   3   1]
 [  2 -10   3]
 [  1   3  10]]
 
原值矩阵b:
[[14]
 [-5]
 [14]]
 
迭代收敛
迭代次数
22

二、MATLAB实现

function [x,h,k] = Jacobiiter(A,b,x0,N,p)  %N代表最大迭代次数
n=length(A);
n1=length(b);
D=zeros(n);L=zeros(n);U=zeros(n);
x=zeros(n,1);
x0=transpose(x0);
b=transpose(b);
if n~=n1
    disp('维数不一致')
    return
end
D(n,n)=A(n,n);
for i=1:n-1
    D(i,i)=A(i,i);
    for j=1:i-1
        L(i,j)=-A(i,j);
    end
    for r=i+1:n-1
        U(i,r)=-A(i,r);
    end
end
k=0;
while k<=N
    x=(eye(n)-inv(D)*A)*x0+inv(D)*b;
    k=k+1;
    h=norm((x-x0),inf);
    if h<=p
        break
    end
    
    x0=x;
    
end
if k>N
    disp(['迭代次数= ',num2str(k),',算法超过最大迭代次数。']);
    
else
    disp(['迭代次数= ',num2str(k)]);
    
end
x=x0
end

输入

 A=[10,3,1;2,-10,3;1,3,10];
b=[14,-5,14];
x0=[0,0,0];
N=100;
p=0.00000001;
[x,h,k] = Jacobiiter(A,b,x0,N,p);

输出

 A=[10,3,1;2,-10,3;1,3,10];
b=[14,-5,14];
x0=[0,0,0];
N=100;
p=0.00000001;
[x,h,k] = Jacobiiter(A,b,x0,N,p);