【JY】推开土木工程振型求解之兰索斯法(Lanczos法)的大门

一、写在文前

【前言】子空间迭代法可同时求解几个极端特征值和相应的特征向量，但它有收敛较慢，运算量较大，积累误差的缺点；随后，人们对其作了进一步的研究，出现了预处理子空间迭代法，这种方法的运算量较之子空间迭代法有所减小，但还是比较大，另外，当要求的特征值与不要求的特征值之间的间隔较大时，预处理子空间迭代法收敛也比较慢。因此需要有一种更高效的方法来求解，今天给大家带来的是兰索斯法(Lanczos法)振型求解思路方法。

关于子空间迭代法可以观看下列文章：

【JY】振型求解之子空间迭代

【JY|STR】求解器之三维结构振型分析

    目前科学计算、理论研究、科学实验并列为当今世界科学活动的三种主要方式。许多科学和工程领域如果没有科学计算不可能有一流的研究成果。而在工程领域，矩阵计算具有举足轻重的地位。

矩阵计算的基本问题有三个：

• 其一，求解线性方程组问题；

• 其二，线性最小二乘问题；

• 其三，矩阵特征值问题。

    在土木工程上，矩阵的特征值具有广泛的应用，如建筑结构中的固有频率分析问题以及钢结构的稳定性，建筑结构的振动问题，大型桥梁的颤振问题等等，都与矩阵的特征值密切相关。一个复杂结构的动力分析和稳定性分析可转化为大型稀疏矩阵的特征值/特征向量问题。

    求解结构振型和频率实质上是对以下式子的特征值和特征周期进行求解。但是，根据伽罗瓦理论，对于次数大于四的多项式求根不存在一种通用的方法。换句话说，对于大型矩阵并无法直接求解得到通法解析解，从而得到特征值和特征向量。于是，人们通过寻求其他的求解途径，提出了许多很好的思想方法和具体算法，促进了该学科的不断发展。

    本期介绍一种针对稀疏矩阵高效的特征值、特征向量计算方法——Lanczos迭代法。Lanczos方法存储量个大，计算速度较快，而且适用于求解矩阵的极端特征值。由于在实际问题中，矩阵的阶数虽然很大，但往往只是需要求其依模最大或最小的特征值，因此Lanczos方法得以广泛的应用，特别是适合求大型稀疏对称矩阵的部分特征值。

    对于Lanczos方法的本质：将实对称矩阵转换成为三对角矩阵（稀疏矩阵），从而便于计算机储存和后续的计算。（这样就把一个求实对称矩阵的特征值问题转化为求一个三对角矩阵的特征值问题。）

    有兴趣的小伙伴也可以了解下另外两种方法将一个实对称矩阵转化为三对角矩阵的方法Givens法及Householdes。

二、方程原理

令结构特征矩阵为：

    假定该结构特征矩阵是大型、稀疏的实对称矩阵，求其几个最大或最小的特征值的问题，可以用Lanczos算法解决，它的原理是先产生一个三对角阵T，于是问题转化为求一个三对角矩阵的特征值，变得相对简便。

    因此，我们的目的是将矩阵A转化为T，并且只要求得T的特征值和特征向量，则可以通过一定关系得到原结构特征矩阵的特征值和特征向量，也即可以得到结构振型和频率。

    首先对于一般结构来说，均需要进行初始向量的预设进行迭代，但是大部分振型都难以提前预估，我们可以听下威尔金森(Wilkinson)的建议预先不必猜测，而统一初始的向量为：

    然后进行预设需求的振型数量为 i≤n (结构特征矩阵维度)进行求解，以下流程图对于自己编写Lanczos方法中的T矩阵集成有较好的理解帮助。

    通过上述的迭代方法，并将求得的α和β带入矩阵T中，即可将结构特征矩阵是大型、稀疏的实对称矩阵化为一个三对角阵T。

    那么结构的频率和振型与该三对角矩阵T的关系是什么呢？

    接下来讨论下T矩阵和结构的频率和振型的关系，继上面的公式推导得到三对角矩阵T，由Gram-Schmidt正交化条件得到：