薛冰 解加芳 張可心
(北方工業(yè)大學(xué) 理學(xué)院, 北京 100144)
非完整系統(tǒng)是一類(lèi)受到不可積微分約束的動(dòng)力學(xué)系統(tǒng)[1],廣泛應(yīng)用于場(chǎng)論、機(jī)電動(dòng)力系統(tǒng)、控制理論、工程科學(xué)等領(lǐng)域[2].20世紀(jì)80年代我國(guó)學(xué)者梅鳳翔研究了帶參數(shù)約束的一類(lèi)可控系統(tǒng)[3]、變質(zhì)量非完整約束系統(tǒng).Maggi在1896年推廣了拉格朗日第二類(lèi)方程[4],對(duì)線(xiàn)性非完整約束系統(tǒng)得到一類(lèi)動(dòng)力學(xué)方程,后人稱(chēng)為Maggi方程[5],這些方程后來(lái)被推廣到非線(xiàn)性非完整系統(tǒng)[6].Maggi方程是力學(xué)系統(tǒng)[7]各大運(yùn)動(dòng)方程的中間產(chǎn)物,對(duì)研究非完整系統(tǒng)的運(yùn)動(dòng)具有重要意義.
Birkhoff動(dòng)力學(xué)理論是Hamilton動(dòng)力學(xué)的自然推廣,它是包括齊次Hamilton系統(tǒng)和非齊次Hamilton系統(tǒng)的更一般動(dòng)力學(xué)理論,是最一般辛結(jié)構(gòu)的局部實(shí)現(xiàn),只有Birkhoff系統(tǒng)與一般辛幾何結(jié)構(gòu)之間才有一一對(duì)應(yīng)關(guān)系.因此 Birkhoff 系統(tǒng)動(dòng)力學(xué)[8]的研究對(duì)于完善和深化分析力學(xué)的理論體系具有重要意義, 尤其是對(duì)于非齊次 Hamilton 動(dòng)力學(xué)系統(tǒng)的幾何結(jié)構(gòu)分析具有重要應(yīng)用價(jià)值[9].本文針對(duì)非完整系統(tǒng)高階Maggi方程,在其滿(mǎn)足一定條件下,將其進(jìn)行Birkhoff化。并通過(guò)一個(gè)算例驗(yàn)證上述理論分析的正確性,再分別采用Runge-Kutta方法和Birkhoff辛算法對(duì)其進(jìn)行數(shù)值計(jì)算[10],并將數(shù)值結(jié)果進(jìn)行比較,給出Birkhoff辛算法在長(zhǎng)時(shí)計(jì)算時(shí)的優(yōu)越性.
設(shè)力學(xué)系統(tǒng)的位形由n個(gè)廣義坐標(biāo)qs(s=1,…,n)確定,系統(tǒng)受有g(shù)個(gè)理想m階非完整約束[11]
(1)
其中
根據(jù)d’Alembert-Lagrange原理可以導(dǎo)出Maggi形式為:
(σ=1,…,ε)
(2)
式中Qσ為廣義力,T為系統(tǒng)動(dòng)能.
令
(ν=1,…,ε;k=1,...,n)
(3)
則方程(2)有形式
(4)
(5)
約束對(duì)初始條件的限制為:
(6)
(7)
約束對(duì)初始條件的限制為:
(8)
將式(7)化為標(biāo)準(zhǔn)一階形式,令
(s=1,…,n)
(9)
則式(7)可以寫(xiě)成形式
(10)
其中
aν=xν,σs=xn+s,…,
σ(m-2)n+s=x(m-1)n+s,
σ(m-1)n+s=hs(s=1,…,n)
(11)
為將式(10)表為Birkhoff形式,其階必為偶數(shù)[9],即mn=2N,如果mn為奇數(shù)2N-1,可增加一個(gè)方程
(12)
使其成為偶階.從而,要使高階非完整系統(tǒng)的Maggi方程可表示成Birkhoff形式
(μ=1,…,2N)
(13)
即要求滿(mǎn)足
(μ=1,…,2N)
(14)
設(shè)式(10)的2n個(gè)第一積分Iμ(t,a)彼此無(wú)關(guān),即
(μ=1,…,2n)
(15)
(16)
根據(jù)Hojman方法,則式(14)的Birkhoff函數(shù)組Rμ由下式確定
(17)
Birkhoff量B為:
(18)
其中Gα需滿(mǎn)足條件
(19)
(20)
其中
(21)
假定方程組(20)的一種離散可以記為:
[?tR(zi,ti)]i
(22)
式(22)稱(chēng)為式(20)的一個(gè)離散格式,如果它決定的格式Φ保持離散的K(z,t)辛格式,即
(23)
(24)
以及
(25)
此映射為:
(26)
式(26)的Jacobi矩陣為:
(27)
映射α要滿(mǎn)足
(28)
其中
當(dāng)K與z無(wú)關(guān)時(shí),根據(jù)式(27)和式(28)可以得到
(29)
且需滿(mǎn)足下面的截面條件
|CαM+Dα|≠0
設(shè)Bα=Cα=0,則式(29)可成
(30)
由式(30)可以得到
(31)
生成函數(shù)φ(ω,t,t0)=α(t,t0),可以構(gòu)造廣義Birkhoff辛差分格式.當(dāng)步長(zhǎng)τ>0足夠小的時(shí)候,取
(m=1,2,…)
(32)
α1(zk+1,zk,tk+1,tk)
(33)
假設(shè)某一力學(xué)系統(tǒng)的位形由2個(gè)廣義坐標(biāo)q1,q2確定,其系統(tǒng)動(dòng)能為
(34)
且該系統(tǒng)受到1個(gè)2階非完整約束:
(35)
(36)
顯然式(36)有解:
q1=1-arctan(t)
q2=-3t·arctan(t)+2ln(t2+1)+t
由此可以構(gòu)造Birkhoff函數(shù)Rμ,B如下[11]:
R2=0
R4=0
從而將該系統(tǒng)的Maggi方程可以Birkhoff化:
結(jié)合式(24)和式(25),我們可以確定生成函數(shù)φ(ω,t,t0),之后采用上述的Birkhoff廣義辛算法式(32), 得到此 Birkhoff 系統(tǒng)的二階K(z,t)離散格式[11-13].
對(duì)該題采用二階K(z,t)算法和二階Runge-Kutta算法進(jìn)行計(jì)算.在計(jì)算過(guò)程中,先取如下初值:q1=1,C1=C2=1,取步長(zhǎng)τ=0.01.并通過(guò)比較兩種數(shù)值方法計(jì)算所得數(shù)值解和解析解q1=1-arctan(t)之間的相對(duì)誤差來(lái)說(shuō)明兩種數(shù)值方法的差別.
對(duì)比圖1和圖2可以看出,Runge-Kutta方法在長(zhǎng)期跟蹤后與解析解有著大幅度的相對(duì)誤差,而B(niǎo)irkhoff辛算法則相對(duì)誤差非常的小.因此,Birkhoff辛算法結(jié)果更加精確.
圖1 二階Runge-Kutta算法相對(duì)誤差Fig.1 Relative error of second order Runge-Kutta algorithm
圖2 Birkhoff辛算法相對(duì)誤差Fig.2 Relative error of Birkhoff symplectic algorithm
本文對(duì)一定條件下的非完整系統(tǒng)的高階Maggi方程(2)先進(jìn)行了Birkhoff化,得到廣義Birkhoff方程(13),并針對(duì)該方程,應(yīng)用Birkhoff廣義辛差分格式與傳統(tǒng)Runge-Kutta算法分別進(jìn)行計(jì)算,比較兩種算法,最后得出Birkhoff廣義辛差分算法在求解非完整系統(tǒng)高階Maggi方程中更加優(yōu)越.
動(dòng)力學(xué)與控制學(xué)報(bào)2024年1期