數(shù)值分析計(jì)算實(shí)習(xí)第一次大作業(yè)

1、【精品文檔】如有侵權(quán)，請(qǐng)聯(lián)系網(wǎng)站刪除，僅供學(xué)習(xí)與交流數(shù)值分析計(jì)算實(shí)習(xí)第一次大作業(yè).精品文檔.學(xué)生：黃小琦 學(xué)號(hào)：S201150110 聯(lián)系方式法：1.編輯反冪法的子程序1）設(shè)置二維數(shù)組a（5，501）來(lái)存儲(chǔ)矩陣A 的帶內(nèi)元素，Aij=ai-j+2,j；2）對(duì)于主對(duì)角線上的元素，即a2j設(shè)置平移量n，并存儲(chǔ)到原來(lái)單元內(nèi)，使得對(duì)于j=0，1，2，500，a2j=a2j-n；3）對(duì)矩陣A進(jìn)行LU分解，相對(duì)應(yīng)于?。?,501）中的元素執(zhí)行如下循環(huán)：對(duì)k=0，1，2，500執(zhí)行 j=k,k+1,，min(k+2,500) i=k+1,k+2,，min(k+2,500)；k<

2、;n4）若平移量x=0，即對(duì)矩陣A進(jìn)行LU分解，所得a(5,501)第二行元素的乘積即為矩陣A的行列式的值，detA=a0,0*a0,1*a0,500;5）進(jìn)行反冪法迭代前，先對(duì)平移后的矩陣的行列式做判斷，看其是否為0，若行列式為0，則絕對(duì)值最小的特征值為0。相對(duì)于原來(lái)的矩陣A，其特征值為平移量x。若行列式不為0，則開(kāi)始進(jìn)行反冪法的迭代求解；6）取非零向量組u(501)作為初始值，取u的所有元素都為1，每迭代一次就將新求得的u(501)存儲(chǔ)于原來(lái)的存儲(chǔ)單元內(nèi)，覆蓋原來(lái)的值，執(zhí)行如下循環(huán)：對(duì)k=0，1，2，500（500為最大迭代次數(shù)，超過(guò)這個(gè)次數(shù)仍不收斂則停止計(jì)算）設(shè)置一個(gè)數(shù)comp來(lái)存儲(chǔ)前一

3、次運(yùn)行所得的特征值，以便與第二次運(yùn)行的值相比較，求出誤差范圍，comp=q；利用3）矩陣的LU分解結(jié)果求解方程組,求解Lz=y和Uu=z的計(jì)算公式是： (i=1，2，500) (i=499，498，0)，判斷，若，結(jié)束循環(huán)轉(zhuǎn)7），否則繼續(xù)循環(huán)7）將q值帶回主函數(shù)。2.編輯冪法的子程序1）設(shè)置二維數(shù)組a（5，501）來(lái)存儲(chǔ)矩陣A 的帶內(nèi)元素，Aij=ai-j+2,j；2）對(duì)于主對(duì)角線上的元素，即a2j設(shè)置平移量n，并存儲(chǔ)到原來(lái)單元內(nèi)，使得對(duì)于j=0，1，2，500，a2j=a2j-n；3）取非零向量組v(501)作為初始值，取u的所有元素都為1，每迭代一次就將新求得的v(501)存儲(chǔ)于原來(lái)的存儲(chǔ)

4、單元內(nèi)，覆蓋原來(lái)的值，執(zhí)行如下循環(huán)：對(duì)k=0，1，2，500（500為最大迭代次數(shù)，超過(guò)這個(gè)次數(shù)仍不收斂則停止計(jì)算）設(shè)置一個(gè)數(shù)comp來(lái)存儲(chǔ)前一次運(yùn)行所得的特征值，以便與第二次運(yùn)行的值相比較，求出誤差范圍，comp=c； ，但是由于A中的元素已被存為a(5,501)，對(duì)于中每一個(gè)分量和中的需對(duì)迭代公式做一下變換： 對(duì)i=0，1，2，500， 判斷，若，結(jié)束循環(huán)轉(zhuǎn)7），否則繼續(xù)循環(huán)7）將c值帶回主函數(shù)。3編輯主函數(shù)1）調(diào)用反冪法子函數(shù)，令平移量為0，得到矩陣A的按模最小特征值，并輸出A的行列式的值；2）調(diào)用冪法子函數(shù)，令其平移量為0，得到按模最大特征值；3）調(diào)用冪法子函數(shù)，令平移量為，得到的，其

5、中c是冪法子函數(shù)的帶回值；4）如果，則，；若，則，；5）設(shè)置數(shù)組s(39)，close(39)，其中close(39)用來(lái)存放A中與數(shù)s(39)最接近的特征值，執(zhí)行如下循環(huán)： 對(duì)k=0，1，38 ； 調(diào)用反冪法子函數(shù)，取平移量； ，其中q為反冪法子函數(shù)的帶回值；6）輸出和cond(A)2，并輸出所有close（39）；全部源程序：#include<stdio.h>#include<math.h>double min(double x) 編輯反冪法子程序，該子程序中也包含的矩陣的LU分解double a5501=0,u501,y501,z501,b=0,n=0;doubl

6、e m,s,sum,deta,comp,q=0;int i,j,k,t,temp,min,max;for(i=0;i<501;i+)if(i+2<501) a0i+2=-0.064;if(i+1<501) a1i+1=0.16;s=sin(0.2*(i+1);m=exp(0.1/(i+1);a2i=(1.64-0.024*(i+1)*s-0.64*m-x;if(i-1>=0) a3i-1=0.16;if(i-2>=0) a4i-2=-0.064;for(k=0;k<501;k+)min=(k+2<500)?k+2:500;for(j=k;j<=m

7、in;j+)sum=0;temp=(0>k-2)?0:k-2;max=(temp>j-2)?temp:j-2;for(t=max;t<=k-1;t+)sum=sum+ak-t+2t*at-j+2j;ak-j+2j=ak-j+2j-sum;for(i=0;i<501;i+) ui=1;if(k<500)for(i=k+1;i<=min;i+)sum=0;temp=(0>i-2)?0:i-2;max=(temp>k-2)?temp:k-2;for(t=max;t<=k-1;t+)sum=sum+ai-t+2t*at-k+2k;ai-k+2k=(

8、ai-k+2k-sum)/a2k;deta=1;if(x=0)for(i=0;i<501;i+)deta=a2i*deta;printf("deta=%.12en",deta);if(deta<=1e-12) return(x);elsefor(k=0;k<500;k+)comp=q;sum=0;for(j=0;j<501;j+)sum=uj*uj+sum;n=sqrt(sum);for(i=0;i<501;i+)yi=ui/n;z0=y0;for(i=1;i<501;i+)max=(0>i-2)?0:i-2;sum=0;for(t

9、=max;t<=i-1;t+)sum=ai-t+2t*zt+sum;zi=yi-sum;u500=z500/a2500;for(i=499;i>=0;i-)min=(i+2<500)?i+2:500;sum=0;for(t=i+1;t<=min;t+)sum=sum+ai-t+2t*ut;ui=(zi-sum)/a2i;b=0;for(i=0;i<501;i+)b=yi*ui+b;q=1/b;if(fabs(comp-q)<=1e-12) break;return(q);double max(double x) 編輯冪法子程序，帶入平移量的數(shù)值，輸出按模最大

10、特征值double a5501=0,v501,w501=0,m=0,c=0;double temp,sum,d,comp,s;int i,j,k;for(i=0;i<501;i+)if(i+2<501) a0i+2=-0.064;if(i+1<501) a1i+1=0.16;s=sin(0.2*(i+1);m=exp(0.1/(i+1);a2i=(1.64-0.024*(i+1)*s-0.64*m-x;if(i-1>=0) a3i-1=0.16;if(i-2>=0) a4i-2=-0.064;for(i=0;i<501;i+) vi=1;for(k=0;k&

11、lt;500;k+)comp=c;sum=0;for(j=0;j<501;j+)sum=vj*vj+sum;m=sqrt(sum);for(i=0;i<501;i+)wi=vi/m;for(i=0;i<501;i+)d=0;for(j=0;j<501;j+)if(i-j+2>=0&&i-j+2<5)d=wj*ai-j+2j+d;vi=d;c=0;for(j=0;j<501;j+)c=wj*vj+c;temp=fabs(c-comp);if(temp<=1e-12) break;return(c);void main() 編輯主函數(shù)

12、，調(diào)用反冪法子程序和冪法子程序，求解特征值int k;double s39,close39,conda,root,roots,root1,root501,comp;roots=min(0);root=max(0);conda=fabs(root/roots);comp=max(root)+root;if(root>0)root501=root;root1=comp;elseroot1=root;root501=comp;printf("root1=%.12e,root501=%.12e,roots=%.12en",root1,root501,roots);printf

13、("conda=%.12en",conda);for(k=0;k<39;k+)sk=root1+(k+1)*(root501-root1)/40;closek=min(sk)+sk;if(k%3=0) printf("n");printf("%.12e ",closek);程序的輸出結(jié)果：detA=2.772786141752e+118；1.070011361514e+001；cond(A)2=1.925204272920e+003；對(duì)于k=1,2,39，依次為如下值：1.018293403315e+001 9.58570742

14、5059e+000 9.172672423928e+0008.652284007898e+000 8.093483808672e+000 7.659405407692e+0007.119684648691e+000 6.611764339397e+000 6.066103226595e+0005.585101052628e+000 5.114083529812e+000 4.578872176865e+0004.096470926271e+000 3.554211215751e+000 3.041090018133e+0002.533974279002e+000 2.003230769562e

15、+000 1.503557611227e+0009.935586060075e-001 4.870426738850e-001 2.231736249575e-002 5.324174742069e-001 1.052898962694e+000 1.589445881881e+000 2.060330460171e+000 2.558075597072e+000 3.080240509307e+000 3.613620867692e+000 4.091378510451e+000 4.603035378279e+000 5.132924283898e+000 5.594906348083e+

16、000 6.080933857026e+000 6.680354092112e+000 7.293877448128e+000 7.717111714236e+000 8.225220014050e+000 8.648666065193e+000 9.254200344575e+000討論：1. 迭代初始向量中所含0的個(gè)數(shù)對(duì)計(jì)算結(jié)果有很大的影響，所含0越多，利用冪法和反冪法求得的數(shù)就越小。因?yàn)槿舫跏紩r(shí)將該分量賦值為零，在計(jì)算中，該分量仍然為0，會(huì)使疊加的結(jié)果變小，對(duì)應(yīng)于冪法得到的特征值小于實(shí)際的特征值，而對(duì)反冪法所得特征值會(huì)大于實(shí)際的特征值。2. 將冪法和反冪法分別編輯成兩個(gè)子程序，并且將LU

17、分解放在反冪法的子程序中，每次運(yùn)行冪法和反冪法之前都對(duì)二維數(shù)組a（5,501）賦值，在賦值時(shí)對(duì)a（5,501）第三行元素，即原矩陣A的對(duì)角線元素設(shè)置平移量，這樣在求解時(shí)可以直接帶入平移量，得到返回值為平移后矩陣的特征值，避免的程序的重復(fù)編寫(xiě)。3. 帶狀矩陣A存儲(chǔ)時(shí)只存儲(chǔ)帶內(nèi)元素，存為a（5,501），Aij=ai-j+2,j，應(yīng)注意在C語(yǔ)言中存儲(chǔ)是下標(biāo)的變化。4. 對(duì)原矩陣A進(jìn)行原點(diǎn)平移時(shí)，應(yīng)該先判斷平移后是否使得矩陣的特征值變?yōu)?，若平移后矩陣特征值變?yōu)?，則平移量就是A的 特征值；若平移后矩陣的特征值不是0，才能進(jìn)行反冪法的計(jì)算，以免產(chǎn)生錯(cuò)誤。5. C語(yǔ)言中利用循環(huán)疊加時(shí)要注意初始變量清零，否則會(huì)在初始變量上一直疊加，使計(jì)算結(jié)果錯(cuò)誤。6. 反冪法中利用LU分解時(shí)，Ly=b，Ux=y，為了節(jié)省存儲(chǔ)空間，會(huì)將計(jì)算所得y值存儲(chǔ)于原來(lái)b