清华大学高等数值分析第一次实验作业Word文档下载推荐.docx
- 文档编号:16658827
- 上传时间:2022-11-25
- 格式:DOCX
- 页数:29
- 大小:341.97KB
清华大学高等数值分析第一次实验作业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'
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 清华大学高等数值分析 第一次实验作业 清华大学 高等 数值 分析 第一次 实验 作业