北航 数值分析大作业 幂法反幂法

#include <iostream> #include <math.h> using namespace std; const int N=501; double a[5][N];//数组A double u[501];//特征向量 const double E=1.0e-12;//精度水平 //////////////////创建矩阵A/////////////////// void create_matrix(){ int i; for(i=0;i<N;i++) if(i<2) a[0][i]=0; else a[0][i]=-0.064; for(i=0;i<N;i++) if(i<1) a[1][i]=0; else a[1][i]=0.16; for(i=0;i<N;i++) a[2][i]=(1.64-0.024*(i+1))*sin(0.2*(i+1))-0.64*exp(0.1/(i+1)); for(i=0;i<N;i++) if(i>=N-1) a[3][i]=0; else a[3][i]=0.16; for(i=0;i<N;i++) if(i>=N-2) a[4][i]=0; else a[4][i]=-0.064; } ////////////////////原点平移///////////////// void move (double d2){ int i; if (d2!=0.0) for(i=0;i<N;i++) a[2][i]-=d2; } /////////////////////两数取大///////////////// int max(int p,int q){ int x; x=p>q?p:q; return x; } //////////////////////两数取小/////////////////// int min(int p,int q){ int x; x=p<q?p:q; return x; } ///////////////////LU分解/////////////////////// void LU(){ int k,i,j,t; double sum; for (k=0;k<N;k++) { for (j=k;j<=min(k+2,N-1);j++) { for(sum=0.0,t=max(0,max(k-2,j-2));t<=k-1;t++) sum+=a[k-t+2][t]*a[t-j+2][j]; a[k-j+2][j]-=sum; } if (k<N-1) { for (i=k+1;i<=min(k+2,N-1);i++) { for(sum=0.0,t=max(0,max(i-2,k-2));t<=k-1;t++) sum+=a[i-t+2][t]*a[t-k+2][k]; a[i-k+2][k]=(a[i-k+2][k]-sum)/a[2][k]; } } } } ///////////////////幂法//////////////////// double mifa(double d1){ double b1,b2=0.0,sum,y[N]; int i; for (i=0;i<N;i++)//确定初始向量 u[i]=1; move(d1);//在做幂法之前确定A的平移量 do{ b1=b2; for(sum=0.0, i=0;i<N;i++) sum+=u[i]*u[i]; sum=sqrt(sum); for(i=0;i<N;i++) y[i]=u[i]/sum; for(i=0;i<N;i++) { u[i]=0; for(int j=0;j<N;j++) if (abs(j-i)<=2) { u[i]+=a[i-j+2][j]*y[j]; } } for(b2=0.0,i=0;i<N;i++) b2+=y[i]*u[i]; } while (fabs(b1-b2)/fabs(b2)>E); move(-d1);//退出还原平移量 return b2; } ///////////////////反幂法///////////////// double fanmifa(double d1){ double b1,b2=0.0,sum,y[N],m[N];; int i,t; for (i=0;i<N;i++)//确定初始向量 u[i]=1; move(d1);//在做反幂法之前确定A的平移量 LU(); //LU分解 do{ b1=b2; for(sum=0.0, i=0;i<N;i++) sum+=u[i]*u[i]; sum=sqrt(sum); for(i=0;i<N;i++) y[i]=u[i]/sum; m[0]=y[0]; for(i=1;i<N;i++) { for(sum=0.0,t=max(0,i-2);t<i;t++) sum+=a[i-t+2][t]*m[t]; m[i]=y[i]-sum; } u[N-1]=m[N-1]/a[2][N-1]; for(i=N-2;i>=0;i--) { for(sum=0.0,t=i+1;t<=min(i+2,N-1);t++) sum+=a[i-t+2][t]*u[t]; u[i]=(m[i]-sum)/a[2][i]; } for(b2=0.0,i=0;i<N;i++) b2+=y[i]*u[i]; } while (fabs(b1-b2)/fabs(b2)>E); create_matrix(); //退出还原数组A return 1.0/b2; } ///////////////////主程序/////////////////// int main(){ double d=0.0,Tzzmax,Tzzmin,Tzz1,Tzz501; // 按模最大特征值Tzzmax，按模最小特征值Tzzmin，第和个特征值Tzz1,Tzz501 double condA2; double u[39],v[39];//与u[i]接近的特征值为v[i] double sum=1.0; int i=0; cout.setf(ios::uppercase |ios::scientific ); cout.precision(12); cout.unsetf(ios::uppercase); create_matrix(); Tzzmax=mifa(d);//求按模最大特征值 //讨论按模最大特征值正负，再通过原点平移求按模最大，确定λ和λ if(Tzzmax>0) { Tzz501=Tzzmax; Tzz1=mifa(Tzzmax)+Tzzmax; } else { Tzz1=Tzzmax; Tzz501=mifa(Tzzmax)+Tzzmax; } //输出特征值λ1和λ501 cout<<"特征值λ1为："<<Tzz1<<endl; cout<<"特征值λ501为："<<Tzz501<<endl; //求按模最小特征值并输出 Tzzmin=fanmifa(d); cout<<"按模最小特征值λs为："<<Tzzmin<<endl; //求与与数μi最接近的特征值λik for (int i=0;i<39;i++) { u[i]=Tzz1+(i+1)*(Tzz501-Tzz1)/40; v[i]=fanmifa(u[i])+u[i]; cout<<"与数μ"<<i+1<<"="<<u[i]<<"最接近的特征值λi"<<i+1<<"为："<<v[i]<<endl; } //求条件数 condA2=fabs(Tzzmax/Tzzmin); cout<<"A的条件数cond(A)2为："<<condA2<<endl; //LU分解，detA=detU LU(); for(i=0;i<N;i++) sum*=a[2][i]; cout<<"A的行列式detA为："<<sum<<endl; //数组A还原 create_matrix(); return 0; }