在自然科学和工程设计中的许多问题

第五章矩阵的特征值与特征向量的MATLAB程序59．高等教育出版社教育电子音像出版社作者：任玉杰在自然科学和工程设计中的许多问题，如电磁振荡、桥梁振动、机械振动等，常归结为求矩阵的特征值和特征向量.求矩阵的特征值和特征向量的问题是代数计算中的重要课题.本章着重介绍直接计算矩阵的特征值和特征向量的MATLAB程序、间接计算矩阵的特征值和特征向量的幂法、反幂法、雅可比方法、豪斯霍尔德方法和QR方法及其它们的MATLAB计算程序.最后我们还讨论广义特征值问题.5.1直接计算特征值和特征向量的MATLAB程序5.1.4计算特征值和特征向量的MATLAB程序从以上的讨论可以看到，有许多问题归结为求矩阵的特征值和特征向量，而用手工计算高阶矩阵的特征值与特征向量的难度较大，但是，计算机软件MATLAB提供了直接计算特征值与特征向量的MATLAB函数（见表5–1），下面介绍这些函数的使用方法.表5–1命令功能b=eig(A)输入方阵A，运行后输出b为由方阵A的全部特征值构成的列向量[V,D]=eig(A)输入对称矩阵A，运行后输出D为由A的全部特征值构成的对角矩阵，V的各列为对应于特征值的特征向量构成的矩阵，使得AV=DV[V,D]=eig(A,'nobalance')输入方阵A，运行后输出D为由A的全部特征值构成的对角矩阵，V的各列为对应于特征值的特征向量构成的矩阵，使得AV=DV；如果A是对称矩阵，则输出的结果与程序[V,D]=eig(A)的运行结果相同5.2幂法及其MATLAB程序幂法是求实矩阵A的主特征值（即实矩阵A按模最大的特征值）及其对应的特征向量的一种迭代方法.5.2.2幂法的MATLAB程序设n阶实矩阵A的n个特征值为n,,,21，且满足021n，A的主特征值1对应的特征向量为1X，则我们可以用下面的MATLAB程序计算1和1X的近似值和近似向量.用幂法计算矩阵A的主特征值和对应的特征向量的MATLAB主程序输入的量：n阶实矩阵A、n维初始实向量V0、计算要求的精度jd、迭代的最大次数max1；输出的量：迭代的次数k、A的主特征值1的近似值lambda、1对应的特征向量1X的近似向量Vk、相邻两次迭代的误差Wc.如果迭代次数已经达到最大的迭代次数max1，则给出提示的相关信息.第五章矩阵的特征值与特征向量的计算第五章矩阵的特征值与特征向量的MATLAB程序60．高等教育出版社教育电子音像出版社作者：任玉杰根据迭代公式（5.20），现提供用幂法计算矩阵A的主特征值和对应的特征向量的MATLAB主程序如下：function[k,lambda,Vk,Wc]=mifa(A,V0,jd,max1)lambda=0;k=1;Wc=1;,jd=jd*0.1;state=1;V=V0;while((k=max1)&(state==1))Vk=A*V;[mj]=max(abs(Vk));mk=m;tzw=abs(lambda-mk);Vk=(1/mk)*Vk;Txw=norm(V-Vk);Wc=max(Txw,tzw);V=Vk;lambda=mk;state=0;if(Wcjd)state=1;endk=k+1;Wc=Wc;endif(Wc=jd)disp('请注意：迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：')elsedisp('请注意：迭代次数k已经达到最大迭代次数max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相邻两次迭代的误差Wc如下：')endVk=V;k=k-1;Wc;例5.2.2用幂法计算下列矩阵的主特征值和对应的特征向量的近似向量，精度510.并把（1）和（2）输出的结果与例5.1.1中的结果进行比较.(1)4211A；(2)633312321B；(3)1124111221C；（4）20101350144D.解（1）输入MATLAB程序A=[1-1;24];V0=[1,1]';[k,lambda,Vk,Wc]=mifa(A,V0,0.00001,100),[V,D]=eig(A),Dzd=max(diag(D)),wuD=abs(Dzd-lambda),wuV=V(:,2)./Vk,运行后屏幕显示结果请注意：迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：k=lambda=Wc=333.000001738368048.691862856124999e-007Vk=V=wuV=-0.49999942054432-0.707106781186550.44721359549996-0.894428227562941.000000000000000.70710678118655-0.89442719099992-0.89442719099992Dzd=wuD=31.738368038406435e-006由输出结果可看出，迭代33次，相邻两次迭代的误差Wc8.6919e-007，矩阵A的主特征值的近似值lambda3.00000和对应的特征向量的近似向量Vk（-0.50000，1.00000T),lambda与例5.1.1中A的最大特征值32近似相等，绝对误差约为1.73837e-006，Vk与特征向量XT22kT)1,21()0(2k的第1个分量的绝对误差约等于0，第2个分量的绝对值相同.由wuV可以看出，2的特征向量V(:,2)与Vk的对应分量的比值近似相等.因此，用程序mifa.m计算的结果达到预先给定的精度510.(2)输入MATLAB程序B=[123;213;336];V0=[1,1,1]';[k,lambda,Vk,Wc]=mifa(B,V0,0.00001,100),[V,D]=eig(B),Dzd=max(diag(D)),wuD=abs(Dzd-lambda),wuV=V(:,3)./Vk,运行后屏幕显示结果请注意：迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的第五章矩阵的特征值与特征向量的MATLAB程序61．高等教育出版社教育电子音像出版社作者：任玉杰误差Wc如下：k=lambda=Wc=Dzd=wuD=39090Vk=wuV=0.500000000000000.816496580927730.500000000000000.816496580927731.000000000000000.81649658092773V=0.707106781186550.577350269189630.40824829046386-0.707106781186550.577350269189630.408248290463860-0.577350269189630.81649658092773由输出结果可看出，迭代3次，相邻两次迭代的误差Wc=0，实对称矩阵B的主特征值的近似值lambda=9和对应的特征向量的近似向量Vk=（0.50000，0.50000，1.00000T),lambda与例5.1.1中B的最大特征值93相同，Vk与特征向量XT33kT)2,1,1()0(3k的对应分量成比例.从wuV的每个分量的值也可以看出，3的特征向量V(:,3)与Vk的对应分量的比值相等.因此，用程序mifa.m计算的结果达到预先给定的精度510.此例说明，幂法对实对称矩阵B的迭代速度快且计算结果精度高，(3)输入MATLAB程序C=[122;1-11;4-121];V0=[1,1,1]';[k,lambda,Vk,Wc]=mifa(C,V0,0.00001,100),[V,D]=eig(C),Dzd=max(diag(D)),wuD=abs(Dzd-lambda),Vzd=V(:,1),wuV=V(:,1)./Vk,运行后屏幕显示请注意：迭代次数k已经达到最大迭代次数max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相邻两次迭代的误差Wc如下：k=lambda=Wc=1000.090909090909102.37758124193119Dzd=wuD=1.000000000000010.90909090909091Vk=Vzd=wuV=0.999999999999930.904534033733290.904534033733350.999999999999950.301511344577760.301511344577781.00000000000000-0.30151134457776-0.30151134457776由输出结果可见，迭代次数k已经达到最大迭代次数max1=100，并且lambda的相邻两次迭代的误差Wc2.377582,由wuV可以看出，lambda的特征向量Vk与真值Dzd的特征向量Vzd对应分量的比值相差较大，所以迭代序列发散.实际上，实数矩阵C的特征值的近似值为i,i,010000000001.000321　，并且对应的特征向量的近似向量分别为XT1=1k（0.90453403373329，0.30151134457776，-0.30151134457776）T，XT22k（-0.72547625011001,-0.21764287503300-0.07254762501100i,0.58038100008801-0.29019050004400i）T，XT33k(-0.72547625011001,-0.21764287503300+0.07254762501100i,0.58038100008801+0.29019050004400i)T0,0(21kk,03k是常数）.此例说明，当n阶实矩阵有复数特征值时，不宜用幂法计算它的主特征值1对应的特征向量1X.（4）输入MATLAB程序D=[-4140;-5130;-102];V0=[1,1,1]';[k,lambda,Vk,Wc]=mifa(D,V0,0.00001,100),[V,Dt]=eig(D),Dtzd=max(diag(Dt)),wuDt=abs(Dtzd-lambda),Vzd=V(:,2),wuV=V(:,2)./Vk,运行后屏幕显示结果请注意：迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：k=lambda=Wc=196.000006539495286.539523793591684e-006第五章矩阵的特征值与特征向量的MATLAB程序62．高等教育出版社教育电子音像出版社作者：任玉杰Dtzd=wuDt=6.000000000000006.539495284840768e-006Vk=Vzd=wuV=0.797400480535640.797400480535640.797400480535640.714285947838860.569571771811170.79740021980618-0.24999918247180-0.199350120133910.79740308813370由输出结果可见，迭代19次，相邻两次迭代的误差Wc6.53952e-006，矩阵D的主特征值的近似值lambda6.00001和对应的特征向量的近似向量为Vk（0.79740，0.71429，-0.25000T).用eig(A)计算，lambda与对应的特征值的真值Dtzd的绝对误差wuDt510，二者的特征向量V(:,2)与Vk的对应分量的比值近似相等（请对比wuV的每个分量的值）.因此，用程序mifa.m计算的结果达到预先给定的精度510.由例5.2.2的计算结果可见，用幂法计算实对称矩阵B的主特征值对应的特征向量时，得到的迭代序列的收敛速度最快且计算结果精度也最高；非实对称矩阵对应的