目录
实验六:矩阵特征值与迭代法
练习二
1. 已知线性方程组Ax=b,其中A是50×50如下形式的矩阵
b为任意的非零向量。按要求完成下面的实验任务:
(1)方程组Ax=b是否有解与右边常向量b是否有关?
无关。不论右边的常量为何值该方程组均有解。(该方阵为满秩方阵)
(2)对此方程组,雅可比和高斯-赛德尔迭代法是否收敛?若收敛,则任意取定一个非零向量b,分别用雅可比和高斯-赛德尔迭代法计算方程组的近似解;
clc;clear;
upp=-ones(50,1);
up=-2*ones(50,1);
d=12*ones(50,1);
a=spdiags([upp,up,d,up,upp],[2,1,0,-1,-2],50,50);
A=full(a);
%将此矩阵分解为d,l和u
d=full(spdiags(d,0,50,50));
l=-full(spdiags([upp,up],[-2,-1],50,50));
u=-full(spdiags([up,upp],[1,2],50,50));
b=ones(50,1);x0=rand(50,1);x1=x0;
%jacobi迭代法
dd1=inv(d)*(l+u);
pbj1=max(abs(eig(dd1)));%0.4981<1,收敛
for i=1:15
x0=dd1*x0+inv(d)*b;
end
x0
%gauss-seidel迭代法
dd2=inv(d-l)*u;
pbj2=max(abs(eig(dd2)));%0.2608<1,收敛
for i=1:15
x1=dd2*x1+inv(d-l)*b;
end
x1
(3) 在相同的精度限制下,查看两种迭代法的迭代次数。
for i=1:30
x0=dd1*x0+inv(d)*b;
if A*x0-b<10^-8
break;
end
end
i
%gauss-seidel迭代法
for i=1:30
x1=dd2*x1+inv(d-l)*b;
if A*x1-b<10^-8
break;
end
end
i
i =
13
i =
3
注:本题的矩阵输入是一个值得学习的地方:学会用spdiags函数进行输入特殊的矩阵。
2. 设有线性方程组Ax=b,其中矩阵
,
建立下面的迭代法解该线性方程组:
要求:
(1)求出使上述迭代法收敛的ω的取值范围;
clc;clear;
A=[-1,1,0;-4,3,0;1,0,2];
b=[1;1;1];
syms w
ddjz=eye(3)-w*A;
tzz=eig(ddjz);
% 1 - w
% 1 - w
% 1 - 2*w
h1=ezplot('max(abs(1-x),abs(1-2*x))');
hold on
h2=ezplot('x','1',[-5,5]);
axis equal
grid on
由图像我们可以得知,若迭代收敛则w属于(0,1)
(2)求出该迭代法收敛速度最快的ω值,并画出该迭代法迭代矩阵的谱半径与ω之间的关系图(答案:最佳ω的取值是2/3);
收敛速度最快即为上图的最小值;
syms x
solve(1-x==2*x-1)
ans =2/3
m=[];
for w=0:0.01:1.4
ddjz=eye(3)-w*A;
m=[m,max(abs(eig(ddjz)))];
end
plot(0:0.01:1.4,m,0:0.01:1.4,m,'.');
grid on
(3) 画出不同ω的取值对应的迭代误差的下降曲线。
clc;clear;
A=[-1,1,0;-4,3,0;1,0,2];
b=[1;1;1];
%选取w为0.1,0.3,2/3,0.8
x=rand(3,1);
w=[0.1,0.3,0.6];
for i=1:3
ddjz=eye(3)-w(i)*A;
m=[];
for n=1:20
x=ddjz*x+w(i)*b;
wz=abs(A*x-b);
m=[m,wz];
end
figure(1)
plot(1:20,m);
end
legend('w=0.1','w=0.3','w=0.6');
推荐下一篇文章:数学实验第三版(主编:李继成 赵小艳)课后练习答案(七)(1)(2)https://blog.csdn.net/2301_80199493/article/details/136043637?spm=1001.2014.3001.5501本文由作者自创,由于时间原因,难免出现一些错误,还请大家批评指正;创作不易,希望大家一键三连!