求矩阵特征值算法及程序简介

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
此结果说明迭代 6 次，求得误差为 0.0000101442 的按模最大的特征值为 44.99999952, 及其对应的一个特征向量:{1.000000000,0.3333333371,-0.6666666704}
2.反幂法
1、 反幂法规范化算法
（1）输入矩阵 A 、初始向量 （2） k 1 ； （3）计算 AV
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”按扭，得如下输出结果：
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;
说明 本程序用于求对称矩阵 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、例题与实验
(k ) ( 0) 注：如上算法中 Ak (aij ) ， A0 (aij ) (aij ) 。
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[mu0,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]]]];
注意到本题按模最小的特征值为 1，因此求解效果较满意。
3.Jacobi 方法
1、Jacobi 旋转法算法
1、输入矩阵 A ，误差 ； 2、 For
k 1,2, ；
( k 1) ( k 1) 记录 p, q max aij i j
（2.1）选取 a pq
（2.2）由 tan 2
,终止；
注：如上算法中解方程 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]," Print[" 特征向量=",N[u,10]]]; If[knmax,Print["迭代超限"]] 误差=",N[t1,10]];
（2.6） a pp a pp
(k )
( k 1)
（2.7） aqq a pp
(k )
( k 1)
（2.8）计算 E ( Ak )
a
i j
(k ) 2 ij
（2.9）如果 E ( Ak ) ，输出对角矩阵 D diag (1 , 2 ,, n ) 和特征向量矩阵 J ，停止
2 2 3 例 3. 用反幂法求矩阵 A 1 1 1 1 3 1
的按模最小的特征值及其相应特征向量，要求误差 eps 10 5 。 解：执行幂法程序后在输入的四个窗口中按提示分别输入: {{2.,-2.,3.},{1,1.,1},{1.,3,-1}}、{1,0,1}、0.00001、100 每次输入后用鼠标点击窗口的“OK”按扭，得如下输出结果：
求矩阵特征值算法及程序简介
1.幂法
1、幂法规范化算法 （1）输入矩阵 A 、初始向量 （2） k 1 ； （3）计算 V
(k )
(0)
，误差 eps ；
A ( k 1) ；
(k )
（4） mk max( V （5）
(k )
)
, mk 1 max( V ( k 1) ) ；
V ( k ) / mk ；
(1) )
（6）如果 mk mk 1 eps ,则显示特征值 1 和对应的特征向量 x （7） k k 1,转（3）
,终止；
注：如上算法中的符号 max( V ) 表示取向量 V 中绝对值最大的分量。本算法使用了数据 规范化处理技术以防止计算过程中出现益出错误。
说明：本程序用于求矩阵 A 按模最小的特征值及其相应特征向量。程序执行后，先通过 键盘输入矩阵 A 、迭代初值向量 ( 0 ) 、精度控制 eps 和迭代允许最大次数 n max ，程序 即可给出每次迭代的次数和对应的迭代特征值、特征向量及误差序列，它们都按 10 位 有效数输出。其中最后输出的结果即为所求的特征值和特征向量序。如果迭代超出 n max 次还没有求出满足精度的根则输出迭代超限提示，此时可以根据输出序列判别收 敛情况。 程序中变量说明 a：存放矩阵 A u：初始向量 ( 0 ) 和迭代过程中的向量 ( k ) 及所求特征向量
j p, q ； j p, q ；
（2.4） aqj a pj （2.5） aij
(k )
( k 1)
( k 1) aij
i, j p, q
( k 1) k 1) cos 2 aqq sin 2 2a (pq sin cos ( k 1) k 1) sin 2 aqq cos 2 2a (pq sin cos
m0=m1; m1=fmax[v]; t=Abs[m1-m0]//N; Print["k=",k," 特征值=",N[m1,10]," Print[" 特征向量=",N[u,10]]]; If[knmax,Print["迭代超限"]] 误差=",N[t,10]];
说明：本程序用于求矩阵 A 按模最大的特征值及其相应特征向量。程序执行后，先通过 键盘输入矩阵 A 、迭代初值向量 ( 0 ) 、精度控制 eps 和迭代允许最大次数 n max ，程序 即可给出每次迭代的次数和对应的迭代特征值、特征向量及误差序列，它们都按 10 位 有效数输出。其中最后输出的结果即为所求的特征值和特征向量序列。如果迭代超出 n max 次还没有求出满足精度的根则输出迭代超限提示，此时可以根据输出序列判别收 敛情况。 程序中变量说明 a:存放矩阵 A ； u:初始向量 ( 0 ) 和迭代过程中的向量 ( k ) 及所求特征向量； v:存放迭代过程中的向量 V ( k ) ； m1:存放所求特征值和迭代过程中的近似特征值； nmax:存放迭代允许的最大次数； eps:存放误差精度； fmax[x]: 给出向量 x 中绝对值最大的分量； k:记录迭代次数； t1:临时变量； 注：迭代最大次数可以修改为其他数字。 3、例题与实验
v:存放迭代过程中的向量 V ( k ) a1：存放逆矩阵 A 1 m1: 存放所求特征值和迭代过程中的近似特征值 nmax:存放迭代允许的最大次数 eps:存放误差精度 fmax[x]: 给出向量 x 中绝对值最大的分量 k:记录迭代次数 t1:临时变量 注：迭代最大次数可以修改为其他数字。
3、例题与实验
(k )
(0)
，误差 eps ；
( k 1) 求出解 V ( k ) ；
(k )
（4） mk max( V （5）
(k )
)
, mk 1 max( V ( k 1) ) ；
V ( k ) / mk ；
(1) )
（6）如果 mk mk 1 eps ,则显示特征值 1 和对应的特征向量 x （7） k k 1,转（3）