反幂法求矩阵特征值.docx

一． 问题描述
用幂法与反幂法求解矩阵特征值
求n阶方阵A的特征值和特征向量，是实际计算中常常碰到的问题，如：机械、结构或电磁振动中的固有值问题等。对于n阶矩阵A，若存在数和n维向量x满足
Ax=x （1）
则称为矩阵A的特征值，x为相应的特征向量。
由线性代数知识可知，特征值是代数方程
|I-A|=+a+…+a+a=0 （2）
的根。从表面上看，矩阵特征值与特征向量的求解问题似乎很简单，只需求解方程（2）的根，就能得到特征值，再解齐次方程组
（I-A）x=0 （3）
的解，就可得到相应的特征向量。
上述方法对于n很小时是可以的。但当n稍大时，计算工作量将以惊人的速度增大，并且由于计算带有误差，方程（2）未必是精确的特征方程，自然就不必说求解方程（2）与（3）的困难了。幂法与反幂法是一种计算矩阵主特征值及对应特征向量的迭代方法, 特别是用于大型稀疏矩阵。
这里用幂法与反幂法求解带状稀疏矩阵A[501][501]的特征值。

二． 算法设计
1． 幂法
（1）取初始向量u（例如取u=(1,1,…1)）,置精度要求，置k=1.
（2）计算
v=Au, m=max(v), u= v/ m
（3）若| m-m|<，则停止计算（m作为绝对值最大特征值，u作为相应的特征向量）否则置k=k+1，转（2）
2. 反幂法
（1）取初始向量u（例如取u=(1,1,…1)）,置精度要求，置k=1.
（2）对A作LU分解，即A=LU
（3）解线性方程组 Ly=u,Uv=y
（4）计算
m=max(v), u= v/ m
（5）若|m-m|<，则停止计算（1/m作为绝对值最小特征值，u作为相应的特征向量）;否则置k=k+1，转（3）.

三． 程序框图
1. 主程序

开始
设置二维数组Q[5][501]存放矩阵A的带内元素
平移求最大特征值和最小特征值
结束
fanmifa(A)反幂法求按模最小特征值
mifa(A)幂法求按模最大特征值
计算u[1]到u[39]

2. 子程序
(1). 幂法迭代程序框图

开始
计算幂法初始向量
fabs(b-c)/fabs(b)<=1e-12？

Return(b)并
printf("幂法成功！")
返回
>
b=sgn(h)*u[l];
找到Max(fabs(u[i]))
迭代求u[i]=u[i]+A[i-j+2][j]*y[j];
y[i]=u[i]/fabs(h);
c=b;

(2). 反幂法迭代程序框图
开始
计算反幂法初始向量

Return(1/b)并
printf("反幂法成功")
返回
>
LU分解前保存A[i][j]和y[i]的值
c=b;
LU分解
fabs(b-c)/fabs(b)<=1e-12？
解线性方程组求出u[i]
迭代求b=b+y[i]*u[i];

四． 结果显示

计算结果如下：
矩阵的按模最大特征值为:-1.070011361487e+001
矩阵的按模最小特征值为:-5.557910794230e-003
矩阵最大的特征值为:9.724634101479e+000
矩阵最小的特征值为:-1.070011361487e+001

与各 最接近的（用表示）的值如下：
v[ 1]=-1.018293403315e+001 u[ 1]=-1.018949492196e+001
v[ 2]=-9.585707425068e+000 u[ 2]=-9.678876229054e+000
v[ 3]=-9.172672423928e+000 u[ 3]=-9.168257536145e+000
v[ 4]=-8.652284007898e+000 u[ 4]=-8.657638843237e+000
v[ 5]=-8.093483808675e+000 u[ 5]=-8.147020150328e+000
v[ 6]=-7.659405407692e+000 u[ 6]=-7.636401457419e+000
v[ 7]=-7.119684648691e+000 u[ 7]=-7.125782764510e+000
v[ 8]=-6.611764339397e+000 u[ 8]=-6.615164071601e+000
v[ 9]=-6.066103226595e+000 u[ 9]=-6.104545378693e+000
v[10]=-5.585101052628e+000 u[10]=-5.593926685784e+000
v[11]=-5.114083529812e+000 u[11]=-5.083307992875e+000
v[12]=-4.578872176865e+000 u[12]=-4.572689299966e+000
v[13]=-4.096470926260e+000 u[13]=-4.062070607058e+000
v[14]=-3.554211215751e+000 u[14]=-3.551451914149e+000
v[15]=-3.041090018133e+000 u[15]=-3.040833221240e+000
v[16]=-2.533970311130e+000 u[16]=-2.530214528331e+000
v[17]=-2.003230769563e+000 u[17]=-2.019595835422e+000
v[18]=-1.503557611227e+000 u[18]=-1.508977142514e+000
v[19]=-9.935586060075e-001 u[19]=-9.983584496049e-001
v[20]=-4.870426738850e-001 u[20]=-4.877397566962e-001
v[21]=2.231736249575e-002 u[21]=2.287893621262e-002
v[22]=5.324174742069e-001 u[22]=5.334976291214e-001
v[23]=1.052898962693e+000 u[23]=1.044116322030e+000
v[24]=1.589445881881e+000 u[24]=1.554735014939e+000
v[25]=2.060330460274e+000 u[25]=2.065353707848e+000
v[26]=2.558075597073e+000 u[26]=2.575972400756e+000
v[27]=3.080240509307e+000 u[27]=3.086591093665e+000
v[28]=3.613620867692e+000 u[28]=3.597209786574e+000
v[29]=4.091378510451e+000 u[29]=4.107828479483e+000
v[30]=4.603035378279e+000 u[30]=4.618447172392e+000
v[31]=5.132924283898e+000 u[31]=5.129065865300e+000
v[32]=5.594906348083e+000 u[32]=5.639684558209e+000
v[33]=6.080933857027e+000 u[33]=6.150303251118e+000
v[34]=6.680354092112e+000 u[34]=6.660921944027e+000
v[35]=7.293877448127e+000 u[35]=7.171540636935e+000
v[36]=7.717111714236e+000 u[36]=7.682159329844e+000
v[37]=8.225220014050e+000 u[37]=8.192778022753e+000
v[38]=8.648666065193e+000 u[38]=8.703396715662e+000
v[39]=9.254200344575e+000 u[39]=9.214015408571e+000

五． 程序

#include<stdio.h>
#include<math.h>
#define N 501
void main()
{
double Q[5][501];
double mifa(double A[5][501]);
double fanmifa(double A[5][501]);
double lm,lmax,lmin,ls,delta,u[39],v[39];
int i,j,k;
double A[5][501];

A[0][0]=A[0][1]=A[1][0]=A[3][500]=A[4][499]=A[4][500]=0.0;//输入*501矩阵
for(i=2;i<N;i++)
A[0][i]=-0.064;
for(i=1;i<N;i++)
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<500;i++)
A[3][i]=0.16;
for(i=0;i<499;i++)
A[4][i]=-0.064;

for(i=0;i<5;i++)//保存A
for(j=0;j<501;j++)
Q[i][j]=A[i][j];

lm=mifa(A);//按模最大特征值,函数mifa()不会改变矩阵A的值，不需还原

for(i=0;i<N;i++) //平移A
{
A[2][i]=A[2][i]-lm;
}
lmax=mifa(A);//平移后A的按模最大特征值

lmax=lmax+lm;//最大特征值或最小特征值
if(lmax<lm)
{
lmin=lmax;
lmax=lm;
}
else
lmin=lm;
for(i=0;i<N;i++)//还原A
for(j=0;j<5;j++)
A[j][i]=Q[j][i];

ls=fanmifa(A);//按模最小特征值
for(i=0;i<N;i++)//还原A
for(j=0;j<5;j++)
A[j][i]=Q[j][i];

for(k=0;k<39;k++)//计算u1-u39
u[k]=lmin+(k+1)*((lmax-lmin)/40);
for(k=0;k<39;k++)
{
for(j=0;j<N;j++)
A[2][j]=A[2][j]-u[k];
v[k]=fanmifa(A)+u[k];
for(i=0;i<N;i++)//还原A
for(j=0;j<5;j++)
A[j][i]=Q[j][i];
}

printf("矩阵的按模最大特征值为:%.12e",lm);
printf("\n");
printf("矩阵的按模最小特征值为:%.12e",ls);
printf("\n");
printf("矩阵最大的特征值为:%.12e",lmax);
printf("\n");
printf("矩阵最小的特征值为:%.12e",lmin);
printf("\n");
for(k=0;k<39;k++)
{
printf("v[%2d]=%.12e ",k+1,v[k]);
printf("u[%2d]=%.12e",k+1,u[k]);
printf("\n");
}
}

double sgn(double a)//符号函数
{
if(a>0)
return 1;
else if(a=0)
return 0;
else return -1;
}

int max2(int a,int b)
{
return a>b?a:b;
}

int max3(int a,int b,int c)
{
return max2(a,b)>c?max2(a,b):c;
}

int min(int a,int b)
{
return a<b?a:b;
}

void LU(double A[5][501],double u[501],double B[501])//LU分解法
{
double X[501];
int i,j,k,t,l;
double m=0,n=0;
for(k=1;k<=N;k++)//求L，U
{
for(j=k;j<=min(N,k+2);j++)//U
{
m=0;
for(t=max3(1,k-2,j-2);t<=k-1;t++)
{
m+=A[k-t+2][t-1]*A[t-j+2][j-1];
}
A[k-j+2][j-1]=A[k-j+2][j-1]-m;
}
for(i=k+1;i<=min(N,k+2);i++)//L
if(k<N)
{
n=0;
for(l=max3(1,i-2,k-2);l<=k-1;l++)
{
n+=A[i-l+2][l-1]*A[l-k+2][k-1];
}
A[i-k+2][k-1]=(A[i-k+2][k-1]-n)/A[2][k-1];
}
}
for(i=2;i<=N;i++)//回代过程
{
m=0;
for(t=max2(1,i-2);t<=i-1;t++)
m+=A[i-t+2][t-1]*B[t-1];
B[i-1]=B[i-1]-m;
}
X[N-1]=B[N-1]/A[2][N-1];//回代过程
for(i=N-1;i>=1;i--)
{
n=0;
for(t=i+1;t<=min(N,i+2);t++)
n+=A[i-t+2][t-1]*X[t-1];
X[i-1]=(B[i-1]-n)/A[2][i-1];
}
for(i=1;i<=N;i++)//输出方程结果
{
u[i-1]=X[i-1];
}
}

double mifa(double A[5][501])//幂法
{
int i,j,l=0;
double u[501],t[501];
double y[501];
double h,b,c;
c=0;
for(i=0;i<N;i++)//幂法初始向量
u[i]=1;
while(1)
{
for(i=0;i<N;i++)
t[i]=0;
h=u[0];
for(i=0;i<N;i++)//无穷范数
{
if(fabs(h)<fabs(u[i]))
{
h=u[i];
l=i;
}
}
for(i=0;i<N;i++)
y[i]=u[i]/fabs(h);
for(i=2;i<499;i++)
{
for(j=i-2;j<=i+2;j++)
{
t[i]=t[i]+A[i-j+2][j]*y[j];
}
u[i]=t[i];
}
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[499]=A[4][497]*y[497]+A[3][498]*y[498]+A[2][499]*y[499]+A[1][N-1]*y[N-1];
u[N-1]=A[4][498]*y[498]+A[3][499]*y[499]+A[2][N-1]*y[N-1];
b=sgn(h)*u[l];
if((fabs(b-c)/fabs(b))<=1e-12)
{
//printf("幂法成功！");
//printf("\n");
break;
}
c=b;
}
return b;
}

double fanmifa(double A[5][501])//反幂法
{
double u[501],y[501];
double P[5][501],Y[501];//LU分解前用于保存A和y的值
double m=0,n=0,b=0,c=0;
int i,j;
for(i=0;i<N;i++)//反幂法初始向量
u[0]=1;
while(1)
{
b=0;
n=0;
for(i=0;i<N;i++)
n=n+u[i]*u[i];
n=sqrt(n);
for(i=0;i<N;i++)
y[i]=u[i]/n;
for(i=0;i<N;i++)//保存A和y
{
Y[i]=y[i];
for(j=0;j<5;j++)
{
P[j][i]=A[j][i];
}
}
LU(A,u,y);//LU分解法,会改变A，y，u的值（目的只需求出u）
for(i=0;i<N;i++)//还原A和y
{
y[i]=Y[i];
for(j=0;j<5;j++)
{
A[j][i]=P[j][i];
}
}
for(i=0;i<N;i++)
b=b+y[i]*u[i];
if((fabs(b-c)/fabs(b))<=1e-12)
{
//printf("反幂法成功！");
//printf("\n");
break;
}
c=b;
}
return 1/b;
}