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

高等数值分析第一次实验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)计算结果如下（绘制出来的随迭代次数的变化如上图所示）：N1000200030004000-81.8505-98.3653-126.3256-115.8889运行时间（s）4205448105.792648289.610550由上表可以看出对于对称正定矩阵A，CG算法还是比较稳定的，但求解步数=阶数时，CG算法的解即为准确解（误差极小）。对于问题2：当A的最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值时，方法的收敛性如何？构造1000阶的对称正定矩阵如下，收敛准则取为绝对ε10(-10)：首先构造一个特征值分别为0.1到1的对称正定矩阵A，代码如下（算例1）：D = diag(linspace(0.1,1,N));U = orth(rand(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;计算结果如右图所示，首先对比可以发现矩阵的收敛速度跟其条件数大小有关，条件数小时，收敛速度快，算例123，同时，A的中间特征值分布对CG的收敛速度有巨大的影响。实际上，在经过几步后，CG的收敛因子将是：而非因此，本题中算例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算法代码同上，在计算时取停机准则为绝对误差e10-12，Lanczos算法的代码如下所示（这里取得矩阵为对称正定矩阵，因此Lanczos过程求解yk时采用cholesky分解）：Xk = 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;if r~=0 bj(1)=norm(r);endq1=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=[q q1];while norm(rk)e 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(