3D刚体摆姿态动力学的离散变分数值方法研究

2.1 平面单摆 连续形式：设平面单摆的位形空间[4]为 S1 q 2 | 阵群，且矩阵的行列式均为 1，即 SO(2) R 中，平面单摆的拉格朗日动力学方程如下
（12）

q 1 。设 SO(2) 为一个 2 2 的特别正交矩

22
ˆ 0 - ，其 | R T R I 22 , det[ R ] 1 。定义 0
1 1 ˆ ˆ ˆ mgle Re L( R, ) ml 2 2 mgle2 Re2 ml 2 , 2 2 2 2
（13）
式中， g 为重力加速度； m 为摆的质量； l 为摆的长度。 e2 [0;1] 2 是指惯性坐标的重力的单位向量。 并且，旋转矩阵 R I 22 SO(2) 是指固定坐标系到惯性坐标系的一个转换矩阵。利用指数函数计算旋转 矩阵 R 及角速度 的变分如下

（34）
对 3D 摆作仿真实验，给出初始条件： m 1 kg ， [0, 0, 0.3] m ， J diag[0.13,0.28,0.17] kgm ，
2
R0 I ， 0 [4.14, 4.14, 4.14] rad/s ， h 0.01 seconds ， t f 100 seconds 。用离散变积分算法仿真
( k
（31）
h T T Fk mg Rk e3 ) h Fk J d J d FkT 2
（32） （33）
Rk 1 R k Fk
由于式(32)在仿真实验中所需的仿真时间比较长，所以利用指数函数及无穷小等价量可得出式（32） 的等价方程为
Fk e f k sin f k 1 cos f k h2 T h k mg Rk e3 Jf k f k Jfk 2 fk fk
3D 刚体摆姿态动力学的离散变分数值方法研究 1)
姜夏英*,2) 戈京 100192）
摘要：本文研究了 3D 刚体摆的姿态动力学数值计算方法，即利用离散变分的数值积分方法针对 3D 刚体摆的运动 学及动力学方程解的稳定性及其能量、角动量与范数的精度进行了数值计算研究。3D 摆是指在重力的作用下由一 个无摩擦的支撑点联结的非对称性刚体，是一个具有三自由度旋转摆系统。本文分别对 3D 摆的悬垂及倒置两种 特殊情况利用由离散变分原理推导出的数值积分方法进行求解，将其结果再和 4 阶龙格-库塔方法求解的结果进 行比较，并在计算机上实现数值仿真。从数值仿真实验中可以看出，离散变分的数值积分方法相比龙格-库塔方 法具有更高的精度。另外，本文还针对龙格-库塔方法提出了改进策略，即增强了动力学方程及运动学方程之间 的耦合性。结果发现其能量、角动量与范数的精度提高了很多，但和离散变分的数值积分方法相比，离散变分的 数值积分方法所得的精度仍然更高。因此可以表明离散变分的数值积分方法具有较高的计算精度和较好的保能 量，保持结构特性。 关键词：离散变分；龙格-库塔；姿态动力学；3D 刚体摆
k 0 j 1
N 1
N
= ( D1 Ld (qk , qk 1 , h) D2 Ld (qk 1 , qk , h)) qk
k 0
N 1
D2 Ld (qN 1 , qN , h) qN D2 Ld (q1 , q0 , h) q0
由离散积分及 q0 qN 0 ，并对任意的 qk ， S 0 可得到离散形式的拉格朗日动力学方程
D1 L(qk , qk 1 , h) D2 L(qk 1 , qk , h) 0
令 pk D1 Ld ( qk , qk 1 ) ，则可得到离散形式的哈密尔顿动力学方程
（11）
pk D1 Ld (qk , qk 1 ) pk 1 D2 Ld (qk , qk 1 ) 2.仿真计算
出 3D 摆的能量图、角动量图、范数图及角速度图。为了说明该算法的保结构，保能量特性，再用比较常 用的 4 阶龙格-库塔方法作比较。得出的角速度、角动量、能量及范数的仿真图如图 3、图 4、图 5、图 6 所示。
1 tr[ J ]I 33 J 。 2 h T h T T Fk mg Rk e3 mg Rk 1e3 2 2

再经过勒让德变换，由式（29）及式（30）可得到 3D 摆的离散形式的哈密尔顿动力学方程组如下
k 1 FkT k
U k 1 U k 1 1 T ( Fk 1 J d J d Fk J d FkT Rk 1 RkT1 ) 1 Fk 1 J d ) h ( h Rk 1 Rk 1
T
（29） （30）
Rk 1 Rk Fk
T 式中， U mge3 R , J d
S L ( R, )dt
对式（24）作变分，再根据变分原理 S 0 可得出 3D 摆的连续拉格朗日动力学方程如下
（24）
J mg RT e3 J
R ˆ R
（25） （26）
离散形式：设一个矩阵 Fk SO(3) 使得 Rk 1 Rk Fk ，且设离散时间步长 h 0 ，则 Rk SO(3) 表示在时 刻 t kh t0 时的旋转矩阵量。则可得到式（23）的离散形式为
（9）
对式（9）作变分得
S ( D1 Ld (qk , qk 1 , h) qk D2 Ld (qk , qk 1 , h) qk 1 )
k 0
N 1
（10）
= ( D1 Ld (qk , qk 1 , h) qk ) ( D2 Ld (q j 1 , q j , h) q j )
[1]
1.离散变分积分子介绍
1.1 连续动力学方程 将质点系的动能 T 与势能 V 之差定义为拉格朗日方程，即
L T V
式中， T
1) 2)
（1）
1 ， V V ( q ) ，由此可得拉格朗日方程为 qMq 2
国家自然科学基金(No.11072038)和北京市自然科学基金（Class B, No. KZ201110772039）重点项目资助 jiangshang21@
R
d ˆ R ˆ | 0 R exp d
（14） （15）
ˆ ˆ ˆ ˆ ˆ ˆ ˆ
则对式（13）作变分再利用哈密尔顿原理得出如下连续动力学方程组
g e Re 0 2 1 l
R ˆ R
（16） （17）
离散形式：上面已给出连续动力学方程，下面讨论其离散动力学形式。定义一个矩阵 Fk 使其满足如下 关系式
qk qk 1 q q k 1 k 。因此可得到离散拉格朗日 ，q 2 h
Ld (qk , qk 1 , h) hL(
对拉格朗日方程作离散积分求和，并记
qk qk 1 qk 1 qk , ) 2 h
（8）
S Ld (qk , qk 1 , h)
k 0
N 1
引言
经典动力学方程主要有三种形式：拉格朗日方程、哈密尔顿方程及牛顿方程，本论文主要讨论拉格朗 日形式及哈密尔顿形式。拉格朗日形式的动力学方程及哈密尔顿形式的动力学方程都有一个特有的几何性 质：密尔顿流形是辛结构的，在非保守力的作用下总能量总是守恒的，并且辛形式的动量也是守恒的。例 如，球面摆的位形空间是二维球，并且刚体姿态动力学的配置空间是李群结构，叫做特殊的正交群。本文 所讨论的离散变分积分子是由离散的哈密尔顿原理推导出来的。 哈密尔顿系统的辛几何算法分别是由冯康提出来的较新的求解动力学方程的求解问题的方法。哈密尔 顿体系是动力系统的重要体系，一切真实的、耗散可忽略不计的物理过程都可表示成哈密尔顿体系 ，它 的应用范围很广，包括结构生物学、半导体、超导、等离子体、天体力学、材料和偏微分方程等。且动力 学方程的应用及其广泛，其求解问题各个领域都会碰到，如航天器、机器人及卫星等。传统的数值积分算 法较好地解决了低阶微分方程的解的问题，但在计算高阶微分方程的数值解，特别是长时间的计算时，经 常有发散现象且伴有能量耗散。本文将离散变分积分子算法引入动力学计算。意在提高求解过程中的收敛 性及解的精确性。
图 1 角动量
图 2 总能量
2.2 3D 摆 连续形式： 3D 摆是指一刚体绕无摩擦的固定支点旋转，刚体受到恒重力(或万有引力)作用且具有三 个转动自由度[5]。3D 摆的连续型拉格朗日函数为
L ( R , )
1 T T R mge3 R 2
（23）
式中， R 为旋转变换矩阵， 为 3D 摆的角速度， m 为 3D 摆的质量， g 为重力加速度， 为支点到质心 的位置向量。对式（23）作积分，并记
Ld ( Rk , Fk )
对式（27）作和，并记
1 T T tr[(Fk I 33 ) J d (Fk I 33 )] mge3 Rk mge3 Rk 1 2h
N 1 k 0
（27）
S Ld ( Rk , Fk )
（28）
对式（28）作变分，再根据离散变分原理 S 0 可得出 3D 摆的离散形式的拉格朗日动力学方程如下
Rk 1 Rk Fk
（18）
ˆ R ( R R ) F I ，因此可得到离散形式的拉格朗日动力学方程 由此得出， h k k 1 k k 22
T
Ld ( Rk , Fk )
1 h h T T ml 2 Fk I 22 , Fk I 22 mgle2 Rk e2 mgle2 Rk 1e2 2h 2 2
（6）
q
1.2 离散变分积分子方法
H H ,p p q
（7）
以上给出了两种形式的连续型动力学方程，但在实际的计算当中，习惯上采用离散的形式。现在考虑 动力学方程的离散形式，设一个离散时间步长 h 0 ，且设位置向量的两个端点 q0 q (0) ， q1 q ( h) ， 且 k 表示第 k 个离散变量 qk q (0 kh) ，并设 q 方程