第四章 结构固有振动特征值问题的数值解

第四章结构固有振动特征值问题的数值解§4.1 概述根据结构振动的数学模型，即振动微分方程所形成的矩阵特征值问题，求解结构的固有振动特性——固有频率与固有振型，是结构振动分析的一个主要任务。

结构的固有振动特性是结构振动的内因。

固有振动特性也是进行结构振动响应分析和结构动力学设计的基础。

对于简单的结构，如均匀直梁、均匀直杆等，可以用解析的方法解得其固有振动特性。

对于一般结构，如果只需获得结构有限阶的固有振动特性，也可以采用试验测试（模态识别）的方法来获得。

但是对于大型复杂结构，不可能用解析分析的方法得到其固有振动特性，而采用试验测试的方法不仅花费高，而且周期长，对于处于设计状态的结构，显然也无法进行试验。

所以对复杂的工程结构，常用的方法是建立结构的数学模型，用数值求解的方法获得结构的固有振动特性。

随着计算技术飞速发展和特征值计算方法的研究进展，通过矩阵特征值问题的求解来获得结构固有振动特性，是已经被振动工程界普遍接受的一个有效和可靠的途径。

从数学理论上也可以证明，许多特征值计算方法具有相当好的精度，并且获得了实践和实验的证明。

由于结构固有振动特性求解与矩阵特征值求解问题的密切关系，在结构振动分析中，矩阵特征值问题已经成为结构固有振动特性分析的一个代名词。

所以在本章中，只要不作说明，一般讲的矩阵特征值问题就是指结构的固有振动特性求解问题。

所谓系统的特征值就代指结构的固有频率，特征向量代指结构的固有振型（固有模态）矩阵特征值问题的数值求解方法可以分为三类：矩阵分解法、迭代法和矩阵变换法。

由于矩阵（代数）特征值问题本身就是一个完整的系统，本章只能根据结构固有振动分析问题的需要，介绍一些常用的求解方法。

详尽的矩阵特征值问题的数值求解方法可以参考威尔金森的名著《代数特征值问题》。

本章的论述是建立在已经用有限元素法建立了结构振动运动数学模型的基础上。

§4.2 结构振动特征值问题的性质根据结构振动方程，可以得到结构固有振动的代数特征值问题：}]{[}]{[2x M x K ω=（4－1）或 }]{[}]{[x M x K λ= （2ωλ=） （4－2）振动特征值问题除了第二章所述的性质外，在特征值问题的数值求解中，还要用到如下一些性质： 1． 移轴特性对特征值问题}]{[}]{[x M x K λ= （4－3）若μ为一已知实数，则有：}]{)[(}]){[]([x M x M K μλμ-=- （4－4）新的特征值问题可写为：}]{[}]{ˆ[x M x Kρ= （4－5） ][][]ˆ[M K Kμ-= （μλρ-=） （4－6） 显然，上面两个特征值问题具有相同的特征向量，而特征值间的关系为：μρλ+=i i （4－7）μ称为移轴量。

在结构振动分析中，移轴特性常用来消除刚度矩阵的奇异性，也可用来加速迭代求解的收敛速度。

2． 特征值对合同变换的不变性［合同变换］：若][A 为一个非奇异矩阵，则变换]][[][A M A T 称为对][M 矩阵的合同变换。

对矩阵特征值问题：}]{[}]{[x M x K λ= （4－8）中的质量阵和刚度阵，用矩阵][A 作如下合同变换：]][[][]ˆ[A M A MT = （4－9） ]][[][]ˆ[A K A KT = （4－10） 则得到一个新的特征值问题}]{ˆ[}]{ˆ[φμφM K= (4－11) 从而有：][][][][]][[][]][[][]ˆ[]ˆ[=-=-=-A M K A A M A A K A M K TT T μμμ (4－12)由于0][≠A ，故有0][][=-M K μ （4－13）显然变换后的特征值不变：μλ= （4－14）且可以容易地证明，变换前后两个特征值问题的特征向量之间具有关系：}]{[}{φA x = （4－15）即合同变换不改变矩阵的特征值，特征向量具有转换关系。

3． 特征值对相似变换的不变性［相似变换］：对非奇异阵][A ，变换]][[][1A K A -称为相似变换。

对矩阵特征值问题：}]{[}]{[x M x K λ= （4－16）中的质量阵和刚度阵，用矩阵][A 作如下相似变换：]][[][]ˆ[1A M A M-= （4－17） ]][[][]ˆ[1A K A K-= （4－18） 则得到一个新的特征值问题：}]{ˆ[}]{ˆ[φμφM K= （4－19） 从而有：][][][][]][[][]][[][]ˆ[]ˆ[111=-=-=----A M K A A M A A K A M K μμμ （4－20）由于0][≠A ，故有0][][=-M K μ （4－21）显然变换后的特征值不变：μλ= （4－22）且可以证明，变换前后的特征向量间具有关系：}]{[}{φA x = （4－23）即相似变换不改变矩阵的特征值，特征向量具有转换关系。

4． 特征值对正交变换的不变性［正交变换］：在相似变换中，若对于非奇异矩阵][A ，有][][][I A A T =，即1][][-=A A T ，则用][A 进行的相似变换称为正交变换。

正交变换不仅不改变矩阵的特征值，而且正交变换后矩阵的对称性不变。

5． 特征值对旋转变换的不变性［旋转变换］：取变换矩阵为：⎥⎦⎤⎢⎣⎡-=θθθθcos sin sin cos ][R （4－24） 可以证明][R 为正交矩阵，用][R 矩阵进行的相似变换，称为旋转变换。

通过一系列的旋转变换，可以使对称阵变为对角阵，使一般矩阵变为三角阵。

利用上述特征值问题的性质，可以得到不同的特征值求解方法。

结构振动特征值的数值求解中常用的方法主要有： 【多项式迭代法】多项式迭代法是利用特征值使特征多项式等于零的性质，即0][][)(=-=M K p λλ （4－25）【利用特征多项式的Sturm 序列性质的方法】通过Sturm 序列的性质，分离出每个特征值所在的区间，通过二分法和其它加速查找方法，逐步缩小特征值所在的分隔区间，最后得到所需的特征值。

【矢量迭代法】矢量迭代法利用的基本关系式为：}]{[}]{[i i i x M x K λ= （4－26）具体的迭代方法有：正迭代法、逆迭代法、子空间迭代（同时迭代）法。

【矩阵变换法】矩阵变换法利用的是“矩阵特征值对合同变换的不变性”这一性质。

利用具有某种性质的变换矩阵，经过一系列变换，得到一个形式简单的、容易进行特征值求解的矩阵特征值问题。

常用的有雅可比方法：它是经过一系列旋转变换，最后使一个对称阵成为一个对角阵，然后很容易地求解该对角阵的特征值。

QR 法：QR 法是经过一系列的相似变换，将一个实数矩阵转化为一个三角阵。

然后求解其特征值。

豪斯霍尔德方法：豪斯霍尔德方法 是通过正交相似变换，将一个实对称矩阵化为三对角阵。

然后用QR 法求其特征值。

还有其它一些方法。

本章重点介绍特征多项式迭代法、基于Sturm 序列性质的方法、矢量迭代法、子空间迭代法，对于矩阵变换方法仅对变换格式作简略介绍。

在求解结构固有振动特性时，有时各种方法也可以联合使用，来加快求解的速度。

§4.2 多项式迭代方法使特征多项式][][)(M K p λλ-=（4－27）等于零的根就是特征值问题}]{[}]{[x M x K λ=（4－28）的特征值解。

当阶数4>n 时，无法用公式求根，只能用数值计算方法求解上方程。

而求解方程0][][)(=-=M K p λλ （4－29）的方法通常有显式和隐式两类。

并且一般都采用“频率扫描”方法求解。

【显式迭代】将)(λp 展开成λ的多项式，然后进行迭代求解，因而，如果][],[M K 阵阶数很大时，将多项式展开就有一个缺点——方程的系数计算误差引起根的稳定性很差，因而在实际结构振动分析中很少采用。

【隐式迭代】先给定λ值，对矩阵（][][M K λ-）进行高斯三角分解，再计算出)(λp 的值。

即做分解]][[][][S L M K =-λ（4－30）][L 和][S 分别是单位下三角阵和上三角阵。

从而：∏==⋅=-=ni ii s S S L M K p ][][][][][)(λλ （4－31）实际计算时，是从λ的初值0λ开始，依次计算)(0λp ，)(0λ∆λ+p ，)2(0λ∆λ+p ， )(0λ∆λn p +的值，根据它们的正负号，确定出根所在的区间（隔根区间），然后对各隔根区间用二分法（或0.618优选法）进行细分直到求出符合给定精度的解，也可以用更快的割线法来加速寻根：)()()()()1()()1()()()()1(--+---=r r r r r r r p p p λλλλλλλ（4－32）但要注意，使用割线法时，迭代的初值很重要，一般与二分法联合使用，以保证迭代的迅速收敛。

【例】计算矩阵特征值问题}]{[}]{[x M x K λ=的第一个特征值1λ。

⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----=210141012][K ，⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=5.000010005.0][M 【解】由于矩阵][],[M K 都是正定阵，故取迭代初值1)0(-=λ，0)1(=λ25.2646/1050015/230015.25.21015115.2)1(=--=----=-p 127/120012/7001221141012)0(=--=----=p 8421.0)]1(0[25.2612120)2(=----=λ7150.4)8421.0(=p ，3871.1)3(=λ，8467.1)3871.1(=p 。

9870.1)6(=λ，026347.0)9870.1(=p ，9993.1)7(=λ迭代法的优缺点：理论简明，方法简单，可以对任一个特征根进行迭代搜索，但计算量大，且容易漏根。

§4.3 基于Sturm 序列性质的方法 【Sturm 序列】对广义特征值问题}]{[}]{[x M x K λ= （4－33）其对应的特征多项式为：][][)(M K p λλ-= （4－34）若将][K ，][M 中的最后r 行和最后r 列去掉，得到相应的一个新广义特征值问题：}]{[}]{[)()()()()(r r r r r x M x K λ= （4－35）称为原系统的“第r 阶伴随约束系统的特征值问题”，这相当于对原系统增加了r 个约束。

定义第r 阶伴随约束系统的特征多项式：][][)()()()()()(r r r r r m k p λλ-= 取)1(,2,1,0-=n r ，得到n 个伴随特征多项式：[]M K p m m m m k k k kp m k p n n n n n n λλλλλλ-=⎥⎦⎤⎢⎣⎡-⎥⎦⎤⎢⎣⎡=-=------][)()()()0()0(22211211)2(22211211)2()2(11)1(11)1()1(（4－36）这n 个特征多项式组成的序列称为Sturm 序列，根据瑞利约束原理，第1+r 阶伴随系统与第r 阶伴随系统特征值)1(+r i λ和)(r i λ之间有关系：rr n r r n r r n r r r r -+----++≤≤≤≤≤≤λλλλλλλ)1(1)(1)1(2)(2)1(1)(1 （4－37）这一性质又叫特征值隔离定理，在第一章已经介绍过。