幂法与反幂法

幂法与反幂法
1 功能 幂法是一种计算矩阵主特征值(矩阵按模最大的特征值)及对应特征向量的迭代方法, 特别是用于大型稀疏矩阵。

反幂法用来计算矩阵按模最小的特征值及其特征向量，也可用来计算对应与一个给定近似特征值的特征向量。

2算法描述
2.1 幂法
（1）取初始向量u
)0(（例如取u )0(=(1,1,…1)T ）,置精度要求ε，置k=1. （2）计算
v )(k =Au )1(-k ,m k =max(v )(k ), u )(k = v )(k / m k
（3）若| m k = m 1-k |<ε，则停止计算（m k 作为绝对值最大特征值1λ，u
)(k 作为相应
的特征向量）否则置k=k+1，转（2） 2.2 反幂法
（1）取初始向量u )0(（例如取u )0(=(1,1,…1)T
）,置精度要求ε，置k=1. （2）对A 作LU 分解，即A=LU
（3）解线性方程组 Ly
)(k =u )1(-k ,Uv )(k =y )(k （4）计算
m k =max(v )(k ), u )(k = v )(k / m k
（5）若|m k =m 1-k |<ε，则停止计算（1/m k 作为绝对值最小特征值n λ，u
)(k 作为相应
的特征向量）;否则置k=k+1，转（3）.
3 Matlab 程序的实现
3.1 幂法
function [m,u]=pow(A,ep,N)
%A 为矩阵；ep 为精度要求；N 为最大迭代次数；m 为绝对值最大的特征值；u 为对应最大特征值的特征向量。

N=100;
ep=1e-6;
n=length(A);u=ones(n,1);
index=0;k=0;m1=0;
while k<=N
v=A*u;[vmax,i]=max(abs(v));
m=v(i);u=v/m;
if abs(m-m1)<ep
index=1;break;
end
m1=m;k=k+1;
end
输入：
A=[7 3 -2;3 4 -1;-2 -1 3];
[m,u]=pow(A,1e-6) Enter
结果：
m = 9.6056
u =1.0000
0.6056
-0.3944
4.2 反幂法
function[m ,u]=pow_inv(A,ep,N)
%A为矩阵；ep为精度要求；N为最大迭代次数；m为绝对值最大的特征值；u为对应最大特征值的特征向量。

N=100;
ep=1e-6;
n=length(A);u=ones(n,1);
index=0;k=0;m1=0;
invA=inv(A);
while k<=N
v=invA*u;[vmax,i]=max(abs(v));
m =v(i);u=v/ m ;
if abs(m-m1)<ep
index=1;break;
end
m1=m ;k=k+1;
end
m =1/ m ;
输入：
A=[6 2 1;2 3 1;1 1 3];
[m,u]=pow_inv(A,1e-6) Enter
结果：
m = 1.7733
u = -0.3473
1.0000
-0.5321
4 Matlab程序调用示例
随机产生一对称矩阵，对不同的原点位移和初值(至少取3个)分别使用幂法求计算矩阵的主特征值及主特征向量，用反幂法求计算矩阵的按模最小特征值及特征向量。

要求1）比较不同的原点位移和初值说明收敛性
2）给出迭代结果，生成DOC 文件。

3）程序清单，生成M 文件。

wo
解答：
>> A=rand(5) %随机产生5*5矩阵 求随机矩阵
A =
0.7094 0.1626 0.5853 0.6991 0.1493
0.7547 0.1190 0.2238 0.8909 0.2575
0.2760 0.4984 0.7513 0.9593 0.8407
0.6797 0.9597 0.2551 0.5472 0.2543
0.6551 0.3404 0.5060 0.1386 0.8143
>> B=A+A' %A 矩阵和A 的转置相加，得到随机对称矩阵B
B =
1.4187 0.9173 0.8613 1.3788 0.8044
0.9173 0.2380 0.7222 1.8506 0.5979
0.8613 0.7222 1.5025 1.2144 1.3467
1.3788 1.8506 1.2144 1.0944 0.3929
0.8044 0.5979 1.3467 0.3929 1.6286
B=⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡6286.13929.03467.15979.08044.03929.00944.12144.18506.13788.13467.12144.15025.17222.08613.05979.08506.17222.02380.09173.08044.03788.18613.09173.04187.1
5 C 语言程序实现
#include "stdio.h"
#include "math.h"
#define M 3
void main()
{
float fan(),max(),e1,e2,r1,r2;
void au(),ex(),print_x(),std();
static float a[M][M]={{1.0,1.0,0.5},{1.0,1.0,0.25},{0.5,0.25,2.0}}; static float u0[M],u1[M],maxn0,maxn1;
int i;
printf("*********************************\n");
printf("****** 幂法 *********\n");
printf("******求特征值与特征向量*********\n");
printf("*********************************\n\n");
printf("input precision e1,e2:");
scanf("%f,%f",&e1,&e2);
printf("\ninput u(%d):",M);
for (i=0;i<M;++i)
{
scanf("%f",&u0[i]);
}
std(u0);
maxn0=max(u0);
i=0;
printf("\n- - - - - - - - - - - - - - - - - -\n");
printf(" ............NMD\n");
do
{
au(a,u0,u1);
maxn1=max(u1);
std(u1);
r1=fan(u0,u1);
r2=(float)fabs(maxn0-maxn1);
maxn0=maxn1;
if (r1>e1 || r2>e2)
{
printf("%4d",i++);
print_x(u0);
printf("\n");
ex(u0,u1);
}
else
break;
} while (1);
}
void au(a,u0,u1)
float a[][M],u0[],u1[];
{
int i,j;
for (i=0;i<M;i++)
{
u1[i]=0;
for (j=0;j<M;j++)
{
u1[i]+=a[i][j]*u0[j];
}
}
}
void std(u)
float u[];
{
int i;
float t,max();
t=max(u);
for (i=0;i<M;i++)
{
u[i]=u[i]/t;
}
}
float fan(u0,u1)
float u0[],u1[];
{
float max();
int i;
float uu[M];
for (i=0;i<M;i++)
{
uu[i]=u0[i]-u1[i];
}
return max(uu);
}
float max(u)
float u[];
{
int i;
float m;
m=u[0];
for (i=0;i<M;i++)
{
if (u[i]>m)
{
m=u[i];
}
}
return m;
}
void ex(u0,u1)
float u0[],u1[];
{
int i;
for (i=0;i<M;i++) {
u0[i]=u1[i]; }
}
void print_x(u)
float u[];
{
int i;
for (i=0;i<M;i++) {
printf("%12.6f",u[i]); }
}
6 C 程序示例调用
用幂法求 A=
⎥⎥⎥⎦⎤⎢⎢⎢⎣
⎡0.225.05.025.00.10.15.00.10.1 的特征值和特征向量。

运行结果如下：。