第五章矩阵的特征值与特征向量的计算

1、第五章 矩阵的特征值与特征向量的计算 5.2 幂法及其MATLAB程序5.2.2 幂法的MATLAB程序用幂法计算矩阵的主特征值和对应的特征向量的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(kjd)state=1;endk=k+1;Wc=Wc;endif(WcA=1 -1;2 4; V0=1,1; k,lambda,Vk,Wc=mifa(A,V0,0.00001,100), V,D = eig (A), Dzd=max(diag(D),

2、 wuD= abs(Dzd- lambda), wuV=V(:,2)./Vk, 运行后屏幕显示结果请注意：迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：k = lambda = Wc = 33 3.804 8.4999e-007Vk = V = wuV =-0.432 -0.655 0.996 -0.2941.000 0.655 -0.992 -0.992Dzd = wuD = 3 1.6435e-006 由输出结果可看出，迭代33次，相邻两次迭代的误差Wc 8.69 19e-007，矩阵的主特征值的近似值lambda3.000 00和对应的特

3、征向量的近似向量Vk （-0.500 00，1.000 00, lambda与例5.1.1中的最大特征值近似相等，绝对误差约为1.738 37e-006，Vk与特征向量 的第1个分量的绝对误差约等于0，第2个分量的绝对值相同.由wuV可以看出，的特征向量V(:,2) 与Vk的对应分量的比值近似相等.因此，用程序mifa.m计算的结果达到预先给定的精度.(2) 输入MATLAB程序 B=1 2 3;2 1 3;3 3 6; V0=1,1,1; k,lambda,Vk,Wc=mifa(B,V0,0.00001,100), V,D = eig (B), Dzd=max(diag(D), wuD= a

4、bs(Dzd- lambda), wuV=V(:,3)./Vk,运行后屏幕显示结果请注意：迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：k = lambda = Wc = Dzd = wuD = 3 9 0 9 0Vk = wuV = 0.000 0.773 0.000 0.773 1.000 0.773V = 0.655 0.963 0.386 -0.655 0.963 0.386 0 -0.963 0.773 (3) 输入MATLAB程序 C=1 2 2;1 -1 1;4 -12 1;V0=1,1,1; k,lambda,Vk,Wc=mif

5、a(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 = 100 0.910 2.119 Dzd = wuD = 1.001 0.091Vk= Vzd = wuV =0.993 0.329 0.3350.995 0.776 0.7781.000 -0.776 -0.776由输出

6、结果可见，迭代次数k已经达到最大迭代次数max1=100，并且lambda的相邻两次迭代的误差Wc2.377 582,由wuV可以看出，lambda的特征向量Vk与真值Dzd的特征向量Vzd对应分量的比值相差较大，所以迭代序列发散.实际上，实数矩阵C的特征值的近似值为，并且对应的特征向量的近似向量分别为=（0.329，0.776，-0.776），（-0.001,-0.300-0.100i, 0.801-0.400i），( -0.001, -0.300 + 0.100i, 0.801 + 0.400i) , 是常数）.（4）输入MATLAB程序 D=-4 14 0;-5 13 0;-1 0 2;

7、 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 = 19 6.528 6.1684e-006Dtzd = wuDt = 6.000 6.0768e-006Vk = Vzd = wuV =0.564 0.564 0.564 0.8

8、86 0.117 0.618 -0.180 -0.391 0.3705.3 反幂法和位移反幂法及其MATLAB程序5.3.3 原点位移反幂法的MATLAB程序（一） 原点位移反幂法的MATLAB主程序1用原点位移反幂法计算矩阵的特征值和对应的特征向量的MATLAB主程序1function 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); if RA1=0disp(请注意：因为A-aE的n阶行列式hl等于零，所以A-aE不能进行LU分解.)retur

9、nendlambda=0;if RA1=0 for p=1:nh(p)=det(A1(1:p, 1:p);endhl=h(1:n);for i=1:nif h(1,i)=0disp(请注意：因为A-aE的r阶主子式等于零，所以A-aE不能进行LU分解.) returnendend if h(1,i)=0 disp(请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.)k=1;Wc =1;state=1; Vk=V0;while(kjd)state=1;endk=k+1;%Vk=Vk2,mk=mk1,endif(Wc A=1 -1 0;-2 4 -2;0 -1 2;V0=1,

10、1,1;k,lambda,Vk,Wc=ydwyfmf(A,V0,0.2,0.0001,10000)运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：k = lambda = Wc = hl = 3 0.2384 1.0213e-007 0.8000 1.0400 0.2720Vk = V = D = 1.0000 -0.2424 -1.0000 -0.5707 5.1249 0 0 0.7616 1.000

11、0 -0.7616 0.3633 0 0.2384 0 0.4323 -0.3200 -0.4323 1.0000 0 0 1.6367（2）输入MATLAB程序 A=1 -1;2 4;V0=20,1;k,lambda,Vk,Wc=ydwyfmf(A,V0,2.001,0.0001,100) 运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：k = lambda = Wc = hl = 2 2.0020 5

12、.1528e-007 -1.0010 -0.0010Vk = V = D = 1.0000 -1.0000 0.5000 2 0 -1.0000 1.0000 -1.0000 0 3（3）输入MATLAB程序 A=-11 2 15;2 58 3;15 3 -3;V0=1,1,-1;k,lambdan,Vk,Wc=ydwyfmf(A,V0,8.26, 0.0001,100) 运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代

13、的误差Wc如下：k = lambdan= Wc = hl = 2 8.2640 6.9304e-008 -19.2600 -961.9924 -6.1256Vk = V = D = -0.7692 0.7928 0.6081 0.0416 -22.5249 0 0 0.0912 0.0030 -0.0721 0.9974 0 8.2640 0 -1.0000 -0.6095 0.7906 0.0590 0 0 58.2609例 5.3.3 用原点位移反幂法的迭代公式（5.28），计算的分别对应于特征值，的特征向量，的近似向量，相邻迭代误差为0.001.将计算结果与精确特征向量比较.解 （1）计

14、算特征值对应的特征向量的近似向量.输入MATLAB程序 A=0 11 -5;-2 17 -7;-4 26 -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,运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代

15、的误差Wc如下：hl = -1.000 5.000 -0.000k = lambda = RA1 = 5 1.000 -0.000Vk = VD = wuV = -0.000 -0.386 0.773 -0.000 -0.386 0.773 -1.000 -0.773 0.773Wc = Dzd = wuD = 1.5562e-009 1.000 0.000 从输出的结果可见，迭代5次，特征向量的近似向量的相邻两次迭代的误差Wc1.379 e-009，由wuV可以看出，= Vk与VD 的对应分量的比值相等.特征值的近似值lambda 1.002与初始值1.001的绝对误差为0.001，而与的绝

16、对误差为0.002，其中,.（2）计算特征值对应特征向量的近似向量.输入MATLAB程序 A=0 11 -5;-2 17 -7;-4 26 -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,运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差W

17、c如下：hl = -2.000 -8.000 0.000k = Wc = lambda = WD = 2 3.2120e-007 2.016 0.016Vk = VD = wuV = -0.999 0.599 -0.401 -0.999 0.198 -0.398 -1.000 0.397 -0.397从输出的结果可见，迭代2次，特征向量的近似向量的相邻两次迭代的误差Wc3.131e-007，与的对应分量的比值近似相等.特征值的近似值lambda2.002与初始值2.001的绝对误差约为0.001，而lambda与的绝对误差约为0.002，其中,.（3）计算特征值对应特征向量的近似向量.输入MA

18、TLAB程序 A=0 11 -5;-2 17 -7;-4 26 -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,运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：hl =-4.000 -30.000 -0.999k = la

19、mbda = Wc = WD = 2 4.990 1.4842e-007 0.990Vk = VD = wuV = 0.001 -0.153 -0.380 0.001 -0.229 -0.381 1.000 -0.381 -0.381从输出的结果可见，迭代2次，特征向量的近似向量的相邻两次迭代的误差Wc1.996e-007，与的对应分量的比值近似相等.特征值的近似值的绝对误差近似为，而lambda与的绝对误差约为0.002，其中-0.400 000 000 000 00，-0.600 000 000 000 00，-1.000 000 000 000 00, .（二）原点位移反幂法的MATLA

20、B主程序2用原点位移反幂法计算矩阵的特征值和对应的特征向量的MATLAB主程序2function 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(kjd)state=1;endk=k+1;endif(Wc A=0 11 -5;-2 17 -7;-4 26 -10;V0=1,1,1; k,lambda,Vk,Wc=wfmifa1(A,V0,1.001,0.001,100)

21、运行后屏幕显示结果请注意迭代次数k,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：k = lambda = Wc =5 1.138 1.6924e-006Vk = -0.000 -0.000 -1.000同理可得，另外与两个特征值对应的特征向量.（2）再用两种原点位移反幂法的MATLAB主程序，求对应的特征向量.输入MATLAB程序 A=0 11 -5;-2 17 -7;-4 26 -10;V0=1,1,1;k,lambda,Vk,Wc=ydwyfmf(A,V0,0.997,0.001,100)运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零

22、，所以A-aE能进行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：hl = -0.997 6.045 0.010Vk = 0.000 0.000 1.000Wc = 4.6759e-013RA1 = 1.9192e-013k = 2lambda = 1.000输入MATLAB程序 A=0 11 -5;-2 17 -7;-4 26 -10;V0=1,1,1;k,lambda,Vk,Wc=wfmifa1(A,V0, 0.997,0.001,100)运行后屏幕显示结果Vk = 0.000

23、 0.000 1.000请注意迭代次数k,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：k = 3lambda = 1.000Wc = 5.7640e-0165.4 雅可比(Jacobi)方法及其MATLAB程序5.4.3 雅可比方法的MATLAB程序用雅可比方法计算对称矩阵的特征值和对应的特征向量的MATLAB主程序function k,Bk,V,D,Wc=jacobite(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);m1 i=m

24、ax(Aij);m2 j=max(m1);i=i(j);while (kjd)state=1;else returnendPk;Vk;Bk=B2;Wc;endif(kmax1)disp(请注意迭代次数k已经达到最大迭代次数max1,迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D如下：)elsedisp(请注意迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D如下：) endWc;k=k; V=Vk;Bk=B2;D=diag(diag(Bk);V1,D1 =eig(A,nobalance)例5.4.2 用雅可比方法的MATLAB程序

25、计算矩阵的特征值和对应的特征向量（）.解 （1）保存名为jacobite.m为M文件；（2）输入MATLAB程序 A=12 -56 3 -1;-56 7 2 0;3 2 5 1;-1 0 1 12;k,B,V,D,Wc=jacobite(A,0.001,100)（3）运行后屏幕显示如下：k = i = j = mk = Wc =1 2 1 -56 56c = t = -0.286 -0.972pii = pij = 0.843 -0.510Pk = 0.843 0.510 0 0 -0.510 0.843 0 0 0 0 1.000 0 0 0 0 1.000Vk = 0.843 0.510

26、0 0 -0.510 0.843 0 0 0 0 1.000 0 0 0 0 1.000Bk = 65.456 0 0.509 -0.843 -0.001 -46.456 3.217 -0.510 0.509 3.217 5.000 1.000 -0.843 -0.510 1.000 12.000k = i = j = mk = Wc = 2 3 2 3.217 3.217c = t = -7.824 -0.129pii = pij = 0.446 -0.592Pk = 1.000 0 0 0 0 0.446 0.592 0 0 -0.592 0.446 0 0 0 0 1.000Vk = 0

27、.843 0.035 0.127 0 -0.510 0.581 0.449 0 0 -0.592 0.446 0 0 0 0 1.000Bk = 65.456 -0.092 0.672 -0.843 -0.093 -46.285 0 -0.627 0.672 0.000 5.829 0.318 -0.843 -0.627 0.318 12.000k = i = j = mk = Wc = 3 4 3 0.318 0.318c = t = -3.213 -0.116pii = pij = 0.324 -0.460Pk = 1.000 0 0 0 0 1.000 0 0 0 0 0.324 0.4

28、60 0 0 -0.460 0.324Vk = 0.843 0.035 0.032 0.403 -0.510 0.581 0.712 0.143 0 -0.592 0.096 0.842 0 0 -0.460 0.324Bk = 65.456 -0.092 0.061 -0.628 -0.093 -46.285 0.562 -0.880 0.061 0.562 5.532 -0.000 -0.628 -0.880 -0.000 12.297k = i = j = mk = Wc = 4 1 3 0.061 0.061c = t = -34.430 -0.914pii = pij = 0.186

29、 -0.877Pk = 0.186 0 -0.877 0 0 1.000 0 0 0.877 0 0.186 0 0 0 0 1.000Vk = 0.465 0.035 0.062 0.403 -0.280 0.581 0.080 0.143 0.712 -0.592 0.105 0.842 -0.214 0 -0.206 0.324Bk = 65.633 -0.808 -0.000 -0.964 -0.809 -46.285 0.177 -0.880 -0.000 0.177 5.356 0.823 -0.964 -0.880 0.823 12.297k = i = j = mk = Wc

30、= 5 4 2 -0.880 0.880c = t = 39.084 0.264pii = pij = 0.114 0.836Pk = 1.000 0 0 0 0 0.114 0 -0.836 0 0 1.000 0 0 0.836 0 0.114Vk = 0.465 0.603 0.062 -0.628 -0.280 0.160 0.080 -0.525 0.712 -0.220 0.105 0.737 -0.214 0.032 -0.206 0.964Bk = 65.633 -0.679 -0.000 -0.674 -0.680 -46.078 0.590 0.000 -0.000 0.5

31、90 5.356 0.860 -0.674 0.000 0.860 12.090k = i = j = mk = Wc = 6 4 1 -0.674 0.674c = t = -43.409 -0.503pii = pij = 0.402 -0.366Pk = 0.402 0 0 0.366 0 1.000 0 0 0 0 1.000 0 -0.366 0 0 0.402Vk = 0.899 0.603 0.062 0.595 -0.777 0.160 0.080 -0.294 0.931 -0.220 0.105 0.404 -0.145 0.032 -0.206 0.500Bk = 65.

32、122 -0.392 -0.377 -0.000 -0.393 -46.078 0.590 -0.883 -0.377 0.590 5.356 0.150 -0.000 -0.883 0.150 12.603k = i = j = mk = Wc = 7 3 2 0.590 0.590c = t = -2.9764e+002 -0.240pii = pij = 0.167 -0.616Pk = 1.000 0 0 0 0 0.167 0.616 0 0 -0.616 0.167 0 0 0 0 1.000请注意迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D如下

33、：V1 = 0.497 -0.484 0.431 -0.450 0.300 -0.737 -0.977 0.337 -0.277 -0.472 0.471 -0.225 0.193 0.200 0.219 0.899D1 = -46.736 0 0 0 0 5.727 0 0 0 0 12.702 0 0 0 0 65.307k = 10B = 65.945 0.175 -0.566 0.836 0.175 -46.739 0.984 0.000 -0.566 0.984 5.090 -0.737 0.836 -0.000 -0.737 12.704V = 0.507 0.617 0.689

34、0.487 -0.201 0.816 0.570 -0.062 0.107 -0.977 0.379 0.122 -0.898 0.755 -0.946 0.520D = 65.945 0 0 0 0 -46.739 0 0 0 0 5.090 0 0 0 0 12.704Wc = 6.7158e-0045.5 豪斯霍尔德(Householder)方法及其MATLAB程序5.5.1 豪斯霍尔德方法及其MATLAB程序求初等反射矩阵，使得的第一个分量以外的其余的分量都为零的MATLAB主程序function xigema,rou,miou,P,PX=Householder(X)n=size(X);nX=norm(X,2); xigema=nX*sign(X(1);rou=xigema*(x