清华大学高等数值分析实验设计及答案

高等数值分析实验一

工物研13 成彬彬2004310559

一．用CG，Lanczos和MINRES方法求解大型稀疏对称正定矩阵Ax=b

作实验中，A是利用A= sprandsym(S,[],rc,3)随机生成的一个对称正定阵，S是1043阶的一个稀疏阵

A= sprandsym(S,[],0.01,3);

检验所生成的矩阵A的特征如下：

rank(A-A')=0 %即A=A’，A是对称的；

rank(A)=1043 %A满秩

cond(A)= 28.5908 %A是一个“好”阵

1．CG方法

利用CG方法解上面的线性方程组

[x,flag,relres,iter,resvec] = pcg(A,b,1e-6,1043);

结果如下：

Iter=35，表示在35步时已经收敛到接近真实x

relres= norm(b-A*x)/norm(b)= 5.8907e-007为最终相对残差

绘出A的特征值分布图和收敛曲线：

S=svd(A); %绘制特征值分布

subplot(211)

plot(S);

title('Distribution of A''s singular values');;

xlabel('n')

ylabel('singular values')

subplot(212); %绘制收敛曲线

semilogy(0:iter,resvec/norm(b),'-o');

title('Convergence curve');

xlabel('iteration number');

ylabel('relative residual');

得到如下图象：

为了观察CG方法的收敛速度和A的特征值分布的关系，需要改变A的特征值：

(1).研究A的最大最小特征值的变化对收敛速度的影响

在A的构造过程中，通过改变A= sprandsym(S,[],rc,3)中的参数rc(1/rc为A的条件数)，可以达到改变A的特征值分布的目的：

通过改变rc=0.1，0.0001得到如下两幅图

以上三种情况下，由收敛定理2.2.2计算得到的至多叠代次数分别为：48，14和486，由于上实验结果可以看出实际叠代次数都比上限值要小较多。

由以上三图比较可以看出，A的条件数越大，即A的最大最小特征值的差别越大，叠代所需要的步骤就越多，收敛越慢。

(2)研究A的中间特征值的分布对于收敛特性的影响：

为了研究A的中间特征值的分布对收敛速度的影响，进行了如下实验：

固定A的条件数，即给定A的最大最小特征值，改变中间特征值得分布，再来生成A，具体的实现方法是，先将原来的生成A进行特征值分解：

[U,S]=svd(A);

靠近最小特征值，得到对角阵C1。再通过如下命令得到经处理后的矩阵A3：

A3=U*C1*U’;

对A3重新利用CG方法求解线性方程组得到如下结果图：

当他们都更靠近最小特征值时，结果如下：

靠近最大特征值，得到对角阵C2。再通过如下命令得到经处理后的矩阵A4：

A4=U*C2*U’;

对A4重新利用CG方法求解线性方程组得到如下结果图：

当他们都更靠近最大特征值时，结果如下：

从以上实验基本可以得出以下结论，A的中间特征的分布对收敛速度的影响是很大的，其影响规律可以总结为：当中间特征值越靠近两头时，也就是说集中到一端时，收敛速度越快，当中间特征值的分布越对称时，收敛速度越慢，这里只能给出实验结果，还没能从数学理论上证明结论的正确性，尚待商讨。

2．Lanczos方法

自己编制Lanczos方法的程序如下：

clear

load A%A is the same as in the CG Method

b=ones(length(A),1);

bn=norm(b);

q(:,1)=zeros(length(b),1);q(:,2)=b/bn;

r(:,1)=b;bata=1;

for k=(1:length(b)) %Realize k-step Lanczos process

r(:,k+1)=A*q(:,k+1)-bata*q(:,k);

alfa=q(:,k+1)'*r(:,k+1);

r(:,k+1)=r(:,k+1)-alfa*q(:,k+1);

rm=norm(r(:,k+1));

if(rm>1e-6)

bata=rm;

q(:,k+2)=r(:,k+1)/bata;

end

e(k)=norm(eye(k)-q(:,2:k+1)'*q(:,2:k+1));

T=q(:,2:k+1)'*A*q(:,2:k+1); %Construct Tk

T1=svd(T);

if (T1(k)>1e-12) %T is not singular

t=inv(T);

rn(k)=t(k,1)*bata; %The relative residual is lie on the

if abs(rn(k))<3e-6 %Last element of inv(T)'s first column

yk=t(:,1)*bn; %Once reach the required precision,stop

break %and calculate yk,then x.

end

end

end

x=q(:,2:k+1)*yk;

semilogy(e);

程序运行及调试情况：一开始将残量rn定义在1e-6时认为收敛，结果发现程序收敛不了，单步运行时发现执行到k=33时残量就开始一直增大，跟踪Tk发现这时它已经和三对角阵误差很大，例如有些本来应改为0的元素已经达到1e-6的量级，再要求最终结果精度在1e-6就不现实了。如果强制将Tk转化为三对角(即，把该为0的地方清0)，则不会出现这种情况，结果收敛。也就是说该程序达不到任意期望的精度。

不过运行时发现程序执行的时间比直接用inv(A)*b命令的执行时间还是要短的，虽然不知道此命令使用什么方法解的，但是说明Lanczos方法的执行效率还是很高。

解得如下结果(收敛曲线)：可看出曲线后面参量已经开始上跳。