实验6反幂法求矩阵按模最小特征值
数值方法课程设计幂法反幂法计算矩阵特征值和特征向量-附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 maximum characteristic 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 相关定理。
实验6反幂法求矩阵按模最小特征值
西华数学与计算机学院上机实践报告课程名称:计算方法A年级: 上机实践成绩: 指导教师:严常龙姓名: 上机实践名称:反幂法求矩阵按模最小特征值 学号:上机实践日期: 上机实践编号:1上机实践时间: 一、目的1.通过本实验加深对反幂法的构造过程的理解;2.能对反幂法提出正确的算法描述编程实现,得到计算结果。
二、内容与设计思想自选方阵,用反幂法求解其按模最小特征值。
可使用实例:⎪⎪⎪⎭⎫ ⎝⎛---=90688465441356133A三、使用环境操作系统:Windows XP软件环境:Microsoft Visual C++四、核心代码及调试过程#include<stdio.h>#include<math.h>#define MAX_N 20 //矩阵最大维数#define MAXREPT 100#define epsilon 0.00001 //求解精度int main(){int n;int i,j,k;double xmax,oxmax;static double a[MAX_N][MAX_N];static double l[MAX_N][MAX_N],u[MAX_N][MAX_N];static double x[MAX_N],nx[MAX_N];printf("\n 请输入矩阵阶数n:"); //输入矩阵维数scanf("%d",&n);if(n>MAX_N){printf("the input n is larger than MAX_N,please redefine the MAX_N.\n");return 1;}if(n<=0){printf("please input a number between 1 and %d.\n",MAX_N);return 1;}//输入A矩阵printf("请输入矩阵的值a[i][j] i,j=0...%d;\n",n-1);for(i=0;i<n;i++)for(j=0;j<n;j++)scanf("%lf",&a[i][j]);for(i=0;i<n;i++)x[i]=1;oxmax=0;for(i=0;i<MAXREPT;i++){for(j=0;j<n;j++) //幂乘{nx[j]=0;for(k=0;k<n;k++)nx[j]+=a[j][k]*x[k];}xmax=0.0;for(j=0;j<n;j++) //规范化if(fabs(nx[j])>xmax)xmax=fabs(nx[j]);for(j=0;j<n;j++)nx[j]/=xmax;for(j=0;j<n;j++)x[j]=nx[j];if(fabs(xmax-oxmax)<epsilon){printf("solve...max lamda=%lf\n",xmax); //输出printf("the vector is:\n");for(i=0;i<n;i++)printf("%f\n",nx[i]);break;//return 0;}oxmax=xmax;}//printf("after %d repeat ,max no result ...\n",MAXREPT); for(i=0;i<n;i++)u[i][i]=1; //U矩阵对角元为for(k=0;k<n;k++){for(i=k;i<n;i++) //计算L矩阵{l[i][k]=a[i][k];for(j=0;j<=k-1;j++)l[i][k]-=(l[i][j]*u[j][k]);}for(j=k+1;j<n;j++) //计算U矩阵{u[k][j]=a[k][j];for(i=0;i<=k-1;i++)u[k][j]-=(l[k][i]*u[i][j]);u[k][j]/=l[k][k];}}for(i=0;i<n;i++)x[i]=1;for(i=0;i<MAXREPT;i++){for(j=0;j<n;j++) //反幂乘{nx[j]=x[j];for(k=0;k<=j-1;k++)nx[j]-=l[j][k]*nx[k];nx[j]/=l[j][j];}for(j=n-1;j>=0;j--){x[j]=nx[j];for(k=j+1;k<n;k++)x[j]-=u[j][k]*x[k];}xmax=0.0;for(j=0;j<n;j++)if(fabs(x[j])>xmax)xmax=fabs(x[j]);for(j=0;j<n;j++)x[j]/=xmax;if(fabs(xmax-oxmax)<epsilon){printf("solve... min lamda=%lf\n",1/xmax); //输出printf("the vector is:\n");for(i=0;i<n;i++)printf("%f\n",x[i]);break;//return 0;}oxmax=xmax;}return 1;}五、总结本次试验利用C语言实现了用反幂法求解其按模最小特征值,通过本实验加深对反幂法的构造过程的理解。
逆幂迭代法求特征值
逆幂迭代法求特征值
逆幂迭代法是一种求解特征值的方法,主要应用于对称矩阵或非奇异矩阵。
它通过不断迭代矩阵的逆幂次来逼近矩阵的最小特征值和对应的特征向量。
具体步骤如下:
1.选择一个初始向量x0,通常为一个随机向量。
2.计算xn+1 = A^-1xn,其中A为给定矩阵,A^-1为A的逆矩阵。
3.对向量xn进行归一化,得到单位向量xn+1。
4.计算新的特征值逼近λn+1= xn+1^T A xn+1。
5.判断λn+1和λn的差值是否满足收敛条件,如果满足则结束迭代,否
则将λn+1作为新的逼近特征值,重复步骤2-5。
逆幂迭代法相比其他数值方法的优点在于,它可以快速地找到近似特征值,并且收敛速度较快。
然而,逆幂迭代法也存在一些局限性,例如对初始特征值逼近的选择较为敏感,不同的初始值可能导致不同的结果。
此外,逆幂迭代法在求解重复特征值时可能会出现困难。
对于非对称矩阵的特征值问题,逆幂迭代法可能无法收敛。
以上信息仅供参考。
矩阵特征值快速求法
矩阵特征值快速求法矩阵特征值是矩阵分析中十分重要的概念。
它在物理、工程、数学等许多领域都有着广泛的应用。
矩阵特征值是指矩阵运动时特殊的运动状态,是一种宏观量度矩阵运动的指标。
求解矩阵特征值是一项复杂的任务,通常需要使用高级算法来完成。
本文将介绍几种常用的求解矩阵特征值的算法,其中包括幂法、反幂法、QR算法、分裂Broyden算法等。
一、幂法幂法是求解矩阵特征值的一种基础算法,其基本思想是通过迭代来逐步逼近矩阵的最大特征值。
幂法的核心公式如下:x_(k+1)=A*x_k/||A*x_k||其中,x_k表示第k次迭代中得到的特征向量,A表示原始矩阵。
幂法通过不断的迭代来逼近A的最大特征值,当迭代次数趋近于无限大时,得到的特征向量就是A的最大特征值所对应的特征向量。
幂法的运算量较小,适用于比较简单的矩阵。
反幂法与幂法类似,不同之处在于每次迭代时采用的是A的逆矩阵来进行计算。
其核心公式如下:x_(k+1)=(A-λI)^(-1)*x_k其中,λ表示要求解的特征值。
反幂法能够求解非常接近于特征值λ的特征向量,并且对于奇异矩阵同样适用。
需要注意的是,在实际计算中,如果A-λI的秩不满,那么反幂法就无法使用。
三、QR算法1. 将原矩阵A进行QR分解,得到A=Q*R。
2. 计算A的近似特征矩阵A1=R*Q。
5. 重复步骤3-4,直到A的对角线元素全部趋近于所求特征值为止。
QR算法的计算量较大,但其具有收敛速度快、精度高等优点,广泛应用于科学计算中。
四、分裂Broyden算法分裂Broyden算法是QR算法的一种改进算法,其基本思想是将矩阵分解成上下三角形式,然后再对其进行QR分解,以减少QR算法中的乘法运算量。
具体实现过程如下:2. 构造一个倒数矩阵B=U^(-1)*L^(-1)。
4. 计算A的近似特征矩阵A1=Q^(-1)*L^(-1)*A*R^(-1)*U^(-1)*Q。
分裂Broyden算法的计算量较小,能够有效地解决QR算法中的乘法运算量过大的问题。
计算方法之计算矩阵的特征值和特征量
计算方法之计算矩阵的特征值和特征量计算矩阵的特征值和特征向量是线性代数中的一个重要问题,它在科学研究和工程应用中有着广泛的应用。
本文将介绍计算矩阵特征值和特征向量的方法,包括特征方程法、幂法、反幂法和QR方法。
一、特征值和特征向量的定义给定一个n阶方阵A,如果存在一个非零向量x和一个标量λ,满足以下方程:Ax=λx其中,x被称为A的特征向量,λ被称为A的特征值。
二、特征方程法特征方程法是计算矩阵特征值和特征向量的一种常用方法,其基本思想是通过求解矩阵的特征方程来求得特征值。
对于一个n阶方阵A,其特征方程为:A-λI,=0其中,I是n阶单位矩阵,A-λI,表示A-λI的行列式。
解特征方程可以得到n个特征值λ₁,λ₂,...,λₙ。
然后,将这些特征值带入原方程组(A-λI)x=0,求解线性方程组得到n个特征向量x₁,x₂,...,xₙ。
三、幂法幂法是一种通过迭代来计算矩阵最大特征值和对应的特征向量的方法。
首先,随机选择一个非零向量b₀,并进行归一化,得到单位向量x₀=b₀/,b₀。
然后,通过迭代的方式,计算xₙ₊₁=Axₙ,其中xₙ为第k次迭代得到的向量。
在迭代过程中,向量xₙ的模长会逐渐趋近于最大特征值对应的特征向量。
当迭代收敛后,xₙ就是矩阵A的最大特征值对应的特征向量。
四、反幂法反幂法是一种通过迭代来计算矩阵最小特征值和对应的特征向量的方法。
首先,随机选择一个非零向量b₀,并进行归一化,得到单位向量x₀=b₀/,b₀。
然后,通过迭代的方式,计算xₙ₊₁=(A-σI)⁻¹xₙ,其中σ为待求的特征值。
在迭代过程中,向量xₙ的模长会逐渐趋近于特征值σ对应的特征向量。
当迭代收敛后,xₙ就是矩阵A的特征值为σ的特征向量。
五、QR方法QR方法是一种通过迭代来计算矩阵特征值和特征向量的方法。
首先,将矩阵A进行QR分解,得到矩阵A=QR,其中Q是正交矩阵,R是上三角矩阵。
然后,计算矩阵B=RQ,重复以上步骤,直到矩阵B收敛。
矩阵特征值求法的十种求法(非常经典)
矩阵特征值求法的十种求法(非常经典)以下是矩阵特征值求法的十种经典求法:1. 幂法(Power Method)幂法(Power Method)幂法是求解特征值的常用方法之一。
它基于一个重要的数学原理:对于一个非零向量$x$,当它连续乘以矩阵$A$的$k$次幂后,$Ax$的方向将趋于特征向量相应的特征值。
这种方法通常需要进行归一化,以防止向量过度增长。
2. 反幂法(Inverse Power Method)反幂法(Inverse Power Method)反幂法是幂法的一种变体。
它通过计算矩阵$A$的逆来求解最小的特征值。
使用反幂法时,我们需要对矩阵$A$进行LU分解,以便更高效地求解线性方程组。
3. QR方法QR方法QR方法是一种迭代方法,可以通过将矩阵$A$分解为$QR$形式来逐步逼近特征值。
这种方法是通过多次应用正交变换来实现的,直到收敛为止。
QR方法不仅可以求解特征值,还可以求解特征向量。
4. Jacobi方法Jacobi方法Jacobi方法是一种迭代方法,通过施加正交相似变换将矩阵逐步变为对角矩阵。
在每个迭代步骤中,Jacobi方法通过旋转矩阵的特定元素来逼近特征值。
这种方法适用于对称矩阵。
5. Givens旋转法Givens旋转法Givens旋转法是一种用于特征值求解的直接方法。
它通过施加Givens旋转矩阵将矩阵逐步变为对角矩阵。
这种方法是通过旋转矩阵的特定元素来实现的。
6. Householder变换法Householder变换法Householder变换法是一种用于特征值求解的直接方法。
它通过施加Householder变换将矩阵逐步变为Hessenberg形式,然后再进一步将其变为上三角形式。
这种方法是通过对矩阵的列向量进行反射来实现的。
7. Lanczos方法Lanczos方法Lanczos方法是一种迭代方法,用于对称矩阵的特征值求解。
该方法创建一个Krylov子空间,并使用正交投影找到最接近特征值的Krylov子空间中的特征值。
反幂法用来计算矩阵按模最小的特征值及其特征向量
(1)
lim
k
uk
xj
,
max (x j )
(2)
lim
k
max(
vk
)
j
1
p
,即
6
p
1 max (vk
)
j
(当k ),
且收敛速度由比值 r j p / mi确inj 定.i p
由该定理知,对 A(其p中I 可用来计算特征向量 x .j
) p应用反j 幂法,
只要选择的 是p 的一j 个较好的近似且特征值分离情
况较好,一般 r很小,常常只要迭代一二次就可完成特征
向量的计算.
反幂法迭代公式中的 是vk通过解方程组
( A pI )vk uk1
求得的. 为了节省工作量,可以先将 A 进pI行三角分解
7
P( A pI ) LU ,
x3 (1, 1 3, 2 3)T (1, 0.73205, 0.26795)T , 由此看出 u是2 的x3相当好的近似.
特征值3 1.2679 1/ 2 1.267949013 , 的真值 3 3 3 1.26794912 .
v1 (12692, 9290.3, 3400.8)T ,
11
u1 (1, 0.73198, 0.26795)T , 由 LUv2 ,P得u1
v2 (20404, 14937, 5467.4)T , u2 (1, 0.73206, 0.26796)T ,
对3 应的特征向量是
如果矩阵 ( A 存p在I ),1 其特征值为
3
幂法,反幂法求解矩阵最大最小特征值及其对应的特征向量(DOC)
数值计算解矩阵的按模最大最小特征值及对应的特征向量一.幂法1. 幂法简介:当矩阵A 满足一定条件时,在工程中可用幂法计算其主特征值(按模最大)及其特征向量。
矩阵A 需要满足的条件为: (1) 的特征值为A i n λλλλ,0||...||||21≥≥≥>(2) 存在n 个线性无关的特征向量,设为n x x x ,...,,21 1.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); % 是原值,而非其绝对值。
反幂法用来计算矩阵按模最小的特征值及其特征向量
a(k) kn
a(k) k 1,n
a(k) n,k 1
a(k) nn
17
其中
ck
(ak( k)1,k
,,
a(k) nk
)T
Rnk ,
A1(1k ) 为
k 阶上海森伯
格阵,A2(2k ) R . (nk )(nk )
设 ck 0,于是可选择初等反射阵 Rk 使 Rkck ke1, 其中, Rk 计算公式为
13
(1) 设
a11 a12 a1n
A
a21
an1
a22
an 2
a2n
ann
a11 c1
A(1) 12
A(1) 22
,
其中 c1 (a21,, an1)T R n1,不妨设 c1 0, 否则这一 步不需要约化.
本节讨论两个问题 (1) 用初等反射阵作正交相似变换约化一般实矩阵A 为上海森伯格阵. (2) 用初等反射阵作正交相似变换约化对称矩阵A 为 对称三对角阵. 于是,求原矩阵特征值问题,就转化为求上海森伯格阵 或对称三对角阵的特征值问题.
8.3.2 用正交相似变换约化一般矩阵为上海森柏格阵
设 A (aij ) R nn . 可选择初等反射阵U1,U2 ,,Un2 使 A 经正交相似变换约化为一个上海森伯格阵.
则
A2
U1 A1U1
a11 R1c1
A(1) 12
R1
R1
A(1) 22
R1
(3.1)
15
a11
实验6反幂法求矩阵按模最小特征值
实验6反幂法求矩阵按模最小特征值反幂法(Inverse Power Method)是求解矩阵按模最小特征值和相应特征向量的一种数值方法。
其基本思想是通过不断迭代来逼近矩阵的按模最小特征值,同时得到相应的特征向量。
反幂法的步骤如下:1.选择一个初始向量x0,并进行归一化处理,即使其模长为12.计算y0=A^(-1)*x0,其中A^(-1)表示矩阵A的逆矩阵,这里的计算可以通过求解线性方程组Ax=y来实现。
3.对y0进行归一化处理,即使其模长为14.迭代计算y(k)=A^(-1)*x(k-1),其中x(k-1)表示第k-1次迭代得到的向量,同样进行归一化处理得到y(k)。
5.计算δ(k)=,y(k)-y(k-1),其中,·,表示向量的模长。
6.如果δ(k)足够小或者达到预定的迭代次数,停止迭代;否则,回到步骤4接下来我们通过一个具体的例子来展示如何使用反幂法求解一个矩阵的按模最小特征值。
假设我们有一个3x3的矩阵A,其具体数值如下:123456789我们希望求解A的按模最小特征值和相应特征向量。
首先,选择一个初始向量x0,比如[1,1,1],并进行归一化处理,得到x0=[1/√3,1/√3,1/√3]。
然后,计算y0=A^(-1)*x0。
由于矩阵A是奇异矩阵,其逆矩阵不存在,所以我们需要对A做一些处理。
我们可以通过将A变成正定矩阵来解决这个问题。
一种常见的方法是计算A的转置矩阵A^T,然后计算A^T*A的最大特征值和相应特征向量,然后取其倒数作为A的逆矩阵的一个估计。
在本例中,我们计算A^T*A得到143250327411650116182然后,计算A^T*A的最大特征值和相应特征向量,即为按模最小特征值的倒数的一个估计。
可以通过使用特征值分解或者其他方法来求解特征值和特征向量。
在本例中,我们假设最大特征值为30,对应的特征向量为[1,2,3]。
将其倒数作为A的逆矩阵的一个估计,即A^(-1)≈[1/30,1/15,1/10]然后,计算y0=A^(-1)*x0,即y0=[1/30,1/15,1/10]*[1/√3,1/√3,1/√3]^T=1/30√3+1/15√3+1/ 10√3≈0.183接下来,对y0进行归一化处理,得到y0=[0.183/,0.183,,0.183/,0.183,,0.183/,0.183,]≈[0.577,0.577,0.577]然后,迭代计算y(k)=A^(-1)*x(k-1),其中x(k-1)表示第k-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)作为初始向量进行迭代,可得出以上结果。
幂法和反幂法求矩阵特征值课程知识讲解
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
幂法-反幂法求解矩阵最大最小特征值与其对应的特征向量
数值计算解矩阵的按模最大最小特征值及对应的特征向量一.幂法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); % 是原值,而非其绝对值。
反幂法求矩阵特征值
一. 问题描述用幂法与反幂法求解矩阵特征值求n 阶方阵A 的特征值和特征向量,是实际计算中常常碰到的问题,如:机械、结构或电磁振动中的固有值问题等。
对于n 阶矩阵A ,若存在数λ和n 维向量x 满足 Ax=λx (1) 则称λ为矩阵A 的特征值,x 为相应的特征向量。
由线性代数知识可知,特征值是代数方程 |λI-A|=λn+a 1λ1-n +…+a 1-n λ+a n =0 (2)的根。
从表面上看,矩阵特征值与特征向量的求解问题似乎很简单,只需求解方程(2)的根,就能得到特征值λ,再解齐次方程组(λI-A )x=0 (3) 的解,就可得到相应的特征向量。
上述方法对于n 很小时是可以的。
但当n 稍大时,计算工作量将以惊人的速度增大,并且由于计算带有误差,方程(2)未必是精确的特征方程,自然就不必说求解方程(2)与(3)的困难了。
幂法与反幂法是一种计算矩阵主特征值及对应特征向量的迭代方法, 特别是用于大型稀疏矩阵。
这里用幂法与反幂法求解带状稀疏矩阵A[501][501]的特征值。
二. 算法设计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. 反幂法 (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)计算mk =max(v)(k), u)(k= v)(k/ mk(5)若|mk -m1-k|<ε,则停止计算(1/m k作为绝对值最小特征值nλ,u)(k作为相应的特征向量);否则置k=k+1,转(3).三.程序框图1.主程序2.子程序(1). 幂法迭代程序框图(2). 反幂法迭代程序框图四. 结果显示计算结果如下:矩阵A 的按模最大特征值为:-1.070011361487e+001 矩阵A 的按模最小特征值为:-5.557910794230e-003 矩阵A 最大的特征值为:9.724634101479e+000 矩阵A 最小的特征值为:-1.070011361487e+001与各k μ(1,2,...,39)k =最接近的ik λ(用[]V k 表示)的值如下: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+000v[ 4]=-8.652284007898e+000 u[ 4]=-8.657638843237e+000 v[ 5]=-8.0934********e+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.0661********e+000 u[ 9]=-6.104545378693e+000 v[10]=-5.585101052628e+000 u[10]=-5.593926685784e+000 v[11]=-5.114083529812e+000 u[11]=-5.0833********e+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.0410********e+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.0913********e+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 501void 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++)//保存Afor(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;}elselmin=lm;for(i=0;i<N;i++)//还原Afor(j=0;j<5;j++)A[j][i]=Q[j][i];ls=fanmifa(A);//按模最小特征值for(i=0;i<N;i++)//还原Afor(j=0;j<5;j++)A[j][i]=Q[j][i];for(k=0;k<39;k++)//计算u1-u39u[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++)//还原Afor(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++)//Lif(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和yy[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;}。
反幂法求特征值例题
反幂法求特征值例题反幂法是求解矩阵特征值中的一种迭代算法。
它通过对矩阵的反幂幂乘,来逼近特征值的倒数,进而求得特征值。
本文将通过一个例题来详细介绍反幂法的原理和应用。
假设有一个方阵A,我们需要求解其最大特征值和对应的特征向量。
我们令x为初始的特征向量估计,即x=(x₁, x₂, ..., xn)T,其中n是方阵A的阶数。
首先,我们需要将x归一化,即使得,x,₂=1,其中,x,₂表示向量x的二范数。
然后,我们通过以下迭代过程来逼近最大特征值:1.计算y=A⁻¹x,得到向量y。
2.根据y的二范数,将y归一化,即使得,y,₂=13. 计算新的特征值的近似值λ=xtAy。
接下来,我们通过多次迭代来计算特征值的近似值,并逐步逼近最大特征值。
每次迭代后,都需要将特征向量归一化。
下面我们通过一个例题来说明反幂法的应用。
例题:假设有一个矩阵A如下:A=[732][341][213]我们要求解矩阵A的最大特征值和对应的特征向量。
首先,我们选择一个初始特征向量估计,假设x₀=[111]T。
然后,我们进行反幂计算。
1.计算y=A⁻¹x₀:A⁻¹=[4/19-5/192/19][-5/1910/19-1/19][2/19-1/193/19]y=A⁻¹x₀=[4/19-5/192/19][1]=[4/19][1][-5/19][1][2/19]2.归一化y,使得,y,₂=1:y,₂=√((4/19)²+(-5/19)²+(2/19)²)=√(56/361)≈0.26归一化y,即得到y₁=(4/19)/0.265≈0.715,y₂=(-5/19)/0.265≈-0.844,y₃=(2/19)/0.265≈0.3013.计算新的特征值近似值:λ=x₀ᵀAy=[111][732][4/19]=[1][7*(4/19)+3*(-5/19)+2*(2/19)][1][7*(4/19)+3*(-5/19)+2*(2/19)][1]利用正交归一性,我们得到近似值λ≈(7*(4/19)+3*(-5/19)+2*(2/19))/(1*(4/19)+1*(-5/19)+1*(2/19))=16.242以上就是第一次迭代的结果。
数值代数中的特征值计算算法
数值代数中的特征值计算算法在数值代数中,特征值计算是一项重要的任务,它在很多领域中都有广泛的应用,如物理学、工程学和计算机科学等。
特征值计算的目标是找到一个方阵的特征值以及对应的特征向量。
在本文中,我们将介绍几种常用的特征值计算算法,并对它们进行比较和评估。
一、幂法幂法是一种最简单且最常用的特征值计算算法之一。
它的基本思想是通过迭代过程逐渐逼近矩阵的最大特征值。
具体步骤如下:1. 初始化一个非零向量x,并对其进行归一化。
2. 计算矩阵A与向量x的乘积Ax。
3. 更新向量x为Ax,并进行归一化。
4. 重复步骤2和3,直到收敛或达到预设的迭代次数。
幂法的收敛条件是向量x的变化趋于稳定,即x的模长变化小于设定的阈值。
该算法的缺点是对于矩阵存在多个特征值的情况,只能收敛到模长最大的特征值对应的特征向量。
二、反幂法反幂法是幂法的一个变种,它用于计算矩阵的最小特征值。
相比于幂法,反幂法的迭代过程中需要对矩阵A的逆进行操作。
具体步骤如下:1. 初始化一个非零向量x,并对其进行归一化。
2. 计算矩阵A的逆与向量x的乘积A^(-1)x。
3. 更新向量x为A^(-1)x,并进行归一化。
4. 重复步骤2和3,直到收敛或达到预设的迭代次数。
与幂法类似,反幂法的收敛条件也是向量x的变化趋于稳定。
反幂法常用于计算矩阵的最小特征值,但对于特征值过接近零的情况,该算法可能会发散。
三、QR算法QR算法是一种迭代算法,用于计算一个方阵的特征值。
其基本思想是通过相似变换将方阵转化为上三角矩阵,从而容易求解特征值。
具体步骤如下:1. 初始化矩阵A为原始方阵。
2. 对矩阵A进行QR分解,得到矩阵Q和上三角矩阵R。
3. 计算矩阵R与Q的乘积QR。
4. 更新矩阵A为QR,并重复步骤2和3。
5. 当矩阵A的对角线元素收敛时,这些元素就是矩阵A的特征值。
QR算法的优点是适用于一般的方阵,并且通常具有较快的收敛速度。
但对于特征值重复且接近的情况,QR算法可能会产生不稳定的结果。
逆幂迭代法求特征值 -回复
逆幂迭代法求特征值-回复标题:逆幂迭代法求特征值逆幂迭代法是求解矩阵特征值的一种有效方法,尤其适用于大规模稀疏矩阵的特征值问题。
本文将详细阐述逆幂迭代法的基本原理、步骤以及应用。
1. 基本原理:逆幂迭代法是一种基于幂迭代法的改进算法。
我们知道,在幂迭代法中,我们通过不断对矩阵A进行乘方操作来逼近最大特征值对应的特征向量。
然而,对于小特征值,这种方法可能会很慢。
为了解决这个问题,逆幂迭代法提出了一种新的策略:通过对矩阵A的逆矩阵进行乘方操作,来逼近最小特征值对应的特征向量。
2. 步骤:逆幂迭代法的具体步骤如下:(1)初始化:选择一个非零初始向量x0,并设置迭代次数k的最大值和收敛精度ε。
(2)计算:对矩阵A的逆矩阵A^-1进行k次乘方操作,得到A^(-k)。
(3)更新:用A^(-k)乘以初始向量x0,得到新的向量xk。
(4)检查:如果xk - x(k-1) < ε,则认为找到了最小特征值对应的特征向量,停止迭代;否则,返回步骤(2),继续迭代。
(5)计算特征值:找到最小特征值对应的特征向量后,可以通过内积的方法计算出该特征值。
具体做法是,令λ= xT A x / xTx,其中x就是最小特征值对应的特征向量。
3. 应用:逆幂迭代法在许多领域都有重要应用。
例如,在化学反应动力学中,逆幂迭代法可以用来求解反应速率常数;在结构力学中,逆幂迭代法可以用来求解结构的固有频率和振型;在网络科学中,逆幂迭代法可以用来求解网络的PageRank值等等。
总结起来,逆幂迭代法是一种有效的求解矩阵特征值的方法,特别是对于那些需要寻找最小特征值的问题。
虽然逆幂迭代法的计算量较大,但是由于其优秀的收敛性能,使其在实际应用中仍然具有很高的价值。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
西华数学与计算机学院上机实践报告
课程名称:计算方法A
年级:2010级 上机实践成绩: 指导教师:严常龙
姓名:李国强 上机实践名称:反幂法求矩阵按模最小特征值
学号:362011********* 上机实践日期:2013.12.18
上机实践编号:6
上机实践时间:14:00 一、目的
1.通过本实验加深对反幂法的构造过程的理解;
2.能对反幂法提出正确的算法描述编程实现,得到计算结果。
二、内容与设计思想 自选方阵,用反幂法求解其按模最小特征值。
可使用实例:
⎪⎪⎪⎭
⎫ ⎝⎛---=90688465441356133A
三、使用环境
操作系统:Win 8
软件平台:Visual C++ 6.0
四、核心代码及调试过程
#include<stdio.h>
#include<math.h>
#define MAX_N 20 //矩阵最大维数
#define MAXREPT 100
#define epsilon 0.00001 //求解精度
int main()
{
int n;
int i,j,k;
double xmax,oxmax;
static double a[MAX_N][MAX_N];
static double l[MAX_N][MAX_N],u[MAX_N][MAX_N];
static double x[MAX_N],nx[MAX_N];
printf("\n 请输入矩阵阶数n:"); //输入矩阵维数
scanf("%d",&n);
if(n>MAX_N)
{
printf("the input n is larger than MAX_N,please redefine the MAX_N.\n");
return 1;
}
if(n<=0)
{
printf("please input a number between 1 and %d.\n",MAX_N);
return 1;
}
//输入A矩阵
printf("请输入矩阵的值a[i][j] i,j=0...%d;\n",n-1);
for(i=0;i<n;i++)
for(j=0;j<n;j++)
scanf("%lf",&a[i][j]);
for(i=0;i<n;i++)
x[i]=1;
oxmax=0;
for(i=0;i<MAXREPT;i++)
{
for(j=0;j<n;j++) //幂乘
{
nx[j]=0;
for(k=0;k<n;k++)
nx[j]+=a[j][k]*x[k];
}
xmax=0.0;
for(j=0;j<n;j++) //规范化
if(fabs(nx[j])>xmax)
xmax=fabs(nx[j]);
for(j=0;j<n;j++)
nx[j]/=xmax;
for(j=0;j<n;j++)
x[j]=nx[j];
if(fabs(xmax-oxmax)<epsilon)
{
printf("solve...max lamda=%lf\n",xmax); //输出
printf("the vector is:\n");
for(i=0;i<n;i++)
printf("%f\n",nx[i]);
break;
//return 0;
}
oxmax=xmax;
}
//printf("after %d repeat ,max no result ...\n",MAXREPT); for(i=0;i<n;i++)
u[i][i]=1; //U矩阵对角元为
for(k=0;k<n;k++)
{
for(i=k;i<n;i++) //计算L矩阵
{
l[i][k]=a[i][k];
for(j=0;j<=k-1;j++)
l[i][k]-=(l[i][j]*u[j][k]);
}
for(j=k+1;j<n;j++) //计算U矩阵
{
u[k][j]=a[k][j];
for(i=0;i<=k-1;i++)
u[k][j]-=(l[k][i]*u[i][j]);
u[k][j]/=l[k][k];
}
}
for(i=0;i<n;i++)
x[i]=1;
for(i=0;i<MAXREPT;i++)
{
for(j=0;j<n;j++) //反幂乘
{
nx[j]=x[j];
for(k=0;k<=j-1;k++)
nx[j]-=l[j][k]*nx[k];
nx[j]/=l[j][j];
}
for(j=n-1;j>=0;j--)
{
x[j]=nx[j];
for(k=j+1;k<n;k++)
x[j]-=u[j][k]*x[k];
}
xmax=0.0;
for(j=0;j<n;j++)
if(fabs(x[j])>xmax)
xmax=fabs(x[j]);
for(j=0;j<n;j++)
x[j]/=xmax;
if(fabs(xmax-oxmax)<epsilon)
{
printf("solve... min lamda=%lf\n",1/xmax); //输出
printf("the vector is:\n");
for(i=0;i<n;i++)
printf("%f\n",x[i]);
break;
//return 0;
}
oxmax=xmax;
}
return 1;
}
显示结果:
五、总结
本次试验利用C语言实现了用反幂法求解其按模最小特征值,通过本实验加深对反幂法的构
造过程的理解。
对于利用C语言实现该算法有了进一步的掌握,为今后的深入学习打下坚实基础。
六、附录
《数值计算方法与算法(第二版)》,张韵华主编,2012.1。