求矩阵特征值算法及程序
普通矩阵特征值的QR算法
普通矩阵特征值的QR 算法摘 要求矩阵的特征值有多种不同的办法，本文主要介绍用QR 法求矩阵的特征值，QR 法是目前求中等大小矩阵全部特征值的最有效方法之一，使用于求实矩阵或复矩阵的特征值，它和雅可比法类似，也是一种变换迭代法。
关键词：QR 分解 迭代序列 特征值 Matlab一 、QR 方法的理论：对任意一个非奇异矩阵（可逆矩阵）A ，可以把它分解成一个正交阵Q 和一个上三角阵的乘积，称为对矩阵A 的QR 分解，即A=QR 。
如果规定R 的对角元取正实数，这种分解是唯一的。
若A 是奇异的，则A 有零特征值。
任取一个不等于A 的特征值的实数μ，则A-μI 是非奇异的。
只要求出A-μI 的特征值和特征向量就容易求出矩阵A 的特征值和特征向量，所以假设A 是非奇异的，不是一般性。
设A=A 1 ，对A 1 作QR 分解，得A 1 = Q 1R 1，，交换该乘积的次序，得A 2 = R 1Q 1=，由于Q 1正交矩阵，A 1到A 2的变换为正交相似变换，于是A 1和A 2就有相同的特征值。
一般的令A 1=A ，对k=1，2，3，…..⎩⎨⎧==+)()(1迭代定义分解kk k k k k Q R A QR R Q A这样，可得到一个迭代序列{A k }，这就是QR 方法的基本过程。
二、QR 方法的实际计算步骤Householder A Hessenberg B -----→=用阵作正交相似变换上第阵一步............*::::***⎛⎫⎪ ⎪ ⎪* ⎪** ⎪ ⎪ ⎪ ⎪*⎝⎭Householder 变换：如果 v 给出为单位向量而 I 是单位矩阵，则描述上述线性变换的是 豪斯霍尔德矩阵 （v * 表示向量 v 的共轭转置）H=I -2VV*1k k kGiven k k k B Q R B B R Q +=⎧-------→→⎨=⎩用变换产生迭代序列第二步12***n λλλ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦Householder A B -----→=用阵作正交相似变换（对称阵）三对角阵*****⎛⎫ ⎪ ⎪⎪*⎪**⎪ ⎪* ⎪ ⎪*⎝⎭三、化一般矩阵为Hessenberg 阵称形如11121112122212323331n n n n n nn nn h h h h h h h h h h h H h h ---⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦的矩阵为上海森堡（Hessenberg ）阵。
第8章 矩阵特征值计算
第八章 矩阵特征值计算1 特征值性质和估计工程实践中有许多种振动问题，如桥梁或建筑物的振动，机械机件的振动，飞机机翼的颤动等，这些问题的求解常常归纳为求矩阵的特征值问题。
另外，一些稳定分析问题及相关问题也可以转化为求矩阵特征值与特征向量的问题。
1.1 特征值问题及性质设矩阵n n ⨯∈A R （或n n ⨯C ），特征值问题是：求C λ∈和非零向量n R ∈x ，使λ=Ax x (1.1)其中x 是矩阵A 属于特征值λ的特征向量。
A 的全体特征值组成的集合记为sp()A 。
求A 的特征值问题(1.1)等价于求A 的特征方程()det()0p I λλ=-=A (1.2)的根。
因为一般不能通过有限次运算准确求解()0p λ=的根，所以特征值问题的数值方法只能是迭代法。
反之，有时为了求多项式111()n n n n q a a a λλλλ--=++++的零点，可以把()q λ看成矩阵123101010n a a a a ----⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦的特征多项式（除(1)n -因子不计）。
这是一个Hessenberg 矩阵，可用QR 方法求特征值，从而求出代数方程()0q λ=的根。
矩阵特征值和特征向量的计算问题可分为两类：一类是求矩阵A 的全部特征值及其对应的向量；另一类是求部分特征值（一个或几个、按模最大或最小）及其对应的特征向量。
本章介绍部分特征值和特征向量的幂法、内积法；求实对称矩阵全部特征值的雅可比法、Given 方法和Householder 方法；求任意矩阵全部特征值的QR 算法。
在第5章已给出特征值的一些重要性质，下面再补充一些基本性质。
定理1 设n n R ⨯∈A ，则(1) 设λ为A 的特征值，则λμ-为μ-A I 的特征值；(2) 设12,,,n λλλ是A 的特征值，()p x 是一多项式，则矩阵()p A 的特征值是12(),(),,()n p p p λλλ。
求矩阵特征值算法及程序
求矩阵特征值算法及程序简介1.幂法1、幂法规范化算法（1）输入矩阵A 、初始向量)0(μ，误差eps ；（2）1⇐k ； （3）计算)1()(-⇐k k A Vμ；（4）)max (,)max ()1(1)(--⇐⇐k k k k V m V m ；（5）k k k m V /)()(⇐μ；（6）如果eps m m k k <--1,则显示特征值1λ和对应的特征向量)1(x ),终止；（7）1+⇐k k ,转（3）注：如上算法中的符号)max(V 表示取向量V 中绝对值最大的分量。
本算法使用了数据规范化处理技术以防止计算过程中出现益出错误。
2、规范化幂法程序 Clear[a,u,x];a=Input["系数矩阵A="];u=Input["初始迭代向量u(0)="]; n=Length[u];eps=Input["误差精度eps ="];nmax=Input["迭代允许最大次数nmax="]; fmax[x_]:=Module[{m=0,m1,m2}, Do[m1=Abs[x[[k]]];If[m1>m,m2=x[[k]];m=m1], {k,1,Length[x]}]; m2] v=a.u;m0=fmax[u]; m1=fmax[v];t=Abs[m1-m0]//N; k=0;While[t>eps&&k<nmax, u=v/m1; v=a.u; k=k+1;m0=m1;m1=fmax[v];t=Abs[m1-m0]//N;Print["k=",k," 特征值=",N[m1,10]," 误差=",N[t,10]]; Print[" 特征向量=",N[u,10]]]; If[k ≥nmax,Print["迭代超限"]]说明：本程序用于求矩阵A 按模最大的特征值及其相应特征向量。
矩阵特征值快速求法
矩阵特征值快速求法矩阵特征值是矩阵分析中十分重要的概念。
它在物理、工程、数学等许多领域都有着广泛的应用。
矩阵特征值是指矩阵运动时特殊的运动状态，是一种宏观量度矩阵运动的指标。
求解矩阵特征值是一项复杂的任务，通常需要使用高级算法来完成。
本文将介绍几种常用的求解矩阵特征值的算法，其中包括幂法、反幂法、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分解求矩阵特征值C语言程序
QR分解求矩阵特征值C语言程序QR分解是一种常用的矩阵分解方法，其中Q是一个正交矩阵，R是一个上三角矩阵。
利用QR分解可以求解矩阵的特征值。
以下是一个使用C语言实现QR分解求矩阵特征值的程序示例：#include <stdio.h>#include <math.h>#define N 3 // 定义矩阵的维度void eigenvalues(double A[N][N], double lambda[N]);int maindouble A[N][N] = {{2, -1, 0}, {-1, 2, -1}, {0, -1, 2}}; // 待求解的矩阵double Q[N][N], R[N][N]; // 存储QR分解结果double lambda[N]; // 存储特征值eigenvalues(A, lambda); // 求解特征值printf("矩阵的特征值为：\n");for (int i = 0; i < N; i++)printf("%.3f ", lambda[i]);}printf("\n");return 0;//初始化Q矩阵为单位阵for (int i = 0; i < N; i++)for (int j = 0; j < N; j++)if (i == j)Q[i][j]=1;elseQ[i][j]=0;}}//进行QR迭代for (int k = 0; k < N - 1; k++) //计算Ak矩阵的第k列的范数double norm = 0;for (int i = k; i < N; i++) norm += A[i][k] * A[i][k];}norm = sqrt(norm);//计算v向量double v[N] = {0};v[k] = A[k][k] + norm;for (int i = k + 1; i < N; i++)v[i]=A[i][k];}//计算P矩阵double P[N][N];for (int i = 0; i < N; i++)for (int j = 0; j < N; j++)if (i == j)P[i][j] = 1 - 2 * v[i] * v[j] / (norm * norm); elseP[i][j] = -2 * v[i] * v[j] / (norm * norm);}}//更新A和Q矩阵double temp[N][N] = {0};for (int i = 0; i < N; i++)for (int j = 0; j < N; j++)for (int m = 0; m < N; m++)temp[i][j] += P[i][m] * A[m][j];}}}for (int i = 0; i < N; i++)for (int j = 0; j < N; j++)A[i][j] = temp[i][j];Q[i][j]=Q[i][j]*P[i][j];}}}//计算矩阵Rfor (int i = 0; i < N; i++)for (int j = i; j < N; j++)R[i][j]=A[i][j];}}void eigenvalues(double A[N][N], double lambda[N])//通过矩阵R的对角线元素即可得到矩阵A的特征值for (int i = 0; i < N; i++)lambda[i] = A[i][i];}请注意，上述程序是一个简化版本的QR分解方法，并不适用于所有的矩阵。
和法求矩阵最大特征值
和法求矩阵最大特征值是一种常用的算法，用于求解大规模矩阵的最大特征值及其对应的特征向量。
这种方法的基本思想是将矩阵分解为若干个子矩阵，并将这些子矩阵的和作为输入，通过求解子矩阵的最大特征值来逐步逼近原矩阵的最大特征值。
具体来说，和法求矩阵最大特征值的基本步骤如下：
1. 将原矩阵分解为若干个子矩阵，通常选择主对角线上的子矩阵作为基本子矩阵。
2. 对于每个基本子矩阵，使用适当的算法（如Jacobi方法、SOR方法等）求解其最大特征值。
3. 将所有基本子矩阵的最大特征值相加得到一个近似值，该值即为原矩阵的最大特征值的近似值。
4. 重复步骤2和3，直到达到预定的精度要求或达到最大迭代次数。
在求解过程中，需要注意以下几点：
* 算法的收敛性：和法求矩阵最大特征值算法需要保证收敛到真实解，因此需要选择合适的算法和参数设置。
* 算法的稳定性：算法需要保证稳定运行，避免出现数值不稳定的情况。
* 矩阵分解的精度：分解的子矩阵大小和数量会影响到求解的精度和速度，需要根据实际情况进行选择。
总体来说，和法求矩阵最大特征值是一种常用的算法，适用于大规模矩阵的特征值求解问题。
通过将矩阵分解为若干个子矩阵并逐步逼近原矩阵的最大特征值，可以获得相对准确的解。
但是，该方法需要较长时间和计算资源，因此在处理大规模问题时需要权衡精度和效率。
第7章 计算矩阵的特征值和特征向量
A (7)
0 0 2.125825 = 0 8.388761 0.000009 0 0.000009 4.485401
从而A的特征值可取为 A
λ1≈2.125825, λ2≈8.388761, λ3≈4.485401
V
下面分析吉文斯变换作用到对称矩阵后正 交相似的变换效果。
注:
bpp = (app cosϑ − aqp sinϑ) cosϑ − (apq cosϑ − aqq sinϑ) sinϑ
注:
1 − s
2 t 1 − t
2
= 0 ⇒
t
2
+ 2 st
− 1 = 0
雅可比方法就是对A连续施行以上变换的方法。 取p,q使 a pq = max aij
由于求解高次多项式的根是件困难的事上述方法一般无法解出阶数略大n4的矩阵特征值的精确解在实际计算中难以按定义计算矩阵特征值
第7章 计算矩阵的特征值和特征向量
在线性代数中，一个n阶矩阵A，若有数λ及非零n维向量 v满足Av=λv，则称λ为A的特征值，v为属于特征值λ的特征 向量。在线性代数中，先计算矩阵A的特征多项式，即计算 det(λ I - A)= λn +…+(-1)ndetA的根，算出A的n个特征值λ i， i=1，2， …，n。然后解线性方程组(A- λiI)v=0，计算出对 应于λ i的特征向量。由于求解高次多项式的根是件困难的事， 上述方法一般无法解出阶数略大(n>4)的矩阵特征值的精确解， 在实际计算中难以按定义计算矩阵特征值。 本章介绍一些简单有效的计算矩阵特征值和特征向量 的近似值的数值方法。
7 ．2
反 幂 法
反幂法
Av = λv ⇒ A v =
矩阵特征值的计算
物理、力学和工程技术中的许多问题在数学上都归结为求矩 阵的特征值和特征向量问题。
� 计算方阵 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 的按模最大的特征值。因此，当给
QR迭代法求矩阵特征值
本文用QR迭代法求解矩阵A的特征值:第一步先用豪斯荷尔德变换将矩阵A化为上海 森伯格矩阵AH ，第二步再对AH 进行QR迭代(使用吉文斯变换)，当迭代满足精度要求时 输出n个特征值。
1
问题描述
对于给定的矩阵A(存在n个不相等的实特征值),采用QR方法求其n个特征值。本文中的矩 阵A由文件”gr9 009 00c rg.mm”给出。并以.mm格式描述。
2 AHi,i 2 + AHi +1,i
si = 于是，第i行和第i+1行的元素变换为：
AHi+1,i
2 + AH 2 AHi,i i+1,i
AHi,j = ci ∗ AHi,j + si ∗ AHi+1,j AHi+1,j = ci ∗ AHi+1,j − si ∗ AHi,j j = 1, 2, · · · · · · , n . 对i从1至n-1,重复以上过程,则等效实现QR分解，并得到上三角矩阵R。 在第二个过程中，R矩阵右乘上Q矩阵，即： AHnew = RQ = RG(1, 2; θ1 )T · · · · · · G(n − 2, n − 1; θn−1 )T G(n − 1, n; θn−1 )T 等效于依次对R矩阵的第i列和第i+1列进行列变换。变换后的元素为： Rj,i = ci ∗ Rj,i + si ∗ Rj,i+1 Rj,i+1 = ci ∗ Rj,i+1 − si ∗ Rj,i j = 1, 2, · · · · · · , n . 对i从1至n-1,重复以上过程,则等效地实现AHnew = RQ过程。 至此，完成了一次QR迭代。 本函数(AHnew=QR-iteration-once(AH,n))的流程图如图三所示。
矩阵特征值求解
矩阵特征值求解的分值算法12组1.1矩阵计算的基本问题(1)求解线性方程组的问题•即给定一个n阶非奇异矩阵A和n维向量b，求一个n维向量X,使得Ax =b (1. 1. 1 )(2)线性最小二乘问题，即给定一个mx n阶矩阵A和m维向量b ,求一个n维向量X,使得|AX -b| =min{ | Ay -比严R n} (1.1.2 )(3)矩阵的特征问题，即给定一个n阶实(复)矩阵A,求它的部分或全部特征值以及对应的特征向量，也就是求解方程(1. 1. 3 )一对解(4 X),其中R(C), x- R n(C n),即A为矩阵A的特征值，X为矩阵Ax = ZxA的属于特征值A的特征向量。
在工程上，矩阵的特征值具有广泛的应用，如大型桥梁或建筑物的振动问题：机械和机件的振动问题；飞机机翼的颤振问题；无线电电子学及光学系统的电磁振动问题;调节系统的自振问题以及声学和超声学系统的振动问题•又如天文、地震、信息系统、经济学中的一些问题都与矩阵的特征值问题密切相关。
在科学上，计算流体力学、统计计算、量子力学、化学工程和网络排队的马尔可夫链模拟等实际问题，最后也都要归结为矩阵的特征值问题.由于特征值问题在许多科学和工程领域中具有广泛的应用，因此对矩阵的特征值问题的求解理论研究算法的开发软件的制作等是当今计算数学和科学与工程计算研究领域的重大课题，国际上这方面的研究工作十分活跃。
1.2矩阵的特征值问题研究现状及算法概述对一个nxn阶实(复)矩阵A,它的特征值问题，即求方程(1.1.3)式的非平凡解，是数值线性代数的一个中心问题•这一问题的内在非线性给计算特征值带来许多计算问题•为了求(1.1.3)式中的A ,—个简单的想法就是显式地求解特征方程det (A 一几I)二0 (121 ) 除非对于个别的特殊矩阵，由于特征方程的系数不能够用稳定的数值方法由行列式的计算来求得，既使能精确计算出特征方程的系数，在有限精度下，其特征多项式f〃)二det(A-ZJ)的根可能对多项式的系数非常敏感能•因此，这个方法只在理论上是有意义的，实际计算中对一般矩阵是不可行的数 _ . _ . 人较大，则行列式det (A -几I)的计算量将非常大；其次，根据•首先，右矩即AfbJ阳数大于四的多项式求根不存在一种通用的方法，基于上述原Galois理论对于次因，人们只能寻求其它途径•因此，如何有效地！精确地求解’矩阵特征值问题，就成为数值线性代数领域的一个中心问题.目前，求解矩阵特征值问题的方法有两大类：一类称为变换方法，另一类称为向量迭代方法•变换方法是直接对原矩阵进行处理，通过一系列相似变换，使之变换成 一个易于求解特征值的形式，如Jacobi 算法,Givens 算法,QR 算法等。
求矩阵特征值算法及程序简介
k 1) 2a (pq
a
( k 1) pp
a
( k 1) qq
,
4
确定旋转角 ，获得旋转矩阵 J p, q, ；
（2.3） a pj a pj
(k ) (k )
( k 1)
( k 1) k) k) cos aqj sin , a (jp a (pj ( k 1) k) (k ) sin aqj cos , a (jp aqj
说明 本程序用于求对称矩阵 A 的所有特征值及其相应特征向量。程序执行后，先通过 键盘输入矩阵 A 、矩阵阶数 n、精度控制 eps 和迭代允许最大次数 n max ，程序即可给出 每次迭代的次数和对应的迭代特征值、特征向量及误差序列。其中最后输出的结果即为 所求的特征值和特征向量。如果迭代超出 n max 次还没有求出满足精度的根则输出迭代 超限提示。此外，输出的特征值矩阵可以不是真正的对角矩阵，但它们的主对角元素就 是满足要求的所有特征值。 程序中变量说明 a：存放矩阵 A 及其相似变换过程中的 Ak bb:存放特征向量矩阵 J 的转置 v:存放迭代过程中的向量 V ( k ) m1:存放所求特征值和迭代过程中的近似特征值 nmax:存放迭代允许的最大次数 eps:存放误差精度 k:记录迭代次数 t1,mu,s,c,ea,p,q,m:临时变量 a1,qq,pp：临时向量 3、例题与实验
133 6 135 例 1. 用幂法求矩阵 A 44 5 46 88 6 90
的按模最大的特征值及其相应特征向量，要求误差 eps 10 4 。
解：执行幂法程序后在输入的四个窗口中按提示分别输入: {{133,6,135},{44,5,46},{-88,-6,-90}}、{1,1,1}、0.0001、20 每次输入后用鼠标点击窗口的“OK”按扭，得如下输出结果：
矩阵特征值和特征向量计算.ppt
j
＝1
1
1
1
j
i
n 2
i
i 1
k
1
j
i
n 2
i
i 1
k
1
i
j
i
j
( 4.2)
lim
uk
j
k
uk1
j
1 ,
故k充
分
大
时
， uk
j
uk1
j
1 ，
(j
1,2,, n)
由（4.1）显然知k充分大时， 0 ,
x 故 uk ( 1k1 1 )就 是1对 应 的 近 似 特 征 向 量 。
v u v u u 如用
m
m
或 m
m
代替 继续迭代, m
u( )m max
(u ) min m
u u u 这里(
m )max 和(
m )min 分 别 表 示 向 量(
)的 绝 对 值
m
最 大 的 分 量 和 最 小 分 量;
4. 由(4.1),乘 幂 法 的 速 度 与 比 值| 2 | 有 关, 1
n
A1
x
1
x
一 定 是A1的
按
模
最
大
的
征值，故对A1用乘幂法— 反幂法，可得1 的近似值
算法（步1）骤：u0 0
n
( 2) (3)
计 算u k
1 A uk 1
(k 1,2,3,)
u 若k充分大后 ( u(
k)j c, ) k 1 j
则n
1 ,
c
uk
就
是n
注：实际相计对算应: A的u特征u向量三。角分解A LU ,
幂法求特征值算法
幂法求特征值算法幂法是一种用于求解矩阵特征值和特征向量的迭代算法。
它是数值线性代数中的重要算法之一，在科学和工程领域有着广泛的应用。
幂法的基本思想是利用一个非零向量的矩阵幂序列逐渐逼近矩阵的主特征向量。
主特征向量对应矩阵的最大特征值，因此通过逼近主特征向量，我们也能够得到矩阵的最大特征值。
算法步骤如下：1.随机选择一个非零向量b作为初始向量。
2.计算矩阵A乘以向量b的结果，得到向量b1，即b1=A*b。
3.对向量b1进行归一化，使其成为单位向量，即b1=b1/||b1||，其中||b1||表示向量b1的模。
4.判断新向量b1与旧向量b之间的误差是否足够小，如果误差小于给定的阈值，则停止计算，否则继续迭代。
5.将新向量b1赋值给旧向量b，即b=b1。
6.重复步骤2-5，直到满足停止条件。
幂法的核心思想是利用矩阵与向量的乘积进行迭代，使得向量收敛到矩阵的特征向量。
由于迭代过程中，向量b的每一次迭代都会趋近于主特征向量，因此迭代的次数越多，结果越接近主特征向量。
最终，我们可以得到矩阵的主特征向量以及对应的特征值。
幂法的收敛性取决于特征向量对应的特征值的相对大小。
如果矩阵的特征值不重复且特征值的模最大，那么幂法能够收敛到主特征向量。
在实际应用中，为了加快收敛速度，通常会对矩阵进行特征值的平移，使得矩阵的主特征值接近于零，然后再进行幂法的迭代计算。
幂法求解特征值的时间复杂度为O(n^2m)，其中n为矩阵的维度，m为迭代的次数。
这是由于每一次迭代需要进行一次矩阵与向量的乘积，而矩阵与向量的乘积的时间复杂度为O(n^2)。
总结起来，幂法是一种简单而有效的求解矩阵特征值和特征向量的方法。
它通过迭代的方式，利用矩阵与向量的乘积逼近特征向量，从而得到矩阵的特征值。
幂法在科学和工程领域有广泛的应用，如电力系统中的潮流计算、结构分析、图像处理等。
虽然幂法存在一些限制，如只能求解特征模最大的特征值和对应的特征向量，但其简单性和高效性使得它成为一种常用的数值线性代数算法。
求解矩阵特征值问题的算法研究
求解矩阵特征值问题的算法研究求解矩阵特征值问题是线性代数和数值计算中的重要问题。
特征值问题的一般形式为Ax = λx，其中A是一个矩阵，x是一个非零向量，λ是一个标量。
求解特征值问题即寻找矩阵A 的特征值λ和对应的特征向量x。
以下是一些常用的求解矩阵特征值问题的算法：1. 幂迭代法：幂迭代法是求解特征值问题的一种迭代方法。
该方法的基本思想是选择一个随机的初始向量x0，通过迭代计算xk+1 = Axk / ||Axk||，其中||.||表示向量的2-范数。
随着迭代的进行，x收敛到矩阵A的主特征向量，而λ则收敛到对应的主特征值。
2. QR迭代法：QR迭代法是求解特征值问题的一种迭代方法。
该方法首先通过QR分解将矩阵A分解为Q和R的乘积，其中Q是一个正交矩阵，R是一个上三角矩阵。
然后，将分解后的矩阵重新相乘得到新的矩阵A'，迭代此过程直到A'的对角线元素收敛到矩阵的特征值。
3. 特征向量分解法：特征向量分解法是求解特征值问题的一种直接方法。
该方法通过对矩阵A 进行特征向量分解得到矩阵V和对角矩阵Λ，其中V的列向量是A的特征向量，Λ的对角线元素是A的特征值。
特征向量分解法在理论上可以得到A的所有特征值和特征向量，但实际计算中对大型矩阵会较为困难。
4. Davidson方法：Davidson方法是求解特征值问题的一种迭代方法。
该方法采用子空间迭代的方式逐步构建特征子空间，并寻找特征值和对应的特征向量。
Davidson方法可以有效地处理大型矩阵，特别适合于求解特征值问题中的一部分特征对。
这些算法都有各自的优点和适用范围，研究者可根据具体问题选择合适的算法进行求解。
幂法反幂法求解矩阵最大最小特征值及其对应的特征向量
幂法反幂法求解矩阵最大最小特征值及其对应的特征向量幂法是一种迭代算法，用于求解矩阵的最大特征值及其对应的特征向量。
反幂法是幂法的一种变体，用于求解矩阵的最小特征值及其对应的特征向量。
两种方法在求解特征值问题时有相似的步骤，但反幂法需要对矩阵进行一定的变换。
幂法的基本思想是通过不断迭代的方式逼近矩阵的最大特征值及其对应的特征向量。
求解的过程如下：1.随机选择一个初始向量x0，并进行归一化，即使其模长为12. 根据公式计算新的向量xk+1 = Axk，其中A为待求解特征值的矩阵。
3. 对xk+1进行归一化。
4. 计算矩阵A关于xk+1的雷神特征值λk+1 = (Axk+1)·xk+1 / xk+1·xk+1，其中·表示向量的内积。
5.重复步骤2至4，直到满足收敛条件。
幂法的收敛条件一般是设置一个精度，当迭代的过程中特征向量的变化小于该精度时，认为结果已经收敛。
最终得到的特征值就是矩阵A的最大特征值，对应的特征向量为收敛时的xk+1反幂法是对幂法的一种改进，用于求解矩阵的最小特征值及其对应的特征向量。
反幂法的基本思想是通过将矩阵A的特征值问题转化为矩阵B=(A-μI)^-1的特征值问题来求解，其中μ为一个非常接近待求解特征值的数。
求解的过程如下：1.随机选择一个初始向量x0，并进行归一化，即使其模长为12. 根据公式求解新的向量xk+1 = (A-μI)^-1xk，其中A为待求解特征值的矩阵，μ为一个非常接近待求解特征值的数。
3. 对xk+1进行归一化。
4. 计算矩阵B关于xk+1的雷神特征值λk+1 = (Bxk+1)·xk+1 / xk+1·xk+1，其中·表示向量的内积。
5.重复步骤2至4，直到满足收敛条件。
反幂法的收敛条件与幂法相似，一般也是设置一个精度。
最终得到的特征值就是矩阵A的最小特征值，对应的特征向量为收敛时的xk+1总结：幂法和反幂法是求解矩阵最大最小特征值的常用迭代算法。
反幂法求矩阵特征值
一． 问题描述用幂法与反幂法求解矩阵特征值求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;}。
计算方法之计算矩阵的特征值和特征量
an1
an2 ann
—— 特征多项式方程。
2
在线性代数中按如下三步计算：
1、计算出A的特征多项式│A- E│； 2、求出特征方程│A-E│=0的全部根i 3、将i代入(A-iE)X=0 求出基础解系，即得A 的对应于i的特征向量，而基础解系的线性组合即 为A的对应于i 的全部特征向量。
例
求矩阵
7 44.99953 14.99983 -29.99968 1 0.33333 -0.66667 44.99953
由表可知，最大特征值为： 1=44.99953
对应特征向量为：( 1 , 0.33333 , -0.66667 )T
26
（二）按模最大特征值1是单实根，但1<0
此时迭代向量序列{V(2k)}和{V(2k+1)}将分别收敛 于互为反号的向量。
1
定义 设A为 n 阶方阵，若存在常数 与 n 维 非零向量X 使 AX=X成立，则称 为方阵A的特 征值，非零向量 X 为A的对应于 的特征向量。
由AX=X (A- E)X=0 此方程有非零解的充要条件是: |A- E|=0 , 即：
a11 a12
a21 a22
a1n a2n 0
k
V(k)
U(k)
max(V(k))
0
1
1
1
1
1
1
1 274
95
-184
1 0.34672 -0.67153
2 44.42377 14.84322 -29.64262 1 0.33413 -0.66727 44.42377
3 44.92333 14.97623 -29.95048 1 0.33337 -0.66670 44.92333
幂迭代法求特征值
幂迭代法求特征值一、引言幂迭代法是求解矩阵特征值和特征向量的常用方法之一，其基本思想是通过不断迭代矩阵的一个向量来逼近该矩阵的最大特征值和对应的特征向量。
本文将从幂迭代法的原理、步骤和算法流程、收敛性及优缺点等方面进行详细介绍。
二、幂迭代法原理幂迭代法是基于以下定理而提出的：定理：设A为n阶方阵，λ1, λ2, …, λk为A的k个不同特征值，且满足|λ1|>|λ2|≥…≥|λk|，则对于任意非零向量x0∈R^n，当k=1时，有lim(k→∞) Akx0/||Akx0||=λ1v1，其中v1为A属于λ1的单位特征向量。
根据该定理可知，在不断迭代矩阵向量时，当k趋近于无穷大时，所得到的向量将趋近于属于最大特征值对应的单位特征向量。
因此，可以通过不断迭代来逼近最大特征值和对应的特征向量。
三、幂迭代法步骤及算法流程1. 初始化向量x0，使其满足||x0||=1；2. 进行迭代计算，即不断将向量xk乘以矩阵A，得到新的向量xk+1=A*xk；3. 对新的向量进行归一化处理，即令xk+1=xk+1/||xk+1||；4. 判断迭代是否收敛，若未收敛，则返回第二步；否则，输出结果。
幂迭代法的算法流程如下所示：```pythondef power_iteration(A, x0, epsilon, max_iter):x = x0for i in range(max_iter):y = A @ xeigenvalue = np.linalg.norm(y)x = y / eigenvalueif np.linalg.norm(A @ x - eigenvalue * x) < epsilon:breakreturn eigenvalue, x```四、幂迭代法收敛性分析幂迭代法并不保证能够收敛到最大特征值和对应的特征向量。
其收敛性与矩阵A的特征值分布有关。
当A的所有特征值都是实数且正交时，幂迭代法一定能够收敛到最大特征值和对应的特征向量。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
求矩阵特征值算法及程序简介1.幂法1、幂法规范化算法（1）输入矩阵A 、初始向量)0(μ，误差eps ； （2）1⇐k ；（3）计算)1()(-⇐k k A V μ；（4）)max (,)max ()1(1)(--⇐⇐k k k k V m V m ； （5）k k k m V /)()(⇐μ；（6）如果eps m m k k <--1,则显示特征值1λ和对应的特征向量)1(x),终止；（7）1+⇐k k ,转（3） 注：如上算法中的符号)max(V 表示取向量V 中绝对值最大的分量。
本算法使用了数据规范化处理技术以防止计算过程中出现益出错误。
2、规范化幂法程序Clear[a,u,x];a=Input["系数矩阵A="];u=Input["初始迭代向量u(0)="];n=Length[u];eps=Input["误差精度eps ="];nmax=Input["迭代允许最大次数nmax="];fmax[x_]:=Module[{m=0,m1,m2},Do[m1=Abs[x[[k]]];If[m1>m,m2=x[[k]];m=m1],{k,1,Length[x]}];m2]v=a.u;m0=fmax[u];m1=fmax[v];t=Abs[m1-m0]//N;k=0;While[t>eps&&k<nmax,u=v/m1;v=a.u;k=k+1;m0=m1;m1=fmax[v];t=Abs[m1-m0]//N;Print["k=",k," 特征值=",N[m1,10]," 误差=",N[t,10]];Print[" 特征向量=",N[u,10]]];If[k ≥nmax,Print["迭代超限"]]说明：本程序用于求矩阵A 按模最大的特征值及其相应特征向量。
程序执行后，先通过键盘输入矩阵A 、迭代初值向量)0(μ、精度控制eps 和迭代允许最大次数max n ，程序即可给出每次迭代的次数和对应的迭代特征值、特征向量及误差序列，它们都按10位有效数输出。
其中最后输出的结果即为所求的特征值和特征向量序列。
如果迭代超出max n 次还没有求出满足精度的根则输出迭代超限提示，此时可以根据输出序列判别收敛情况。
程序中变量说明a:存放矩阵A ；u:初始向量)0(μ和迭代过程中的向量)(k μ及所求特征向量；v:存放迭代过程中的向量)(k V ；m1:存放所求特征值和迭代过程中的近似特征值；nmax:存放迭代允许的最大次数；eps:存放误差精度；fmax[x]: 给出向量x 中绝对值最大的分量；k:记录迭代次数；t1:临时变量；注：迭代最大次数可以修改为其他数字。
3、例题与实验例1. 用幂法求矩阵⎪⎪⎪⎭⎫ ⎝⎛---=90688465441356133A 的按模最大的特征值及其相应特征向量，要求误差410-<eps 。
解：执行幂法程序后在输入的四个窗口中按提示分别输入:{{133,6,135},{44,5,46},{-88,-6,-90}}、{1,1,1}、0.0001、20每次输入后用鼠标点击窗口的“OK ”按扭，得如下输出结果：此结果说明迭代6次，求得误差为0.0000101442的按模最大的特征值为44.99999952,及其对应的一个特征向量:{1.000000000,0.3333333371,-0.6666666704}2.反幂法1、 反幂法规范化算法（1）输入矩阵A 、初始向量)0(μ，误差eps ； （2）1⇐k ；（3）计算)1()(-⇐k k AV μ求出解)(k V ；（4）)max (,)max ()1(1)(--⇐⇐k k k k V m V m ； （5）k k k m V /)()(⇐μ；（6）如果eps m m k k <--1,则显示特征值1λ和对应的特征向量)1(x),终止；（7）1+⇐k k ,转（3）注：如上算法中解方程)1()(-⇐k k AV μ可以使用Dololittle 分解法。
本算法使用了数据规范化处理技术以防止计算过程中出现益出错误。
2、规范化反幂法程序Clear[a,u,x];a=Input["系数矩阵A="];u=Input["初始迭代向量u(0)="];n=Length[u];eps=Input["误差精度eps ="];nmax=Input["迭代允许最大次数nmax="];fmax[x_]:=Module[{m=0,m1,m2},Do[m1=Abs[x[[k]]];If[m1>m,m2=x[[k]];m=m1],{k,1,Length[x]}];m2];v=a.u;a1=Inverse[a];m0=fmax[u];m1=fmax[v];t=Abs[m1-m0]//N;k=0;While[t>eps&&k<nmax,u=v/m1;v=a1.u;k=k+1;m0=m1;m1=fmax[v];t=Abs[m1-m0]//N;t1=Abs[1/m1-1/m0]//N;Print["k=",k," 特征值=",N[1/m1,10]," 误差=",N[t1,10]];Print[" 特征向量=",N[u,10]]];If[k ≥nmax,Print["迭代超限"]]说明：本程序用于求矩阵A 按模最小的特征值及其相应特征向量。
程序执行后，先通过键盘输入矩阵A 、迭代初值向量)0(μ、精度控制eps 和迭代允许最大次数max n ，程序即可给出每次迭代的次数和对应的迭代特征值、特征向量及误差序列，它们都按10位有效数输出。
其中最后输出的结果即为所求的特征值和特征向量序。
如果迭代超出max n 次还没有求出满足精度的根则输出迭代超限提示，此时可以根据输出序列判别收敛情况。
程序中变量说明a ：存放矩阵Au ：初始向量)0(μ和迭代过程中的向量)(k μ及所求特征向量v:存放迭代过程中的向量)(k Va1：存放逆矩阵1-Am1: 存放所求特征值和迭代过程中的近似特征值nmax:存放迭代允许的最大次数eps:存放误差精度fmax[x]: 给出向量x 中绝对值最大的分量k:记录迭代次数t1:临时变量注：迭代最大次数可以修改为其他数字。
3、例题与实验例3. 用反幂法求矩阵⎪⎪⎪⎭⎫ ⎝⎛--=131111322A 的按模最小的特征值及其相应特征向量，要求误差510-<eps 。
解：执行幂法程序后在输入的四个窗口中按提示分别输入:{{2.,-2.,3.},{1,1.,1},{1.,3,-1}}、{1,0,1}、0.00001、100每次输入后用鼠标点击窗口的“OK ”按扭，得如下输出结果：注意到本题按模最小的特征值为1，因此求解效果较满意。
3.Jacobi 方法1、Jacobi 旋转法算法1、输入矩阵A ，误差ε；2、Λ,2,1=k For ；（2.1）选取)1()1(m ax -≠-=k ij j i k pqa a 记录q p , （2.2）由4,22tan )1()1()1(πθθ≤-=---k qq k pp k pqa a a 确定旋转角θ，获得旋转矩阵()θ,,q p J ；（2.3）q p j a a a a a k pj k jp k qj k pjk pj ,,sin cos )()()1()1()(≠⇐+⇐--θθ； （2.4）q p j a a a a a k qj k jp k qj k pjk qj ,,cos sin )()()1()1()(≠⇐+-⇐--θθ； （2.5）q p j i a a k ijk ij ,,)1()(≠⇐- （2.6）θθθθcos sin 2sin cos )1(2)1(2)1()(---++⇐k pq k qq k pp k pp a a a a（2.7）θθθθcos sin 2cos sin )1(2)1(2)1()(----+⇐k pq k qq k pp k qq a a a a （2.8）计算()∑≠=j i k ij k a A E 2)()(（2.9）如果ε<)(k A E ，输出对角矩阵),,,(21n diag D λλλΛ=和特征向量矩阵J ，停止注：如上算法中)()(k ij k a A =，)()()0(0ij ij a a A ==。
2、Jacobi 算法程序Clear[a,bb];a=Input["矩阵A="];n=Input["矩阵阶数n="];eps=Input["误差精度eps ="];nmax=Input["迭代允许最大次数nmax="];k=0;bb=IdentityMatrix[n];ea=Sum[a[[i,j]]^2,{i,1,n},{j,1,n}]-Sum[a[[i,i]]^2,{i,1,n}]//N;While[ea>eps&&k<nmax,m=0;Print["迭代次数k=",k];Do[If[Abs[a[[i,j]]]>m,m=Abs[a[[i,j]]];p=i;q=j],{i,1,n},{j,i+1,n}]; mu=a[[p,p]]-a[[q,q]];If[mu 0,thi=Pi/4,thi=ArcTan[2*a[[p,q]]/mu]/2];s=Sin[thi]//N;c=Sqrt[1-s^2];a1=bb[[p]];bb[[p]]=c*bb[[p]]+s*bb[[q]];bb[[q]]=-s*a1+c*bb[[q]];pp=a[[p,p]]*c*c+a[[q,q]]*s*s+2a[[p,q]]*s*c;qq=a[[p,p]]*s*s+a[[q,q]]*c*c-2a[[p,q]]*s*c;Do[a1=a[[p,j]];a[[p,j]]=c*a[[p,j]]+s*a[[q,j]];a[[j,p]]=a[[p,j]];a[[q,j]]=c*a[[q,j]]-s*a1;a[[j,q]]=a[[q,j]],{j,1,n}];a[[p,p]]=pp;a[[q,q]]=qq;a[[p,q]]=0;a[[q,p]]=0;ea=Sum[a[[i,j]]^2,{i,1,n},{j,1,n}]-Sum[a[[i,i]]^2,{i,1,n}]//N; k=k+1;Print["误差=",ea];Print["相似矩阵A="];Print[MatrixForm[a]];Print["特征向量J"];Print[MatrixForm[Transpose[bb]]]];说明 本程序用于求对称矩阵A 的所有特征值及其相应特征向量。