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

高等数值分析实验一
工物研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方法的执行效率还是很高。

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

3．MINRES方法：
程序如下：
[x,flag,relres,iter,resvec] = minres(A,b,1e-6,1043);
semilogy(0:iter,resvec/norm(b),'-o');
title('Convergence curve');
xlabel('iteration number');
ylabel('relative residual');
绘得的收敛曲线如下：
二．当b是由A的部分特征向量生成时，三种方法的收敛情况：[U S]=svd(A); %对A进行特征值分解
n=10;
C=rand(n,1); %随机生成向量C
b=U(:,1:n)*C; %构造b由A的前n个特征相量生成
(1).CG的结果如下：三幅图分别对应n=10，100，1000的情况
(2)Lanczos的结果为：
以下三幅图分别对应n=10，100，1000的收敛曲线
下面的四幅图是念0，100，1000和原来的一般情况下norm(I-Q(k)’*Q(k)随k的增长情况，由图中可以看出，当n较大和在原来的情况下，Q(k)已经和一个正交阵相差较远，即k 较大时，它的各列已不能看成是正交的。

(3).MINRES方法的结果为：
结论：
由各种方法的结果均可以看出，当b由A的n个特征相量生成时，CG收敛性明显较好，并且n越小越好，当n接近A的阶数时，叠代次数趋向于一般情况.
三．利用三种方法解对称不定A的线性方程组，给出收敛曲线。

还是利用上述的1043阶的矩阵A，对它进行特征值分解后，修改它的特征值为有正有负后，再构造新的A。

[U,S]=svd(A); %对A特征值分解
g=10; %通过g可以改变负特征值的位置
for(i=1:5) %同过igai变特征值的个数
S(g+i,g+i)=-S(g+i,g+i);
end
A=U*S*U'; %重新构建A
(1). CG运行结果：
b=ones(length(A),1)时，
方法很快中断，不能够得到方程的解，如下图所示：
当b由A的10个特征相量构成时，方法能够得到方程的解：如下图
用CG方法解对称不定方程时，可能发生中断，而得不到方程的解。

所以CG方法对于对称不定方程不能普遍适用。

(2). Lanczos方法：
用Lanczos方法解对称不定方程，残量会产生振荡收敛的情况，振荡点代表Tk接近奇异，图中有5个振荡峰，这是由于A有5个负的特征值造成的，这说明Lanczos方法在解对称不定方程时不时最优的，没有使残量单调下降.
(3).MINRES方法
可以看出，解对称不定方程的MINRES方法的残量是单调下降的。

下面构造A只有5个不同的特征值，其中有2个是负的，以验证结论：当A只有k个不同的特征值时，MINRES方法至多k步找到准确解：。