清华大学高等数值分析 第一次实验作业

1、高等数值分析第二次作业 第四题姓名：- 学号：- 时间：11月20日1、构造例子说明 CG 的数值性态. 当步数 = 阶数时 CG 的解如何? 当 A的最大特征值远大于第二个最大特征值, 最小特征值远小于第二个最小特征值时方法的收敛性如何?答：1）构造矩阵如下首先，构造一个非常病态的矩阵。构造矩阵的条件数K(A)为1015，可以得到收敛曲线如下图所示第二种方法构造矩阵（良态）：条件数不大的矩阵也能够使迭代步数大于等于矩阵阶数，构造一个良态的矩阵如下：使其特征值为等比数列，即由此得到的矩阵A最大特征值为105，最小为1，的条件数是105，特征值分布如下图收敛曲线如下图，发现使用CG法仍然可以收敛

2、，但是收敛速度非常缓慢，即使到1000步也难以达到理想的残差，说明即使问题的条件数较小，即良态的，CG方法也可能收敛非常慢。从上面的结果中可以看出，当矩阵的迭代步数等于阶数时a. 往往收敛曲线不平滑，震荡非常严重，残差的2范数没有最优性；b. 但是方法的相对残差大体上仍然具有下降趋势，最终仍然能够收敛。矩阵的病态性延迟了方法的收敛。c. 条件数较小的矩阵也会出现收敛缓慢的现象，这与特征值的分布有关。2）构造矩阵1002阶的A，其中，A的最大的特征值为109，最小特征值为1，其具有如下的特征值分布图 1 A的特征值分布该矩阵的条件数k=1011，是个病态问题，结果如下图 2 A的收敛性态图代入C

3、G法，得到的收敛曲线如Error! Reference source not found.所示。从图中可以看出，一开始，方法残差下降很慢，等到后面，在150步左右，以另一种收敛速度下降。实验结果与理论结果相同。3） 结论：a. 对于非常病态的问题，CG法的迭代次数可能大于等于阶数，残差震荡非常严重，但是大体趋势仍然是下降的，最终也能够收敛；对于极其病态的问题，CG法也可能无能为力，无法收敛。b. 对于最大特征值远远大于第二大特征值的情况，收敛曲线出现了两种下降速度，说明特征值的分布对于CG方法的收敛速度影响很大。c.2、对于同样的例子, 比较 CG 和 Lanczos 的计算结果.在相同的A下

4、，取b=(1,1,1,1.1)T,x0=(0,0,0.0)T，分别用CG法和Lanczos法求解。停机准则为相对误差小于10-8。将A取以下几种：良态正定矩阵、近似奇异正定矩阵、病态正定矩阵、良态不定矩阵进行实验，结果分别如下： 良态的正定问题取矩阵A的特征值分布如下：lamda=1001, 1000:-1:1, 1从特征值的分布上可以看出，A的条件数为103，是个良态正定的问题。CG结果如下图从上图可以看出，对称正定良态问题，CG和Lanczos方法收敛性态几乎一模一样。 近似奇异的问题上一个A矩阵中的一个特征值减小，取为10-3，特征值分布为从收敛特性上看，对于近似奇异的正定矩阵（数值上不

5、奇异），CG和Lanczos方法得到的结果一致。但是CG的求解速度比Lanczos小很多。 病态的正定问题在以上相同的A（n=1002）的基础上，将最大特征值提高为1e15，则条件数为1e15，是个非常病态的问题，同时使用CG法和Lanczos方法得到的收敛曲线下从结果中可以看出，对于这个问题，Lanczos可能比CG的迭代次数小，但是，程序的运行时间来看，CG法为2.13s，而Lanczos为7.14s，因此对于正定问题CG法的效率要明显高于Lanczos。 良态的不定问题在以上相同的A的基础上，增加两个负特征值-1和-10，在此基础上，使用CG和Lanczos计算，得到的结果如下从图中可以

6、看出，虽然问题比较良性，但是也出现了尖峰，即近似中断现象。其中，CG法也能够在没有中断的情况下求解不定问题，但是如果一旦中断，那么无法得到结果。但是Lanczos中断之后，却能够继续进行下去。即使用Lanczos解不定问题的风险更小。总结：1) 对于正定问题，不论良态还是病态，Lanczos和CG方法求解的结果非常相近，收敛曲线也近似一致，收敛步数也相差不大，这说明两个方法的原理是一致的。但是CG方法的速度明显快于Lanczos法。2) 对于不定问题，CG方法也可能收敛，如果收敛，Lanczos方法和CG方法的到的结果也类似（步数和收敛性态）。但是，如果CG方法一旦中断，前面的所有运算都白费了

7、，是恶性中断，而Lanczos仍然能够继续。因此，再不定问题中，Lanczos比CG方法的风险小。3、当A 只有m个不同特征值时, 对于大的m和小的m, 观察有限精度下 Lanczos方法如何收敛.答：对于不同的m值进行实验，得到的结果如下：当m=20时，取特征值分别为1,2,3,4.20，每个特征值的重数均为60重，得到结果为两个方法均在20步的时候收敛，即nitm。当m=50时，取特征值分别为1,2,3,4.50，每个特征值的重数均为24重，得到的结果为在38步左右收敛，CG法与Lanczos方法的收敛速度一致。同样也满足nitm。当m=500时，取特征值分别为1,2,3,4.500，每个

8、特征值的重数均为2重，得到的结果为从结果中可以看出，迭代120多步就收敛了，而且CG法与Lanczos法的结果一致。同样满足nitm。当m=1000时，取特征值分别为1,2,3,4.100，每个特征值的重数均为1重，得到的结果为经过多次实验，将得到的迭代次数与不同特征值个数的关系做图如下，从图中可以看出，当m很小时，迭代次数等于不同特征值个数，随着m增大，迭代次数也不断增大，但是增大的趋势变缓，最终，几乎不随着m增大而增大。同时，发现，迭代次数始终小于不同特征值的个数，即与书上讲的“若A有m个不同的特征值，则至多m步收敛”一致！注：图中蓝色曲线为y=x曲线。总结：1) 特征值的分布对于Lanc

9、zos方法的收敛性态影响非常大；2) 如果A有m个不同的特征值，那么就会至多m步收敛。4、取初始近似解为零向量, 右端项 b仅由 A 的 m个不同个特征向量的线性组合表示时, Lanczos 方法的收敛性如何? 数值计算中方法的收敛性和 m的大小关系如何?答：首先，确定一个比较良态的矩阵A，取其特征值为lamda=3100:10:2500, 1000:-1:1分布如下图分别取b为bm=Qm(1,1,1.1)mT其中Qm的列为A的特征向量。变化分别取m的值为5,10,15,30,40,50,80,90,100,200,300,400,500,600,800,900,1000分别得到收敛曲线，下面

10、仅列出m=5、15、40、50、90、100、400、600、900时的收敛曲线。从图中可以看出，各个收敛性态都很好，曲线比较平滑。通过实验，可以得到m与收敛迭代步长的曲线如下图注：图中蓝色曲线为y=x曲线。从图中可以看出，收敛性态与特征向量有关！随着m的增大，迭代次数有上升趋势，即若b仅有少数几个特征向量的线性组合而成，那么可以得到更好的收敛效果。注：图中蓝色曲线为y=x曲线。总结：1) 收敛特性与特征向量有关；2) 若b由m个特征向量线性组合而成，m越小，收敛速度越快；3) 理论上，至多m步收敛，但是实验中并不是严格按照这个关系，说明特征向量问题比特征值问题复杂一些。5、构造对称不定的矩阵

11、, 验证 Lanczos 方法的近似中断, 观察收敛曲线中的峰点个数和特征值的分布关系; 观察当出现峰点时, MINRES 方法的收敛性态怎样.答：取b=(1,1,1,.,1,1)T，x0=(0,0,0,.,0,0)T，停机准则为10-8。首先，取特征值为lamda=3100:10:2500, 1000:-1:1,-1 ;其中有一个负特征值-1，分别使用MINRES和Lanczos方法得到的曲线如下图从图中可以看出，在40步左右Lanczos发生近似中断现象，然而MINRES方法的相对误差单调不增，没有出现尖峰。我们再多实验几组数据。在前面的特征值分布基础上增加一个负特征值-10得到的收敛曲线

12、如下图可以看出，MIRES方法的相对误差仍然具有单调不增的特性。再将特征值取得更多，再增加10个负特征值，-10：-10：-100，得到的曲线如下通过对比可以发现，负特征值越多，Lanczos方法的纹波越大，而MINRES方法却仍然维持单调不减的特性。总结：1) 对于Lanczos方法，随着负特征值的增多，振荡越来越严重，发生近似中断的次数越来越多；2) 然而，对于相同的A，Minres方法的相对残差没有出现峰值，随着迭代数增加而单调不减。收敛性态比Lanczos好。程序代码CG法 CGllb11.m% CG method for solving Ax=b% By Liu libin 2011

13、-11-06% Email: function x,Error,i,flag=CGllb11(A,b,x,ErrSet,uplimit) %deal with input datam,n=size(b);if mn b=b;endm,n=size(x);if mn x=x;end%CG method beginsr=b-A*x;p=r;%Loop begini=1; temp_rkrkplus=r*r;Error=sqrt(temp_rkrkplus)/norm(b,2);while 1 temp_AP=A*p; temp_rkrk=temp

14、_rkrkplus; temp_pAP=p*temp_AP; if abs(temp_pAP)1e-12 disp(恶性中断！) break; end a=temp_rkrk/(temp_pAP); x=x+a*p; r=r-a*temp_AP; temp_rkrkplus=r*r; beta=temp_rkrkplus/temp_rkrk; p=r+beta*p; % decide the error Err=sqrt(temp_rkrkplus)/norm(b,2); %/(norm(b)+temp_AP); if Erruplimit flag=0; break end endEndLa

15、nczos法 lanczosllb.m% Lanczos method for solving Ax=b% By Liu libin 2011-11-15% Email: function T,Q,x,k,Errs,tiao=Lanczosllb(A,b,x0,Err) %deal with input datax=;tiao=;m,n=size(b);if mn b=b;endm,n=size(x0);if mn m=n; x0=x0;end size(b)size(A*x0)r=b-A*x0;r_zeros=r;q=r/norm(r);q

16、0=0;beta0=0;T=0;k=1;Q=q;Errs=norm(r,2)/norm(b,2);%solve for tri diag matrix Twhile 1 %clc disp(k=);disp(k); T=T,zeros(k,1);zeros(1,k),0; %Add a row and a collom to T r=A*q-beta0*q0; T(k,k)=q*r; r=r-T(k,k)*q; T(k+1,k)=norm(r,2); beta0=T(k+1,k); T(k,k+1)=beta0; q0=q; q=r/beta0; Tk=T(1:k,1:k) ; %solve

17、for ykif k=1 disp(跳过此步) tiao=tiao,k; else L,U=LanczosLU(Tk); bk=norm(r_zeros,2)*1;zeros(k-1,1); %求Ux Ux=zeros(k,1); Ux(1)=bk(1); for i=2:k Ux(i)=bk(i)-L(i,i-1)*Ux(i-1); end %求ym ym=zeros(k,1); ym(k)=Ux(k)/U(k,k); for i=k-1:-1:1 ym(i)=(Ux(i)-ym(i+1)*U(i,i+1)/U(i,i); end Errs=Errs,abs(ym(k)*beta0)/nor

18、m(b,2); if abs(ym(k)*beta0)/norm(b,2)Err x=x0+Q*ym; disp(method converge, got the accurate result!) break; else if beta01e-10 x=x0+Q*ym; disp(Good Luck!) end end %x=x0+Q*ym; end if k=1.1*m %x=x0+Q*ym; disp(Not Converge!) break; end Q=Q,q; k=k+1;endendMINRES法 MINRES.m% Minres method for solving Ax=b%

19、 By Liu libin 2011-11-15% Email: function T,Q,x,k,Errs=minresllb(A,b,x0,Err) %deal with input datam,n=size(b);if mn b=b;endm,n=size(x0);if mn m=n; x0=x0;endr=b-A*x0;r_zeros=r;q=r/norm(r);q0=0;beta0=0;T=0;k=1;Q=q;Errs=norm(r_zeros,2);%solve for tri diag matrix Twhile 1 clc d

20、isp(k=);disp(k); T=T,zeros(k,1);zeros(1,k),0; %Add a row and a collom to T r=A*q-beta0*q0; T(k,k)=q*r; r=r-T(k,k)*q; T(k+1,k)=norm(r,2); beta0=T(k+1,k); T(k,k+1)=beta0; q0=q; q=r/beta0; Q=Q,q;%*Lanzcos process completed! %Interruption for Uplimited times if k=2 else%*start to devide Tk by QR method

21、Rk=T(1:k+1,1:k); Qk=eye(k+1,k); for i=1:k Ci=Rk( i ,i)/(Rk(i,i)2+Rk(i+1,i)2)0.5); Si=Rk(i+1,i)/(Rk(i,i)2+Rk(i+1,i)2)0.5); Rk_temp=Rk; Rk_temp( i ,:)= Ci*Rk(i,:)+Si*Rk(i+1,:); Rk_temp(i+1,:)=-Si*Rk(i,:)+Ci*Rk(i+1,:); Rk=Rk_temp; %compute Qk Qk_temp=Qk; Qk_temp( i ,:)= Ci*Qk(i,:)+Si*Qk(i+1,:); Qk_temp

22、(i+1,:)=-Si*Qk(i,:)+Ci*Qk(i+1,:); Qk=Qk_temp; end %*QR devision completed Errsk=abs(Qk(k+1,1)*norm(r_zeros,2)/norm(b,2); Errs=Errs;Errsk; if ErrskErr | k=m gk=Qk(1:k,1)*norm(r_zeros,2); yk=zeros(k,1); yk( k )=gk( k )/Rk(k,k); yk(k-1)=(gk(k-1)-yk(k)*Rk(k-1,k)/Rk(k-1,k-1); for i=k-2:-1:1 yk(i)=(gk(i)-

23、yk(i+1)*Rk(i,i+1)-yk(i+2)*Rk(i,i+2)/Rk(i,i); end if ErrskErr disp(Method Converge！Congratulations！) else disp(Failed!) end x=x0+Q(:,1:k)*yk; break; end end k=k+1;end end其他附属过程函数A和b的生成矩阵 MakeMatrix.mclcclear%build Matrix A %lamdaglobal Aglobal bbeta=0.5;%lamda=1e9, 1000:-1:1, 0.01;%正定的病态问题lamda=3100:

24、10:2500, 1000:-1:1,-1,-10:-10:-100;%正定的良态问题%lamda=1001, 1000:-1:1, 0.001;%近似奇异问题%lamda=1e15, 1000:-1:1,1;%病态正定问题%lamda=1001, 1000:-1:1,-1,-10,-100;%病态正定问题%lamda=1001, 1000:-1:1, -0.1, -10;%不定的良态问题%lamda= 100:-1:1 , -100:0.5: -1,0 ;%近似不定的良态问题%lamda=100*ones(1,210),-10*ones(1,210),1*ones(1,210),1e3*on

25、es(1,210)%lamda=1:1100, 1e12,0.001;% lamda=1e9,1,2,3,4,5,6,5,8,9,10,11,4,78;%lamda=100,50,20,10,1;n=length(lamda) semilogx(lamda,ones(n,1),bo)Q=orth(randn(n,n);title(A的特征值分布情况)xlabel(特征值)set(gca, ytick,) b=ones(n,1);B=mdiagllb(lamda); A=Q*B*Q; disp(矩阵生成完毕！ 可以进行后续计算！)disp(A(1:5) 对称三对角矩阵的LU分解 LanczosLU.mfunction L,U=LanczosLU(A)m,n=size(A);L=eye(n,n);U=eye(n,n); U(1,1)=A(1,1);%L(2,1)=A(2,1)/A(1,1);U(1,2)=A(1,2); for k=2:m-1 L(k,k-1)=A(k,k-1)/U(k-1,k-1); U(k,k)=A(k,k)-L(k,k-1)*U(k-1,k)