3矩阵特征值与特征向量的计算

3.2
幂法是求方阵的最大特征值及对应特征向量的一种迭代法。
3.2.1
设An有n个线性相关的特征向量v1，v2，…，vn，对应的特征值1，2，…，n，满足
|1| > |2|…|n| (3.2.1)
1.
因为{v1，v2，…，vn}为Cn的一组基，所以任给x(0)0， ——线性表示
所以有
(3.2-2)
若a10，则因 知，当k充分大时A(k)x(0)1ka1v1=cv1属1的特征向量
另一方面，记max(x) =xi，其中|xi| = ||x||，则当k充分大时，
若a1= 0，则因舍入误差的影响，会有某次迭代向量在v1方向上的分量不为0，迭代下去可求得1及对应特征向量的近似值。
2.
在实际计算中，若|1| > 1则|1ka1|，若|1| < 1则|1ka1|0都将停机。
须采用“规范化”的方法
输出0
例2用幂法求方阵A的最大模特征值，并用Aitkem加速法
解：见（P314）
x0=[1;1;1];y0=x0/maxa(x0)
x1=A*y0;y1=x1/maxa(x1)
x2=A*y1;y2=x2/maxa(x2)
l0=maxa(x2)-(maxa(x2)-maxa(x1))^2/(maxa(x2)-2*maxa(x1)+maxa(x0))
A()的盖氏圆为：
因为A(0) =D的n个特征值a11，a12，…，ann，恰为A的盖氏圆圆心，当由0增大到1时，i()画出一条以i(0) =aii为始点，i(1) =i为终点的连续曲线，且始终不会越过Gi；
不失一般性，设A开头的k个圆盘是连通的，其并集为S，它与后n–k个圆盘严格分离，显然，A()的前k个盖氏圆盘与后n–k个圆盘严格分离。
方法：设A1=Q1R1（QR分解），令A2=R1Q1，设A2=Q2R2，令A3=R2Q2，即
k= 1，2，…(3.3-1)
{Ak}的性质：
①Ak~A：Ak+1=RkQk= (Qk-1Ak)Qk=Rk=Qk-1Ak
3)因为例1中A为实方阵，所以若为A的特征值，则 也是A的特征值，所以G2，G4中各有一个实特征值。
3.1.3
由于特征值是相似不变量，所以代数上常用相似变换将矩阵化简以得到特征向量，这里也可用相似变换将盖氏圆的半径变小，以得到更好的估计。
原理：取对角阵作相似变换阵：P= diag(b1，b2，…，bn)其中bi> 0，i= 1，2，…，n
若1–a为A–aI的最大模特征值，且 。（k–a是A–aI的次最大模特征值），则对A–aI计算1–a及对应的特征向量比对 计算收敛得快，此即为原点平移法。
计算1–a及特征向量的迭代公式
特征向量： ，max(x(k))1–a，a+ max(x(k))1。
注：a的选取较为困难。
例3设 ， ，求最大模特征值及特征向量。
迭代公式：x(k+1)=A–1x(k)，k= 0，1，2，…，
但A–1不易求，通常可解方程组Ax(k+1)=x(k)来求x(k+1)
即有 (3.2.12)
当k时有
注：为解(3.2-12)中的方程组。对 作LR分解（带行交模）PA=LR
则有
2.
——反幂法结合原点平移法
思想：若已知 为j的近似值，则 的特征值是
y
maxa(x1)-4
3.
定义设A对称，x0，则称 为 关于 的Rayleigh商
思想：A对称 特征值1，2，…，n均为实数，且存在特征向量v1，v2，…，vn为标准正交基。
设 ，a10，则
当k充分大时，M'与k无关)
注；此比Aitken加速中的(3.2-6)更快
公式 称为Rayleigh商加速法。
证明：设Ax=x(x0)，若k使得
因为
例1估计方阵 特征值的范围
解：
G1= {z：|z–1|0.6}；G2= {z：|z–3|0.8}；
G3= {z：|z+ 1|1.8}；G4= {z：|z+ 4|0.6}。
注：定理称A的n个特征值全落在n个盖氏圆上，但未说明每个圆盘内都有一个特征值。
3.1.2
称相交盖氏圆之并构成的连通部分为连通部分。
x0=[1;1;1];
B=A+6.42*eye(3);
C=lu(B);
R = triu(C,0);
L =eye(3)+ tril(C,-1);
y=x0/maxa(x0);
z=[1,1,1]';
x1=inv(R)*z
while(abs(maxa(x1)-maxa(x0)))>0.001
x0=x1;
y=x0/maxa(x0)
而显然 非常大（最大）比值 很小
迭代公式：
当k时有
注：(1)若有分LR解
则迭代公式
(3.2-16)
(2)在(3.2-16)中直接取z(1)= (1，…，1)T作初值开始迭代称为半次迭代法
例5设 的一个特征值的的近似值 ，用带原点平移的反幂法求及相应的特征向量
见[P320]
format long;
A=[-1,2,1;2,-4,1;1,1,-6];
解：（P315）
幂法：
A=[-3,1,0;1,-3,-3;0,-3,4];
x0=[0;0;1];k=1;
y=x0/maxa(x0)
x1=A*y
while(abs(maxa(x1)-maxa(x0)))>0.001
x0=x1;
y=x0/maxa(x0)
k=k+1
x1=A*y
end;
y
maxa(x1)
原点平移法：
其中
注：有了R(x(k))，R(x(k+1))，R(x(k+2))，的值，可再用Aitken加速法得到 的一个更好的近似值：
因为
所以
例4设 ，用Rayleigh商加速法求 的最大模特征值及特征向量，并与幂法相比较。
解：（P317）
幂法：
A=[6,2,1;2,3,1;1,1,1];
x0=[1;1;1];k=1;
y=x0/maxa(x0)
x1=A*y
while(abs(r1-r))>0.001
x0=x1;r1=r;
y=x0/maxa(x0)
x1=A*y
r = y'*x1/(y'*y)
k=k+1
end;
y
maxa(x1)
r
3.2.3
——用 代替 作幂法，即反幂法
1.
若 可逆，|1| > |2|…|n|为其特征值，则 为A-1的最大模特征值。
当= 0时，A(0) =D的前k个特征值刚好落在前k个圆盘G1，…，Gk中，而另n–k个特征值则在区域S之外，从0变到1时， 与 始终分离（严格）。连续曲线始终在S中，所以S中有且仅有A的k个特征值。
注：1)每个孤立圆中恰有一个特征值。
2)例1中G2，G4为仅由一个盖氏圆构成的连通部分，故它们各有一个特征值，而G1，G3构成的连通部分应含有两个特征值。
G3= {z：|z–0.4|0.03}
3G3，较好。
为了更好地估计另外两个特征值可取b3最小：
取b1=b2= 1，b3= 0.1即 ，则
所以G1'= {z：|z–0.9|0.022}；G2'= {z：|z–0.8|0.023}；G3'= {z：|z–0.4|0.3}
三个盖氏圆分离，故有1G1'，2G2'，3G3。
A=[-3,1,0;1,-3,-3;0,-3,4];
x0=[0;0;1];k=1;
y=x0/maxa(x0)
x1=(A+4*eye(3))*y
while(abs(maxa(x1)-maxa(x0)))>0.001
x0=x1;
y=x0/maxa(x0)
k=k+1
x1=(A+4*eye(3))*y
end;
z=inv(L)*y
x1=inv(R)*z
end;
-6.42+1/maxa(x1)
预备知识：矩阵论
1.
定理设 可逆，则存在正交阵Q与上三角阵R使A=QR
注：方法1)使用史密斯正交变换
2)使用Householder变换（反）
3)使Givens变换
2.
定理设 ，则存在正交阵Q使 ——实Schur型
其中 至多2阶。若1阶，其元素即A的特征值若2阶其特征值为A的一对共轭复特征值。
y=x0/maxa(x0)
x1=A*y
while(abs(maxa(x1)-maxa(x0)))>0.001
x0=x1;
y=x0/maxa(x0)
x1=A*y
end;
y
maxa(x1)
3.2.2
幂法的迭代公式：
当k时 ，max(x(k))1，其中|1| > |2|…|n|
注：幂法的收敛速度取决于比值|2| / |1|，考虑收敛加速
第
一些工程技术问题需要用数值方法求得矩阵的全部或部分特征值及相关的特征向量。
3.1
较粗估计(A)||A||
欲将复平面上的特征值一个个用圆盘围起来。
3.1.1
定义3.1-1设A= [aij]nn，称由不等式 所确定的复区域为A的第i个盖氏图，记为Gi，i= 1，2，…，n。
定理3.1-1若为A的特征值，则
孤立的盖氏圆本身也为一个连通部分。
定理3.1-2若由A的k个盖氏圆组成的连通部分，含且仅含A的k个特征值。
证明:令D= diag(a11，a12，…，ann)，M=A–D，记
则显然有A(1) =A，A(0) =D，易知A()的特征多项式的系数是的多项式，从而A()的特征值1()，2()，…，n()为的连续函数。
计算x1=A y0，max(x1)，y1= x1/max(x1)
x2 = A y1，1 =0
计算max(x2)
y2= x2/max(x2)
0=max(x2)-[max(x2)-max(x1)]^2/[max(x2)-2max(x1)+max(x0)]
x0 = x1，x1 = x2
|1 -0| > eps
数组x = [n]
kΒιβλιοθήκη Baidu= 1
for(i = 2 to n, i++)
若|x[i]| > |x[k]|
T
k = i
max = x[k]
幂法流程：
输入数组x0, eps, A
x1 = x0
y = x1/maxa(x1)
x0 = Ay
|maxa(x1)–maxa(x0)| < eps
输出y, maxa(x0)
while (abs(l1-l0))>0.01
x0=x1;x1=x2;l1=l0;
x2=A*y2
maxk=maxa(x2)
y2=x2/maxk
l0=maxa(x2)-(maxa(x2)-maxa(x1))^2/(maxa(x2)-2*maxa(x1)+maxa(x0))
end;
2.
思想：由矩阵论知，若为A的特征值则–a为A–aI的特征值，且特征向量相同。
例1，用幂法求 的最大模特征值及对应特征向量
见P312
function y = maxa(x)
k=1;n=length(x);
for i=2:n
if (abs(x(i))>abs(x(k)),k=i;
end;
end;
y=x(k);
A=[2,4,6;3,9,15;4,16,36];
x0=[1;1;1];
注：想加快迭代速度通常先将A化为上Hessenberg阵
3.
正交相似于一个n阶上Hessenberg矩阵
（ ）
证明：见（P125）
§
QR方法即使用QR分解构造迭代序列，是目前求一般矩阵全部特征值的最有效并广泛使用的方法之一。
3.3.1
思想：从A1=A出发用正交相似变换得序列{Ak}使当k时，Ak本质收敛到块上三角阵
y=x0/maxa(x0)
x1=A*y
while(abs(maxa(x1)-maxa(x0)))>0.001
x0=x1;
y=x0/maxa(x0)
x1=A*y
k=k+1
end;
y
maxa(x1)
Rayleigh商加速法：
A=[6,2,1;2,3,1;1,1,1];
x0=[1;1;1];k=1;r=0;
1.
(1)思想：由定理3.2.1的证明知
(3.2.6)
解之得
(3.2.7)
使用1(k+2)作为1的近似值的算法称为Aitken加速法。
(2) Aitken加速法
设{xk}线性收敛到x*，即存在c，|c| < 1，满足
xk+1–x*= (c–k)(xk–x*)，其中
令 则
算法：
计算
流程图
输入x0
计算max(x0)，y0 = x0/max(x0)
则 与A有相同特征值.
而B的第i个盖氏圆为： ，
适当选取b1，b2，…，bn就有可能使B的某些盖氏圆的半径比A的相应盖氏圆的半径小。
1)欲缩小Gi，可取bi最大。
2)欲缩小除Gi外的圆，可取bi最小。
例2，估计 的特征值范围。
解：A的三个盖氏圆分别为：
G1= {z：|z–0.9|0.13}；G2= {z：|z–0.8|0.14}；
，k= 0，1，2，…(3.2-4)
定理3.2-1任给初始向量 有， (3.2-5)
证明：
而
注：若 的特征值不满足条件(3.2.1)，幂法收敛性的分析较复杂，但若1=2=…=r且|1| > |r+1|…|n|则定理结论仍成立。
此时不同初始向量的迭代向量序列一般趋向于1的不同特征向量。
3.
求maxa(x)的流程，设数组x(n)数向量x的n个分量