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

《高等数值分析作业-第一次实验》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、高等数值分析第一次实验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=

2、rk'*rk;fork=1:Nk=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;endplot(r,m,'r-');Ek=m(k)计算结果如下（绘制出来的log10rk随迭代次数的变化如上图所示）：N1000200030004000log10rk-81.8505-98.3653-126.3256-115.8889运

3、行时间（s）4.30985530.205448105.792648289.610550由上表可以看出对于对称正定矩阵A，CG算法还是比较稳定的，但求解步数=阶数时，CG算法的解即为准确解（误差极小）。第

4、页8对于问题2：当A的最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值时，方法的收敛性如何？Ø构造1000阶的对称正定矩阵如下，收敛准则取为绝对ε<10(-10)：首先构造一个特征值分别为0.1到1的对称正定矩阵A，代码如下（算例1）：D=diag(linspace(0.1,1,N));U=orth(r

5、and(N,N));A=U'*D*U;在之前的基础上，将最大特征值调为107，最小特征值为10-5，代码如下（算例2）：DIA=linspace(0.1,1,N);DIA(1)=10^(-5);DIA(N)=10^7;D=diag(DIA);U=orth(rand(N,N));A=U'*D*U;最后生成一个特征值在10-5到107均匀分布的矩阵（算例3）:DIA=linspace(10^(-5),10^7,N);D=diag(DIA);U=orth(rand(N,N));A=U'*D*U;计算结果如右图所示，首先对比可以

6、发现矩阵的收敛速度跟其条件数大小有关，条件数小时，收敛速度快，算例1>2>3，同时，A的中间特征值分布对CG的收敛速度有巨大的影响。实际上，在经过几步后，CG的收敛因子将是：λ2λn-1-1λ2λn-1+1而非λ1λn-1λ1λn+1因此，本题中算例2的矩阵的收敛速度较算例3快很多，而与算例1较为接近。T2.对于同样的例子，比较CG和Lanczos的计算结果。Answer:首先构造一个1000阶的对称正定矩阵，代码如下：D=diag(linspace(1,1000,N));U=orth(rand(N,N));A=U'*D

7、*U;第

8、页8在计算时，取X0=zeros(N,1);b=ones(N,1);CG算法代码同上，在计算时取停机准则为绝对误差e<10-12，Lanczos算法的代码如下所示（这里取得矩阵为对称正定矩阵，因此Lanczos过程求解yk时采用cholesky分解）：第

9、页8Xk=X0;n=size(A,1);aj=zeros(n,1);bj=zeros(n,1);r0=b-A*Xk;rk=r0;q=rk/norm(rk);r=A*q;aj(1)=q'*r;r=r-aj(1)*q;ifr~=0bj(1)=norm(r);end

10、q1=r/bj(1);qk=q;k=1;TK(1,1)=aj(1);L=chol(TK);lyk=YK(L',norm(r0)*eye(k,1));yk=YKN(L,lyk);Xk=X0+qk*yk;rk=b-A*Xk;m(k)=log10(norm(rk));num(k)=k;TK(1,2)=bj(1);TK(2,1)=bj(1);qk=[qq1];whilenorm(rk)>e&k

11、;endq=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);T