第六章计算方法

λ1 = 45 , λ 2 = 2, λ 3 = 1
x1 ≈
x1 max( x1 )
uk
1 -0.67153 -0.66727 -0.66670 -0.66667 -0.66667 -0.66667 -0.66667
x1 = (3,1, − 2 ) T , x 2 = (3, 2 , − 3) T , x3 = ( 2 ,1, − 2 ) T
Di = z − aii ≤ Λi , i = 1, 2, , n
Di = z − aii ≤ Λi , i = 1, 2,
,n
P128
例 1用Gerschgorin定理估计矩阵 的特征值的范围。 解
⎡4 1 1 ⎤ ⎥ A=⎢ ⎢1 10 1 ⎥ ⎢ ⎥ 1 1 5 ⎣ ⎦
C 3 = { z : z − 5 < 2}
a11 − λ det( A − λE ) = a21 an1
a12 a22 − λ an2
a1n a2n ann − λ
建立特征方程：
det ( A − λ E ) = 0
求解特征方程，所得根λ0 即为矩阵A的特征值，然 后求解方程组﹙Ａ﹣λ0E﹚Ｘ﹦0，就可得出矩阵A对应 于特征值λ0的特征向量X。
A的三个特征值与特征向量分别是：
解:
k 0 1 2 3 4 5 6 7
取初始向量V0= u0=(1,1,1)T ，计算出Vk,uk和mk,迭代 7次的结果列于下表
Vk
1 274 44.43277 44.92333 44.99572 44.99959 44.99953 44.99953 1 95 14.84322 14.97623 14.99865 14.99988 14.99983 14.99983 1 -184 -29.64262 -29.95048 -29.99722 -29.99974 -29.99968 -29.99968 1 1 1 1 1 1 1 1 1 0.34672 0.33413 0.33337 0.33334 0.33333 0.33333 0.33333
2
2.反幂法（求矩阵的按模最小特征值） 基本思路： 设A没有零特征值，则A非奇异，即A的逆阵存 在。设A的特征值为 λ ≥ λ ≥ > λ > 0
1 2 n
反幂法规范化算法 (1) 任取一非零向量u0=V0∈Rn,一般可取V0=(1,1,.…,1)T (2) 计算Vk=A-1 uk-1 求解AVk=uk-1 (3) mk =max(Vk ), uk= Vk / mk 当k足够大时, 即可得到： λn≈1/ mk 有时为了加速反幂法，将(2)中
2
=1
平面π 方程 W T x = 0 ∀x ∈ π
T
Householder变换阵的一个重要性质是可在一个 向量中引入零分量. 定理6-5 对任意非零向量, x=(x1,x2,…,xn)T,存在一个 Householder变换阵 H,使 Hx=σe1 ⎛ x1 ⎞ ⎛1⎞ 其中,e1=(1,0,…,0)T, σ=±‖x ‖2 ⎜ ⎟ ⎜ ⎟ x 0 证明 H ⎜ 2 ⎟ =σ ⎜ ⎟ 令 U=x-σe1,W=U/‖U‖2 ⎜ ⎟ ⎜ ⎟ 构造Householder变换阵 ⎜ ⎟ ⎜ ⎟
Vk = A−1uk −1
改为
Vk = ( A − qE )−1 uk −1
这时求得特征值(λ-q)-1 ，这里λ为所有λi中与q最接近者。
注意： −1 由于求逆非常费时。故在求Vk时， Vk = A uk −1 可采用解方程组 AV k = u k − 1 的办法。 由于每次解方程组的系数矩阵都相同，可预先作三角分解， 这样每次迭代仅仅求解两个三角方程组就可以了。特别当 n 较大时，将大大地节省计算量。 反幂法的适用范围是求矩阵的按模最小特征值及对应的特征向 量。
但众所周知，高次代数方程求根是相当困难的，而且重 根的计算精度较低。同时，矩阵A求特征多项式系数的过程对 舍入误差十分敏感，这对最后计算结果影响很大。因此，从 数值计算角度来看，上述方法缺乏实用价值。 目前，求矩阵特征值问题实际采用的是迭代法。这里 将介绍两种方法：幂法、反幂法以及QR方法。 引入Gerishgorin圆盘
(V k ) j
α 1 ( x1 ) j + [ ∑ α i (
i=2 n
mk = Vk 幂法(计算矩阵A的按模最大特征值) (1)任取一非零向量V0∈Rn,一般可取V0=(1,1,.…,1)T (2)计算Vk=AVk-1 (3)当k足够大时,即可得到： λ 1 ≈ ( V k ) j
(V k − 1 )
i=2
n
i
n
(
λi k ) xi ] λ1
i
Αxi = λi xi
同理可得 假定
( x1 ) j ≠ 0
V k −1 = λ
,因为
k −1 1
[α 1 x1 +
∑α
i=2
(
λ i k −1 ) xi ] λ1
, n ) ,故得
2 2 2 V2 = AV +αnλ2 n xn 1 = A V0 = α1λ 1 x1 +α2λ2 x2 +
ρ U 2 1 1 2 ⎛ x1 − σ ⎞ 2 2 2 其中 ρ = U 2 = (( x1 − σ ) + x 2 + ... + x n ) ⎜ ⎟ 2 2 x2 ⎟ U = x − σ e1 = ⎜ 1 ⎜ ⎟ = (2σ 2 − 2 x1σ ) = σ (σ − x1 ) ⎜ ⎟ 2 x n ⎝ ⎠ 取 σ = − sign ( x1 ) x 2
H = HT
H = E − 2 WW
T
, W
2
=1
HH T = H 2 = ( E − 2WW T )( E − 2WW T ) = E − 4WW T + 4WW T WW T = E .
H −1 = H T = H
⎡ 1 ⎤ ⎢ 2⎥ 1 ⎞ ⎢ ⎥⎛ 1 0 ⎜ H = E − 2WW T = E − 2 ⎢ 0 ⎥ ⎜ ⎟ ⎟ 2 2 ⎠ 1 ⎢ ⎥⎝ ⎢ 2⎥ ⎣ ⎦ 0 − 1⎤ ⎡ 0 1 0 ⎥ = ⎢ H为Householder变换阵 ⎢ 0 ⎥ ⎢ 1 0 0 ⎥ − ⎣ ⎦
当
λ λ2 << 1 时 幂法收敛快，当 2 ≈ 1 时 幂法收敛慢。 λ1 λ1
m 2 = 44 . 42377 , m 3 = 44 . 92333 , m 4 = 44 . 99572 m 5 = 44 .99959 , m 6 = 44 . 99953 , m 7 = 44 .99953
小结： 幂法适用求矩阵的按模最大特征值及相应的特征向量， 其优点是算法简单，容易编写程序在计算机上实现，缺点是收 敛速度慢，其有效性依赖与矩阵特征值的分布情况. 为了加速 幂法，人们提出若干措施: （1）原点平移法，（2）Aitken加 速法，Rayleigh商加速，见 P155。
6.1 引言
工程技术的许多实际问题，例如振动问题，稳定问题的求 解，有时会归结为求矩阵的特征值λ和对应的特征向量χ。 学过线性代数后，我们已知求矩阵A的特征值λ和特征向量χ 的解法，即先求出A的特征多项式：
6.1 引言 6.2 幂法与反幂法 6.3 矩阵的正交分解 6.4 QR方法 6.5 雅可比方法
x1 max( x 1 )
lim m k = λ 1
k→∞
幂法规范化算法 (1) 任取一非零向量u0=V0∈Rn,一般可取V0=(1,1,.…,1)T (2) 计算Vk=Auk-1 (3) mk =max(Vk ), uk= Vk / mk 当k足够大时, 即可得到： λ1≈ mk
⎛ 1 ⎞ ⎛−1/2⎞ ⎜ ⎟ ⎜ ⎟ V V 1 = ⎜ − 2 ⎟ , m 1 = max( V 1 ) = − 2 , u 1 = 1 = ⎜ 1 ⎟ m 1 ⎜ 0 ⎟ ⎜ ⎟ 0 ⎝ ⎠ ⎝ ⎠
j
∞
定理6-3：在定理1 的条件下,规范化向量序列{uk}收敛于 矩阵 A 按模最大的特征值 λ1 对应的特征向量 , 而向量序列 {Vk}的绝对值最大的分量mk收敛于λ1,即
k→∞
lim u k =
按上述计算过程，有一严重缺点，当|λ1|>>1 （或| λ1 |<<1 时） {Vk} 中不为零的分量将随 k 的增大而无限增大，计算机就 可能出现上溢（或随 k 的增大而很快出现下溢），因此，在实 际计算时，须按规范法计算，每步先对向量Vk进行“规范化”， 即取Vk中绝对值最大的一个分量记作mk =max(Vk )，用mk遍除 的向量Vk 的所有分量，得到规范化向量uk (uk = Vk / mk)。
uk ≈
x1 max( x1 )
例1 用幂法求矩阵
⎡ 133 6 135 ⎤ A = ⎢ 44 5 46 ⎥ ⎥ ⎢ ⎥ ⎢ − 88 − 6 − 90⎦ ⎣
由上可见经过7次迭代, m7的值已稳定到小数后5位，故所 求的按模最大特征值和对应的特征向量可取作：
按模最大特征值λ1和对应的特征向量x1
λ1 ≈ 44 .9995 , x1 ≈ (1,0 .333 , − 0 .6667 ) T
2
平面π称为。
若 x ∈ π , Hx = ( E − 2WW ) x = x − 2WW T x = x
（3）若 Hx = y, 则 x 2 = y
(4)镜面反射性质
2
H 变换下，向量的长度保持不变。
（i）Hw = − w −1为H的一个特征值，w为对应的特征向量;
（ii）H的其余n-1个特征值为1。
对应于1的特征向量是与w正交的非零向量。
3
（4）镜映射 − 几何意义
H = E − 2 WW
T
, W
, x n (它们线性无
关)，则对任意一个非零向量V0∈Rn 所构造的向量序 列 Vk = AV k−1 ,有