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

迭代次数

l g （|r k |）

高等数值分析第一次实验

T1. 构造例子说明CG 的数值形态。当步数 = 阶数时CG 的解如何？当A 的最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值时，方法的收敛性如何？ Answer:

对于问题1：当步数 = 阶数时CG 的解如何？ ➢ 在MA TLAB 中构造N 阶对称正定矩阵代码如下：

N=1000

D = 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)

计算结果如下（绘制出来的log 由上表可以看出对于对称正定矩阵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;

计算结果如右图所示，首先对比可以发现矩阵的收敛速度跟其条件数大小有关，条件数小时，收敛速度快，算例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*U;

迭代次数

l g （|r k |）

在计算时，取X0=zeros(N,1);b=ones(N,1);

CG算法代码同上，在计算时取停机准则为绝对误差e<10-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);

end

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=[q q1];

while norm(rk)>e & k

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，LL T y=‖r0‖e1的代码如下：

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);

end

迭代次数

l

g

（

|

r

k

|

）