数值分析5-用SOR法求解方程组Ax=b的matlab程序

题目如图所示：

要求程序中不存系数矩阵A,分别对不同的阶数取w=1.1, 1.2, …,1.9进行迭代；
首先编写根据矩阵的阶数生成系数矩阵的函数，然后再SOR迭代程序中调用此函数；观察系数矩阵的规律，编写函数SOR_A如下：
代码块1：根据阶数生成系数矩阵：

function [A,b]=SOR_A(n)
A1=(-4)*eye(n);
for i=1:nfor j=1:nif i==j+1L(i,j)=1;else L(i,j)=0;end%生成下三角矩阵Lif j==i+1U(i,j)=1;else U(i,j)=0;end%生成上三角矩阵Uend
end
A=A1+L+U;%生成系数矩阵A
b1(1)=-3;
for i=2:n-1b1(i)=-2;
end
b1(n)=-3;
b=b1;%生成结果向量b

SOR迭代法和Jacobi迭代和Gauss-Seidel迭代法构造思路一致，只是迭代矩阵的构造不同，外加了一个松弛因子w，当w=1时，此迭代法和高斯塞德尔迭代法一致，SOR迭代法代码块如下：
代码块2：SOR迭代法：

%函数SOR_A根据输入的方程阶数自动生成矩阵A
clc;clear;
n=input('请输入系数矩阵A的阶数:n=');
%n=15;
[A,b1]=SOR_A(n);%上机实习题5中的矩阵A；b为行向量，计算时需要转置
%A=input('请输入系数矩阵A:A=');%把上行注释掉后自由度更大，自己定义矩阵A；
%b=input('请输入结果向量b:b=');%把上行注释掉后自由度更大，自己定义矩阵b；
x1=input('请输入初始向量x1:x1=');%x1为行向量，计算时需要转置
%x1=[3 3 3 3 3 3 3 3 3 3 3 3 3 3 3];
%eps=input('请输入停止精度要求:eps=');
eps=1e-6;
w=input('请输入松弛因子:w=');
for i=1:nfor j=1:nif i==jD(i,j)=A(i,j);else D(i,j)=0;end%生成矩阵Dif i>jL(i,j)=-A(i,j);else L(i,j)=0;end%生成矩阵Lif i<jU(i,j)=-A(i,j);else U(i,j)=0;end%生成矩阵Uend
end
D1=inv(D-w*L);%矩阵(D-w*L)的逆矩阵记作D1
D2=((1-w)*D+w*U);
Lw=D1*D2;
f=w*D1*b1';%f为行向量
k=2;
xk=Lw*x1'+f;%xk为列向量，储存时需要转置
x=[x1;xk'];
bk=A*xk;
b=[b1;bk'];
while max(abs(x(k,:)-x(k-1,:)))>eps%if k>500||max(abs(eig(Lw)))>1%如果迭代矩阵B的谱半径大于1，则迭代不收敛error('迭代失败!可能是迭代精度过高或迭代矩阵谱半径大于1的缘故,迭代矩阵谱半径为:%0.8f\n',max(abs(eig(Lw))))elsek=k+1;%第一次为k=4xk=Lw*x(k-1,:)'+f;%x(k-1,:)为X矩阵的第k-1行的行元素x=[x;xk'];%第一次后X矩阵为4行3列；bk=A*xk;b=[b;bk'];end
end
fprintf('经过%d次迭代，方程组的根的近似解为:\n',k-1)
disp(vpa(x(k,:),8))
fprintf('SOR迭代法迭代矩阵的谱半径为：%0.8f\n',max(abs(eig(Lw))))

运行并输入相应变量，结果如下图：

分别输入不同的阶数n和不同的松弛因子w，得到不同的结果，并对比迭代次数，观察迭代因子w与迭代矩阵谱半径和迭代速度的关系

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！