[精文优选]北航研究生数值分析编程大作业.doc

数值分析大作业一、算法设计方案1、矩阵初始化矩阵[]501501⨯=ij a A 的下半带宽r=2,上半带宽s=2，设置矩阵[][]5011++s r C ，在矩阵C 中检索矩阵A 中的带内元素ij a 的方法是：j s j i ij c a ,1++-=。

这样所需要的存储单元数大大减少，从而极大提高了运算效率。

2、利用幂法求出5011λλ，幂法迭代格式：0111111nk k k k kk T k k k u R y u u Ay y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止。

首先对于矩阵A 利用幂法迭代求出一个λ，然后求出矩阵B ，其中I A B λ-=(I 为单位矩阵），对矩阵B 进行幂法迭代，求出λ'，之后令λλλ+'=''，比较的大小与λλ''，大者为501λ，小者为1λ。

3、利用反幂法求出ik s λλ,反幂法迭代格式：0111111nk k k k kk T k k k u R y u Au y y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止，1s k λβ=。

每迭代一次都要求解一次线性方程组1-=k k y Au ，求解过程为：（1）作分解LU A =对于n k ,...,2,1=执行[][]s k n r k k k i c c c c c n s k k k j c cc c k s ks k t k s k r i t t s t i k s k i k s k i js j t k s j r k t t s t k j s j k j s j k <+++=-=++=-=+++----=++-++-++-++----=++-++-++-∑∑);,min(,...,2,1/)(:),min(,...,1,:,1,11),,1max(,1,1,1,11),,1max(,1,1,1（2）求解y Ux b Ly ==,（数组b 先是存放原方程组右端向量，后来存放中间向量P))1,...,2,1(/)(:/:),...,3,2(:,1),min(1.1.11),1max(,1--=-===-=+++-++-+--=++-∑∑n n i c x c b x c b x n i b c b b i s t n s i i t t s t i i i ns n n ti r i t t s t i i i使用反幂法，直接可以求得矩阵按模最小的特征值s λ。

求与数)39,...,2,1(4015011=-+=k k k λλλμ最接近的特征值ik λ，对矩阵I A k μ-实行反幂法，即可求出对应的k k ik k k μλλβλ+==,/1。

4、求出A 的条件数和行列式 根据max 2()scond A λλ=，其中分子分母分别对应按模最大和最小的特征值。

det()A 的计算：由于A LU =,其中L 为下三角矩阵，且对角线元素为1，故det()1L =，所以有A LU U ==，又U 为上三角矩阵，故det()U 为对其对角线上各元素的乘积，最后可得det()det()A U =。

二、程序源代码（1）定义所需要的函数：#include<stdio.h>#include<conio.h>#include<math.h>#defineN501#defineR2#defineS2intmin(inta,intb);//求最小值intmaG(inta,intb,intc);//求最大值doubleFan_two(doubleG[N]);//计算二范数voidFenjieLU(double(GC)[N]);//解线性方程组的LU 分解过程voidSolve(double(GC)[N],doubleGb,doubleGG);//解线性方程组的求解过程 doublePowerMethod(doubleC[][N],doubleu[N],doubleP[N],doublebta,doubleD);//幂法doubleInversePowerMethod(doubleC[][N],doubleu[N],doubleP[N],double bta,doubleD);//反幂法};（2）程序的主函数，Main.cpp 代码如下：voidmain(){doubleC[R+S+1][N];doubleu[N];doubleP[N];doublemiu[39];doubleC1[R+S+1][N];doublebta=1.0;doubleNamda1,Namda501,NamdaS;doubleNamda[39];doubleCondA2;doubledetA=1.0;doubleD=1.0e-12;inti,j,k;FILEGfp;fp=fopen("Namda.tGt","w");//对数组进行初始化//inti,j;for(i=0;i<N;i++){u[i]=1;}for(i=0;i<R+S+1;i++){for(j=0;j<N;j++){if(i==0||i==4){C[i][j]=-0.064;}elseif(i==1||i==3){C[i][j]=0.16;}elseif(i==2){C[i][j]=(1.64-0.024G(j+1))Gsin(0.2G(j+1))-0.64GeGp(0.1/(j+1));}}}//幂法求Namda1//Namda1=PowerMethod(C,u,P,bta,D);printf("\n================================== ==============\n");printf("Namda1=%12.11e",Namda1);printf("\n================================================\n");//幂法求Namda501//bta=1.0;for(i=0;i<R+S+1;i++){for(j=0;j<N;j++){if(i==2)C1[i][j]= C[i][j]-Namda1;elseC1[i][j]=C[i][j];}}Namda501=algorism.PowerMethod(C1,u,P,bta,D)+Namda1;printf("\n================================== ==============\n");printf("Namda501=%12.11e",Namda501);printf("\n================================== ==============\n");//反幂法求NamdaS//bta=1.0;NamdaS=InversePowerMethod(C,u,P,bta,D);printf("\n================================== ==============\n");printf("NamdaS=%12.11e",NamdaS);printf("\n================================== ==============\n");//反幂法求Namda[k]//printf("\n================================== ==============\n");for(k=0;k<39;k++){miu[k]=Namda1+(k+1)G(Namda501-Namda1)/40.0;bta=1.0;for(i=0;i<R+S+1;i++){for(j=0;j<N;j++){if(i==2)C1[i][j]= C[i][j]-miu[k];elseC1[i][j]=C[i][j];}}Namda[k]=InversePowerMethod(C1,u,P,bta,D)+miu[k];fprintf(fp,"与%12.11e最接近的特征值为:%12.11e\n",miu[k],Namda[k]);}printf("求与miu[k]最接近的Namda[k]的计算结果已经输出到文件Namda.tGt中");printf("\n================================== ==============\n");//求A的谱范数//printf("\n================================== ==============\n");printf("A的谱范数为:%12.11e",sqrt(Namda501));printf("\n================================== ==============\n");//求A的条件数//CondA2=fabs(Namda1/NamdaS);printf("\n================================== ==============\n");printf("A的谱范数的条件数Cond(A)2为:%12.11e",CondA2);printf("\n================================== ==============\n");//求det(A)2的值//for(j=0;j<N;j++)detAG=C[2][j];printf("\n================================== ==============\n");printf("行列式A的值为:%12.11e",detA);printf("\n================================== ==============\n");fclose(fp);_getch();return;}（3）成员函数的实现intmin(inta,intb){returna<b?a:b;}intmaG(inta,intb,intc){inttemp;temp=a>b?a:b;returntemp>c?temp:c;doubleFan_two(doubleG[N]){doublesum=0.0;inti;for(i=0;i<N;i++){sum+=pow(G[i],2);}returnsqrt(sum);}voidFenjieLU(double(GC)[N]){doublesum=0;inti,j,k,t;for(k=0;k<N;k++){j=k;i=k+1;while(1){if(j==min(k+S+1,N))break;for(t=maG(0,k-R,j-S);t<=k-1;t++){sum+=C[k-t+S][t]GC[t-j+S][j];}C[k-j+S][j]=C[k-j+S][j]-sum;sum=0.0;j++;if(k==N-1)break;if(i==min(k+R+1,N))break;for(t=maG(0,i-R,k-S);t<=k-1;t++){sum+=C[i-t+S][t]GC[t-k+S][k];}C[i-k+S][k]=(C[i-k+S][k]-sum)/C[S][k];sum=0;i++;}}voidSolve(double(GC)[N],doubleGb,doubleGG){doublesum=0;inti,t;sum=0;for(i=1;i<N;i++){for(t=maG(0,i-R);t<=i-1;t++){sum+=C[i-t+S][t]Gb[t];}b[i]=b[i]-sum;sum=0;}G[N-1]=b[N-1]/C[S][N-1];for(i=N-2;i>=0;i--){for(t=i+1;t<=min(i+S,N-1);t++){sum+=C[i-t+S][t]GG[t];}G[i]=(b[i]-sum)/C[S][i];sum=0;}}doublePowerMethod(doubleC[][N],doubleu[N],doubleP[N],doublebta,double D){doubleita;doublesum=0;doubletemp=0.0;inti,j,k=0;while(fabs(bta-temp)/fabs(bta)>D){temp=bta;ita=Fan_two(u);for(i=0;i<N;i++){P[i]=u[i]/ita;}for(i=0;i<N;i++){for(j=maG(0,i-R);j<min(i+S+1,N);j++){sum+=C[i-j+S][j]GP[j];}u[i]=sum;sum=0;}for(i=0;i<N;i++){sum+=P[i]Gu[i];}bta=sum;sum=0;k++;}returnbta;}doubleInversePowerMethod(doubleC[][N],doubleu[N],doubleP[N],doublebta, doubleD){doubleTC[R+S+1][N];doubletP[N];doubleita;doublesum=0;doubletemp=0.0;inti,j,k=0;FenjieLU(C);while(abs(1/bta-1/temp)/abs(1/bta)>D){temp=bta;ita=Fan_two(u);for(i=0;i<N;i++){P[i]=u[i]/ita;}//用到临时存储数组TC[][]和tP[][]是因为函数Solve执行过程中会改变A[][]和P[][]for(i=0;i<R+S+1;i++){for(j=0;j<N;j++)TC[i][j]=C[i][j];}for(i=0;i<N;i++)tP[i]=P[i];Solve(C,P,u);for(i=0;i<R+S+1;i++){for(j=0;j<N;j++)C[i][j]=TC[i][j];}for(i=0;i<N;i++)P[i]=tP[i];for(i=0;i<N;i++){sum+=P[i]Gu[i];}bta=sum;sum=0;k++;}bta=1.0/bta;returnbta;}三、程序运行结果下图为主程序运行结果的结果输出在Namda.tGt文件中，结果如下：其中ik四、分析迭代初始向量对计算结果的影响选择不同的初始向量[]u N可能会得到不同的特征值。