第五章矩陣的特征值與特征向量的計算第五章矩陣的特征值與特征向量的計算在自然科學和工程設計中的許多問題，如電磁振蕩、橋梁振動、機械振動等，常歸結(jié)為求矩陣的特征值和特征向量.求矩陣的特征值和特征向量的問題是代數(shù)計算中的重要課題.本章著重介紹直接計算矩陣的特征值和特征向量的MATLAB程序、間接計算矩陣的特征值和特征向量的冪法、反冪法、雅可比方法、豪斯霍爾德方法和QR方法及其它們的MATLAB計算程序.最后我們還討論廣義特征值問題.5.1直接計算特征值和特征向量的MATLAB程序5.1.4計算特征值和特征向量的MATLAB程序從以上的討論可以看到，有許多問題歸結(jié)為求矩陣的特征值和特征向量，而用手工計算高階矩陣的特征值與特征向量的難度較大，但是，計算機軟件MATLAB提供了直接計算特征值與特征向量的MATLAB函數(shù)（見表5–1），下面介紹這些函數(shù)的使用方法.表5–1命令功能b=eig(A)輸入方陣A，運行后輸出b為由方陣A的全部特征值構(gòu)成的列向量[V,D]=eig(A)輸入對稱矩陣A，運行后輸出D為由A的全部特征值構(gòu)成的對角矩陣，V的各列為對應于特征值的特征向量構(gòu)成的矩陣，使得AV=DV[V,D]=eig(A,'nobalance')輸入方陣A，運行后輸出D為由A的全部特征值構(gòu)成的對角矩陣，V的各列為對應于特征值的特征向量構(gòu)成的矩陣，使得AV=DV；如果A是對稱矩陣，則輸出的結(jié)果與程序[V,D]=eig(A)的運行結(jié)果相同5.2冪法及其MATLAB程序冪法是求實矩陣的主特征值（即實矩陣按模最大的特征值）及其對應的特征向量的一種迭代方法.5.2.2冪法的MATLAB程序設階實矩陣的個特征值為，且滿足，的主特征值對應的特征向量為，則我們可以用下面的MATLAB程序計算和的近似值和近似向量.用冪法計算矩陣的主特征值和對應的特征向量的MATLAB主程序輸入的量：階實矩陣、維初始實向量V0、計算要求的精度jd、迭代的最大次數(shù)max1；輸出的量：迭代的次數(shù)k、的主特征值的近似值lambda、對應的特征向量的近似向量Vk、相鄰兩次迭代的誤差Wc.如果迭代次數(shù)已經(jīng)達到最大的迭代次數(shù)max1，則給出提示的相關(guān)信息.根據(jù)迭代公式（5.20），現(xiàn)提供用冪法計算矩陣的主特征值和對應的特征向量的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(Wc>jd)state=1;endk=k+1;Wc=Wc;endif(Wc<=jd)disp('請注意：迭代次數(shù)k,主特征值的近似值lambda,主特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：')elsedisp('請注意：迭代次數(shù)k已經(jīng)達到最大迭代次數(shù)max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相鄰兩次迭代的誤差Wc如下：')endVk=V;k=k-1;Wc;例5.2.2用冪法計算下列矩陣的主特征值和對應的特征向量的近似向量，精度.并把（1）和（2）輸出的結(jié)果與例5.1.1中的結(jié)果進行比較.(1)；(2)；(3)；（4）.解（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,運行后屏幕顯示結(jié)果請注意：迭代次數(shù)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由輸出結(jié)果可看出，迭代33次，相鄰兩次迭代的誤差Wc8.6919e-007，矩陣的主特征值的近似值lambda3.00000和對應的特征向量的近似向量Vk（-0.50000，1.00000,lambda與例5.1.1中的最大特征值近似相等，絕對誤差約為1.73837e-006，Vk與特征向量的第1個分量的絕對誤差約等于0，第2個分量的絕對值相同.由wuV可以看出，的特征向量V(:,2)與Vk的對應分量的比值近似相等.因此，用程序mifa.m計算的結(jié)果達到預先給定的精度.(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,運行后屏幕顯示結(jié)果請注意：迭代次數(shù)k,主特征值的近似值lambda,主特征向量的近似向量Vk,相鄰兩次迭代的誤差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由輸出結(jié)果可看出，迭代3次，相鄰兩次迭代的誤差Wc=0，實對稱矩陣B的主特征值的近似值lambda=9和對應的特征向量的近似向量Vk=（0.50000，0.50000，1.00000,lambda與例5.1.1中的最大特征值相同，Vk與特征向量的對應分量成比例.從wuV的每個分量的值也可以看出，的特征向量V(:,3)與Vk的對應分量的比值相等.因此，用程序mifa.m計算的結(jié)果達到預先給定的精度.此例說明，冪法對實對稱矩陣的迭代速度快且計算結(jié)果精度高，(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,運行后屏幕顯示請注意：迭代次數(shù)k已經(jīng)達到最大迭代次數(shù)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由輸出結(jié)果可見，迭代次數(shù)k已經(jīng)達到最大迭代次數(shù)max1=100，并且lambda的相鄰兩次迭代的誤差Wc2.37758>2,由wuV可以看出，lambda的特征向量Vk與真值Dzd的特征向量Vzd對應分量的比值相差較大，所以迭代序列發(fā)散.實際上，實數(shù)矩陣C的特征值的近似值為，并且對應的特征向量的近似向量分別為=（0.90453403373329，0.30151134457776，-0.30151134457776），（-0.72547625011001,-0.21764287503300-0.07254762501100i,0.58038100008801-0.29019050004400i），(-0.72547625011001,-0.21764287503300+0.07254762501100i,0.58038100008801+0.29019050004400i),是常數(shù)）.此例說明，當階實矩陣有復數(shù)特征值時，不宜用冪法計算它的主特征值對應的特征向量.（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,運行后屏幕顯示結(jié)果請注意：迭代次數(shù)k,主特征值的近似值lambda,主特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：k=lambda=Wc=196.000006539495286.539523793591684e-006Dtzd=wuDt=6.000000000000006.539495284840768e-006Vk=Vzd=wuV=0.797400480535640.797400480535640.797400480535640.714285947838860.569571771811170.79740021980618-0.24999918247180-0.199350120133910.79740308813370由輸出結(jié)果可見，迭代19次，相鄰兩次迭代的誤差Wc6.53952e-006，矩陣D的主特征值的近似值lambda6.00001和對應的特征向量的近似向量為Vk（0.79740，0.71429，-0.25000.用eig(A)計算，lambda與對應的特征值的真值Dtzd的絕對誤差wuDt，二者的特征向量V(:,2)與Vk的對應分量的比值近似相等（請對比wuV的每個分量的值）.因此，用程序mifa.m計算的結(jié)果達到預先給定的精度.由例5.2.2的計算結(jié)果可見，用冪法計算實對稱矩陣的主特征值對應的特征向量時，得到的迭代序列的收斂速度最快且計算結(jié)果精度也最高；非實對稱矩陣對應的迭代序列的斂散性不定，有時發(fā)散（例如矩陣C）,有時收斂（例如矩陣A和D），且收斂速度較慢，比較（1）和（4）的計算結(jié)果可知，用冪法得到的迭代序列的收斂速度與矩陣的階數(shù)無關(guān)；當實矩陣有復數(shù)特征值時，不宜用冪法計算它的主特征值對應的特征向量.5.3反冪法和位移反冪法及其MATLAB程序反冪法是求非奇異矩陣的按模最小特征值及其對應的特征向量的一種迭代方法.5.3.3原點位移反冪法的MATLAB程序設階實矩陣的個特征值都不相同，且滿足，且的按模最小特征值對應的特征向量為，則我們可以根據(jù)迭代公式（5.27）和（5.28）編寫兩種MATLAB程序分別計算和的近似值和近似向量.（一）原點位移反冪法的MATLAB主程序1根據(jù)原點位移反冪法的迭代公式（5.28），現(xiàn)提供用原點位移反冪法計算矩陣的按模最小特征值和對應的特征向量的MATLAB主程序如下：用原點位移反冪法計算矩陣的特征值和對應的特征向量的MATLAB主程序1輸入的量：階實矩陣、維初始實向量V0、特征值的近似值jlamb、計算的精度jd、迭代的最大次數(shù)max1；輸出的量：迭代的次數(shù)k、的特征值的近似值lambdan、與對應的特征向量的近似向量Vk、相鄰兩次迭代的誤差Wc.如果迭代次數(shù)已經(jīng)達到最大的迭代次數(shù)max1，則給出提示的相關(guān)信息.function[k,lambdan,Vk,Wc]=ydwyfmf(A,V0,jlamb,jd,max1)[n,n]=size(A);A1=A-jlamb*eye(n);jd=jd*0.1;RA1=det(A1);ifRA1==0disp('請注意：因為A-aE的n階行列式hl等于零，所以A-aE不能進行Lu分解.')returnendlambda=0;ifRA1~=0forp=1:nh(p)=det(A1(1:p,1:p));endhl=h(1:n);fori=1:nifh(1,i)==0disp('請注意：因為A-aE的r階主子式等于零，所以A-aE不能進行Lu分解.')returnendendifh(1,i)~=0disp('請注意：因為A-aE的各階主子式都不等于零，所以A-aE能進行Lu分解.')k=1;Wc=1;state=1;Vk=V0;while((k<=max1)&(state==1))[LU]=lu(A1);Yk=L\Vk;Vk=U\Yk;[mj]=max(abs(Vk));mk=m;Vk1=Vk/mk;Yk1=L\Vk1;Vk1=U\Yk1;[mj]=max(abs(Vk1));mk1=m;Vk2=(1/mk1)*Vk1;tzw1=abs((mk-mk1)/mk1);tzw2=abs(mk1-mk);Txw1=norm(Vk)-norm(Vk1);Txw2=(norm(Vk)-norm(Vk1))/norm(Vk1);Txw=min(Txw1,Txw2);tzw=min(tzw1,tzw2);Vk=Vk2;mk=mk1;Wc=max(Txw,tzw);Vk=Vk2;mk=mk1;state=0;if(Wc>jd)state=1;endk=k+1;%Vk=Vk2,mk=mk1,endif(Wc<=jd)disp('A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：')elsedisp('A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k已經(jīng)達到最大迭代次數(shù)max1,按模最小特征值的迭代值lambda,特征向量的迭代向量Vk,相鄰兩次迭代的誤差Wc如下：')endhl,RA1endend[V,D]=eig(A,'nobalance'),Vk;k=k-1;Wc;lambdan=jlamb+1/mk1;例5.3.2用原點位移反冪法的迭代公式（5.28），根據(jù)給定的下列矩陣的特征值的初始值，計算與對應的特征向量的近似向量,精確到0.0001.（1）,；（2），；（3），.解（1）輸入MATLAB程序>>A=[1-10;-24-2;0-12];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,0.2,0.0001,10000)運行后屏幕顯示結(jié)果請注意：因為A-aE的各階主子式都不等于零，所以A-aE能進行Lu分解.A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：k=lambda=Wc=hl=30.23841.0213e-0070.80001.04000.2720Vk=V=D=1.0000-0.2424-1.0000-0.57075.1249000.76161.0000-0.76160.363300.238400.4323-0.3200-0.43231.0000001.6367（2）輸入MATLAB程序>>A=[1-1;24];V0=[20,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,2.001,0.0001,100)運行后屏幕顯示結(jié)果請注意：因為A-aE的各階主子式都不等于零，所以A-aE能進行Lu分解.A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：k=lambda=Wc=hl=22.00205.1528e-007-1.0010-0.0010Vk=V=D=1.0000-1.00000.500020-1.00001.0000-1.000003（3）輸入MATLAB程序>>A=[-11215;2583;153-3];V0=[1,1,-1]';[k,lambdan,Vk,Wc]=ydwyfmf(A,V0,8.26,0.0001,100)運行后屏幕顯示結(jié)果請注意：因為A-aE的各階主子式都不等于零，所以A-aE能進行Lu分解.A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：k=lambdan=Wc=hl=28.26406.9304e-008-19.2600-961.9924-6.1256Vk=V=D=-0.76920.79280.60810.0416-22.5249000.09120.0030-0.07210.997408.26400-1.0000-0.60950.79060.05900058.2609例5.3.3用原點位移反冪法的迭代公式（5.28），計算的分別對應于特征值，，的特征向量，，的近似向量，相鄰迭代誤差為0.001.將計算結(jié)果與精確特征向量比較，其中0.99999999999997…，1.99999999999999…，4.00000000000002…，分別對應特征向量為0.50000000000000，0.50000000000000，1.00000000000000，-0.24999999999999，-0.50000000000000，-1.00000000000000，-0.40000000000000，-0.60000000000000，-1.00000000000000.解（1）計算特征值對應的特征向量的近似向量.輸入MATLAB程序>>A=[011-5;-217-7;-426-10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,1.001,0.001,100),[V,D]=eig(A);Dzd=min(diag(D)),wuD=abs(Dzd-lambda),VD=V(:,1),wuV=V(:,1)./Vk,運行后屏幕顯示結(jié)果請注意：因為A-aE的各階主子式都不等于零，所以A-aE能進行Lu分解.A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：hl=-1.001000000000005.98500100000000-0.00299600100000k=lambda=RA1=51.00200000000000-0.00299600100000Vk=VD=wuV=-0.50000000000000-0.408248290463860.81649658092773-0.50000000000000-0.408248290463860.81649658092773-1.00000000000000-0.816496580927730.81649658092773Wc=Dzd=wuD=1.378794763695562e-0091.000000000000000.00200000000000從輸出的結(jié)果可見，迭代5次，特征向量的近似向量的相鄰兩次迭代的誤差Wc1.379e-009，由wuV可以看出，=Vk與VD的對應分量的比值相等.特征值的近似值lambda1.002與初始值1.001的絕對誤差為0.001，而與的絕對誤差為0.002，其中,.（2）計算特征值對應特征向量的近似向量.輸入MATLAB程序>>A=[011-5;-217-7;-426-10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,2.001,0.001,100),[V,D]=eig(A);WD=lambda-D(2,2),VD=V(:,2),wuV=V(:,2)./Vk,運行后屏幕顯示結(jié)果請注意：因為A-aE的各階主子式都不等于零，所以A-aE能進行Lu分解.A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：hl=-2.00100000000000-8.012999000000000.00200099900000k=Wc=lambda=WD=23.131363162302120e-0072.002000000000160.00200000000016Vk=VD=wuV=-0.249999999999990.21821789023599-0.87287156094401-0.499999999999990.43643578047198-0.87287156094398-1.000000000000000.87287156094397-0.87287156094397從輸出的結(jié)果可見，迭代2次，特征向量的近似向量的相鄰兩次迭代的誤差Wc3.131e-007，與的對應分量的比值近似相等.特征值的近似值lambda2.002與初始值2.001的絕對誤差約為0.001，而lambda與的絕對誤差約為0.002，其中,.（3）計算特征值對應特征向量的近似向量.輸入MATLAB程序>>A=[011-5;-217-7;-426-10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,4.001,0.001,100)[V,D]=eig(A);WD=lambda-max(diag(D)),VD=V(:,3),wuV=V(:,3)./Vk,運行后屏幕顯示結(jié)果請注意：因為A-aE的各階主子式都不等于零，所以A-aE能進行Lu分解.A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：hl=-4.00100000000000-30.00899900000000-0.00600500099999k=lambda=Wc=WD=24.001999999999901.996084182914842e-0070.00199999999990Vk=VD=wuV=0.40000000000001-0.32444284226153-0.811107105653800.60000000000001-0.48666426339229-0.811107105653811.00000000000000-0.81110710565381-0.81110710565381從輸出的結(jié)果可見，迭代2次，特征向量的近似向量的相鄰兩次迭代的誤差Wc1.996e-007，與的對應分量的比值近似相等.特征值的近似值的絕對誤差近似為，而lambda與的絕對誤差約為0.002，其中-0.40000000000000，-0.60000000000000，-1.00000000000000,.綜上所述，用原點位移反冪法的迭代公式（5.28）求矩陣的全部特征值分別對應特征向量的收斂速度快，且精確度高.但是求矩陣的全部特征值的結(jié)果不理想.（二）、原點位移反冪法的MATLAB主程序2根據(jù)迭代公式（5.27），我編寫用原點位移反冪法計算矩陣的按模最小特征值和對應的特征向量的MATLAB主程序如下：用原點位移反冪法計算矩陣的特征值和對應的特征向量的MATLAB主程序2輸入的量：階實矩陣、維初始實向量V0、特征值的近似值jlambn、計算的精度jd、迭代的最大次數(shù)max1.輸出的量：迭代的次數(shù)k、的特征值的近似值lambdan、對應的特征向量的近似向量Vk、相鄰兩次迭代的誤差Wc.如果迭代次數(shù)已經(jīng)達到最大的迭代次數(shù)max1，則給出提示的相關(guān)信息.function[k,lambdan,Vk,Wc]=wfmifa1(A,V0,jlamb,jd,max1)[n,n]=size(A);jd=jd*0.1;A1=A-jlamb*eye(n);nA1=inv(A1);lambda1=0;k=1;Wc=1;state=1;U=V0;while((k<=max1)&(state==1))Vk=A1\U;[mj]=max(abs(Vk));mk=m;Vk=(1/mk)*Vk;Vk1=A1\Vk;[m1j]=max(abs(Vk1));mk1=m1,Vk1=(1/mk1)*Vk1;U=Vk1,Txw=(norm(Vk1)-norm(Vk))/norm(Vk1);tzw=abs((lambda1-mk1)/mk1);Wc=max(Txw,tzw);lambda1=mk1;state=0;if(Wc>jd)state=1;endk=k+1;endif(Wc<=jd)disp('請注意迭代次數(shù)k,特征值的近似值lambda,對應的特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：')elsedisp('請注意迭代次數(shù)k已經(jīng)達到最大迭代次數(shù)max1,特征值的近似值lambda,對應的特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：')end[V,D]=eig(A,'nobalance'),Vk=U;k=k-1;Wc;lambdan=jlamb+1/mk;例5.3.4用原點位移反冪法的迭代公式（5.27），計算例題5.3.3，并且將這兩個例題的計算結(jié)果進行比較.再用兩種原點位移反冪法的MATLAB主程序，求對應的特征向量.解（1）計算特征值對應特征向量的近似向量.輸入MATLAB程序>>A=[011-5;-217-7;-426-10];V0=[1,1,1]';[k,lambda,Vk,Wc]=wfmifa1(A,V0,1.001,0.001,100)運行后屏幕顯示結(jié)果請注意迭代次數(shù)k,特征值的近似值lambda,對應的特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：k=lambda=Wc=51.002000000001381.376344154436924e-006Vk’=-0.50000000000000-0.50000000000000-1.00000000000000同理可得，另外與兩個特征值對應的特征向量，將計算結(jié)果列表5–3.表5–3例題5.3.3和例題5.3.4計算結(jié)果比較由（5.27）式計算結(jié)果由（5.28）式計算結(jié)果迭代次數(shù)k和其近似向量Vk的相鄰迭代誤差Wc和其近似值lambdak=5=(0.500,0.500,1.000)Vk=(-0.500,-0.500,-1.000)Wc=1.376344154436924e-006lambda=1.00200000000138k=5=(0.500,0.500,1.000)Vk=(-0.500,-0.500,-1.000)Wc=1.378794763695562e-009lambda=1.00200000000000迭代次數(shù)k和其近似向量Vk的相鄰迭代誤差Wc和其近似值lambdak=2=(-0.250,-0.500,-1.000)Vk=(-0.250,-0.500,-1.000)Wc=6.252343455491521e-004lambda=2.00199999968703k=2=(-0.250,-0.500,-1.000)Vk=(-0.250,-0.500,-1.000)Wc=3.131258787820493e-007lambda=2.00200000000016迭代次數(shù)k和其近似向量Vk的相鄰迭代誤差Wc和其近似值lambdak=2=(-0.400,-0.600,-1.000)Vk=(0.400,0.600,1.000)Wc=3.997600240971801e-004lambda=4.00199999980030k=2=(-0.400,-0.600,-1.000)Vk=(0.400,0.600,1.000)Wc=1.996005395108913e-007lambda=4.00199999999990通過比較表5–3中列出的例題5.3.3和例題5.3.4的計算結(jié)果可見，由原點位移反冪法的迭代公式（5.27）的迭代序列的收斂速度比迭代公式（5.28）稍慢些.（2）再用兩種原點位移反冪法的MATLAB主程序，求對應的特征向量.輸入MATLAB程序>>A=[011-5;-217-7;-426-10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,0.99999999999997,0.001,100)運行后屏幕顯示結(jié)果請注意：因為A-aE的各階主子式都不等于零，所以A-aE能進行Lu分解.A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：hl=-0.999999999999976.000000000000450.00000000000010Vk=0.50000000000000Vk=0.500000000000000.500000000000001.00000000000000Wc=4.317692037236759e-0131.039168751049192e-013k=2lambda=1.00000000000000輸入MATLAB程序>>A=[011-5;-217-7;-426-10];V0=[1,1,1]';[k,lambda,Vk,Wc]=wfmifa1(A,V0,0.99999999999997,0.001,100)運行后屏幕顯示結(jié)果Vk=0.50000000000000Vk=0.500000000000000.500000000000001.00000000000000k=3lambda=1.00000000000000Wc=5.412337245047640e-0165.4雅可比(Jacobi)方法及其MATLAB程序雅可比方法是用來計算實對稱矩陣的全部特征值和對應的特征向量的一種迭代的方法，最早由雅可比給出.自從計算機出現(xiàn)以后，古典的雅可比方法已有了不少的改進和推廣.5.4.3雅可比方法的MATLAB程序設階實對稱矩陣的個特征值為，特征值對應的特征向量為，則我們可以用下面的MATLAB程序計算的全部特征值和對應的特征向量的近似值和向量.用雅可比方法計算對稱矩陣的特征值和對應的特征向量的MATLAB主程序輸入的量：階實對稱矩陣、計算的精度jd、迭代的最大次數(shù)max1；輸出的量：計算全過程中每次迭代的序號k,選取的旋轉(zhuǎn)主元mk及其所在的行數(shù)和列數(shù)，c，t，pii和pij的值，旋轉(zhuǎn)矩陣Pk，正交矩陣Vk=P1P2…Pk，對稱矩陣Bk,判斷是否滿足控制迭代終止的條件Wc，以特征向量為列向量的矩陣V,特征值為對角元的對角矩陣D.如果迭代次數(shù)已經(jīng)達到最大的迭代次數(shù)max1，則給出提示的相關(guān)信息.根據(jù)雅可比迭代求階實對稱矩陣A的特征值和對應的特征向量的近似向量，精度為的一般步驟，現(xiàn)提供MATLAB主程序如下：function[k,Bk,V,D,Wc]=renyujjacobite(A,jd,max1)[n,n]=size(A);Vk=eye(n);Bk=A;state=1;k=0;P0=eye(n);Aij=abs(Bk-diag(diag(Bk)));[m1i]=max(Aij);[m2j]=max(m1);i=i(j);while((k<=max1)&(state==1))k=k+1,aij=abs(Bk-diag(diag(Bk)));[m1i]=max(abs(aij));[m2j]=max(m1);i=i(j),j,Aij=(Bk-diag(diag(Bk)));mk=m2*sign(Aij(i,j)),Wc=m2,Dk=diag(diag(Bk));Pk=P0;c=(Bk(j,j)-Bk(i,i))/(2*Bk(i,j)),t=sign(c)/(abs(c)+sqrt(1+c^2)),pii=1/(sqrt(1+t^2)),pij=t/(sqrt(1+t^2)),Pk(i,i)=pii;Pk(i,j)=pij;Pk(j,j)=pii;Pk(j,i)=-pij;Pk,B1=Pk'*Bk;B2=B1*Pk;Vk=Vk*Pk,Bk=B2,if(Wc>jd)state=1;elsereturnendPk;Vk;Bk=B2;Wc;endif(k>max1)disp('請注意迭代次數(shù)k已經(jīng)達到最大迭代次數(shù)max1,迭代次數(shù)k,對稱矩陣Bk,以特征向量為列向量的矩陣V,特征值為對角元的對角矩陣D如下：')elsedisp('請注意迭代次數(shù)k,對稱矩陣Bk,以特征向量為列向量的矩陣V,特征值為對角元的對角矩陣D如下：')endWc;k=k;V=Vk;Bk=B2;D=diag(diag(Bk));[V1,D1]=eig(A,'nobalance')例5.4.2用雅可比方法的MATLAB程序計算例5.4.1中矩陣的特征值和對應的特征向量（）.解（1）保存名為renjacobite.m為M文件；（2）輸入MATLAB程序>>A=[12-563-1;-56720;3251;-10112];[k,B,V,D,Wc]=renyujjacobite(A,0.001,100)（3）運行后屏幕顯示計算的全過程中每次迭代的序號k,選取的旋轉(zhuǎn)主元mk及其所在的行數(shù)和列數(shù)，c，t，pii和pij的值，旋轉(zhuǎn)矩陣Pk，正交矩陣Vk=P1P2…Pk，對稱矩陣Bk,判斷是否滿足控制迭代終止的條件Wc，以特征向量為列向量的矩陣V,特征值為對角元的對角矩陣D和相關(guān)信息如下：k=i=j=mk=Wc=121-5656c=t=-0.04464285714286-0.95635313919972pii=pij=0.72270271801843-0.69115901308510Pk=0.722702718018430.6911590130851000-0.691159013085100.7227027180184300001.0000000000000000001.00000000000000Vk=0.722702718018430.6911590130851000-0.691159013085100.7227027180184300001.0000000000000000001.00000000000000Bk=65.5557757951845600.78579012788509-0.72270271801843-0.00000000000001-46.555775795184563.51888247529217-0.691159013085100.785790127885093.518882475292175.000000000000001.00000000000000-0.72270271801843-0.691159013085101.0000000000000012.00000000000000k=i=j=mk=Wc=2323.518882475292173.51888247529217c=t=-7.32558932518824-0.06793885568129pii=pij=0.99770011455446-0.06778260409592Pk=1.0000000000000000000.997700114554460.0677826040959200-0.067782604095920.9977001145544600001.00000000000000Vk=0.722702718018430.689569426530350.046848557751270-0.691159013085100.721040584555810.0489866722144900-0.067782604095920.9977001145544600001.00000000000000Bk=65.55577579518456-0.053262901140920.78398290060672-0.72270271801843-0.05326290114093-46.794844643832850-0.757352030626270.783982900606720.000000000000005.239068848648290.95085155680318-0.72270271801843-0.757352030626270.9508515568031812.00000000000000k=i=j=mk=Wc=3430.950851556803180.95085155680318c=t=-3.55519802380213-0.13796227443116pii=pij=0.99061693994324-0.13666776612460Pk=1.0000000000000000001.0000000000000000000.990616939943240.1366677661246000-0.136667766124600.99061693994324Vk=0.722702718018430.689569426530350.046408974920320.00640268773403-0.691159013085100.721040584555810.048527027327120.006694899061430-0.067782604095920.988338634460960.1363534459184200-0.136667766124600.99061693994324Bk=65.55577579518456-0.053262901140920.87539690801061-0.60877636330628-0.05326290114093-46.794844643832850.10350561019562-0.750245751038800.875396908010610.103505610195625.10788720522532-0.00000000000000-0.60877636330628-0.75024575103880-0.0000000000000012.13118164342297k=i=j=mk=Wc=4130.875396908010610.87539690801061c=t=-34.52598931799430-0.01447880833914pii=pij=0.99989519853186-0.01447729093877Pk=0.999895198531860-0.01447729093877001.00000000000000000.0144772909387700.9998951985318600001.00000000000000Vk=0.723298853944650.689569426530350.035941333680620.00640268773403-0.690384038712800.721040584555810.058528051740800.006694899061430.01430846595712-0.067782604095920.988235055121050.13635344591842-0.001978579012140-0.136653443142060.99061693994324Bk=65.56845049923633-0.05175883827808-0.00000000000000-0.60871256264964-0.05175883827809-46.794844643832850.10426586517177-0.75024575103880-0.000000000000000.104265865171775.095212501173560.00881343252823-0.60871256264964-0.750245751038800.0088134325282312.13118164342297k=i=j=mk=Wc=542-0.750245751038800.75024575103880c=t=39.271149623750840.01272992971264pii=pij=0.999918984291140.01272889838836Pk=1.0000000000000000000.999918984291140-0.01272889838836001.00000000000000000.0127288983883600.99991898429114Vk=0.723298853944650.689595059736030.03594133368062-0.00237529014628-0.690384038712800.721067387631600.05852805174080-0.002483695665250.01430846595712-0.066041483482200.988235055121050.13720519702737-0.001978579012140.01260946237032-0.136653443142060.99053668440964Bk=65.56845049923633-0.05950288535679-0.00000000000000-0.60800441437674-0.05950288535680-46.804395219510780.104369603285900.00000000000000-0.000000000000000.104369603285905.095212501173560.00748552889860-0.608004414376740.000000000000000.0074855288986012.14073221910090k=i=j=mk=Wc=641-0.608004414376740.60800441437674c=t=-43.93694931878409-0.01137847012503pii=pij=0.99993527149402-0.01137773361366Pk=0.99993527149402000.0113777336136601.0000000000000000001.000000000000000-0.01137773361366000.99993527149402Vk=0.723279061308990.689595059736030.035941333680620.00585436528595-0.690311092357770.721067387631600.05852805174080-0.010338540582940.01274645560931-0.066041483482200.988235055121050.13735911385404-0.013248513471450.01260946237032-0.136653443142060.99045005670500Bk=65.57536865930122-0.05949903382392-0.00008516835377-0.00000000000000-0.05949903382393-46.804395219510780.10436960328590-0.00067700797883-0.000085168353770.104369603285905.095212501173560.00748504437150-0.00000000000000-0.000677007978830.0074850443715012.13381405903603k=i=j=mk=Wc=7320.104369603285900.10436960328590c=t=-2.486337309269764e+002-0.00201098208240pii=pij=0.99999797798167-0.00201097801616Pk=1.0000000000000000000.999997977981670.0020109780161600-0.002010978016160.9999979779816700001.00000000000000…………請注意迭代次數(shù)k,對稱矩陣Bk,以特征向量為列向量的矩陣V,特征值為對角元的對角矩陣D如下：V1=0.68990429476497-0.037324232224840.00588594854431-0.722913771734500.72058252860300-0.05998661236737-0.010283221619770.69069289931337-0.06802029759277-0.987953684104720.13841044442471-0.012779125692250.012888857681930.137680884982000.990304074432190.01325486405899D1=-46.8046366141973600005.09541442877727000012.13382202426702000065.57540016115307k=10B=65.575400160459450.00000000000175-0.000204819675660.000000148628360.00000000000175-46.804636614197390.000000627399840.00000000000000-0.000204819675660.000000627399845.09541442947090-0.000000000007370.00000014862836-0.00000000000000-0.0000000000073712.13382202426704V=0.722913898115070.689904295216170.037321775686890.00588595055487-0.690692696132010.720582529328160.05998894273570-0.010283223540620.01278247108107-0.068020285649770.987953641643790.13841044446122-0.013255333078980.01288885601755-0.137680840249460.99030407439520D=65.575400160459450000-46.8046366141973900005.09541442947090000012.13382202426704Wc=6.920584967017158e-0045.5豪斯霍爾德(Householder)方法及其MATLAB程序在雅可比方法中，每一次迭代產(chǎn)生兩個值為零的非對角元，但接下來的迭代使它們變成非零.因此需要許多次迭代才能使得非對角元足夠接近零.豪斯霍爾德(Householder)方法可以在一次迭代中產(chǎn)生多個值為零的非對角元，而且在接下來的迭代中使它們保持為零.5.5.1豪斯霍爾德方法及其MATLAB程序設向量，確定一個n階初等反射矩陣，使得的后分量個分量都為零，則我們可以用下面的MATLAB程序計算.求初等反射矩陣，使得的第一個分量以外的其余的分量都為零的MATLAB主程序輸入的量：維向量；輸出量：迭代次數(shù)k,對稱矩陣B,以特征向量為列向量的矩陣V,特征值為對角元的對角矩陣D.如果迭代次數(shù)已經(jīng)達到最大的迭代次數(shù)max1，則給出提示的相關(guān)信息.根據(jù)（5.49）至（5.52）式，現(xiàn)提供MATLAB主程序如下：function[xigema,rou,miou,P,PX]=renHouseholder(X)n=size(X);nX=norm(X,2);xigema=nX*sign(X(1));rou=xigema*(xigema+X(1));miou=[xigema,zeros(1,n-1)]'+X,E=eye(n,n);C=2*miou*(miou)';P=E-C/(norm(miou,2)^2);PX=P*X;例5.5.1設向量，確定一個初等反射矩陣，使得的后兩個分量為零.解輸入MATLAB程序>>X=[221]';[xigema,rou,miou,P,PX]=renHouseholder(X)運行后屏幕顯示結(jié)果P=PX=-0.6667-0.6667-0.3333-3.0000-0.66670.7333-0.13330.0000-0.3333-0.13330.93330.00005.5.2矩陣約化為上豪斯霍爾德（Hou