北航數值分析大作業(yè)(一)

北京航空航天大學《數值分析A》課程計算實習《數值分析A》計算實習題目一姓名：學號：學院：2014年11月算法設計方案：A矩陣的存儲與檢索：根據題目的設定，A矩陣為上半帶寬與下半帶寬均為2的帶狀矩陣，且除主對角線元素為變量ai外，其余元素均為常數b或c，故可以將A矩陣轉存為新矩陣A[5][501]，在存入時原A矩陣的主對角元素存入新A矩陣的第三行，且各元素的列號保持不變。求解λ1、λ501、λs邏輯關系圖：求解A與數μk=λ可對A矩陣進行平移變換，對矩陣B=A+μkΙ采用反冪法求解按模最小的特征值βk，再經過反平移得到βk求矩陣A的（譜范數）條件數cond(A)2和行列式在采用反冪法求解λs的過程中要對矩陣A進行Doolittle分解（LU分解），det?(A)為U矩陣對角線元素的乘積。而根據公式cond(A)2=λmaxλmin可以得到A的條件數，其中λ程序源代碼#include<stdio.h>#include<math.h>#include<stdlib.h>#defineN501#defineM5#definee1.0e-12doublea[M][N]={0.0},u[N]={0.0};voidread_A()/*設定A矩陣函數*/{ doubleb=0.16,c=-0.064; inti; doublea_0[N]; for(i=1;i<=N;i++) a_0[i-1]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); for(i=2;i<N;i++) a[0][i]=c; for(i=1;i<N;i++) a[1][i]=b; for(i=0;i<N;i++) a[2][i]=a_0[i]; for(i=0;i<N-1;i++) a[3][i]=b; for(i=0;i<N-2;i++) a[4][i]=c;}voidread_u()/*讀取初始迭代向量*/{ inti; for(i=0;i<N;i++) u[i]=double(1231.0*rand()/1991.0);}double*unitization(doubleu[N],doubletotal)/*二范數向量單位化函數*/{ inti; doubley[N]; for(i=0;i<N;i++)y[i]=u[i]/total; return(y);}doublesecond_num(doubleu[N])/*求矩陣二范數函數*/{ inti;doublesum=0.0,total=0.0;for(i=0;i<N;i++) sum+=u[i]*u[i]; total=sqrt(sum); return(total);}doublemin(doublea,doubleb)/*最小值函數*/{ if(a<b) returna; else returnb;}doublemax(doublea,doubleb)/*最大值函數*/{ if(a<b) returnb; else returna;}double*Get_u(doublea[M][N],doubley[N])/*冪法中獲得新迭代向量*/ { inti; doubleu[N]={0.0}; u[0]=a[2][0]*y[0]+a[1][1]*y[1]+a[0][2]*y[2]; u[1]=a[3][0]*y[0]+a[2][1]*y[1]+a[1][2]*y[2]+a[0][3]*y[3]; u[N-2]=a[4][N-4]*y[N-4]+a[3][N-3]*y[N-3]+a[2][N-2]*y[N-2]+a[1][N-1]*y[N-1]; u[N-1]=a[4][N-3]*y[N-3]+a[3][N-2]*y[N-2]+a[2][N-1]*y[N-1]; for(i=2;i<N-2;i++) u[i]=a[4][i-2]*y[i-2]+a[3][i-1]*y[i-1]+a[2][i]*y[i]+a[1][i+1]*y[i+1]+a[0][i+2]*y[i+2];return(u);}voidResolve_LU(doublea[M][N])/*LU分解函數*/{ intk,i,j,t; doubletemp; for(k=1;k<=N;k++){ for(j=k;j<=min(k+2,N);j++) { temp=0; for(t=max(max(1,k-2),j-2);t<=k-1;t++) temp+=a[k-t+2][t-1]*a[t-j+2][j-1]; a[k-j+2][j-1]=a[k-j+2][j-1]-temp; }for(i=k+1;i<=min(k+2,N);i++) { temp=0; for(t=max(max(1,i-2),k-2);t<=k-1;t++) temp+=a[i-t+2][t-1]*a[t-k+2][k-1]; a[i-k+2][k-1]=(a[i-k+2][k-1]-temp)/a[2][k-1]; }}}double*back(doublea[M][N],doubleu[N])/*方程組回代函數*/{ inti,t; doubletemp=0.0; for(i=2;i<=N;i++) {temp=0.0; for(t=max(1,i-2);t<i;t++) temp+=a[i-t+2][t-1]*u[t-1]; u[i-1]-=temp; } u[N-1]=u[N-1]/a[2][N-1]; for(i=N-1;i>0;i--) { temp=0.0; for(t=i+1;t<=min(i+2,N);t++) temp+=a[i-t+2][t-1]*u[t-1]; u[i-1]=(u[i-1]-temp)/a[2][i-1]; }return(u);}doublefanmifa(doublea[M][N],doubleu[N])/*反冪法求特征值*/{inti,j,k;doubley[N]={0.0},Y[N]={0.0},x,b,B,total,E; doubleea[M][N]={0.0}; doubleex[N]={0.0}; double*p,*q,*n; B=0.0;for(i=0;;i++){b=1/B;B=0.0; total=second_num(u); p=unitization(u,total); for(j=0;j<N;j++) y[j]=*(p+j); for(j=0;j<N;j++)ex[j]=y[j]; for(j=0;j<M;j++) for(k=0;k<N;k++)ea[j][k]=a[j][k]; Resolve_LU(ea); q=back(ea,ex); for(j=0;j<N;j++) u[j]=*(q+j); total=second_num(u); n=unitization(u,total); for(j=0;j<N;j++) Y[j]=*(n+j); for(j=0;j<N;j++) B+=y[j]*u[j];x=1/B; E=abs((1/B)-b)/abs(1/B); if(E<=e)break; }x=1/B;return(x);}doublemifa(doublea[M][N],doubleu[N])/*冪法求特征值*/{inti,j;doubley[N]={0.0},Y[N]={0.0},x,b,B,total,E; double*p,*q,*n; B=0.0;for(i=0;;i++){b=B;B=0.0; total=second_num(u); p=unitization(u,total); for(j=0;j<N;j++) y[j]=*(p+j); q=Get_u(a,y); for(j=0;j<N;j++) u[j]=*(q+j); total=second_num(u); n=unitization(u,total); for(j=0;j<N;j++) Y[j]=*(n+j); for(j=0;j<N;j++) B+=y[j]*u[j];x=B; E=abs(B-b)/abs(B); if(E<=e)break; }x=B;return(x);}voiddisplace(doublea[M][N],doublex){ inti; for(i=0;i<N;i++) a[2][i]=a[2][i]+x;}voidmain()/*主程序*/{/*第一部分求解三個特征值*/ inti; doubleEig_1,Eig_s,Eig_501;doubleEig_01,Eig_02,Eig_03; doubleA[M][N],U[39],Eig[39],Eig_temp[39]; doublecond_A=0.0,DetA=1.0; read_A();read_u(); Eig_s=fanmifa(a,u);/*采用反冪法直接得到按模最小的特征值*/ read_u(); Eig_01=mifa(a,u);/*采用冪法獲得按模最大的特征值作為平移量*/ displace(a,Eig_01); read_u(); Eig_02=fanmifa(a,u); Eig_03=Eig_02-Eig_01; Eig_1=min(Eig_01,Eig_03); Eig_501=max(Eig_01,Eig_03); printf("******************數值分析作業(yè)一**********************\n"); printf("結果為：\n"); printf("（1）\n"); printf("特征值λ1=%19.11e\n",Eig_1); printf("特征值λ501=%19.11e\n",Eig_501);printf("特征值λs=%19.10e\n",Eig_s); /*第二部分求解指定的最接近特征值*/ printf("（2）\n"); for(i=0;i<39;i++) { U[i]=Eig_1+((i+1)*(Eig_501-Eig_1)/40); read_A(); displace(a,-U[i]); read_u(); Eig_temp[i]=fanmifa(a,u); Eig[i]=Eig_temp[i]+U[i]; printf("k=%d,μ%d=%19.11e,最接近的特征值為：λ%d=%19.11e\n",i+1,i+1,U[i],i+1,Eig[i]); } /*第三部分求A的譜范數條件數cond(A)和行列式detA*/printf("（3）\n"); read_A(); Resolve_LU(a); for(i=0;i<N;i++) DetA=DetA*a[2][i]; cond_A=abs(Eig_01/Eig_s); printf("cond(A)=%19.11e\n",cond_A); printf("DetA=%19.11e\n",DetA); system("pause");}程序運行結果：結果截圖如下圖1程序結果圖關于計算實習題目的討論和思考：迭代初始向量對結果的影響：在求解矩陣特征值時使用的方法為冪法和反冪法，這兩種方法都是迭代方法的一種形式，所以初始向量的選擇對計算結果產生一定的影響。當初始迭代向量中的0元素個數增大時（即x1在本程序中每一次使用冪法或者反冪法的初始迭代向量均為隨機生成的。為了進一步說明初始迭代向量的影響，這里單獨提取出冪法迭代函數部分進行測試。在這部分測試中，兩次迭代循環(huán)的初始迭代向量為隨機生成（完全不同），迭代結果如下圖顯示：圖2第一次迭代結果從兩個結果圖的對比中可以發(fā)現，在精度要求（e-12）的要求下，因為結果顯示即為12位有效數字，故雖兩次迭代循環(huán)的初始迭代向量不同，但最終收斂的結果一致。但在循環(huán)的過程中，可以發(fā)現第二次循環(huán)的初始向量中元素值為0的比例大于第一次循環(huán)，故第二次迭代的循環(huán)圖3第二次迭代結果次數比第一次迭代多，耗時也較久。接下來進一步討論初始迭代向量u中各個分量的值對迭代的影響。從迭代的結果和步數來看，當u[i]=u0(i=1,2,···,501)且u0的絕對值相對較大時，結果收斂但速度較慢，原因是u0與A矩陣元素數量級差距較大，導致了運算量增大；當u[i]=u0(i=1,2,···,501)且u0的絕對值相對