北航數(shù)值分析大作業(yè)一(純?cè)瓌?chuàng),高分版)

數(shù)值分析上機(jī)實(shí)習(xí)作業(yè)一題目：設(shè)有501×501的實(shí)對(duì)稱(chēng)矩陣A，A其中ai=1.64-0.024isin0.2λ求λ1,和λ501的值。求A的與數(shù)μk=λ求A的（普范數(shù)）條件數(shù)cond(A)和行列式detA。說(shuō)明在所用的算法中，凡是要給出精度水平的ε，都取ε=10選擇算法的時(shí)候應(yīng)使矩陣A的所有零元素都不存儲(chǔ)。打印以下內(nèi)容：(1)算法設(shè)計(jì)方案和思路。(2)全部源程序。(3)特征值λ1,λ501,λs,λik(4)發(fā)現(xiàn)的現(xiàn)象與問(wèn)題。算法設(shè)計(jì)思路和方案因?yàn)轭}目只要求求部分特征值，及最大的特征值和最小的特征值以及距離某些給定實(shí)數(shù)最近的特征值，而矩陣A是實(shí)對(duì)稱(chēng)矩陣，所有特征值均為實(shí)數(shù)，故我們可以用冪法和反冪法結(jié)合適當(dāng)?shù)钠揭屏壳蠼?。其中求解最大的特征值和最小的特征值，采用冪法并結(jié)合適當(dāng)?shù)钠揭屏壳蠼?，而求解λs，和距離μk最近的特征值采用帶平移的反冪法求解。算法分析一、冪法求解λ我們知道冪法適用于求解矩陣A的絕對(duì)值最大的特征值，設(shè)矩陣A的特征值為λ1≥λ2≥?≥|λn|，對(duì)于情況λ1>λ2，冪法是適用的，并可以求出特征值λ1，然而對(duì)于λ此時(shí)，由于λ1≤λ2≤?≤λ501B=A-λI然后對(duì)矩陣B用冪法。由線性代數(shù)知識(shí)知，矩陣B的特征值為λ1-λ,λ2-λ,?,λ501-λ，若λ=λ1，則有0=λ1-λ≤λ2-λ≤?≤λ501-λ，對(duì)B綜上，對(duì)于λ1≠-λ501，只需要對(duì)A用冪法即可得到A的絕對(duì)值最大的特征值λ，然后做平移B=A-λI，對(duì)B用冪法即可求λ由于這種情況，不能直接對(duì)矩陣A使用冪法。此時(shí)有兩種解決辦法，一種是做適當(dāng)平移B=A-pI，對(duì)用冪法迭代求出其絕對(duì)值最大的特征值λ后，則-P+λ，P+λ={λ綜合1，2，求解矩陣A的特征值λ1,λ501算法思路為設(shè)定最大迭代次數(shù)和收斂精度，然后對(duì)矩陣A使用冪法，如果收斂，說(shuō)明λ1≠-λ501，則可以得到一個(gè)特征值λ，然后做平移B=A-λI，對(duì)矩陣B使用冪法得到B的一個(gè)特征值λ'，則有{λ1,λ501}=λ,λ'+λ；若達(dá)到最大迭代次數(shù)還沒(méi)有達(dá)到收斂要求的精度，說(shuō)明λ二、反冪法求解λ設(shè)Bμ=A-μI-1，由于Bμ的特征值為1λ1-μ,1λ2-μ,?,1λ501-μ，可知以次取μ=0,u1,μ2,?,μ39時(shí)，矩陣Bμ三、求cond由定義知，condA2=A2A-12=λmλs,λm=max1≤i≤501λi，λm,λs在前面都已經(jīng)求出，可以直接用來(lái)算出condA2。而對(duì)于最后對(duì)于矩陣A的存儲(chǔ)，考慮到矩陣A的特殊性，只需定義一個(gè)一維數(shù)組用來(lái)存儲(chǔ)矩陣A的對(duì)角線上的元素，和兩個(gè)浮點(diǎn)數(shù)存儲(chǔ)b，c即可。算法實(shí)現(xiàn)流程圖冪法求解λ1反冪法求λs計(jì)算結(jié)果λ-1.070011361502E+01λ-9.935586060075E-01λ9.724634099619E+00λ-4.870426738850E-01λ-5.557910794229E-03λ2.231736249575E-02λ-1.018293403315E+01λ5.324174742069E-01λ-9.585707425068E+00λ1.052898962693E+00λ-9.172672423928E+00λ1.589445881881E+00λ-8.652284007898E+00λ2.060330460274E+00λ-8.093483808675E+00λ2.558075597073E+00λ-7.659405407692E+00λ3.080240509307E+00λ-7.119684648691E+00λ3.613620867692E+00λ-6.611764339397E+00λ4.091378510451E+00λ-6.066103226595E+00λ4.603035378279E+00λ-5.585101052628E+00λ5.132924283898E+00λ-5.114083529812E+00λ5.594906348083E+00λ-4.578872176865E+00λ6.080933857027E+00λ-4.096470926260E+00λ6.680354092112E+00λ-3.554211215751E+00λ7.293877448127E+00λ-3.041090018133E+00λ7.717111714236E+00λ-2.533970311130E+00λ8.225220014050E+00λ-2.003230769563E+00λ8.648666065193E+00λ-1.503557611227E+00λ9.254200344575E+00cond1.925204273902e+003det2.772786141752e+118初始迭代向量為：u0=(1,1,1,?,1)。分析初始向量對(duì)結(jié)果的影響分析初始向量對(duì)計(jì)算結(jié)果的影響，以計(jì)算λ1,λ501為例(初始迭代向量λ1迭代次數(shù)正確與否是否收斂u-1.070011361502e+001341是是u-2.080981085336e+000158否是μ-2.080981085338e+000184否是u-1.070011361509e+001412是是初始迭代向量λ50迭代次數(shù)正確與否是否收斂μ9.723340141234e+00010000是否u9.724634098771e+000641是是u9.724634099653e+00081否是u9.724634099652e+00073是是u9.724634098776e+0001023是是分析上表的結(jié)果可以看到，使用冪法求矩陣特征值的時(shí)候，初始向量的選取對(duì)結(jié)果的影響：一方面，會(huì)影響冪法的收斂速度，使得迭代收斂緩慢，如取初始向量為u0=(1,1,1,?,1)計(jì)算λ501時(shí)，迭代次數(shù)達(dá)到給定的最大值時(shí)仍未收斂；另一方面可能會(huì)得不到想要的結(jié)果，如取u0程序源代碼冪法#include<stdio.h>#include<math.h>#include<stdlib.h>doubleerro(doublelamd0,doublelamd1);//每步迭代結(jié)束后計(jì)算誤差的函數(shù)voiddiedai(doublea[],doubleb,doublec,double*x0,double*x1);//冪法迭代函數(shù)doubledianji(doublex[],doubley[]);//計(jì)算兩個(gè)向量?jī)?nèi)積的函數(shù)intmain(){ doublea[501]={0},x[505]={0.0},y[505]={0.0},b=0.16,c=-0.064,lamd=0,lamd0=0,lamd1=0,e1=0,d=0,*x1=y,*x0=x,*t; inti,k1,k2;//變量定義即初始化，其中一維數(shù)組x，y用于存儲(chǔ)每步迭代后的向量，指針變量x0，x1用于存取x，y的地址用于迭代計(jì)算。為了迭代格式的統(tǒng)一，對(duì)初始向量首尾各擴(kuò)充了兩個(gè)零元素 for(inti=1;i<502;++i) { a[i-1]=(1.64-i*0.024)*sin(0.2*i)-0.64*exp(0.1/i); x[i+1]=1;//賦值初始迭代向量 } for(i=0;i<10000;++i) { d=sqrt(dianji(x0,x0));//計(jì)算初始向量也即uk-1的2范數(shù) for(intj=2;j<503;j++) { x0[j]=x0[j]/d;//計(jì)算yk-1 } diedai(a,b,c,x0,x1);//根據(jù)向量yk-1計(jì)算x1，uk lamd1=dianji(x0,x1);//計(jì)算yk-1，uk的內(nèi)積，即βk e1=erro(lamd0,lamd1);//計(jì)算誤差 if(e1<1e-12)//判斷是否收斂 { lamd=lamd1; break; } else//未收斂則將uk，βk用于下次迭代 t=x0; x0=x1; x1=t; lamd0=lamd1; } k1=i-1; for(inti=0;i<501;++i) { a[i]=a[i]-lamd;//以第一次計(jì)算結(jié)果為平移量進(jìn)行冪法迭代 x[501+1]=i; } for(i=0;i<100000;++i) { d=sqrt(dianji(x0,x0)); for(intj=2;j<503;j++) { x0[j]=x0[j]/d; } diedai(a,b,c,x0,x1); lamd1=dianji(x0,x1); e1=erro(lamd0,lamd1); if(e1<1e-12) { lamd1=lamd1+lamd; k2=i;//記錄迭代步數(shù) break; } else t=x0; x0=x1; x1=t; lamd0=lamd1; } lamd1=lamd1+lamd; k2=i; printf("lamd1=%20.12e,K1=%d\nlamd501=%20.12e,K2=%d\n",lamd<lamd1?lamd:lamd1,lamd<lamd1?k1:k2,lamd>lamd1?lamd:lamd1,lamd>lamd1?k1:k2);//輸出結(jié)果 system("pause"); return0;}doubledianji(doublex[],doubley[])//定義向量?jī)?nèi)積{ doublee=0; for(inti=2;i<503;++i) { e=e+x[i]*y[i]; } return(e);}doubleerro(doublelamd0,doublelamd1)//計(jì)算兩個(gè)數(shù)之間的相對(duì)誤差{ doublee,d,l; e=fabs(lamd1-lamd0); d=fabs(lamd1); l=e/d; return(l);}voiddiedai(doublea[],doubleb,doublec,double*x0,double*x1){ for(inti=0;i<501;++i) { x1[i+2]=c*x0[i]+b*x0[i+1]+a[i]*x0[i+2]+b*x0[i+3]+c*x0[i+4];//計(jì)算uk=Ayk-1 }}反冪法#include<stdio.h>#include<math.h>#include<stdlib.h>doubleerro(doublelamd0,doublelamd1);//每步迭代結(jié)束后計(jì)算誤差的函數(shù)voiddiedai(double(*L)[501],double(*U)[3],double*x0,double*x1);//反冪法迭代函數(shù)，根據(jù)LU分解函數(shù)結(jié)果反解ukdoubledianji(doublex[],doubley[]);//計(jì)算兩個(gè)向量?jī)?nèi)積的函數(shù)voidLU(doublea[],doubleb,doublec,double(*L)[501],double(*U)[3]);//對(duì)A進(jìn)行LU分解intmain(){ doublea[501]={0},aa[501],x[501]={0.0},y[501]={0.0},l[3][501]={0},u[501][3]={0},(*L)[501],(*U)[3],b=0.16,c=-0.064,lamd0=0,lamd1=0,e1=0,d=0,*x1=y,*x0=x,*t; for(inti=1;i<502;++i)//變量定義即初始化，其中一維數(shù)組x，y用于存儲(chǔ)每步迭代后的向量，指針變量x0，x1用于存取x，y的地址用于迭代計(jì)算。其中3X501維數(shù)組l,u用于存放對(duì)a分解后的矩陣L，U，指針數(shù)組L，U用于存放l，u的元素地址 { a[i-1]=(1.64-i*0.024)*sin(0.2*i)-0.64*exp(0.1/i); x[i-1]=1;//賦值初始迭代向量 } L=l; U=u; x0=x; x1=y;LU(a,b,c,L,U);//對(duì)A進(jìn)行LU分解 for(inti=0;i<1000;++i) { d=sqrt(dianji(x0,x0));//計(jì)算初始向量也即uk-1的2范數(shù) for(intj=0;j<501;j++) { x0[j]=x0[j]/d;//計(jì)算yk-1 } diedai(L,U,x0,x1);/根據(jù)向量yk-1反解x1即uk lamd1=dianji(x0,x1);//計(jì)算yk-1，uk的內(nèi)積，即βk e1=erro(lamd0,lamd1);/計(jì)算誤差 if(e1<1e-12)//判斷是否收斂 { printf("lamds=%.12e,K=%d\n",1/lamd1,i); break; } else t=x0;//未收斂則將uk，βk用于下次迭代 x0=x1; x1=t; lamd0=lamd1; } doubledetA=1,condA,lamds=1/lamd1; for(inti=0;i<501;++i) { detA=detA*U[i][0];//根據(jù)LU分解所得的矩陣U計(jì)算detA } printf("detA=%.12e\n",detA); doublelamd=-10.700113615016,lamd501=9.724634099619,delta,miu; condA=fabs(lamd/lamds); printf("condA=%.12e\n",condA); delta=(lamd501-lamd)/40; for(intk=1;k<40;++k)//計(jì)算距離μk最近的特征值 { miu=lamd+delta*k;//計(jì)算μk for(inti=0;i<501;++i) { aa[i]=a[i]-miu;//對(duì)矩陣A進(jìn)行平移量為μk的平移 x[i]=1; y[i]=2; } LU(aa,b,c,L,U);//計(jì)算平移后的矩陣A的LU分解，并存儲(chǔ)到數(shù)組l，u中用于下面的迭代計(jì)算 lamd0=0; for(inti=0;i<1000;++i) { d=sqrt(dianji(x0,x0)); for(intj=0;j<501;j++) { x0[j]=x0[j]/d; } diedai(L,U,x0,x1); lamd1=dianji(x0,x1); e1=erro(lamd0,lamd1); if(e1<1e-12) { printf("lamdi_%d=%.12e,K=%d\n",k,1/lamd1+miu,i); break; } else t=x0; x0=x1; x1=t; lamd0=lamd1; } } system("pause"); return0;}doubledianji(doublex[],doubley[])//定義向量?jī)?nèi)積{ doublee=0; for(inti=0;i<501;++i) { e=e+x[i]*y[i]; } return(e);}doubleerro(doublelamd0,doublelamd1){ doublee,d,l; e=fabs(lamd1-lamd0); d=fabs(lamd1); l=e/d; return(l);}voidLU(doublea[],doubleb,doublec,double(*L)[501],double(*U)[3])//計(jì)算矩陣A類(lèi)型的矩陣的LU分解{ U[0][0]=a[0]; U[0][1]=b; U[0][2]=c; L[1][0]=b/U[0][0]; L[2][0]=c/U[0][0]; U[1][0]=a[1]-L[1][0]*U[0][1]; U[1][1]=b-L[1][0]*U[0][2]; U[1][2]=c;L[1][1]=(b-L[2][0]*U[0][1])/U[1][0];L[2][1]=c/U[1][0]; for(intk=2;k<499;++k) { U[k][0]=a[k]-L[2][k-2]*U[k-2][2]-L[1][k-1]*U[k-1][1]; U[k][1]