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

高等数值分析第一次实验T1. 构造例子说明CG的数值形态。当步数 = 阶数时CG的解如何？当A的最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值时，方法的收敛性如何？Answer:对于问题1：当步数 = 阶数时CG的解如何？ 在MATLAB中构造N阶对称正定矩阵代码如下：N=1000D = diag(rand(N,1);U = orth(rand(N,N);A = U*D*U;在计算时，取X0=zeros(N,1);b=ones(N,1);自己编写CG算法，如下：Xk = X0;rk=b-A*Xk;pk=rk;crk_1=rk*rk;for k=1:N k=k+1; apk=A*pk; ak=crk_1/(pk*apk); Xk=Xk+ak*pk; rk=rk-ak*apk; crk=rk*rk; bk_1=crk/crk_1; crk_1=crk; pk=rk+bk_1*pk; m(k)=norm(rk); r(k)=k;end plot(r,m,r-); Ek=m(k)计算结果如下（绘制出来的log10rk随迭代次数的变化如上图所示）：N1000200030004000log10rk-81.8505-98.3653-126.3256-115.8889运行时间（s）4.30985530.205448105.792648289.610550由上表可以看出对于对称正定矩阵A，CG算法还是比较稳定的，但求解步数=阶数时，CG算法的解即为准确解（误差极小）。对于问题2：当A的最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值时，方法的收敛性如何？ 构造1000阶的对称正定矩阵如下，收敛准则取为绝对23，同时，A的中间特征值分布对CG的收敛速度有巨大的影响。实际上，在经过几步后，CG的收敛因子将是：2n-1-12n-1+1而非1n-11n+1因此，本题中算例2的矩阵的收敛速度较算例3快很多，而与算例1较为接近。T2. 对于同样的例子，比较CG和Lanczos的计算结果。Answer:首先构造一个1000阶的对称正定矩阵，代码如下：D = diag(linspace(1,1000,N);U = orth(rand(N,N);A = U*D*U;在计算时，取X0=zeros(N,1);b=ones(N,1);CG算法代码同上，在计算时取停机准则为绝对误差ee & kn k=k+1; r=A*q1-bj(k-1)*q; aj(k)=q1*r; r=r-aj(k)*q1; if r=0 bj(k)=norm(r); end q=q1; q1=r/bj(k); TK(k,k)=aj(k); L=chol(TK); lyk=YK(L,norm(r0)*eye(k+1,1); yk=YKN(L,lyk); Xk=X0+qk*yk; rk=b-A*Xk; m(k)=log10(norm(rk); num(k)=k; TK(k,k+1)=bj(k); TK(k+1,k)=bj(k); qk=qk q1;end在计算时直接采用cholesky分解（直接采用matlab自带的函数）计算yk，LLTy=r0e1的代码如下： L=chol(TK); lyk=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);for i=2:n yk(i)=(b(i)-TK(i,i-1)*yk(i-1)/TK(i,i);endYKN函数：n=size(TK,1);yk=zeros(n,1);yk(n)=b(n)/TK(n,n);for i=1:n-1 j=n-i; yk(j)=(b(j)-TK(j,j+1)*yk(j+1)/TK(j,j);end计算结果如上图所示，计算时间与迭代步数如下：计算方法CG方法Lanczos方法迭代步数234235运行时间（s）0.1870810.513452可以得到如下结论： 当矩阵A为对称正定时，两种方法效果相当，每一步误差也基本相同，但收敛速度基本一样。 由于Lanczos方法的每一步迭代中都有一个Lanczos过程，其中需要构造Tk和Qk，以及计算yk，故该方法需要耗费更长的计算时间。 同时计算过程中发现，当取得收敛准则比较严格时，CG算法较Lanczos方法稳定。T3. 当A只有m个不同特征值时，对于大的m和小的m，观察有限精度下Lanczos方法如何收敛。Answer：首先构造1000阶有m个特征值的矩阵，构造是保证矩阵的条件数相同（这里均取为1000），构造矩阵的代码如下：N=1000m=10; DIA=linspace(1,1000,m); VEC=zeros(N,1);k=N/m;i=1;ii=0;while i=m for j=1:k VEC(ii+j)=DIA(i); end ii=ii+k; i=i+1;endD=diag(VEC);U = orth(rand(N,N); A1 = U*D*U;对构造好的A1进行计算，在计算时，取X0=zeros(N,1);b=ones(N,1); 在计算时取停机准则为绝对误差e10-10，Lanczos方法的代码同上，得出的结果如上图所示（分别绘制出log10rk和rk随迭代次数的变化图），各种情况下的迭代次数及计算时间如下表所示（这里取得矩阵为对称正定矩阵，因此Lanczos过程求解yk时采用cholesky分解）：特征值个数1020501005001000运行时间0.0157270.0285230.0722000.1075780.2839860.453770迭代次数10204973159215可以看到，几个问题的条件数虽然相同，如果A只有m个不同的特征值，则Lanczos方法至多m步就可以找到精确解。实验中，在m较大的时候，算法收敛较快，远小于m。当m较小时，可能需要接近于m步才能找到准确解。另外，由上面计算可以看到在m值较小时，其第k步的误差rk会有一个在计算时会有先增大后减小的过程，最后收敛。T4. 取初始值近似解为零向量，右端项b仅由A的m个不同特征向量的线性组合表示时，Lanczos方法的收敛性如何？数值计算中方法的收敛性和m的大小关系如何？Answer：在前面的基础上，我们知道之前构造对称正定矩阵时，U的每一个行向量均为对应A的特征值（D的某一个元素）的特征向量，因此构造m不同时的又端项b代码如下所示：N=1000m=10;DIA=linspace(1,N,N); D=diag(DIA);U = orth(rand(N,N);A1 = U*D*U;b1=zeros(N,1);for i=1:m b1=b1+rand(1)*(U(i,:);endX0=zeros(N,1);e=10(-10);ticS1,ek,num=LanczosT3(A1,X0,b1,e);toc计算时采用的Lanczos算法同之前的题目，计算出来的结果如上所示（这里取得矩阵为对称正定矩阵，因此Lanczos过程求解yk时采用cholesky分解）。各种情况下的迭代次数及计算时间如下表所示：m1020501005007501000运行时间0.10740.19500.29670.35460.35460.43250.4275迭代次数68120165179202205212我们可以得到如下结论： M较小时，计算过程中收敛误差会产生波动，m较大时，波动减小，最后均收敛； 收敛步数随着 m 值的增大而增加，但是由上表我们可以看出收敛步数的增幅逐渐减小，可以推测当 m 达到一定值后，收敛步数可能趋近于某一定值。T5. 构造对称不定矩阵，验证Lanczos方法的近似中断，观察收敛曲线中的峰点个数和特征值的分布关系；观测当出现峰点时，MINRES方法的收敛形态怎样。Answer：对于对称不定矩阵，求解yk时采用追赶法，追赶法的代码如下所示：function s = 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);for i=1:n-1 d(1)=b(1); u(i)=c(i)/d(i); d(i+1)=b(i+1)-a(i)*u(i);end%-追的过程-y=zeros(n,1);y(1)=bb(1)/d(1);for i=2:n y(i)=(bb(i)-a(i-1)*y(i-1)/d(i);end%-赶的过程-s(n)=y(n);for i=n-1:-1:1 s(i)=y(i)-u(i)*s(i+1);endend而对于MINRES方法，考虑求解yk时用MATLAB自带的qr分解进行计算，其中与Lanczos算法不同之处在于，求解yk部分代码如下： ek=zeros(k,1); ek(k)=1; TKT=TK bj(k)*ek; Q,R=qr(TKT); gk=Q*norm(r0)*eye(k+1,1); yk=minresYK(R,gk);其中，minresYK为自己编写的有Rgk=yk求解yk的函数，代码如下：function yk = minresYK(TK,b)n=size(TK,2);yk=zeros(n,1);yk(n)=b(n)/TK(n,n);for i=1:n-1 k=n-i; yk(k)=b(k); for j=k+1:n yk(k)=yk(k)-TK(k,j)*yk(j); end yk(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); 在计算时取停机准则为绝对误差e10-10，计算结果如下：负特征值个数：0Lanczos 步数：212 计算时间：0.368697MINRES 步数：210 计算时间：0.544282负特征值个数：10Lanczos 步数：369 计算时间：0.841989MINRES 步数：367 计算时间：1.776133负特征值个数：50Lanczos 步数：697 计算时间：2.673012MINRES 步数：689 计算时间：10.111938负特征值个数：100Lanczos 步数：944 计算时间：5.155608MINRES 步数：939 计算时间：24.470218负特征值个数：150Lanczos 未收敛 计算时间：5.897715MINRES 未收敛 计算时间：28.467015负特征值个数：200Lanczos 未收敛 计算时间：5.897715MINRES 未收敛 计算时间：28.467015由上述结果得到的结论如下： Lanczos方法有峰点说明