迭代法 求解线性方程组

迭代法 求解线性方程组 （MATLAB)

统筹了 李庆扬《数值分析》第五版中关于求解Ax=b的四种常用迭代法
一码多用
Jacobi、Gauss-Seidel、SOR、SSOR四种迭代法
可以自行选择迭代方法，自定义精度，选择收敛判定方案
交互式软件般的体验

代码如下

clear
%% 程序说明
% 迭代法求解Ax=b. 可选求解器模型
% A为n阶系数矩阵;b为常数项，err指定精度
% R为迭代矩阵B的谱半径，k为迭代步数，x为最终解

%% 输入端
while 1
A=input('请输入方程组系数矩阵A: \n');
b=input('请输入方程组系数b（行向量）:\n');b=b';
% 精度选择
disp('选择默认精度，请输入0；自定义请输入1');
chr=input("");
if (chr ~= 0)
err=input('请输入指定精度:');
else
    err=1e-7;
end
%迭代次数选择
disp('选择默认最大迭代次数，请输入0;自定义请输入1');
chr=input("");
if (chr ~= 0)
N=input('请指定最大迭代次数:');
else
    N=1e6;
end
[m,n]=size(A);
if m~=n
    disp('error,矩阵A的行数和列数必须相同!'); 
    return;
end
if m~=size(b)
    disp('error,b的大小必须和A的行数和列数相同');
    return;
end
if rank(A)~=rank([A,b])
    disp('error,A的矩阵的秩和增广矩阵的秩不相同，方程不存在唯一的解');
    return;
end
if m==n
    break;
end
end

%% 初始化
x=zeros(n,1);x0=zeros(n,1);
k=0;

%% 迭代模型选择
fprintf('\n');
disp("请输入迭代模型对应的数字，以选择该模型");
fprintf('1、Jacobi 迭代法 \n');
fprintf('2、Gauss-Seidel 迭代法 \n');
fprintf('3、SOR 迭代法 \n');
fprintf('4、SSOR 迭代法 \n');
while 1
choose=input('请输入模型对应数字:');
if choose==1||choose==2
    break;
else
    if choose==4||choose==3
         w=input('请预先指定一个松弛因子:');
         break;
    else
        disp('error，所选数字不在模型范围内');
        choose=input('请重新输入模型对应数字:');
    end
end
end

%% 判敛方案选择
fprintf('\n');
disp('请选择判敛方案');
fprintf('0、默认方案 \n');
fprintf('1、与精确解作比较 \n');
fprintf('2、与前一步近似解作比较 \n');
while 1
example=input('请输入方案对应数字:');
if example==2||example==0
    break;
else
    if example==1
    xstar=input('请输入精确解（行向量)：\n');
    xstar=xstar';
    break;
    else
        disp('error，所选数字不在方案范围内');
        choose=input('请重新输入方案对应数字:');
    end
end
end

%% 判敛范数选择
fprintf('\n');
disp("请指定判敛准则的向量范数p（p>0)");
while 1
disp('如果为无穷范数请在下栏输入:inf');
p=input('请输入范数值(正数）:');
if p<=0
    disp('error,p必须是一个正数');
else
    break;
end
end

%% 具体模型应用
% 将A进行DLU分裂
% B为迭代矩阵；f迭代式的次要结构，R为谱半径
D=diag(diag(A));L=tril(-A,-1);U=triu(-A,1);
switch choose
    case 1
        %% Jacobi 迭代法
B=D\(L+U);
f=D\b;
 R=max(abs(eig(B)));
  if R>=1
     disp('error，迭代不收敛');
  end
while 1
    x0=x;
    x=B*x0+f;
    k=k+1;
    if example==1
      if norm(x-xstar,p)<err
        break;
      end
    else
      if norm(x-x0,p)<err
        break;
      end
    end
     if k>N
        break;
    end
end

case 2
    %% Gauss-Seidel 迭代法 
B=(D-L)\U;
f=(D-L)\b;
R=max(abs(eig(B)));
   if R>=1
      disp('error，迭代不收敛');   
   end
while 1
    x0=x;
    x=B*x0+f;
    k=k+1;
    if example==1
      if norm(x-xstar,p)<err
        break;
      end
    else
      if norm(x-x0,p)<err
        break;
      end
    end
     if k>N
        break;
    end
end

    case 3
        %% SOR 迭代法 
% w的输入控制
while  1 
    if w>0 && w<2
       % 用谱半径判断收敛性 不满足敛散性退出程序并报错
       B=(D-w*L)\((1-w)*D+w*U);
       f=(D-w*L)\(w*b);
       R=max(abs(eig(B)));
       if R<1
          break;
       else
          disp('错误，R>1迭代不收敛');
          disp('请尝试重新输入一个w:');
          w=input('请输入w:');
       end
    else 
       disp('请尝试重新输入一个w∈(0,2}:');
       w=input('再次输入w:');
    end  
end
% 迭代主体
while 1
    x0=x;
    x=B*x0+f;
    k=k+1;
    if example==1
      if norm(x-xstar,p)<err
        break;
      end
    else
      if norm(x-x0,p)<err
        break;
      end
    end
    if k>N
        break;
    end
end  

    case 4
      %% SSOR 迭代法 
% w的输入控制
while  1 
    if w>0 && w<2
       % 用谱半径判断收敛性 不满足敛散性退出程序并报错
       B1=(D-w*L)\((1-w)*D+w*U);
       B2=(D-w*U)\((1-w)*D+w*L);
       B=B2*B1;
       f1=(D-w*L)\(w*b);
       f2=(D-w*U)\(w*b);
       R=max(abs(eig(B)));
       if R<1
          break;
       else
          disp('错误，R>1迭代不收敛');
          disp('请尝试重新输入一个w:');
          w=input('请输入w:');
       end
    else 
       disp('请尝试重新输入一个w∈(0,2}:');
       w=input('再次输入w:');
    end  
end
% 迭代主体
while 1
    x0=x;
    x1=B1*x0+f1;
    x=B2*x1+f2;
    k=k+1;
    if example==1
      if norm(x-xstar,p)<err
        break;
      end
    else
      if norm(x-x0,p)<err
        break;
      end
    end
     if k>N
        break;
    end
end        
end

%% 输出端
fprintf('\n');
fprintf(' 所选迭代法的迭代矩阵B：\n');
disp(B);
fprintf('\n 迭代矩阵B的谱半径R: %6.2f \n',R);
fprintf('\n 迭代步数k: %d ',k);
fprintf('\n');
fprintf('\n 方程的解x为：\n');
disp(x');

感兴趣的可以自行测试