基于Lagrange插值多項(xiàng)式擬合的力學(xué)系統(tǒng)的變分積分子*

孔新雷 楊雪

（北方工業(yè)大學(xué)理學(xué)院，北京，100144）

引言

Lagrange系統(tǒng)是由Euler-Lagrange方程

所描述的一類力學(xué)系統(tǒng)，其中x=(x1，x1，…，xn)代表位 形 流形M的局部坐標(biāo) ，L(x，x?)是系統(tǒng)的Lagrange函數(shù)，一般為系統(tǒng)的動(dòng)能和勢(shì)能之差.如果進(jìn)一步考慮非保守力，那么Euler-Lagrange方程的形式會(huì)發(fā)生些許變化［1］，但這種情況在本文中暫不研究.當(dāng)系統(tǒng)（1）是一組非線性微分方程時(shí)，對(duì)其解析求解是非常困難甚至是不可能的，因此需要借助于數(shù)值手段，通常是采用經(jīng)典差分格式，用差分方程代替微分方程實(shí)現(xiàn)數(shù)值求解.

然而，人為地采用經(jīng)典差分格式離散系統(tǒng)的運(yùn)動(dòng)方程，不可避免地會(huì)破壞連續(xù)系統(tǒng)具有的幾何結(jié)構(gòu)，引入不必要的數(shù)值耗散.因此，在“數(shù)值算法應(yīng)盡可能地保持原問(wèn)題的本質(zhì)特征”［2］的指導(dǎo)原則下，針對(duì)系統(tǒng)（1），人們提出并設(shè)計(jì)了保結(jié)構(gòu)算法［2，3］.較早的保結(jié)構(gòu)算法主要是保辛差分格式 .借助于Legendre變換，Euler-Lagrange方程可以轉(zhuǎn)化為Hamilton正則方程，而基于Hamilton系統(tǒng)的辛結(jié)構(gòu)，利用生成函數(shù)法可以構(gòu)造系統(tǒng)的保辛差分格式［2］.變分積分子是另外一種構(gòu)造思路的保結(jié)構(gòu)算法.考慮到Euler-Lagrange方程可以由Hamilton原理誘導(dǎo)，仿照這一連續(xù)情形，在對(duì)系統(tǒng)（1）構(gòu)造數(shù)值算法時(shí)可以先離散Hamilton原理，由離散后的Hamilton原理可以導(dǎo)出離散Euler-Lagrange方程.作為離散變分原理的產(chǎn)物，離散Euler-Lagrange方程在作為數(shù)值差分格式時(shí)，真實(shí)地繼承了連續(xù)系統(tǒng)具有的幾何特性，不僅是自然保辛的，而且滿足Noether定理進(jìn)而保持系統(tǒng)的動(dòng)量［4，5］.自然簡(jiǎn)便的構(gòu)造方式和優(yōu)越顯著的計(jì)算性能，使得變分積分子迅速發(fā)展，不僅局限于系統(tǒng)（1），對(duì)于非保守系統(tǒng)［5，6］、非光滑系統(tǒng)［7］、完整和非完整約束系統(tǒng)［8-10］、Birkhoff系 統(tǒng)［11-13］、連續(xù)介質(zhì)力學(xué)［14，15］、隨機(jī)系統(tǒng)［16，17］等都可以構(gòu)造系統(tǒng)的變分積分子，進(jìn)行高效的數(shù)值模擬仿真.

本文結(jié)合局部路徑擬合［18，19］方法提出一種新的變分積分子的設(shè)計(jì)構(gòu)造方式.考慮到變分積分子的經(jīng)典構(gòu)造過(guò)程最終歸結(jié)為計(jì)算離散Lagrange函數(shù)的偏導(dǎo)數(shù)，本文先對(duì)相關(guān)涉及到的偏導(dǎo)數(shù)進(jìn)行解析計(jì)算，發(fā)現(xiàn)一個(gè)新的依賴于Euler-Lagrange方程的積分隨之產(chǎn)生.在對(duì)該積分實(shí)現(xiàn)科學(xué)有效的近似之后，離散Lagrange函數(shù)的偏導(dǎo)數(shù)就轉(zhuǎn)化為連續(xù)Lagrange函數(shù)關(guān)于速度變量的偏導(dǎo)數(shù)，這樣就省略了利用經(jīng)典求積公式計(jì)算離散Lagrange函數(shù)的過(guò)程.因此，這種新的構(gòu)造方式不僅讓最終得到的差分格式仍然繼承了變分積分子特有的優(yōu)越計(jì)算性能，而且顯著地簡(jiǎn)化了其構(gòu)造過(guò)程，特別是對(duì)于高階精度的變分積分子.

1 變分積分子

如引言中所述，Lagrange力學(xué)系統(tǒng)是由Euler-Lagrange方程

變分積分子是一類特殊的數(shù)值差分格式.與經(jīng)典差分格式不同，構(gòu)造力學(xué)系統(tǒng)的變分積分子，出發(fā)點(diǎn)不是直接離散運(yùn)動(dòng)方程，而是離散能夠誘導(dǎo)運(yùn)動(dòng)方程的變分原理.

具體到Lagrange力學(xué)系統(tǒng)，首先考慮離散Lagrange函數(shù)

其中，xk是x(tk)的近似，τ=tk+1-tk是選定的時(shí)間步長(zhǎng).這樣

進(jìn) 一 步 注 意 到δx0=δx(0)=0 以 及δxN=δx(T)=0，可得離散Euler-Lagrange方程

寫(xiě)成分量形式為

當(dāng)?2Ld滿足非退化條件時(shí)，離散Euler-Lagrange方程（3）就確定了差分格式

這種由離散變分原理所得到的數(shù)值算法Φ就被稱作變分積分子.獨(dú)特的構(gòu)造方式避免了人為的離散運(yùn)動(dòng)方程所帶來(lái)的數(shù)值耗散，使得變分積分子能夠真實(shí)地繼承連續(xù)系統(tǒng)的解所具有的幾何性質(zhì)，例如保辛性，即算法Φ滿足Φ?Ω=Ω，其中微分形式

另外，與連續(xù)情形類似，當(dāng)離散Lagrange函數(shù)Ld進(jìn)一步滿足李群作用不變性時(shí)，算法Φ還保持動(dòng)量（映射）不變.根據(jù)后誤差分析理論，變分積分子所表現(xiàn)出的優(yōu)越計(jì)算性能，如高精度、長(zhǎng)時(shí)間穩(wěn)定、保持系統(tǒng)守恒量等，都與其保結(jié)構(gòu)特性有著直接關(guān)系.當(dāng)力學(xué)系統(tǒng)受到若干完整約束，即位形坐標(biāo)x=(x1，x1，…，xn)滿足代數(shù)方程組

時(shí)，借鑒上述先離散后變分的思想，并結(jié)合Lagrange乘子法，同樣可以構(gòu)造完整約束系統(tǒng)

的變分積分子

其中，

值得注意的是，如果引入擴(kuò)展的離散Lagrange函數(shù)

或

那么，方程組（5）與以Lˉd(xˉk，xˉk+1)為離散 Lagrange函數(shù)的離散Euler-Lagrange方程

等價(jià)［5］.

2 基于Lagrange插值多項(xiàng)式擬合的力學(xué)系統(tǒng)的變分積分子

根據(jù)離散Euler-Lagrange方程（3）的形式可知，構(gòu)造力學(xué)系統(tǒng)的變分積分子最后歸結(jié)為計(jì)算離散Lagrange函數(shù)的偏導(dǎo)數(shù)——?1Ld和?2Ld.為了計(jì)算這些偏導(dǎo)數(shù)，首先引入自由變量μ，在實(shí)際計(jì)算偏導(dǎo)數(shù)時(shí)μ可以代表或.由于離散Lagrange函數(shù)

那么，

進(jìn)一步地，利用分部積分公式可得

在上式中，記

同理，當(dāng)考慮時(shí)間區(qū)間[tk-1，tk]時(shí)，

假設(shè)在時(shí)間區(qū)間[tk，tk+1]內(nèi)，

其中，

而ηα=(，，…)T是擬合軌跡x(t)在離散時(shí)間結(jié)點(diǎn)處的值，并且滿足Euler-Lagrange方程，那么根據(jù)x(t)的表達(dá)式可以計(jì)算x?(t)、x?(t)以及

這樣，由（7）式可得

同理，當(dāng)考慮時(shí)間區(qū)間[tk-1，tk]時(shí)，

利用（9）式和（10）式，聯(lián)立方程組

其中，i=1，2，…，n，α=1，2，…，m-1，則首先由方程（12）可求得ηα=ηα(xk，xk+1)，平行地，當(dāng)時(shí)間區(qū)間為[tk-1，tk]時(shí)，ηα=ηα(xk-1，xk)；然后將求得的ηα代入x(tk)和x?(tk)，并進(jìn)一步帶入（9）式和（10）式，那么方程（11）就轉(zhuǎn)化關(guān)于xk-1、xk和xk+1的代數(shù)方程，即變分積分子Φ：(xk-1，xk)→ (xk，xk+1).

注1：?jiǎn)蝹€(gè)步長(zhǎng)區(qū)間[tk，tk+1]內(nèi)，在對(duì)局部軌跡x(t)進(jìn)行插值擬合時(shí)，除了Lagrange插值多項(xiàng)式（8）之外，還可以采用其他的插值方式，例如二項(xiàng)插值

式可得

代入偏導(dǎo)數(shù)?1Ld和?2Ld，同樣可以得到

對(duì)應(yīng)地

記

平行地

在時(shí)間區(qū)間[tk，tk+1]內(nèi)，如果局部軌跡x(t)和λ(t)都按照Lagrange插值多項(xiàng)式進(jìn)行擬合，則最終可得

而方程（14）則隨之變?yōu)?/p>

3 數(shù)值算例

3.1 算例1—簡(jiǎn)諧振子

考慮簡(jiǎn)諧振子方程

在時(shí)間區(qū)間[tk，tk+1]內(nèi)，令

則求導(dǎo)可得

在t=tk+τ/2處，由Euler-Lagrange方程，即系統(tǒng)運(yùn)動(dòng)方程（15）成立可得

將η1帶入x?(t)相應(yīng)可得

那么，由離散Euler-Lagrange方程（11）可得

迭代算法（16）就是利用新的構(gòu)造方法針對(duì)簡(jiǎn)諧振子方程（15）設(shè)計(jì)的變分積分子.假定系統(tǒng)（15）的初值為x(0)=x?(0)=1，以τ=0.1為時(shí)間步長(zhǎng)，利用差分格式（16）對(duì)系統(tǒng)進(jìn)行數(shù)值求解，則對(duì)應(yīng)的數(shù)值解如圖1所示.顯然，該差分格式非常精確地?cái)?shù)值模擬了簡(jiǎn)諧振動(dòng)，在圖示的尺度下，方程的數(shù)值解和精確解所對(duì)應(yīng)的曲線完全重合，無(wú)任何差別.

圖1 方程（15）的精確解與通過(guò)差分格式（16）算得的數(shù)值解對(duì)照Fig.1 Comparison between the exact solution and the numerical solution computed by algorithm（16）for equation（15）

除了有效性得以驗(yàn)證之外，由圖2可以非常直觀地看到，差分格式（16）也非常真實(shí)地模擬了系統(tǒng)的能量，所計(jì)算出的數(shù)值能量曲線呈現(xiàn)出長(zhǎng)時(shí)間的守恒行為，且能量誤差始終保持在有界的范圍內(nèi).保持保守系統(tǒng)的守恒量，這是變分積分子特有的計(jì)算優(yōu)越性，因此，圖2中的數(shù)值結(jié)果也驗(yàn)證了通過(guò)局部路徑擬合方法構(gòu)造的變分積分子繼承了經(jīng)典變分積分子優(yōu)越的計(jì)算性能.實(shí)際上不僅如此，在同階精度的情況下，借助新方法得到的變分積分子具有更小的誤差.

圖2 由差分格式（16）算得的系統(tǒng)（15）的數(shù)值能量曲線Fig.2 The numerical energy curve of system（15）computed by algorithm（16）

為了驗(yàn)證這一事實(shí)，采用經(jīng)典求積公式—中點(diǎn)格式—逼近（2）式中的積分可得與系統(tǒng)（15）對(duì)應(yīng)的離散Lagrange函數(shù)

根據(jù)離散Euler-Lagrange方程（3）可得經(jīng)典變分積分子

差分格式（16）和差分格式（17）盡管都是具有二階精度的數(shù)值算法，然而由表1中數(shù)據(jù)可以直觀地對(duì)比出差分格式（16）具有更小的局部誤差和整體誤差.

表1 不同時(shí)刻差分格式(16)和差分格式(17)的誤差對(duì)比Table 1 Local errors of numerical solutions and energies at different times computed with algorithms(16)and(17)respectively

當(dāng)然，如果單個(gè)步長(zhǎng)區(qū)間[tk，tk+1]內(nèi)的離散結(jié)點(diǎn)取得更加密集，那么可以得到更加精確的差分格式.例如，在時(shí)間區(qū)間[tk，tk+1]內(nèi)，如果設(shè)

對(duì)應(yīng)地

代入離散Euler-Lagrange方程（11）可得差分格式

又例如，如果在[tk，tk+1]內(nèi)再多取一個(gè)結(jié)點(diǎn)，即令

則最終可得變分積分子

如表2中數(shù)據(jù)所示，差分格式（18）和差分格式（19）較之差分格式（16）具有更小的誤差.內(nèi)部結(jié)點(diǎn)選取得越密集，得到的數(shù)值算法越精確，這一規(guī)律與經(jīng)典構(gòu)造過(guò)程中愈加復(fù)雜的積分逼近格式誘導(dǎo)愈加精確的變分積分子相一致.

表2 不同時(shí)刻差分格式(16-19)的數(shù)值解誤差對(duì)比Table 2 Local errors of numerical solutions at different times computed with algorithms(16-19)respectively

3.2 算例2—單擺

單擺運(yùn)動(dòng)方程

也可以等價(jià)地表述為Euler-Lagrange方程，對(duì)應(yīng)的Lagrange函數(shù)為

在單個(gè)步長(zhǎng)區(qū)間內(nèi)如果只選取一個(gè)結(jié)點(diǎn)，那么對(duì)于系統(tǒng)（20）可得如下方程組

或

其中方程組（21）是基于Lagrange插值（8）誘導(dǎo)得來(lái)，而方程組（22）是基于二項(xiàng)插值（13）得來(lái).隱式地求解方程組（21）或（22）均可得變分積分子(θk，θk+1)= Φ(θk-1，θk).利用兩種差分格式算得系統(tǒng)（20）在相空間中的數(shù)值軌跡和數(shù)值能量如圖3所示.顯然，兩種差分格式不僅有效地模擬了系統(tǒng)的運(yùn)動(dòng)，而且非常真實(shí)地反映了系統(tǒng)能量守恒.

圖3 由差分格式（21）和（22）分別算得的系統(tǒng)（20）在相空間中的數(shù)值軌跡（a）和數(shù)值能量曲線（b）Fig.3 Numerical orbits in the phase space（a）and numerical energy curves（b）of system（20）computed by algorithms（21）and（22）respectively

3.3 算例3—?jiǎng)蛩賵A周運(yùn)動(dòng)

單位質(zhì)量的質(zhì)點(diǎn)在二維平面內(nèi)做勻速圓周運(yùn)動(dòng)，其運(yùn)動(dòng)方程可以表述為

其中位形坐標(biāo)(x1，x2)滿足約束條件x21+x22=1.基于Lagrange插值多項(xiàng)式（8），可得該完整約束系統(tǒng)的變分積分子

式（23）構(gòu)成了一個(gè)具有二階精度的隱式差分格式.

4 結(jié)論

本文結(jié)合局部路徑擬合方法提出了一種新的力學(xué)系統(tǒng)變分積分子的設(shè)計(jì)方式.數(shù)值結(jié)果表明，這種基于局部路徑擬合的構(gòu)造方法不僅讓最終得到的數(shù)值算法仍然繼承了變分積分子特有的優(yōu)越計(jì)算性能，而且顯著地提高了算法的精度.另外，這種構(gòu)造方法不僅適用于自由無(wú)約束系統(tǒng)，而且適用于完整約束系統(tǒng).