吳志剛, 徐小明

(中山大學(xué)·深圳 航空航天學(xué)院,深圳 518107)

1 引 言

多體系統(tǒng)的動(dòng)力學(xué)模型通常采用兩種坐標(biāo)描述方法,即最小坐標(biāo)與非最小坐標(biāo)[1]。與最小坐標(biāo)相比,非最小坐標(biāo)具有數(shù)學(xué)表達(dá)方便、簡(jiǎn)潔以及方法通用性強(qiáng)等優(yōu)勢(shì)。但是,非最小坐標(biāo)描述不可避免地引入了代數(shù)約束條件,導(dǎo)致其數(shù)學(xué)表達(dá)為微分代數(shù)方程組。與常微分方程相比,微分代數(shù)方程組的數(shù)值求解相對(duì)復(fù)雜,面臨著約束違約和剛性方程問(wèn)題等,其數(shù)值計(jì)算的精度與穩(wěn)定性受到嚴(yán)峻考驗(yàn)[2]。因此,相關(guān)問(wèn)題的數(shù)值計(jì)算方法研究一直是多體系統(tǒng)動(dòng)力學(xué)研究的重要課題。

近幾十年,在多體系統(tǒng)動(dòng)力學(xué)的應(yīng)用背景下,微分代數(shù)方程組的計(jì)算方法得到了快速發(fā)展[3]。對(duì)于約束哈密頓系統(tǒng),保辛算法可以較好地處理約束違約與剛性方程問(wèn)題,取得了良好的計(jì)算穩(wěn)定性與精度[4,5]。然而,經(jīng)典的數(shù)值計(jì)算格式存在累加的能量耗散和相位誤差問(wèn)題。雖然辛算法通過(guò)保證離散過(guò)程的辛對(duì)稱(chēng)性,進(jìn)而在數(shù)值計(jì)算過(guò)程中保持了系統(tǒng)的守恒律,解決了能量耗散問(wèn)題,但相位誤差隨時(shí)間快速積累導(dǎo)致算法精度仍受到極大的影響[6]。Reich[7]很早就通過(guò)誤差估計(jì)分析了單步辛算法相位誤差產(chǎn)生機(jī)制。邢譽(yù)峰等[8]較早地闡述了算法的累積相位誤差,并針對(duì)單自由度線性系統(tǒng)提出了單步保辛算法相位誤差的修正公式。隨后,多位學(xué)者提出了保辛算法的相位誤差補(bǔ)償方法[9-14]。但這些方法沒(méi)有推廣應(yīng)用于約束哈密頓系統(tǒng)。最近,Noh等[15-17]分析了算法相位誤差隨時(shí)間累積的問(wèn)題,稱(chēng)這種現(xiàn)象為算法的彌散性(Dispersion),并指出能量的耗散性和相位的彌散性都可以通過(guò)設(shè)定隱式積分格式的待定參數(shù)來(lái)改變,恰當(dāng)?shù)膮?shù)取值可以使積分算法在仿真分析中給出更為精確結(jié)果。然而,該方法的待定參數(shù)與離散方式相關(guān),改變其取值將破壞積分格式的原有結(jié)構(gòu)(如辛結(jié)構(gòu)),因此不能直接應(yīng)用。

本文針對(duì)非最小坐標(biāo)描述的平面剛體系統(tǒng),展開(kāi)具有低相位誤差積累特點(diǎn)的保辛積分構(gòu)造方法研究。首先,本文采用笛卡爾坐標(biāo)描述剛體的平移與旋轉(zhuǎn)運(yùn)動(dòng)。這是一種非最小坐標(biāo)描述,導(dǎo)致相關(guān)動(dòng)力學(xué)模型的質(zhì)量矩陣表達(dá)形式的不唯一性[18]。本文通過(guò)引入正交投影矩陣(Orthogonal Projection),將系統(tǒng)的質(zhì)量矩陣描述為關(guān)于待定參數(shù)的等價(jià)表達(dá)式,進(jìn)而推導(dǎo)出了改進(jìn)的拉格朗日方程[19]。研究發(fā)現(xiàn),離散過(guò)程中關(guān)于待定參數(shù)的項(xiàng)將產(chǎn)生額外的數(shù)值誤差。通過(guò)調(diào)整參數(shù)的大小和正負(fù),進(jìn)而調(diào)整相關(guān)誤差的大小和方向,可以達(dá)到抵消系統(tǒng)其他項(xiàng)誤差的效果,從而大幅降低數(shù)值方法相位誤差隨時(shí)間累積的速度。由于該方法中的待定參數(shù)是在推導(dǎo)動(dòng)力學(xué)方程過(guò)程中得到的,與離散過(guò)程無(wú)關(guān),所以待定參數(shù)取值的改變并不會(huì)破壞數(shù)值積分格式原有的離散性質(zhì)。利用這一思想,本文針對(duì)平面剛體系統(tǒng),提出了預(yù)參數(shù)調(diào)節(jié)保辛積分格式的構(gòu)造方法。數(shù)值結(jié)果表明,該方法在保持保辛積分格式無(wú)能量耗散、長(zhǎng)時(shí)間積分穩(wěn)定的情況下,同時(shí)大幅降低了離散系統(tǒng)的相位誤差,極大程度上改善了數(shù)值方法的軌跡精度。由于非最小坐標(biāo)描述中質(zhì)量矩陣的不唯一性普遍存在,所以該方法不僅可以用于平面剛體系統(tǒng),還可以推廣到其他約束多體系統(tǒng)。

2 改進(jìn)的拉格朗日方程

2.1 Lagrange方程描述

當(dāng)采用非最小坐標(biāo)描述平面剛體運(yùn)動(dòng)時(shí),選取的廣義坐標(biāo)維度大于獨(dú)立運(yùn)動(dòng)自由度數(shù),其動(dòng)力學(xué)方程可以表達(dá)為L(zhǎng)agrange方程形式

(1)

Φ(q)=0

(2)

式中q∈n為n維廣義坐標(biāo),p∈n為廣義動(dòng)量,T為動(dòng)能,Fex∈n為外力,Φ∈m為m個(gè)完整約束,A=?Φ/?q為約束的雅克比矩陣,λ∈m為拉格朗日乘子。假設(shè)動(dòng)能的具體表達(dá)式為

(3)

式中M(q)∈n×n為n維質(zhì)量矩陣。對(duì)于約束系統(tǒng),質(zhì)量矩陣一般為半正定矩陣。由解的存在唯一性條件,其矩陣的秩滿足[18]

rank(M)≥n-m

(4)

根據(jù)動(dòng)能表達(dá)式(3),系統(tǒng)的廣義動(dòng)量可以表達(dá)為

(5)

實(shí)際上,上述動(dòng)能表達(dá)式與唯一性條件表明,質(zhì)量矩陣中只有投影到n-m維約束零空間中的分量對(duì)動(dòng)能產(chǎn)生正的作用,而投影到m維約束空間中的分量只能產(chǎn)生零動(dòng)能。在物理上,任何可能的運(yùn)動(dòng)(即滿足約束方程的運(yùn)動(dòng))都不可能與零動(dòng)能相關(guān)聯(lián)。這暗示了約束系統(tǒng)質(zhì)量矩陣的不唯一性。

2.2 基于正交投影的等價(jià)表達(dá)式

為將質(zhì)量矩陣表達(dá)為關(guān)于待定參數(shù)的等價(jià)形式,定義約束的零空間矩陣

N(A)={z∈n:Az=0}

(6)

則可以定義S∈n×n為N(A)上的正交投影矩陣,滿足R(S)=N(A),S2=S,以及S2=S,其中R(S)為矩陣的值域(Range of the matrix)。根據(jù)正交矩陣定義,可以推導(dǎo)得到其正交補(bǔ)矩陣為I-S,對(duì)任意滿足其中I為單位矩陣。這樣就可以定義質(zhì)量矩陣的等價(jià)表達(dá)式

Mγ(q)=M(q)+γ(I-S)

(7)

式中γ為任意參數(shù)。將式(7)代入式(3)替換質(zhì)量矩陣M,可得具有待定參數(shù)的等價(jià)動(dòng)能表達(dá)式

(8)

與式(3)的質(zhì)量矩陣相比,式(8)增加了與參數(shù)γ相關(guān)的項(xiàng),其只能產(chǎn)生零動(dòng)能,不會(huì)對(duì)系統(tǒng)的真實(shí)運(yùn)動(dòng)狀態(tài)產(chǎn)生影響。將式(8)代入式(1),可以得到Lagrange方程的等價(jià)表達(dá)式

(9)

式中 廣義動(dòng)量p已經(jīng)重新定義為

(10)

式(9)稱(chēng)為改進(jìn)的拉格朗日方程MLEs(Modified Lagrange’s equations)。

3 參數(shù)預(yù)調(diào)節(jié)保辛算法思路

文獻(xiàn)[19]闡述了改進(jìn)的拉格朗日方程中與參數(shù)γ相關(guān)的廣義約束力項(xiàng)的離散誤差產(chǎn)生機(jī)制。誤差分析與數(shù)值結(jié)果表明,當(dāng)滿足一定的離散條件,參數(shù)γ相關(guān)項(xiàng)將產(chǎn)生額外的數(shù)值誤差。通過(guò)預(yù)先調(diào)整參數(shù)γ的取值,可以抵消數(shù)值積分產(chǎn)生的誤差。受此啟發(fā),本文提出利用改進(jìn)的拉格朗日方程構(gòu)造具有低相位誤差累積的保辛積分方法,具體思路如下。

(1) 在Lagrange框架或Hamilton框架下推導(dǎo)系統(tǒng)的動(dòng)力學(xué)方程,得到原始的質(zhì)量矩陣,如式(1～5)所示。

(2) 根據(jù)約束方程的具體形式,推導(dǎo)約束零空間上的正交投影矩陣S,得到其正交補(bǔ)映射I-S,進(jìn)而構(gòu)造具有待定參數(shù)的質(zhì)量矩陣與動(dòng)能表達(dá)式,如式(7,8)所示。

(3) 在上述公式基礎(chǔ)上,推導(dǎo)改進(jìn)的拉格朗日方程,如式(9,10)所示。

(4) 聯(lián)合改進(jìn)的拉格朗日方程與約束方程,構(gòu)造保辛算法。

(5) 通過(guò)先驗(yàn)誤差估計(jì)或后驗(yàn)誤差估計(jì),得到使得相位誤差顯著降低的待定參數(shù)取值。

在進(jìn)行參數(shù)預(yù)調(diào)節(jié)算法構(gòu)造時(shí),需要注意以下幾點(diǎn),離散格式可以采用針對(duì)指標(biāo)-3微分代數(shù)方程組的任意一種保辛格式,如辛龍格庫(kù)塔法或變分的算法等;待定參數(shù)是在仿真前選取的,仿真中不再改變。當(dāng)采用后驗(yàn)誤差估計(jì)方法確定待定參數(shù)γ取值時(shí),本文首先假設(shè)軌跡誤差δW是參數(shù)γ的線性函數(shù),計(jì)算兩個(gè)參考點(diǎn)(γ1,δW1)與(γ2,δW2)的值,然后令δW=0,通過(guò)線性插值就可以得到γ的估計(jì)值

γopt=(γ1δW2-γ2δW1)/(δW2-δW1)

(11)

即為使得軌跡誤差為零的最優(yōu)參數(shù)。在選取參考點(diǎn)時(shí),γ1一般可選為零,γ2可選為平面剛體系統(tǒng)轉(zhuǎn)動(dòng)慣量或質(zhì)量矩陣相關(guān)的常值,然后通過(guò)短時(shí)間仿真計(jì)算得到相應(yīng)的后驗(yàn)誤差值δW1與δW2,然后再將參考點(diǎn)的值代入式(11),求取最優(yōu)參數(shù)值γopt。文獻(xiàn)[20]更為詳細(xì)地介紹了確定待定參數(shù)γ的后驗(yàn)誤差估計(jì)方法,限于本文篇幅不再介紹。

4 平面剛體系統(tǒng)建模

4.1 單剛體系統(tǒng)

圖1 笛卡爾坐標(biāo)描述的平面剛體

xi=[xiyixx+1yi+1]T

(12)

與單位正交矢量

ui=Di1xi,vi=Di2xi

(13)

(14)

并且需要滿足約束條件

(15)

根據(jù)式(13),可以定義平面運(yùn)動(dòng)的剛體旋轉(zhuǎn)矩陣

此外，對(duì)微博的良好管理，離不開(kāi)對(duì)網(wǎng)絡(luò)技術(shù)的掌握。地方政府可以向官員組織開(kāi)設(shè)微博的相關(guān)課程，使地方政府官員熟知微博各項(xiàng)操作的具體方法，做到準(zhǔn)確地使用微博，避免因操作失誤引起誤會(huì)，影響政府官員微博的公信力。

Ri=[uivi]

(16)

這樣全局坐標(biāo)系下剛體的質(zhì)心坐標(biāo)可表示為

(17)

(18)

據(jù)此可以推導(dǎo)得到單剛體的動(dòng)能表達(dá)式

(19)

(20)

(21)

(22)

(23)

將式(13)代入式(23),可以得到單剛體的等價(jià)動(dòng)能表達(dá)式

(24)

(25)

4.2 重力場(chǎng)中的多體平面擺

為簡(jiǎn)化推導(dǎo)過(guò)程,本文只考慮如圖2所示平面多體擺系統(tǒng),其中Oi是平面鉸連接,坐標(biāo)為(xi,yi),點(diǎn)O0的坐標(biāo)為(0,0)。定義系統(tǒng)的廣義坐標(biāo)向量為

圖2 多體平面擺

q=[x1y1…xNyN]T

(26)

其滿足約束方程

(27)

剛體系統(tǒng)的總動(dòng)能表達(dá)式為

(28)

(29)

(30)

重力加速度g=9.81 m/s2。然后可以推導(dǎo)得到系統(tǒng)的外力Fex=-?V/?qT。

為獲得式(28)中參數(shù)γ的合理取值,可以定義參數(shù)[21]

(31)

為所有平面剛體轉(zhuǎn)動(dòng)慣量的均值,離散誤差估計(jì)表明當(dāng)γ=γR時(shí),數(shù)值積分的離散旋轉(zhuǎn)動(dòng)能誤差較小。進(jìn)一步,還可以考慮減小平移動(dòng)能的離散誤差,此時(shí)可以定義

(32)

5 數(shù)值算例

考慮由兩個(gè)剛體組成的兩自由度平面擺,其構(gòu)型參數(shù)為l1=2,l2=3,m1=m2=1以及

(33)

式中i=0,1。給定初始條件

(34)

本文采用Gauss-Lobatto SPARK算法來(lái)構(gòu)造保辛積分格式,具體格式參考文獻(xiàn)[22]。需要注意,由于式(28)的引入,保辛積分格式中含有待定系數(shù)γ。首先取待定系數(shù)γ=0與γ=γT=0.916667進(jìn)行數(shù)值結(jié)果比較,其中γT的定義見(jiàn)式(32)。為此,定義均值誤差

(35)

(36)

(37)

式中Nh為總的時(shí)間步數(shù),N=2為剛體數(shù)目,h為時(shí)間步長(zhǎng),δ(#)表示數(shù)值解和參考解之間的誤差,H=T+V為系統(tǒng)的能量,上標(biāo)k表示tk=hk時(shí)刻的值。本文參考解選取為比待比較數(shù)值解步長(zhǎng)小100倍的數(shù)值解。

表1給出了當(dāng)待定參數(shù)γ取0與取γT時(shí)的數(shù)值比較結(jié)果,其中h=0.04,Nh=100,(s,s)代表SPARK方法選取了s個(gè)Gauss-Lobatto積分點(diǎn)與s個(gè)Lobatto積分點(diǎn)進(jìn)行離散。若數(shù)值積分選取s個(gè)積分點(diǎn),則收斂階數(shù)為2s階。數(shù)值結(jié)果顯示,對(duì)于2階、4階、6階和8階SPARK方法,參數(shù)γ=γT時(shí)的軌跡精度與能量精度均優(yōu)于γ=0。證實(shí)了參數(shù)預(yù)調(diào)節(jié)保辛算法可以通過(guò)調(diào)整γ的取值改變數(shù)值積分的精度。圖3給出了能量誤差隨步長(zhǎng)收斂曲線,其中s=1?？梢钥闯?無(wú)論參數(shù)取何值,取1個(gè)積分點(diǎn)的SPARK方法都保證了2階收斂性。實(shí)際上,當(dāng)s=2,3和4時(shí),數(shù)值結(jié)果也顯示了算法與待定參數(shù)γ取值無(wú)關(guān)的2s階收斂性,限于篇幅本文不列出其相關(guān)數(shù)值結(jié)果。

表1 數(shù)值結(jié)果比較

圖3 能量誤差隨步長(zhǎng)收斂曲線

圖4呈現(xiàn)了數(shù)值誤差隨待定參數(shù)γ變化曲線,其中h=0.04,Nh=100,數(shù)值誤差由式(35～37)定義。如圖4所示,無(wú)論是s=1還是s=2,用逐點(diǎn)誤差絕對(duì)值評(píng)估的平均軌跡誤差與待定參數(shù)γ呈現(xiàn)了V字型的曲線。這些數(shù)值誤差在γ的正半軸達(dá)到最小值。這意味著,如果通過(guò)優(yōu)化γ來(lái)最小化積分的軌跡誤差,算法的數(shù)值精度將得到極大的提升。然而,上述曲線中軌跡誤差和能量誤差曲線的最小值對(duì)應(yīng)的橫坐標(biāo)(γ的取值)不同,因此應(yīng)當(dāng)考慮能量和軌跡的綜合精度來(lái)確定γ的最佳值。本文以軌跡精度最佳為目標(biāo)來(lái)求參數(shù)γ的最佳取值γopt。

圖4 數(shù)值誤差隨待定參數(shù)γ變化曲線

圖5顯示了s=1時(shí)算法數(shù)值誤差隨參數(shù)γ變化曲線,其中h=0.01,Nh=1000。盡管本文選擇了不同的時(shí)間步長(zhǎng)和更長(zhǎng)的總仿真時(shí)間,但是圖5與圖4(a)的誤差曲線是精確一致的。這意味著雖然時(shí)間步長(zhǎng)變化了,但最優(yōu)軌跡精度的參數(shù)γopt幾乎不受時(shí)間步長(zhǎng)與仿真時(shí)間長(zhǎng)短的影響。為進(jìn)一步分析軌跡誤差隨參數(shù)γ的變化規(guī)律,進(jìn)一步定義如下均值軌跡誤差

(38)

(39)

通過(guò)上述定義,軌跡的均值誤差變成可正可負(fù)。圖6給出了采用式(38,39)定義的誤差隨參數(shù)γ變化趨勢(shì),其中h=0.01,Nh=10。數(shù)值結(jié)果顯示軌跡誤差隨參數(shù)γ呈現(xiàn)線性變化趨勢(shì),并且使誤差最小的參數(shù)γopt與圖5的值一致。因此,可以采用較少的仿真時(shí)長(zhǎng)預(yù)估最優(yōu)參數(shù)γopt。

為預(yù)估γopt的最優(yōu)參數(shù)值,本文定義軌跡誤差均值的和為δW=E[δx]+E[δy]。如圖6所示,誤差δW隨參數(shù)γ呈現(xiàn)線性變化趨勢(shì)。進(jìn)一步令γ1=0,γ2=γT,δW1=δW(γ1),δW2=δW(γ2),然后將這些取值代入式(11),計(jì)算得到γopt。

圖6 軌跡均值誤差-γ變化曲線趨勢(shì)比較

表2給出了γ1=0和γ2=γT時(shí)的誤差δW1與δW2的取值,以及通過(guò)式(11)計(jì)算出的γopt取值,其中h=0.01,Nh=10。表3給出了參數(shù)γ=γopt以及h=0.04時(shí)軌跡誤差?？梢钥闯?s=1和s=2時(shí),x方向和y方向的軌跡誤差分別減小了兩個(gè)和一個(gè)數(shù)量級(jí)。圖7給出了s=1和h=0.04時(shí)多體平面擺x2分量與y2分量的數(shù)值結(jié)果,其中黑色線為解析解,紅色○代表γ=0的數(shù)值解,藍(lán)色×代表γ=γopt的數(shù)值解。數(shù)值結(jié)果表明,通過(guò)預(yù)先調(diào)節(jié)待定參數(shù)γ的取值,數(shù)值積分的相位誤差得到了大幅降低。實(shí)際上對(duì)于s≥2的高階情況,通過(guò)預(yù)參數(shù)調(diào)節(jié)同樣可以大幅降低,本文限于篇幅略去相關(guān)討論。需要指出,通過(guò)式(11)估計(jì)最優(yōu)參數(shù)γopt需要一定的額外計(jì)算量(不隨仿真時(shí)間而增長(zhǎng))。所以只有在長(zhǎng)時(shí)間仿真和額外的計(jì)算量可以忽略不計(jì)時(shí),該方法才具有較大優(yōu)勢(shì)。

表2 最優(yōu)參數(shù)γopt計(jì)算

表3 優(yōu)化精度后的數(shù)值結(jié)果

圖7 兩體平面擺數(shù)值結(jié)果比較

6 結(jié) 論

針對(duì)笛卡爾坐標(biāo)描述的平面多剛體系統(tǒng),本文提出了一種參數(shù)預(yù)調(diào)節(jié)保辛積分方法。在該方法構(gòu)造過(guò)程中,首先引入約束零空間的正交投影矩陣,進(jìn)而推導(dǎo)了平面多剛體系統(tǒng)的改進(jìn)的拉格朗日方程,將動(dòng)力學(xué)方程描述成包含待定參數(shù)γ的微分代數(shù)方程組,然后參數(shù)預(yù)調(diào)節(jié)的保辛積分方法通過(guò)離散化改進(jìn)的拉格朗日方程得到。為了減小數(shù)值方法的相位誤差累積,本文利用數(shù)值積分的軌跡誤差與參數(shù)γ的近似線性關(guān)系,通過(guò)對(duì)誤差線性插值,得到最優(yōu)的參數(shù)值γopt。數(shù)值結(jié)果表明通過(guò)預(yù)先估計(jì)最優(yōu)的待定參數(shù)γ,包含兩個(gè)剛體平面擺的數(shù)值積分精度提升了1～2個(gè)數(shù)量級(jí),較大幅度地減小了保辛算法隨時(shí)間累積的相位誤差。