gaussseidel.m

线行方程组的Gauss Jacobi 及一个简单的自适应算法

function [x,flag] = gaussseidel(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);
   %1er   = [];
% 主循环   
   while ~flag & it<maxit
       for k = 1:n,                      % Jacobi公式
           p = [1:k-1 k+1:n];            % Gauss-Seidel
           x(k) = ( b(k)-A(k,p)*x(p) ) / A(k,k);
       end
       it   = it + 1;                    % 统计迭代步
       %x    = xnew;                      % 更新
       flag = [ norm(A*x-b)<=ep ] ;      % 设置收敛指标
       %1er   = [ er norm(A*x-b) ];
       %1semilogy(er);
   end