基于符號計算的對角隱式辛龍格-庫塔算法結構研究

余英 盧圳

摘 要：龍格-庫塔算法是數值計算非線性微分方程常用的方法之一。本文通過在計算代數系統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）不存在實系數的對角隱式辛龍格庫塔算法?？偟膩碚f，正如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.