吳志剛, 徐小明

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

1 引 言

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

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

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

2 改進的拉格朗日方程

2.1 Lagrange方程描述

當采用非最小坐標描述平面剛體運動時,選取的廣義坐標維度大于獨立運動自由度數(shù),其動力學方程可以表達為Lagrange方程形式

(1)

Φ(q)=0

(2)

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

(3)

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

rank(M)≥n-m

(4)

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

(5)

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

2.2 基于正交投影的等價表達式

為將質量矩陣表達為關于待定參數(shù)的等價形式,定義約束的零空間矩陣

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ù)正交矩陣定義,可以推導得到其正交補矩陣為I-S,對任意滿足其中I為單位矩陣。這樣就可以定義質量矩陣的等價表達式

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

(7)

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

(8)

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

(9)

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

(10)

式(9)稱為改進的拉格朗日方程MLEs(Modified Lagrange’s equations)。

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

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

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

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

(3) 在上述公式基礎上,推導改進的拉格朗日方程,如式(9,10)所示。

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

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

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

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

(11)

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

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

4.1 單剛體系統(tǒng)

圖1 笛卡爾坐標描述的平面剛體

xi=[xiyixx+1yi+1]T

(12)

與單位正交矢量

ui=Di1xi,vi=Di2xi

(13)

(14)

并且需要滿足約束條件

(15)

根據(jù)式(13),可以定義平面運動的剛體旋轉矩陣

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

Ri=[uivi]

(16)

這樣全局坐標系下剛體的質心坐標可表示為

(17)

(18)

據(jù)此可以推導得到單剛體的動能表達式

(19)

(20)

(21)

(22)

(23)

將式(13)代入式(23),可以得到單剛體的等價動能表達式

(24)

(25)

4.2 重力場中的多體平面擺

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

圖2 多體平面擺

q=[x1y1…xNyN]T

(26)

其滿足約束方程

(27)

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

(28)

(29)

(30)

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

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

(31)

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

(32)

5 數(shù)值算例

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

(33)

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

(34)

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

(35)

(36)

(37)

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

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

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

圖3 能量誤差隨步長收斂曲線

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

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

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

(38)

(39)

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

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

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

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

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

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

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

6 結 論

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