高等数值分析实验一
工物研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方法的执行效率还是很高。
解得如下结果(收敛曲线)：可看出曲线后面参量已经开始上跳。
