jacobi.m

来自「线行方程组的Gauss Jacobi 及一个简单的自适应算法」· M 代码 · 共 30 行

function [x,flag] = jacobi(A,b,x0,ep,maxit)
% 参数个数判断
   if nargin<5,
       maxit = 1000;
       if nargin<4,
           ep = 1e-8;
           if nargin<3,
               x0 = zeros(size(b));
           end
       end
   end
% 初始化
   x    = x0;
   n    = length(b);
   flag = 0;
   it   = 0;
   xnew = zeros(n,1);
   er   = [];
% 主循环   
   while ~flag & it<maxit
       for k = 1:n,                      % Jacobi公式
           p = [1:k-1 k+1:n];            % Gauss-Seidel
           xnew(k) = ( b(k)-A(k,p)*x(p) ) / A(k,k);
       end
       it   = it + 1;                    % 统计迭代步
       x    = xnew;                      % 更新
       flag = [ norm(A*x-b)<=ep ] ;      % 设置收敛指标
       er   = [ er norm(A*x-b) ];
       semilogy(er);
   end