高等数值分析作业-第一次实验

高等数值分析作业-第一次实验（总

8页）本页仅作为预览文档封面，使用时请删除本页-#YKN函数：n=size(TK,1);yk=zeros(n,1);yk(n)=b(n)/TK(n,n);fori=1:n-1j=n-i；yk(j)=(b(j)-TK(j,j+1)*yk(j+1))/TK(j,j)end首先构造1000阶有m个特征值的矩阵，1000YKN函数：n=size(TK,1);yk=zeros(n,1);yk(n)=b(n)/TK(n,n);fori=1:n-1j=n-i；yk(j)=(b(j)-TK(j,j+1)*yk(j+1))/TK(j,j)end首先构造1000阶有m个特征值的矩阵，1000)，构造矩阵的代码如下：N=1000m=10;DIA=linspace(1,1000,m);VEC=zeros(N,1);k=N/m;i=1；ii=0;whilei<=mforj=1:kVEC(ii+j)=DIA(i);endii=ii+k;i=i+1;构造是保证矩阵的条件数相同（这里均取为kg 10个特征值 20个特征值 50个特征值 100个特征值 500个特征值1000个特征值-12100 150 200迭代次数2502010100908070605040 10个特征值 20个特征值 50个特征值 100个特征值 500个特征值1000个特征值00 50100 150迭代次数200 250lyk=YK(L',norm(r0)*eye(k,1));yk=YKN(L,lyk);YK函数：n=size(TK,1);yk=zeros(n,1);yk(1)=b(1)/TK(1,1);fori=2:nyk(i)=(b(i)-TK(i,i-1)*yk(i-1))/TK(i,i);end计算结果如上图所示，计算时间与迭代步数如下:计算方法CG方法Lanczos方法迭代步数234235运行时间（s）可以得到如下结论：当矩阵A为对称正定时，两种方法效果相当，每一步误差也基本相同，但收敛速度基本一样。由于Lanczos方法的每一步迭代中都有一个Lanczos过程，其中需要构造Tk和Qk，以及计算yk，故该方法需要耗费更长的计算时间。同时计算过程中发现，当取得收敛准则比较严格时，CG算法较Lanczos方法稳定。T3.当A只有m个不同特征值时，对于大的m和小的m，观察有限精度下Lanczos方法如何收敛。Answer：

endD=diag(VEC);U=orth(rand(N,N));A1=U'*D*U;对构造好的A1进行计算，在计算时，取X0=zeros(N,1);b=ones(N,1);在计算时取停机准则为绝对误差e<10-i0,Lanczos方法的代码同上，得出的结果如上图所示(分别绘制出log10ll*ll和llrMI随迭代次数的变化图)，各种情况下的迭代次数及计算时间如下表所示(这里取得矩阵为对称正定矩阵，因此Lanczos过程求解yk时采用cholesky分解)：特征值个数1020501005001000运行时间迭代次数10204973159215可以看到，几个问题的条件数虽然相同，如果A只有m个不同的特征值，则Lanczos方法至多m步就可以找到精确解。实验中，在m较大的时候，算法收敛较快，远小于m。当m较小时，可能需要接近于m步才能找到准确解。另外，由上面计算可以看到在m值较小时，其第k步的误差llrkll会有一个在计算时会有先增大后减小的过程，最后收敛。T4.取初始值近似解为零向量，右端项b仅由A的m个不同特征向量的线性组合表示时，Lanczos方法的收敛性如何？数值计算中方法的收敛性和m的大小关系如何？Answer：在前面的基础上，我们知道之前构造对称正定矩阵时，U的每一个行向量均为对应A的特征值(D的某一个元素)的特征向量，因此构造m不同时的又端项b代码如下所示：N=1000m=10;k|gg-6-8-10k|gg-6-8-100 50 100 150 200 250迭代次数-12D=diag(DIA);U=orth(rand(N,N));A1=U'*D*U;b1=zeros(N,1);fori=1:mb1=b1+rand(1)*(U(i,:))';endX0=zeros(N,1);e=10八(-10);tic[S1,ek,num]=LanczosT3(A1,X0,b1,e);toc计算时采用的Lanczos算法同之前的题目，计算出来的结果如上所示(这里取得矩阵为对称正定矩阵，因此Lanczos过程求解yk时采用cholesky分解)。各种情况下的迭代次数及计算时间如下表所示：m1020501005007501000运行时间

迭代次数68120165179202205212我们可以得到如下结论：M较小时，计算过程中收敛误差会产生波动，m较大时，波动减小，最后均收敛；收敛步数随着m值的增大而增加，但是由上表我们可以看出收敛步数的增幅逐渐减小，可以推测当m达到一定值后，收敛步数可能趋近于某一定值。T5.构造对称不定矩阵，验证Lanczos方法的近似中断，观察收敛曲线中的峰点个数和特征值的分布关系；观测当出现峰点时，MINRES方法的收敛形态怎样。Answer：对于对称不定矩阵，求解yk时采用追赶法，追赶法的代码如下所示:functions=Zhuigan(A,bbfunctions=Zhuigan(A,bb)n=size(A,1);s=zeros(n,1);b=diag(A);a=diag(A,-1);c=diag(A,1);d=zeros(n,1);u=zeros(n-1,1);fori=1:n-1d(1)=b(1);u(i)=c(i)/d(i);d(i+1)=b(i+1)-a(i)*u(i);end而对于MINRES方法，% 追的过程 y=zeros(n,1);y(1)=bb(1)/d(1);fori=2:ny(i)=(bb(i)-a(i-1)*y(i-1))/d(i);end% 赶的过程 s(n)=y(n);fori=n-1:-1:1s(i)=y(i)-u(i)*s(i+1);Lanczos算法不同之处在于，ek=zeros(k,1);ek(k)=1;endend考虑求解yk时用MATLAB自带的qr分解进行计算，其中与求解yk部分代码如下：TKT=[TKbj(k)*ek]';[Q,R]=qr(TKT);gk=Q'*norm(r0)*eye(k+1,1);yk=minresYK(R,gk);其中，：^口侬丫长为自己编写的有R-gk=yk求解yk的函数，代码如下:functionyk=minresYK(TK,b)n=size(TK,2);yk=zeros(n,1);yk(n)=b(n)/TK(n,n);fori=1:n-1k=n-i;yk(k)=b(k);forj=k+1:nyk(k)=yk(k)-TK(k,j)*yk(j);endyk(k)=yk(k)/TK(k,k);

endend对于有m个特征值的对称不定矩阵，本文中采用的生成代码如下(基本思想同之前相同，生成矩阵的特征值分别为-m,-(m-1),……,-1,1,2,……N-m，其中N为矩阵阶数，计算过程中按1000考虑)：N=1000%10个负特征值m=10;DIA=linspace(-m+1,N-m,N);DIA(m)=-m;D=diag(DIA);U=orth(rand(N,N));A1=U'*D*U;取X0=zeros(N,1);b=ones(N,1);在计算时取停机准则为绝对误差e<10-i0，计算结果420)|k|r(-2-4lg-6-8420)|k|r(-2-4lg-6-8-100 50 100150迭代20次0数250 300 35010-12400负特征值个数：Lanczos步数：Lanczos步数：212计算时间：MINRES步数：210计算时间：Lanczos步数：369计算时间：MINRES步数：367计算时间：420-4-6-8-10|r(lg) -20 100 200 300 400 500 600 700迭代次数420-4-6-8-10|r(lg) -20 100 200 300 400 500 600 700迭代次数-12负特征值个数：50Lanczos步数：697计算时间：MINRES步数：689计算时间：Lanczos步数：944计算时间：MINRES步数：939计算时间：负特征值个数：200Lanczos未收敛计算时间：MINRES未收敛计算时间：Lanczos未收敛计算时间：负特征值个数：200Lanczos未收敛计算时间：MINRES未收敛计算时间：Lanczos方法有峰点说明，Lanczos方法可能会近似中断；负特征值越多，峰点个数越多。实际上，更一般的结论是，当正