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

北京航空航天大學(xué)《數(shù)值分析A》課程計(jì)算實(shí)習(xí)《數(shù)值分析A》計(jì)算實(shí)習(xí)題目一姓名：學(xué)號(hào)：學(xué)院：2014年11月算法設(shè)計(jì)方案：A矩陣的存儲(chǔ)與檢索：根據(jù)題目的設(shè)定，A矩陣為上半帶寬與下半帶寬均為2的帶狀矩陣，且除主對(duì)角線元素為變量ai外，其余元素均為常數(shù)b或c，故可以將A矩陣轉(zhuǎn)存為新矩陣A[5][501]，在存入時(shí)原A矩陣的主對(duì)角元素存入新A矩陣的第三行，且各元素的列號(hào)保持不變。求解λ1、λ501、λs邏輯關(guān)系圖：求解A與數(shù)μk=λ可對(duì)A矩陣進(jìn)行平移變換，對(duì)矩陣B=A+μkΙ采用反冪法求解按模最小的特征值βk，再經(jīng)過(guò)反平移得到βk求矩陣A的（譜范數(shù)）條件數(shù)cond(A)2和行列式在采用反冪法求解λs的過(guò)程中要對(duì)矩陣A進(jìn)行Doolittle分解（LU分解），det?(A)為U矩陣對(duì)角線元素的乘積。而根據(jù)公式cond(A)2=λmaxλmin可以得到A的條件數(shù)，其中λ程序源代碼#include<stdio.h>#include<math.h>#include<stdlib.h>#defineN501#defineM5#definee1.0e-12doublea[M][N]={0.0},u[N]={0.0};voidread_A()/*設(shè)定A矩陣函數(shù)*/{ 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)/*二范數(shù)向量單位化函數(shù)*/{ inti; doubley[N]; for(i=0;i<N;i++)y[i]=u[i]/total; return(y);}doublesecond_num(doubleu[N])/*求矩陣二范數(shù)函數(shù)*/{ 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)/*最小值函數(shù)*/{ if(a<b) returna; else returnb;}doublemax(doublea,doubleb)/*最大值函數(shù)*/{ 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分解函數(shù)*/{ 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])/*方程組回代函數(shù)*/{ 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()/*主程序*/{/*第一部分求解三個(gè)特征值*/ 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("******************數(shù)值分析作業(yè)一**********************\n"); printf("結(jié)果為：\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的譜范數(shù)條件數(shù)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");}程序運(yùn)行結(jié)果：結(jié)果截圖如下圖1程序結(jié)果圖關(guān)于計(jì)算實(shí)習(xí)題目的討論和思考：迭代初始向量對(duì)結(jié)果的影響：在求解矩陣特征值時(shí)使用的方法為冪法和反冪法，這兩種方法都是迭代方法的一種形式，所以初始向量的選擇對(duì)計(jì)算結(jié)果產(chǎn)生一定的影響。當(dāng)初始迭代向量中的0元素個(gè)數(shù)增大時(shí)（即x1在本程序中每一次使用冪法或者反冪法的初始迭代向量均為隨機(jī)生成的。為了進(jìn)一步說(shuō)明初始迭代向量的影響，這里單獨(dú)提取出冪法迭代函數(shù)部分進(jìn)行測(cè)試。在這部分測(cè)試中，兩次迭代循環(huán)的初始迭代向量為隨機(jī)生成（完全不同），迭代結(jié)果如下圖顯示：圖2第一次迭代結(jié)果從兩個(gè)結(jié)果圖的對(duì)比中可以發(fā)現(xiàn)，在精度要求（e-12）的要求下，因?yàn)榻Y(jié)果顯示即為12位有效數(shù)字，故雖兩次迭代循環(huán)的初始迭代向量不同，但最終收斂的結(jié)果一致。但在循環(huán)的過(guò)程中，可以發(fā)現(xiàn)第二次循環(huán)的初始向量中元素值為0的比例大于第一次循環(huán)，故第二次迭代的循環(huán)圖3第二次迭代結(jié)果次數(shù)比第一次迭代多，耗時(shí)也較久。接下來(lái)進(jìn)一步討論初始迭代向量u中各個(gè)分量的值對(duì)迭代的影響。從迭代的結(jié)果和步數(shù)來(lái)看，當(dāng)u[i]=u0(i=1,2,···,501)且u0的絕對(duì)值相對(duì)較大時(shí)，結(jié)果收斂但速度較慢，原因是u0與A矩陣元素?cái)?shù)量級(jí)差距較大，導(dǎo)致了運(yùn)算量增大；當(dāng)u[i]=u0(i=1,2,···,501)且u0的絕對(duì)值相對(duì)