空间迭代法

¾ 自适应多重逆迭代(adaptive multiple inverse iteration)[8]，是在迭代前对每个试向量 xi 进行对应的
多次逆迭代，从而加速高阶特征向量的收敛。
80
(k +1)
¾ 超松弛幂迭代(over-relaxed and matrix power accelerated subspace iteration)[9]，是用 X 和 X (k +1)
理想，而移频能很好的解决高阶重频特征向量的收敛速度慢的情况。因此，可以把超松弛幂迭代改进 方案用于移频的子空间迭代法。
ε |
Kϕ
(k i
)
−
λi(
k
)
M
ϕ
(k i
)
|2
<
|
Kϕ
(k i
)
|2Biblioteka ϕ(3)来控制收敛。 εϕ 可取 10-3 或 10-4，对于矩阵阶数比较大且求解模态数比较多时，最好取 10-4。
81
2.1 超松弛幂迭代 为了方便介绍超松弛幂迭代[9]改进方案，下面先给出矩阵幂迭代(matrix power accelerated
的增多抵消了一部分用移频所减少的工作量。
移频子空间迭代法中，子空间维数可取 q = max( s , 4) ，其中 s 为 L 中一行的平均非零元个数，
Imax 由第 II 与 III 步计算量之比确定[19]
I max
=
⎛ max ⎜
⎝
0.5Ns2 2Nsq + 3Nq2 + 2Nq
+ 10q3
⎞ ,4⎟
％。
文章分为五部分，第一部分为引言，简单介绍了近十年来子空间迭代法的进展；第二章介绍了
超松弛幂迭代和大移频；第三章把移频与超松弛幂迭代结合在一起，并从子空间维数选取等方面进行
了研究；第四部分为算例和讨论；最后一部分为结论。
2 子空间迭代算法
子空间迭代法是反幂法的推广。反幂法的矩阵迭代每次只计算一个特征值和特征向量，从最低 阶特征值开始，逐渐推进到高阶特征值。后来开始用多个向量同时进行迭代(simultaneous iteration)， 并引入了 Gram-Schmidt 过程来正交化向量，能同时求出几个特征值和特征向量。上世纪七十年代初， 在其中加入了子空间上的 Rayleigh-Ritz 过程[1]，明显改善了收敛速度。
1. 将 X (k ) 进行 M-正交归一化
2. 解试向量矩阵 (LDLT ) X (k+1) = MX (k ) ；
3. 计算 K 和 M 在 X (k+1) 上的投影
K (k +1)
=
T (k +1)
XK
X M (k +1) ， (k +1)
=
T (k +1)
XM
(k +1)
X
83
4.
求解 q 阶广义特征值问题
⎠
以下是移频子空间迭代算法的主要步骤： I．初始化 1. 确定子空间的维数 q
2. 选取初始向量矩阵 X 0
3. 设定每次移轴的最大迭代次数 Imax II．移轴与 Sturm 序列校核 1. 计算移轴 μ，应设法保证它不是特征值
2. 分解移轴刚度矩阵 K − μM＝LDLT
3. Sturm 序列校核 III．迭代 Imax 次，完成后转向 I
1 引言
在结构动力学、电磁学、声学、量子化学等计算领域中经常会求解广义特征值问题。广义特征值 问题作为科学计算和工程应用中的一个重要研究课题，已经有许多比较成熟的算法。在现今应用中经 常需要求解几十到几百个的特征值，矩阵规模也不断增加，因此需要寻找有效、稳定和快速的算法来 提高单机上的解题规模和速度。
广义特征值问题求解方法的选择取决于系统自由度大小、矩阵的稀疏性、所求解特征值的个数和 其在特征值谱上的位置。在结构动力学中需要求解
Kϕ − λMϕ = 0
(1)
的 p 个低阶特征值与特征向量( λi , ϕi )，i=1,2,…, p。其中，刚度矩阵 K 为稀疏的对称正定对称
（Symmetrical Positive Definite, SPD）矩阵，质量矩阵 M 一般为半正定的对角矩阵。但矩阵束正定， 即对任意的正数 μ，矩阵 K+μM 正定。在实际工程计算中，矩阵 K 和 M 的阶数可以从几十到几百万。
¾ 矩 阵 幂 迭 代 (matrix power accelerated subspace iteration)[7] ， 是 用 X (k+1) = (K −1M )2 X (k ) 代 替
X (k+1) = K −1MX (k ) ，也就是每进行两次逆迭代，才有一次 Rayleigh-Ritz 过程。
μ = λs+1 + λs+2
(5)
2
这时需要放松 Sturm 序列校核，即需要检查已收敛的特征值个数加 1 与因子矩阵 D 中负对角元的个 数是否相符。大移频移位方案比传统的移位方案平均快 17%[11]。
3 一种新的改进方案
超松弛幂迭代[9]改进方案的加速效果很有成效，缺点是在 λp 附近的高阶重频特征向量收敛速度不
当 K 和 M 规模比较大，例如矩阵阶数超过 1000 的大型特征值问题，子空间迭代法无疑是最受
青睐的方法之一。许多有限元软件, 例如 ABAQUS、ADINA、ANSYS 和 NASTRAN，早已把子空间 迭代法作为它们的广义特征值问题求解器。与迭代 Lanczos 法和迭代 Ritz 向量法相比，子空间迭代法 的速度慢一些，但稳定性要好很多，并且算法易于实现。
子空间迭代法自从上世纪 70 年代提出就得到广泛的应用，并在相当长的时间内成为有限元分析
的首选特征值问题解法。当时，需要求解的特征值个数较少(少于 20 个) ，但在现今应用中，抗震规
范要求 90％的质量参与系数经常导致求解几十甚至几百个振型。
当求解的模态数较多(大于 40)时，用传统上两次迭代之间的特征值的相对误差来控制收敛不能很 好控制高阶特征向量的收敛，需要用模态误差[10]
其中 λ (k +1) =
λ (k ) q
1 − δ (k +1)
82
1
δ (k +1) 的经验取值为 δ (k +1)
= 0.01δq(k )
=
0.01
⎡ ⎢
(λq(k
)
)
2
(k
xq
+1) T
(k +1)
pq
−
2λq(k )
(k +1)T
xq
p(k +1) q
⎤2 + 1⎥
⎢⎣
⎥⎦
具体的证明过程可参考[9]。
X = X Q (k +1)
(k +1) (k +1)
9)计算特征值的相对误差
tol
(k i
+1)
，并判断是否全部收敛。若全部收敛，则退出循环；否则，跳
到 2)继续进行迭代。 矩阵幂迭代加速收敛的效果明显，但稳定性有些差。在有些情况下，在一次循环内进行两次逆
迭代之后， X (k +1) 会变为几乎线形相关的，使投影后的刚度矩阵和质量矩阵变成奇异的，导致数值不
子空间迭代法假设 q 个线形无关的初始向量同时进行迭代，求得前 p 个特征值和特征向量。传 统上，q = min(2p, p+8)，用两次迭代之间的特征值的相对误差
( k +1)
tol i
=
λ λ − (k +1)
(k)
i
i
λ (k +1) i
< tol , i =1, 2, …, p
(2)
来控制收敛，tol 通常取 10-6。
6)计算 K 和 M 在 X (k ) 上的投影
K (k +1)
=
T (k +1)
XK
X M (k +1) ， (k +1)
=
T (k +1)
XM
(k +1)
X
7)求解投影后的 q 阶广义特征值问题
K Q = M Q Λ (k +1) (k +1)
(k +1) (k +1) (k +1)
8)形成新的近似特征向量
K Q = M Q Λ (k +1) (k +1)
(k +1) (k +1) (k +1)
5. 形成新的近似特征向量 X (k +1) = X Q (k +1) (k +1)
6. 按模态误差判断特征值和特征向量的收敛，移出已收敛的特征向量,并在 X 中加入随机向量 传统的移位方案[2,20]是相当保守的，它只用来加速数轴上其右侧模态的收敛速度，因为它左侧的 模态已经收敛。大移频(large shifting)[10]是将移位移到下两个待收敛特征值之间，即
subspace iteration)[7]的步骤：
1)选取初始向量矩阵 X (0)
2)计算 P(k +1) = MX (k ) , k =0,1,2,…
(k +1)
3)解出 X
(k +1)
KX
= P(k+1)
(k +1)
(k)
4) 计算 P = M X
5) 解出 X (k +1)
( k +1)
K X (k+1) = P
2.2 大移频
移频，又称谱变换(Spectrum Transformation)或移位，相当于在频率谱上对频率作了一维坐标变换， 从而加速待求的特征值和特征向量的收敛速度。移频[2]最早由 Bathe 在上世纪 70 年代提出，近十年来 随着求解特征模态个数的增加才广泛应用到特征值问题解法中。