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 方法尤为严重，原因是采用逆迭代技术时，高阶的、密集的特征值不易分离。

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

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

||||||||k k k k kk k k k K M K M K M ϕϕλϕϕλϕεϕλϕ−−≈≤%%%%%%%%%(5)模态误差有明显的物理意义，k K ϕ是振型k ϕ的最大弹性节点力，而kk M λϕ%%是振型k ϕ的最大惯性节点力。

式(5)的最左端是非平衡节点力与最大弹性节点力之比，中间是非平衡节点力与最大惯性力之比[3]。

大量算例表明，在模态误差的意义之下，收敛过程平稳，三种解法效率相差不是太大。

迭代Ritz 向量法平均最快，子空间迭代法最慢。

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

1 算法描述1.1 子空间迭代法 最初是由Clint 和Jennings 提出，是反幂法的推广[4]。

稍后，Bathe 和Wilson 在其中加入了子空间上的Rayleigh-Ritz 过程。

它可以明显地改善收敛速度[4,5]。

以下是一个子空间迭代算法的主要步骤。

I ．初始化1. 确定子空间的维数q2. 选取初始向量矩阵 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 与生成步数r2. 选取初始向量矩阵0N q Q R ×∈3. 设定每次移轴的最大迭代次数max I II ．移轴1. 计算移轴µ，应设法保证它不是特征值2. 分解移轴刚度矩阵T K M LDL µ−=3. Sturm 序列校核III ．迭代max I 次，完成后转向II1. 对0,1,,k r =L 解 1T k k LDL Q MQ +=，然后用将 1k 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 次，完成后转向II1. 对0,1,,k r =L 解 1T k k LDL Q MQ +=，然后用将 1k Q +对已收敛的特征向量以及12,,,k Q Q Q L 作M-正交归一化，并形成1k Q +。

在此过程中依次形成*T2. 求解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 校核均取后验方式，移轴的首选为下两个待收敛特征值之中点，即120.5()k k µλλ++=+%% (8)若失败，则选为10.98()k k k µλλλ+=+−% (9)其中kλ是最后一个已收敛的特征值，1k λ+%，2k λ+%是下两个待收敛的特征值。