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

病态：高等数值分析实验作业一高等数值分析实验作业一1.构造例子说明CG的数值性态。当步数=阶数时CG的解如何？当A的最大特征值远大于第二个最大特征值,最小特征值远小于第二个最小特征值时方法的收敛性如何?解：（1）构造1002阶正定矩阵A：A的特征值向量为，分为良态和病态情形;（2）CG法求解Ax=b，A良态:利用matlab编程实现CG算法。b=ones(1002,1)，x0=zeros(1002,1)。计算每一步迭代的残差rk相对于初始残差的2范数。相对残差2范数的对数值与迭代步数的关系曲线如图1所示：图1CG法求解良态问题的收敛曲线（3）CG法求解Ax=b，A利用matlab编程实现CG算法。b=ones(1002,1)，x0=zeros(1002,1)。计算每一步迭代的残差rk相对于初始残差的2范数。相对残差2范数的对数值与迭代步数的关系曲线如图2所示：02040608010012014016018010-810-710-610-510-410-310-210-1100CG法收敛曲线(阶数n=1002)迭代次数||rk||/||b||良态：高等数值分析实验作业一图2CG法求解病态问题的收敛曲线结论：对于良态问题，CG法收敛较快，相对残差的2范数随迭代步数的增加而减小；对于病态问题CG法收敛的很慢，而且收敛曲线不平滑，震荡非常严重，相对残差的2范数没有最优性。2.对于同样的例子，比较CG和Lanczos的计算结果。解：（1）构造1002阶正定矩阵A：（2）Lanczos法求解Ax=b，A利用matlab编程实现CG算法。b=ones(1002,1)，x0=zeros(1002,1)。计算每一步迭代的残差rk相对于初始残差的2范数。相对残差2范数的对数值与迭代步数的关系曲线如图3所示：05010015020025030035040010-810-610-410-2100102104CG法收敛曲线(阶数n=1002)迭代次数||rk||/||b||图3Lanczos法求解良态问题的收敛曲线02040608010012014016018010-910-810-710-610-510-410-310-210-1100Lanzcos算法的收敛曲线(阶数n=1002)迭代次数||rk||/||b||病态：高等数值分析实验作业一（3）Lanczos法求解Ax=b，A利用matlab编程实现CG算法。b=ones(1002,1)，x0=zeros(1002,1)。计算每一步迭代的残差rk相对于初始残差的2范数。相对残差2范数的对数值与迭代步数的关系曲线如图4所示：图4Lanczos法求解病态问题的收敛曲线结论：对于良态正定问题，Lanczos法和CG法收敛性一样，相对残差的2范数随迭代步数的增加而减小；对于病态正定问题Lanczos法较CG法收敛的要慢一点，而且收敛曲线更加不平滑，震荡较严重，相对残差的2范数同样没有最优性。05010015020025030035040010-810-610-410-2100102104Lanzcos算法的收敛曲线(阶数n=1002)迭代次数||rk||/||b||事实上，Lanczos法主要是解对称不定问题，针对这类问题才显现出它的优越性，这也从侧面表明，对称正定问题不适合用Lanzcos法求解。3.当A只有m个不同特征值时，对于大的m和小的m，观察有限精度下Lanczos方法如何收敛。解：分别构建m=10、50、100、400、800五个矩阵A，分别求解Ax=b，收敛曲线如图5所示：图5对不同的m，Lanczos法求解Ax=b的收敛曲线02040608010012014016010-1510-1010-5100Lanzcos算法的收敛曲线(阶数n=1002)迭代次数||rk||/||b||m=10m=50m=100m=400m=800高等数值分析实验作业一结论：如果A只有m个不同的特征值，则Lanczos方法至多m步就可以找到精确解。实验中，在m较大的时候，算法收敛较快，迭代次数远小于m。当m较小时，可能需要接近于m步才能找到准确解。4.取初始值近似解为零向量，右端项b仅由A的m个不同特征向量的线性组合表示时，Lanczos方法的收敛性如何？数值计算中方法的收敛性和m的大小关系如何？解：设A的特征值向量为：1=[1002:-1:1]，分别构建m=10、50、100、400、800时的b，即b分别由10、50、100、400、800个A的特征向量线性表示，分别求解Ax=b，收敛曲线如图6所示：图6对于b由m个特征向量线性表示时，Lanczos法求解Ax=b的收敛曲线0246810121416182010-1210-1010-810-610-410-2100Lanzcos算法的收敛曲线(阶数n=1002)迭代次数||rk||/||b||m=10m=50m=100m=400m=800结论：理论上，b仅由A的m个不同特征向量的线性组合表示时，Lanczos方法必然m步收敛。但由于A的阶数为1002，是比较大的，精度方面的限制导致计算得到的m个特征向量并不都线性无关，所以，m=800时只需20步迭代。高等数值分析实验作业一5.构造对称不定矩阵，验证Lanczos方法的近似中断，观察收敛曲线中的峰点个数和特征值的分布关系；观测当出现峰点时，MINRES方法的收敛形态怎样。解：分别构建负特征值个数为m=10、50、100、200、500的矩阵A，分别计算Ax=b的解，Lanzcos方法的收敛曲线如图7-图11所示，Minres的收敛曲线如图12-图16所示：图7m=10时，Lanczos法求解Ax=b的收敛曲线05010015020025030035010-1010-810-610-410-2100102Lanzcos算法的收敛曲线(阶数n=1002）迭代次数||rk||/||b||图8m=50时，Lanczos法求解Ax=b的收敛曲线010020030040050060070010-1010-810-610-410-2100102104Lanzcos算法的收敛曲线(阶数n=1002）迭代次数||rk||/||b||图9m=100时，Lanczos法求解Ax=b的收敛曲线010020030040050060070080090010-1010-810-610-410-2100102Lanzcos算法的收敛曲线(阶数n=1002）迭代次数||rk||/||b||高等数值分析实验作业一图10m=200时，Lanczos法求解Ax=b的收敛曲线02004006008001000120010-1010-810-610-410-2100102104Lanzcos算法的收敛曲线(阶数n=1002）迭代次数||rk||/||b||图11m=500时，Lanczos法求解Ax=b的收敛曲线结论：负特征值越多，Lanczos方法的收敛曲线的峰点个数越多，纹波越大，发生近似中断的次数越来越多。05001000150010-1010-810-610-410-2100102Lanzcos算法的收敛曲线(阶数n=1002）迭代次数||rk||/||b||高等数值分析实验作业一图12m=10时，Minres法求解Ax=b的收敛曲线05010015020025030035010-1010-810-610-410-2100102Minres算法的收敛曲线(阶数n=1002）迭代次数||rk||/||b||图13m=50时，Minres法求解Ax=b的收敛曲线图14m=100时，Minres法求解Ax=b的收敛曲线010020030040050060070010-1010-810-610-410-2100102Minres算法的收敛曲线(阶数n=1002）迭代次数||rk||/||b||010020030040050060070080090010-1010-810-610-410-2100102Minres算法的收敛曲线(阶数n=1002）迭代次数||rk||/||b||高等数值分析实验作业一图15m=200时，Minres法求解Ax=b的收敛曲线02004006008001000120010-310-210-1100101102Minres算法的收敛曲线(阶数n=1002）迭代次数||rk||/||b||图16m=500时，Minres法求解Ax=b的收敛曲线结论：对于Lanczos方法，随着负特征值的增多，收敛曲线的峰点个数增多，振荡越来越严重，发生近似中断的次数越来越多；然而，对于相同的A，Minres方法的相对残差没有出现峰值，随着迭代数增加而单调下降。收敛性态比Lanczos好，但代价是计算时间的急剧增大。02004006008001000120010-310-210-1100101102Minres算法的收敛曲线(阶数n=1002）迭代次数||rk||/||b||高等数值分析实验作业一附件：主要算法代码CG法CG.mfunction[x,Error,i,flag]=CG(A,b,x,ErrSet,uplimit)[m,n]=size(b);ifmnb=b';end[m,n]=size(x);ifmnx=x';endr=b-A*x;p=r;i=1;temp_rkrkplus=r'*r;Error=sqrt(temp_rkrkplus)/norm(b,2);while1temp_AP=A*p;temp_rkrk=temp_rkrkplus;temp_pAP=p'*temp_AP;ifabs(temp_pAP)1e-12disp('¶ñÐÔÖжϣ¡')break;enda=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;Err=sqrt(temp_rkrkplus)/norm(b,2);%/(norm(b)+temp_AP);ifErrErrSetdisp('Methodconverge£¡')disp(i)flag=1;break;endError=[Error,sqrt(temp_rkrkplus)/norm(b,2)];ifiuplimitflag=0;breakendendend高等数值分析实验作业一endLanczos法lanczos.mfunction[T,Q,x,k,Errs,tiao]=Lanczos(A,b,x0,Err)x=[];tiao=[];[m,n]=size(b);ifmnb=b';end[m,n]=size(x0);ifmnm=n;x0=x0';endsize(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)];while1disp('k=');disp(k);T=[T,zeros(k,1);zeros(1,k),0];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);ifk==1disp('Ìø¹ý´Ë²½')tiao=[tiao,k];else[L,U]=LanczosLU(Tk);bk=norm(r_zeros,2)*[1;zeros(k-1,1)];高等数值分析实验作业一Ux=zeros(k,1);Ux(1)=bk(1);fori=2:kUx(i)=bk(i)-L(i,i-1)*Ux(i-1);endym=zeros(k,1);ym(k)=Ux(k)/U(k,k);fori=k-1:-1:1ym(i)=(Ux(i)-ym(i+1)*U(i,i+1))/U(i,i);endErrs=[Err