北航数值分析大作业第一题幂法与反幂法.doc

数 值 分 析（B ） 大 作 业（一）姓名： 学号： 电话：1、算法设计：①求1λ、501λ和s λ的值：s λ：s λ表示矩阵的按模最小特征值，为求得s λ直接对待求矩阵A 应用反幂法即可。

1λ、501λ：若矩阵A 的特征值满足关系 1n λλ<<且1n λλ≠，要求1λ、及501λ时，可按如下方法求解： a ． 对矩阵A 用幂法，求得按模最大的特征值1m λ。

b ． 按平移量1m λ对矩阵A 进行原点平移得矩阵1m BA I λ=+，对矩阵B 用反幂法求得B 的按模最小特征值2m λ。

c ． 321m m m λλλ=-则：113min(,)m m λλλ=，13max(,)n m m λλλ=即为所求。

②求和A 的与数5011140k k λλμλ-=+最接近的特征值ik λ（k=0，1，…39）：求矩阵A 的特征值中与P 最接近的特征值的大小，采用原点平移的方法：先求矩阵 B=A-PI 对应的按模最小特征值k β，则k β+P 即为矩阵A 与P 最接近的特征值。

在本次计算实习中则是先求平移矩阵k B A I μ=-，对该矩阵应用反幂法求得s λ，则与k μ最接近的A 的特征值为：s P λ+重复以上过程39次即可求得ik λ（k=0，1，…39）的值。

③求A 的（谱范数）条件数2cond()A 和行列式det A ：在（1）中用反幂法求矩阵A 的按模最小特征值时，要用到Doolittle 分解方法，在Doolittle 分解完成后得到的两个矩阵分别为L 和U ，则A 的行列式可由U 阵求出，即：det(A)=det(U)。

求得det(A)不为0，因此A 为非奇异的实对称矩阵，则： max 2()scond A λλ=，max λ和s λ分别为模最大特征值与模最小特征值。

2、程序源代码：#include "Stdio.h"#include "Conio.h"#include "math.h"//****************************************************************************// // 在存储带状矩阵时，下面的几个量在程序中反复用到，为方便编程故把它们定义成宏.// // M ：转换后的矩阵的行数，M=R+S+1。

北航研究⽣数值分析作业第⼀题北航研究⽣数值分析作业第⼀题：⼀、算法设计⽅案1.要求计算矩阵的最⼤最⼩特征值，通过幂法求得模最⼤的特征值，进⾏⼀定判断即得所求结果；2.求解与给定数值接近的特征值，可以该数做漂移量，新数组特征值倒数的绝对值满⾜反幂法的要求，故通过反幂法即可求得；3.反幂法计算时需要⽅程求解中间过渡向量，需设计Doolite分解求解；4.|A|=|B||C|，故要求解矩阵的秩，只需将Doolite分解后的U矩阵的对⾓线相乘即为矩阵的Det。

算法编译环境：vlsual c++6.0需要编译函数：幂法，反幂法，Doolite分解及⽅程的求解⼆、源程序如下：#include#include#include#includeint Max(int value1,int value2);int Min(int value1,int value2);void Transform(double A[5][501]);double mifa(double A[5][501]);void daizhuangdoolite(double A[5][501],double x[501],double b[501]); double fanmifa(double A[5][501]); double Det(double A[5][501]);/***定义2个判断⼤⼩的函数，便于以后调⽤***/int Max(int value1,int value2){return((value1>value2)?value1:value2);}int Min(int value1,int value2){return ((value1}/*****************************************//***将矩阵值转存在⼀个数组⾥，节省空间***/void Transform(double A[5][501],double b,double c){int i=0,j=0;A[i][j]=0,A[i][j+1]=0;for(j=2;j<=500;j++)A[i][j]=c;i++;j=0;A[i][j]=0;for(j=1;j<=500;j++)A[i][j]=b;i++;for(j=0;j<=500;j++)A[i][j]=(1.64-0.024*(j+1))*sin(0.2*(j+1))-0.64*exp(0.1/(j+1)); i++;for(j=0;j<=499;j++)A[i][j]=b;A[i][j]=0;i++;for(j=0;j<=498;j++)A[i][j]=c;A[i][j]=0,A[i][j+1]=0;}/***转存结束***///⽤于求解模最⼤的特征值，幂法double mifa(double A[5][501]){int s=2,r=2,m=0,i,j;double b2,b1=0,sum,u[501],y[501];for (i=0;i<=500;i++){u[i] = 1.0;}do{sum=0;if(m!=0)b1=b2;m++;for(i=0;i<=500;i++)sum+=u[i]*u[i];for(i=0;i<=500;i++)y[i]=u[i]/sqrt(sum);for(i=0;i<=500;i++){u[i]=0;for(j=Max(i-r,0);j<=Min(i+s,500);j++)u[i]=u[i]+A[i-j+s][j]*y[j];}b2=0;for(i=0;i<=500;i++)b2=b2+y[i]*u[i];}while(fabs(b2-b1)/fabs(b2)>=exp(-12));return b2;}//带状DOOLITE分解，并且求解出⽅程组的解void daizhuangdoolite(double A[5][501],double x[501],double b[501]) { int i,j,k,t,s=2,r=2;double B[5][501],c[501];for(i=0;i<=4;i++){for(j=0;j<=500;j++)B[i][j]=A[i][j];}for(i=0;i<=500;i++)c[i]=b[i];for(k=0;k<=500;k++){for(j=k;j<=Min(k+s,500);j++){for(t=Max(0,Max(k-r,j-s));t<=k-1;t++)B[k-j+s][j]=B[k-j+s][j]-B[k-t+s][t]*B[t-j+s][j]; }for(i=k+1;i<=Min(k+r,500);i++){for(t=Max(0,Max(i-r,k-s));t<=k-1;t++)B[i-k+s][k]=B[i-k+s][k]-B[i-t+s][t]*B[t-k+s][k]; B[i-k+s][k]=B[i-k+s][k]/B[s][k];}}for(i=1;i<=500;i++)for(t=Max(0,i-r);t<=i-1;t++)c[i]=c[i]-B[i-t+s][t]*c[t];x[500]=c[500]/B[s][500];for(i=499;i>=0;i--){x[i]=c[i];for(t=i+1;t<=Min(i+s,500);t++)x[i]=x[i]-B[i-t+s][t]*x[t];x[i]=x[i]/B[s][i];}}//⽤于求解模最⼤的特征值，反幂法double fanmifa(double A[5][501]){int s=2,r=2,m=0,i;double b2,b1=0,sum=0,u[501],y[501];for (i=0;i<=500;i++){u[i] = 1.0;}do{if(m!=0)b1=b2;m++;sum=0;for(i=0;i<=500;i++)sum+=u[i]*u[i];for(i=0;i<=500;i++)y[i]=u[i]/sqrt(sum);daizhuangdoolite(A,u,y);b2=0;for(i=0;i<=500;i++)b2+=y[i]*u[i];}while(fabs(b2-b1)>=fabs(b1)*exp(-12));return 1/b2;}//⾏列式的LU分解，U的主线乘积即位矩阵的DET double Det(double A[5][501]) {int i,j,k,t,s=2,r=2;for(k=0;k<=500;k++){for(j=k;j<=Min(k+s,500);j++){for(t=Max(0,Max(k-r,j-s));t<=k-1;t++)A[k-j+s][j]=A[k-j+s][j]-A[k-t+s][t]*A[t-j+s][j];}for(i=k+1;i<=Min(k+r,500);i++){for(t=Max(0,Max(i-r,k-s));t<=k-1;t++)A[i-k+s][k]=A[i-k+s][k]-A[i-t+s][t]*A[t-k+s][k];A[i-k+s][k]=A[i-k+s][k]/A[s][k];}}double det=1;for(i=0;i<=500;i++)det*=A[s][i];return det;}void main(){double b=0.16,c=-0.064,p,q;int i,j;double A[5][501];Transform(A,b,c); //进⾏A的赋值cout.precision(12); //定义输出精度double lamda1,lamda501,lamdas;double k=mifa(A);if(k>0) //判断求得最⼤以及最⼩的特征值.如果K>0，则它为最⼤特征值值，//并以它为偏移量再⽤⼀次幂法求得新矩阵最⼤特征值，即为最⼤ //与最⼩的特征值的差{lamda501=k;for(i=0;i<=500;i++)A[2][i]=A[2][i]-k;lamda1=mifa(A)+lamda501;for(i=0;i<=500;i++)A[2][i]=A[2][i]+k;}else //如果K<=0，则它为最⼩特征值值，并以它为偏移量再⽤⼀次幂法//求得新矩阵最⼤特征值，即为最⼤与最⼩的特征值的差{lamda1=k;for(i=0;i<=500;i++)A[2][i]=A[2][i]-k;lamda501=mifa(A)+lamda1;for(i=0;i<=500;i++)A[2][i]=A[2][i]+k;}lamdas=fanmifa(A);FILE *fp=fopen("result.txt","w");fprintf(fp,"λ1=%.12e\n",lamda1);fprintf(fp,"λ501=%.12e\n",lamda501);fprintf(fp,"λs=%.12e\n\n",lamdas);fprintf(fp,"\t要求接近的值\t\t\t实际求得的特征值\n");for(i=1;i<=39;i++) //反幂法求得与给定值接近的特征值{p=lamda1+(i+1)*(lamda501-lamda1)/40;for(j=0;j<=500;j++)A[2][j]=A[2][j]-p;q=fanmifa(A)+p;for(j=0;j<=500;j++)A[2][j]=A[2][j]+p;fprintf(fp,"µ%d: %.12e λi%d: %.12e\n",i,p,i,q);}double cond=fabs(mifa(A)/fanmifa(A));double det=Det(A);fprintf(fp,"\ncond(A)=%.12e\n",cond);fprintf(fp,"\ndetA=%.12e\n",det);}三、程序运⾏结果λ1=-1.069936345952e+001λ501=9.722283648681e+000λs=-5.557989086521e-003要求接近的值实际求得的特征值µ1: -9.678281104107e+000 λi1: -9.585702058251e+000µ2: -9.167739926402e+000 λi2: -9.172672423948e+000µ3: -8.657198748697e+000 λi3: -8.652284007885e+000µ4: -8.146657570993e+000 λi4: -8.0934********e+000µ5: -7.636116393288e+000 λi5: -7.659405420574e+000µ6: -7.125575215583e+000 λi6: -7.119684646576e+000µ7: -6.615034037878e+000 λi7: -6.611764337314e+000µ8: -6.104492860173e+000 λi8: -6.0661********e+000µ9: -5.593951682468e+000 λi9: -5.585101045269e+000µ10: -5.0834********e+000 λi10: -5.114083539196e+000µ11: -4.572869327058e+000 λi11: -4.578872177367e+000µ12: -4.062328149353e+000 λi12: -4.096473385708e+000µ13: -3.551786971648e+000 λi13: -3.554211216942e+000µ14: -3.0412********e+000 λi14: -3.0410********e+000µ15: -2.530704616238e+000 λi15: -2.533970334136e+000µ16: -2.020*********e+000 λi16: -2.003230401311e+000µ17: -1.509622260828e+000 λi17: -1.503557606947e+000µ18: -9.990810831232e-001 λi18: -9.935585987809e-001µ19: -4.885399054182e-001 λi19: -4.870426734583e-001µ20: 2.200127228676e-002 λi20: 2.231736249587e-002µ21: 5.325424499917e-001 λi21: 5.324174742068e-001µ22: 1.043083627697e+000 λi22: 1.052898964020e+000µ23: 1.553624805402e+000 λi23: 1.589445977158e+000µ24: 2.064165983107e+000 λi24: 2.060330427561e+000µ25: 2.574707160812e+000 λi25: 2.558075576223e+000µ26: 3.0852********e+000 λi26: 3.080240508465e+000µ27: 3.595789516221e+000 λi27: 3.613620874136e+000µ28: 4.106330693926e+000 λi28: 4.0913********e+000µ29: 4.616871871631e+000 λi29: 4.603035354280e+000µ30: 5.127413049336e+000 λi30: 5.132924284378e+000µ31: 5.637954227041e+000 λi31: 5.594906275501e+000µ32: 6.148495404746e+000 λi32: 6.080933498348e+000µ33: 6.659036582451e+000 λi33: 6.680354121496e+000µ34: 7.169577760156e+000 λi34: 7.293878467852e+000µ35: 7.680118937861e+000 λi35: 7.717111851857e+000µ36: 8.190660115566e+000 λi36: 8.225220016407e+000µ37: 8.701201293271e+000 λi37: 8.648665837870e+000µ38: 9.211742470976e+000 λi38: 9.254200347303e+000µ39: 9.722283648681e+000 λi39: 9.724634099672e+000cond(A)=1.925042185755e+003detA=2.772786141752e+118四、分析如果初始向量选择不当,将导致迭代中X1的系数等于零.但是，由于舍⼊误差的影响,经若⼲步迭代后,.按照基向量展开时,x1的系数可能不等于零。

数值分析幂法和反幂法数值分析中的幂法和反幂法是求解矩阵最大特征值和最小特征值的常用方法。

这两种方法在许多数值计算问题中都有着广泛的应用，包括图像压缩、数据降维、谱聚类等。

幂法（Power Method）是一种迭代算法，通过不断迭代矩阵与一个向量的乘积，来逼近原矩阵的最大特征值和对应的特征向量。

其基本思想是，对于一个矩阵A和一维向量x，可以通过不断迭代计算Ax，Ax，Ax...，来使得向量x逼近最大特征值对应的特征向量。

具体的迭代过程如下：1.初始化一个向量x0（可以是单位向量或任意非零向量）2.令x1=Ax0，对向量进行归一化（即除以向量的范数）得到x13.重复步骤2，即令x2=Ax1，x3=Ax2...，直到收敛（即相邻迭代向量的差的范数小于一些阈值）为止4. 最终得到的向量xn就是A的最大特征值对应的特征向量在实际求解时，我们可以将迭代过程中的向量进行归一化，以防止数值溢出或下溢。

此外，为了提高迭代速度，我们可以选择使得xn与xn-1的内积大于0的方向作为迭代方向，这样可以使得特征值的模快速收敛到最大特征值。

幂法的收敛性是保证的，但收敛速度可能较慢，尤其是当最大特征值与其他特征值非常接近时。

此时可能需要使用一些改进的方法来加速收敛，例如Rayleigh商或位移策略。

相反，反幂法（Inverse Power Method）是求解矩阵的最小特征值和对应的特征向量的方法。

它的基本思想和幂法类似，但在每次迭代中，需要计算A和依其逆矩阵A-1的乘积。

迭代过程如下：1.初始化一个向量x0（可以是单位向量或任意非零向量）2.令x1=A-1x0，对向量进行归一化（即除以向量的范数）得到x13.重复步骤2，即令x2=A-1x1，x3=A-1x2...4. 最终得到的向量xn就是A的最小特征值对应的特征向量反幂法和幂法的区别在于迭代过程中乘以了A的逆矩阵，从而可以利用矩阵的特殊结构或性质来提高迭代速度。

同时，在实际求解时，可能需要将矩阵进行一些变换，以确保A-1存在或数值稳定性。

数值分析大作业一一、算法设计方案1、求λ1和λ501的值：思路：采用幂法求出按模最大特征值λmax，该值必为λ1或λ501，若λmax小于0，则λmax=λ1；否则λmax=λ501。

再经过原点平移，使用幂法迭代出矩阵A-λmax I的特征值，此时求出的按模最大特征值即为λ1和λ501的另一个值。

2、求λs的值：采用反幂法求出按模最小的特征值λmin即为λs，其中的方程组采用LU分解法进行求解。

3、求与μk最接近的特征值：对矩阵A采用带原点平移的反幂法求解最小特征值，其中平移量为：μk。

4、A的条件数cond(A)=| λmax/λmin|；5、A的行列式的值：先将A进行LU分解，再求U矩阵对角元素的乘积即为A 行列式的值。

二、源程序#include<iostream>#include<iomanip>#include<math.h>#define N 501#define E 1.0e-12 //定义精度常量#define r 2#define s 2using namespace std;double a[N];double cc[5][N];void init();double mifa();double fmifa();int max(int aa,int bb);int min(int aa,int bb);int max_3(int aa,int bb,int cc);void LU();void main(){double a1,a2,d1,d501=0,ds,det=1,miu[39],lamta,cond;int i,k;init();/*************求λ1和λ501********************/a1=mifa();if(a1<0)d1=a1; //若小于0则表示λ1的值elsed501=a1; //若大于0则表示λ501的值for(i=0;i<N;i++)a[i]=a[i]-a1;a2=mifa()+a1;if(a2<0)d1=a2; //若小于0则表示λ1的值elsed501=a2; //若大于0则表示λ501的值cout<<"λ1="<<setiosflags(ios::scientific)<<setprecision(12)<<d1<<"\t";cout<<"λ501="<<setiosflags(ios::scientific)<<setprecision(12)<<d501<<endl;/**************求λs*****************/init();ds=fmifa();cout<<"λs="<<setiosflags(ios::scientific)<<setprecision(12)<<ds<<endl;/**************求与μk最接近的特征值λik**************/cout<<"与μk最接近的特征值λik:"<<endl;for(k=0;k<39;k++){miu[k]=d1+(k+1)*(d501-d1)/40;init();for(i=0;i<N;i++)a[i]=a[i]-miu[k];lamta=fmifa()+miu[k];cout<<"λi"<<k+1<<"\t\t"<<setiosflags(ios::scientific)<<setprecision(12)<<lamta<<en dl;}/**************求A的条件数**************/cout<<"矩阵A的条件式";cond=abs(max(abs(d1),abs(d501))/ds);cout<<"cond="<<setiosflags(ios::scientific)<<setprecision(12)<<cond<<endl;/**************求A的行列式**************/cout<<"矩阵A的行列式";init();LU();for(i=0;i<N;i++){det*=cc[2][i];}cout<<"det="<<setiosflags(ios::scientific)<<setprecision(12)<<det<<endl;system("pause");}/**************初始化函数，给a[N]赋值*************/void init(){int i;for(i=1;i<=501;i++)a[i-1]=(1.64-0.024*i)*sin((double)(0.2*i))-0.64*exp((double)(0.1/i)); }/**************幂法求最大绝对特征值**************/double mifa(){int i,k=0;double u[N],y[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++) //控制最大迭代次数为2000{/***求y(k-1)***/double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;}/****求新的uk****/u[0]=a[0]*y[0]+b*y[1]+c*y[2];u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3]; //前两列和最后两列单独拿出来求中D间的循环求for(i=2;i<N-2;i++){u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];}u[N-2]=c*y[N-4]+b*y[N-3]+a[N-2]*y[N-2]+b*y[N-1];u[N-1]=c*y[N-3]+b*y[N-2]+a[N-1]*y[N-1];/***求beta***/double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}//cout<<"Beta"<<k<<"="<<Beta<<"\t"; 输出每次迭代的beta /***求误差***/error=abs(Beta-Beta_)/abs(Beta);if(error<=E) //若迭代误差在精度水平内则可以停止迭代{return Beta;} //控制显示位数Beta_=Beta; //第个eta的值都要保存下来，为了与后个值进行误差计算 }if(k==2000){cout<<"error"<<endl;return 0;} //若在最大迭代次数范围内都不能满足精度要求说明不收敛}/**************反幂法求最小绝对特¬征值**************/double fmifa(){int i,k,t;double u[N],y[N]={0},yy[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++){double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;yy[i]=y[i]; //用重新赋值，避免求解方程组的时候改变y的值}/****LU分解法解方程组Au=y，求新的***/LU();for(i=2;i<=N;i++){double temp_b=0;for(t=max(1,i-r);t<=i-1;t++)temp_b+=cc[i-t+s][t-1]*yy[t-1];yy[i-1]=yy[i-1]-temp_b;}u[N-1]=yy[N-1]/cc[s][N-1];for(i=N-1;i>=1;i--){double temp_u=0;for(t=i+1;t<=min(i+s,N);t++)temp_u+=cc[i-t+s][t-1]*u[t-1];u[i-1]=(yy[i-1]-temp_u)/cc[s][i-1];}double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}error=abs(Beta-Beta_)/abs(Beta);if(error<=E){return (1/Beta);}Beta_=Beta;}if(k==2000){cout<<"error"<<endl;return 0;} }/**************求两数最大值的子程序**************/int max(int aa,int bb){return(aa>bb?aa:bb);}/**************求两数最小值的子程序**************/int min(int aa,int bb){return(aa<bb?aa:bb);}/**************求三数最大值的子程序**************/int max_3(int aa,int bb,int cc){ int tt;if(aa>bb)tt=aa;else tt=bb;if(tt<cc) tt=cc;return(tt);}/**************LU分解**************/void LU(){int i,j,k,t;double b=0.16,c=-0.064;/**赋值压缩后矩阵cc[5][501]**/for(i=2;i<N;i++)cc[0][i]=c;for(i=1;i<N;i++)cc[1][i]=b;for(i=0;i<N;i++)cc[2][i]=a[i];for(i=0;i<N-1;i++)cc[3][i]=b;for(i=0;i<N-2;i++)cc[4][i]=c;for(k=1;k<=N;k++){for(j=k;j<=min(k+s,N);j++){double temp=0;for(t=max_3(1,k-r,j-s);t<=k-1;t++)temp+=cc[k-t+s][t-1]*cc[t-j+s][j-1];cc[k-j+s][j-1]=cc[k-j+s][j-1]-temp;}//if(k<500){for(i=k+1;i<=min(k+r,N);i++){double temp2=0;for(t=max_3(1,i-r,k-s);t<=k-1;t++)temp2+=cc[i-t+s][t-1]*cc[t-k+s][k-1];cc[i-k+s][k-1]=(cc[i-k+s][k-1]-temp2)/cc[s][k-1];}}}}三、程序结果。

数值分析上机作业（一）一、算法的设计方案1、幂法求解λ1、λ501幂法主要用于计算矩阵的按模最大的特征值和相应的特征向量，即对于|λ1|≥|λ2|≥.....≥|λn|可以采用幂法直接求出λ1，但在本题中λ1≤λ2≤……≤λ501，我们无法判断按模最大的特征值。

但是由矩阵A的特征值条件可知|λ1|和|λ501|之间必然有一个是最大的，通过对矩阵A使用幂法迭代一定次数后得到满足精度ε=10−12的特征值λ0，然后在对矩阵A做如下的平移：B=A-λ0I由线性代数(A-PI)x=(λ-p)x可得矩阵B的特征值为：λ1-λ0、λ2-λ0…….λ501-λ0。

对B矩阵采用幂法求出B矩阵按模最大的特征值为λ∗=λ501-λ0，所以λ501=λ∗+λ0，比较λ0与λ501的大小，若λ0>λ501则λ1=λ501，λ501=λ0；若λ0<λ501，则令t=λ501，λ1=λ0，λ501=t。

求矩阵M按模最大的特征值λ的具体算法如下：任取非零向量u0∈R nηk−1=u T(k−1)∗u k−1y k−1=u k−1ηk−1u k=Ay k−1βk=y Tk−1u k(k=1,2,3……)当|βk−βk−1||βk|≤ε=10−12时，迭终终止，并且令λ1=βk2、反幂法计算λs和λik由已知条件可知λs是矩阵A 按模最小的特征值，可以应用反幂法直接求解出λs。

使用带偏移量的反幂法求解λik，其中偏移量为μk=λ1+kλ501−λ140(k=1,2,3…39)，构造矩阵C=A-μk I，矩阵C的特征值为λik−μk，对矩阵C使用反幂法求得按模最小特征值λ0，则有λik=1λ0+μk。

求解矩阵M按模最小特征值的具体算法如下：任取非零向量u 0∈R n ηk−1= u T (k−1)∗u k−1y k−1=u k−1ηk−1 Au k =y k−1βk =y T k−1u k (k=1,2,3……)在反幂法中每一次迭代都要求解线性方程组Au k =y k−1，当K 足够大时，取λn =1βk 。

《数值分析》计算实习题目第一题：1. 算法设计方案（1）1λ，501λ和s λ的值。

1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。

2)使用反幂法求λs ，其中需要解线性方程组。

因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。

（2）与140k λλμλ-5011=+k 最接近的特征值λik 。

通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。

（3）2cond(A)和det A 。

1）1=nλλ2cond(A)，其中1λ和n λ分别是按模最大和最小特征值。

2）利用步骤（1）中分解矩阵A 得出的LU 矩阵，L 为单位下三角阵,U 为上三角阵，其中U 矩阵的主对角线元素之积即为det A 。

由于A 的元素零元素较多，为节省储存量，将A 的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A 中的元素。

2.全部源程序#include <stdio.h>#include <math.h>void init_a();//初始化Adouble get_an_element(int,int);//取A 中的元素函数double powermethod(double);//原点平移的幂法double inversepowermethod(double);//原点平移的反幂法int presolve(double);//三角LU 分解int solve(double [],double []);//解方程组int max(int,int);int min(int,int);double (*u)[502]=new double[502][502];//上三角U 数组double (*l)[502]=new double[502][502];//单位下三角L 数组double a[6][502];//矩阵Aint main(){int i,k;double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det;init_a();//初始化Alambdat1=powermethod(0);lambdat2=powermethod(lambdat1);lambda1=lambdat1<lambdat2?lambdat1:lambdat2;lambda501=lambdat1>lambdat2?lambdat1:lambdat2;presolve(0);lambdas=inversepowermethod(0);det=1;for(i=1;i<=501;i++)det=det*u[i][i];for (k=1;k<=39;k++){mu[k]=lambda1+k*(lambda501-lambda1)/40;presolve(mu[k]);lambda[k]=inversepowermethod(mu[k]);}printf("------------所有特征值如下------------\n");printf("λ=%1.11e λ=%1.11e\n",lambda1,lambda501);printf("λs=%1.11e\n",lambdas);printf("cond(A)=%1.11e\n",fabs(lambdat1/lambdas));printf("detA=%1.11e \n",det);for (k=1;k<=39;k++){printf("λi%d=%1.11e ",k,lambda[k]);if(k % 3==0) printf("\n");} delete []u;delete []l;//释放堆内存return 0;}void init_a()//初始化A{int i;for (i=3;i<=501;i++) a[1][i]=a[5][502-i]=-0.064;for (i=2;i<=501;i++) a[2][i]=a[4][502-i]=0.16;for (i=1;i<=501;i++) a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); }double get_an_element(int i,int j)//从A中节省存储量的提取元素方法{if (fabs(i-j)<=2) return a[i-j+3][j];else return 0;}double powermethod(double offset)//幂法{int i,x1;double beta=0,prebeta=-1000,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0;//设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;for (x1=1;x1<=501;x1++){u[x1]=0;for (int x2=1;x2<=501;x2++)u[x1]=u[x1]+((x1==x2)?(get_an_element(x1,x2)-offset):get_an_element(x1,x2))*y[x2];} prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);}double inversepowermethod(double offset)//反幂法{int i;double u[502],y[502];double beta=0,prebeta=0,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0; //设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;solve(u,y);prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];beta=1/beta;if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);}int presolve(double offset)//三角LU分解{int i,k,j,t;double sum;for (k=1;k<=501;k++)for (j=1;j<=501;j++){u[k][j]=l[k][j]=0;if (k==j) l[k][j]=1;} //初始化LU矩阵for (k=1;k<=501;k++){for (j=k;j<=min(k+2,501);j++){sum=0;for (t=max(1,max(k-2,j-2)) ; t<=(k-1) ; t++)sum=sum+l[k][t]*u[t][j];u[k][j]=((k==j)?(get_an_element(k,j)-offset):get_an_element(k,j))-sum;}if (k==501) continue;for (i=k+1;i<=min(k+2,501);i++){sum=0;for (t=max(1,max(i-2,k-2));t<=(k-1);t++)sum=sum+l[i][t]*u[t][k];l[i][k]=(((i==k)?(get_an_element(i,k)-offset):get_an_element(i,k))-sum)/u[k][k];}}return 0;}int solve(double x[],double b[])//解方程组{int i,t;double y[502];double sum;y[1]=b[1];for (i=2;i<=501;i++){sum=0;for (t=max(1,i-2);t<=i-1;t++)sum=sum+l[i][t]*y[t];y[i]=b[i]-sum;}x[501]=y[501]/u[501][501];for (i=500;i>=1;i--){sum=0;for (t=i+1;t<=min(i+2,501);t++)sum=sum+u[i][t]*x[t];x[i]=(y[i]-sum)/u[i][i];}return 0;}int max(int x,int y){return (x>y?x:y);}int min(int x,int y){return (x<y?x:y);}3.计算结果结果如下图所示:部分中间结果：给出了偏移量(offset)，误差(err)，迭代次数(k)4.讨论迭代初始向量的选取对计算结果的影响,并说明原因使用u[i]=1（i=1,2,...,501）作为初始向量进行迭代，可得出以上结果。

北京航空航天大学数值分析大作业一学院名称自动化专业方向控制工程学号ZY*******学生姓名许阳教师孙玉泉日期2021 年11月26 日设有501501⨯的实对称矩阵A ，⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=5011A a b c b c c b c b a其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1.0-==⋅⋅⋅=--=c b i e i i a ii 。

矩阵A 的特征值为)501,,2,1(⋅⋅⋅=i i λ，并且有||min ||,501150121i i s λλλλλ≤≤=≤⋅⋅⋅≤≤1λ，501λ和s λ的值。

A 的与数4015011λλλμ-+=kk 最接近的特征值)39,,2,1(⋅⋅⋅=k k i λ。

A 的(谱范数)条件数2)A (cond 和行列式detA 。

一 方案设计1 求1λ，501λ和s λ的值。

s λ为按模最小特征值，||min ||5011i i s λλ≤≤=。

可使用反幂法求得。

1λ，501λ分别为最大特征值及最小特征值。

可使用幂法求出按模最大特征值，如结果为正，即为501λ，结果为负，那么为1λ。

使用位移的方式求得另一特征值即可。

2 求A 的与数4015011λλλμ-+=kk 最接近的特征值)39,...,2,1(=k k i λ。

题目可看成求以k μ为偏移量后，按模最小的特征值。

即以k μ为偏移量做位移，使用反幂法求出按模最小特征值后，加上k μ，即为所求。

3 求A 的(谱范数)条件数2)(A cond 和行列式detA 。

矩阵A 为非奇异对称矩阵，可知，||)(min max2λλ=A cond(1-1)其中m ax λ为按模最大特征值，min λ为按模最小特征值。

detA 可由LU 分解得到。

因LU 均为三角阵，那么其主对角线乘积即为A 的行列式。

二 算法实现1 幂法使用如下迭代格式：⎪⎪⎩⎪⎪⎨⎧⋅===⋅⋅⋅=------||max |)|sgn(max ||max /),,(111111)0()0(10k k k k k k k k Tn u u Ay u u u y u u u β任取非零向量 (2-1)终止迭代的控制理论使用εβββ≤--||/||1k k k ， 实际使用εβββ≤--||/||||||1k k k(2-2)由于不保存A 矩阵中的零元素，只保存主对角元素a[501]及b,c 值。

数值分析—计算实习作业一学院：机械工程学院专业：材料加工工程姓名：暴一品学号：SY12071342012-10-29一、算法设计方案观察矩阵A ，结构为带状，且与主对角线相邻的两个带的值b 和c 都是常数。

从而可以用带原点平移的幂法或反幂法计算λ1和λ501。

所以算法的设计方案如下：1.求按模最大的特征值，并记为max_eigenvalue ，算法如下所示⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=======------≤≤-),2,1()sgn(),,(/max ),,()(1)()(11)1(11)1(1)1()0()0(10ΛΛΛk h h h h Ay u h u y h h h h u k r k r k Tk nk k kk r k k k j nj k rTn β任取非零向量2.平移矩阵得到A ’=A-max_eigenvalueI ，再次用幂法，这次求出的A ’的按模大的特征值pymax_eigenvalue 就是与步骤1求出的特征值相差最大的特征值。

即两者一个为最大的特征值，另一个为最小的特征值。

3.根据max_eigenvalue 和pymax_eigenvalue 的正负性，直接确定λ1，和λ501。

4.对原矩阵A 用反幂法，求出其按模最小的特征值，记为s_eigenvalue ，此即λs 。

⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=====∈--------),2,1(/111111110Λk u y y Au u y u u R u k T k k k k k k k k Tk k n βηη任取非零向量在反幂法的求解过程中，每迭代一次都要求满足解线性方程组Auk=yk-1。

本题中矩阵A 上半带宽为2，下半带宽也为2 。

故选择采用三角分解法求解方程组：先将原矩阵改写成5行501列的矩阵C （不存储A 的0元素） A 的带内元素aij=c 中的元素ci-j+3。

再对C 矩阵做LU 分解。

对于k=1,2，…，n ，执行∑---=+-+-+-+--=1)2,2,1max(,3,3,3,3:k j k t jj t t t k j j k j j k ccc c [j=k,k+1,…,min(k+2,n)]kk s k r i t k k t t t i k k i k k i c ccc c ,31),,1max(,3,3,3,3/)(:∑---=+-+-+-+--=[i=k+1,k+2,…,min(k+2,n);k<n]求解Lx=b ，Uuk=x （数组b 先是存放原方程组右端向量yk-1，后来存放中间向量x ）∑--=+--=1),1max(,3:i r i t tt t i i i bcb b （i=2,3,…,n ）nn kn c b u ,3/:=in i i t kt tt i i ki c u cb u ,3),2min(1,3/)(:∑++=+--= (i=n-1,n-2, (1)5.对k=1，2，……39执行：先根据题中给出的公式算出μk ，再将矩阵平移A ”=A-μk ，对矩阵A ”运用反幂法（线性方程组的解法同上），就可以求出与μk 最接近的特征值λik ，保存在数组py_eigenvalue 中。

数值分析幂法和反幂法数值分析中，幂法（Power method）和反幂法（Inverse Power method）是求解矩阵的特征值和特征向量的两种常用方法。

它们都是通过迭代过程逼近特征值和特征向量。

1.幂法：幂法是求解矩阵的最大特征值和对应的特征向量的一种迭代方法。

幂法的原理是通过迭代过程，将一个任意选择的初始向量不断与矩阵相乘，使其逼近对应最大特征值的特征向量。

幂法的迭代公式为：$x^{(k+1)} = \frac{Ax^{(k)}}{\，Ax^{(k)}\，}$幂法的迭代过程是不断对向量进行归一化，使其逐渐逼近最大特征值对应的特征向量。

当迭代次数足够多时，可以得到非常接近最大特征值的估计。

2.反幂法：反幂法是幂法的一种变形，用于求解矩阵的最小特征值和对应的特征向量。

反幂法的原理是通过迭代过程，将一个任意选择的初始向量不断与矩阵的逆相乘，使其逼近对应最小特征值的特征向量。

反幂法的迭代公式为：$x^{(k+1)} = \frac{A^{-1}x^{(k)}}{\，A^{-1}x^{(k)}\，}$反幂法的迭代过程同样是不断对向量进行归一化，使其逐渐逼近最小特征值对应的特征向量。

当迭代次数足够多时，可以得到非常接近最小特征值的估计。

3.收敛性分析：幂法和反幂法的收敛性分析与矩阵的特征值分布有关。

对于幂法而言，如果矩阵$A$的最大特征值是唯一的，并且其他特征值的绝对值小于最大特征值的绝对值，那么幂法是收敛的，而且收敛速度是指数级的。

对于反幂法而言，如果矩阵$A$的最小特征值是唯一的，并且其他特征值的绝对值大于最小特征值的绝对值，那么反幂法是收敛的，而且同样是指数级的收敛速度。

4.实际应用：幂法和反幂法在实际中广泛应用于各个领域，例如物理、工程、计算机科学等。

比如在结构力学中，幂法可以用来求解结构的自振频率和相应的振型；在电力系统中，反幂法可以用来求解电力系统决定性特征值，例如功率稳定性的最小特征值。

数值分析大作业一、算法设计方案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 先是存放原方程组右端向量，后来存放中间向量y))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 λ。

一、算法的设计方案：（一）各所求值得计算方法1、最大特征值λ501，最小特征值λ1，按模最小特征值λs的计算方法首先使用一次幂法运算可以得到矩阵的按模最大的特征值λ，λ必为矩阵A的最大或最小特征值，先不做判断。

对原矩阵A进行一次移项，即（A-λI），在进行一次幂法运算，可以得到另一个按模最大特征值λ0。

比较λ和λ的大小，较大的即为λ501，较小的即为λ1。

对矩阵A进行一次反幂法运算，即可得到按模最小特征值λs。

2、A与μk 值最接近的特征值λik的计算方法首先计算出k所对应的μk 值，对原矩阵A进行一次移项，即（A-μkI），得到一个新的矩阵，对新矩阵进行一次反幂法运算，即可得到一个按模最小特征值λi 。

则原矩阵A与μk值最接近的特征值λik=λi+μk。

3、A的（谱范数）条件数cond（A）2的计算方法其中错误!未找到引用源。

矩阵A的按模最大和按模最小特征值。

（二）程序编写思路。

由于算法要求A的零元素不存储，矩阵A本身为带状矩阵，所以本题的赋值，LU分解，反幂法运算过程中，均应采用Doolittle分解法求解带状线性方程组的算法思路。

幂法、反幂法和LU分解均是多次使用，应编写子程序进行反复调用。

二、源程序：#include<stdio.h>#include<iostream>#include<stdlib.h>#include<math.h>#include<float.h>#include<iomanip> /*头文件*//*定义全局变量*/#define N 502 /*取N为502，可实现从1到501的存储，省去角标变换的麻烦*/ #define epsilon 1.0e-12 /*定义精度*/#define r 2 /*r，s为带状矩阵的半带宽，本题所给矩阵二者都是2*/ #define s 2double c[6][N]; /*定义矩阵存储压缩后的带状矩阵*/double fuzhi(); /*赋值函数*/void LUfenjie(); /*LU分解程序*/int max(int a,int b); /*求两个数字中较大值*/int min(int a,int b); /*求两个数字中较小值*/double mifa(); /*幂法计算程序*/double fanmifa(); /*反幂法计算程序*/double fuzhi() /*赋值程序，按行赋值，行从1到5，列从1到501，存储空间的第一行第一列不使用，角标可以与矩阵一一对应，方便书写程序*/{int i,j;i=1;for(j=3;j<N;j++){c[i][j]=-0.064;}i=2;for(j=2;j<N;j++){c[i][j]=0.16;}i=3;for(j=1;j<N;j++){c[i][j]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);}i=4;for(j=1;j<N-1;j++){c[i][j]=0.16;}i=5;for(j=1;j<N-2;j++){c[i][j]=-0.064;}return(c[i][j]);}int max(int a,int b){ return((a>b)?a:b);}int min(int a,int b){ return((a<b)?a:b);}void LUfenjie() /*LU分解程序，采用的是带状矩阵压缩存储后的LU分解法*/{double temp;int i,j,k,t;for(k=1;k<N;k++){ for(j=k;j<=min(k+s,N-1);j++){temp=0;for(t=max(1,max(k-r,j-s));t<=(k-1);t++){temp=temp+c[k-t+s+1][t]*c[t-j+s+1][j];}c[k-j+s+1][j]=c[k-j+s+1][j]-temp;}for(i=k+1;i<=min(k+r,N-1);i++){temp=0;for(t=max(1,max(i-r,k-s));t<=(k-1);t++){temp=temp+c[i-t+s+1][t]*c[t-k+s+1][k];}c[i-k+s+1][k]=(c[i-k+s+1][k]-temp)/c[s+1][k];}}}double mifa() /*幂法计算程序*/ {double u0[N],u1[N];double temp,Lu,beta=0,beta0;int i,j;for(i=1;i<N;i++) /*选取迭代初始向量*/{u0[i]=1;}do{beta0=beta;temp=0;for(i=1;i<N;i++){temp=temp+u0[i]*u0[i]; }Lu=sqrt(temp);for(i=1;i<N;i++){u1[i]=u0[i]/Lu;}for(i=1;i<N;i++){temp=0;for(j=max(i-1,1);j<=min(i+2,N-1);j++){temp=temp+c[i-j+s+1][j]*u1[j]; }u0[i]=temp;} //新的u0temp=0;for(i=1;i<N;i++){temp=temp+u1[i]*u0[i]; }beta=temp;}while(fabs(beta-beta0)/fabs(beta)>=epsilon); /*迭代运行条件判断*/return(beta);}double fanmifa() /*反幂法计算程序*/{double u0[N],u1[N],u2[N];double temp,Lu,beta=0,beta0;int i,t;for(i=1;i<N;i++) /*选取迭代初始向量*/{u0[i]=1;}do{beta0=beta;temp=0;for(i=1;i<N;i++){temp=temp+u0[i]*u0[i]; }Lu=sqrt(temp);for(i=1;i<N;i++){u1[i]=u0[i]/Lu;u2[i]=u1[i];}fuzhi();LUfenjie();/*带状矩阵压缩存储并进行LU分解后，求解线性方程组得到迭代向量u k，即程序中的u0*/for(i=2;i<N;i++){ temp=0;for(t=max(1,i-r);t<=(i-1);t++){temp=temp+c[i-t+s+1][t]*u2[t];}u2[i]=u2[i]-temp;}u0[N-1]=u2[N-1]/c[s+1][N-1];for(i=N-2;i>=1;i--){ temp=0;for(t=i+1;t<=min(i+s,N-1);t++){temp=temp+c[i-t+s+1][t]*u0[t];}u0[i]=(u2[i]-temp)/c[s+1][i];}temp=0;for(i=1;i<N;i++){temp=temp+u1[i]*u0[i]; }beta=temp;beta=1/beta; /*beta即为所求特征值，可直接返回*/}while(fabs(beta-beta0)/fabs(beta)>=epsilon); /*迭代运行条件判断*/return(beta);}void main(){double u[40]; /*定义数组，存放k值运算得到的μk值*/double lambda1,lambda501,lambdak,a,b,d,cond,det;int i,j,k;fuzhi();a=mifa(); /*幂法计算按模最大值*/fuzhi();d=fanmifa(); /*反幂法计算按模最小值*/fuzhi();for(j=1;j<N;j++){c[3][j]=c[3][j]-a;}b=mifa()+a; /*移项后幂法计算按模最大值*/if(a>b) /*比较两个按模最大值大小，并相应输出最大特征值λ501和最小特征值λ1*/ {lambda1=b;lambda501=a;printf("矩阵A最小的特征值lambda1=%13.11e\n",lambda1);printf("矩阵A最大的特征值lambda501=%13.11e\n",lambda501);}else{lambda1=a;lambda501=b;printf("矩阵A最小的特征值lambda1=%13.11e\n",lambda1);printf("矩阵A最大的特征值lambda501=%13.11e\n",lambda501);}printf("矩阵A按模最小特征值lambdas=%13.11e\n",d); /*输出按模最小特征值λs*/for(k=1;k<40;k++) /*对每一个进行移项反幂法运算，求出最接近μk的特征值并输出*/ {u[k]=(lambda501-lambda1)*k/40+lambda1;fuzhi();for(j=1;j<N;j++){c[3][j]=c[3][j]-u[k];}lambdak=fanmifa()+u[k];i=k;printf("矩阵A最接近uk特征值lambdak%d=%13.11e\n",i,lambdak);}cond=fabs(a/d);printf("A的条件数=%13.11e\n",cond); /*计算A条件数并输出*/fuzhi(); /*计算A的行列式值并输出*/LUfenjie();det=1;for(i=1;i<N;i++){det=det*c[3][i];}printf("行列式的值detA=%13.11e\n",det);}三、程序的运行结果：四、初始向量的选取对计算结果的影响：（一）选取形式不变，数值变换1、取u0为[0.5,0.5………..0.5]，运行结果如下：2、取u0为[50,50………..50]，运行结果如下：从运行结果来看，此类初始向量的选取对结果不会产生影响，即使选成0，结果也不变化。

数值分析—计算实习作业一学院：17系专业：精密仪器及机械姓名：张大军学号：DY14171142014-11-11数值分析计算实现第一题报告一、算法方案算法方案如图1所示。

（此算法设计实现完全由本人独立完成）图1算法方案流程图二、全部源程序全部源程序如下所示#include <iostream.h>#include <iomanip.h>#include <math.h>int main(){double a[501];double vv[5][501];double d=0;double r[3];double uu;int i,k;double mifayunsuan(double *a,double weiyi);double fanmifayunsuan(double *a,double weiyi);void yasuo(double *A,double (*C)[501]);void LUfenjie(double (*C)[501]);//赋值语句for(i=1;i<=501;i++){a[i-1]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);}//程序一：使用幂方法求绝对值最大的特征值r[0]=mifayunsuan(a,d);//程序二：使用幂方法求求平移λ[0]后绝对值最大的λ，得到原矩阵中与最大特征值相距最远的特征值d=r[0];r[1]=mifayunsuan(a,d);//比较λ与λ-λ[0]的大小，由已知得if(r[0]>r[1]){d=r[0];r[0]=r[1];r[1]=d;}//程序三：使用反幂法求λr[2]=fanmifayunsuan(a,0);cout<<setiosflags(ios::right);cout<<"λ["<<1<<"]="<<setiosflags(ios::scientific)<<setprecision(12)<<r[0]<<endl;cout<<"λ["<<501<<"]="<<setiosflags(ios::scientific)<<setprecision(12)<<r[1]<<endl;cout<<"λ[s]="<<setiosflags(ios::scientific)<<setprecision(12)<<r[2]<<endl;//程序四：求A的与数u最接近的特征值for(k=1;k<40;k++){uu=r[0]+k*(r[1]-r[0])/40;cout<<"最接近u["<<k<<"]"<<"的特征值为"<<setiosflags(ios::scientific)<<setprecision(12)<<fanmifayunsuan(a,uu)<<endl;}//程序五：谱范数的条件数是绝对值最大的特征值除以绝对值最小的特征值的绝对值cout<<"cond(A)2="<<fabs(r[0]/r[2])<<endl;//程序六：A的行列式的值就是A分解成LU之U的对角线的乘积yasuo(a,vv);LUfenjie(vv);uu=1;for(i=0;i<501;i++){uu=uu*vv[2][i];}cout<<"Det(A)="<<uu<<endl;return 1;}double mifayunsuan(double *a,double weiyi){int i,k;double b=0.16;double c=-0.064;double ee,w,v1,v2,mm,sum;double u[501];double y[505]={0};for(i=0;i<501;i++)u[i]=1;//给u赋初值if (weiyi!=0){for (i=0;i<501;i++)a[i]-=weiyi;}ee=1;k=0;//使得初始计算时进入循环语句while(ee>1e-12){mm=0;for(i=0;i<501;i++){mm=mm+u[i]*u[i];}w=sqrt(mm);for(i=0;i<501;i++){y[i+2]=u[i]/w;//注意此处编程与书上不同，之后会解释它的巧妙之处1 }for(i=0;i<501;i++){u[i]=c*y[i]+b*y[i+1]+a[i]*y[i+2]+b*y[i+3]+c*y[i+4];//1显然巧妙之处凸显出来}sum=0;for(i=0;i<501;i++){sum+=y[i+2]*u[i];}v1=v2;v2=sum;//去除特殊情况，减少漏洞if(k==0){k++;}else{ee=fabs(v2-v1)/fabs(v2);}}if (weiyi!=0){for (i=0;i<501;i++)a[i]+=weiyi;}//还原A矩阵return (v2+weiyi);}double fanmifayunsuan(double *a,double weiyi){int i,k;double b=0.16;double c=-0.064;double ee,w,v1,v2,mm,sum;double u[501];double y[501];double C[5][501];void yasuo(double *A,double (*C)[501]);void LUfenjie(double (*C)[501]);void qiuU(double (*C)[501],double *y,double *u);//把A阵压缩到C阵中for(i=0;i<501;i++)u[i]=1;//给u赋初值if (weiyi!=0){for (i=0;i<501;i++)a[i]-=weiyi;}yasuo(a,C);LUfenjie(C);ee=1;k=0; //使得初始计算时进入循环语句while(ee>1e-12){mm=0;for(i=0;i<501;i++){mm=mm+u[i]*u[i];}w=sqrt(mm);for(i=0;i<501;i++){y[i]=u[i]/w;}qiuU(C,y,u);sum=0;for(i=0;i<501;i++){sum+=y[i]*u[i];}v1=v2;v2=sum;//去除特殊情况，减少漏洞if(k==0){k++;}else{ee=fabs(1/v2-1/v1)/fabs(1/v2);}}if (weiyi!=0){for (i=0;i<501;i++)a[i]+=weiyi;}//还原A矩阵return (1/v2+weiyi);}void yasuo(double *A,double (*C)[501]){double b=0.16;double c=-0.064;int i;for(i=0;i<501;i++){C[0][i]=c;C[1][i]=b;C[2][i]=A[i];C[3][i]=b;C[4][i]=c;}}void LUfenjie(double (*C)[501]){int k,t,j;int r=2,s=2;double sum;int minn(int ,int );int maxx(int ,int );for(k=0;k<501;k++){for(j=k;j<=minn(k+s,501-1);j++){if(k==0)sum=0;else{sum=0;for(t=maxx(k-r,j-s);t<k;t++){sum=sum+C[k-t+s][t]*C[t-j+s][j];}}C[k-j+s][j]=C[k-j+s][j]-sum;}for(j=k+1;j<=minn(k+r,501-1);j++){if(k<501-1){if(k==0)sum=0;else{sum=0;for(t=maxx(j-r,k-s);t<k;t++){sum=sum+C[j-t+s][t]*C[t-k+s][k];}}C[j-k+s][k]=(C[j-k+s][k]-sum)/C[s][k];}}}}void qiuU(double (*C)[501],double *y,double *u){int i,t;double b[501];double sum;int r=2,s=2;int minn(int ,int );int maxx(int ,int );for(i=0;i<501;i++){b[i]=y[i];}for(i=1;i<501;i++){sum=0;for(t=maxx(0,i-r);t<i;t++){sum=sum+C[i-t+s][t]*b[t];}b[i]=b[i]-sum;}u[500]=b[500]/C[s][500];for(i=501-2;i>=0;i--){sum=0;for(t=i+1;t<=minn(i+s,500);t++){sum=sum+C[i-t+s][t]*u[t];}u[i]=(b[i]-sum)/C[s][i];}}int minn(int x,int y){int min;if(x>y)min=y;elsemin=x;return min;}int maxx(int b,int c){int max;if(b>c){if(b>0)max=b;elsemax=0;}else{if(c>0)max=c;elsemax=0;}return max;}三、特征值以及的值λ[1]=-1.070011361502e+001 λ[501]=9.724634098777e+000λ[s]=-5.557910794230e-003最接近u[1]的特征值为-1.018293403315e+001最接近u[2]的特征值为-9.585707425068e+000最接近u[3]的特征值为-9.172672423928e+000最接近u[4]的特征值为-8.652284007898e+000最接近u[5]的特征值为-8.0934********e+000最接近u[6]的特征值为-7.659405407692e+000最接近u[7]的特征值为-7.119684648691e+000最接近u[8]的特征值为-6.611764339397e+000最接近u[9]的特征值为-6.0661********e+000最接近u[10]的特征值为-5.585101052628e+000最接近u[11]的特征值为-5.114083529812e+000最接近u[12]的特征值为-4.578872176865e+000最接近u[13]的特征值为-4.096470926260e+000最接近u[14]的特征值为-3.554211215751e+000最接近u[15]的特征值为-3.0410********e+000最接近u[16]的特征值为-2.533970311130e+000最接近u[17]的特征值为-2.003230769563e+000最接近u[18]的特征值为-1.503557611227e+000最接近u[19]的特征值为-9.935586060075e-001最接近u[20]的特征值为-4.870426738850e-001最接近u[21]的特征值为2.231736249575e-002最接近u[22]的特征值为5.324174742069e-001最接近u[23]的特征值为1.052898962693e+000最接近u[24]的特征值为1.589445881881e+000最接近u[25]的特征值为2.060330460274e+000最接近u[26]的特征值为2.558075597073e+000最接近u[27]的特征值为3.080240509307e+000最接近u[28]的特征值为3.613620867692e+000最接近u[29]的特征值为4.0913********e+000最接近u[30]的特征值为4.603035378279e+000最接近u[31]的特征值为5.132924283898e+000最接近u[32]的特征值为5.594906348083e+000最接近u[33]的特征值为6.080933857027e+000最接近u[34]的特征值为6.680354092112e+000最接近u[35]的特征值为7.293877448127e+000最接近u[36]的特征值为7.717111714236e+000最接近u[37]的特征值为8.225220014050e+000最接近u[38]的特征值为8.648666065193e+000最接近u[39]的特征值为9.254200344575e+000cond(A)2=1.925204273902e+003 Det(A)=2.772786141752e+118四、现象讨论在大作业的程序设计过程当中，初始向量的赋值我顺其自然的设为第一个分量为1，其它分量为0的向量，计算结果与参考答案存在很大差别，计算结果对比如下图2所示（左侧为正确结果，右侧为错误结果），导致了我花了很多的时间去检查程序算法。

一、算法的设计方案1、求矩阵最大特征值，最小特征值与按模最小特征值的方法首先用幂法求出矩阵A 的一个特征值λ，则其必为最大特征值与最小特征值二者其一，之后对矩阵A 进行一次移项，即A-λI ，然后再次用幂法求出另一个按模最大特征值，再比较这两个值的大小，则较大的为矩阵A 的最大特征值，较小的为矩阵A 的最小特征值。

用反幂法可以求得矩阵的按模最小特征值λs 2、求矩阵A 与k μ最接近的特征值k i λ可以先对矩阵A 进行移项，即A-k μI ，对这个移项后的矩阵用反幂法求出按模最小的特征值，然后再加上k μ，就求出所要求的k i λ。

3、求矩阵A 的条件数cond(A)2和行列式detA由于矩阵A 是非奇异的实对称矩阵，所以可以用以下公式方便地求出矩阵A 的条件数cond(A)2=sλλ501对于矩阵A 行列式的求法也比较简单。

由于在用反幂法的过程中对A 进行了Doolittle LU 分解，所以detA=detL*detU ，而detL=1，detU 可以用对角线元素相乘方便地算出，所以detA 就是U 阵对角线元素的乘积。

4、几点说明由于A 中的零元素都不存储，所以在存储矩阵的时候采用书上26页的压缩存储方式。

在反幂法中采用LU 分解求解带状线性方程组的算法来求解每一次迭代的方程组，由于每一次方程左边的系数都相同，所以只要进行一次LU 分解即可。

因为幂法，反幂法，LU 分解，求最大值与最小值在程序编写的过程中多次用到，所以这几项作为子函数单独进行编写。

二、源程序如下：#include "stdio.h"#include "math.h"# define s 2# define r 2# define N 501double c[5][N]={0};double lameda[40];double max(double x,double y);double min(double x,double y);double mifa();double fanmifa();void LUfenjie();void main(){int i=0,j=0;/*============对数组进行赋值==============*/ for(j=3;j<=N;j++)c[0][j-1]=-0.064;for(j=2;j<=N;j++)c[1][j-1]=0.16;for(j=1;j<=N;j++)c[2][j-1]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);for(j=1;j<=N-1;j++)c[3][j-1]=0.16;for(j=1;j<=N-2;j++)c[4][j-1]=-0.064;/*========幂法求最大和最小特征值==============*/ double a=mifa();for(j=1;j<=N;j++)c[2][j-1]-=a;double b=mifa()+a;double lameda501=max(a,b);double lameda1=min(a,b);printf("矩阵A最大的特征值=%13.11e\n",lameda501);printf("矩阵A最小的特征值=%13.11e\n",lameda1);/*========反幂法求绝对值最小特征值===========*/for(j=1;j<=N;j++)c[2][j-1]+=a;double lamedas=fanmifa();printf("矩阵A按模最小的特征值=%13.11e\n",lamedas);/*========求条件数和行列式的值===========*/double detA=1;for(j=1;j<=N;j++)detA*=c[2][j-1];printf("矩阵A的行列式=%13.11e\n",detA);double condA=fabs(lameda501/lamedas);printf("矩阵A的条件数=%13.11e\n",condA);/*========反幂法求与uk最接近的特征值========*/for(int k=1;k<40;k++){ for(j=3;j<=N;j++)c[0][j-1]=-0.064;for(j=2;j<=N;j++)c[1][j-1]=0.16;for(j=1;j<=N;j++)c[2][j-1]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);for(j=1;j<=N-1;j++)c[3][j-1]=0.16;for(j=1;j<=N-2;j++)c[4][j-1]=-0.064;for(j=1;j<=N;j++)c[2][j-1]-=(lameda1+k*(lameda501-lameda1)/40);lameda[k]=fanmifa()+(lameda1+k*(lameda501-lameda1)/40);printf("矩阵A最接近u%d的特征值=%13.11e\n",k,lameda[k]);}}double max(double x,double y) //求两数中的最大值{double z;z=x>y ? x:y;return(z);}double min(double x,double y) //求两数中的最小值{double z;z=x<y ? x:y;return(z);}double mifa() //幂法求按模最大特征值{double u[N]={0};double sum=0;double zero=0;double y[N]={0};double b1=0,b2=0;int i=0,j=0;for(i=0;i<N;i++)u[i]=1;do{ sum=0;for(i=0;i<N;i++)sum+=u[i]*u[i];for(i=0;i<N;i++)y[i]=u[i]/sqrt(sum);for(i=0;i<N;i++){ u[i]=0;for(j=max(0,i-2);j<=min(i+2,N-1);j++){u[i]+=c[i-j+s][j]*y[j];}}b2=b1;zero=0;for(i=0;i<N;i++)zero+=y[i]*u[i];b1=zero;}while((fabs(b1-b2)/fabs(b1))>1e-12);return(b1);}double fanmifa() //反幂法求按模最小特征值{double u[N]={0};double sum=0;double zero=0;double y[N]={0};double b[N]={0};double b1=0,b2=0;int i=0,j=0,t=0;for(i=0;i<N;i++)u[i]=1;LUfenjie();do{ sum=0;for(i=0;i<N;i++)sum+=u[i]*u[i];for(i=0;i<N;i++){b[i]=u[i]/sqrt(sum);y[i]=b[i];}for(i=1;i<N;i++)for(t=max(0,i-r);t<i;t++)y[i]-=c[i-t+s][t]*y[t];u[N-1]=y[N-1]/c[s][N-1];for(i=N-2;i>=0;i--){ for(t=i+1;t<=min(i+s,N-1);t++)y[i]-=c[i-t+s][t]*u[t];u[i]=y[i]/c[s][i];}b2=b1;zero=0;for(i=0;i<N;i++)zero+=b[i]*u[i];b1=1/zero;}while((fabs(b1-b2)/fabs(b1))>1e-12);return(b1);}void LUfenjie() //对矩阵做LU分解{ int k=0;int j=0;int i=0;int t=0;for(k=0;k<N;k++){ for(j=k;j<=min(k+s,N-1);j++)for(t=max(max(0,k-r),j-s);t<k;t++)c[k-j+s][j]-=c[k-t+s][t]*c[t-j+s][j];for(i=k+1;i<=min(k+r,N-1);i++){ for(t=max(max(0,i-r),k-s);t<k;t++)c[i-k+s][k]-=c[i-t+s][t]*c[t-k+s][k];c[i-k+s][k]/=c[s][k];}}}三、运行结果如下四、初始向量对计算结果的影响在本程序的编写中，为了方便起见，所以迭代向量的初值选为u=[1 1 1 ...1]。

一、问题分析及算法描述1. 问题的提出：（1）用幂法、反幂法求矩阵A =[a ij ]20×20的按摸最大和最小特征值，并求出相应的特征向量。

其中 a ij ={sin (0.5i +0.2j ) i ≠j 1.5cos (i +1.2j ) i =j要求：迭代精度达到10−12。

（2）用带双步位移的QR 法求上述的全部特征值，并求出每一个实特征值相应的特征向量。

2. 算法的描述：(1) 幂法幂法主要用于计算矩阵的按摸为最大的特征值和相应的特征向量。

其迭代格式为：{ 任取非零向量u 0=(h 1(0)，⋯，h n (0))T|h r (k−1)|=max 1≤j≤n |h r (k−1)| y ⃑ k−1=u ⃑ k−1|h r (k−1)| u ⃑ k =Ay ⃑ k−1=(h 1(k )，⋯，h n (k ))T βk =sgn (h r (k−1))h r (k ) (k =1，2，⋯) 终止迭代的控制选用≤ε。

幂法的使用条件为n ×n 实矩阵A 具有n 个线性无关的特征向量x 1，x 2，⋯，x n ，其相应的特征值λ1，λ2，⋯，λn 满足不等式|λ1|>|λ2|≥|λ3|≥⋯≥|λn |或λ1=λ2=⋯=λm|λ1|>|λm+1|≥|λm+2|≥⋯≥|λn |幂法收敛速度与比值|λ2λ1|或|λm+1λ1|有关，比值越小，收敛速度越快。

(2) 反幂法反幂法用于计算n ×n 实矩阵A 按摸最小的特征值，其迭代格式为：{任取非零向量u 0∈R nηk−1=√u ⃑ k−1T u ⃑ k−1 y ⃑ k−1=u ⃑ k−1ηk−1⁄ Au ⃑ k =y ⃑ k−1 βk =y ⃑ k−1u ⃑ k (k =1，2，⋯) 每迭代一次都要求解一次线性方程组Au ⃑ k =y ⃑ k−1。

当k 足够大时，λn ≈1βk ，y ⃑ k−1可近似的作为矩阵A 的属于λn 的特征向量。

数值计算解矩阵的按模最大最小特征值及对应的特征向量一.幂法1. 幂法简介：当矩阵A 满足一定条件时，在工程中可用幂法计算其主特征值(按模最大)及其特征向量。

矩阵A 需要满足的条件为： (1) 的特征值为A i n λλλλ,0||...||||21≥≥≥>(2) 存在n 个线性无关的特征向量，设为n x x x ,...,,211.1计算过程：i ni i i u xx αα,1)0()0(∑==，有对任意向量不全为0，则有1111112211211111111011)()(...u u a u a u λu λαu αA x A Ax x k n n k n k k ni ik i i ni i i k )(k (k))(k αλλλλλα++++=+=+++≈⎥⎦⎤⎢⎣⎡+++======∑∑ 可见，当||12λλ越小时，收敛越快；且当k 充分大时，有1)1111)11111λαλαλ=⇒⎪⎩⎪⎨⎧==+++(k )(k k(k k )(k x x u x u x ，对应的特征向量即是)(k x 1+。

2 算法实现.,, 3,,1 , ).5()5(,,,,||).4();max(,).3()(max(;0,1).2(,).1()()()(停机否则输出失败信息转置若转否则输出若计算最大迭代次数，误差限，初始向量输入矩阵βλβεβλβλε←+←<<-←←=←←k k N k y x Ay x x abs x y k N x A k k k3 matlab 程序代码function [t,y]=lpowerA,x0,eps,N) % t 为所求特征值，y是对应特征向量k=1;z=0; % z 相当于λy=x0./max(abs(x0)); % 规范化初始向量x=A*y; % 迭代格式b=max(x); % b 相当于βif abs(z-b)<eps % 判断第一次迭代后是否满足要求t=max(x);return;endwhile abs(z-b)>eps && k<Nk=k+1;z=b;y=x./max(abs(x));x=A*y;b=max(x);end[m,index]=max(abs(x)); % 这两步保证取出来的按模最大特征值t=x(index); % 是原值，而非其绝对值。