带原点平移的反幂法解特征值
数值分析幂法和反幂法
k=k+1;
end
在matlab输入面板，输入
A=rand（4）；%产生一个4维随机矩阵
B=A+A’；
u=[1 1 1 1]’;%设立初始向量
[m,u,index,k]=pow(B，u,ep,it_max)%最多可省略2个参数
程序结束。
在M文件中可以通过改变m0的值改变原点位移，从而达到原点位移加速。
m =
2.6814
u =
0.8576
0.6934
0.5623
1.0000
index =
0
k =
1001
修改M0=0
m =
2.6820
u =
0.8577
0.6937
0.5624
1.0000
index =
1
k =
7
总结以上,幂法如下：
U
m0
m
u
index
k
[1 1 1 1]
0.0001
2.6813
[0.8576 0.6934 0.5623 1.0000]
2、对于幂法的定理
按式（1）计算出m 和u 满足
m = , u =
（二）反幂法算法的理论依据及推导
反幂法是用来计算绝对值最小的特征值忽然相应的特征向量的方法。是对幂法的修改，可以给出更快的收敛性。
1、反幂法的迭代格式与收敛性质
设A是非奇异矩阵，则零不是特征值，并设特征值为
| |≥| |≥…≥| |>| |
反幂法程序设计代码：
在matlab中建立一个M文件并保存。
%pow_inv.m
function[m,u,index,k]=pow_inv(A,u,ep,it_max)
七、矩阵特征值的乘幂方法和反乘幂方法
1、用幂法计算矩阵A的主特征值和对应的特征向量。
function [k,lambda,Vk,Wc]=mifa(A,V0,jd,max1)lambda=0;k=1;Wc=1;jd=jd*0.1;state=1;V=V0;while((k<=max1)&(state==1))Vk=A*V;[m j]=max(abs(Vk));mk=m*sign(Vk(j));tzw=abs(lambda-mk);Vk=(1/mk)*Vk;txw=norm(V-Vk);Wc=max(txw,tzw);V=Vk;lambda=mk;state=0;if(Wc>jd)state=1;endk=k+1;endif(Wc<=jd)disp('迭代次数k，主特征值的近似值lambda，主特征向量的近似向量Vk，相邻两次迭代的误差Wc如下：')elsedisp('迭代次数k已经达到最大迭代次数max1，主特征值的迭代值lambda，主特征向量的迭代向量Vk，相邻两次迭代的误差Wc如下：')endVk=V;k=k-1;Wc;>> A=[1 -1;2 4];>> V0=[1,1]';>> [k,lambda,Vk,Wc]=mifa(A,V0,0.00001,100)迭代次数k，主特征值的近似值lambda，主特征向量的近似向量Vk，相邻两次迭代的误差Wc如下：k =33lambda =3.0000Vk =-0.50001.0000Wc =8.6919e-007>> [V,D]=eig(A)V =-0.7071 0.44720.7071 -0.8944D =2 00 3>> Dzd=max(diag(D))Dzd =3>> wuD=abs(Dzd-lambda)wuD =1.7384e-006>> wuV=V(:,2)./VkwuV =-0.8944-0.8944>> A=[1 2 3;2 1 3;3 3 6];>> V0=[1 1 1]';>> [k,lambda,Vk,Wc]=mifa(A,V0,0.00001,100)迭代次数k，主特征值的近似值lambda，主特征向量的近似向量Vk，相邻两次迭代的误差Wc如下：k =3lambda =9Vk =0.50000.50001.0000Wc =>> [V,D]=eig(A)V =0.7071 0.5774 0.4082-0.7071 0.5774 0.40820 -0.5774 0.8165D =-1.0000 0 00 -0.0000 00 0 9.0000>> Dzd=max(diag(D))Dzd =9>> wuD=abs(Dzd-lambda)wuD =>> wuV=V(:,2)./VkwuV =1.15471.1547-0.5774>> A=[1 2 2;1 -1 1;4 -12 1];>> V0=[1 1 1]';>> [k,lambda,Vk,Wc]=mifa(A,V0,0.00001,100)迭代次数k已经达到最大迭代次数max1，主特征值的迭代值lambda，主特征向量的迭代向量Vk，相邻两次迭代的误差Wc如下：k =100lambda =-0.0909Vk =1.00001.00001.0000Wc =1.9582>> [V,D]=eig(A)V =0.9045 -0.7255 -0.72550.3015 -0.2176 - 0.0725i -0.2176 + 0.0725i-0.3015 0.5804 - 0.2902i 0.5804 + 0.2902iD =1.0000 0 00 -0.0000 + 1.0000i 00 0 -0.0000 - 1.0000i>> Dzd=max(diag(D))Dzd =1.0000>> wuD=abs(Dzd-lambda)wuD =1.0909>> wuV=V(:,2)./VkwuV =-0.7255-0.2176 - 0.0725i0.5804 - 0.2902i（4）>> A=[-4 14 0;-5 13 0;-1 0 2];>> V0=[1 1 1]';>> [k,lambda,Vk,Wc]=mifa(A,V0,0.00001,100)迭代次数k，主特征值的近似值lambda，主特征向量的近似向量Vk，相邻两次迭代的误差Wc如下：k =22lambda =6.0000Vk =1.00000.7143-0.2500Wc =8.1744e-007>> [V,D]=eig(A)V =0 0.7974 0.66670 0.5696 0.33331.0000 -0.1994 -0.6667D =2.0000 0 00 6.0000 00 0 3.0000>> Dzd=max(diag(D))Dzd =6.0000>> wuD=abs(Dzd-lambda)wuD =8.1744e-007>> wuV=V(:,2)./VkwuV =0.79740.79740.79742、用原点位移反幂法计算矩阵A的特征值和对应的特征向量。
幂法和反幂法
此例中比值为 2 2 . 1 3
例2：用幂法计算下面矩阵的主特征值及对应的特征向量。
解： 取初始向量 01
2 4 A 3 9
4 16
v u 1 1 1 ,按(3.7)迭代5次得到数据T如下 表: 00
1 11
11
ukT
6 15 36
k
vkT
(规范化向量)
5 0.1859 0.4460 1 8.156 19.57 43.88
v (i) k 1 v (i) k
1?
即两相邻迭代向量的对应非零分量的比值一定收敛到主特征值？
不一定. 先讨论以下情况：
情形1: 设n n阶实矩阵A的特征值i (i 1, 2, , n) 满足 1 2 n 且与i (i 1, 2, , n)相应的特征
向量x1 , x2 , , xn 线性无关。
v (1) 2
v (1) 1
0.41 ,
v (2) 2
v (2) 1
0.41666,
v (1) 3
0.41260,
v (2) 3
0.41249,
v (1) 2
v (2) 2
v (1) 4
v (1) 3
0.41263,
v (2) 4
v (2) 3
0.41263,
问题:是否任何矩阵的幂法，当k比较大时，一定有
故按模特征值为:
1 43.88 对应的特征向量为:
u1 0.1859 0.4460 1.0000T
例3 用幂法求矩阵 的主特征值和主特征向量.
1 1 0.5 A 1 1 0.25
0.5 0.25 2
解 : 取初始向量u0 (1,1,1)T , 按(3.2)的计算结果如表9 1。
41第一节-幂法和反幂法
=a11v1 a22v2 annvn x2 Ax1 a11 Av1 a22 Av2 ann Avn
=a112v1 a222v2 ann2vn
即对 x0 a1v1 a2v2 anvn 用公式 xk Axk1, k 1, 2,
幂法的迭代公式为
xk Axk1 k 1, 2,
当k充分大时，有
xk
1ka1v1
1
xk1 i xk i
收敛速度取决于比值 2 ,比值越小,收敛越快. 1
3. 误差分析
幂法的迭代公式为 xk Axk1 k 1, 2,
当k充分大时, 有
xk 1ka1v1
称1为A的按摸最大特征值(也称主特征值).
任取非零向量
x0
( x1(0) ,
x(0) 2
,
,
x(0) n
)T,
则
x0 a1v1 a2v2 anvn
设 a10, 由A构造向量序列{xk}
xk Axk1, k 1, 2,
其中 x1 Ax0 a1 Av1 a2 Av2 an Avn
a1n a2n 0
an1
an2
ann
的根；求A的属于特征值的特征向量等价于求
非零解.
( E A)x 0
设为A∈Rn×n的特征值, x 称为A的与特 征值 相对应的一个特征向量,即Ax= x, (x≠0)
则有
(1) cx (c≠0为常数)也是 A的与特征值 相对 应的一个特征向量，即A(cx)=(cx)；
4. 实用计算公式
yk Axk1
mk max yk
矩阵特征值的计算
物理、力学和工程技术中的许多问题在数学上都归结为求矩 阵的特征值和特征向量问题。
� 计算方阵 A 的特征值，就是求特征多项式方程：
| A − λI |= 0 即 λn + p1λn−1 + ⋅ ⋅ ⋅ + pn−1λ + pn = 0
的根。求出特征值 λ 后，再求相应的齐次线性方程组：
(13)
为了防止溢出，计算公式为
⎧ Ay k = xk −1
⎪ ⎨
m
k
=
max(
yk )
( k = 1, 2, ⋅ ⋅⋅)
⎪ ⎩
x
k
=
yk
/ mk
(14)
相应地取
⎧ ⎪
λ
n
⎨
≈
1 mk
⎪⎩ v n ≈ y k ( 或 x k )
(15)
9
(13)式中方程组有相同的系数矩阵 A ，为了节省工作量，可先对
11
11
≤ ≤ ⋅⋅⋅ ≤
<
λ1 λ2
λn −1
λn
对应的特征向量仍然为 v1, v2 ,⋅⋅⋅, vn 。因此，计算矩阵 A 的按模
最小的特征值，就是计算 A−1 的按模最大的特征值。
� 反幂法的基本思想：把幂法用到 A−1 上。
任取一个非零的初始向量 x0 ，由矩阵 A−1 构造向量序列：
xk = A−1xk−1 , k = 1, 2, ⋅⋅⋅
如果 p 是矩阵 A 的特征值 λi 的一个近似值，且
| λi − p |<| λ j − p | , i ≠ j
1 则 λ i − p 是矩阵 ( A − pI )−1 的按模最大的特征值。因此，当给
北航数值分析计算实习题目一 幂法 反幂法 求矩阵特征值
《数值分析》计算实习题目第一题：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;double lambda[40];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 u[502],y[502];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）作为初始向量进行迭代，可得出以上结果。
数值方法课程设计幂法反幂法计算矩阵特征值和特征向量-附Matlab程序
矩阵的特征值与特征向量的计算摘要物理，力学，工程技术中的很多问题在数学上都归结于求矩阵特征值的问题，例如振动问题（桥梁的振动，机械的振动，电磁振动等）、物理学中某些临界值的确定问题以及理论物理中的一些问题。
矩阵特征值的计算在矩阵计算中是一个很重要的部分，本文使用幂法和反幂法分别求矩阵的按模最大，按模最小特征向量及对应的特征值。
幂法是一种计算矩阵主特征值的一种迭代法，它最大的优点是方法简单，对于稀疏矩阵比较合适，但有时收敛速度很慢。
其基本思想是任取一个非零的初始向量。
由所求矩阵构造一向量序列。
再通过所构造的向量序列求出特征值和特征向量。
反幂法用来计算矩阵按模最小特征向量及其特征值，及计算对应于一个给定近似特征值的特征向量。
本文中主要使用反幂法计算一个矩阵的按模最小特征向量及其对应的特征值。
计算矩阵按模最小特征向量的基本思想是将其转化为求逆矩阵的按模最大特征向量。
然后通过这个按模最大的特征向量反推出原矩阵的按模最小特征向量。
关键词：矩阵；特征值；特征向量；冥法；反冥法THE CALCULATIONS OF EIGENVALUE AND EIGENVECTOR OF MATRIXABSTRACTPhysics, mechanics, engineering technology in a lot of problems in mathematics are attributed to matrix eigenvalue problem, such as vibration (vibration of the bridge, mechanical vibration, electromagnetic vibration, etc.) in physics, some critical values determine problems and theoretical physics in some of the problems. Matrix eigenvalue calculation is a very important part in matrix computation. In this paper, we use the power method and inverse power method to calculate the maximum of the matrix, according to the minimum characteristic vector and the corresponding characteristic value.Power method is an iterative method to calculate the eigenvalues of a matrix. It has the advantage that the method is simple and suitable for sparse matrices, but sometimes the convergence rate is very slow. The basic idea is to take a non - zero initial vector. Construct a vector sequence from the matrix of the matrix. Then the eigenvalues and eigenvectors are obtained by using the constructed vector sequence.The inverse power method is used to calculate the minimum feature vectors and their eigenvalues of the matrix, and to calculate the eigenvalues of the matrix. In this paper, we use the inverse power method to calculate the minimum eigenvalue of a matrix and its corresponding eigenvalues. The basic idea of calculating the minimum characteristic vector of a matrix is to transform it to the maximumc haracteristic vector of the modulus of the inverse matrix. Then, according to the model, the minimum feature vector of the original matrix is introduced.Key words: Matrix ；Eigenvalue ；Eigenvector ；Iteration methods;目录1引言 (1)2相关定理。
幂法和反幂法求矩阵特征值课程知识讲解
v (k ) =Au (k 1) ,m =max(v (k ) ), u (k ) = v (k ) / m
k
k
（3）若|
m= k
m k 1 |<
，则停止计算（m k 作为绝对值最大特征值 1 ，u (k) 作为
相应的特征向量）否则置 k=k+1，转（2）
2、反幂法算法
（1）取初始向量 u (0) （例如取 u (0) =(1,1,…1) T ）,置精度要求 ，置 k=1.
要
2.选择合适问题求解的数值计算方法；
求
3.设计程序并进行计算；
4.对结果进行解释说明；
对于幂法和反幂法求解矩阵特征值和特征向量的问题将从问题分析，算 法设计和流程图，理论依据，程序及结果进行阐述该问题。
一．问题的分析：
求 n 阶方阵 A 的特征值和特征向量，是实际计算中常常碰到的问题，如：
采
机械、结构或电磁振动中的固有值问题等。对于 n 阶矩阵 A，若存在数 和
按式（1）计算出 m 和 u (k ) 满足 k
lim
k
m
k
=
1
,
lim u (k ) = x1
k
max( x1 )
（二）反幂法算法的理论依据及推导
反幂法是用来计算绝对值最小的特征值忽然相应的特征向量的方法。是对 幂法的修改，可以给出更快的收敛性。 1、反幂法的迭代格式与收敛性质
设 A 是非奇异矩阵，则零不是特征值，并设特征值为 | 1 |≥| 2 |≥…≥| n1|>| n |
则按 A 1 的特征值绝对值的大小排序，有
| 1 |>| 1 |≥…≥| 1 |
n
n 1
1