三种常用固有振动特征值解法的比较

2005全国结构动力学学术研讨会

海南省海口市，2005.12.19-20

中国振动工程学会结构动力学专业委员会

三种常用固有振动特征值解法的比较

宫玉才1

周洪伟 陈 璞 袁明武

（北京大学力学与工程科学系 北京，100871）

Email ：yuanmw@

摘要： 本文以高效的细胞稀疏直接快速解法为核心步骤，实现了快速的固有振动广义特征值问题解法, 并在相同的允许模态误差的意义下检验了三种结构动力学中常用的大型矩阵特征模态算法——子空间迭代法、迭代Ritz 向量法和迭代Lanczos 法的计算效率。迭代Ritz 向量法平均而言最快，子空间迭代法最慢，三种解法效率相差不是太大。与ANSYS 的子空间迭代和Lanczos 法相比，本文的子空间迭代比ANSYS 的效率高很多，Lanczos 法和ANSYS 的差不多 。大量较大规模的例题显示，本文对特征值算法的改进是十分有效的，算法的健壮性，通用性都达到了高水平。

关键词：特征值，结构振动，迭代法，高效能计算

1

高等学校博士学科点专项科研基金资助项目 （编号：20030001112）

引言

在工程有限元分析中常常要求解广义代数特征值问题

0K M ϕλϕ−= (1)的部分低阶特征值与特征向量。对于矩阵阶数超过1000的大型问题，子空间迭代法、Ritz 向量法和Lanczos 法被公认为求解部分低阶极端特征值和特征向量的有效方法。尽管国内外的有限元软件都提供广义代数特征值问题(1)的多种解法，但结果仍然不能令人完全满意，漏根与多根、自由模态误判都时有发生。

传统上，低端特征值问题求解过程极度依赖于谱变换的线性方程组

()T K M x LDL x My µ−==

(2)

的解法，移轴矩阵K M µ−的LDLT 三角分解是计算量最大的主要步骤。在以变带宽解法为核心步骤的特征值解法中，它常常占到特征值问题计算时间的70%到90%。本文采用了文[1]提出的一个效率非常高的有限元解法-细胞稀疏直接快速解法（简称细胞解法）替换变带宽解法，极大地提高了三角分解的效率。

如果要求不太多的特征模态，例如10个，通常认为Ritz 向量法和Lanczos 法具有比子空间迭代法更高的计算效率，Ritz 向量法和Lanczos 法比子空间迭代法平均快4～10倍[2]。但是，标准的Ritz 向量法和Lanczos 方法对收敛的判定是相对含糊的，在实用的工程计算中可能造成漏根或多根。

传统上，子空间迭代用特征值的两次迭代之相对误差不等式

()(1)(1)

||

||

l l k k l k λλλελ++−≤ (3)

控制收敛，而Lanczos 法用其过程中的不等式 1||||k i q qi s νλνβε−−≤<

(4)

控制收敛。在大量的工程计算中，发现在允许误差510λνεε−==的情形下，除最低的十几阶模态之外，子空间迭代与Lanczos 法所得到的特征向量精度都可能不令人满意。这一现象对Lanczos 方法尤为严重，原因是采用逆迭代技术时，高阶的、密集的特征值不易分离。

关于特征模态的收敛，不同的算法往往采用不同的标准，相对速度的比较不是很客观。对于特

征模态的近似(,)k

k λϕ%%，在各种算法中可以统一用模态误差(5) 代替特征值误差作为收敛判据，来衡量算法的效率。

||||||||k k k k k

k k k k K M K M K M ϕϕλϕϕλϕεϕλϕ

−−≈≤%%%%%%%%%

(5)

模态误差有明显的物理意义，k K ϕ是振型k ϕ的最大弹性节点力，而k

k M λϕ%%是振型k ϕ的最大惯性节点力。式(5)的最左端是非平衡节点力与最大弹性节点力之比，中间是非平衡节点力与最大惯性力之比

[3]。

大量算例表明，在模态误差的意义之下，收敛过程平稳，三种解法效率相差不是太大。迭代Ritz 向量法平均最快，子空间迭代法最慢。

本文最后还与 ANSYS V.8.1的子空间迭代法和块Lanczos 法进行了比较。

1 算法描述

1.1 子空间迭代法 最初是由Clint 和Jennings 提出，是反幂法的推广[4]。稍后，Bathe 和Wilson 在其中加入了子空间上的Rayleigh-Ritz 过程。它可以明显地改善收敛速度[4,5]。以下是一个子空间迭代算法的主要步骤。

I ．初始化

1. 确定子空间的维数q

2. 选取初始向量矩阵 N q X R ×∈

3. 设定每次移轴的最大迭代次数max I II ．移轴与Sturm 序列校核

1. 计算移轴µ，应设法保证它不是特征值

2. 分解移轴刚度矩阵T K M LDL µ−＝

3. Sturm 序列校核

III ．迭代max I 次，完成后转向步骤I 1. 将X 进行M-正交归一化 2. 向量矩阵 11()T X LDL MX −=；

3. 计算K 和M 在X 1上的投影，*11T K X KX =，*11T M X MX =

4. 求解q 阶广义特征值问题 ***K M ∗∗Φ=ΦΛ

5. 形成新的近似特征向量 1X Y ∗=Φ

6. 按模态误差判断特征值和特征向量的收敛，移出已收敛的特征向量，并在X 中加入随机向

量或减缩子空间的大小。

子空间迭代法假设q 个初始向量同时进行迭代，求得前p 个特征值和特征向量。传统上，

min(2,8)q p p =+，但当模态数目需求较多时，这种取法显然是不现实的。经验表明，子空间维数

可取4)q =，其中s 为L 中一行的平均非零元个数，由第II 与III 步计算量之比确定。

1.2 迭代Ritz 向量法 Ritz 向量法由Wilson ，Yuan （袁明武）和Dickens 在1982年提出的[6], 也称为

WYD-Ritz 向量法，最初用来求解地震的动力响应问题。后来，袁明武等将其用于大型特征值问题的计算，使它成为一种极为有效的特征值算法[7]。引入迭代可以改善特征值与特征向量的精度，具体步骤如下：

I ．初始化

1. 确定块Ritz 向量法块宽q 与生成步数r

2. 选取初始向量矩阵0N q Q R ×∈

3. 设定每次移轴的最大迭代次数max I II ．移轴

1. 计算移轴µ，应设法保证它不是特征值

2. 分解移轴刚度矩阵T K M LDL µ−=

3. Sturm 序列校核

III ．迭代max I 次，完成后转向II

1. 对0,1,,k r =L 解 1T k k LDL Q MQ +=，然后用将 1

k Q +对已收敛的特征向量以及12,,,k Q Q Q L 作M-正交归一化，并形成1k Q +。

2. 计算K 在12(,,,)r Q Q Q Q =L 上的投影，*T K Q KQ =

3. 求解q ×r 阶标准特征值问题**

K ∗∗Φ=ΦΛ 4. 形成新的近似特征向量X Q ∗=Φ

5. 按模态误差判断特征值和特征向量的收敛，移出已收敛的特征向量

6. 如果达到了预期的特征值个数，退出；否则将未收敛的前q 个近似向量作为初始向量进行下一次迭代

1.3 迭代Lanczos 方法 Lanczos 方法是在五十年代初提出的，它用正交向量组约化对称矩阵为三对角矩阵。七十年代以前它被认为不稳定，用于实际计算不多。1972年Paige 证明了失去了正交性的充分必要条件是其投影矩阵的特征值收敛到原矩阵的特征值。此后，Wilkinson 等建议了重正交化方案，Golub, Cullum 和Donath, Underwood 等建议了块Lanczos 方法，Underwood 建议了迭代Lanczos 方法[8,9]。