一种用于简易Lawson积分的显式六阶Runge–Kutta方法

《Journal of Computational Mathematics and Data Science》：An explicit sixth order Runge–Kutta method for simple Lawson integration

　　

当动力学系统中存在刚性线性算子时，显式Runge–Kutta（RK）格式会变得不实用，这种失效模式在流体与等离子体数值模拟中十分常见。Lawson于1967年提出针对刚性问题的广义Runge–Kutta过程（Generalized Runge–Kutta Processes），其中刚性线性算子通过矩阵指数被完全隐式处理。任意RK格式都可诱导出合法的Lawson积分，但若节点横坐标ci有序且等间距，则该格式实现异常简单。经典RK4满足此要求，但在该约束下推导高效高阶格式十分困难。研究人员通过Newton–Raphson迭代识别出一种显式六阶方法，可提供简易Lawson积分。

研究背景方面，自主常微分方程组? = g(u) + ?u广泛来源于流体动力学与等离子体物理中偏微分方程的空间离散，其中?为具有大负实部特征值的线性算子，导致刚性问题。现有处理途径包括隐式–显式（IMEX）RK格式与指数Runge–Kutta（ExpRK）格式，前者需反复求解(I? ? a_iih?)线性系统，后者需谨慎计算大量解析函数。Lawson广义RK过程将刚性部分用exp(c_ih?)精确演化，但任意RK格式套用需计算O(s)个矩阵指数；Lawson从RK4实验中发现，若c_i在[0,1]有序等距且c₁=0，则整个算法仅需一个矩阵指数exp(Δc h?)，研究人员将此特例命名为简易LawsonRunge–Kutta（SLRK）。三阶段以上RK中Heun三阶Δc=1/3、RK4 Δc=1/2均满足，但已有文献中大量六阶RK（Huta、Butcher、Fehlberg、Prince–Dormand、Verner等）均无有序等距c_i。因此研究人员旨在构造满足该约束的六阶显式RK。

研究人员采用的关键技术方法为：以七阶段试算Δc=1/6与1/5做数万次随机初值Newton–Raphson迭代（x为b_i⊕a_ij系数集，F为37个六阶条件加c_i偏离等距网格的残差，J⁺为Jacobian伪逆，γ∈(0,1]为阻尼），七阶段无收敛，八阶段迅速得根；将双精度根映射为小分母有理数并符号验证37条件全满足，得八阶段六阶RK6-uniform（c_i=0,1/6,…,1）。配套SLRK6仅穿插exp(h?/6)。测试样本含R?ssler方程（初值[1,2,3]积到t=5）、二维周期域1024×1024谱方法解不可压Navier–Stokes（ν=10^?2，?=ν?²，强迫sin4y?）、van der Pol（μ=100，初值[2,0]积到t=200，?取初值Jacobian）。

研究结果部分，小标题2. Main results中，研究人员给出RK6-uniform的Butcher表（八阶段六阶，Δc=1/6），在R?ssler流误差对比中其误差低于Prince–Dormand RK6(5)^?与Verner RK6(5)9b^?等多数文献六阶格式；由该表诱导的SLRK6在Navier–Stokes涡量演化和van der Pol振荡器中均达六阶收敛，van der Pol在μΔt～1处出现12阶超收敛再回落六阶，且稳定步长受μΔt?1限制。小标题稳定性函数中，RK4与RK6-uniform的稳定性多项式Φ_RK4(z)、Φ_RK6(z)分别截至z⁴、z⁶加406/495·z⁷/7!项；SLRK4/SLRK6的二变量稳定性函数Φ(z₁,z₂)=exp(z₂)·Φ_RK(z₁)，对Re(z₂)?0刚度情形稳定域远大于标准RK，但因最高次z₁^p（p=4或7）使SLRK4在强负z₂下比SLRK6稳定域更宽，此为Lawson格式普遍特征。

讨论与结论部分，研究人员指出Lawson积分借显式RK适应刚性系统，c_i有序等距令实现仅需单个矩阵指数。RK6-uniform是首例八阶段有理系数、有序等距显式六阶RK，其自身及诱导SLRK6在R?ssler流、不可压Navier–Stokes与van der Pol振荡器上经符号代数严格校验37阶条件与六阶收敛。该格式由Newton–Raphson在阶条件与等距约束上搜得，表明现代计算可继续用此法探更高阶SLRK友好格式。论文发表于《Journal of Computational Mathematics and Data Science》。结论是：研究人员提出RK6-uniform并证实其提供简易Lawson六阶积分能力，刚性耗散算子下数值稳定域显著扩大，为流体与等离子体高维刚性ODE/PDE求解提供低实现成本高阶工具。