清华大学高等数值分析第一次实验作业Word文档下载推荐.docx
《清华大学高等数值分析第一次实验作业Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《清华大学高等数值分析第一次实验作业Word文档下载推荐.docx（29页珍藏版）》请在冰豆网上搜索。
lamda=[1001,1000:
-1:
1,1]
从特征值的分布上可以看出，A的条件数为103，是个良态正定的问题。
CG结果如下图
从上图可以看出，对称正定良态问题，CG和Lanczos方法收敛性态几乎一模一样。
近似奇异的问题
上一个A矩阵中的一个特征值减小，取为10-3，特征值分布为
从收敛特性上看，对于近似奇异的正定矩阵（数值上不奇异），CG和Lanczos方法得到的结果一致。
但是CG的求解速度比Lanczos小很多。
病态的正定问题
在以上相同的A（n=1002）的基础上，将最大特征值提高为1e15，则条件数为1e15，是个非常病态的问题，同时使用CG法和Lanczos方法得到的收敛曲线下
从结果中可以看出，对于这个问题，Lanczos可能比CG的迭代次数小，但是，程序的运行时间来看，CG法为2.13s，而Lanczos为7.14s，因此对于正定问题CG法的效率要明显高于Lanczos。
良态的不定问题
在以上相同的A的基础上，增加两个负特征值-1和-10，在此基础上，使用CG和Lanczos计算，得到的结果如下
从图中可以看出，虽然问题比较良性，但是也出现了尖峰，即近似中断现象。
其中，CG法也能够在没有中断的情况下求解不定问题，但是如果一旦中断，那么无法得到结果。
但是Lanczos中断之后，却能够继续进行下去。
即使用Lanczos解不定问题的风险更小。
总结：
1）对于正定问题，不论良态还是病态，Lanczos和CG方法求解的结果非常相近，收敛曲线也近似一致，收敛步数也相差不大，这说明两个方法的原理是一致的。
但是CG方法的速度明显快于Lanczos法。
2）对于不定问题，CG方法也可能收敛，如果收敛，Lanczos方法和CG方法的到的结果也类似（步数和收敛性态）。
但是，如果CG方法一旦中断，前面的所有运算都白费了，是恶性中断，而Lanczos仍然能够继续。
因此，再不定问题中，Lanczos比CG方法的风险小。
3、当A只有m个不同特征值时,对于大的m和小的m,观察有限精度下Lanczos方法如何收敛.
对于不同的m值进行实验，得到的结果如下：
当m=20时，取特征值分别为1,2,3,4...20，每个特征值的重数均为60重，得到结果为
两个方法均在20步的时候收敛，即nit≤m。
当m=50时，取特征值分别为1,2,3,4...50，每个特征值的重数均为24重，得到的结果为
在38步左右收敛，CG法与Lanczos方法的收敛速度一致。
同样也满足nit≤m。
当m=500时，取特征值分别为1,2,3,4...500，每个特征值的重数均为2重，得到的结果为
从结果中可以看出，迭代120多步就收敛了，而且CG法与Lanczos法的结果一致。
同样满足nit≤m。
当m=1000时，取特征值分别为1,2,3,4...100，每个特征值的重数均为1重，得到的结果为
经过多次实验，将得到的迭代次数与不同特征值个数的关系做图如下，从图中可以看出，当m很小时，迭代次数等于不同特征值个数，随着m增大，迭代次数也不断增大，但是增大的趋势变缓，最终，几乎不随着m增大而增大。
同时，发现，迭代次数始终小于不同特征值的个数，即与书上讲的“若A有m个不同的特征值，则至多m步收敛”一致！
注：
图中蓝色曲线为y=x曲线。
1）特征值的分布对于Lanczos方法的收敛性态影响非常大；
2）如果A有m个不同的特征值，那么就会至多m步收敛。
4、取初始近似解为零向量,右端项b仅由A的m个不同个特征向量的线性组合表示时,Lanczos方法的收敛性如何?
数值计算中方法的收敛性和m的大小关系如何?
首先，确定一个比较良态的矩阵A，取其特征值为
lamda=[3100:
10:
2500,1000:
1]
分布如下图
分别取b为
bm=Qm（1,1,1...1）mT
其中Qm的列为A的特征向量。
变化分别取m的值为5,10,15,30,40,50,80,90,100,200,300,400,500,600,800,900,1000分别得到收敛曲线，下面仅列出m=5、15、40、50、90、100、400、600、900时的收敛曲线。
从图中可以看出，各个收敛性态都很好，曲线比较平滑。
通过实验，可以得到m与收敛迭代步长的曲线如下图
从图中可以看出，收敛性态与特征向量有关！
随着m的增大，迭代次数有上升趋势，即若b仅有少数几个特征向量的线性组合而成，那么可以得到更好的收敛效果。
1）收敛特性与特征向量有关；
2）若b由m个特征向量线性组合而成，m越小，收敛速度越快；
3）理论上，至多m步收敛，但是实验中并不是严格按照这个关系，说明特征向量问题比特征值问题复杂一些。
5、构造对称不定的矩阵,验证Lanczos方法的近似中断,观察收敛曲线中的峰点个数和特征值的分布关系;
观察当出现峰点时,MINRES方法的收敛性态怎样.
取b=（1,1,1,...,1,1）T，x0=（0,0,0,...,0,0）T，停机准则为10-8。
首先，取特征值为
1,-1];
其中有一个负特征值-1，分别使用MINRES和Lanczos方法得到的曲线如下图
从图中可以看出，在40步左右Lanczos发生近似中断现象，然而MINRES方法的相对误差单调不增，没有出现尖峰。
我们再多实验几组数据。
在前面的特征值分布基础上增加一个负特征值-10得到的收敛曲线如下图
可以看出，MIRES方法的相对误差仍然具有单调不增的特性。
再将特征值取得更多，再增加10个负特征值，-10：
-10：
-100，得到的曲线如下
通过对比可以发现，负特征值越多，Lanczos方法的纹波越大，而MINRES方法却仍然维持单调不减的特性。
1）对于Lanczos方法，随着负特征值的增多，振荡越来越严重，发生近似中断的次数越来越多；
2）然而，对于相同的A，Minres方法的相对残差没有出现峰值，随着迭代数增加而单调不减。
收敛性态比Lanczos好。
程序代码
CG法CGllb11.m
%CGmethodforsolvingAx=b
%ByLiulibin@2011-11-06
%Email:
********************
function[x,Error,i,flag]=CGllb11（A,b,x,ErrSet,uplimit）
%dealwithinputdata
[m,n]=size（b）;
ifm<
n
b=b'
;
end
[m,n]=size（x）;
x=x'
%CGmethodbegins
r=b-A*x;
p=r;
%Loopbegin
i=1;
temp_rkrkplus=r'
*r;
Error=[sqrt（temp_rkrkplus）/norm（b,2）];
while1
temp_AP=A*p;
temp_rkrk=temp_rkrkplus;
temp_pAP=p'
*temp_AP;
ifabs（temp_pAP）<
1e-12
disp（'
恶性中断！
'
）
break;
end
a=temp_rkrk/（temp_pAP）;
x=x+a*p;
r=r-a*temp_AP;
temp_rkrkplus=r'
beta=temp_rkrkplus/temp_rkrk;
p=r+beta*p;
%decidetheerror
Err=sqrt（temp_rkrkplus）/norm（b,2）;
%/（norm（b）+temp_AP）;
ifErr<
ErrSet
Methodconverge！
disp（i）
%disp（x）
flag=1;
Error=[Error,sqrt（temp_rkrkplus）/norm（b,2）];
ifi>
uplimit
flag=0;
break
End
Lanczos法lanczosllb.m
%LanczosmethodforsolvingAx=b
%ByLiulibin@2011-11-15
function[T,Q,x,k,Errs,tiao]=Lanczosllb（A,b,x0,Err）
x=[];
tiao=[];
[m,n]=size（x0）;
n
m=n;
x0=x0'
size（b）
size（A*x0）
r=b-A*x0;
r_zeros=r;
q=r/norm（r）;
q0=0;
beta0=0;
T=[0];
k=1;
Q=[q];
Errs=[norm（r,2）/norm（b,2）];
%solvefortridiagmatrixT
%clc
k='
）;
disp（k）;
T=[T,zeros（k,1）;
zeros（1,k）,0];
%AddarowandacollomtoT
r=A*q-beta0*q0;
T（k,k）=q'
r=r-T（k,k）*q;
T（k+1,k）=norm（r,2）;
beta0=T（k+1,k）;
T（k,k+1）=beta0;
q0=q;
q=r/beta0;
Tk=T（1:
k,1:
k）;
%solveforyk
ifk==1
跳过此步'
tiao=[tiao,k];
else
[L,U]=LanczosLU（Tk）;
bk=norm（r_zeros,2）*[1;
zeros（k-1,1）];
%求Ux
Ux=zeros（k,1）;
Ux
（1）=bk
（1）;
fori=2:
k
Ux（i）=bk（i）-L（i,i-1）*Ux（i-1）;
%求ym
ym=zeros（k,1）;
ym（k）=Ux（k）/U（k,k）;
fori=k-1:
1
ym（i）=（Ux（i）-ym（i+1）*U（i,i+1））/U（i,i）;
Errs=[Errs,abs（ym（k）*beta0）/norm（b,2）];
ifabs（ym（k）*beta0）/norm（b,2）<
Err
x=x0+Q*ym;
methodconverge,gottheaccurateresult!
elseifbeta0<
1e-10
GoodLuck!
%x=x0+Q*ym;
ifk==1.1*m
NotConverge!
Q=[Q,q];
k=k+1;
MINRES法MINRES.m
%MinresmethodforsolvingAx=b
function[T,Q,x,k,Errs]=minresllb（A,b,x0,Err）
Errs=[norm（r_zeros,2）];
clc
%*********Lanzcosprocesscompleted!
%InterruptionforUplimitedtimes
ifk<
=2
%*******starttodevideTk^byQRmethod
Rk=T（1:
k+1,1:
k）;
Qk=eye（k+1,k）;
fori=1:
Ci=Rk（i,i）/（（Rk（i,i）^2+Rk（i+1,i）^2）^0.5）;
Si=Rk（i+1,i）/（（Rk（i,i）^2+Rk（i+1,i）^2）^0.5）;
Rk_temp=Rk;
Rk_temp（i,:
）=Ci*Rk（i,:
）+Si*Rk（i+1,:
Rk_temp（i+1,:
）=-Si*Rk（i,:
）+Ci*Rk（i+1,:
Rk=Rk_temp;
%computeQk
Qk_temp=Qk;
Qk_temp（i,:
）=Ci*Qk（i,:
）+Si*Qk（i+1,:
Qk_temp（i+1,:
）=-Si*Qk（i,:
）+Ci*Qk（i+1,:
Qk=Qk_temp;
%****QRdevisioncompleted
Errsk=abs（Qk（k+1,1））*norm（r_zeros,2）/norm（b,2）;
Errs=[Errs;
Errsk];
ifErrsk<
Err||k==m
gk=Qk（1:
k,1）*norm（r_zeros,2）;
yk=zeros（k,1）;
yk（k）=gk（k）/Rk（k,k）;
yk（k-1）=（gk（k-1）-yk（k）*Rk（k-1,k））/Rk（k-1,k-1）;
fori=k-2:
yk（i）=（gk（i）-yk（i+1）*Rk（i,i+1）-yk（i+2）*Rk（i,i+2））/Rk（i,i）;
Err
MethodConverge！
Congratulations！
Failed!
x=x0+Q（:
1:
k）*yk;
其他附属过程函数
A和b的生成矩阵MakeMatrix.m
clc
clear
%buildMatrixA
%lamda
globalA
globalb
beta=0.5;
%lamda=[1e9,1000:
1,0.01];
%正定的病态问题
1,-1,-10:
-10:
-100];
%正定的良态问题
%lamda=[1001,1000:
1,0.001];
%近似奇异问题
%lamda=[1e15,1000:
1,1];
%病态正定问题
1,-1,-10,-100];
1,-0.1,-10];
%不定的良态问题
%lamda=[100:
1,-100:
0.5:
-1,0];
%近似不定的良态问题
%lamda=[100*ones（1,210）,-10*ones（1,210）,1*ones（1,210）,1e3*ones（1,210）]
%lamda=[1:
1100,1e12,0.001];
%lamda=[1e9,1,2,3,4,5,6,5,8,9,10,11,4,78];
%lamda=[100,50,20,10,1];
n=length（lamda）
semilogx（lamda,ones（n,1）,'
bo'
Q=orth（randn（n,n））;
title（'
A的特征值分布情况'
xlabel（'
特征值'
set（gca,'
ytick'
[]）
b=ones（n,1）;
B=mdiagllb（lamda）;
A=Q'
*B*Q;
disp（'
矩阵生成完毕！
可以进行后续计算！
disp（A（1:
5））
对称三对角矩阵的LU分解LanczosLU.m
function[L,U]=LanczosLU（A）
[m,n]=size（A）;
L=eye（n,n）;
U=eye（n,n）;
U（1,1）=A（1,1）;
%L（2,1）=A（2,1）/A（1,1）;
U（1,2）=A（1,2）;
fork=2:
m-1
L（k,k-1）=A（k,k-1）/U（k-1,k-1）;
U（k,k）=A（k,k）-L（k,k-1）*U（k-1,k）;
U（k,k+1）=A（k,k+1）-L（k,k-1）*U（k-1,k+1）;
k=m;
L（m,m-1）=A（m,m-1）/U（m-1,m-1）;
实验控制文件（调用这一系列函数）Fortest.m
clc
globalb
n_Matrix=m;
%调用CG
uplimit=2000;
tic
[y,error,n]=CGllb11（A,b,zeros（n_Matrix,1）,1e-8,uplimit）;
toc
figure
（2）
semilogy（[1:
length（error）],error,'
b.-'
norm（A*y-b）
stringsCG='
算法的收敛曲线（阶数n='
stringsCG=[stringsCG,num2str（m）,'
）'
];
title（stringsCG）
迭代次数'
ylabel（'
||r_k||/||b||'
gridon
%调用Lanczos
[~,~,x,k,Errs]=Lanczosllb（A,b,zeros（n_Matrix,1）,1e-8）;
最大的误差为'
disp（max（A*x-b））
holdon
semilogy（1:
length（Errs）,Errs,'
r.-'
%调用系统函数symmlq
[x,flag,relres,iter,rv]=symmlq（A,b,1e-8,500）%系统函数
semilogy（1:
length（rv）,rv/norm（b,2）,'
g.-'
%调用MINRES
[T,Q,x,k,Errs]=minresllb（A,b,zeros（n_Matrix,1）,1e-8）;
holdon
k-1,Errs,'
legend（'
CG'
'
Lanczos'
Symmlq@Matlab'
Minres'
升级会员 