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

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

海南省海口市，2005.12.19-20

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

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

宫玉才1

周洪伟 陈 璞 袁明武

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

Email ：yuanmw@https://www.360docs.net/doc/3d10942404.html,

摘要： 本文以高效的细胞稀疏直接快速解法为核心步骤，实现了快速的固有振动广义特征值问题解法, 并在相同的允许模态误差的意义下检验了三种结构动力学中常用的大型矩阵特征模态算法——子空间迭代法、迭代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]。

本文采用一个带重正交的迭代块Lanczos 方法，向量生成步骤与上面的迭代Ritz 向量法一致，其差别是Lanczos 方法利用了Ritz 向量生成过程中的正交归一化系数。迭代Lanczos 方法的第I ，II 步与迭代Ritz 向量法完全一致，为了节省篇幅，我们仅给出第III 大步中的第1，2步：

III ．迭代 I max 次，完成后转向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 +。在此过程中依次形成*T

2. 求解q ×r 阶标准特征值问题**T ??Φ=ΦΛ。

文献上一般都将Krylov 空间12(,,,)r Span Q Q Q L 的维数取得较大，例如待求特征值个数的2倍，以期一次完整的Lanczos 过程得到全部想要的模态。本文是在较小的Krylov 空间上完成Lanczos 过程，然后用最好的q 个近似特征向量作为下一次Lanczos 过程的初始向量。这样的方案是一个子空间迭代与Krylov 空间结合的算法，具有与子空间迭代法同样的可靠性，在文献上还较为少见。

1.4 程序实现 在程序实现中，移轴三角分解，向前消元和向后回代采用了细胞解法[1]，它的综合效率是变带宽解法的数倍至数百倍。在其它方面，循环展开也广泛地应用于各种计算加速中，例如正交归一化等。算法中需要多次计算的乘积KX 采用稀疏总体刚度矩阵与向量乘积的方案计算，它的计算量与计算时间比求解方程(2)小一个量级以上，因而不影响整体的计算效率。

如果在正交化过程中，一个向量正交化以前的模与正交化以后的模之比超过了某一阈值，我们将对此向量实施双正交化，即对已正交化的向量再实施正交化。在下面的数值试验中，这一阈值取为1210。

为节省计算量，Ritz 向量法与Lanczos 法的III.4在形成新的特征向量时，未计算全部Rayleigh-Ritz 特征值相对应的近似特征向量。

在子空间迭代法中，投影特征值问题的解法选用了广义Jacobi 方法[3]，允许误差取为242?。在Lanczos 方法和Ritz 向量法中，求解投影特征值问题采用了Householder 变换与QL 方法的组合[10]。 根据经验，子空间迭代法中与一次移轴三角分解相应的最大迭代次数取为

223max max(0.5/(3210),6)I Ns Nq Nq q =++ (6)

迭代Ritz 向量法与迭代Lanczos 法块的大小一般取为4q =，步数为6r =；这相当于Krylov 子空间的维数是24。类比于子空间迭代法，与一次移轴三角分解相应的最大迭代次数取为

22233max max(0.5/(3210),4)I Ns Nq r Nqr q r =++ (7)

三个算法的Sturm 校核均取后验方式，移轴的首选为下两个待收敛特征值之中点，即

12

0.5()k k μλλ++=+%% (8)若失败，则选为

10.98()k k k μλλλ+=+?% (9)其中k

λ是最后一个已收敛的特征值，1k λ+%，2

k λ+%是下两个待收敛的特征值。

2 算例与讨论

大量的实际工程问题被用来检验本文的三种方法，限于篇幅仅在表1和表2中列出其中的一小部分，其中PKUSTK 系列在以前的研究中已多次使用[1,11]。必须强调，基于半带宽解法或变带宽解法的特征值算法来求解这些阶数的问题时，在时间和空间上都是十分困难的，特别是求解数十个特征模态时。

算例与讨论的第一部分是三个解法自身的比较，第二部分是本文的方法与Windows NT ANSYS

8.1的速度比较。测试平台是Window 2000或Windows XP 的机器，CPU 均为Pentium IV 系列，编译器

选用的是Compaq Visual Fortran 6.5，

编译中未对特定处理器优化。测试程序接受Harwell-Boeing RSA 矩阵交换格式[12]。

2.1 三种广义特征值解法的比较 这一部分的数值试验都在一台带IDE 磁盘，操作系统为中文Windows2000的Pentium IV 2.0机器上进行的，内存为512 MB 。目标是在相同的精度要求下，比较三种解法在模态误差意义下的求解效率，即计算量与求解时间。三种方法的数据区大小均控制在384 MB ，所有计算可在内存中完成。表1给出了例题的工程背景以及大小，表2则是计算时间、移轴次数和迭代次数。这部分试验例题是由SAP84生成的，并间接转换成Harwell-Boeing RSA 格式[13]。

表中neq 是方程组的大小，而K 与L 分别是刚度矩阵与它的三角因子的存储量。表中的八个算例中，PKUSTK03、PKUSTK13 、MM-08特征值分布密集。

计算过程中，除PKUSTK03的允许模态误差选为4

10?ε?=之外，其它各例题的允许模态误差均取

为3

10?ε?=。三种算法求得的固有振动特征值的前6位有效数字完全一致。各算例的各阶特征值最大相对误差为52.5610?×，平均为61.5610?×。这间接说明了三种迭代方法的可靠性。平均意义下，迭代Ritz 向量法最快，它的计算速度是子空间迭代法的2倍左右。迭代Lanczos 法比迭代Ritz 向量法稍慢一点，但仍然比子空间迭代法快。

算例PKUSTK03是在一大方钢仓内有24个互相连接的筒仓，它的固有振动特征值十分密集。在70.49到130.82的区间内分布了79个特征值，两个相邻特征值之间的相对误差有时在410?的量级。如

采用式(3)的特征值相对误差控制子空间迭代的收敛，在取8

10λε?=时，仍不能很好地控制收敛。用310?ε?=在计算中也出现了收敛到不正确特征值与特征向量的例子，3种方法的同阶特征值计算结果

的最大相对误差达31.4710?×。在取4

10?ε?=后，3个方法的特征值结果的最后相对误差降为

61.1210?×。一般来说，取更小的?ε需更长的求解时间。测试表明310?ε?=与4

10?ε?=的求解时间的

相对差为10%到17%，但也有例外，MM-08的子空间迭代法用4

10?

ε?=作为收敛判据的求解时间为885.43秒，而3

10?ε?=时的求解时间为1098.28秒，这可能是由于低阶模态的精度提高所引起的。

表1 试验例题以及其说明

矩阵大小(MB)算例 说明 neq |K| |L| PKUSTK01植物园温室22,044 3.86 15.97PKUSTK03群仓 63,33612.37 65.94PKUSTK11围堰 87,80420.54 214.53PKUSTK12吉建大厦 94,65330.07 127.06PKUSTK13机械零件 94,89326.09 234.21MM-08 体育场 160,07425.68 175.56GUIZHEN 世纪坛晷针6,840 1.48 6.82PALACE

少年宫

10,860

2.28 12.82

表2 计算10个与80个特征模态所需的CPU 时间（单位：秒）与计算次数

子空间 Ritz Lanczos

算例 特征模态个数

CPU 时间 移轴次数迭代次数CPU 时间 移轴 次数 迭代次数CPU

时间 移轴次数迭代次数10 6.41 2 22 4.27 1 4 3.89 1 4 PKUSTK01 80 70.83 9 16935.787 21 39.037 25 10 53.53 3 41 28.95

2 7 28.88 2 8 PKUSTK0

3 80 330.84

11 201

154.42

8 26 133.859 30 10 51.14 1 6 48.52 1 3 47.67 1 3 PKUSTK11 80 945.62

8 145

399.32

7 26 330.17 6 21 10 39.70 1 14 32.66 1 4 40.92 2 5 PKUSTK12

80 651.3511 201

237.76

6 23 257.73

7 2

8 10 73.20 1 11 50.91

1 3 50.86 1 4 PKUSTK13 80 504.18 4 75 375.18 6 21 416.457 25 10 129.64

2 30 78.64 2 6 79.64 2 6 MM-08 80 1098.28

10 185

418.78

7 27 523.98

10 35

10 1.30 1 11 1.05 1 3 0.95 1 3 GUIZHEN 80 19.90 7 129

7.25 5 20 7.12 6 22

10 3.59 2 15 2.70 1 4 3.61 2 5 PALACE

80

39.01 7 126

19.93

7 24 19.84

7 25

子空间迭代法的实际测试中发现，310?ε?=大致相当于7

10λε?=或更小，410?ε?=则相当于

910λε?=。另外，本文的Lanczos 方法与WYD-Ritz 向量法采用的方案基本上是一致的，所以它们的计算时间大致相同。

2.2 与ANSYS 的比较 为了说明本文方法与国际先进水平的相对速度，与Window NT ANSYS 8.1的子空间迭代法和块Lanczos 方法进行了对比。表3所列的有限元模型是在ANSYS 中形成的，并通过DUMP 命令导出Harwell-Boeing RSA 格式的总体刚度矩阵与聚集总体质量矩阵。

所有的比较是在Pentium IV 2.4，512MB RAM 的机器上进行的，操作系统为Windows XP 。在测

试中，本文的三种算法仍采用3

10?ε?=作为收敛判据，

ANSYS 的子空间迭代法的收敛判据应为式(3)，510λε?=。ANSYS 与NASTRAN 的块Lanczos 法均采用了文献[14]上的一个精巧方法。表3与表4分别给出了试验例题以及计算所用的CPU 时间。

表3 ANSYS 试验例题以及它们的说明

算例

说明

neq |K|(MB)

TABT 高层建筑，8800个SHELL63单元， 1000个BEAM188单元60,108 11.05 BUAA1直升机齿轮，2840个

SOLID95单元 46,080 24.91 BUAA2

直升机齿轮，5360个

SOLID95单元

82,800 55.26

表4 ANSYS 与本文特征特征值计算方法的CPU 时间比较（秒）

ANSYS

本文 算例

模态数

子空间 Lanczos 子空间 Lanczos Ritz

10 68.11 15.87 18.16 16.30 12.66 20 162.53 23.20 68.83 37.91 25.77 40 254.92 29.44 146.23 49.19 41.66 TABT

80 533.22 53.45 289.31 93.81 84.03 10 822.16 72.42 42.14 23.92 24.16 20 1738.2092.67 69.91 47.73 34.19 40 4752.39132.28 166.83 84.00 87.41 BUAA1

80 5496.28221.50 378.72 185.28 169.36

10 405.44 282.78 83.09 53.47 55.78 20 >9000 345.89 174.27 95.55 102.33

40 484.58 317.70 220.81 222.25

BUAA2

80 649.20 586.47 448.75 481.95

算例TABT 的特征值分布不密集；算例BUAA1为1刚体自由度结构，特征值分布较密集；算例

BUAA2为6刚体自由度结构，重根，特征值分布较密集。 除BUAA2之外，ANSYS 与本文几种算法的周期结果在6位有效数字的范围内完全一致。BUAA2为一具有6个刚体自由度的轴对称结构，它有多对二重根。计算中，ANSYS 的Lanczos 算法未能正确识别全部6个刚体自由度。ANSYS 的子空间迭代则不能发现全部重根；在计算10个特征值时遗漏了4个重根从而导致迭代失败，计算20个特征值时出现不收敛的现象。

与ANSYS 的子空间迭代比较，本文的子空间迭代法不仅具有明显的速度优势，而且在计算控制上也更胜一筹。与ANSYS 的块Lanczos 法相比，本文的迭代Ritz 向量法与迭代Lanczos 法是相当保守的和安全的，算法上还不够精巧。本文的方案的迭代次数与三角分解次数明显多于ANSYS 的，但由于采用的细胞解法的三角分解速度以及消元回代速度比ANSYS 的相应算法快一些，才使得本文的方案在整体上略有优势。在计算少量特征模态时，三角分解的计算时间起控制作用，本文的方法比ANSYS 快，但求较大数量的特征值时，迭代占去较大部分的计算时间，与ANSYS 之间的差距缩小。

ANSYS 的块Lanczos 方法采用了多Lanczos 步的方案[14]，它对特征值不密集的结构相当有效，但当特征值密集时，本文采用的少Lanczos 步的方案似乎更为有利。

如只采用类似式(3)的特征值不等式控制迭代Lanczos 方法收敛，本文的方法速度可以提高30%到70%，但求解大量特征值时出现数个“幽灵”特征值。

3 结论

本文在广义特征值计算中实现了以下几个改进： 1) 用细胞解法替代了变带宽解法;

2) 用模态误差替代了特征值相对误差； 3) 循环展开用于特征值解法的各个环节。

本文的改进使得固有振动广义特征值计算的速度提高1个量级以上。细胞解法的引入改变了特征值解法中三角分解与迭代的计算量之比，多次移轴在大中型实用计算上成为可能；模态误差收敛判据使固有振动特征模态计算的过程变得平稳。

经验表明，模态误差比特征值误差更能反映特征值问题计算的精度。在计算较多模态时，模态误差应该作为首选的收敛判据。在实际工程计算中可将式(5)与式(3)的逻辑“或”作为收敛判据，并取

42

10λ?εε?=。这一做法保证了绝大部分特征向量是按模态误差收敛的。

在没有移轴的情形下，许多研究都认为Lanczos方法与Ritz向量法比子空间迭代快5到10倍[9,13]，但在允许以较小的代价实施移位时，本文的Lanczos方法与Ritz向量法仅比子空间迭代法快2倍左右。

我们期待在近期将稀疏解法与更精巧的固有振动广义特征值方法结合，进一步缩短计算时间，提高计算效率。

致谢作者特别感谢美国Virginia大学的秦钧博士关于Lanczos方法的通信以及部分源代码，作者感谢北京航天航空大学504教研室的硕士生陈展同学提供了ANSYS的算例BUAA1和BUAA2，并帮助进行计算。作者也感谢审稿人的建议。

参考文献

1. Chen P, Zheng D, Sun SL, Yuan MW. High Performance Sparse Static Solver in Finite Element Analysis with

Loop-unrolling. Advances in Engineering Software, 2003, 34:203-215

2. Yuan MW, Chen P, Sun SL, Meng ZL. Fast Sparse Static and Eigenvalue Solvers for Finite Element Analysis.

见:袁明武等主编.中国计算力学大会 ’2001论文集, 工程与科学中的计算力学. 北京: 北京大学出版社, 2001;101-111

3. Bathe KJ. Finite Element Procedure. Prentice Hall International, Inc. 1996

4. Clint M, Jennings A. The Evaluation of Eigenvalue and Eigenvectors of Real Symmetric Matrices by

Simultaneous Iteration. The Computer Journal，V ol.13, No.1, 1970; 10: 76-80

5. Bathe KJ, Wilson EL. Large Eigenvalue Problems in Dynamic Analysis. ASCE J Energ Meth Div,

1972;98:1471-1485

6. Wilson EL, Yuan MW and Dickens JM. Dynamic analysis by Direct superposition of Ritz vectors. Earthquake

Eng Struct Dyn, 1982; 813-832

7. Yuan MW et al. The WYD method in large eigenvalue problems. Eng Comput, 1989;6:9-57

8. Golub GH. Some Uses of Lanczos Algorithm in Numerical Linear Algabra. Topics in Numerical Analysis (ed.

by Miller JJH),1972:173-194

9. Underwood R. An Iterative Block Lanczos Method for the Solution of Large Sparse Symmetric Eigenproblems.

Ph.D. Dissertation. Stanford University, 1975

10. Press WH et al. Numerical Recipes in FORTRAN 90: the Art of Parallel Scientific Computing (2nd ed).

Cambridge University Press, 1996

11. 陈璞, 孙树立,袁明武. 有限元快速解法. 力学学报, 2002;34(2):216-222.

12. Duff IS, Grimes R, Lewis J. The User's Guide for the Harwell-Boeing Sparse Matrix Collection (Release I).

Technical Report TR/PA/92/86, CERFACS, Lyon, France, 1992.

13. Nguyen DT, Bunting CF, Moeller KJ, Runesha H, Qin J. Subspace and Lanczos sparse eigen-solvers for finite

element structural and electromagnetic applications. J Advances in Engineering Software, 2000; 31:599-606

14. Grimes RG, Lewis JG, Simon HD. A Shifted Block Lanczos Algorithm for Solving Sparse Symmetric

Generalized Eigenproblems. SIAM J Matrix Anal. Appl. 1994;15(1)