建筑结构丨惊！结构人搞了个比牛顿法快60多倍的算法，快准稳！|容差|牛顿|算子|算法

来源：哲设计。

https//doi.org/10.1002/nme.6456

这个算法其实只有三个特点：

快、准、稳

1快

先解释一下标题里的“60多倍”是怎么测算的。

为了便于大家验证，我们用了陆老大放在个人主页上的一栋42层钢筋混凝土框架-核心筒结构的OpenSees模型[下载链接见文末]。它有2万多个结点和将近4万个单元。

在OpenSees中，分别采用NewtonLineSearch和我们的新算法ExpressNewton计算结构在调幅至PGA=5.1m/s2的1940 El Centro NS记录[下载链接见文末]作用下的地震反应。为公平起见，二者均采用Newmark-β积分算子和0.01s的固定时间步长，且均未使用并行计算。此外，在NewtonLineSearch中使用0.001m的位移容差。

模型的顶点位移、速度和加速度反应如下图所示。可见二者几乎完全吻合。但是为了得到这样的结果，NewtonLineSearch耗时272.26小时（11.3天），而ExpressNewton只用了4.36小时，是前者的1/62.4。

这就是标题里“60多倍”的由来。

当然，计算效率提升的幅度也和模型的复杂程度和反应的非线性程度等诸多因素有关。欢迎大家用自己的算例来测试。那么，这么残暴的效率，是怎么实现的呢？

先来复习一下迭代求解非线性方程的牛顿（Newton）法。

通常采用逐步积分法求解结构的非线性动力反应，也就是把整个时间历程离散成很多小的时间段，一段一段逐步求解。这里以Newmark-β法作为时域积分算子，比如我们已经知道了第i时刻的力和位移以及第i+1时刻的力，那么牛顿法在一个时间步内不断利用当前的广义刚度矩阵J去逼近第i+1时刻的真实位移，直到迭代计算的结果满足某一容差要求，比如下图中ε。

好像说得有点儿复杂不太符合本号的调性了，但我实在想不出更直观的解释方法了。

重点是，牛顿法在每个时间步内都需要多次更新广义刚度矩阵J，很花时间。为了省时间，人们提出修正牛顿法：不在每次迭代的时候都更新广义刚度矩阵，而在一个时间步内始终采用时间步开始时刻的广义刚度矩阵。

即使如此，修正牛顿法仍然需要在每个时间步开始的时候更新一下广义刚度矩阵J，还是挺花时间的，并且为了满足和牛顿法相同的容差，在一个时间步内往往需要进行更多次迭代计算。

下面隆重推出简单粗暴的极速牛顿法：

（1）在所有时间步内都只使用零时刻的初始广义刚度矩阵J0进行迭代；

（2）每个时间步内只迭代两次就直接进入下一时间步，不检查容差。

更新广义刚度矩阵往往是非线性动力分析中最耗时的工作，并且随着矩阵规模的增大，计算成本指数上升。极速牛顿法只在最开始的时候计算一次广义刚度矩阵，而在随后的整个时程分析中不再更新广义刚度矩阵，想不快都难。

江见鲸老师有个比方，牛顿法是走一步绕一个圈，但是步子大；修正牛顿法是走几步绕一个圈。这么说来，我们的算法就是不用绕圈直接跑。

不过，这也太粗暴了吧？能算得准吗？不怕发散吗？

2准

上面的例子其实已经说明了，4个小时算出来的结果和11天的结果相差无几。

怎么做到的？你以为在每个时间步里面只迭代两次是拍脑袋拍出来的吗？为什么不是一次，为什么不是三次，而偏偏是两次？

对，确实是拍脑袋拍出来的。然而，作者在拍脑袋之余还给出了证明：迭代两次就够了。

怎么就算够了呢？说好的精益求精呢？直接把容差取消了，这样真的好吗？

非线性动力反应的求解至少有两步近似。

首先，逐步积分法把一个连续的反应时程切成许多小段，然后在一系列离散的时间点上求解结构的反应状态（位移、速度、加速度）。从已知的反应状态推算下一个时刻的未知状态，是时域积分算子的任务。这时需要引入一些假设，比如常用的Newmark-β法假设加速度反应在时间步内为常数或线性变化。引入假设的同时也引入了误差，即时域积分算子的误差，姑且称之为积分误差。

其次，对于非线性体系，还需要在时间步内求解一个非线性方程。上面说的牛顿法、修正牛顿法和极速牛顿法等都是在做这件事。这时还会引入一定的误差，即非线性求解算法的误差，姑且称之为算法误差。

总的误差取决于积分误差和算法误差里面较大的那一个，有点儿像短板效应。仍以Newmark-β法为例，它具有3阶精度。可以证明，极速牛顿法具有2N阶精度，N是迭代次数。因此，当与Newmark-β法配合使用时，迭代两次（N=2）就够了。如果与其他更高阶的积分算子配合使用，可以通过增加迭代次数得到够用的精度阶。

积分误差和算法误差就像两个人在森林里遇到熊，你不需要比熊跑得快，只要比另一个人跑得快就够了。

初学有限元的同学经常恨不得把算法的容差设置得无穷小，似乎容差越小结果就越精确似的，其实不然。总要被“无法收敛”折磨千百遍之后才慢慢学会放过容差。

不过现在好了，极速牛顿法不需要设置容差了。

3稳

众所周知，Newmark-β法是无条件稳定的。简单粗暴的极速牛顿法会不会破坏这个美好的性质呢？

可以证明，当下式成立时，极速牛顿法不影响Newmark-β法的无条件稳定性。

式中，K0是结构的初始刚度矩阵，KT是切线刚度矩阵，σ>=1是针对强化型结构打的一个补丁。

对于如上文图示的常见的软化型结构，即使σ=1，上式也总是成立的。当结构表现出强化行为时，可以首先选用恰当的σ保证无条件稳定性，然后在极速牛顿法的迭代计算中将初始广义刚度矩阵J0中的刚度矩阵K0放大σ倍。如下图所示。

如果仅仅是办公室算法，再怎么吹也没有用。为了方便大家评测，我们把极速牛顿法放进了OpenSees，预计将于近期随OpenSees3.2.2发布。届时，只要在你的动力分析步中将algorithm改成ExpressNewton并开启-initial开关，就可以立即体验极速牛顿法带来的快准稳。具体语法如下：

algorithm ExpressNewton <$iter><-initial><$kMultiplier>

$iter：迭代次数

-initial：使用初始广义刚度矩阵进行迭代的开关

$kMultiplier：刚度矩阵放大系数σ，默认为1.0

例如，当与Newmark积分算子配合使用时，可以使用以下命令行。其中，ExpressNewton之后的第一个参数是迭代次数N。

algorithm ExpressNewton 2 -initial 1.0integrator Newmark 0.5 0.25analysis Transient

这篇之所以吹得比较没有顾忌，主要原因是我只是打个酱油。下面两位才是本文的技术担当。

第一作者：徐俊杰，中国地震局工程力学研究所副研究员，毕业于大连理工大学和清华大学。

通讯作者：黄羽立，见之前，毕业于清华大学和加州大学伯克利分校。

RC框筒结构的OpenSees模型（陆新征个人主页）：http://www.luxinzheng.net/download/OpenSeesTallBuildings.zip

1940 El Centro记录：

http://www.vibrationdata.com/elcentro.htm