数值计算方法

数值计算方法 第三章矩阵特征值与特征向量的计算

一些工程技术问题需要用数值方法求得矩阵的全部或部分特征值及相关的特征向量。一些工程技术问题需要用数值方法求得矩阵的全部或部分特征值及相关的特征向量。

3.1 特征值的估计 • 较粗估计(A)  ||A|| • 欲将复平面上的特征值一个个用圆盘围起来。 • 3.1.1 盖氏圆 • 定义3.1-1 设A = [aij]nn，称由不等式 • 所确定的复区域为A的第i个盖氏圆，记为Gi： • i = 1，2，…，n。 • 定理3.1-1 若为A的特征值，则

定理3.1-1 若为A的特征值，则 • 证明：设Ax = x (x 0)，若k使得 • 因为 •  •  • 

G4 G1 G2 G3 • 例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 盖氏圆的连通部分 • 称相交盖氏圆之并构成的连通部分为连通部分。 • 孤立的盖氏圆本身也为一个连通部分。 • 定理3.1-2 若由A的k个盖氏圆组成的连通部分，含且仅含A的k个特征值。

定理3.1-2 若由A的k个盖氏圆组成的连通部分，含且仅含A的k个特征值。 • 证明: 令D = diag(a11，a22，…，ann)，M = A–D，记 • 则显然有A(1) = A，A(0) = D， • 易知A()的特征多项式的系数是的多项式， • 从而A()的特征值1()，2()，…，n()为的连续函数。

i aii • A()的盖氏圆为： • 因为A(0) = D的n个特征值a11，a22，…，ann，恰为A的盖氏圆圆心，当由0增大到1时，i()画出一条以i(0) = aii为始点，i(1) = i为终点的连续曲线，且始终不会越过Gi；

不失一般性，设A开头的k个圆盘是连通的，其并集为S，它与后n–k个圆盘严格分离，显然，A()的前k个盖氏圆盘与后n–k个圆盘严格分离。不失一般性，设A开头的k个圆盘是连通的，其并集为S，它与后n–k个圆盘严格分离，显然，A()的前k个盖氏圆盘与后n–k个圆盘严格分离。 • 当= 0时，A(0) = D的前k个特征值刚好落在前k个圆盘G1，…，Gk中，而另n–k个特征值则在区域S之外，从0变到1时，与始终分离（严格）。连续曲线始终在S中，所以S中有且仅有A的k个特征值。

注：1) 每个孤立圆中恰有一个特征值。 • 2) 例1中G2，G4为仅由一个盖氏圆构成的连通部分，故它们各有一个特征值，而G1，G3构成的连通部分应含有两个特征值。 • 3) 因为例1中A为实方阵，所以若为A的特征值，则也是A的特征值，所以G2，G4中各有一个实特征值。

3.1.3 盖氏圆与相似变换 • 由于特征值是相似不变量，所以代数上常用相似变换将矩阵化简以得到特征向量，这里也可用相似变换将盖氏圆的半径变小，以得到更好的估计。 • 原理：取对角阵作相似变换阵：P = diag(b1，b2，…，bn)其中bi > 0，i = 1，2，…，n • 则 • 与A有相同特征值. • 而B的第i个盖氏圆为：，

而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}； • G3 = {z：| z– 0.4| 0.03} • 3G3，较好。

为了更好地估计另外两个特征值可取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} • 三个盖氏圆分离，故有1G1'，2G2'，3G3。

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) • 若a1 0，则因知， • 当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 • 都将停机。须采用“规范化”的方法 • ， k = 0，1，2，… (3.2-4)

规范化： • ， k = 0，1，2，… (3.2-4) • 定理3.2-1 任给初始向量x(0)  0 ，有 • (3.2-5)

证明：

而 • 注：若A的特征值不满足条件(3.2-1)，幂法收敛性的分析较复杂，但若1 = 2 = … =  r，且|1| > | r +1| … |n|，则定理结论仍成立。 • 此时不同初始向量的迭代向量序列一般趋向于1的不同特征向量。

3. 算法 • 求maxa(x)的流程，设数组x(n)存放向量x的n个分量 • 幂法流程：

Matlab幂法流程：

例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]; • 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|，考虑收敛加速

1. 特征值的Aitken加速法 • (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*)，其中 • 令 • 则 • 步骤： • 计算

流程图

Matlab流程图

例2 用幂法求方阵A的最大模特征值，并用Aitkem加速法 • 解：见（P314）

幂法 • A=[-4,14,0;-5,13,0;-1,0,2]; • x0=[1;1;1];k=1 • y=x0/maxa(x0) • x1=A*y • while(abs(maxa(x1)-maxa(x0)))>0.01 • x0=x1;k=k+1 • maxa(x0) • y=x0/maxa(x0) • x1=A*y • end;

Aitkem加速 • A=[-4,14,0;-5,13,0;-1,0,2];l1=0;k=1 • 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)) • while (abs(l1-l0))>0.01 • x0=x1;x1=x2;l1=l0;k=k+1 • 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–a为A–aI的最大模特征值，且 • （k–a是A–aI的次最大模特征值） • 则对A–aI计算1–a及对应的特征向量比对A计算收敛得快，此即为原点平移法。 • 计算1–a及特征向量的迭代公式 • 特征向量： • 特征值：max(x(k)) 1–a，a + max(x(k)) 1。

注：a的选取较为困难。 • 例3 设，， • 求最大模特征值及特征向量。 • 解：（P315）

幂法： • A=[-3,1,0;1,-3,-3;0,-3,4]; • x0=[0;0;1];k=1 • y=x0/maxa(x0) • maxa(x0) • x1=A*y • while(abs(maxa(x1)-maxa(x0)))>0.01 • x0=x1; k=k+1 • y=x0/maxa(x0) • maxa(x0) • x1=A*y • end

原点平移法： • A=[-3,1,0;1,-3,-3;0,-3,4]; • x0=[0;0;1];k=1 • y=x0/maxa(x0) • maxa(x0) • x1=(A+4*eye(3))*y • while(abs(maxa(x1)-maxa(x0)))>0.01 • x0=x1; k=k+1 • y=x0/maxa(x0) • maxa(x0) • x1=(A+4*eye(3))*y • end; • y • maxa(x1)-4

3. 对称矩阵的Rayleigh商加速法 • 定义设A对称，x 0，则称 • 为x关于A的Rayleigh商 • 思想：A对称，特征值1，2，…，n均为实数，且存在特征向量v1，v2，…，vn为标准正交基。 • 设，a1 0，则

当k充分大时，M'与k无关)

注: 此比Aitken加速中的(3.2-6)更快 • 公式 • 称为Rayleigh商加速法。 • 其中

注：有了R(x(k))，R(x(k+1))，R(x(k+2))，的值，可再用Aitken加速法得到的一个更好的近似值：注：有了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) • maxa(x0) • x1=A*y • while(abs(maxa(x1)-maxa(x0)))>0.001 • x0=x1; k=k+1 • y=x0/maxa(x0) • maxa(x0) • x1=A*y • end;

Rayleigh商加速法： • A=[6,2,1;2,3,1;1,1,1]; • x0=[1;1;1];r=0;k=1 • y=x0/maxa(x0) • maxa(x0) • x1=A*y • while(abs(r1-r))>0.001 • x0=x1;r1=r; k=k+1 • y=x0/maxa(x0) • maxa(x0) • x1=A*y • r = y'*x1/(y'*y) • end

3.2.3 反幂法 • ——用A–1代替A作幂法，即反幂法 • 1. 求最小模特征值及相应的特征向量 • 若A可逆，|1| > |2| … |n|为其特征值，则为A-1的最大模特征值。 • 迭代公式：x(k+1) = A–1x(k)，k = 0，1，2，…， • 但A–1不易求，通常可解方程组Ax(k+1) = x(k)来求x(k+1) • 即有 • (3.2-12)

(3.2-12) • 当k时有 • 注：为解(3.2-12)中的方程组。对作LR分解（带行交换）PA = LR则有

2. 求任一特征值及相应特征向量 • ——反幂法结合原点平移法 • 思想：若已知为j的近似值，则的特征值是 • 而显然非常大（最大），比值很小 • 迭代公式：

迭代公式： • 当k时有 • 注：(1) 若有LR分解，则迭代公式 • (3.2-16)

(2) 在(3.2-16)中直接取 z(1) = (1，…，1)T作初值开始迭代称为半次迭代法

数值计算方法

数值计算方法

Presentation Transcript

数值计算方法