四元數(shù)與歐拉角剛體動力學(xué)數(shù)值積分算法及其比較

徐小明　鐘萬勰

摘要：為推廣四元數(shù)保辛積分在工程中的應(yīng)用，對歐拉角表示的狀態(tài)方程數(shù)值積分與四元數(shù)的保辛積分進(jìn)行比較.重陀螺的數(shù)值仿真結(jié)果表明四元數(shù)保辛積分的數(shù)值結(jié)果明顯優(yōu)于歐拉角狀態(tài)方程積分.與歐拉角狀態(tài)方程積分相比，四元數(shù)保辛積分在剛體動力學(xué)的數(shù)值仿真中更具優(yōu)勢.

關(guān)鍵詞：四元數(shù)； 歐拉角； 剛體動力學(xué)； 保辛； 重陀螺

中圖分類號： TP391.9； O313.3

文獻(xiàn)標(biāo)志碼： A

0引言

四元數(shù)、歐拉角和方向余弦[1]是描述剛體旋轉(zhuǎn)最主要的3種坐標(biāo)形式.方向余弦法需要9個參量，應(yīng)用較少；而另外2種則應(yīng)用廣泛，如航天器姿態(tài)控制和捷聯(lián)式慣性導(dǎo)航系統(tǒng)[1]等，對于兩者的研究也卷帙浩繁，但對兩者優(yōu)劣的評價卻褒貶不一.

描述剛體在三維空間中的運動姿態(tài)可采用2類12種歐拉角系統(tǒng)，分別對應(yīng)于不同的旋轉(zhuǎn)軸先后次序.目前公認(rèn)的用歐拉角描述旋轉(zhuǎn)的固有缺陷是奇異性問題[2]，即：無論采用哪種歐拉角系統(tǒng)，都不可避免地會含有奇異點，使得在該點附近區(qū)域進(jìn)行的數(shù)值積分精度不高.對于小角度旋轉(zhuǎn)，只要采用適當(dāng)?shù)臍W拉角系統(tǒng)便可避開奇異點；而在大角度旋轉(zhuǎn)時，若想避開奇異點，必須提供2套歐拉角系統(tǒng)交替進(jìn)行計算.據(jù)此，黃雪樵[3]提出一種“雙歐法”的計算方法；近幾年仍有學(xué)者[4]在繼續(xù)研究該方法.雙歐法雖然解決奇異性問題，但是計算過程較復(fù)雜.

四元數(shù)用于描述剛體旋轉(zhuǎn)，沒有奇異性問題，可很好地描述剛體的全角度旋轉(zhuǎn).然而，四元數(shù)需要滿足長度等于1的單位約束，這成為制約其應(yīng)用的限制.在實際應(yīng)用中，經(jīng)常采用的正交化修正等方法只能緩解長度的偏移，無法從根本上解決問題；黃雪樵[3]在其雙歐法中也有所討論.目前，對于單位約束最佳的解決方案是將四元數(shù)表示的剛體動力學(xué)方程導(dǎo)入微分代數(shù)方程范疇，近年來逐漸有學(xué)者[57]展開相關(guān)問題的研究.徐小明等[8]提出一種基于四元數(shù)理論描述剛體旋轉(zhuǎn)的保辛數(shù)值積分方法.該方法先將問題導(dǎo)入微分代數(shù)方程系統(tǒng)，然后利用分析結(jié)構(gòu)力學(xué)理論[9]進(jìn)行逐步積分，該積分保辛.該方法在積分點上嚴(yán)格滿足四元數(shù)長度等于1的約束條件，而在區(qū)段內(nèi)部則由插值近似，屬于祖沖之類方法[10].數(shù)值算例表明仿真效果優(yōu)異.

本文簡要介紹四元數(shù)和歐拉角的基本理論，以重陀螺為例對2種表示形式的數(shù)值積分進(jìn)行比較.對于歐拉角表述，應(yīng)用比較普遍的狀態(tài)方程表述.研究表明，以歐拉角和角速度為狀態(tài)變量形成的1階微分方程在使用差分近似積分時，精度損失很大，能量不守恒；該現(xiàn)象被周江華等[11]稱為“睡陀螺”.這是一個值得注意的問題，卻未得到廣泛關(guān)注；而采用文獻(xiàn)[8]給出的保辛格式，四元數(shù)單位長度約束條件得到滿足，仿真結(jié)果優(yōu)異，能量也達(dá)到守恒.

1剛體旋轉(zhuǎn)及其運動學(xué)表示

剛體運動由質(zhì)心平動和繞質(zhì)心轉(zhuǎn)動組成.如果剛體上有一固定點，則稱為剛體定點轉(zhuǎn)動問題.假設(shè)Oxyz為系統(tǒng)的慣性坐標(biāo)系，O為原點，亦為固定點.Ox′y′z′為隨體坐標(biāo)系，固定于剛體上.若將剛體的隨體坐標(biāo)軸取為與固定點有關(guān)的主軸，則剛體定點轉(zhuǎn)動可由式（1）描述.

2種數(shù)值積分的系統(tǒng)能量隨時間變化情況見圖4.在歐拉角表示中，雖然采用辛歐拉格式進(jìn)行數(shù)值積分，但是能量卻保持得不好，這驗證對狀態(tài)方程應(yīng)用辛歐拉格式并不能保辛.與之相反，四元數(shù)的數(shù)值積分保辛，其能量保持得很好，這也是保辛積分的優(yōu)勢所在.四元數(shù)保辛積分的約束誤差見圖5，表明在時間積分過程中單位長度約束條件滿足得很好.

4結(jié)束語

介紹2種剛體旋轉(zhuǎn)的數(shù)值積分，一種基于歐拉角表示，另一種基于四元數(shù)表示.以重陀螺的高速旋轉(zhuǎn)為例，對2種數(shù)值積分進(jìn)行比較發(fā)現(xiàn)：以歐拉角、角速度組成狀態(tài)變量，然后直接使用辛歐拉格式不

能保辛，能量衰減很快，數(shù)值積分存在缺陷；與之相

反，采用四元數(shù)表示，根據(jù)分析結(jié)構(gòu)力學(xué)的保辛積分方法，并結(jié)合祖沖之類方法的思想，可以很好地避免約束違約，仿真效果令人滿意，可作為陀螺等仿真分析的有力工具.

本文僅對以歐拉角、角速度組成狀態(tài)變量的數(shù)值積分進(jìn)行研究，對其他形式并未涉及，對其積分效果不佳的成因亦未研究，還有很多工作有待展開.

參考文獻(xiàn)：

[1]張樹俠， 孫靜. 捷聯(lián)式慣性導(dǎo)航系統(tǒng)[M]. 北京： 國防工業(yè)出版社， 1992： 4880.

[2]趙曉穎， 溫立書， 么彩蓮. 歐拉角參數(shù)表示下姿態(tài)的2階運動奇異性[J]. 科學(xué)技術(shù)與工程， 2012， 12（3）： 634637.

[3]黃雪樵. 克服歐拉方程奇異性的雙歐法[J]. 飛行力學(xué)， 1994， 12（4）： 2837.

[4]李躍軍， 閻超. 飛行器姿態(tài)角解算的全角度雙歐法[J]. 北京航空航天大學(xué)學(xué)報， 2007， 33（5）： 505508.

[5]NIKRAVESH P E， WEHAGE R A， KWON K. Euler parameters in computational kinematics and dynamics： Part 1[J]. J Mechanisms， Transmissions & Automation Des， 1985， 107（3）： 358365.

[6]BETSCH P， SIEBERT R. Rigid body dynamics in terms of quaternions： Hamiltonian formulation and conserving numerical integration[J]. Int J Numer Methods Eng， 2009， 79（4）： 444473.

[7]UDWADIA F E， SCHUTTE A D. An alternative derivation of the quaternion equations of motion for rigidbody rotational dynamics[J]. J Appl Mech， 2010， 77（4）： 44505.

[8]徐小明， 鐘萬勰. 剛體動力學(xué)的四元數(shù)表示及保辛積分[J]. 應(yīng)用數(shù)學(xué)和力學(xué)， 2014， 35（1）： 111.

[9]鐘萬勰， 高強(qiáng). 約束動力系統(tǒng)的分析結(jié)構(gòu)力學(xué)積分[J]. 動力學(xué)與控制學(xué)報， 2006， 4（3）： 193200.

[10]鐘萬勰， 高強(qiáng)， 彭海軍. 經(jīng)典力學(xué)辛講[M]. 大連： 大連理工大學(xué)出版社， 2013： 202241.

[11]周江華， 苗育紅， 李宏， 等. 四元數(shù)在剛體姿態(tài)仿真中的應(yīng)用研究[J].飛行力學(xué)， 2000， 18（4）： 2832.

[12]HAIRER E， LUBICH C， WANNER G. Geometric numerical integration： structurepreserving algorithms for ordinary differential equations[M]. Berlin： Springer， 2006： 189.

（編輯于杰）

摘要：為推廣四元數(shù)保辛積分在工程中的應(yīng)用，對歐拉角表示的狀態(tài)方程數(shù)值積分與四元數(shù)的保辛積分進(jìn)行比較.重陀螺的數(shù)值仿真結(jié)果表明四元數(shù)保辛積分的數(shù)值結(jié)果明顯優(yōu)于歐拉角狀態(tài)方程積分.與歐拉角狀態(tài)方程積分相比，四元數(shù)保辛積分在剛體動力學(xué)的數(shù)值仿真中更具優(yōu)勢.

關(guān)鍵詞：四元數(shù)； 歐拉角； 剛體動力學(xué)； 保辛； 重陀螺

中圖分類號： TP391.9； O313.3

文獻(xiàn)標(biāo)志碼： A

0引言

四元數(shù)、歐拉角和方向余弦[1]是描述剛體旋轉(zhuǎn)最主要的3種坐標(biāo)形式.方向余弦法需要9個參量，應(yīng)用較少；而另外2種則應(yīng)用廣泛，如航天器姿態(tài)控制和捷聯(lián)式慣性導(dǎo)航系統(tǒng)[1]等，對于兩者的研究也卷帙浩繁，但對兩者優(yōu)劣的評價卻褒貶不一.

描述剛體在三維空間中的運動姿態(tài)可采用2類12種歐拉角系統(tǒng)，分別對應(yīng)于不同的旋轉(zhuǎn)軸先后次序.目前公認(rèn)的用歐拉角描述旋轉(zhuǎn)的固有缺陷是奇異性問題[2]，即：無論采用哪種歐拉角系統(tǒng)，都不可避免地會含有奇異點，使得在該點附近區(qū)域進(jìn)行的數(shù)值積分精度不高.對于小角度旋轉(zhuǎn)，只要采用適當(dāng)?shù)臍W拉角系統(tǒng)便可避開奇異點；而在大角度旋轉(zhuǎn)時，若想避開奇異點，必須提供2套歐拉角系統(tǒng)交替進(jìn)行計算.據(jù)此，黃雪樵[3]提出一種“雙歐法”的計算方法；近幾年仍有學(xué)者[4]在繼續(xù)研究該方法.雙歐法雖然解決奇異性問題，但是計算過程較復(fù)雜.

四元數(shù)用于描述剛體旋轉(zhuǎn)，沒有奇異性問題，可很好地描述剛體的全角度旋轉(zhuǎn).然而，四元數(shù)需要滿足長度等于1的單位約束，這成為制約其應(yīng)用的限制.在實際應(yīng)用中，經(jīng)常采用的正交化修正等方法只能緩解長度的偏移，無法從根本上解決問題；黃雪樵[3]在其雙歐法中也有所討論.目前，對于單位約束最佳的解決方案是將四元數(shù)表示的剛體動力學(xué)方程導(dǎo)入微分代數(shù)方程范疇，近年來逐漸有學(xué)者[57]展開相關(guān)問題的研究.徐小明等[8]提出一種基于四元數(shù)理論描述剛體旋轉(zhuǎn)的保辛數(shù)值積分方法.該方法先將問題導(dǎo)入微分代數(shù)方程系統(tǒng)，然后利用分析結(jié)構(gòu)力學(xué)理論[9]進(jìn)行逐步積分，該積分保辛.該方法在積分點上嚴(yán)格滿足四元數(shù)長度等于1的約束條件，而在區(qū)段內(nèi)部則由插值近似，屬于祖沖之類方法[10].數(shù)值算例表明仿真效果優(yōu)異.

本文簡要介紹四元數(shù)和歐拉角的基本理論，以重陀螺為例對2種表示形式的數(shù)值積分進(jìn)行比較.對于歐拉角表述，應(yīng)用比較普遍的狀態(tài)方程表述.研究表明，以歐拉角和角速度為狀態(tài)變量形成的1階微分方程在使用差分近似積分時，精度損失很大，能量不守恒；該現(xiàn)象被周江華等[11]稱為“睡陀螺”.這是一個值得注意的問題，卻未得到廣泛關(guān)注；而采用文獻(xiàn)[8]給出的保辛格式，四元數(shù)單位長度約束條件得到滿足，仿真結(jié)果優(yōu)異，能量也達(dá)到守恒.

1剛體旋轉(zhuǎn)及其運動學(xué)表示

剛體運動由質(zhì)心平動和繞質(zhì)心轉(zhuǎn)動組成.如果剛體上有一固定點，則稱為剛體定點轉(zhuǎn)動問題.假設(shè)Oxyz為系統(tǒng)的慣性坐標(biāo)系，O為原點，亦為固定點.Ox′y′z′為隨體坐標(biāo)系，固定于剛體上.若將剛體的隨體坐標(biāo)軸取為與固定點有關(guān)的主軸，則剛體定點轉(zhuǎn)動可由式（1）描述.

2種數(shù)值積分的系統(tǒng)能量隨時間變化情況見圖4.在歐拉角表示中，雖然采用辛歐拉格式進(jìn)行數(shù)值積分，但是能量卻保持得不好，這驗證對狀態(tài)方程應(yīng)用辛歐拉格式并不能保辛.與之相反，四元數(shù)的數(shù)值積分保辛，其能量保持得很好，這也是保辛積分的優(yōu)勢所在.四元數(shù)保辛積分的約束誤差見圖5，表明在時間積分過程中單位長度約束條件滿足得很好.

4結(jié)束語

介紹2種剛體旋轉(zhuǎn)的數(shù)值積分，一種基于歐拉角表示，另一種基于四元數(shù)表示.以重陀螺的高速旋轉(zhuǎn)為例，對2種數(shù)值積分進(jìn)行比較發(fā)現(xiàn)：以歐拉角、角速度組成狀態(tài)變量，然后直接使用辛歐拉格式不

能保辛，能量衰減很快，數(shù)值積分存在缺陷；與之相

反，采用四元數(shù)表示，根據(jù)分析結(jié)構(gòu)力學(xué)的保辛積分方法，并結(jié)合祖沖之類方法的思想，可以很好地避免約束違約，仿真效果令人滿意，可作為陀螺等仿真分析的有力工具.

本文僅對以歐拉角、角速度組成狀態(tài)變量的數(shù)值積分進(jìn)行研究，對其他形式并未涉及，對其積分效果不佳的成因亦未研究，還有很多工作有待展開.

參考文獻(xiàn)：

[1]張樹俠， 孫靜. 捷聯(lián)式慣性導(dǎo)航系統(tǒng)[M]. 北京： 國防工業(yè)出版社， 1992： 4880.

[2]趙曉穎， 溫立書， 么彩蓮. 歐拉角參數(shù)表示下姿態(tài)的2階運動奇異性[J]. 科學(xué)技術(shù)與工程， 2012， 12（3）： 634637.

[3]黃雪樵. 克服歐拉方程奇異性的雙歐法[J]. 飛行力學(xué)， 1994， 12（4）： 2837.

[4]李躍軍， 閻超. 飛行器姿態(tài)角解算的全角度雙歐法[J]. 北京航空航天大學(xué)學(xué)報， 2007， 33（5）： 505508.

[5]NIKRAVESH P E， WEHAGE R A， KWON K. Euler parameters in computational kinematics and dynamics： Part 1[J]. J Mechanisms， Transmissions & Automation Des， 1985， 107（3）： 358365.

[6]BETSCH P， SIEBERT R. Rigid body dynamics in terms of quaternions： Hamiltonian formulation and conserving numerical integration[J]. Int J Numer Methods Eng， 2009， 79（4）： 444473.

[7]UDWADIA F E， SCHUTTE A D. An alternative derivation of the quaternion equations of motion for rigidbody rotational dynamics[J]. J Appl Mech， 2010， 77（4）： 44505.

[8]徐小明， 鐘萬勰. 剛體動力學(xué)的四元數(shù)表示及保辛積分[J]. 應(yīng)用數(shù)學(xué)和力學(xué)， 2014， 35（1）： 111.

[9]鐘萬勰， 高強(qiáng). 約束動力系統(tǒng)的分析結(jié)構(gòu)力學(xué)積分[J]. 動力學(xué)與控制學(xué)報， 2006， 4（3）： 193200.

[10]鐘萬勰， 高強(qiáng)， 彭海軍. 經(jīng)典力學(xué)辛講[M]. 大連： 大連理工大學(xué)出版社， 2013： 202241.

[11]周江華， 苗育紅， 李宏， 等. 四元數(shù)在剛體姿態(tài)仿真中的應(yīng)用研究[J].飛行力學(xué)， 2000， 18（4）： 2832.

[12]HAIRER E， LUBICH C， WANNER G. Geometric numerical integration： structurepreserving algorithms for ordinary differential equations[M]. Berlin： Springer， 2006： 189.

（編輯于杰）

摘要：為推廣四元數(shù)保辛積分在工程中的應(yīng)用，對歐拉角表示的狀態(tài)方程數(shù)值積分與四元數(shù)的保辛積分進(jìn)行比較.重陀螺的數(shù)值仿真結(jié)果表明四元數(shù)保辛積分的數(shù)值結(jié)果明顯優(yōu)于歐拉角狀態(tài)方程積分.與歐拉角狀態(tài)方程積分相比，四元數(shù)保辛積分在剛體動力學(xué)的數(shù)值仿真中更具優(yōu)勢.

關(guān)鍵詞：四元數(shù)； 歐拉角； 剛體動力學(xué)； 保辛； 重陀螺

中圖分類號： TP391.9； O313.3

文獻(xiàn)標(biāo)志碼： A

0引言

四元數(shù)、歐拉角和方向余弦[1]是描述剛體旋轉(zhuǎn)最主要的3種坐標(biāo)形式.方向余弦法需要9個參量，應(yīng)用較少；而另外2種則應(yīng)用廣泛，如航天器姿態(tài)控制和捷聯(lián)式慣性導(dǎo)航系統(tǒng)[1]等，對于兩者的研究也卷帙浩繁，但對兩者優(yōu)劣的評價卻褒貶不一.

描述剛體在三維空間中的運動姿態(tài)可采用2類12種歐拉角系統(tǒng)，分別對應(yīng)于不同的旋轉(zhuǎn)軸先后次序.目前公認(rèn)的用歐拉角描述旋轉(zhuǎn)的固有缺陷是奇異性問題[2]，即：無論采用哪種歐拉角系統(tǒng)，都不可避免地會含有奇異點，使得在該點附近區(qū)域進(jìn)行的數(shù)值積分精度不高.對于小角度旋轉(zhuǎn)，只要采用適當(dāng)?shù)臍W拉角系統(tǒng)便可避開奇異點；而在大角度旋轉(zhuǎn)時，若想避開奇異點，必須提供2套歐拉角系統(tǒng)交替進(jìn)行計算.據(jù)此，黃雪樵[3]提出一種“雙歐法”的計算方法；近幾年仍有學(xué)者[4]在繼續(xù)研究該方法.雙歐法雖然解決奇異性問題，但是計算過程較復(fù)雜.

四元數(shù)用于描述剛體旋轉(zhuǎn)，沒有奇異性問題，可很好地描述剛體的全角度旋轉(zhuǎn).然而，四元數(shù)需要滿足長度等于1的單位約束，這成為制約其應(yīng)用的限制.在實際應(yīng)用中，經(jīng)常采用的正交化修正等方法只能緩解長度的偏移，無法從根本上解決問題；黃雪樵[3]在其雙歐法中也有所討論.目前，對于單位約束最佳的解決方案是將四元數(shù)表示的剛體動力學(xué)方程導(dǎo)入微分代數(shù)方程范疇，近年來逐漸有學(xué)者[57]展開相關(guān)問題的研究.徐小明等[8]提出一種基于四元數(shù)理論描述剛體旋轉(zhuǎn)的保辛數(shù)值積分方法.該方法先將問題導(dǎo)入微分代數(shù)方程系統(tǒng)，然后利用分析結(jié)構(gòu)力學(xué)理論[9]進(jìn)行逐步積分，該積分保辛.該方法在積分點上嚴(yán)格滿足四元數(shù)長度等于1的約束條件，而在區(qū)段內(nèi)部則由插值近似，屬于祖沖之類方法[10].數(shù)值算例表明仿真效果優(yōu)異.

本文簡要介紹四元數(shù)和歐拉角的基本理論，以重陀螺為例對2種表示形式的數(shù)值積分進(jìn)行比較.對于歐拉角表述，應(yīng)用比較普遍的狀態(tài)方程表述.研究表明，以歐拉角和角速度為狀態(tài)變量形成的1階微分方程在使用差分近似積分時，精度損失很大，能量不守恒；該現(xiàn)象被周江華等[11]稱為“睡陀螺”.這是一個值得注意的問題，卻未得到廣泛關(guān)注；而采用文獻(xiàn)[8]給出的保辛格式，四元數(shù)單位長度約束條件得到滿足，仿真結(jié)果優(yōu)異，能量也達(dá)到守恒.

1剛體旋轉(zhuǎn)及其運動學(xué)表示

剛體運動由質(zhì)心平動和繞質(zhì)心轉(zhuǎn)動組成.如果剛體上有一固定點，則稱為剛體定點轉(zhuǎn)動問題.假設(shè)Oxyz為系統(tǒng)的慣性坐標(biāo)系，O為原點，亦為固定點.Ox′y′z′為隨體坐標(biāo)系，固定于剛體上.若將剛體的隨體坐標(biāo)軸取為與固定點有關(guān)的主軸，則剛體定點轉(zhuǎn)動可由式（1）描述.

2種數(shù)值積分的系統(tǒng)能量隨時間變化情況見圖4.在歐拉角表示中，雖然采用辛歐拉格式進(jìn)行數(shù)值積分，但是能量卻保持得不好，這驗證對狀態(tài)方程應(yīng)用辛歐拉格式并不能保辛.與之相反，四元數(shù)的數(shù)值積分保辛，其能量保持得很好，這也是保辛積分的優(yōu)勢所在.四元數(shù)保辛積分的約束誤差見圖5，表明在時間積分過程中單位長度約束條件滿足得很好.

4結(jié)束語

介紹2種剛體旋轉(zhuǎn)的數(shù)值積分，一種基于歐拉角表示，另一種基于四元數(shù)表示.以重陀螺的高速旋轉(zhuǎn)為例，對2種數(shù)值積分進(jìn)行比較發(fā)現(xiàn)：以歐拉角、角速度組成狀態(tài)變量，然后直接使用辛歐拉格式不

能保辛，能量衰減很快，數(shù)值積分存在缺陷；與之相

反，采用四元數(shù)表示，根據(jù)分析結(jié)構(gòu)力學(xué)的保辛積分方法，并結(jié)合祖沖之類方法的思想，可以很好地避免約束違約，仿真效果令人滿意，可作為陀螺等仿真分析的有力工具.

本文僅對以歐拉角、角速度組成狀態(tài)變量的數(shù)值積分進(jìn)行研究，對其他形式并未涉及，對其積分效果不佳的成因亦未研究，還有很多工作有待展開.

參考文獻(xiàn)：

[1]張樹俠， 孫靜. 捷聯(lián)式慣性導(dǎo)航系統(tǒng)[M]. 北京： 國防工業(yè)出版社， 1992： 4880.

[2]趙曉穎， 溫立書， 么彩蓮. 歐拉角參數(shù)表示下姿態(tài)的2階運動奇異性[J]. 科學(xué)技術(shù)與工程， 2012， 12（3）： 634637.

[3]黃雪樵. 克服歐拉方程奇異性的雙歐法[J]. 飛行力學(xué)， 1994， 12（4）： 2837.

[4]李躍軍， 閻超. 飛行器姿態(tài)角解算的全角度雙歐法[J]. 北京航空航天大學(xué)學(xué)報， 2007， 33（5）： 505508.

[5]NIKRAVESH P E， WEHAGE R A， KWON K. Euler parameters in computational kinematics and dynamics： Part 1[J]. J Mechanisms， Transmissions & Automation Des， 1985， 107（3）： 358365.

[6]BETSCH P， SIEBERT R. Rigid body dynamics in terms of quaternions： Hamiltonian formulation and conserving numerical integration[J]. Int J Numer Methods Eng， 2009， 79（4）： 444473.

[7]UDWADIA F E， SCHUTTE A D. An alternative derivation of the quaternion equations of motion for rigidbody rotational dynamics[J]. J Appl Mech， 2010， 77（4）： 44505.

[8]徐小明， 鐘萬勰. 剛體動力學(xué)的四元數(shù)表示及保辛積分[J]. 應(yīng)用數(shù)學(xué)和力學(xué)， 2014， 35（1）： 111.

[9]鐘萬勰， 高強(qiáng). 約束動力系統(tǒng)的分析結(jié)構(gòu)力學(xué)積分[J]. 動力學(xué)與控制學(xué)報， 2006， 4（3）： 193200.

[10]鐘萬勰， 高強(qiáng)， 彭海軍. 經(jīng)典力學(xué)辛講[M]. 大連： 大連理工大學(xué)出版社， 2013： 202241.

[11]周江華， 苗育紅， 李宏， 等. 四元數(shù)在剛體姿態(tài)仿真中的應(yīng)用研究[J].飛行力學(xué)， 2000， 18（4）： 2832.

[12]HAIRER E， LUBICH C， WANNER G. Geometric numerical integration： structurepreserving algorithms for ordinary differential equations[M]. Berlin： Springer， 2006： 189.

（編輯于杰）