余英 卢圳

摘 要：龙格-库塔算法是数值计算非线性微分方程常用的方法之一。本文通过在计算代数系统Sagemath中建立软件包dis_rk.sage，探讨了对角隐式辛龙格-库塔算法的系数所构成的仿射簇具有的结构性质。将所得到的研究结果与已有结果进行比较，验证所设计软件包的可行性与正确性，同时也给出了在软件包下进行的数学实验所得到的新结论。

关键词：对角隐式辛龙格-库塔算法；仿射簇；积分不变量；Sagemath

中图分类号：O242.1  文献标识码：A  文章编号：1673-260X（2023）09-0031-04

1 引言

随着科技的进步，非线性现象在生活和生产领域有着广泛的应用。很明显，非线性现象较线性现象更能反映实际模型。鉴于目前许多非线性微分方程的解析解很难求得，因而学者们转向数值研究这类模型。对非线性微分方程的数值计算目前没有统一的算法，龙格-库塔算法（Runge-Kutta method）是非线性方程数值计算最常用的一种方法。在对非线性方程的数值模拟中，由于隐式的龙格-库塔算法在每一步的计算中，都要求一个非线性方程组的解。这不仅增加了计算成本，而且求非线性方程组的解至今没有一致且有效的算法[1]，因而给计算带来不便。在这一点上，显式的龙格库塔法就比较受欢迎。它被设计成专门的程序，并在各种计算代数系统中都有实现它的软件包。

近三十年来，构造守恒的算法受到了广大学者的高度重视。守恒的要求在物理现象、化学现象、天体运动和其他自然现象中处处可见。比如，C. Marchal指出三体问题对应的系统具有10个经典的守恒量，即：总能量不变量（the invariant of the energy）、动量不变量（the invariant of impulse moment，它包括三个标量不变量）、角动量不变量（the invariant of angular moment，它包括三个标量不变量）和质量中心不变量（the invariant of the center of mass，它包括三个标量不变量）[2]。又比如钱旭指出变系数的一维非线性Schr？觟dinger方程的解满足电荷守恒、全局动量守恒和全局能量守恒[3]。数值解要满足守恒定律比一般的过程更重要有许多理由。如果得到的数值解不具有物理意义，如能量守恒，是我们不期望的。换句话说，可以认为这样的数值解存在巨大的误差。早在1928年，Courant， Friedrichs和lewy就在构造算法时渗透了保持能量守恒的思想，只不过当时的这个保持的量是一个正定的量，还没有能量守恒这种提法。在国内，据了解，学者郭本瑜是比较早构造具有这种结构算法的学者之一[4-6]。再后来，冯康先生提出了辛几何算法思想[7]，并利用生成函数与辛方法的联系，证明了任意的高阶辛算法是存在的这个结论。该思想在国内引起了极大的反响，直接推动了保结构算法的研究。构造保结构算法的原则是问题的模型在离散化后其基本特征应当得到最大限度地保持[8-11]。随后，保结构算法的构造成为世界前沿的问题之一。近年来，保结构算法已在量子力学、流体力学、结构力学、水动力学，电磁学、地球物理学等科学和工程领域得到了广泛的应用。

然而，遗憾的是E. Hairer等人指出显式的龙格-库塔法对高于一次的多项式不变量不守恒[12]。J. M. Sanz-Serna指出辛龙格-库塔算法对二次量守恒的结论[13]。近年，H. Zhang等提出了能量不变量二次化法（IEQ method），该方法通过引入辅助变量，将高于二次的多项式能量不变量转化为二次守恒量[14]。因而，设计对二次不变量的守恒算法具有重要的研究意义。Y. Yu研究了一般的辛隐式龙格-库塔算法[15]，对其设计了软件包rk.sage，借助于该软件包，给出了一些相关的研究结果。考虑到一般的隐式辛龙格-库塔算法在实际的数值运算过程中的计算复杂性。隐式龙格-库塔算法中一类较简单的是对角隐式辛龙格-库塔算法。由于它除了左下三角和对角线元素不为零外，其余都为0，就大大减少了大型计算过程中的计算复杂性，因而应用该类算法会大大地节约计算成本。

本文在计算代数系统Sagemath中建立软件包dis_rk.sage。借助于该软件包，研究了对角隐式辛龙格-库塔算法的系数所构成的仿射簇（Affine Variety）的性质。将所得到的研究结果与已有结果进行比较，验证所设计软件包的可行性与正确性。同时，借助于该软件包，还发现了一些新结果。

2 预备知識

本文考虑以下自治的微分系统

其中t是自变量，通常代表时间；x是n维向量 （x1，…xn），通常代表坐标；f是一个向量函数（f1，…，fn），其中fi是坐标x的函数。

定义1 如果存在一个关于x的度为2的多项式I（x），使得条件

满足，则称I（x）是关于自治系统（1）的二次积分不变量。

对于动态系统（1）s级数的龙格-库塔算法写为

以及

定义2 如果龙格-库塔的系数满足

biaij+bjaji-bibj=0，i，j=1，…，s，

则称其为辛的。当对动态系统（1）数值积分时，只有辛龙格-库塔算法才能对二次积分不变量守恒[13]。

定义3 如果由式（2）、（3）迭代得到的近似值与精确值的前p项相等，则称该系数对应的龙格-库塔算法是p阶近似的，即截断误差为O（hp+1）[8]。

以定义2，定义3确定的对龙格-库塔系数的限制条件作为多项式环中的方程组，得到了由龙格-库塔算法的系数构建的多项式环。用类似于Y. Yu的研究方法[15]，在计算代数系统Sagemath中建立软件包dis_rk.sage，通过多项式环的信息来研究龙格-库塔算法的系数的性质，由此得到龙格-库塔算法的各种格式。

3 算法设计

在本部分，给出计算S级p阶近似的对角隐式辛龙格-库塔算法的系数满足的方程组的算法。下令dt=t-t0，t0为初值起点。

步骤1：对于问题（1），计算满足初值条件x（t0）=0的精确解x（t）在点t=t0处的泰勒展开式。它是关于函数f的任意阶导数以及f′（t0），f"（t0），……的多项式表达式，将该表达式记为

g1（f′（t0），[f′（t0）]2，……，f′（t0），[f"（t0）]2，…，f（p-1）（t0），[f（p-1）（t0）]2， …，dt），简记为g1;

步骤2：对于问题（1），通过将S级p阶近似的对角隐式辛龙格-库塔算法的步骤一表达式进行泰勒展开后代入步骤二，计算满足初值条件x（t0）=0的近似解的表达式。它是关于函数f的任意阶导数f′（t0），f"（t0），……，a11，a12，……，ass，b1，……，bs以及dt的多项式表达式，将该表达式记为

g2（f′（t0），[f′（t0）]2，……，f（p-1）（t0），[f（p-1）（t0）]2，……，dt，a11，a12，……，ass，b1，……，bs），简记为g2;

步骤3：若上述算法的精确度为p阶，则表达式g1与表达式g2的前p项应该一致，即dt，d2t，……，d（p-1）t处前的系数应该相等，由此得到p个关于龙格-库塔系数a11，a12，……，ass，b1，……，bs以及函数f的任意阶导数f′（t0），f"（t0），……和它的多项式的等式。再由f的任意性可知，单项式，f′（t0），f"（t0），……的系数必须为0，得到若干个关于龙格-库塔系数a11，a12，……，ass，b1，……，bs的等式组，从而形成一个Q[a11，a12，……，ass，b1，……，bs]的一个理想[15]。

4 软件包应用举例

在仿射空间（a11，a12，……，ass，b1，……，bs）中描述仿射簇，该流体上的点的坐标可以看作是一个满足辛要求和s级数且具有p阶近似的对角隐式辛龙格库系数的一组取值。以下将该仿射簇记为I（s，p）。软件包dis_rk.sage中函数rk_var_dis（s）返回s级数且具有p阶近似的对角隐式辛龙格库系数的所有需要的符号变量和环Kab=Q[a11，a12，……，ass，b1，……，bs]。按照算法1，p阶近似的条件产生了一系列的方程，函数rk_series（s，p）返回这些方程的左边。让我们来看一个I（2，2）例子：

sage ： var （’x，t，dt ’）

（x ， t ， dt ）

sage ： load （’sage /dis_rk. sage ’）

None

sage ： rk_var_dis （2）

None

sage ： a

[a00 0]

[a10 a11]

sage ： b

（b0 ， b1）

sage ： c

[a00 ， a10 + a11]

sage ： eqs = rk_symplectic （2，2）

sage ： eqs

[2* a00 * b0 - b0 ^2 ， a10 * b1 - b0 *b1 ， a10 * b1 - b0 * b1，

2* a11 * b1 - b1 ^2 ， - b0 - b1 + 1 ， -2* a00 * b0 - 2* a10 *

b1 - 2* a11 * b1 + 1]

eqs返回的是理想J在仿射空间（a11，a12，……，ass，b1，……，bs）中生成的仿射簇I（2，2）

sage ： J = Kab * eqs

sage ： J.is_prime （）

False

sage ： eq = [J . primary_decomposition （） [ i ]. gens （） for i in range （len （ J.primary_

decomposition （） ））]

sage ： eq

[[ b1， b0 - 1， 2* a00 - 1]， [ b1 - 1， b0， 2* a11 - 1， a10 ]， [ b0 + b1 - 1， 2* a11 - b1， a10 + b1 - 1， 2* a00 + b1 - 1]]

sage ： Jp = Kab . ideal （ eq [2]）

sage ： Jp

Ideal （ b0 + b1 - 1， 2* a11 - b1， a10 + b1 - 1， 2* a00 + b1 - 1） of Multivariate Polynomial Ring in a00， a10， a11， b0， b1 over Rational Field

sage ： Jp.dimension （）

1

sage ： Jp.genus （）

0

上述數学实验结果dim（I（2，2））=1且genus（I（2，2））=0。

5 仿射簇I（s，p）的性质

在这一部分，借助于软件包dis_rk.sage，对对角隐式辛龙格-库塔法的系数所构成的仿射簇展开研究，给出其一些性质。由这些系数得到的龙格-库塔算法的格式不仅具有计算方便的优势，而且也对二次积分不变量守恒。不失一般性，我们假设bi≠0（i=1，…，s）。若bk=0，则辛条件导致了biaik=0（i=1，…，s），则该方法等价于低阶的对应算法。通过在计算代数系统Sagemath中的数学实验，得到了以下关于仿射簇I（s，p）的维数（dimension）的结论：

维度（I（2，4））=维度（I（2，5））=维度（I（2，6））=维度（I（3，5））=维度（I（3，6））=维度（I（4，5））==维度（I（4，6））=维度（I（5，6））=-1；

维度（I（2，3））=维度（I（3，4））=维度（I（5，0））=0；

维度（I（2，2））=维度（I（3，3））=维度（I（4，4））=1；

维度（I（3，2））=维度（I（4，3））=维度（I（5，4））=2；

维度（I（4，2））=维度（I（5，3））=3；

维度（I（5，2））=4。

其中0维代表该集是有限集，-1维代表该集是空集。當该理想的维数为1时，借助于软件包可以计算其相应的亏格（genus），借助数学实验，得到以下结论：亏格（I（3，3））=1，亏格（I（2，2））=0。下面对5种情形的I（s，p）结构展开具体的研究。

（1）情形I（2，2）

数学实验I（2，2）对应的仿射簇所代表的曲线的维数为1，亏格为0，由代数几何知识可知仿射簇I（2，2）存在一个有理的参数化方程，即：

（2）情形I（2，3）

由上述结论可知，仿射簇I（2，3）是由有限个点组成的，其相应的龙格—库塔系数为：

注记2：I（2，3）不存在实系数元素，也就是说，不存在2级3阶的实系数对角辛隐式龙格库塔算法。

（3）情形I（3，3）

仿射簇I（3，3）对应的龙格—库塔系数为

其中b0，b1，b2满足

3b12b2-3b12+3b22b1-6b1b2+3b1-3b22+3b2-1=0b0+b1+b2=1

若令b1=b0，计算代数系统Sagemath中的数学实验给出了以下结论：

该结论与文献[8]中第279页的结论一致。其中a=1.351207，它是多项式6x3-12x2+6x-1=0的一个实根。

（4）情形I（3，4）

仿射簇I（3，4）是由有限个点组成的，计算代数系统Sagemath中的数学实验显示它包含以下两类对角隐式辛龙格-库塔算法。

类型1：

注记3：数学实验显示I（3，4）只含有一个实系数的点，也就是说，只有一个相应的实系数对角隐式辛龙格库塔算法。

（5）情形I（5，5）

很明显，I（5，5）包含有限个点，因此相应的对角隐式辛龙格库塔算法是有限的。然而，数学实验告诉我们I（5，5）不存在实系数的对角隐式辛龙格库塔算法。总的来说，正如G. Sun所指出的那样高于4阶的实系数的对角辛龙格库塔算法是不存在的[18]。

6 结论

本文首先在计算代数系统Sagemath中建立了软件包dis_rk.sage。基于该软件包，对由s级数且具有p阶近似的对角隐式辛龙格库塔法的系数所构成的仿射簇I（s，p）在该数学软件中进行了一系列的数学实验。数学实验证明我们的软件包是可行且正确有效的。对于I（s，p）（s，p≤4）给出了相应的的结构性质。而且，数学实验结果显示高于4阶的实系数的对角辛龙格库塔算法不存在。

——————————

参考文献：

〔1〕H. Press， B. P. Flannery， S. A. Teukolsky， W. T. Vetterling， et al.. Numerical recipes[M]. Cambridge University Press Cambridge， 1989.

〔2〕C.Marchal， The Three-Body Problem[M]. Elsevier Science publishing company，1990.

〔3〕钱旭.几类偏微分方程的保结构算法研究[D].北京：国防科学技术大学，2014.

〔4〕郭本瑜.偏微分方程的差分方法[M].北京：科学出版社，1988.

〔5〕Guo B Y，Vázquez L. A numerical scheme for nonlinear Klein-Gordon equation. J Appl Sci，1983，1（01）：25-32.

〔6〕Guo B Y， Pascual P J. Numerical solution of the Sine-Gordon equation. Appl Math Comput， 1986：18：1-14.

〔7〕Feng K.On difference schemes and symplectic geometry[C]proceedings of the 1984 Beijing symposium on Differential Geometry and Differential Equations. Beijing：Science Press，1985，01：42-58.

〔8〕冯康，秦孟兆.哈密尔顿系统的辛几何算法[M].杭州：浙江科学技术出版社，2003.

〔9〕Feng K， Qin M. Symplectic Geometric Algorithms for Hamiltonian Systems[M]. Berlin： Springer， 2010.

〔10〕Feng K， Wu H， Qin M Z， et al. Construction of canonical difference schemes for Hamiltonian formalism via generating functions[J]. J. Comput. Math.， 1989， 7：71-96.

〔11〕秦孟兆，王雨順.偏微分方程中的保结构算法[M].杭州：浙江科学技术出版社，2011.

〔12〕E. Hairer， G. Wanner， and S. P. N？覬rsett， Solving ordinary differential equations I： Nonstiff Problems[M]. 3rd ed. New York： Springer， 1993.

〔13〕J. M. Sanz-Serna. Symplectic Runge–Kutta schemes for adjoint equations， automatic differentiation， optimal control， and more[J]. SIAM REVIEW. 2016，58 （01）： 3–33.

〔14〕H. Zhang， X. Qian， J. Yan， and S. Song. Highly efficient invariant- conserving explicit Runge–Kutta schemes for nonlinear Hamiltonian differential equations[J]. Journal of Computational Physics， 2020，48 ：1-16.

〔15〕Y.Yu.The symbolic problems associated with Runge-Kutta methods and their solving in Sage，Discrete and ontinuous models and applied computational science4[J].2019，27（01）：33-41.

〔16〕J. C. Butcher， On Runge-Kutta processes of high order， Journal of the Australian Mathematical Society [J].1964， 4 （02）： 179-194.

〔17〕G. Cooper.Stability of Runge–Kutta methods for trajectory problems[J]. IMA journal of numerical analysis. 1987，7（01）： 1-13.

〔18〕G．Sun. A simple way constructing symplectic Runge-Kutta methods[J]. J. Comput. Math. ，2000（18）：61-68.