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

高等数值分析第一次实验
T1. 构造例子说明CG的数值形态。当步数 = 阶数时CG的解如何？当A的最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值时，方法的收敛性如何？
Answer:
对于问题1：当步数 = 阶数时CG的解如何？
Ø 在MATLAB中构造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)
计算结果如下（绘制出来的log10rk随迭代次数的变化如上图所示）：
N
1000
2000
3000
4000
log10rk
-81.8505
-98.3653
-126.3256
-115.8889
运行时间（s）
4.309855
30.205448
105.792648
289.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;
计算结果如右图所示，首先对比可以发现矩阵的收敛速度跟其条件数大小有关，条件数小时，收敛速度快，算例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;
在计算时，取X0=zeros(N,1);b=ones(N,1);
CG算法代码同上，在计算时取停机准则为绝对误差e<10-12，Lanczos算法的代码如下所示（这里取得矩阵为对称正定矩阵，因此Lanczos过程求解yk时采用cholesky分解）：
第 | 页 8

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 & kk+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);
end

YKN函数：
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方法
迭代步数
234
235
运行时间（s）
0.187081
0.513452
可以得到如下结论：
Ø 当矩阵A为对称正定时，两种方法效果相当，每一步误差也基本相同，但收敛速度基本一样。
Ø 由于Lanczos方法的每一步迭代中都有一个Lanczos过程，其中需要构造Tk和Qk，以及计算yk，故该方法需要耗费更长的计算时间。
Ø 同时计算过程中发现，当取得收敛准则比较严格时，CG算法较Lanczos方法稳定。
T3. 当A只有m个不同特征值时，对于大的m和小的m，观察有限精度下Lanczos方法如何收敛。
Answer：
首先构造1000阶有m个特征值的矩阵，构造是保证矩阵的条件数相同（这里均取为1000），构造矩阵的代码如下：
N=1000
m=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;
end
D=diag(VEC);
U = orth(rand(N,N));
A1 = U'*D*U;
对构造好的A1进行计算，在计算时，取X0=zeros(N,1);b=ones(N,1); 在计算时取停机准则为绝对误差e<10-10，Lanczos方法的代码同上，得出的结果如上图所示（分别绘制出log10rk和rk随迭代次数的变化图），各种情况下的迭代次数及计算时间如下表所示（这里取得矩阵为对称正定矩阵，因此Lanczos过程求解yk时采用cholesky分解）：
特征值个数
10
20
50
100
500
1000
运行时间
0.015727
0.028523
0.072200
0.107578
0.283986
0.453770
迭代次数
10
20
49
73
159
215
可以看到，几个问题的条件数虽然相同，如果A只有m个不同的特征值，则Lanczos方法至多m步就可以找到精确解。实验中，在m较大的时候，算法收敛较快，远小于m。当m较小时，可能需要接近于m步才能找到准确解。另外，由上面计算可以看到在m值较小时，其第k步的误差rk会有一个在计算时会有先增大后减小的过程，最后收敛。
T4. 取初始值近似解为零向量，右端项b仅由A的m个不同特征向量的线性组合表示时，Lanczos方法的收敛性如何？数值计算中方法的收敛性和m的大小关系如何？
Answer：
在前面的基础上，我们知道之前构造对称正定矩阵时，U的每一个行向量均为对应A的特征值（D的某一个元素）的特征向量，因此构造m不同时的又端项b代码如下所示：
N=1000
m=10;
DIA=linspace(1,N,N);
D=diag(DIA);
U = orth(randN,N));
A1 = U'*D*U;
b1=zeros(N,1);
for i=1:m
b1=b1+rand(1)*(U(i,:))';
end
X0=zeros(N,1);
e=10^(-10);
tic
[S1,ek,num]=LanczosT3(A1,X0,b1,e);
toc
计算时采用的Lanczos算法同之前的题目，计算出来的结果如上所示（这里取得矩阵为对称正定矩阵，因此Lanczos过程求解yk时采用cholesky分解）。
各种情况下的迭代次数及计算时间如下表所示：
m
10
20
50
100
500
750
1000
运行时间
0.1074
0.1950
0.2967
0.3546
0.3546
0.4325
0.4275
迭代次数
68
120
165
179
202
205
212
我们可以得到如下结论：
Ø 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(n1);
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);
end
end

而对于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为自己编写的有R·gk=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);
end
end
对于有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-10，计算结果如下：

负特征值个数：0
Lanczos 步数：212 计算时间：0.368697
MINRES 步数：210 计算时间：0.544282
负特征值个数：10
Lanczos 步数：369 计算时间：0.841989
MINRES 步数：367 计算时间：1.776133

负特征值个数：50
Lanczos 步数：697 计算时间：2.673012
MINRES 步数：689 计算时间：10.111938
负特征值个数：100
Lanczos 步数：944 计算时间：5.155608
MINRES 步数：939 计算时间：24.470218

负特征值个数：150
Lanczos 未收敛 计算时间：5.897715
MINRES 未收敛 计算时间：28.467015
负特征值个数：200
Lanczos 未收敛 计算时间：5.897715
MINRES 未收敛 计算时间：28.467015
由上述结果得到的结论如下：
Ø Lanczos方法有峰点说明，Lanczos方法可能会近似中断；
Ø 负特征值越多，峰点个数越多。实际上，更一般的结论是，当正负特征值个数相差较多时，峰点个数较少，当正负特征值个数相当时，峰点个数最多；
Ø 通过观察说明在Lanczos出现峰点时，MINRES方法收敛过程很稳定；
Ø 达到相同精度时，MINRES迭代次数较Lanczos略少，但是差别很小，而MINRES方法要比Lanczos方法花费更多的计算时间；
Ø 当负特征值个数过多时，两种方法均无法在1000步内达到收敛进度，但是MINRES的精度较Lanczos方法高。