清华大学高等数值分析第一次实验作业.docx
《清华大学高等数值分析第一次实验作业.docx》由会员分享,可在线阅读,更多相关《清华大学高等数值分析第一次实验作业.docx(29页珍藏版)》请在冰豆网上搜索。
![清华大学高等数值分析第一次实验作业.docx](https://file1.bdocx.com/fileroot1/2022-11/25/92de13ee-1d50-47ac-b124-de0b712e2717/92de13ee-1d50-47ac-b124-de0b712e27171.gif)
清华大学高等数值分析第一次实验作业
高等数值分析第二次作业第四题
姓名:
----学号:
-------时间:
11月20日
1、构造例子说明CG的数值性态.当步数=阶数时CG的解如何?
当A的最大特征值远大于第二个最大特征值,最小特征值远小于第二个最小特征值时方法的收敛性如何?
答:
1)构造矩阵如下
首先,构造一个非常病态的矩阵。
构造矩阵的条件数K(A)为1015,可以得到收敛曲线如下图所示
第二种方法构造矩阵(良态):
条件数不大的矩阵也能够使迭代步数大于等于矩阵阶数,构造一个良态的矩阵如下:
使其特征值为等比数列,即
由此得到的矩阵A最大特征值为105,最小为1,的条件数是105,特征值分布如下图
收敛曲线如下图,发现使用CG法仍然可以收敛,但是收敛速度非常缓慢,即使到1000步也难以达到理想的残差,说明即使问题的条件数较小,即良态的,CG方法也可能收敛非常慢。
从上面的结果中可以看出,当矩阵的迭代步数等于阶数时
a.往往收敛曲线不平滑,震荡非常严重,残差的2范数没有最优性;
b.但是方法的相对残差大体上仍然具有下降趋势,最终仍然能够收敛。
矩阵的病态性延迟了方法的收敛。
c.条件数较小的矩阵也会出现收敛缓慢的现象,这与特征值的分布有关。
2)构造矩阵1002阶的A,其中,A的最大的特征值为109,最小特征值为1,其具有如下的特征值分布
图1A的特征值分布
该矩阵的条件数k=1011,是个病态问题,结果如下
图2A的收敛性态图
代入CG法,得到的收敛曲线如图2所示。
从图中可以看出,一开始,方法残差下降很慢,等到后面,在150步左右,以另一种收敛速度下降。
实验结果与理论结果相同。
3)结论:
a.对于非常病态的问题,CG法的迭代次数可能大于等于阶数,残差震荡非常严重,但是大体趋势仍然是下降的,最终也能够收敛;对于极其病态的问题,CG法也可能无能为力,无法收敛。
b.对于最大特征值远远大于第二大特征值的情况,收敛曲线出现了两种下降速度,说明特征值的分布对于CG方法的收敛速度影响很大。
c.
2、对于同样的例子,比较CG和Lanczos的计算结果.
在相同的A下,取b=(1,1,1,1...1)T,x0=(0,0,0...0)T,分别用CG法和Lanczos法求解。
停机准则为相对误差小于10-8。
将A取以下几种:
良态正定矩阵、近似奇异正定矩阵、病态正定矩阵、良态不定矩阵进行实验,结果分别如下:
Ø良态的正定问题
取矩阵A的特征值分布如下:
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:
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与收敛迭代步长的曲线如下图
注:
图中蓝色曲线为y=x曲线。
从图中可以看出,收敛性态与特征向量有关!
随着m的增大,迭代次数有上升趋势,即若b仅有少数几个特征向量的线性组合而成,那么可以得到更好的收敛效果。
注:
图中蓝色曲线为y=x曲线。
总结:
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。
首先,取特征值为
lamda=[3100:
10:
2500,1000:
-1:
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);
ifmb=b';
end
[m,n]=size(x);
ifmx=x';
end
%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'*r;
beta=temp_rkrkplus/temp_rkrk;
p=r+beta*p;
%decidetheerror
Err=sqrt(temp_rkrkplus)/norm(b,2);%/(norm(b)+temp_AP);
ifErrdisp('Methodconverge!
')
disp(i)
%disp(x)
flag=1;
break;
end
Error=[Error,sqrt(temp_rkrkplus)/norm(b,2)];
ifi>uplimit
flag=0;
break
end
end
End
Lanczos法lanczosllb.m
%LanczosmethodforsolvingAx=b
%ByLiulibin@2011-11-15
%Email:
********************
function[T,Q,x,k,Errs,tiao]=Lanczosllb(A,b,x0,Err)
%dealwithinputdata
x=[];
tiao=[];
[m,n]=size(b);
ifmb=b';
end
[m,n]=size(x0);
ifmm=n;
x0=x0';
end
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
while1
%clc
disp('k=');disp(k);
T=[T,zeros(k,1);zeros(1,k),0];
%AddarowandacollomtoT
r=A*q-beta0*q0;
T(k,k)=q'*r;
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
disp('跳过此步')
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);
end
%求ym
ym=zeros(k,1);
ym(k)=Ux(k)/U(k,k);
fori=k-1:
-1:
1
ym(i)=(Ux(i)-ym(i+1)*U(i,i+1))/U(i,i);
end
Errs=[Errs,abs(ym(k)*beta0)/norm(b,2)];
ifabs(ym(k)*beta0)/norm(b,2)x=x0+Q*ym;
disp('methodconverge,gottheaccurateresult!
')
break;
elseifbeta0<1e-10
x=x0+Q*ym;
disp('GoodLuck!
')
end
end
%x=x0+Q*ym;
end
ifk==1.1*m
%x=x0+Q*ym;
disp('NotConverge!
')
break;
end
Q=[Q,q];
k=k+1;
end
end
MINRES法MINRES.m
%MinresmethodforsolvingAx=b
%ByLiulibin@2011-11-15
%Email:
********************
function[T,Q,x,k,Errs]=minresllb(A,b,x0,Err)
%dealwithinputdata
[m,n]=size(b);
ifmb=b';
end
[m,n]=size(x0);
ifmm=n;
x0=x0';
end
r=b-A*x0;
r_zeros=r;
q=r/norm(r);
q0=0;
beta0=0;
T=[0];
k=1;
Q=[q];
Errs=[norm(r_zeros,2)];
%solvefortridiagmatrixT
while1
clc
disp('k=');disp(k);
T=[T,zeros(k,1);zeros(1,k),0];
%AddarowandacollomtoT
r=A*q-beta0*q0;
T(k,k)=q'*r;
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;
Q=[Q,q];
%*********Lanzcosprocesscompleted!
%InterruptionforUplimitedtimes
ifk<=2
else
%*******starttodevideTk^byQRmethod
Rk=T(1:
k+1,1:
k);
Qk=eye(k+1,k);
fori=1:
k
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;
end
%****QRdevisioncompleted
Errsk=abs(Qk(k+1,1))*norm(r_zeros,2)/norm(b,2);
Errs=[Errs;Errsk];
ifErrskgk=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:
-1:
1
yk(i)=(gk(i)-yk(i+1)*Rk(i,i+1)-yk(i+2)*Rk(i,i+2))/Rk(i,i);
end
ifErrskdisp('MethodConverge!
Congratulations!
')
else
disp('Failed!
')
end
x=x0+Q(:
1:
k)*yk;
break;
end
end
k=k+1;
end
end
其他附属过程函数
A和b的生成矩阵MakeMatrix.m
clc
clear
%buildMatrixA
%lamda
globalA
globalb
beta=0.5;
%lamda=[1e9,1000:
-1:
1,0.01];%正定的病态问题
lamda=[3100:
10:
2500,1000:
-1:
1,-1,-10:
-10:
-100];%正定的良态问题
%lamda=[1001,1000:
-1:
1,0.001];%近似奇异问题
%lamda=[1e15,1000:
-1:
1,1];%病态正定问题
%lamda=[1001,1000:
-1:
1,-1,-10,-100];%病态正定问题
%lamda=[1001,1000:
-1:
1,-0.1,-10];%不定的良态问题
%lamda=[100:
-1:
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);
end
k=m;
L(m,m-1)=A(m,m-1)/U(m-1,m-1);
U(k,k)=A(k,k)-L(k,k-1)*U(k-1,k);
End
实验控制文件(调用这一系列函数)Fortest.m
clc
globalA
globalb
[m,n]=size(A);
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)
xlabel('迭代次数')
ylabel('||r_k||/||b||')
gridon
%调用Lanczos
tic
[~,~,x,k,Errs]=Lanczosllb(A,b,zeros(n_Matrix,1),1e-8);
toc
disp('最大的误差为')
disp(max(A*x-b))
figure
(2)
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);
figure
(2)
holdon
semilogy(1:
k-1,Errs,'g.-');
disp('最大的误差为')
disp(max(A*x-b))
legend('CG','Lanczos','Symmlq@Matlab','Minres')