高等数值计算实践题目一1. 实践目的本次计算实践主要是在掌握共轭梯度法,Lanczos 算法与MINRES 算法的基础上,进一步探讨这3种算法的数值性质,主要研究特征值特征向量对算法收敛性的影响。
2. 实践过程(一)生成矩阵(1)作5个100阶对角阵i D 如下:1D 对角元:1,1,...,20,1+0.1(-20),21,...,100j j d j d j j ====2D 对角元:1,1,...,20,1+(-20),21,...,100j j d j d j j ==== 3D 对角元:,1,...,80,81,81,...,100j j d j j d j ====4D 对角元:,1,...,40,41,41,...,60,41+(60),61,...,100j j j d j j d j d j j =====-= 5D 对角元:,1,...,100j d j j ==记i D 的最大模特征值和最小模特征值分别为1i λ和i n λ,则i D 特征值分布有如下特点:1D 的特征值有较多接近于i n λ,并且1/i i n λλ较小,2D 的特征值有较多接近于i n λ,并且1/i i n λλ较大, 3D 的特征值有较多接近于1i λ,并且1/i i n λλ较大,4D 的特征值有较多接近于中间模特征值,并且1/i i n λλ较大, 5D 的特征值均匀分布,并且1/i i n λλ较大(2)随机生成10个100阶矩阵j M :(100(100))j M fix rand =并作它们的QR 分解,得j Q 和j R ,这样可得50个对称的矩阵Tij j i j A Q DQ =,其中i D 的对角元就是ij A 的特征值,若它们都大于0,则ij A 正定,j Q 的列就是相应的特征向量。
结合(1)可知,ij A 都是对称正定阵。
(二)计算结果以下计算,均选定精确解(100,1)exact x ones =,初值0(100,1)x zeros =由ij exact kA x b =计算得到k b (算法中要求解的精度为10e -)。
Lanczos 算法和MINRES 算法均分别使用不带m 循环和有带m 循环进行计算,其中带m 循环时,取15m =进行计算。
(1) 特征值分布对3种算法的影响特征值分布对共轭梯度法的影响特征值分布对()Lanczos m 算法的影响特征值分布对()MINRES m 算法的影响(2) 特征向量矩阵对3种算法的影响特征向量对共轭梯度法的影响(取定2D )特征向量对Lanczos 算法的影响(取定2D )特征向量对()Lanczos m 算法的影响(取定2D )特征向量对MINRES 算法的影响(取定2D )特征向量对()MINRES m 算法的影响(取定2D )(3)从ij A 选取5个画出收敛率曲线 1.共轭梯度法的1k A k Ae k e -曲线以及kAe k 曲线(当迭代步数()j j k n n >表示在步收敛,1k A kAk Ae e e -和由于初始化均为零)nczos 算法的1k A k Ae k e -曲线以及kAe k 曲线(当迭代步数()j j k n n >表示在步收敛,1k A kAk Ae e e -和由于初始化均为零)3.()Lanczos m 算法的1k A k Ak e -曲线以及kAe k 曲线(当迭代步数()j j k n n >表示在步收敛,1k A kAk Ae e e -和由于初始化均为零)4.MINRES 算法的1k A k Ak e -曲线以及kAe k 曲线(当迭代步数()j j k n n >表示在步收敛,1k A kAk Ae e e -和由于初始化均为零)5.()MINRES m 算法的1k A k Ae k e -曲线以及kAe k 曲线(当迭代步数()j j k n n >表示在步收敛,1k A kAk Ae e e -和由于初始化均为零)(三)相关讨论(1)由矩阵i D 的生成方法,可以知道:1D 的特征值有较多接近于i n λ,并且1/i i n λλ较小,2D 的特征值有较多接近于i n λ,并且1/i i n λλ较大,3D 的特征值有较多接近于1i λ,并且1/i i n λλ较大,4D 的特征值有较多接近于中间模特征值,并且1/i i n λλ较大,5D 的特征值均匀分布,并且1/i i n λλ较大。
观察分析计算结果,有以下几点结论:1.当1/i i n λλ较小时,3种算法都收敛的特别快;当1/i i n λλ较大时,3种算法的收敛速度都要慢一些,这与理论分析是一致的。
2.当特征值有较多接近i n λ与特征值有较多接近于1i λ以及特征值有较多接近中间值时,几乎不影响3种算法的收敛速度,而且特征向量也几乎不影响3种算法的收敛速度。
(2)1.在填特征值分布对3种算法的影响的表中,对于ij A 的选择原则:1.1首先由于要研究特征值分布对3种算法的影响,因此选取的ij A 的特征值应该照顾的不同的特征值分布;1.2对于固定的特征值分布i D 对应的ij A 的生成方法可知,由于j M 是随机的,j Q 也是随机的,因此,选取了1j =的情况,即1121314151,,,,A A A A A 。
2.在填特征值向量对3种算法的影响的表中,取定一个i D 的选择原则:考虑矩阵i D 的特征值分布以及对计算结果的观察(特征值向量对3种算法的影响不大),于是取定2D 。
3.在画收敛率曲线图时,对于ij A 的选择原则:同1的选取原则(3)当矩阵阶数n (本次计算100n =)不是很大时,对于对称正定阵,3种算法的收敛速度都比较快,尤其是当1/i i n λλ较小时,3种算法都收敛的特别快,这与课程上的理论分析是一致的。
其次,若选取合适的m (比如本次计算所选取的=15m ),使用Lanczos m ()算法和MINRES m ()算法都有相当令人满意的收敛性质,而且相比共轭梯度法的收敛速度要快得多。
另外,也可以看到共轭梯度法数值上是十分不稳定的,相比较而言,使用Lanczos m ()算法和MINRES m ()算法要稳定的多。
(4)stationary 迭代法中的SOR 以及使用对称的SOR 公式构造的预处理CG 法的迭代结果如下所示:(算法中要求解的精度为5e -)观察分析计算结果:stationary 迭代法中的SOR 的收敛速度要慢很多,而且解的精度不是很高,同时数值上也不是很稳定;使用预处理CG 法,可以加快收敛速度,并且可以增加数值稳定性。
因此,相比较而言,使用预处理CG 法与Lanczos m ()算法和MINRES m ()算法要好得多。
stationary 迭代法中的SOR 算法SOR CG1.stationary 迭代法中的SOR 算法(当迭代步数()j j k n n >表示在步收敛,1k A kAk Ae e e -和由于初始化均为零)2.使用对称的SOR 公式构造的预处理CG 法(当迭代步数()j j k n n >表示在步收敛,1k A kAk Ae e e -和由于初始化均为零)(5)对于非正定的对称正定系数矩阵,比较Lanczos m ()算法和MINRES m ()算法 1.构造非正定的对称系数矩阵,类似正定对称矩阵的构造方式,选取对角阵6D ,对角元:1,1,...,20,(1+(-20)),21,...,100j j d j d j j ===-=,随机矩阵11(100(100))M fix rand =,作它的QR 分解,得11Q ,令11611T AA Q D Q =,则AA 为非正定的对称阵。
2.选定精确解(100,1)exact x ones =,初值0(100,1)x zeros =由ij exact k A x b =计算得到k b (算法中要求解的精度为10e -)。
L a n c z o s 算法和MINRES 算法均分别使用不带m 循环和有带m 循环进行计算,其中带m 循环时,取15m =进行计算。
3.由于系数矩阵非正定,,Axx Ax =<>会出现负值,故不能作为范数,此时考虑使用2-范数来代替,画出Lanczos 算法和MINRES 算法的212k k e k e -曲线以及2ke k 曲线。
4.结合图表可以看出,对于此非正定的对称系数矩阵AA ,2种算法均有很好的收敛性质,不仅收敛速度特别快,而且解的精度可以很高。
同时,可以看出,Lanczos 算法的收敛曲线图在开始时出现了一个较大的峰值,相应的误差突然增大,表现出数值不稳定。
因此,对于本例,MINRES 算法要略优于Lanczos 算法,MINRES 算法要相对稳定得多。
(当迭代步数()j j k n n >表示在步收敛,2212k k k e e -和由于初始化均为零)(当迭代步数()j j k n n >表示在步收敛,2212k k k e e -和由于初始化均为零)DMAQ.m 用于生成矩阵D,M,A,QCG.m 共轭梯度法Lanczos1.m Lanczos算法Lanczos.m Lanczos(m)算法MINRES1.m MINRES算法MINRES.m MINRES(m)算法SOR.m SOR算法SOR_CG.m 对称SOR的预处理共轭梯度法main.m 主函数,包括画图代码NPDS.m 用于非正定的对称系数矩阵的算例Lanczos1_AA.m 用于非正定的对称系数矩阵的Lanczos算法Lanczos_AA.m 用于非正定的对称系数矩阵的Lanczos(m)算法MINRES1_AA.m 用于非正定的对称系数矩阵的MINRES算法MINRES_AA.m 用于非正定的对称系数矩阵的MINRES(m)算法function [D,Q,A]=DMAQa12=linspace(1.1,9,80);d1=[a11,a12];D1=diag(d1);a21=ones(1,20);a22=linspace(2,81,80);d2=[a21,a22];D2=diag(d2);a31=linspace(1,80,80);a32=81*ones(1,20);d3=[a31,a32];D3=diag(d3);a41=linspace(1,40,40); a42=41*ones(1,20);a43=linspace(42,81,40); d4=[a41,a42,a43];D4=diag(d4);d5=1:100;D5=diag(d5);N=100;M1=fix(100*rand(N)); [Q1,R1]=qr(M1);M2=fix(100*rand(N)); [Q2,R2]=qr(M2);M3=fix(100*rand(N)); [Q3,R3]=qr(M3);M4=fix(100*rand(N)); [Q4,R4]=qr(M4);M5=fix(100*rand(N)); [Q5,R5]=qr(M5);M6=fix(100*rand(N)); [Q6,R6]=qr(M6);M7=fix(100*rand(N)); [Q7,R7]=qr(M7);M8=fix(100*rand(N)); [Q8,R8]=qr(M8);M9=fix(100*rand(N)); [Q9,R9]=qr(M9);M10=fix(100*rand(N)); [Q10,R10]=qr(M10);D=cat(3,D1,D2,D3,D4,D5);Q=cat(3,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,Q10);k=0;for i=1:5for j=1:10k=k+1;A(:,:,k)=Q(:,:,j)*D(:,:,i)*Q(:,:,j)';endendendfunction [x,iter,w1,w2]=CG(A,b,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);x=x0;r=b-A*x;p=r;for i=1:N1t1=x-xx;alpha=(r'*r)/(p'*(A*p));r1=r;x=x+alpha*p;r=r-alpha*A*p;belta=(r'*r)/(r1'*r1);p=r+belta*p;t2=x-xx;w1(i,1)=sqrt((t2'*(A*t2))/(t1'*(A*t1))); w2(i,1)=sqrt(t2'*(A*t2));iter=iter+1;if norm(t2)<tolbreak;endendfunction [x,iter,w1,w2]=Lanczos1(A,bb,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);m=N1;w1=zeros(N1,1);w2=zeros(N1,1);x=x0;q=zeros(N1);a=zeros(1,m);b=zeros(1,m);r0=bb-A*x0;b0=norm(r0);q(:,1)=r0/b0;a(1)=q(:,1)'*A*q(:,1);r1=A*q(:,1)-a(1)*q(:,1);b(1)=norm(r1);q(:,2)=r1/b(1);for i=2:N1a(i)=q(:,i)'*A*q(:,i);r1=A*q(:,i)-a(i)*q(:,i)-b(i-1)*q(:,i-1); b(i)=norm(r1);b(i)=b(i);q(:,i+1)=r1/b(i);m=i;Q=q(:,1:m);H=Q'*A*Q;y=H\(Q'*r0);z=Q*y;t1=x-xx;x=x0+z;t2=x-xx;w1(i,1)=sqrt((t2'*(A*t2))/(t1'*(A*t1))); w2(i,1)=sqrt(t2'*(A*t2));iter=iter+1;if norm(t2)<tolbreak;endendendfunction [x,iter,w1,w2]=Lanczos(A,bb,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);x=x0;for j=1:N1m=15;q=zeros(N1);a=zeros(1,m);b=zeros(1,m);r0=bb-A*x;b0=norm(r0);q(:,1)=r0/b0;a(1)=q(:,1)'*A*q(:,1);r1=A*q(:,1)-a(1)*q(:,1);b(1)=norm(r1);q(:,2)=r1/b(1);for i=2:ma(i)=q(:,i)'*A*q(:,i);r1=A*q(:,i)-a(i)*q(:,i)-b(i-1)*q(:,i-1); b(i)=norm(r1);b(i)=b(i);q(:,i+1)=r1/b(i);endQ=q(:,1:m);H=Q'*A*Q;y=H\(Q'*r0);z=Q*y;t1=x-xx;x=x0+z;t2=x-xx;w1(j,1)=sqrt((t2'*(A*t2))/(t1'*(A*t1)));w2(j,1)=sqrt(t2'*(A*t2));iter=iter+1;if norm(t2)<tolbreak;endx0=x;endendfunction [x,iter,w1,w2]=MINRES1(A,bb,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);x=x0;m=N1;q=zeros(N1);a=zeros(1,m);b=zeros(1,m);r0=bb-A*x0;b0=norm(r0);q(:,1)=r0/b0;a(1)=q(:,1)'*A*q(:,1);r1=A*q(:,1)-a(1)*q(:,1);b(1)=norm(r1);q(:,2)=r1/b(1);for i=2:N1a(i)=q(:,i)'*A*q(:,i);r1=A*q(:,i)-a(i)*q(:,i)-b(i-1)*q(:,i-1); b(i)=norm(r1);b(i)=b(i);q(:,i+1)=r1/b(i);m=i;Q=q(:,1:m);HH=zeros(m+1,m);for j=1:mHH(j,j)=a(j);endfor j=1:m-1HH(j,j+1)=b(j);endfor j=2:mHH(j,j-1)=b(j-1);endHH(m+1,m)=b(m);H=HH;[Q1,R1]=qr(H);RR=R1(1:m,:);e1=zeros(m+1,1);e1(1,1)=1;bb1=b0*Q1'*e1;bb2=bb1(1:m,1);y=RR\(bb2);z=Q*y;t1=x-xx;x=x0+z;t2=x-xx;w1(i,1)=sqrt((t2'*(A*t2))/(t1'*(A*t1))); w2(i,1)=sqrt(t2'*(A*t2));iter=iter+1;if norm(t2)<tolbreak;endendendfunction [x,iter,w1,w2]=MINRES(A,bb,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);x=x0;for j=1:N1m=15;q=zeros(N1);a=zeros(1,m);b=zeros(1,m);r0=bb-A*x;b0=norm(r0);q(:,1)=r0/b0;a(1)=q(:,1)'*A*q(:,1);r1=A*q(:,1)-a(1)*q(:,1);b(1)=norm(r1);q(:,2)=r1/b(1);for i=2:ma(i)=q(:,i)'*A*q(:,i);r1=A*q(:,i)-a(i)*q(:,i)-b(i-1)*q(:,i-1); b(i)=norm(r1);b(i)=b(i);q(:,i+1)=r1/b(i);endQ=q(:,1:m);QQ=q(:,1:(m+1));H=QQ'*A*Q;[Q1,R1]=qr(H);RR=R1(1:m,:);e1=zeros(m+1,1);e1(1,1)=1;bb1=(b0*Q1'*e1)';bb2=bb1(1,1:m);y=RR\(bb2');z=Q*y;t1=x-xx;x=x0+z;t2=x-xx;w1(j,1)=sqrt((t2'*(A*t2))/(t1'*(A*t1)));w2(j,1)=sqrt(t2'*(A*t2));iter=iter+1;if norm(t2)<tolbreak;endx0=x;endendfunction [x,iter,w1,w2]=SOR(A,b,x0,xx)iter=0;tol=1.0e-4;w=1.5;[N1,N2]=size(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\b*w;x=x0;w1=zeros(N1,1);w2=zeros(N1,1);for i=1:N1t1=x-xx;x=B*x+f;t2=x-xx;iter=iter+1;if norm(t2)<tolbreak;endw1(i,1)=sqrt((t2'*(A*t2))/(t1'*(A*t1))); w2(i,1)=sqrt(t2'*(A*t2));endendfunction [x,iter,w1,w2]=SOR_CG(A,b,x0,xx)iter=0;%µü´ú²½Êýtol=1.0e-4;%Îó²î¿ØÖÆ[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);w=1.2;D=diag(diag(A));L=-tril(A,-1);M=(1/(w*(2-w)))*(D-w*L)*(D^(-1))*(D-w*L');x=x0;r=b-A*x;z=M\r;p=z;for i=1:N1alpha=(r'*z)/(p'*A*p);t1=x-xx;x=x+alpha*p;r1=r-alpha*A*p;z1=M\r1;belta=(z1'*z1)/(r'*z);p=z1+belta*p;r=r1;z=z1;t2=x-xx;iter=iter+1;if norm(t2)<tolbreak;endw1(i,1)=sqrt((t2'*(A*t2))/(t1'*(A*t1))); w2(i,1)=sqrt(t2'*(A*t2));endendfunction [x,iter,w1,w2]=SOR_CG(A,b,x0,xx)iter=0;tol=1.0e-4;[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);w=1.2;D=diag(diag(A));L=-tril(A,-1);M=(1/(w*(2-w)))*(D-w*L)*(D^(-1))*(D-w*L');x=x0;r=b-A*x;z=M\r;p=z;for i=1:N1alpha=(r'*z)/(p'*A*p);t1=x-xx;x=x+alpha*p;r1=r-alpha*A*p;z1=M\r1;belta=(z1'*z1)/(r'*z);p=z1+belta*p;r=r1;z=z1;t2=x-xx;iter=iter+1;if norm(t2)<tolbreak;endw1(i,1)=sqrt((t2'*(A*t2))/(t1'*(A*t1))); w2(i,1)=sqrt(t2'*(A*t2));endendclc;clear;[D,Q,A]=DMAQ;xx=ones(100,1);x0=zeros(100,1);B=zeros(100,50);KK=zeros(50,2);X_CG=zeros(100,50);Iter_CG=zeros(50,1);W1_CG=zeros(100,50);W2_CG=zeros(100,50);EE_CG=zeros(50,2);X_Lanczos1=zeros(100,50);Iter_Lanczos1=zeros(50,1);W1_Lanczos1=zeros(100,50);W2_Lanczos1=zeros(100,50);EE_Lanczos1=zeros(50,2);X_Lanczos=zeros(100,50);Iter_Lanczos=zeros(50,1);W1_Lanczos=zeros(100,50);W2_Lanczos=zeros(100,50);EE_Lanczos=zeros(50,2);X_MINRES1=zeros(100,50);Iter_MINRES1=zeros(50,1);W1_MINRES1=zeros(100,50);W2_MINRES1=zeros(100,50);EE_MINRES1=zeros(50,2);X_MINRES=zeros(100,50);Iter_MINRES=zeros(50,1);W1_MINRES=zeros(100,50);W2_MINRES=zeros(100,50);EE_MINRES=zeros(50,2);X_SOR=zeros(100,50);Iter_SOR=zeros(50,1);W1_SOR=zeros(100,50);W2_SOR=zeros(100,50);EE_SOR=zeros(50,2);X_SOR_CG=zeros(100,50);Iter_SOR_CG=zeros(50,1);W1_SOR_CG=zeros(100,50);W2_SOR_CG=zeros(100,50);EE_SOR_CG=zeros(50,2);for k=1:50K=A(:,:,k);temp=cond(K);KK(k,1)=temp;KK(k,2)=(temp-1)/(temp+1);B(:,k)=A(:,:,k)*xx;[X_CG(:,k),Iter_CG(k,1),W1_CG(:,k),W2_CG(:,k)]=CG(K,B(:,k),x0,xx); WW=W1_CG(:,k);EE_CG(k,1)=mean(WW(1:Iter_CG(k,1)));EE_CG(k,2)=norm(X_CG(:,k)-xx);[X_Lanczos(:,k),Iter_Lanczos(k,1),W1_Lanczos(:,k),W2_Lanczos(:,k)]=La nczos(K,B(:,k),x0,xx);WW=W1_Lanczos(:,k);EE_Lanczos(k,1)=mean(WW(1:Iter_Lanczos(k,1)));EE_Lanczos(k,2)=norm(X_Lanczos(:,k)-xx);[X_Lanczos1(:,k),Iter_Lanczos1(k,1),W1_Lanczos1(:,k),W2_Lanczos1(:,k) ]=Lanczos1(K,B(:,k),x0,xx);WW=W1_Lanczos1(:,k);EE_Lanczos1(k,1)=mean(WW(1:Iter_Lanczos1(k,1)));EE_Lanczos1(k,2)=norm(X_Lanczos1(:,k)-xx);[X_MINRES(:,k),Iter_MINRES(k,1),W1_MINRES(:,k),W2_MINRES(:,k)]=MINRES (K,B(:,k),x0,xx);WW=W1_MINRES(:,k);EE_MINRES(k,1)=mean(WW(1:Iter_MINRES(k,1)));EE_MINRES(k,2)=norm(X_MINRES(:,k)-xx);[X_MINRES1(:,k),Iter_MINRES1(k,1),W1_MINRES1(:,k),W2_MINRES1(:,k)]=MI NRES1(K,B(:,k),x0,xx);WW=W1_MINRES1(:,k);EE_MINRES1(k,1)=mean(WW(1:Iter_MINRES1(k,1)));EE_MINRES1(k,2)=norm(X_MINRES1(:,k)-xx);[X_SOR(:,k),Iter_SOR(k,1),W1_SOR(:,k),W2_SOR(:,k)]=SOR(K,B(:,k),x0,xx );WW=W1_SOR(:,k);EE_SOR(k,1)=mean(WW(1:Iter_SOR(k,1)));EE_SOR(k,2)=norm(X_SOR(:,k)-xx);[X_SOR_CG(:,k),Iter_SOR_CG(k,1),W1_SOR_CG(:,k),W2_SOR_CG(:,k)]=SOR_CG(K,B(:,k),x0,xx);WW=W1_SOR_CG(:,k);EE_SOR_CG(k,1)=mean(WW(1:Iter_SOR_CG(k,1)));EE_SOR_CG(k,2)=norm(X_SOR_CG(:,k)-xx);end¨subplot(211);MN=100;TT=1:MN;plot(TT,W1_CG(1:MN,1),TT,W1_CG(1:MN,11),TT,W1_CG(1:MN,21),TT,W1_CG(1: MN,31),TT,W1_CG(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{A}}/{||e_{k-1}||_{A}}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('¹²éîÌݶȷ¨');subplot(212);MN=100;TT=1:MN;plot(TT,W2_CG(1:MN,1),TT,W2_CG(1:MN,11),TT,W2_CG(1:MN,21),TT,W2_CG(1: MN,31),TT,W2_CG(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{A}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('¹²éîÌݶȷ¨');subplot(211);MN=100;TT=1:MN;plot(TT,W1_Lanczos1(1:MN,1),TT,W1_Lanczos1(1:MN,11),TT,W1_Lanczos1(1: MN,21),TT,W1_Lanczos1(1:MN,31),TT,W1_Lanczos1(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{A}}/{||e_{k-1}||_{A}}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('LanczosËã·¨');subplot(212);MN=100;TT=1:MN;plot(TT,W2_Lanczos1(1:MN,1),TT,W2_Lanczos1(1:MN,11),TT,W2_Lanczos1(1: MN,21),TT,W2_Lanczos1(1:MN,31),TT,W2_Lanczos1(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{A}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('LanczosËã·¨');subplot(211);MN=100;TT=1:MN;plot(TT,W1_Lanczos(1:MN,1),TT,W1_Lanczos(1:MN,11),TT,W1_Lanczos(1:MN, 21),TT,W1_Lanczos(1:MN,31),TT,W1_Lanczos(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{A}}/{||e_{k-1}||_{A}}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('Lanczos(m)Ëã·¨');subplot(212);MN=100;TT=1:MN;plot(TT,W2_Lanczos(1:MN,1),TT,W2_Lanczos(1:MN,11),TT,W2_Lanczos(1:MN, 21),TT,W2_Lanczos(1:MN,31),TT,W2_Lanczos(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{A}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('Lanczos(m)Ëã·¨');subplot(211);MN=100;TT=1:MN;plot(TT,W1_MINRES1(1:MN,1),TT,W1_MINRES1(1:MN,11),TT,W1_MINRES1(1:MN, 21),TT,W1_MINRES1(1:MN,31),TT,W1_MINRES1(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{A}}/{||e_{k-1}||_{A}}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('MINRESËã·¨');subplot(212);MN=100;TT=1:MN;plot(TT,W2_MINRES1(1:MN,1),TT,W2_MINRES1(1:MN,11),TT,W2_MINRES1(1:MN, 21),TT,W2_MINRES1(1:MN,31),TT,W2_MINRES1(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{A}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('MINRESËã·¨');subplot(211);MN=100;TT=1:MN;plot(TT,W1_MINRES(1:MN,1),TT,W1_MINRES(1:MN,11),TT,W1_MINRES(1:MN,21) ,TT,W1_MINRES(1:MN,31),TT,W1_MINRES(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{A}}/{||e_{k-1}||_{A}}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('MINRES(m)Ëã·¨');subplot(212);MN=100;TT=1:MN;plot(TT,W2_MINRES(1:MN,1),TT,W2_MINRES(1:MN,11),TT,W2_MINRES(1:MN,21) ,TT,W2_MINRES(1:MN,31),TT,W2_MINRES(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{A}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('MINRES(m)Ëã·¨');subplot(211);MN=100;TT=1:MN;plot(TT,W1_SOR(1:MN,1),TT,W1_SOR(1:MN,11),TT,W1_SOR(1:MN,21),TT,W1_SO R(1:MN,31),TT,W1_SOR(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{A}}/{||e_{k-1}||_{A}}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('SORËã·¨');subplot(212);MN=100;TT=1:MN;plot(TT,W2_SOR(1:MN,1),TT,W2_SOR(1:MN,11),TT,W2_SOR(1:MN,21),TT,W2_SO R(1:MN,31),TT,W2_SOR(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{A}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('SORËã·¨');subplot(211);MN=100;TT=1:MN;plot(TT,W1_SOR_CG(1:MN,1),TT,W1_SOR_CG(1:MN,11),TT,W1_SOR_CG(1:MN,21),TT,W1_SOR_CG(1:MN,31),TT,W1_SOR_CG(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{A}}/{||e_{k-1}||_{A}}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('¶Ô³ÆSORµÄÔ¤´¦Àí¹²éîÌݶȷ¨');subplot(212);MN=100;TT=1:MN;plot(TT,W2_SOR_CG(1:MN,1),TT,W2_SOR_CG(1:MN,11),TT,W2_SOR_CG(1:MN,21) ,TT,W2_SOR_CG(1:MN,31),TT,W2_SOR_CG(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{A}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('¶Ô³ÆSORµÄÔ¤´¦Àí¹²éîÌݶȷ¨');clc;clear;N=100;M11=fix(100*rand(N));[Q11,R11]=qr(M11);a61=ones(1,20);a62=-linspace(2,81,80);d6=[a61,a62];D6=diag(d6);AA=Q11*D6*Q11';xx=ones(100,1);x0=zeros(100,1);bb=AA*xx;norm_AA=zeros(2,1);norm1_AA=zeros(2,1);[X_Lanczos_AA,Iter_Lanczos_AA,W1_Lanczos_AA,W2_Lanczos_AA]=Lanczos_AA (AA,bb,x0,xx);norm_AA(1,1)=norm(X_Lanczos_AA-xx);[X_Lanczos1_AA,Iter_Lanczos1_AA,W1_Lanczos1_AA,W2_Lanczos1_AA]=Lanczo s1_AA(AA,bb,x0,xx);norm1_AA(1,1)=norm(X_Lanczos1_AA-xx);[X_MINRES_AA,Iter_MINRES_AA,W1_MINRES_AA,W2_MINRES_AA]=MINRES_AA(AA,b b,x0,xx);norm_AA(2,1)=norm(X_MINRES_AA-xx);[X_MINRES1_AA,Iter_MINRES1_AA,W1_MINRES1_AA,W2_MINRES1_AA]=MINRES1_AA (AA,bb,x0,xx);norm1_AA(2,1)=norm(X_MINRES1_AA-xx);%»-ͼMN=100;TT=1:MN;subplot(211);plot(TT,W1_Lanczos1_AA,TT,W1_MINRES1_AA)xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{2}}/{||e_{k-1}||_{2}}');legend('LanczosËã·¨','MINRESËã·¨');title('LanczosËã·¨ÓëMINRESËã·¨');subplot(212);plot(TT,W2_Lanczos1_AA,TT,W2_MINRES1_AA)xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{2}');legend('LanczosËã·¨','MINRESËã·¨');title('LanczosËã·¨ÓëMINRESËã·¨');subplot(211);plot(TT,W1_Lanczos_AA,TT,W1_MINRES_AA)xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{2}}/{||e_{k-1}||_{2}}');legend('Lanczos(m)Ëã·¨','MINRES(m)Ëã·¨');title('Lanczos(m)Ëã·¨ÓëMINRES(m)Ëã·¨');subplot(212);plot(TT,W2_Lanczos_AA,TT,W2_MINRES_AA)xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{2}');legend('Lanczos(m)Ëã·¨','MINRES(m)Ëã·¨');title('Lanczos(m)Ëã·¨ÓëMINRES(m)Ëã·¨');function [x,iter,w1,w2]=Lanczos1_AA(A,bb,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);m=N1;w1=zeros(N1,1);w2=zeros(N1,1);x=x0;q=zeros(N1);a=zeros(1,m);b=zeros(1,m);r0=bb-A*x0;b0=norm(r0);q(:,1)=r0/b0;a(1)=q(:,1)'*A*q(:,1);r1=A*q(:,1)-a(1)*q(:,1);b(1)=norm(r1);q(:,2)=r1/b(1);for i=2:N1a(i)=q(:,i)'*A*q(:,i);r1=A*q(:,i)-a(i)*q(:,i)-b(i-1)*q(:,i-1);b(i)=norm(r1);b(i)=b(i);q(:,i+1)=r1/b(i);m=i;Q=q(:,1:m);H=Q'*A*Q;y=H\(Q'*r0);z=Q*y;t1=x-xx;x=x0+z;t2=x-xx;w1(i,1)=norm(t2)/norm(t1);w2(i,1)=norm(t2);iter=iter+1;if norm(t2)<tolbreak;endendendfunction [x,iter,w1,w2]=Lanczos_AA(A,bb,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);x=x0;for j=1:N1m=15;q=zeros(N1);a=zeros(1,m);b=zeros(1,m);r0=bb-A*x;b0=norm(r0);q(:,1)=r0/b0;a(1)=q(:,1)'*A*q(:,1);r1=A*q(:,1)-a(1)*q(:,1);b(1)=norm(r1);q(:,2)=r1/b(1);for i=2:ma(i)=q(:,i)'*A*q(:,i);r1=A*q(:,i)-a(i)*q(:,i)-b(i-1)*q(:,i-1);b(i)=norm(r1);b(i)=b(i);q(:,i+1)=r1/b(i);endQ=q(:,1:m);H=Q'*A*Q;y=H\(Q'*r0);z=Q*y;t1=x-xx;x=x0+z;t2=x-xx;w1(j,1)=norm(t2)/norm(t1);w2(j,1)=norm(t2);iter=iter+1;if norm(t2)<tolbreak;endx0=x;endendfunction [x,iter,w1,w2]=MINRES1_AA(A,bb,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);x=x0;m=N1;q=zeros(N1);a=zeros(1,m);b=zeros(1,m);r0=bb-A*x0;b0=norm(r0);q(:,1)=r0/b0;a(1)=q(:,1)'*A*q(:,1);r1=A*q(:,1)-a(1)*q(:,1);b(1)=norm(r1);q(:,2)=r1/b(1);for i=2:N1a(i)=q(:,i)'*A*q(:,i);r1=A*q(:,i)-a(i)*q(:,i)-b(i-1)*q(:,i-1);b(i)=norm(r1);b(i)=b(i);q(:,i+1)=r1/b(i);m=i;Q=q(:,1:m);HH=zeros(m+1,m);for j=1:mHH(j,j)=a(j);endfor j=1:m-1HH(j,j+1)=b(j);endfor j=2:mHH(j,j-1)=b(j-1);endHH(m+1,m)=b(m);H=HH;[Q1,R1]=qr(H);RR=R1(1:m,:);e1=zeros(m+1,1);e1(1,1)=1;bb1=b0*Q1'*e1;bb2=bb1(1:m,1);y=RR\(bb2);z=Q*y;t1=x-xx;x=x0+z;t2=x-xx;w1(i,1)=norm(t2)/norm(t1);w2(i,1)=norm(t2);iter=iter+1;if norm(t2)<tolbreak;endendendfunction [x,iter,w1,w2]=MINRES_AA(A,bb,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);x=x0;for j=1:N1m=15;q=zeros(N1);a=zeros(1,m);b=zeros(1,m);r0=bb-A*x;b0=norm(r0);q(:,1)=r0/b0;a(1)=q(:,1)'*A*q(:,1);r1=A*q(:,1)-a(1)*q(:,1);b(1)=norm(r1);q(:,2)=r1/b(1);for i=2:ma(i)=q(:,i)'*A*q(:,i);r1=A*q(:,i)-a(i)*q(:,i)-b(i-1)*q(:,i-1); b(i)=norm(r1);b(i)=b(i);q(:,i+1)=r1/b(i);endQ=q(:,1:m);QQ=q(:,1:(m+1));H=QQ'*A*Q;[Q1,R1]=qr(H);RR=R1(1:m,:);e1=zeros(m+1,1);e1(1,1)=1;bb1=(b0*Q1'*e1)';bb2=bb1(1,1:m);y=RR\(bb2');z=Q*y;t1=x-xx;x=x0+z;t2=x-xx;w1(j,1)=norm(t2)/norm(t1);w2(j,1)=norm(t2);iter=iter+1;if norm(t2)<tolbreak;endx0=x;endend。