数值分析算法 MATLAB 实践 线性方程组 雅可比迭代法

数值分析算法 MATLAB 实践 线性方程组迭代法

Jacobi迭代法

雅可比迭代法保证收敛的条件是矩阵A(Ax=b)为严格的行对角占优矩阵，对于每一行，对角线上的元素之绝对值大于其余元素绝对值的和。需要说明的是：即使不满足此条件，雅可比法有时仍可以收敛。

%% 雅可比迭代法 [x,k,index] = Jacobimethod(A,b,ep)
% A为方程组的系数矩阵；
% b为方程组的右端项;
% ep为精度要求，缺省值为1e-5;
% it_max为最大迭代次数，缺省值为100;
% x为方程组的解;
% k为迭代次数;
% index为指标变量，index=0表示迭代失败，index=1表示收敛到指定要求
A = [10 3 1; 2 -10 3 ;1 3 10];
b = [14 -5 14 ]'; %b = [14; -5; 14 ];
eps = 0.005;
[x_0,k0_cnt,index] = Jacobimethod(A,b,eps);
disp('迭代次数:k0_cnt=')
disp(k0_cnt)
disp(['方程组的解：x_0 = '])
disp(x_0)

%% 求线性方程组的Jacobi迭代法，调用格式为[x, k] = JacobiFunc(A,b,x0,eps,it_max)
%  其中, A 为线性方程组的系数矩阵，b 为常数项，eps 为精度要求，默认为1e-6,x0迭代初始值
%  it_max 为最大迭代次数，默认为1000
%  x 为线性方程组的解，k迭代次数
x0=[0,0,0]';%[x1;x2;x3]列向量
it_max = 1000;eps=1e-6;
[x1, k1_cnt] = JacobiFunc(A,b,x0,eps,it_max);
disp('迭代次数:k1_cnt=');
disp(k1_cnt)
disp(['方程组的解：x1 = ']);
disp(x1)

 %% [x,k]=JacobiFunmethod(A,b,x0,N,emg)
 % A:线性方程组左端矩阵，b:线性方程组右端向量，x0:迭代初值
 % N:迭代次数上界，若迭代次数大于n，则迭代失败， emg:精度指标
 % k:迭代次数，
 % x:用迭代法求得的线性方程组的近似解
x0=[0,0,0]';%[x1;x2;x3]列向量
it_max = 1000;eps=1e-6;
[x2,k2_cnt]=JacobiFunmethod(A,b,x0,it_max,eps);
disp('迭代次数:k2_cnt=');
disp(k2_cnt)
disp(['方程组的解：x2 = ']);
disp(x2)

function [x,k] = JacobiFunc(A,b,x0,eps,it_max)
%  求线性方程组的Jacobi迭代法，调用格式为[x, k] = JacobiFunc(A,b,x0,eps,it_max)
%  其中, A 为线性方程组的系数矩阵，b 为常数项，eps 为精度要求，默认为1e-6,x0迭代初始值
%  it_max 为最大迭代次数，默认为200
%  x 为线性方程组的解，k迭代次数
  if nargin == 3
      eps = 1.0e-6;
      M = 200;
  elseif nargin<3
      disp('输入参数数目不足3个');
      return
  elseif nargin ==5
      M = it_max;
  end
  D = diag(diag(A));%求A的对角矩阵
  L = -tril(A,-1);%求A的下三角矩阵
  U = -triu(A,1);%求A的上三角矩阵
  B = D\(L+U);
  f = D\b;
  x = B*x0+f;
  k = 1;%迭代次数
  while norm(x-x0)>=eps
      x0 = x;
      x = B*x0+f;
      k = k+1;
      if(k>=M)
          disp('Warning:迭代次数太多,可能不收敛!');
          return;
      end
  end
end

function [ x,k,index]=Jacobimethod(A,b,ep,it_max)
% 求线性方程组的雅可比迭代法，其中，
% A为方程组的系数矩阵；
% b为方程组的右端项;
% ep为精度要求，缺省值为1e-5;
% it_max为最大迭代次数，缺省值为100;
% x为方程组的解;
% k为迭代次数;
% index为指标变量，index=0表示迭代失败，index=1表示收敛到指定要求,
    [n,m] = size(A);nb = length(b);
%当方程组行与列的维数不相等时，停止计算，并输出出错信息。
    if n ~=m
            error('The rows and columns of matrix A must be equal! ');
            return;
    end
    % 当方程组与右端项的维数不匹配时，停止计算，并输出出错信息。
    if m~=nb
            error ('The columns of A must be equal the length of b! ');
            return;
    end

    if nargin<4
        it_max =100;
    end
    if nargin<3 
        ep = 1e-5;
    end
    
    k=0;x = zeros (n,1);y=zeros (n,1);index=1;
    while 1
        for i=1 :n
                y(i) =b(i) ;
                for j=1:n
                    if j~=i
                            y(i) =y(i) -A(i,j)*x(j);
                    end
                end
       if abs(A(i,i))<1e-10  &&k==it_max       % abs绝对值函数
                    index =0 ;return;
                end
                y(i) =y(i)/A(i,i);
        end
        k = k +1;
        if norm(y-x,inf) <ep
            break;
        end
        x = y;
end

function [x,k]=JacobiFunmethod(A,b,x0,N,emg)
 % A:线性方程组左端矩阵，b:线性方程组右端向量，x0:迭代初值
 % N:迭代次数上界，若迭代次数大于n，则迭代失败， emg:精度指标
 % k:迭代次数，
 % x:用迭代法求得的线性方程组的近似解
     n=length(A);      
     x=zeros(n,1); %设置变量
     X=zeros(n,1); % X
     x=x0; k=0; 
     r=max(abs(b-A*x));
     while (r>emg)          % 迭代循环过程
       for i=1:n
           sum=0;
           for j=1:n
                if i~=j
                 sum=sum+A(i,j)*x(j);
                end
           end
           X(i)=(b(i)-sum)/A(i,i);
       end
       r=max(abs(X-x)); 
       x=X; 
       k=k+1;
       if k>N
          disp('迭代失败，返回');
          return;
       end
end