特征值解法

《结构动力学》大作业结构大型特征值问题的求解0810020035 吴亮秦1振动系统的特征值问题1.1实特征值问题n 自由度无阻尼线性振动系统的运动微分方程可表示为：[]{}[]{}()M u K u F t += （1.1）其中，{}u 是位移向量，[]M 和[]K 分别是系统的质量矩阵和刚度矩阵，都是n 阶正定矩阵，()F t 是激励向量。

此系统的自由振动微分方程为[]{}[]{}0M u K u += （1.2）设其主振型为： {}{}sin()u v t ωϕ=+ （1.3） 其中，{}v 为振幅向量，ω为圆频率，ϕ为初相位。

将(1.3)代入自由振动微分方程(1.2)， 得：[]{}[]{}K v M v λ= （1.4） 其中2λω=，(1.4)具有非零解的条件是()[][]det 0M K λ-= (1.5)式(1.4)称为系统的特征方程，由此可以确定方程的n 个正实根1{}n i i λ=，称为系统的特征值，1{}n i i ω=称为系统的固有频率，{}i v （i=1,2,…..n ）为对应于特征值的特征向量或称为系统的振型或模态。

因为[]M 矩阵正定，则[]M 有Cholesky 分解：[][][]TM L L = (1.6)其中，[]L 是下三角矩阵。

引入向量{}x 满足：{}[]{}Tx L v =，则：1{}([]){}T v L x -= (1.7) 代入(1.4)，得：([][]){}0I P x λ-= (1.8)其中，()11[][][][]TP L K L --=，式(1.8)称为标准实特征值问题。

1.2复特征值问题多自由度阻尼自由振动系统的运动方程为如下二阶常系数微分方程组：[]{()}[]{()}[]{()}0 M x t C x t K x t ++= (1.9) 其中 []M ，[]C ，[]K 分别是n 阶的质量、阻尼和刚度矩阵，{()}q t 是n 维可微向量函数。

用分离变量法，设{()}{}tx t e λφ=，其中{}φ是与时间t 无关的常向量，λ为待定参数。

将{()}{}tx t e λφ=代入上述齐次方程，得确定参数λ,{}φ的特征方程：()2[][][]{}0M C K λλφ++= (1.10)(1.10)具有非零解的条件是()2det [][][]0M C K λλ++= (1.11)式(1.11)的求解就是复特征值问题。

2 实特征值求解方法 2.1特征方程法求解实特征值最直接的方法就是特征方程法，即把式(1.5)展开得到特征值多项式：11100n n n C C C λλλ--++++= (2.1)求解(1.10)即得特征值，特征方程法仅适用于3n ≤的低阶情况特征值求解，并不是求解特征值的一般方法，实际求解大型结构的实特征值问题的方法很多，归纳起来大致分为两类，即基于矩阵相似变换原理的相似变换法和基于瑞利—李兹（Rayleigh-Ritz ）法、迭代方法的瑞利—李兹类方法。

2.2相似变换法相似变换法基于矩阵的相似变换原理，即相似矩阵与原矩阵有相同的特征值，主要有以下几种方法。

2.2.1雅克比法（Jacobi ）雅克比法是求解实对称矩阵全部特征值的简单有效的方法。

它的基本思想是，通过正交 相似变换使矩阵对角化，从而求出矩阵的全部特征值，进行步骤如下：（1） 将系统的自由振动微分方程：[]{}[]{}0M u K u += 化为标准实特征值形式：([][]){}0I P x λ-=；(2)找矩阵[P]中绝对值最大的非对角线元素0uv p ，构造正交矩阵[]1S ，对[P]矩阵作正交变换[][][][]T111S P S =P 使元素0uv p 化成0；（3）找矩阵[]1P 中的绝对值最大非对角线元素1uv p ，构造正交矩阵[]2S ，对[]1P 矩阵作正交变换[][][][]T2122S P S =P 使元素1uv p 化成0；（4）依次进行下去，每次找矩阵[]i P 中的最大非对角线元素i uv p ，构造正交矩阵[]i+1S ，对[]i P 矩阵作正交变换[][][][]Ti+1i i+1i+1S P S =P 使元素i uv p 化成0。

重复使用该变换，每次变换可使矩阵[]i P 更接近于对角矩阵，若干次变换厚，原矩阵[P]化成对角矩阵，对角线元素即是原矩阵的特征值。

每步变换的关键在于构造正交矩阵[]i+1S ，实际采用吉文斯（Givens ）旋转矩阵，通过多次坐标系的旋转来实现原矩阵的对角化，[]i+1S 取如下形式：[]11cos sin 11sin cos 11θθθθ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦i+1S 其中，u 、v 为[]i P 中最大非对角线元素所在的行列序号，θ为旋转角，由正交变换[][][][]Ti+1i i+1i+1S P S =P 使矩阵[]i+1P 的元素uv p =vu p =0来确定，计算得到：arctan 2i i vv uu iuvp p p θ-1=2.2.2 豪斯厚德三对角化法（Householder ）当矩阵[]P 的阶数较高，非对角元素较稠密，且数值较大时，使用雅克比法的收敛速度就不快，为了谋求更快更有效的算法，吉文斯（Givens ）首先提出了一个将实对角矩阵三对角化的方法，将矩阵三对角线外的元素利用正交变换逐一化成0.而豪斯厚德（Householder ）则改进了这一方法，即将三对角线外的元素逐行变换成0，最后将矩阵三对角化，在三对角化的基础上便可结合别的方法迅速地求出全部或部分特征值。

豪斯厚德方法基本思路如下：（1）寻找一个对称的正交矩阵[]1S ，通过正交变换[][][][]T=111S P S P ，将矩阵[]P 的最后一行实现三对角化，因为[]P 是对称矩阵，这样同时也实现了[]P 的最后一列三对角化；（2）寻找一个对称的正交矩阵[]2S ，通过正交变换[][][][]T=2122S P S P ，将矩阵[]1P 的倒数第二行（列）三对角化；（3）寻找一个对称的正交矩阵[]i S ，通过正交变换[][][][]T =i i-1i i S P S P ，将矩阵[]i-1P 的第n-i+1行（列）三对角化；（4）依次进行下去，经过n-2次正交变换后，原矩阵[]P 即变换成了三对角矩阵[]n-2P 。

与雅可比法一样，上述过程的关键就是如何构造每步变换的正交矩阵[]i S ，豪斯厚德的主要贡献就是构造出了满足这一要求的矩阵，表达如下：[][]{}{}/Ti i i u u H =-i S I其中，{}{}12Ti i i H u u =，[]I 为单位矩阵，{}i u 是列向量，取值为： {}(),1,2,1,,0Ti i i i l l l l u pp p-=1l n i =-+是当前处理的行列号，,i l j p 是矩阵i P 的第l 行第j 列元素。

=当,10i l l p -≥时，前面取正号，当,10il l p -<时，取负号。

2.2.3斯特姆排序法（Sturm ） 将矩阵三对角化后，运用斯特姆排序法计算三对角矩阵的全部或部分特征值时非常有效的，该方法基于如下思想：经三对角化后，矩阵P 的特征方程方程为：[]111222311000[]000n n n c b b c b b c P I b b c λλλλλ------==-定义它的零阶主子式为：1=0p以后各阶为首的主子式为：111221112212()()()(r 2,3...,)r r r r c c b c c b b c c p b p n λλλλλλ--⎫=-⎪-⎪==---⎬-⎪⎪=--=⎭12r p p p则0p ，1p ，2p ,……n p 构成一个多项式序列（斯特姆序列）。

任何一个斯特姆序列都与一个―整数函数‖ λS()相联系。

当λ取任意一个实数值，斯特姆序列的各项多项式都可分别算出具体数值：λ0p ()，λ1p ()……λn p ()，而从λ0p ()到λn p ()所发生的符号改变数就是λS()之值。

在数学上可证明，λS()值即为矩阵[P]的小于λ的特征值的个数。

例如，由(2.5)2S =，即说明矩阵[P]有两个特征值是小于2.5的；若(0)0S =，说明矩阵[P]没有负特征值。

还可以证明，仅当自变量λ的值是矩阵[P]的某一特征值时λS()值才会发生改变，例如，若有(2.1)2S =，(2.2)3S =，则矩阵[P]肯定有一特征值介于2.1和2.2之间，只要将λ值分得足够细，总可以根据λS()的变值特性按任意精度计算处矩阵[P]的某一特征值。

2.2.4 QL(或QR)法实对称矩阵P 也可借助于QL 算法通过多次正交变换而按指定的精度逼近对角矩阵。

如果事先将矩阵P 变换成三对角化矩阵，然后再施行QL 变换则是特别有效的，进行QL 变换变换的基本步骤如下：（1）寻找正交矩阵[]1Q ，使：[][][]L =11P Q[]L 1为一个下三对角矩阵（若[]L 1取上三角矩阵则是QR 法）。

（2）由[][][][][]L -1T111=Q P =Q P 作正交变换：[][][][][][]L ==T11111P Q Q P Q（2） 矩阵[]1P 作相同分解：[][][]L =122P Q再做变换：[][][][][][]L ==T222212P Q Q P Q（4）这样反复按下式进行：[][][]L =s-1s s P Q[][][][][][]L ==Ts s s s s-1s P Q Q P Q当→∞s 时，矩阵[]s P 可按任意精度逼近对角矩阵。

与前面介绍的方法一样，QL 法的关键也是在每步变换中寻找到满足条件的正交矩阵[]s Q ，[]s Q 可由n-1个正交矩阵[]si Q (i=1,2,……n -1)组成。

每个[]si Q 的功能是消去矩阵[]s P中一个非对角线元素，实际[]si Q 取为吉文斯旋转矩阵。

2.3 瑞利—李兹类方法上述求解实特征值的第一大类方法-矩阵变换法可求出系统的全部特征值，而第二大类方法—瑞利—李兹法及其相关方法主要用于求解大型系统的部分特征值的近似值。

2.3.1瑞利—李兹法（R-R ）Rayleigh-Ritz 法（简称R-R 法）是一种缩减自由度的方法，用于求解大型系统部分特征值的近似值，R-R 法在理论上的根据是，各阶特征值是Rayleigh 商的极值或驻值。

其基本思想是：选择m 个线性无关的Ritz 基向量{}i q （i=1，2，……m ），它们张成一个m 维子空间m V ，而Rayleigh 商在子空间m V 中存在m 个极值点。

这m 个极值点就是系统前m 阶特征值在m V 子空间的最佳近似值，同时也求出相应的特征向量近似值，R-R 法的具体做法如下：（1）令2ωλ=，把特征方程2([][]){}0K M x ω-=化成[]{}{}[]K x x M λ=； （2）把结构体系自由度折减为m 个，选取m 个线性无关的Ritz 基向量，则向量{}x 可表示成m 个基向量的线性组合：{}[]{}1{}ni i x z Q z ===∑i q其中，[]{}{}{}Q =⎡⎤⎣⎦ 12m q ，q ，q ，{}[]12m z z z z = ，，分别称为Ritz 基向量和Ritz 坐标向量；（3）给出向量{}x 的Rayleigh 商：**{}[]{}{}[]{}{}[]{}{}[]{}T T T T x K x z K z x M x z M z ρ= = 式中，[][]*[][]TK Q K Q =、[][]*[][]TM Q M Q =，分别是矩阵[]K 和[]M 在Ritz 基向量所张子空间的投影，都是m m ⨯阶矩阵； （4）由（3）式可见，Rayleigh 商是Ritz 坐标向量的函数，因此它的极值条件可表示为：{}0ρ∂=∂z 将（3）式代入，得到：{}{}**[][]i i i K z M z ρ=这是m 阶广义特征值问题，它的特解为{}i z 和i ρ（i=1，2，……m ），则前m 阶特征值的近似值：i i λρ=。