首页/文章/ 详情

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

详细信息

作者优秀: 优秀教师/意见领袖/博士学历/特邀专家/独家讲师
平台推荐: 主编推荐

仿真圈

3年前浏览2703

点击图片报名

一、写在文前

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

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

其一，求解线性方程组问题；
其二，线性最小二乘问题；
其三，矩阵特征值问题。

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

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

本期介绍一种针对稀疏矩阵高效的特征值、特征向量计算方法——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迭代法。

四、写在文后

本文课程相关模型和资料，我将分享仿真秀土木工程行业学习包（持续更新），感兴趣的朋友可以自行下载。即日起，仿真秀平台将开启土木工程有限元分析交流月，将给理工科学子和研发工程师带来多场精彩直播和社群讨论，欢迎报名。如下图：

在文章末尾点赞、在看和分享到朋友圈，点击附件下载仿真秀平台土木工程行业学习包（长期永久反复更新），另外诚邀你加入仿真秀土木工程行业交流群，不定期邀请你参加土木工程相关答疑直播，且抱团一起学习理论、软件和行业应用哦。

（完）

作者：建源之光，工学博士，结构工程-结构振动控制方向研究。高级技术经理，自主研发软著 4 项，专利 21 项。主要从事超高层建筑结构设计分析，减隔震设计分析，舒适度及其他振动控制分析，减隔震(振)元件设计开发。建源学堂，土木与有限元优质内容创作者，擅长超高层建筑结构设计分析，减隔震设计分析，舒适度及其他振动控制分析，Sap2000，Etabas，Abaqus，Matlab，OpenSeespy。