Jacobi迭代法求解方程组Ax=b

原理

A = D - L - U

迭代公式：x^(k+1) = D^-1( L + U ) x^(k) + D^-1 b

MATLAB实现

function [result,i,errory] = Jacobi(A,b,x0)
% Jacobi(A,b,x0) Jacobi迭代法MATLAB实现
% 输入参数：矩阵A b 初值x0
% 输出参数：计算结果x，迭代次数cnt
sizeA = size(A);
if sizeA(1) ~= sizeA(2)
    error('A不是方阵');
end
e = 1e-4;
i = 1;
D = diag(diag(A));
L = D - tril(A);
U = D - triu(A);
x1 = x0;
x2 = D \ (L + U) * x1 + invD * b;
errory(i) = norm(x2 - x1,inf);
while (norm(x2 - x1,inf) > e)
    x1 = x2;
    x2 = invD * (L + U) * x1 + D \ b;
    i = i + 1;
    errory(i) = norm(x2 - x1,inf);
end
result = x2;
errory(i + 1) = norm(x2 - x1,1);
%% 收敛速度
errorx = 0:1:i;
errory = exp(errory);
figure
semilogy(errorx,errory);
end