迭代法
求解线性方程组 (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');