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、迭代初值向量u(0)、精度控制eps和迭代允许最大次数nmax,程序即可给出每次迭代的次数和对应的迭代特征值、特征向量及误差序列,它们都按10位有效数输出。
其中最后输出的结果即为所求的特征值和特征向量序。
如果迭代超出nmax次还没有求出满足精度的根则输出迭代超限提示,此时可以根据输出序列判别收敛情况。
程序中变量说明
a:
存放矩阵A
u:
初始向量u(0)和迭代过程中的向量u(k)及所求特征向量
v:
存放迭代过程中的向量V(k)
m1:
存放所求特征值和迭代过程中的近似特征值
nmax:
存放迭代允许的最大次数
eps:
存放误差精度
fmax[x]:
给出向量x中绝对值最大的分量
k:
记录迭代次数
t1:
临时变量
注:
迭代最大次数可以修改为其他数字。
例题与实验
例1.
用幂法求矩阵
的按模最大的特征值及其相应特征向量,要求误差<10-4。
解:
执行幂法程序后在输入的4个窗口中按提示分别输入
{{133,6,135},{44,5,46},{-88,-6,-90}},{1,1,1},0.0001,20
每次输入后用鼠标点击窗口的“OK”按扭,得如下输出结果。
k=1特征值=44.42335766误差=229.5766423
特征向量={1.,0.3467153285,-0.6715328467}
k=2特征值=44.92343082误差=0.5000731606
特征向量={1.,0.3341275058,-0.6672691423}
k=3特征值=44.99546459误差=0.07203376236
特征向量={1.,0.3333729572,-0.6667020234}
k=4特征值=44.99977337误差=0.004308781874
特征向量={1.,0.3333351894,-0.6666684279}
k=5特征值=44.99998937误差=0.0002160020115
特征向量={1.,0.3333334179,-0.6666667492}
k=6特征值=44.99999952误差=0.0000101441501
特征向量={1.,0.3333333371,-0.6666666704}
此结果说明迭代6次,求得误差为err=0.0000101441501的按模最大的特征值=44.99999952及其对应的一个特征向量={1.,0.3333333371,-0.6666666704}。
本题矩阵A的3个特征值为{45.,2.,1.},可见所求结果很好。
但如果执行幂法程序后在输入的4个窗口中按提示分别输入
{{133,6,135},{44,5,46},{-88,-6,-90}},{1,1,-1},0.0001,20
每次输入后用鼠标点击窗口的“OK”按扭,得如下输出结果。
k=1特征值=2.5误差=1.5
特征向量={1.,0.75,-1.}
k=2特征值=2.2误差=0.3
特征向量={1.,0.7,-1.}
k=3特征值=2.090909091误差=0.1090909091
特征向量={1.,0.6818181818,-1.}
k=4特征值=2.043478261误差.0474*******
特征向量={1.,0.6739130435,-1.}
k=5特征值=2.021276596误差.022********
特征向量={1.,0.670212766,-1.}
k=6特征值=2.010526316误差.010********
特征向量={1.,0.6684210526,-1.}
k=7特征值=2.005235602误差=0.005290713695
特征向量={1.,0.667539267,-1.}
k=8特征值=2.002610966误差=0.002624636037
特征向量={1.,0.6671018277,-1.}
k=9特征值=2.001303781误差=0.001307185093
特征向量={1.,0.6668839635,-1.}
k=10特征值=2.000651466误差=0.0006523151668
特征向量={1.,0.6667752443,-1.}
k=11特征值=2.000325627误差=0.0003258389664
特征向量={1.,0.6667209378,-1.}
k=12特征值=2.000162787误差=0.0001628399197
特征向量={1.,0.6666937978,-1.}
k=13特征值=2.000081387误差=0.00008140008032
特征向量={1.,0.6666802311,-1.}
此结果说明迭代13次,求得误差为err=0.00008140008032的按模最大的特征值=2.000081387及其对应的一个特征向量={1.,0.6666802311,-1.}。
选用不同的迭代初值获得两个不同结果,显然第二个特征值=2.000081387不是模最大的特征值。
上面实验说明使用幂法依赖于迭代初值的选取且有时得到的结果不是模最大的特征值(知道是什么原因吗?
)。
不过一般情况下,幂法是可以求出按模最大的特征值的。
如果不放心,可以选用两个不同的初值迭代计算,通过计算结果可以马上确定按模最大的特征值。
例2.用幂法求矩阵
的按模最大的特征值及其相应特征向量,要求误差<10-5。
解:
执行幂法程序后在输入的4个窗口中按提示分别输入
{{2.,-2.,3.},{1,1.,1},{1.,3,-1}},{1,0,1},0.00001,20
每次输入后用鼠标点击窗口的“OK”按扭,得如下输出结果。
…………………………………………………….
k=9特征值=2.990381958误差.0454*******
特征向量={1.,0.9810402963,0.952738931}
k=10特征值=2.974823934误差=0.01555802398
特征向量={0.9684837058,0.9810717388,1.}
k=11特征值=2.99573741误差.020********
特征向量={1.,0.9915058875,0.9787802529}
k=12特征值=2.988679144误差=0.007058265982
特征向量={0.985843744,0.9915041722,1.}
k=13特征值=2.998102573误差=0.009423429932
特征向量={1.,0.996208617,0.9905232775}
k=14特征值=2.994943938误差=0.003158635286
特征向量={0.993679344,0.9962073749,1.}
k=15特征值=2.999155514误差=0.004211575789
特征向量={1.,0.9983114144,0.9957787292}
k=16特征值=2.997748176误差=0.001407337687
特征向量={0.9971851559,0.9983110678,1.}
k=17特征值=2.999624371误差=0.001876194866
特征向量={1.,0.9992487853,0.9981219847}
k=18特征值=2.998998284误差=0.0006260875486
特征向量={0.9987478473,0.9992487055,1.}
k=19特征值=2.999832987误差=0.0008347031267
特征向量={1.,0.9996659782,0.9991649479}
k=20特征值=2.999554616误差=0.0002783707007
特征向量={0.9994432692,0.9996659612,1.}
迭代超限
此结果说明迭代20次后还没有得到满足要求的解,但观察特征值序列发现其是收敛的,因此可以增大迭代次数以求得满足要求的解。
本题将最大迭代次数设定为100后得出在迭代第30次时的满足要求的解为
k=30特征值=2.999992274误差=4.82876832510-6
特征向量={0.9999903425,0.9999942055,1.}
注意到本题按模最大的特征值为3,因此求解效果较满意。
反幂法
反幂法规范化算法
1.输入矩阵A、初始向量u(0),误差eps
2.k1
3.解方程AV(k)=u(k-1)求出解V(k)
4.mkmax(V(k)),mk-1max(V(k-1))
5.ukV(k)/mk
6.如果|mk-mk-1|(1),终止
7.kk+1,转3
注:
如上算法中解方程AV(k)=u(k-1)可以使用Dololittle分解法。
本算法使用了数据规范化处理技术以防止计算过程中出现益出错误。
规范化反幂法程序
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&&ku=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、迭代初值向量u(0)、精度控制eps和迭代允许最大次数nmax,程序即可给出每次迭代的次数和对应的迭代特征值、特征向量及误差序列,它们都按10位有效数输出。
其中最后输出的结果即为所求的特征值和特征向量序。
如果迭代超出nmax次还没有求出满足精度的根则输出迭代超限提示,此时可以根据输出序列判别收敛情况。
程序中变量说明:
a:
存放矩阵A
u初始向量u(0)和迭代过程中的向量u(k)及所求特征向量
v:
存放迭代过程中的向量V(k)
a1:
存放逆矩阵A-1
m1:
存放所求特征值和迭代过程中的近似特征值
nmax:
存放迭代允许的最大次数
eps:
存放误差精度
fmax[x]:
给出向量x中绝对值最大的分量
k:
记录迭代次数
t1:
临时变量
注:
迭代最大次数可以修改为其他数字。
例题与实验
例3.用反幂法求矩阵
的按模最小的特征值及其相应特征向量,要求误差<10-5。
解:
执行幂法程序后在输入的4个窗口中按提示分别输入
{{2.,-2.,3.},{1,1.,1},{1.,3,-1}},{1,0,1},0.00001,100
每次输入后用鼠标点击窗口的“OK”按扭,得如下输出结果。
…………………………………………………….
k=12特征值=0.9995193167误差=0.002553354512
特征向量={-0.9990991337,1.,0.9988281951}
k=13特征值=1.00052327误差=0.001003953492
特征向量={-0.9998950625,0.9994143792,1.}
k=14特征值=0.9998814738误差=0.0006417963548
特征向量={-0.9997696935,1.,0.9997070364}
k=15特征值=1.000131384误差=0.0002499099059
特征向量={-0.9999720617,0.9998535356,1.}
k=16特征值=0.9999705536误差=0.0001608301282
特征向量={-0.9999418581,1.,0.9999267582}
k=17特征值=1.000032909误差=0.00006235517185
特征向量={-0.9999928266,0.9999633802,1.}
k=18特征值=0.9999926591误差=0.00004024964909
特征向量={-0.9999854018,1.,0.9999816895}
k=19特征值=1.000008234误差=0.00001557505036
特征向量={-0.9999981857,0.9999908448,1.}
k=20特征值=0.9999981671误差=0.00001006707749
特征向量={-0.9999963435,1.,0.9999954224}
k=21特征值=1.000002059误差=3.89222610510-6
特征向量={-0.9999995441,0.9999977112,1.}
此结果说明迭代21次后得到满足要求的解:
k=21特征值=1.000002059误差=3.89222610510-6
特征向量={-0.9999995441,0.9999977112,1.}
注意到本题按模最小的特征值为1,因此求解效果较满意。
如果将如上输入的迭代初值改为常用的{1,1,1},则得到如下结果。
k=1特征值=3.误差=2.666666667
特征向量={1.,1.,1.}
k=2特征值=3.误差=2.22044604910-15
特征向量={1.,1.,1.}
这个结果是按模最大的特征值,而不是按模最小的特征值,因此选用初值要小心。
Jacobi方法
Jacobi旋转法算法
注:
如上算法中Ak=(aij(k)),A0=(aij(0))=(aij)。
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&&km=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]]]
];
If[k>=nmax,Print[“迭代超限”]]
说明:
本程序用于求对称矩阵A的所有特征值及其相应特征向量。
程序执行后,先通过键盘输入矩阵A、矩阵阶数n、精度控制eps和迭代允许最大次数nmax,程序即可给出每次迭代的次数和对应的迭代特征值、特征向量及误差序列。
其中最后输出的结果即为所求的特征值和特征向量。
如果迭代超出nmax次还没有求出满足精度的根则输出迭代超限提示。
此外,输出的特征值矩阵可以不是真正的对角矩阵,但它们的主对角元素就是满足要求的所有特征值。
程序中变量说明
a:
存放矩阵A及其相似变换过程中的Ak
bb:
存放特征向量矩阵J的转置
v:
存放迭代过程中的向量V(k)
m1:
存放所求特征值和迭代过程中的近似特征值
nmax:
存放迭代允许的最大次数
eps:
存放误差精度
k:
记录迭代次数
t1,mu,s,c,ea,p,q,m:
临时变量
a1,qq,pp:
临时向量
例题与实验
例4.用Jacobi方法求矩阵
的所有特征值及其相应特征向量,要求误差<10-6。
解:
执行Jacobi方法程序后在输入的4个窗口中按提示分别输入
{{2,-1,0},{-1,2,-1},{0,-1,2}},3,0.000001,30
每次输入后用鼠标点击窗口的“OK”按扭,得如下输出结。
迭代次数k=0
误差=2.
相似矩阵A=
1.0-0.707107
03.-0.707107
-0.707107-0.7071072
特征向量J
0.707107-0.7071070
0.7071070.7071070
001
迭代次数k=1
误差=1.
相似矩阵A=
0.633975-0.3250580
-0.3250583.-0.627963
0-0.6279632.36603
特征向量J
0.627963-0.707107-0.325058
0.6279630.707107-0.325058
0.45970100.888074
迭代次数k=2
误差=0.211325
相似矩阵A=
0.633975-0.276837-0.170364
-0.2768373.386450
-0.17036401.97958
特征向量J
0.627963-0.431846-0.647434
0.6279630.7725740.0937614
0.459701-0.4654440.756332
…………..
迭代次数k=5
误差=8.3089110-6
相似矩阵A=
0.5857880.00203811-0.0000241542
0.002038113.414210
-0.000024154202.
特征向量J
0.499652-0.50036-0.707098
0.7076160.7065970.0000120684
0.499628-0.500360.707115
迭代次数k=6
误差=1.1668510-9
相似矩阵A=
0.5857860-0.0000241542
03.41421-1.740510-8
-0.0000241542-1.740510-82.
特征向量J
0.500012-0.5-0.707098
0.7071070.7071070.0000120684
0.499988-0.50.707115
此结果说明迭代6次,求得误差为1.1668510-9的矩阵A的所有特征值1=0.585786,2=3.41421,3=2及对应的特征向量
{0.500012,0.707107,0.499988},
{-0.5,0.707107,-0.5},
{-0.707098,0.0000120684,0.707115}
本题矩阵A的3个特征值为{2-sqrt[2],2+sqrt[2],2}。
线性代数告诉我们,不同特征值对应的特征向量是正交的,用如上求出的特征向量验证,可以得到
{0.500012,0.707107,0.499988}.{-0.5,0.707107,-0.5}=-5.5511210-17,
{0.500012,0.707107,0.499988}.{-0.707098,0.0000120684,0.707115}=-5.5511210-17
{-0.5,0.707107,-0.5