高等數(shù)值分析作業(yè)-第一次實(shí)驗(yàn)（總

8頁）本頁僅作為預(yù)覽文檔封面，使用時請刪除本頁-#YKN函數(shù)：n=size(TK,1);yk=zeros(n,1);yk(n)=b(n)/TK(n,n);fori=1:n-1j=n-i；yk(j)=(b(j)-TK(j,j+1)*yk(j+1))/TK(j,j)end首先構(gòu)造1000階有m個特征值的矩陣，1000YKN函數(shù)：n=size(TK,1);yk=zeros(n,1);yk(n)=b(n)/TK(n,n);fori=1:n-1j=n-i；yk(j)=(b(j)-TK(j,j+1)*yk(j+1))/TK(j,j)end首先構(gòu)造1000階有m個特征值的矩陣，1000)，構(gòu)造矩陣的代碼如下：N=1000m=10;DIA=linspace(1,1000,m);VEC=zeros(N,1);k=N/m;i=1；ii=0;whilei<=mforj=1:kVEC(ii+j)=DIA(i);endii=ii+k;i=i+1;構(gòu)造是保證矩陣的條件數(shù)相同（這里均取為kg 10個特征值 20個特征值 50個特征值 100個特征值 500個特征值1000個特征值-12100 150 200迭代次數(shù)2502010100908070605040 10個特征值 20個特征值 50個特征值 100個特征值 500個特征值1000個特征值00 50100 150迭代次數(shù)200 250lyk=YK(L',norm(r0)*eye(k,1));yk=YKN(L,lyk);YK函數(shù)：n=size(TK,1);yk=zeros(n,1);yk(1)=b(1)/TK(1,1);fori=2:nyk(i)=(b(i)-TK(i,i-1)*yk(i-1))/TK(i,i);end計算結(jié)果如上圖所示，計算時間與迭代步數(shù)如下:計算方法CG方法Lanczos方法迭代步數(shù)234235運(yùn)行時間（s）可以得到如下結(jié)論：當(dāng)矩陣A為對稱正定時，兩種方法效果相當(dāng)，每一步誤差也基本相同，但收斂速度基本一樣。由于Lanczos方法的每一步迭代中都有一個Lanczos過程，其中需要構(gòu)造Tk和Qk，以及計算yk，故該方法需要耗費(fèi)更長的計算時間。同時計算過程中發(fā)現(xiàn)，當(dāng)取得收斂準(zhǔn)則比較嚴(yán)格時，CG算法較Lanczos方法穩(wěn)定。T3.當(dāng)A只有m個不同特征值時，對于大的m和小的m，觀察有限精度下Lanczos方法如何收斂。Answer：

endD=diag(VEC);U=orth(rand(N,N));A1=U'*D*U;對構(gòu)造好的A1進(jìn)行計算，在計算時，取X0=zeros(N,1);b=ones(N,1);在計算時取停機(jī)準(zhǔn)則為絕對誤差e<10-i0,Lanczos方法的代碼同上，得出的結(jié)果如上圖所示(分別繪制出log10ll*ll和llrMI隨迭代次數(shù)的變化圖)，各種情況下的迭代次數(shù)及計算時間如下表所示(這里取得矩陣為對稱正定矩陣，因此Lanczos過程求解yk時采用cholesky分解)：特征值個數(shù)1020501005001000運(yùn)行時間迭代次數(shù)10204973159215可以看到，幾個問題的條件數(shù)雖然相同，如果A只有m個不同的特征值，則Lanczos方法至多m步就可以找到精確解。實(shí)驗(yàn)中，在m較大的時候，算法收斂較快，遠(yuǎn)小于m。當(dāng)m較小時，可能需要接近于m步才能找到準(zhǔn)確解。另外，由上面計算可以看到在m值較小時，其第k步的誤差llrkll會有一個在計算時會有先增大后減小的過程，最后收斂。T4.取初始值近似解為零向量，右端項(xiàng)b僅由A的m個不同特征向量的線性組合表示時，Lanczos方法的收斂性如何？數(shù)值計算中方法的收斂性和m的大小關(guān)系如何？Answer：在前面的基礎(chǔ)上，我們知道之前構(gòu)造對稱正定矩陣時，U的每一個行向量均為對應(yīng)A的特征值(D的某一個元素)的特征向量，因此構(gòu)造m不同時的又端項(xiàng)b代碼如下所示：N=1000m=10;k|gg-6-8-10k|gg-6-8-100 50 100 150 200 250迭代次數(shù)-12D=diag(DIA);U=orth(rand(N,N));A1=U'*D*U;b1=zeros(N,1);fori=1:mb1=b1+rand(1)*(U(i,:))';endX0=zeros(N,1);e=10八(-10);tic[S1,ek,num]=LanczosT3(A1,X0,b1,e);toc計算時采用的Lanczos算法同之前的題目，計算出來的結(jié)果如上所示(這里取得矩陣為對稱正定矩陣，因此Lanczos過程求解yk時采用cholesky分解)。各種情況下的迭代次數(shù)及計算時間如下表所示：m1020501005007501000運(yùn)行時間

迭代次數(shù)68120165179202205212我們可以得到如下結(jié)論：M較小時，計算過程中收斂誤差會產(chǎn)生波動，m較大時，波動減小，最后均收斂；收斂步數(shù)隨著m值的增大而增加，但是由上表我們可以看出收斂步數(shù)的增幅逐漸減小，可以推測當(dāng)m達(dá)到一定值后，收斂步數(shù)可能趨近于某一定值。T5.構(gòu)造對稱不定矩陣，驗(yàn)證Lanczos方法的近似中斷，觀察收斂曲線中的峰點(diǎn)個數(shù)和特征值的分布關(guān)系；觀測當(dāng)出現(xiàn)峰點(diǎn)時，MINRES方法的收斂形態(tài)怎樣。Answer：對于對稱不定矩陣，求解yk時采用追趕法，追趕法的代碼如下所示:functions=Zhuigan(A,bbfunctions=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);fori=1:n-1d(1)=b(1);u(i)=c(i)/d(i);d(i+1)=b(i+1)-a(i)*u(i);end而對于MINRES方法，% 追的過程 y=zeros(n,1);y(1)=bb(1)/d(1);fori=2:ny(i)=(bb(i)-a(i-1)*y(i-1))/d(i);end% 趕的過程 s(n)=y(n);fori=n-1:-1:1s(i)=y(i)-u(i)*s(i+1);Lanczos算法不同之處在于，ek=zeros(k,1);ek(k)=1;endend考慮求解yk時用MATLAB自帶的qr分解進(jìn)行計算，其中與求解yk部分代碼如下：TKT=[TKbj(k)*ek]';[Q,R]=qr(TKT);gk=Q'*norm(r0)*eye(k+1,1);yk=minresYK(R,gk);其中，：^口儂丫長為自己編寫的有R-gk=yk求解yk的函數(shù)，代碼如下:functionyk=minresYK(TK,b)n=size(TK,2);yk=zeros(n,1);yk(n)=b(n)/TK(n,n);fori=1:n-1k=n-i;yk(k)=b(k);forj=k+1:nyk(k)=yk(k)-TK(k,j)*yk(j);endyk(k)=yk(k)/TK(k,k);

endend對于有m個特征值的對稱不定矩陣，本文中采用的生成代碼如下(基本思想同之前相同，生成矩陣的特征值分別為-m,-(m-1),……,-1,1,2,……N-m，其中N為矩陣階數(shù)，計算過程中按1000考慮)：N=1000%10個負(fù)特征值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);在計算時取停機(jī)準(zhǔn)則為絕對誤差e<10-i0，計算結(jié)果420)|k|r(-2-4lg-6-8420)|k|r(-2-4lg-6-8-100 50 100150迭代20次0數(shù)250 300 35010-12400負(fù)特征值個數(shù)：Lanczos步數(shù)：Lanczos步數(shù)：212計算時間：MINRES步數(shù)：210計算時間：Lanczos步數(shù)：369計算時間：MINRES步數(shù)：367計算時間：420-4-6-8-10|r(lg) -20 100 200 300 400 500 600 700迭代次數(shù)420-4-6-8-10|r(lg) -20 100 200 300 400 500 600 700迭代次數(shù)-12負(fù)特征值個數(shù)：50Lanczos步數(shù)：697計算時間：MINRES步數(shù)：689計算時間：Lanczos步數(shù)：944計算時間：MINRES步數(shù)：939計算時間：負(fù)特征值個數(shù)：200Lanczos未收斂計算時間：MINRES未收斂計算時間：Lanczos未收斂計算時間：負(fù)特征值個數(shù)：200Lanczos未收斂計算時間：MINRES未收斂計算時間：Lanczos方法有峰點(diǎn)說明，Lanczos方法可能會近似中斷；負(fù)特征值越多，峰點(diǎn)個數(shù)越多。實(shí)際上，更一般的結(jié)論是，當(dāng)正