一、写在文前

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

往期精彩

#性能分析

【JY】基于性能的抗震设计浅析（一）

【JY】基于性能的抗震设计浅析（二）

【JY】浅析消能附加阻尼比

【JY】近断层结构设计策略分析与讨论

【JY】浅析各动力求解算法及其算法数值阻尼(人工阻尼)

理念

【JY|体系】结构概念设计之(结构体系概念)

【JY|理念】结构概念设计之(设计理念进展)

【JY】有限单元分析的常见问题及单元选择

【JY】结构动力学之显隐式

【JY】浅谈结构设计

【JY】浅谈混凝土损伤模型及Abaqus中CDP的应用

#概念机理

【JY】基于Ramberg-Osgood本构模型的双线性计算分析

【JY】结构动力学初步-单质点结构的瞬态动力学分析

【JY】从一根悬臂梁说起

【JY】反应谱的详解与介绍

【JY】结构瑞利阻尼与经济订货模型

【JY】主成分分析与振型分解

【JY】浅谈结构多点激励之概念机理（上）

【JY】浅谈结构多点激励之分析方法（下）

【JY】板壳单元的分析详解

【JY】橡胶支座的简述和其力学性能计算

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

【JY】橡胶支座精细化模拟与有限元分析注意要点

#软件讨论

【JY】复合材料分析利器—内聚力单元

【JY】SDOF计算教学软件开发应用分享

【JY】Abaqus案例—天然橡胶隔震支座竖(轴)向力学性能

【JY】Abaqus6.14-4如何关联fortran？

【JY】如何利用python来编写GUI？

【JY】如何解决MATLAB GUI编程软件移植运行问题？

【JY】浅谈结构分析与设计软件

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

【JY】SignalData软件开发应用分享

【JY】基于Matlab的双线性滞回代码编写教程

#其他

【JY】位移角还是有害位移角？

【JY】如何利用python来编写GUI？

【JY】今日科普之BIM

~关注未来更精彩~

欢迎关注微信公众号: 建源学堂