求矩阵特征值算法及程序简介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 的所有特征值及其相应特征向量。