清华大学高等数值分析第一次实验作业Word文档下载推荐.docx

1、lamda=1001, 1000:-1:1, 1从特征值的分布上可以看出，A的条件数为103，是个良态正定的问题。CG结果如下图从上图可以看出，对称正定良态问题，CG和Lanczos方法收敛性态几乎一模一样。近似奇异的问题上一个A矩阵中的一个特征值减小，取为10-3，特征值分布为从收敛特性上看，对于近似奇异的正定矩阵（数值上不奇异），CG和Lanczos方法得到的结果一致。但是CG的求解速度比Lanczos小很多。病态的正定问题在以上相同的A（n=1002）的基础上，将最大特征值提高为1e15，则条件数为1e15，是个非常病态的问题，同时使用CG法和Lanczos方法得到的收敛曲线下从结果中可

2、以看出，对于这个问题，Lanczos可能比CG的迭代次数小，但是，程序的运行时间来看，CG法为2.13s，而Lanczos为7.14s，因此对于正定问题CG法的效率要明显高于Lanczos。良态的不定问题在以上相同的A的基础上，增加两个负特征值-1和-10，在此基础上，使用CG和Lanczos计算，得到的结果如下从图中可以看出，虽然问题比较良性，但是也出现了尖峰，即近似中断现象。其中，CG法也能够在没有中断的情况下求解不定问题，但是如果一旦中断，那么无法得到结果。但是Lanczos中断之后，却能够继续进行下去。即使用Lanczos解不定问题的风险更小。总结：1）对于正定问题，不论良态还是病态，

3、Lanczos和CG方法求解的结果非常相近，收敛曲线也近似一致，收敛步数也相差不大，这说明两个方法的原理是一致的。但是CG方法的速度明显快于Lanczos法。2）对于不定问题，CG方法也可能收敛，如果收敛，Lanczos方法和CG方法的到的结果也类似（步数和收敛性态）。但是，如果CG方法一旦中断，前面的所有运算都白费了，是恶性中断，而Lanczos仍然能够继续。因此，再不定问题中，Lanczos比CG方法的风险小。3、当A 只有m个不同特征值时, 对于大的m和小的m, 观察有限精度下 Lanczos方法如何收敛.对于不同的m值进行实验，得到的结果如下：当m=20时，取特征值分别为1,2,3,4

4、.20，每个特征值的重数均为60重，得到结果为两个方法均在20步的时候收敛，即nitm。当m=50时，取特征值分别为1,2,3,4.50，每个特征值的重数均为24重，得到的结果为在38步左右收敛，CG法与Lanczos方法的收敛速度一致。同样也满足nitm。当m=500时，取特征值分别为1,2,3,4.500，每个特征值的重数均为2重，得到的结果为从结果中可以看出，迭代120多步就收敛了，而且CG法与Lanczos法的结果一致。同样满足nitm。当m=1000时，取特征值分别为1,2,3,4.100，每个特征值的重数均为1重，得到的结果为经过多次实验，将得到的迭代次数与不同特征值个数的关系做图

5、如下，从图中可以看出，当m很小时，迭代次数等于不同特征值个数，随着m增大，迭代次数也不断增大，但是增大的趋势变缓，最终，几乎不随着m增大而增大。同时，发现，迭代次数始终小于不同特征值的个数，即与书上讲的“若A有m个不同的特征值，则至多m步收敛”一致！注：图中蓝色曲线为y=x曲线。1）特征值的分布对于Lanczos方法的收敛性态影响非常大；2）如果A有m个不同的特征值，那么就会至多m步收敛。4、取初始近似解为零向量, 右端项 b仅由 A 的 m个不同个特征向量的线性组合表示时, Lanczos 方法的收敛性如何? 数值计算中方法的收敛性和 m的大小关系如何?首先，确定一个比较良态的矩阵A，取其特

6、征值为lamda=3100:10:2500, 1000: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分别得到收敛曲线，下面仅列出m=5、15、40、50、90、100、400、600、900时的收敛曲线。从图中可以看出，各个收敛性态都很好，曲线比较平滑。通过实验，可以得到m与收敛迭代步长的曲线如下图从图中可以看出，收敛性态与特征向量有关！随着m的增大，迭代次数有上升趋势，即若b仅有少数几个特征向量的线性组合而成，那么可以得

7、到更好的收敛效果。1）收敛特性与特征向量有关；2）若b由m个特征向量线性组合而成，m越小，收敛速度越快；3）理论上，至多m步收敛，但是实验中并不是严格按照这个关系，说明特征向量问题比特征值问题复杂一些。5、构造对称不定的矩阵, 验证 Lanczos 方法的近似中断, 观察收敛曲线中的峰点个数和特征值的分布关系; 观察当出现峰点时, MINRES 方法的收敛性态怎样.取b=（1,1,1,.,1,1）T，x0=（0,0,0,.,0,0）T，停机准则为10-8。首先，取特征值为1,-1 ;其中有一个负特征值-1，分别使用MINRES和Lanczos方法得到的曲线如下图从图中可以看出，在40步左右La

8、nczos发生近似中断现象，然而MINRES方法的相对误差单调不增，没有出现尖峰。我们再多实验几组数据。在前面的特征值分布基础上增加一个负特征值-10得到的收敛曲线如下图可以看出，MIRES方法的相对误差仍然具有单调不增的特性。再将特征值取得更多，再增加10个负特征值，-10：-10：-100，得到的曲线如下通过对比可以发现，负特征值越多，Lanczos方法的纹波越大，而MINRES方法却仍然维持单调不减的特性。1）对于Lanczos方法，随着负特征值的增多，振荡越来越严重，发生近似中断的次数越来越多；2）然而，对于相同的A，Minres方法的相对残差没有出现峰值，随着迭代数增加而单调不减。收

9、敛性态比Lanczos好。程序代码CG法 CGllb11.m% CG method for solving Ax=b% By Liu libin 2011-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）; x=x%CG method beginsr=b-A*x;p=r;%Loop begini=1;temp_rkrkplus=r*r;Error=sqrt（temp_rkrkplus）/norm

10、（b,2）;while 1 temp_AP=A*p; temp_rkrk=temp_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 beta=temp_rkrkplus/temp_rkrk; p=r+beta*p; % decide the error Err=sqrt（temp_rkrkplus）/norm（b,2）; %/（norm（b）+temp_AP）; if

11、 Erruplimit flag=0; breakEndLanczos法 lanczosllb.m% Lanczos method for solving Ax=b% By Liu libin 2011-11-15function T,Q,x,k,Errs,tiao=Lanczosllb（A,b,x0,Err）x=;tiao=;m,n=size（x0）;n m=n; x0=x0 size（b）size（A*x0）r=b-A*x0;r_zeros=r;q=r/norm（r）;q0=0;beta0=0;T=0;k=1;Q=q;Errs=norm（r,2）/norm（b,2）;%solve for

12、tri diag matrix T %clck=）;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-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 for ykif k=1 跳过此步 tiao=tiao,k; else L,U=LanczosLU（Tk）; bk=norm（r_zeros,2）*1;

13、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）; %求ym ym=zeros（k,1）; ym（k）=Ux（k）/U（k,k）; for i=k-1:1 ym（i）=（Ux（i）-ym（i+1）*U（i,i+1）/U（i,i）; Errs=Errs,abs（ym（k）*beta0）/norm（b,2）; if abs（ym（k）*beta0）/norm（b,2）Err x=x0+Q*ym;method converge, got the accurate result!

14、else if beta01e-10Good Luck! %x=x0+Q*ym; if k=1.1*mNot Converge! Q=Q,q; k=k+1;MINRES法 MINRES.m% Minres method for solving Ax=bfunction T,Q,x,k,Errs=minresllb（A,b,x0,Err）Errs=norm（r_zeros,2）; clc%*Lanzcos process completed! %Interruption for Uplimited times if k=2%*start to devide Tk by QR method Rk=

15、T（1:k+1,1:k）; Qk=eye（k+1,k）; for i=1: 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（i+1,:）=-S

16、i*Qk（i,:）+Ci*Qk（i+1,: Qk=Qk_temp;%*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: yk（i）=（gk（i）-yk（i+1）*Rk（i,i+1）-yk（

17、i+2）*Rk（i,i+2）/Rk（i,i）;ErrMethod Converge！Congratulations！Failed! x=x0+Q（:,1:k）*yk;其他附属过程函数A和b的生成矩阵 MakeMatrix.mclcclear%build Matrix A %lamdaglobal Aglobal bbeta=0.5;%lamda=1e9, 1000:1, 0.01;%正定的病态问题1,-1,-10:-10:-100;%正定的良态问题%lamda=1001, 1000:1, 0.001;%近似奇异问题%lamda=1e15, 1000:1,1;%病态正定问题1,-1,-10,-1

18、00;1, -0.1, -10;%不定的良态问题%lamda= 100:1 , -100:0.5: -1,0 ;%近似不定的良态问题%lamda=100*ones（1,210）,-10*ones（1,210）,1*ones（1,210）,1e3*ones（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）,boQ=orth（randn（n,n）;title（A的特征值分布情况x

19、label（特征值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,

20、k）; U（k,k+1）=A（k,k+1）-L（k,k-1）*U（k-1,k+1）; k=m; L（m,m-1）=A（m,m-1）/U（m-1,m-1）; 实验控制文件 （调用这一系列函数） Fortest.mclc global b n_Matrix=m;%调用CGuplimit=2000;ticy,error,n=CGllb11（A,b,zeros（n_Matrix,1）,1e-8,uplimit）;tocfigure（2）semilogy（1:length（error）,error,b.-norm（A*y-b）stringsCG= 算法的收敛曲线 （阶数n=stringsCG=strin

21、gsCG,num2str（m）,）;title（stringsCG）迭代次数ylabel（|r_k|/|b|grid on%调用Lanczos ,x,k,Errs=Lanczosllb（A,b,zeros（n_Matrix,1）,1e-8） ;最大的误差为disp（max（A*x-b） hold onsemilogy（ 1:length（Errs）,Errs,r.-%调用系统函数symmlqx,flag,relres,iter,rv = symmlq（A,b,1e-8,500） %系统函数 semilogy（ 1:length（rv）,rv/norm（b,2）,g.-%调用 MINREST,Q,x,k,Errs=minresllb（A,b,zeros（n_Matrix,1）,1e-8）;hold onk-1,Errs,legend（CG,LanczosSymmlqMatlabMinres