基于符號計算的對角隱式辛龍格-庫塔算法結(jié)構(gòu)研究

余英 盧圳

摘 要：龍格-庫塔算法是數(shù)值計算非線性微分方程常用的方法之一。本文通過在計算代數(shù)系統(tǒng)Sagemath中建立軟件包dis_rk.sage，探討了對角隱式辛龍格-庫塔算法的系數(shù)所構(gòu)成的仿射簇具有的結(jié)構(gòu)性質(zhì)。將所得到的研究結(jié)果與已有結(jié)果進(jìn)行比較，驗證所設(shè)計軟件包的可行性與正確性，同時也給出了在軟件包下進(jìn)行的數(shù)學(xué)實驗所得到的新結(jié)論。

關(guān)鍵詞：對角隱式辛龍格-庫塔算法；仿射簇；積分不變量；Sagemath

中圖分類號：O242.1? 文獻(xiàn)標(biāo)識碼：A? 文章編號：1673-260X（2023）09-0031-04

1 引言

隨著科技的進(jìn)步，非線性現(xiàn)象在生活和生產(chǎn)領(lǐng)域有著廣泛的應(yīng)用。很明顯，非線性現(xiàn)象較線性現(xiàn)象更能反映實際模型。鑒于目前許多非線性微分方程的解析解很難求得，因而學(xué)者們轉(zhuǎn)向數(shù)值研究這類模型。對非線性微分方程的數(shù)值計算目前沒有統(tǒng)一的算法，龍格-庫塔算法（Runge-Kutta method）是非線性方程數(shù)值計算最常用的一種方法。在對非線性方程的數(shù)值模擬中，由于隱式的龍格-庫塔算法在每一步的計算中，都要求一個非線性方程組的解。這不僅增加了計算成本，而且求非線性方程組的解至今沒有一致且有效的算法[1]，因而給計算帶來不便。在這一點上，顯式的龍格庫塔法就比較受歡迎。它被設(shè)計成專門的程序，并在各種計算代數(shù)系統(tǒng)中都有實現(xiàn)它的軟件包。

近三十年來，構(gòu)造守恒的算法受到了廣大學(xué)者的高度重視。守恒的要求在物理現(xiàn)象、化學(xué)現(xiàn)象、天體運動和其他自然現(xiàn)象中處處可見。比如，C. Marchal指出三體問題對應(yīng)的系統(tǒng)具有10個經(jīng)典的守恒量，即：總能量不變量（the invariant of the energy）、動量不變量（the invariant of impulse moment，它包括三個標(biāo)量不變量）、角動量不變量（the invariant of angular moment，它包括三個標(biāo)量不變量）和質(zhì)量中心不變量（the invariant of the center of mass，它包括三個標(biāo)量不變量）[2]。又比如錢旭指出變系數(shù)的一維非線性Schr？觟dinger方程的解滿足電荷守恒、全局動量守恒和全局能量守恒[3]。數(shù)值解要滿足守恒定律比一般的過程更重要有許多理由。如果得到的數(shù)值解不具有物理意義，如能量守恒，是我們不期望的。換句話說，可以認(rèn)為這樣的數(shù)值解存在巨大的誤差。早在1928年，Courant， Friedrichs和lewy就在構(gòu)造算法時滲透了保持能量守恒的思想，只不過當(dāng)時的這個保持的量是一個正定的量，還沒有能量守恒這種提法。在國內(nèi)，據(jù)了解，學(xué)者郭本瑜是比較早構(gòu)造具有這種結(jié)構(gòu)算法的學(xué)者之一[4-6]。再后來，馮康先生提出了辛幾何算法思想[7]，并利用生成函數(shù)與辛方法的聯(lián)系，證明了任意的高階辛算法是存在的這個結(jié)論。該思想在國內(nèi)引起了極大的反響，直接推動了保結(jié)構(gòu)算法的研究。構(gòu)造保結(jié)構(gòu)算法的原則是問題的模型在離散化后其基本特征應(yīng)當(dāng)?shù)玫阶畲笙薅鹊乇３諿8-11]。隨后，保結(jié)構(gòu)算法的構(gòu)造成為世界前沿的問題之一。近年來，保結(jié)構(gòu)算法已在量子力學(xué)、流體力學(xué)、結(jié)構(gòu)力學(xué)、水動力學(xué)，電磁學(xué)、地球物理學(xué)等科學(xué)和工程領(lǐng)域得到了廣泛的應(yīng)用。

然而，遺憾的是E. Hairer等人指出顯式的龍格-庫塔法對高于一次的多項式不變量不守恒[12]。J. M. Sanz-Serna指出辛龍格-庫塔算法對二次量守恒的結(jié)論[13]。近年，H. Zhang等提出了能量不變量二次化法（IEQ method），該方法通過引入輔助變量，將高于二次的多項式能量不變量轉(zhuǎn)化為二次守恒量[14]。因而，設(shè)計對二次不變量的守恒算法具有重要的研究意義。Y. Yu研究了一般的辛隱式龍格-庫塔算法[15]，對其設(shè)計了軟件包rk.sage，借助于該軟件包，給出了一些相關(guān)的研究結(jié)果。考慮到一般的隱式辛龍格-庫塔算法在實際的數(shù)值運算過程中的計算復(fù)雜性。隱式龍格-庫塔算法中一類較簡單的是對角隱式辛龍格-庫塔算法。由于它除了左下三角和對角線元素不為零外，其余都為0，就大大減少了大型計算過程中的計算復(fù)雜性，因而應(yīng)用該類算法會大大地節(jié)約計算成本。

本文在計算代數(shù)系統(tǒng)Sagemath中建立軟件包dis_rk.sage。借助于該軟件包，研究了對角隱式辛龍格-庫塔算法的系數(shù)所構(gòu)成的仿射簇（Affine Variety）的性質(zhì)。將所得到的研究結(jié)果與已有結(jié)果進(jìn)行比較，驗證所設(shè)計軟件包的可行性與正確性。同時，借助于該軟件包，還發(fā)現(xiàn)了一些新結(jié)果。

2 預(yù)備知識

本文考慮以下自治的微分系統(tǒng)

其中t是自變量，通常代表時間；x是n維向量 （x1，…xn），通常代表坐標(biāo)；f是一個向量函數(shù)（f1，…，fn），其中fi是坐標(biāo)x的函數(shù)。

定義1 如果存在一個關(guān)于x的度為2的多項式I（x），使得條件

滿足，則稱I（x）是關(guān)于自治系統(tǒng)（1）的二次積分不變量。

對于動態(tài)系統(tǒng)（1）s級數(shù)的龍格-庫塔算法寫為

以及

定義2 如果龍格-庫塔的系數(shù)滿足

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

則稱其為辛的。當(dāng)對動態(tài)系統(tǒng)（1）數(shù)值積分時，只有辛龍格-庫塔算法才能對二次積分不變量守恒[13]。

定義3 如果由式（2）、（3）迭代得到的近似值與精確值的前p項相等，則稱該系數(shù)對應(yīng)的龍格-庫塔算法是p階近似的，即截斷誤差為O（hp+1）[8]。

以定義2，定義3確定的對龍格-庫塔系數(shù)的限制條件作為多項式環(huán)中的方程組，得到了由龍格-庫塔算法的系數(shù)構(gòu)建的多項式環(huán)。用類似于Y. Yu的研究方法[15]，在計算代數(shù)系統(tǒng)Sagemath中建立軟件包dis_rk.sage，通過多項式環(huán)的信息來研究龍格-庫塔算法的系數(shù)的性質(zhì)，由此得到龍格-庫塔算法的各種格式。

3 算法設(shè)計

在本部分，給出計算S級p階近似的對角隱式辛龍格-庫塔算法的系數(shù)滿足的方程組的算法。下令dt=t-t0，t0為初值起點。

步驟1：對于問題（1），計算滿足初值條件x（t0）=0的精確解x（t）在點t=t0處的泰勒展開式。它是關(guān)于函數(shù)f的任意階導(dǎo)數(shù)以及f′（t0），f"（t0），……的多項式表達(dá)式，將該表達(dá)式記為

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階近似的對角隱式辛龍格-庫塔算法的步驟一表達(dá)式進(jìn)行泰勒展開后代入步驟二，計算滿足初值條件x（t0）=0的近似解的表達(dá)式。它是關(guān)于函數(shù)f的任意階導(dǎo)數(shù)f′（t0），f"（t0），……，a11，a12，……，ass，b1，……，bs以及dt的多項式表達(dá)式，將該表達(dá)式記為

g2（f′（t0），[f′（t0）]2，……，f（p-1）（t0），[f（p-1）（t0）]2，……，dt，a11，a12，……，ass，b1，……，bs），簡記為g2;

步驟3：若上述算法的精確度為p階，則表達(dá)式g1與表達(dá)式g2的前p項應(yīng)該一致，即dt，d2t，……，d（p-1）t處前的系數(shù)應(yīng)該相等，由此得到p個關(guān)于龍格-庫塔系數(shù)a11，a12，……，ass，b1，……，bs以及函數(shù)f的任意階導(dǎo)數(shù)f′（t0），f"（t0），……和它的多項式的等式。再由f的任意性可知，單項式，f′（t0），f"（t0），……的系數(shù)必須為0，得到若干個關(guān)于龍格-庫塔系數(shù)a11，a12，……，ass，b1，……，bs的等式組，從而形成一個Q[a11，a12，……，ass，b1，……，bs]的一個理想[15]。

4 軟件包應(yīng)用舉例

在仿射空間（a11，a12，……，ass，b1，……，bs）中描述仿射簇，該流體上的點的坐標(biāo)可以看作是一個滿足辛要求和s級數(shù)且具有p階近似的對角隱式辛龍格庫系數(shù)的一組取值。以下將該仿射簇記為I（s，p）。軟件包dis_rk.sage中函數(shù)rk_var_dis（s）返回s級數(shù)且具有p階近似的對角隱式辛龍格庫系數(shù)的所有需要的符號變量和環(huán)Kab=Q[a11，a12，……，ass，b1，……，bs]。按照算法1，p階近似的條件產(chǎn)生了一系列的方程，函數(shù)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

上述數(shù)學(xué)實驗結(jié)果dim（I（2，2））=1且genus（I（2，2））=0。

5 仿射簇I（s，p）的性質(zhì)

在這一部分，借助于軟件包dis_rk.sage，對對角隱式辛龍格-庫塔法的系數(shù)所構(gòu)成的仿射簇展開研究，給出其一些性質(zhì)。由這些系數(shù)得到的龍格-庫塔算法的格式不僅具有計算方便的優(yōu)勢，而且也對二次積分不變量守恒。不失一般性，我們假設(shè)bi≠0（i=1，…，s）。若bk=0，則辛條件導(dǎo)致了biaik=0（i=1，…，s），則該方法等價于低階的對應(yīng)算法。通過在計算代數(shù)系統(tǒng)Sagemath中的數(shù)學(xué)實驗，得到了以下關(guān)于仿射簇I（s，p）的維數(shù)（dimension）的結(jié)論：

維度（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維代表該集是空集。當(dāng)該理想的維數(shù)為1時，借助于軟件包可以計算其相應(yīng)的虧格（genus），借助數(shù)學(xué)實驗，得到以下結(jié)論：虧格（I（3，3））=1，虧格（I（2，2））=0。下面對5種情形的I（s，p）結(jié)構(gòu)展開具體的研究。

（1）情形I（2，2）

數(shù)學(xué)實驗I（2，2）對應(yīng)的仿射簇所代表的曲線的維數(shù)為1，虧格為0，由代數(shù)幾何知識可知仿射簇I（2，2）存在一個有理的參數(shù)化方程，即：

（2）情形I（2，3）

由上述結(jié)論可知，仿射簇I（2，3）是由有限個點組成的，其相應(yīng)的龍格—庫塔系數(shù)為：

注記2：I（2，3）不存在實系數(shù)元素，也就是說，不存在2級3階的實系數(shù)對角辛隱式龍格庫塔算法。

（3）情形I（3，3）

仿射簇I（3，3）對應(yīng)的龍格—庫塔系數(shù)為

其中b0，b1，b2滿足

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

若令b1=b0，計算代數(shù)系統(tǒng)Sagemath中的數(shù)學(xué)實驗給出了以下結(jié)論：

該結(jié)論與文獻(xiàn)[8]中第279頁的結(jié)論一致。其中a=1.351207，它是多項式6x3-12x2+6x-1=0的一個實根。

（4）情形I（3，4）

仿射簇I（3，4）是由有限個點組成的，計算代數(shù)系統(tǒng)Sagemath中的數(shù)學(xué)實驗顯示它包含以下兩類對角隱式辛龍格-庫塔算法。

類型1：

注記3：數(shù)學(xué)實驗顯示I（3，4）只含有一個實系數(shù)的點，也就是說，只有一個相應(yīng)的實系數(shù)對角隱式辛龍格庫塔算法。

（5）情形I（5，5）

很明顯，I（5，5）包含有限個點，因此相應(yīng)的對角隱式辛龍格庫塔算法是有限的。然而，數(shù)學(xué)實驗告訴我們I（5，5）不存在實系數(shù)的對角隱式辛龍格庫塔算法。總的來說，正如G. Sun所指出的那樣高于4階的實系數(shù)的對角辛龍格庫塔算法是不存在的[18]。

6 結(jié)論

本文首先在計算代數(shù)系統(tǒng)Sagemath中建立了軟件包dis_rk.sage。基于該軟件包，對由s級數(shù)且具有p階近似的對角隱式辛龍格庫塔法的系數(shù)所構(gòu)成的仿射簇I（s，p）在該數(shù)學(xué)軟件中進(jìn)行了一系列的數(shù)學(xué)實驗。數(shù)學(xué)實驗證明我們的軟件包是可行且正確有效的。對于I（s，p）（s，p≤4）給出了相應(yīng)的的結(jié)構(gòu)性質(zhì)。而且，數(shù)學(xué)實驗結(jié)果顯示高于4階的實系數(shù)的對角辛龍格庫塔算法不存在。

——————————

參考文獻(xiàn)：

〔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〕錢旭.幾類偏微分方程的保結(jié)構(gòu)算法研究[D].北京：國防科學(xué)技術(shù)大學(xué)，2014.

〔4〕郭本瑜.偏微分方程的差分方法[M].北京：科學(xué)出版社，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〕馮康，秦孟兆.哈密爾頓系統(tǒng)的辛幾何算法[M].杭州：浙江科學(xué)技術(shù)出版社，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〕秦孟兆，王雨順.偏微分方程中的保結(jié)構(gòu)算法[M].杭州：浙江科學(xué)技術(shù)出版社，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.