賈哥高等數(shù)值分析第一次實驗(共18頁)

1、高等(godng)數(shù)值分析第一次實驗第一(dy)題：構(gòu)造例子說明CG的數(shù)值(shz)形態(tài)。當(dāng)步數(shù) = 階數(shù)時CG的解如何？當(dāng)A的最大特征值遠(yuǎn)大于第二個最大特征值，最小特征值遠(yuǎn)小于第二個最小特征值時，方法的收斂性如何？解：用Housholder變換和對角陣構(gòu)造1000階正定對稱矩陣A：構(gòu)造對角陣D = diag( linspace(1, 1000, 1000) )；構(gòu)造Householder陣H。取單位向量w=1,0,0,.0T，I為1000階單位矩陣，H = I wTw。構(gòu)造對稱正定矩陣A。A = HTDH。由于D是對角陣，H是對稱的，所以A對稱；且A與D具有相同的特征值linspace(1,

2、 1000, 1000) 0，因此A對陣正定。b=rand(1000,1);取初始解x0=zeros(1000,1);1.計算Ax = b利用matlab編程實現(xiàn)CG算法。由于實際計算存在機(jī)器誤差，因此迭代1000步后的殘差不等于0，因此不能用rk=0作為停機(jī)準(zhǔn)則，否則matlab會無休止地計算下去。本例采用停機(jī)準(zhǔn)則為：迭代步數(shù)=1000步。當(dāng)D = diag( linspace(1, 1000, 1000) )時，條件數(shù)k=1000；當(dāng)D = diag( linspace(1, 100, 1000) )時，條件數(shù)k=100；當(dāng)D = diag( linspace(1, 10, 1000) )

3、時，條件數(shù)k=10；分別計算上述三種條件數(shù)下的CG算法，得到迭代步數(shù)與殘差的曲線圖。圖1：log(rk)與步數(shù)關(guān)系曲線。橫坐標(biāo)是迭代步數(shù)，縱坐標(biāo)是殘差的對數(shù)值。圖 SEQ 圖 * ARABIC 1如圖1所示，矩陣(j zhn)A的條件數(shù)k越小，CG法收斂(shulin)速度越快。附matlab程序(chngx)1-1：clear allclc%條件數(shù)k=1000D=diag(linspace(1,1000,1000);w=eye(1,1000);I=eye(1000);H=I-w*w;A=H*D*H;%生成1000階的對稱矩陣b=rand(1000,1);x=zeros(1000,1);%初始

4、解r=b-A*x;%初始?xì)埩縫=r;%初始搜索方向k=0;semilogy(0,norm(r),r-);hold on;while k1000 alpha = r*p/(p*A*p); x = x+alpha*p; rold = r; r = rold-alpha*A*p; beta = r*r/(rold*rold); p = r+beta*p; semilogy(k,norm(r),r-); hold on; k=k+1;end%條件數(shù)k=100clear allD=diag(linspace(1,100,1000);w=eye(1,1000);I=eye(1000);H=I-w*w;A=

5、H*D*H;%生成1000階的對稱矩陣b=rand(1000,1);x=zeros(1000,1);%初始解r=b-A*x;%初始?xì)埩縫=r;%初始搜索方向k=0;semilogy(0,norm(r),b-);hold on;while k1000 alpha = r*p/(p*A*p); x = x+alpha*p; rold = r; r = rold-alpha*A*p; beta = r*r/(rold*rold); p = r+beta*p; semilogy(k,norm(r),b-); hold on; k=k+1;end%條件(tiojin)數(shù)k=10clear allD=di

6、ag(linspace(1,10,1000);w=eye(1,1000);I=eye(1000);H=I-w*w;A=H*D*H;%生成1000階的對稱(duchn)矩陣b=rand(1000,1);x=zeros(1000,1);%初始(ch sh)解r=b-A*x;%初始?xì)埩縫=r;%初始搜索方向k=0;semilogy(0,norm(r),black-);hold on;while k1000 alpha = r*p/(p*A*p); x = x+alpha*p; rold = r; r = rold-alpha*A*p; beta = r*r/(rold*rold); p = r+be

7、ta*p; semilogy(k,norm(r),black-); hold on; k=k+1;endtitle(條件數(shù)的大小對CG法收斂特性的影響);xlabel(迭代步數(shù))ylabel(殘差對數(shù)log(|rk|) 2.構(gòu)造特殊特征值分布構(gòu)造對稱正定矩陣A1和A2。D1=diag( linspace(1, 1000, 1000) )時，條件數(shù)k=1000，特征值均勻分布；D2=diag(1,linspace(500,600,998),1000)時，條件(tiojin)數(shù)仍為k=1000，最大特征值1000遠(yuǎn)大于第二個最大特征值600，最小特征值1遠(yuǎn)小于第二個最小特征值500。圖2：矩陣特征

8、值分布對CG算法(sun f)收斂性的影響圖 SEQ 圖 * ARABIC 2如圖2所示，A1和A2的條件(tiojin)數(shù)均為1000，但A2的收斂速度遠(yuǎn)高于A1。這是因為，在CG算法中，系數(shù)矩陣的中間特征值分布對CG的收斂速度有巨大的影響。經(jīng)過幾步后，CG的收斂因子將是：=0.046而非：=0.939因此，A2矩陣的收斂速度快得多。附matlab程序1-2：clear allclc%條件數(shù)k=1000,特征值均勻分布D=diag(linspace(1,1000,1000);w=eye(1,1000);I=eye(1000);H=I-w*w;A=H*D*H;%生成(shn chn)1000階

9、的對稱矩陣b=rand(1000,1);x=zeros(1000,1);%初始(ch sh)解r=b-A*x;%初始(ch sh)殘量p=r;%初始搜索方向k=0;semilogy(0,norm(r),r-);hold on;while k1000 alpha = r*p/(p*A*p); x = x+alpha*p; rold = r; r = rold-alpha*A*p; beta = r*r/(rold*rold); p = r+beta*p; semilogy(k,norm(r),r-); hold on; k=k+1;end%條件數(shù)k=1000，最大特征值遠(yuǎn)大于第二個最大特征值，最

10、小特征值遠(yuǎn)小于第二個最小特征值clear allD=diag(1,linspace(500,600,998),1000);w=eye(1,1000);I=eye(1000);H=I-w*w;A=H*D*H;%生成1000階的對稱矩陣b=rand(1000,1);x=zeros(1000,1);%初始解r=b-A*x;%初始?xì)埩縫=r;%初始搜索方向k=0;semilogy(0,norm(r),b-);hold on;while k1000 alpha = r*p/(p*A*p); x = x+alpha*p; rold = r; r = rold-alpha*A*p; beta = r*r/(

11、rold*rold); p = r+beta*p; semilogy(k,norm(r),b-); hold on; k=k+1;endtitle(特征值分布對CG法收斂(shulin)特性的影響);xlabel(迭代(di di)步數(shù))ylabel(殘差對數(shù)(du sh)log(|rk|) 第二題對于同樣的例子，比較CG和Lanczos的計算結(jié)果解：采用和第一題相同的構(gòu)造方法，構(gòu)造三個1000階正定對稱矩陣，使條件數(shù)k分別為：1000,100,10。分別采用CG和Lanczos方法計算Ax=b，且都設(shè)置停機(jī)準(zhǔn)則為：norm(rk)1 bita(j)=norm(r0); Q(:,j+1)=r0

12、/bita(j); end for k=1:j T(k,k)=alpha(k); end for k=1:j-1 T(k+1,k)=bita(k); T(k,k+1)=bita(k); %生成三對角陣T end e(1)=1; e(2:j)=0; y1=T(y*e); W=Q(:,1:j); X=X0+W*y1; r=norm(b-A*X); %求解(qi ji)第k步生成的X及r F(j)=r; if r/norm(b)1e-8 break; end %r的精度達(dá)到要求后停止迭代，得到(d do)最終X end semilogy(F,r);hold on;j %迭代(di di)步數(shù)toc%

13、=CG算法=ticp=R0;for i=1:n R=R0; a=(R*R)/(p*(A*p); x=x+a*p; R=R-a*(A*p); G(i)=norm(R); if G(i)/norm(b)=1e-8 break; end beta=(R*R)/(R0*R0); p=R+beta*p; R0=R; endsemilogy(G,b);toci %迭代步數(shù)legend(Lanczos,CG)title(Lanczos與CG算法收斂性對比,條件數(shù)k=10)xlabel(迭代步數(shù))ylabel(殘差對數(shù)log(|rk|) 第三題當(dāng)A只有m個不同特征值時，對于大的m和小的m，觀察有限精度下Lan

14、czos方法如何收斂。解： 分別(fnbi)構(gòu)建m = 10、50、100、1000四個矩陣A，設(shè)置停機(jī)準(zhǔn)則為：norm(rk)1 bita(j)=norm(r0); Q(:,j+1)=r0/bita(j); end for k=1:j T(k,k)=alpha(k); end for k=1:j-1 T(k+1,k)=bita(k); T(k,k+1)=bita(k); %生成三對角陣T end e(1)=1; e(2:j)=0; y1=T(y*e); W=Q(:,1:j); X=X0+W*y1; r=norm(b-A*X); %求解第k步生成的X及r F(j,i)=r; if r/norm

15、(b)300后，迭代步數(shù)遠(yuǎn)小于m。而同樣由于計算精度限制，在m較小時，迭代步數(shù)可能溢出m步，比如本例m=65和80的情形。另外，m值大于一定值后，隨著m的增加，收斂步數(shù)逐漸趨于定值，可能的原因是：隨著m的增大，特征向量線性相關(guān)性逐漸增強(qiáng)。附matlab程序4：clear allclcn=1000; D=diag(linspace(1,1000,n);P1=rand(n);Q1,R=qr(P1);A=Q1*D*Q1;%生成1000階的對稱矩陣V,D1=eig(A);for i=1:7M=10,50,65,80,300,800,1000;m=M(i);b=m*mean(V(:,1:m),2);X0

16、=zeros(n,1);%初始解x=X0; r0=b-A*X0; R0=r0; %=Lanczos解法 =ticy=norm(r0);F(1)=norm(r0); Q(:,1)=r0/norm(r0); r0=A*Q(:,1); alpha(1)=Q(:,1)*r0; r0=r0-alpha(1)*Q(:,1); bita(1)=norm(r0); Q(:,2)=r0/bita(1); %給各變量賦初始值 for j=2:n r0=A*Q(:,j)-bita(j-1)*Q(:,j-1); alpha(j)=Q(:,j)*r0; r0=r0-alpha(j)*Q(:,j); if norm(r0

17、)1 bita(j)=norm(r0); Q(:,j+1)=r0/bita(j); end for k=1:j T(k,k)=alpha(k); end for k=1:j-1 T(k+1,k)=bita(k); T(k,k+1)=bita(k); %生成(shn chn)三對角陣T end e(1)=1; e(2:j)=0; y1=T(y*e); W=Q(:,1:j); X=X0+W*y1; r=norm(b-A*X); %求解(qi ji)第k步生成的X及r F(j,i)=r; if rsqrt(eps) break; end %r的精度(jn d)達(dá)到要求后停止迭代，得到最終X end

18、tocjclear e y1 X r Tendsemilogy(F(:,1),r);hold on;semilogy(F(:,2),b);hold on;semilogy(F(:,3),y);hold on;semilogy(F(:,4),black);hold on;semilogy(F(:,5),g);hold on;semilogy(F(:,6),m);hold on;semilogy(F(:,7),c);hold on;legend(m=10,m=50,m=65,m=80,m=300,m=800,m=1000)title(b由不同特征向量個數(shù)組成時Lanczos的收斂性)xlabel(迭代步數(shù))ylabel(殘差對數(shù)log(|rk|)第五題構(gòu)造對稱不定矩陣，驗證Lanczos方