求矩阵特征值算法及程序简介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 V a1：存放逆矩阵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 ji k pqa a 记录q p , （2.2）由4,22tan )1()1()1(πθθ≤-=---k qqk pp k pq a 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 qjk jp k qj k pj k 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 ppk qq a a a a（2.8）计算()∑≠=ji k ijk 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 的所有特征值及其相应特征向量。