【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正交化条件得到： 将Gram-Schmidt正交化条件带入上式中可得到： 考虑列向量间的正交性，并将得到的上述公式往会带入可得到： 得到这个标注了三星的重要级公式！其中标记红框部分是左乘部分。若将红框这么一画，就变成了： 再通过上述这样一个化解，原结构和变换后的T的频率和振型的关系一目了然。若令原结构的特征向量（振型）为：
也就是求解T矩阵的完全解为，则该完全解的特征向量与原结构的特征向量(振型)的关系式之间相差一个系数矩阵X，而特征值(频率)是原结构的倒数。因此取Lanczos向量的个数i等于系统的阶数n，则Lanczos变换法给出了原特征值问题的精确解。 一般情况下我们取 i ≤ n 。由于组装结构特征矩阵时候，刚度求逆变换了一次，相当于进行了一次逆迭代，因此Lanczos变换法给出了系统低阶特征值的很好的近似解。值得注意的是，为了保证求解精度，在迭代过程中可以用Gram-Schmidt对迭代向量进行重新正交化，并采用移轴法可提高其效率。 至此，该问题变为求一个i阶三对角阵T的特征值和特征向量。
对于求该降阶的特征值、特征向量，有许多方法如：比较高效的QR分解法等等，这里就不再多赘述，小伙伴们可以自行编程尝试。 建立一个无楼板的简易框架结构，体验一把Abaqus的兰索斯算法，并对比下自编的Lanczos。模型如下：
求解前4阶模态
利用Abaqus 直接提取该结构的质量矩阵、刚度矩阵、阻尼矩阵(非必要)，一个简单的结构密密麻麻的数，感兴趣的小伙伴可以试一试：
并将这三个矩阵放入自行编写的动力计算器（由于该动力计算器尚未完善，所以暂不公开使用，小伙伴们按这个思路也可以自己编写得到属于自己的求解器。）
可以看得到结果非常接近，误差来源可能来源于计算机计算的舍入误差。而且计算速度非常高效，因此针对稀疏矩阵高效的特征值、特征向量的计算方法可以采用可存储量个大，计算速度较快的Lanczos迭代法。
