用于求解刚性微分方程的、时间精确度高且稳定性强的显式两步龙格-库塔方法

《Journal of Computational and Applied Mathematics》：Time-accurate and highly-stable explicit two-step Runge-Kutta methods for stiff differential equations

　　**摘要**
在本文中，我们引入了一种时间精确且高度稳定的显式两步Runge-Kutta方法（TASE-TSRK），用于求解刚性常微分方程。这些方法的灵感来源于Calvo等人[1]提出的技术。我们的目标是通过预先用适当的TASE算子乘以向量场来提高显式TSRK方法的稳定性。在这项工作中，我们分析了相关的阶条件，并对阶数k≤4的TASE-TSRK方法进行了全面的稳定性研究。通过适当选择TASE算子族的自由参数，我们证明了二阶方法是强A稳定的，而三阶和四阶方法几乎可以达到强A稳定。一些数值实验表明，新的TASE-TSRK方案与TSRK方法相比具有竞争力，并且对于刚性问题可能更有优势。

**引言**
我们研究了刚性常微分方程（ODEs）初值问题（IVPs）的数值解。最初假设这些方程是自治形式的：{y′(t)=f(y(t)), t∈[t0,T], f:Rd→Rd, y(t0)=y0}。这些方程也可以来自一维或多维中抛物型偏微分方程（PDEs）的空间离散化。为了高效准确地求解刚性IVPs，数值方法必须具有足够的稳定性。著名的例子包括基于BDF的多步方案[1][2]、隐式Runge-Kutta族（如DIRK和SIRK[3]）或具有较少内部阶段的高度隐式RK公式（例如三阶段Radau IIA）、通用线性方法[4]，以及线性隐式方法（如Rosenbrock和W方法[5]）。更多有趣的数值技术，请参考[6][7][8][9]。
最近，Bassenne、Fu和Mani[10]提出了一类基于显式方法的时间推进方案，用于求解刚性IVPs。在这种方法中，原始的刚性IVP（1）通过用适当的线性算子Tp∈Rd×d预乘向量场f来进行修改，这种算子被称为时间精确且高度稳定的显式（TASE）算子。然后，修改后的问题{u′(t)=Tpf(u(t)), t∈[t0,T], u(t0)=y0}通过显式的s阶Runge-Kutta（RK）方法来求解。TASE算子Tp的选择取决于步长h>0和向量场f，必须确保得到的数值方案至少具有与显式RK方法相同的精度阶数，并具有良好的稳定性。注意，如果Tp是单位算子Id的p阶近似，即Tp=Id+O(hp)，根据常数变分公式，(1)和(2)的精确解满足∥y(t)?u(t)∥=O(hp)。同时，如果显式RK方法的阶数为p，并且TASE算子满足Tq=Id+O(hq)，可以证明稳定化系统(2)的数值解至少具有min(p, q)的阶数（见[11]）。在[10]中提出的TASE算子Tp=Tp(α,hJ)是递归定义的，即Tp(hJ;α)={(Id?αhJ)?1}^{p?1}=(2p?1)^{p?1?1}Tp?1(hJ;α)^{2p?1?1}Tp?1(hJ;α)^{2p?1?1}，其中J=fy(t,y(t))表示时间t处f的雅可比矩阵，α是一个可以适当选择的正自由参数，以获得具有良好稳定性的TASE-RK方法。在[10]中已经证明，当用阶数p≤s的TASE算子进行预处理时，s阶sth阶显式RK方法是A稳定的。
随后，Calvo、Montijano和Rández[1]引入了另一种类型的TASE算子Tp(hJ;α)=∑j=1pγj(Id?αjhJ)?1，其中α∈Rp, αj>0, αj≠αk（对于j≠k），并且γj=(βj)^{p?1}∏k≠j(βj?βk)，βj=1α^jj=1,…,p。公式(6)在零附近的泰勒展开得到Tp=I+KphpJ^{p/2}＋O(h^(p/2+)，其中Kp=∑j=1pγjαjp。这类TASE算子包括p个自由实数参数αj（j=1,…,p），可以选择这些参数以满足稳定性和精度要求。对阶数p≤4的p阶显式RK方案在TASE算子预处理下的A稳定性进行了全面研究。基于[1][10]中TASE-RK方法与W方法类的联系，González-Pinto等人[12]提出了一类新的所谓TASE-W方法，这些方法每个积分步骤所需的线性系统数量减少，且具有一致性阶数和良好的稳定性。Conte等人[13]研究了将近似矩阵分解方法与TASE W方法结合用于抛物型PDEs的时间积分。最近，Acetoa、Conte和Pagano[14]通过允许向量α中出现复数参数对来推广了(5)中定义的TASE算子，从而在不增加计算成本的情况下提高了相应TASE-RK方法的精度和稳定性。
为了降低RK-TASE方法的计算成本，Calvo、Fu、Montijano和Rández[11]开发了一类不同的单TASE算子STp(hJ;α)=∑j=1pβj(Id?αhJ)?^j，其中α>0，实数参数βj（j=1,…,p）作为α的函数唯一确定，以满足STp(hJ;α)=Id+O(hp)。对阶数p≤4的p阶显式RK方案与相同阶数的单TASE算子结合的线性稳定性进行了全面研究。此后，Calvo、Montijano和Rández[15][16]提出了对Runge-Kutta方法每个阶段使用不同TASE算子的改进，以提高效率和精度。Wang和Cong[17]将TASE-RK和TASE-W方法扩展到求解中性延迟微分方程，并研究了TASE-RK和TASE-W与拉格朗日插值结合时的延迟依赖稳定性。
值得一提的是，Conte、Pagano和Paternoster[18]以及Pagano[19]通过将TASE技术应用于显式两步同伴方法，推导出了一类具有良好稳定性的线性隐式并行化数值方案。同伴方法的引入旨在结合多步和Runge-Kutta方法的优点（见[20][21]）。由于结构特性，对于某些系数选择，同伴方法可以并行化。此外，所有阶段都具有相同的精度和稳定性，这使得得到的数值方案在应用于高度刚性微分问题时能够避免阶数降低。结合TASE算子的优势，它们成为半离散化PDEs产生的大型ODE系统数值积分的有吸引力的候选方法。
在本文中，我们旨在采用TASE技术设计一类新的TASE两步方法，与同阶的TSRK和TASE-RK方法相比，这些方法具有更低的计算成本和更好的稳定性。众所周知，两步Runge-Kutta方法通过使用当前和前一积分步骤获得的变量和导数值来推广一步方法。这使得能够构建具有传统Runge-Kutta方法无法实现特性的方法。例如，显式两步Runge-Kutta（TSRK）方法的一个关键优势是，至少对于p≤5的情况，p阶的TSRK方法需要与p?1阶的RK方法相同数量的场函数f的评估次数（见[22]）。TSRK方法的另一个优势是它们的阶段阶数可以高于显式一步RK方法，从而确保内部阶段的更高精度，从而提高中等刚性和高度非线性问题的解的精度。当与TASE算子结合时，这些属性允许TASE-TSRK方法在相对较低的计算成本下实现更高的精度或更好的稳定性。
本文的组织结构如下：第2节中，我们构建了TASE-TSRK方法并得到了阶条件。第3节中，我们展示了适当选择TASE算子（5）的自由系数，使得一阶二阶、两阶三阶和三阶四阶TASE-TSRK方法具有良好的稳定性。第4节进行了三个数值实验，以证明我们提出方法的有效性，强调了它们在求解线性和非线性刚性系统时的潜力及改进的性能。最后，在第5节中讨论了结论性意见。

**TASE-TSRK方法**
对于常微分方程（ODE）(1)，s阶显式TSRK方法定义为：
{Yi[n?1]=yn?2+h∑j=1i?1aijf(Yj[n?1]), Yi[n]=yn?1+h∑j=1i?1aijf(Yj[n]), i=1,…,s,
yn=(1??)yn?1+?yn?2+h∑j=1s(vjf(Yj[n?1])+wjf(Yj[n])),
其中h是常数步长，yn是对y(tn)的近似，Yi[n]是对y(tn?1+cih)的近似，tn=t0+nh（n=0,1,…,N，Nh=T?t0），aij、?、vj、wj是方法的系数。为了方便表示（9），可以使用以下系数表：
其中节点ci=∑j=1i?1aij（i=1,…,s）。

**TASE-TSRK方法的稳定性**
在本节中，我们分析了使用(5)中定义的Tp的TASE-TSRK方法的稳定性。将TASE-TSRK方法应用于线性测试方程y′=λy（λ∈C），我们得到：
y[n]=M(z)y[n?1]，
其中M(z)=[1??+zTp(z;α)wTd?vT100zTp(z;α)d00]，
v=[v1?vs]，w=[w1?ws]，e=[1?1]∈Rs，
d=(I?zTp(z;α)A)^?1e。由于A是一个严格下三角矩阵，d的每个分量都是z的多项式。

**数值示例**
在本节中，我们通过求解一些刚性问题来说明新TASE-TSRK方法的性能。为了比较，我们使用了以下数值方法：
• TSRK2：二阶单步显式TSRK方法，ρ=0.1；
• TSRK3：三阶两步显式TSRK方法，ρ=?0.1，c2=1；
• TSRK4：第四阶三步显式TSRK方法（案例3），ρ=?0.625，v3=0.2；
• TASE2-TSRK2：二阶单步TASE2-TSRK方法（26），ρ=0.1，α1=2，α2=8；
• TASE3-TSRK3：三阶两步TASE3-TSRK方法。

**结论**
在本文中，我们提出了一种新的数值方案，通过结合显式TSRK方法和[1]中介绍的TASE技术来求解刚性ODEs。研究了阶数k≤4的提出的方法的阶条件和稳定性。我们提出了一些强A稳定的二阶TASE-TSRK方案，并证明了三阶和四阶TASE-TSRK方法不能是A稳定的，但可以是近似A稳定的。进行了数值实验来展示它们的稳定性。

**未引用的参考文献**
[26] 5Hongjiong Tian, Junyu Chen, Zhongli Liu
上海师范大学数学系，中国上海桂林路100号，200234