myjacobi.m

解线性方程组的迭代法 C语言

function [x,iter] = myjacobi(A, b, x0, epsilon, iter_max)
%  Jacobi 迭代解线性方程组
%  A 为系数矩阵， b为右端项，x0为初值，epsilon为精度，iter_max为最大迭代步数
%

if prod(diag(A))==0
    disp('Error: diag(A) is not invertible!')
    return
end
n = length(b);
x_old = reshape(x0,n,1);
x_new = zeros(n,1);
for k = 1 : iter_max
    for i = 1 : n
        tmp = b(i) - A(i,[1:i-1])*x_old([1:i-1]) - A(i,[i+1:n])*x_old([i+1:n]);
        x_new(i) = tmp/A(i,i); 
        % x_new(i) = x_old(i) + (b(i) - A(i,:)*x_old)/A(i,i);
    end
    % x_new = x_old + (b - A*x_old)./diag(A);
    error1 = norm(x_new - x_old);
    if error1 < epsilon, break, end
    x_old = x_new;
end
iter = k;    
x = x_new;

% END