动力学中非线性常微分方程的一种新型二阶显式积分方法

《Mathematics》：A Novel Second-Order Explicit Integration Method for Nonlinear Ordinary Differential Equations in Dynamics Gorka Urkullu, Ibai Coria, Igor Fernández de Bustos and Haritz Uriarte

　　

本研究提出了一种面向工程应用中常见的二阶常微分方程（Ordinary Differential Equations, ODEs）的新型显式积分方法。传统上，这类问题通常通过将其改写为一阶方程组以应用单步法（如龙格-库塔（Runge–Kutta）格式），或采用在线性动力学中广泛应用的直接二阶方法（包括广义-_α法、中心差分法（Central Difference, CD）和纽马克法（Newmark））来求解。所提出的方法基于截断至三阶导数的泰勒级数（Taylor series）展开推导而来，形成一种完全显式算法，每个时间步仅需一次函数求值。与纽马克表述类似，该方法包含可调参数，允许研究人员在精度与稳定性之间进行权衡。对于特定的参数选择，该方法展现出与中心差分方案相当的收敛性和稳定性。一个重要优势在于，即使非线性项依赖于一阶导数项，该方法仍保持显式特性。论文从理论上分析了稳定性、局部截断误差、谱性质、数值阻尼和周期伸长等特性。该方法通过四个多体动力学中的测试用例（包括线性和非线性问题）进行了验证。结果表明，显式积分三阶法（Explicit Integration Grade 3, EIG-3）在达到与现有显式二阶积分器相当精度的同时，显著降低了计算成本，尤其在非线性应用中表现突出。

## 研究背景与问题提出

数值积分是计算数学的基石，在结构力学和多体动力学领域，二阶常微分方程源于连续域的空间半离散化，有限元法（Finite Element Method, FEM）是最常见的来源。对于以低频模态主导的刚性系统，隐式积分器如纽马克法、希尔伯-休斯-泰勒法（Hilber–Hughes-Taylor, HHT）和广义-_α法是首选，因其具有无条件稳定性和可控的数值耗散特性。然而薄的反，显式方法在涉及高频现象的场景中不可或缺，如波传播、冲击动力学和复杂非线性多体系统。在这些高维应用中，隐式求解器为确保收敛所需的迭代开销计算上难以承受。

现有显式积分器面临严重的算法权衡：标准中心差分法对线性系统稳健，但当存在速度相关非线性（如非线性阻尼或科里奥利力）时，会变得隐式耦合，需要迭代求解器，破坏其严格显式优势；二阶显式龙格-库塔法（RK2）虽能有效处理非线性，但每步需要多次函数求值，在大规模仿真中成本翻倍；显式线性多步法（如亚当斯-巴什福斯（Adams–Bashforth）族）随精度阶数增加，绝对稳定区域缩小，在对鲁棒稳定性要求高的高频应用中常常失效。为克服这些局限，Gorka Urkullu、Ibai Coria、Igor Fernández de Bustos和Haritz Uriarte开发了EIG-3，一种基于诺斯ieck（Nordsieck）多项式近似和泰勒高阶展开框架的新型显式积分方法。

## 核心技术与方法

研究人员采用的关键技术方法主要包括：基于三阶截断的泰勒级数展开作为数学基础，通过引入可调参数₀和₁构建参数化框架以调节稳定性和谱特性，利用 amplification matrix分析稳定性与收敛性，通过Faddeev-LeVerrier法计算特征多项式，以及采用Horner方法进行多项式求值以提高计算效率和数值鲁棒性。数值验证采用两个评估指标：相对L₂误差范数（用于具有解析解的问题）和相对机械能漂移（用于保守系统）。弹簧摆问题中速度非线性项的处理采用牛顿迭代法配合LDU分解与Rook主元策略。

## 研究结果

**显式积分算法的数学推导**：研究人员从二阶常微分方程的一般形式出发，对未知变量在时间点附近进行泰勒级数展开并截断至二阶导数项。为控制方法的收敛和稳定性质，最高阶导数通过参数_i进行线性插值，得到位移和速度的修正更新公式。算法实施时，方程被重组以隔离_n⁽²⁾，确保数值稳定性并最小化舍入误差。对于初始时间步，采用简化平衡方程直接计算初始二阶导数，再利用泰勒展开估计后续值。

**新显式算法的分析——稳定性**：基于单自由度无阻尼系统的自由振动方程，研究人员推导了放大矩阵（amplification matrix），并通过分析其特征值得出稳定性条件。研究发现，当参数向量=[₀,₁]=[2,3/2]时，该方法获得最优稳定性范围，即_c[0,2]，对应最大允许时间步长为=2/。特征多项式的系数由₀、₁和=决定，稳定性边界可通过不等式(16)确定。

**新显式算法的分析——一致性与精度**：通过将变量展开为截断至三阶导数的泰勒级数，研究人员确定了方法的局部截断误差。速度局部误差为_n⁽¹⁾=(3/2-₁)_n⁽³⁾²+((³))，位移局部误差为_n=(2/3-₀/2)_n⁽³⁾³+((⁴))。为确保方法至少具有二阶精度，需满足一致性条件₁=3/2。该方法因此被命名为显式积分三阶法（EIG-3），因其基于三阶导数信息构建。

**算法比较——谱特性与精度比较**：谱分析表明，对于=[2,3/2]，EIG-3在整个稳定范围内谱半径恒等于1，表现为非耗散行为，算法阻尼比=0，能理想地保持无阻尼系统的能量。此时特征值为复共轭对，确保稳定性。图4显示该配置存在轻微的周期收缩，其他参数组合则允许可控的周期伸长，为用户提供根据问题需求调谐方法的灵活性。对于=[4/3,3/2]，稳定性仅在极限情况下实现，缺乏连续的稳定范围。

**算法比较——与其他显式方法的比较**：与经典显式积分器的基准测试表明，高阶方法如RK3和AB3的稳定区域显著更小。EIG-3与CD的稳定性边界均扩展至_c[0,2]，而RK3仅限_c[0,1.732]。EIG-3与CD是仅有的在整个稳定区域内引入零人工阻尼的方法，而其他方法特别是RK2、AB2和AB3会引入数值耗散。在相对周期误差方面，EIG-3与CD具有相同的色散特性，伴随轻微周期收缩。

**关键优势——相对于中心差分法**：尽管在线性系统中EIG-3与CD具有等价的谱特性，但EIG-3具有三项关键优势：其一，对于含一阶导数非线性项的方程，标准CD方案需要迭代求解以保持二阶精度，而EIG-3对任意二阶常微分方程保持严格显式，避免迭代过程；其二，CD若强制显式处理会降阶至(() )，EIG-3则保持状态变量及其导数的二阶收敛；其三，标准CD为刚性非耗散格式，EIG-3可通过参数引入可控数值耗散，这对稳定仿真或滤除复杂动力系统中的高频噪声至关重要。

**数值验证——单自由度模型**：该模型由弹簧、质点和阻尼元件组成，运动方程为 ⁽²⁾+2+²=。取=4 rad/s，初始条件(0)=1，⁽¹⁾(0)=0，阻尼率=0，外力=0。理论分析确定最大稳定时间步长为_max=0.5 s，数值实验证实此阈值：小于等于0.5 s的步长保持稳定，超过则失稳。与CD、AB2、RK2的比较显示，EIG-3与CD产生相同精度和计算成本，AB2和RK2的稳定区域更小，在相同步长下精度较低，且RK2因每步两次函数求值而计算成本更高。

**数值验证——弹性硬化弹簧**：该系统由方程 ⁽²⁾+100(1+10²)=0 描述，初始条件(0)=1.5，⁽¹⁾(0)=0，振荡周期=0.15153 s。以100为总仿真时间，通过机械能漂移评估精度，能量表达式为=0.5(⁽¹⁾)²+50²+500⁴。结果表明，所有方法在足够小步长下精度相似，但在=110³时AB2和RK2显著失稳且误差更大。计算成本方面EIG-3与CD相近，AB2略高，RK2因两次函数求值而 noticeably更昂贵。

**数值验证——弹性软化弹簧**：由方程 ⁽²⁾+100tanh()=0 描述，初始条件(0)=4，⁽¹⁾(0)=0，周期=1.14 s。总能量为=0.5(⁽¹⁾)²+100ln[cosh()]。结果显示在=110³时所有方法稳定，但仅EIG-3和CD在=110²保持稳定。虽然所有方法均显示二阶收敛，EIG-3和CD在大步长下更精确，即使在小步长下也优于AB2和RK2。

**数值验证——弹簧摆**：该两自由度系统由质量点、无质量杆和线弹簧组成，在重力作用下运动，=9.81 m/s²。质量1 kg，杆长0.5 m，弹簧刚度98.1 N/m。运动方程包含速度非线性项，使CD变为隐式，需用牛顿法迭代求解（采用前一时间步值为初始猜测，LDU分解配合Rook主元策略）。EIG-3保持严格显式，无需解任何方程组。结果显示RK2和CD在较小步长下精度略高于EIG-3和AB2，但EIG-3在计算效率上显著优于CD——即使采用更精细的步长（如210⁶ s），EIG-3仍能以低于CD（步长110⁵ s）的成本获得更好精度。

## 讨论与总结

本研究开发了一种基于泰勒级数展开的新型二阶常微分方程显式积分器。该算法通过可调参数集实现灵活性，这种参数化框架允许用户自定义调谐稳定性和精度，提供了通常与隐式结构积分器相关的灵活性，但处于严格显式的框架内。在多种配置中，特定参数集已被识别出能产生与中心差分法相当的稳定性和收敛阶性质。该新方法每个时间步仅需一次函数求值，实现简单且计算高效。

为验证方法的理论基础，四个选自多体动力学的数值算例已被求解。EIG-3的性能与三种标准二阶显式方法（RK2、AB2和CD）进行了比较。在线性速度系统中，EIG-3与CD计算成本相似，效率优于RK2（因RK2每步需两次函数求值）。此外，EIG-3的单次求值促进了变步长，避免了龙格-库塔方法中间阶段的 overhead。EIG-3还表现出比AB2多步法更鲁棒的稳定范围，验证了使用局部高阶导数而非远处历史点的优势。

尽管所分析的方法形式上均为二阶，但精度存在差异，CD和RK2在某些线性情况下显示略好的精度。然而，在涉及速度非线性项的问题中，CD需要使用牛顿法等迭代求解器以保持精度，这增加了计算成本。相比之下，EIG-3即使面对此类非线性仍保持严格显式，从而保持其低计算开销。在这些场景中，EIG-3表现出相对于CD的优越性能。与RK2相比，EIG-3通常以大约一半的时间完成仿真，具体取决于所求值函数的复杂度。

总体而言，EIG-3作为现有二阶常微分方程显式积分方法的一种竞争性新替代方案而涌现。它提供了与CD方法相当的稳定区域，同时避免了非线性问题中对迭代求解器的需求。此外，在其稳定极限内，EIG-3不引入数值阻尼，保持解的保真度。

此处呈现的基于泰勒的推导作为新一族积分器的模块化基础。未来工作将聚焦于通过纳入泰勒展开的高阶导数项来开发该积分器的高阶版本。进一步研究将探索替代参数组合以增强稳定性和收敛性。另一个有前景的方向是将此方法扩展应用于微分-代数方程（Differential-Algebraic Equations, DAEs）的求解。