求矩阵特征值算法及程序
求矩阵特征值算法及程序
求矩阵特征值算法及程序简介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 按模最大的特征值及其相应特征向量。
计算方法(5)第四章 矩阵特征值和特征向量的计算
n
使得u 0
i xi
i 1
n
n
uk Auk1 Aku0 Ak (i xi ) iik xi
i 1
i 1
1k [1x1
n i2
( i 1
)k i xi ]
由1 0, 1 i (i 2, 3,L , n) 得
lim(
对矩阵A1用乘幂法得 uk
A-1u
k
，
1
因为A1 的计算
比较麻烦，而且往往不能保持矩阵A 的一些好性质
（如稀疏性），因此，反幂法在实际计算时以求解
方程组
Auk
u
k
，代替迭代
1
uk
A-1uk1求得uk，每
迭代一次要解一线性方程组。 由于矩阵在迭代过
程中不变，故可对A 先进行三角分解，每次迭代只 要解两个三角形方程组。
且
2 p 2 n
2 n
2 n 2
1 p 21 2 n 1 n 1 2 1 n 1
因此，用原点平移法求1可使收敛速度加快。
三、反幂法
反幂法是计算矩阵按模最小的特征值及特征向 量的方法，也是修正特征值、求相应特征向量的最 有效的方法。
0
0.226
0.975
做正交相似变换后得到
3.366
A3 =R2 AR2T
0.0735
0.317
0.0735 1.780
0
0.317
0
1.145
雅可比方法是一个迭代过程，它生成的是一个矩阵的
序列 Ak，当k越大时Ak就越接近于对角矩阵，从而
矩阵特征值快速求法
矩阵特征值快速求法矩阵特征值是矩阵分析中十分重要的概念。
它在物理、工程、数学等许多领域都有着广泛的应用。
矩阵特征值是指矩阵运动时特殊的运动状态，是一种宏观量度矩阵运动的指标。
求解矩阵特征值是一项复杂的任务，通常需要使用高级算法来完成。
本文将介绍几种常用的求解矩阵特征值的算法，其中包括幂法、反幂法、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算法中的乘法运算量过大的问题。
和法求矩阵最大特征值
和法求矩阵最大特征值是一种常用的算法，用于求解大规模矩阵的最大特征值及其对应的特征向量。
这种方法的基本思想是将矩阵分解为若干个子矩阵，并将这些子矩阵的和作为输入，通过求解子矩阵的最大特征值来逐步逼近原矩阵的最大特征值。
具体来说，和法求矩阵最大特征值的基本步骤如下：
1. 将原矩阵分解为若干个子矩阵，通常选择主对角线上的子矩阵作为基本子矩阵。
2. 对于每个基本子矩阵，使用适当的算法（如Jacobi方法、SOR方法等）求解其最大特征值。
3. 将所有基本子矩阵的最大特征值相加得到一个近似值，该值即为原矩阵的最大特征值的近似值。
4. 重复步骤2和3，直到达到预定的精度要求或达到最大迭代次数。
在求解过程中，需要注意以下几点：
* 算法的收敛性：和法求矩阵最大特征值算法需要保证收敛到真实解，因此需要选择合适的算法和参数设置。
* 算法的稳定性：算法需要保证稳定运行，避免出现数值不稳定的情况。
* 矩阵分解的精度：分解的子矩阵大小和数量会影响到求解的精度和速度，需要根据实际情况进行选择。
总体来说，和法求矩阵最大特征值是一种常用的算法，适用于大规模矩阵的特征值求解问题。
通过将矩阵分解为若干个子矩阵并逐步逼近原矩阵的最大特征值，可以获得相对准确的解。
但是，该方法需要较长时间和计算资源，因此在处理大规模问题时需要权衡精度和效率。
矩阵特征值的计算
物理、力学和工程技术中的许多问题在数学上都归结为求矩 阵的特征值和特征向量问题。
� 计算方阵 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）输入矩阵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 按模最大的特征值及其相应特征向量。
特征值分解 算法流程
特征值分解算法流程特征值分解算法流程简介特征值分解是一种重要的矩阵分解方法，广泛应用于数学、物理、工程和计算机科学等领域。
它能将一个方阵分解为特征向量和特征值的乘积形式，揭示出矩阵的结构和性质。
算法流程特征值分解的算法流程可以概括为以下几个步骤：1.输入方阵：将待分解的方阵记为A。
2.求解特征方程：通过求解特征方程det(A - λI) =0，得到特征值λ。
其中，A是输入方阵，λ是待求的特征值，I是单位矩阵。
3.求解特征向量：对于每一个特征值，利用A - λI，解出特征向量。
4.归一化特征向量：对每一个特征向量进行归一化处理，使其模长为1。
5.组成特征值矩阵和特征向量矩阵：将求得的特征值按对角线排列形成特征值矩阵Λ，将求得的特征向量按列形成特征向量矩阵P。
6.完成分解：得到特征值矩阵Λ和特征向量矩阵P后，即完成了对输入方阵A的特征值分解。
示例为了更好地理解上述算法流程，以下给出一个简单的示例：输入方阵A为：[ 4 2 ][ 1 3 ]1.求解特征方程：det(A - λI) = 0=> det([4-λ 2 | 1 0 ])[ 1 3-λ | 0 1 ]= (4-λ)(3-λ) - 2·1 = λ^2 - 7λ + 10 = (λ - 2)(λ - 5) = 0得到特征值λ1 = 2 和λ2 = 5。
2.求解特征向量：对于λ1 = 2：[2-2 2 | 1 0 ][ 1 3-2 | 0 1 ]解得特征向量v1 = [1, -1]。
对于λ2 = 5：[4-5 2 | 1 0 ][ 1 3-5 | 0 1 ]解得特征向量v2 = [1, 1]。
3.归一化特征向量：归一化v1 = [1/√2, -1/√2]，归一化v2 = [1/√2, 1/√2]。
4.组成特征值矩阵和特征向量矩阵：特征值矩阵Λ为：[ 2 0 ][ 0 5 ]特征向量矩阵P为：[ 1/√2 1/√2 ][ -1/√2 1/√2 ]5.完成分解：得到特征值矩阵Λ和特征向量矩阵P后，即完成了对输入方阵A的特征值分解。
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))的流程图如图三所示。
北航数值分析计算实习题目一 幂法 反幂法 求矩阵特征值
《数值分析》计算实习题目第一题：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）作为初始向量进行迭代，可得出以上结果。
矩阵特征值求解
矩阵特征值求解的分值算法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 算法等。
特征值的求法
特征值的求法
特征值是线性代数中的一个重要概念，主要用于描述矩阵的某些重要性质。
对于方阵，特征值可以通过求解特征多项式得到。
以下是特征值的基本求法：
1.写出矩阵A的特征多项式f(λ)。
对于n阶矩阵A，其特征多项式为f(λ)=|λE-A|，其中E是n阶单位矩阵。
2.求解特征多项式f(λ)=0的根，这些根就是矩阵A的特征值。
这个方程的解可能是一个或多个实数，也可能是复数。
3.对于每个解出的特征值λ，求解齐次线性方程组(λE-A)x=0的非零解x，这个解x就是对应于特征值λ的特征向量。
以上步骤是求解特征值和特征向量的基本方法。
需要注意的是，对于具体的矩阵，可能需要根据其特点选择合适的求解方法，例如对于大型稀疏矩阵，可能需要使用迭代法等数值方法求解特征值和特征向量。
此外，对于一些特殊的矩阵，如对称矩阵、反对称矩阵、正交矩阵等，其特征值和特征向量具有一些特殊的性质，可以利用这些性质简化求解过程。
以上信息仅供参考，如需更多信息，建议查阅线性代数相关书籍或咨询专业教师。
用幂法求解矩阵特征值和特征向量
x= -0.3930 -0.9774 0.2921 1.0000 第五题 A=[-1 2 1; 2 -4 1; 1 1 -6 ]; v0=[1 1 1]'; tol=1e-4; [lda,x]=mifa(A,v0,tol) lda = -6.4209
第4页
数值分析实验指导
x= -0.0463 -0.3746 000
( 1, 0, 1, 0, 0, 1 )T 105
1 21.30525 6 1.62139 x1 0.8724,0.5401,0.9973,0.5644,0.4972,1.0 T
第1页
数值分析实验指导
2 1 1 2 1 （3） A= 1 2 1 1 2 1 1 2 T 0 104 取 =( 1, 1, 1, 1, 1 ) 参考结果： 3.7321 3 4 2 1 1 3 1 5 （4） A= 3 1 6 2 4 5 2 1 T 2 取 0 =( 1, 1, 1, 1 ) ， 10 。
第3页
数值分析实验指导
x= 0.5000 -0.8660 1.0000 -0.8660 0.5000 第四题 A=[2 1 3 4; 1 -3 1 5; 3 1 6 -2; 4 5 -2 -1 ]; v0=[1 1 1 1]'; tol=1e-2; [lda,x]=mifa(A,v0,tol) lda = -8.0136
下面再考虑主特征值 1 的的计算，用 (vk )i 表示 vk 的第 i 个分量，则
( x ) ( k 1 )i (vk 1 )i 1 1 1 i , (vk )i 1 ( x1 )i ( k )i
故
求解矩阵特征值问题的算法研究
求解矩阵特征值问题的算法研究求解矩阵特征值问题是线性代数和数值计算中的重要问题。
特征值问题的一般形式为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总结：幂法和反幂法是求解矩阵最大最小特征值的常用迭代算法。
幂法,反幂法求解矩阵最大最小特征值及其对应的特征向量
数值计算解矩阵的按模最大最小特征值及对应的特征向量一.幂法1. 幂法简介：当矩阵A 满足一定条件时，在工程中可用幂法计算其主特征值(按模最大)及其特征向量。
矩阵A 需要满足的条件为：(1) 存在n 个线性无关的特征向量，设为n x x x ,...,,211.1计算过程：i n i i i u xx αα,1)0()0(∑==，有对任意向量不全为0，则有 可见，当||12λλ越小时，收敛越快；且当k 充分大时，有1)1111)11111λαλαλ=⇒⎪⎩⎪⎨⎧==+++(k )(k k (k k )(k x x u x u x ，对应的特征向量即是)(k x 1+。
2 算法实现3 matlab 程序代码function [t,y]=lpowerA,*0,eps,N) % t 为所求特征值，y 是对应特征向量k=1;z=0; % z 相当于λy=*0./ma*(abs(*0)); % 规化初始向量*=A*y; % 迭代格式b=ma*(*); % b 相当于 βif abs(z-b)<eps % 判断第一次迭代后是否满足要求t=ma*(*);return ;endwhile abs(z-b)>eps && k<Nk=k+1;z=b;y=*./ma*(abs(*));*=A*y;b=ma*(*);end[m,inde*]=ma*(abs(*)); % 这两步保证取出来的按模最大特征值t=*(inde*); % 是原值，而非其绝对值。
end4 举例验证选取一个矩阵A ，代入程序，得到结果，并与eig(A)的得到结果比拟，再计算A*y-t*y ，验证y 是否是对应的特征向量。
结果如下：结果正确，说明算法和代码正确，然后利用此程序计算15阶Hilb 矩阵，与eig(A)的得到结果比拟，再计算 A*y-t*y ，验证y 是否是对应的特征向量。
设置初始向量为*0=ones(15,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;}。
幂迭代法求特征值
幂迭代法求特征值一、引言幂迭代法是求解矩阵特征值和特征向量的常用方法之一，其基本思想是通过不断迭代矩阵的一个向量来逼近该矩阵的最大特征值和对应的特征向量。
本文将从幂迭代法的原理、步骤和算法流程、收敛性及优缺点等方面进行详细介绍。
二、幂迭代法原理幂迭代法是基于以下定理而提出的：定理：设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. 初始化一个非零向量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、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
求矩阵特征值算法及程序简介1.幂法1、幂法规范化算法(1)输入矩阵A、初始向量( 0)，误差eps；(2) k 1；(3)计算V(k)A(k 1)；(4)m k max(V(k)) ,m k1max( V ( k 1))；(5) (k)V(k)/m k；(6)如果m k m k 1eps,则显示特征值1和对应的特征向量x(1) ),终止；(7)k k 1, 转(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和迭代允许最大次数n max ，程序即可给出每次迭代的次数和对应的迭代特征值、特征向量及误差序列，它们都按10 位有效数输出。
其中最后输出的结果即为所求的特征值和特征向量序列。
如果迭代超出nmax 次还没有求出满足精度的根则输出迭代超限提示，此时可以根据输出序列判别收敛情况。
程序中变量说明a: 存放矩阵A ；u: 初始向量(0)和迭代过程中的向量(k)及所求特征向量；v: 存放迭代过程中的向量V (k )；m1:存放所求特征值和迭代过程中的近似特征值；nmax:存放迭代允许的最大次数；eps: 存放误差精度；fmax[x]: 给出向量x 中绝对值最大的分量；k: 记录迭代次数；t1: 临时变量；注：迭代最大次数可以修改为其他数字。
3、例题与实验133 6 135例1. 用幂法求矩阵A 44 5 4688 6 90的按模最大的特征值及其相应特征向量，要求误差eps 10解：执行幂法程序后在输入的四个窗口中按提示分别输入:{{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) k 1；(3)计算AV(k) (k 1)求出解V(k)；(4)m k max(V (k)) ,m k 1max(V (k 1))；(5) (k)V(k)/m k；(6)如果m k m k 1eps,则显示特征值1和对应的特征向量x(1) ),终止；7)k k 1, 转(3) 注：如上算法中解方程AV(k) (k 1)可以使用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和迭代允许最大次数n max ，程序即可给出每次迭代的次数和对应的迭代特征值、特征向量及误差序列，它们都按10 位有效数输出。
其中最后输出的结果即为所求的特征值和特征向量序。
如果迭代超出nmax 次还没有求出满足精度的根则输出迭代超限提示，此时可以根据输出序列判别收敛情况。
程序中变量说明a：存放矩阵Au：初始向量(0)和迭代过程中的向量(k)及所求特征向量v: 存放迭代过程中的向量V(k )a1：存放逆矩阵A 1m1: 存放所求特征值和迭代过程中的近似特征值nmax:存放迭代允许的最大次数eps: 存放误差精度fmax[x]给出向量x 中绝对值最大k: 记录迭代次数t1: 临时变量注：迭代最大次数可以修改为其他数字3、例题与实验2 2 3例3. 用反幂法求矩阵A 1 1 11 3 1的按模最小的特征值及其相应特征向量，要求误差eps 10 解：执行幂法程序后在输入的四个窗口中按提示分别输入:{{2.,-2.,3.},{1,1.,1},{1.,3,-1}} 、{1,0,1} 、0.00001 、100每次输入后用鼠标点击窗口的“ OK”按扭，得如下输出结果：注意到本题按模最小的特征值为1，因此求解效果较满意。
3.J acobi 方法1、Jacobi 旋转法算法1、输入矩阵A ，误差；2、For k 1,2, ；(2.1 )选取a(p k q1)max ai(j k 1)记录p, q2a(p k q1)(2.2 )由tan 2 (k21)a pq(k1), 确定旋转角，获得旋转矩阵J p,q,a pp a qq 42.3 ) a (p k j) pj a (p k j1) pj cos a q(k j1)sin , a(jp k) aqj jp (k) pjj p,q；2.4 )a(k) a qj a(k a pj 1)sin(k 1) a qj c os ,a(jp k)a(k) a qj jp,q；2.5 )a i(j k) ij a i(j k 1)iji, j p, q2.6 )a(k) a pp a(k 1)a ppcos2a q(k q1)sin22a(p kq 1)sin cos2.7 ) (k) a qq (k 1)a pp2 sin a(qq k 1)cos2(k2a pq 1)sin cos2.8 ) 计算E(A k ) ( k) 2 a ij2.9 )如果E(A k ) ，输出对角矩阵D diag( 1, 2注：如上算法中A k (a i(j k))，A0 (a i(j0)) (a ij )n )和特征向量矩阵J ，停止ij2、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 的所有特征值及其相应特征向量。