数值分析幂法和反幂法
- 格式:doc
- 大小:357.50 KB
- 文档页数:24
数值分析幂法和反幂法数值分析中的幂法和反幂法是求解矩阵最大特征值和最小特征值的常用方法。
这两种方法在许多数值计算问题中都有着广泛的应用,包括图像压缩、数据降维、谱聚类等。
幂法(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存在或数值稳定性。
幂法反幂法求解矩阵最大最小特征值及其对应的特征向量幂法和反幂法是求解矩阵最大最小特征值及其对应特征向量的常用方法。
在本文中,我们将详细介绍这两种方法的原理和具体实现。
一、幂法(Power Method)幂法是一种迭代算法,用于求解矩阵的最大特征值及其对应的特征向量。
其基本思想是通过多次迭代得到矩阵的一个特征值和特征向量的近似值,并使其逼近真实值。
幂法的原理如下:1.初始化一个非零向量b0作为初始特征向量;2.计算b0的归一化向量b0/,b0,得到新的向量b1;3.计算矩阵A和向量b1的乘积Ab1,得到新的向量b2;4.对b2进行归一化,得到新的向量b3;5.重复步骤3和步骤4,直到b的变化趋于稳定;6.计算矩阵A和向量b的乘积Ab,得到新的向量b;7.特征值的近似值λ=,Ab,/,b。
具体实现如下:1.初始化一个非零向量b0;2.迭代n次进行如下操作:a. 计算bn=A*bn-1;b. 将bn进行归一化,得到bn=bn/,bn;3. 计算特征值的近似值lambda=,A*bn,/,bn;4. 特征向量的近似值vbn=bn。
幂法的优点是计算简单、迭代次数少,但对于含有多个特征值接近的矩阵,可能会收敛到次大特征值。
二、反幂法(Inverse Power Method)反幂法是幂法的拓展,用于求解矩阵的最小特征值及其对应的特征向量。
其基本思想是通过多次迭代得到矩阵的一个特征值和特征向量的近似值,并使其逼近真实值。
反幂法的原理如下:1.初始化一个非零向量b0作为初始特征向量;2.计算b0的归一化向量b0/,b0,得到新的向量b1;3.计算矩阵A的逆矩阵Ai和向量b1的乘积Ai*b1,得到新的向量b2;4.对b2进行归一化,得到新的向量b3;5.重复步骤3和步骤4,直到b的变化趋于稳定;6.计算矩阵A的逆矩阵Ai和向量b的乘积Ai*b,得到新的向量b;7.特征值的近似值λ=,Ai*b,/,b。
具体实现如下:1.初始化一个非零向量b0;2.迭代n次进行如下操作:a. 计算bn=inv(A)*bn-1;b. 将bn进行归一化,得到bn=bn/,bn;3. 计算特征值的近似值lambda=,inv(A)*bn,/,bn;4. 特征向量的近似值vbn=bn。
幂法-反幂法求解矩阵最大最小特征值及其对应的特征向量数值计算解矩阵的按模最大最小特征值及对应的特征向量一.幂法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)<="">t=max(x);return;endwhile abs(z-b)>eps && k<n< p="">k=k+1;z=b;y=x./max(abs(x));x=A*y;b=max(x);end[m,index]=max(abs(x)); % 这两步保证取出来的按模最大特征值t=x(index); % 是原值,而非其绝对值。
数值分析幂法和反幂法数值分析中,幂法(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.实际应用:幂法和反幂法在实际中广泛应用于各个领域,例如物理、工程、计算机科学等。
比如在结构力学中,幂法可以用来求解结构的自振频率和相应的振型;在电力系统中,反幂法可以用来求解电力系统决定性特征值,例如功率稳定性的最小特征值。
矩阵的特征值与特征向量的计算摘要物理,力学,工程技术中的很多问题在数学上都归结于求矩阵特征值的问题,例如振动问题(桥梁的振动,机械的振动,电磁振动等)、物理学中某些临界值的确定问题以及理论物理中的一些问题。
矩阵特征值的计算在矩阵计算中是一个很重要的部分,本文使用幂法和反幂法分别求矩阵的按模最大,按模最小特征向量及对应的特征值。
幂法是一种计算矩阵主特征值的一种迭代法,它最大的优点是方法简单,对于稀疏矩阵比较合适,但有时收敛速度很慢。
其基本思想是任取一个非零的初始向量。
由所求矩阵构造一向量序列。
再通过所构造的向量序列求出特征值和特征向量。
反幂法用来计算矩阵按模最小特征向量及其特征值,及计算对应于一个给定近似特征值的特征向量。
本文中主要使用反幂法计算一个矩阵的按模最小特征向量及其对应的特征值。
计算矩阵按模最小特征向量的基本思想是将其转化为求逆矩阵的按模最大特征向量。
然后通过这个按模最大的特征向量反推出原矩阵的按模最小特征向量。
关键词:矩阵;特征值;特征向量;冥法;反冥法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相关定理。
《数值分析》计算实习题目第一题: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)作为初始向量进行迭代,可得出以上结果。
一、问题分析与算法描述1. 问题的提出:〔1〕用幂法、反幂法求矩阵的按摸最大和最小特征值,并求出相应的特征向量。
其中要求:迭代精度达到。
〔2〕用带双步位移的QR法求上述的全部特征值,并求出每一个实特征值相应的特征向量。
2. 算法的描述:(1) 幂法幂法主要用于计算矩阵的按摸为最大的特征值和相应的特征向量。
其迭代格式为:终止迭代的控制选用。
幂法的使用条件为实矩阵A具有n个线性无关的特征向量,其相应的特征值满足不等式或幂法收敛速度与比值或有关,比值越小,收敛速度越快。
(2) 反幂法反幂法用于计算实矩阵A按摸最小的特征值,其迭代格式为:每迭代一次都要求解一次线性方程组。
当k足够大时,,可近似的作为矩阵A的属于的特征向量。
比值越小,收敛的越快。
反幂法要求矩阵A非奇异。
(3) 带双步位移的QR分解法QR方法适用于计算一般实矩阵的全部特征值,尤其适用于计算中小型实矩阵的全部特征值。
本算例中采用带双步位移的QR方法,可加速收敛,其迭代格式为:二、计算结果与分析1. 计算结果:(1) 幂法:初始条件:最大迭代次数L=1000;向量计算结果:第1次迭代结果:最大特征值:0.00000e+000第2次迭代结果:最大特征值:2.48910e+000 相对误差:1.00000e+000 第3次迭代结果:最大特征值:1.67719e+000 相对误差:第4次迭代结果:最大特征值:-2.10960e+000 相对误差:1.79503e+000 第5次迭代结果:最大特征值:-6.13203e-001 相对误差:2.44030e+000 ……第794次迭代结果:最大特征值:-1.97638e+000 相对误差:最大特征值:-1.97638e+000 相对误差:********************最终迭代结果***************特征值:-1.97638e+000 相对误差:迭代次数:795(2) 反幂法:初始条件:最大迭代次数L=1000;向量运行结果:第1次迭代结果:最大特征值:1.07542e+000第2次迭代结果:最大特征值:-3.66550e+000 相对误差:1.29339e+000 第3次迭代结果:最大特征值:1.22709e+001 相对误差:1.29871e+000 第4次迭代结果:最大特征值:-1.03421e+000 相对误差:1.28650e+001 第5次迭代结果:最大特征值:相对误差:……第995次迭代结果:最大特征值:相对误差:第996次迭代结果:最大特征值:相对误差:最大特征值:相对误差:第998次迭代结果:最大特征值:相对误差:第999次迭代结果:最大特征值:相对误差:第1000次迭代结果:最大特征值:相对误差:******************************超过最大设定迭代次数,迭代失败!(3) 带双步位移的QR法:初始条件:最大迭代次数L=1000;向量运行结果:全部特征值:特征向量〔经谱X数归一化〕:实特征值对应特征向量:-0.062705 -0.022368 0.304372 0.064466 0.521833 -0.157024 0.136942 -0.218108 0.250264 -0.043064 -0.228688 -0.184632 -0.072871 0.124721 0.029070 0.102566 -0.136358 0.167727 0.085747 0.546165 实特征值对应特征向量:-0.018001 0.019652 0.273447 0.070528 0.274896 -0.144015 0.048385 0.376439 -0.583051 -0.054008 -0.168682 -0.113430 -0.034709 0.009204 0.472291 0.125664 -0.190617 0.113145 0.046278 0.059871 实特征值对应特征向量:0.106861 0.087709 -0.024967 -0.020897 0.064302 0.034047 0.535143 0.046383 0.028832 0.003479-0.097276 -0.383801 0.089445 -0.039560 -0.036928 -0.021330 0.014811 0.705836 -0.108904 0.082022 实特征值对应特征向量:-0.055201 0.003399 0.242191 0.102847 0.372470 -0.372826 0.113953 0.240659 -0.310401 -0.076590 -0.244632 -0.192549 -0.077259 0.263328 0.201662 0.154166 -0.407814 0.186782 0.094649 0.173302 实特征值对应特征向量:0.427828 -0.546801 0.007822 -0.382580 0.025199 0.012788 0.033241 0.005389 -0.004065 0.043524 -0.032112 -0.044233 0.135395 -0.006564 0.001214 0.020165 0.011678 0.050001 -0.585765 0.013115 实特征值对应特征向量:0.236032 -0.139250 -0.008143 0.638527 -0.009049 -0.002911 -0.001307 0.003054 0.006515 -0.030134 0.012712 0.011368 -0.018792 -0.001753 -0.005749 -0.014290 -0.005292 -0.014591 0.717590 0.001369 实特征值对应特征向量:-0.227404 -0.048154 0.022615 0.297305 0.070372 0.039927 0.078503 0.015822 -0.012182 0.605334 -0.083616 -0.106270 -0.573963 -0.019907 0.003839 0.051362 0.036567 0.115613 0.332707 0.036954 实特征值对应特征向量:-0.027768 -0.051081 -0.159642 -0.054573 -0.084441 0.118378 0.029553 0.211088 0.203867 0.0486272. 结果分析以上三种方法中,幂法计算共进展了795次迭代才达到收敛,计算量较大,收敛性不好;反幂法计算结果未能收敛,通过进一步分析发现,这是因为反幂法迭代程序未考虑按模最小特征值为复数的情况,造成迭代失败。
幂法,反幂法求解矩阵最大最小特征值及其对应的特征向量数值计算解矩阵的按模最大最小特征值及对应的特征向量一.幂法1. 幂法简介:当矩阵A满足一定条件时,在工程中可用幂法计算其主特征值(按模最大) 及其特征向量。
矩阵A需要满足的条件为:|,|,|,|,...,|,|,0,,为A的特征值(1) 12nix,x,...,x(2) 存在n个线性无关的特征向量,设为 12n1.1计算过程:n(0)(0)对任意向量x,有x,,u,,不全为0,则有 ,iii,1i(k,1)(k)k,1(0),,,xAx...Axnn11k,k,,,Aαuαλu,,iiiii11i,i,,,,, k,1k,1k,1n2,,,,,λu()au?()au11122nn,,,,11,,k,1,,,u111,2||可见,当越小时,收敛越快;且当k充分大时,有,111(k,)k,1(k,),,,xu,x,111(k,1),,,,x1,对应的特征向量即是。
)(k)(kkx,xu,,,111,2 算法实现,(1).输入矩阵A,初始向量x,误差限,最大迭代次数N(k)x(k),(2).k,1,,0;y,(k)max(abs(x),(3).计算x,Ay,,max(x);,,,,(4).若|,|,,输出,y,否则,转(5)(5).若 k,N, 置k,k,1,,,,,转3 ,否则输出失败信息,停机.3 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); % 是原值,而非其绝对值。