-
旋转矩阵R的微分
对于线性位移$x(t)$和线性速度$v(t)$，很容易得出$v(t)=\frac{d}{dt}x(t)$。那么旋转矩阵$R(t)$和角速度$\omega(t)$之间的关系是怎样的呢？显然$\dot{R}(t)$不等于$\omega(t)$，因为$R(t)$是一个矩阵，而$\omega(t)$是向量。为了回答这个问题，回忆一下矩阵$R(t)$的实际意义。我们知道$R(t)$的列向量就是其坐标系在世界坐标系中的方向向量，因此$\dot{R}(t)$的列向量描述了刚体旋转时$x$，$y$，$z$轴的“速度”。下图描述了一个刚体以角速度$\omega(t)$进行旋转，考虑刚体上某一位置向量$r(t)$，其速度很容易推导出为：$$\dot{r}(t)=\omega(t) \times r(t)$$
$t$时刻，刚体$x$轴的方向向量在世界坐标系下的坐标为矩阵$R(t)$的第一列：$(r_{xx},r_{xy},r_{xz})^T$，则在$t$时刻$R(t)$第一列的微分就是该向量的变化率。根据上面的推论，我们可得到该微分为：$\omega(t) \times (r_{xx},r_{xy},r_{xz})^T$
对$R(t)$的另外两个列向量情况也一样。因此，我们可以写出：
$$\dot{R(t)}=\begin{pmatrix}
\omega(t) \times\begin{pmatrix} r_{xx}\\ r_{xy}\\ r_{xz} \end{pmatrix}
& \omega(t) \times\begin{pmatrix} r_{yx}\\ r_{yy}\\ r_{yz} \end{pmatrix}
& \omega(t) \times\begin{pmatrix} r_{zx}\\ r_{zy}\\ r_{zz} \end{pmatrix}
\end{pmatrix}$$
根据角速度$\omega$定义斜对称矩阵（skew-symmetric matrix）$\Omega$，斜对称矩阵满足$\Omega^T=-\Omega$：
$$\Omega=\begin{pmatrix}
0 & -\omega_z & \omega_y\\
\omega_z & 0 & -\omega_x\\
-\omega_y & \omega_x & 0
\end{pmatrix}$$
可以证明，对任意向量$P$均有：$$\omega \times P = \Omega P$$
则旋转矩阵$R(t)$的微分可以用角速度的斜对称矩阵表示为：$$\dot{R}(t)=\Omega R(t)$$
-
惯性参考系中的Euler方程
通常Euler方程中各量均为在与刚体固连的坐标系中表达，但这并不意味着此方程只在与刚体固连的坐标系中才成立。事实上很容易证明：对任一坐标系均有$$I \dot{\omega}+\omega \times I \omega =N$$
式中各量均为在上述给定坐标系中的表达式。当然，将Euler方程表示在与刚体固连坐标系中尤其特殊的优点，即惯性张量矩阵是常值矩阵，它可以预先计算或辨识出来，这给Euler方程的应用带来很大方便。如果欧拉方程在惯性坐标系中描述，那么惯性张量是个随着刚体运动而变化的变量，需要实时计算该值。
下面从动量矩定理开始推导惯性坐标系下的欧拉方程。等式中物理量下标$I$表示该量在惯性参考系（Inertia frame）中描述，这里刚体的惯性张量$I$不再恒定，而是随时间变化的量，因此在求导时要注意：
$$\boldsymbol\tau_I = \frac d{dt} (\mathrm I_I\, \boldsymbol \omega_I) =
\mathrm I_I \, \dot {\boldsymbol \omega}_I
+ \dot {\mathrm I}_I \,\boldsymbol \omega_I$$
惯性参考系下的惯性张量$I_I$和刚体局部参考系下的惯性张量$I_b$（$I_b$是常量）之间的关系如下：$$I_I=RI_bR^T$$
其中，矩阵$R$是刚体局部坐标系相对于惯性坐标系（世界坐标系）的变换，即$^I_bR$。考虑惯性张量$I_I$对时间的导数：$\dot{I}_I=\dot{R}I_bR^T+RI_b \dot{R}^T$
变换矩阵$R$的导数可以表示为：$\dot{R}=Sk(\omega_I)R$，$Sk(\omega_I)$是根据角速度$\omega_I$构造的斜对称矩阵。则有：$\dot{I}_I=Sk(\omega_I)RI_bR^T+RI_bR^TSk^T(\omega_I)=Sk(\omega_I)I_I-I_ISk(\omega_I)$
于是$\dot{I}_I\omega_I=Sk(\omega_I)I_I\omega_I-I_ISk(\omega_I)\omega_I=\omega_I \times (I_I\omega_I)-I_I(\omega_I \times \omega_I)=\omega_I \times (I_I\omega_I)$
因此，惯性参考系中的欧拉公式为：
$$\boldsymbol\tau_I =
\mathrm I_I \, \dot {\boldsymbol \omega}_I
+ \boldsymbol \omega_I \times (\mathrm I_I \,\boldsymbol \omega_I)$$
可以看出，惯性参考系下欧拉公式与刚体局部坐标系下欧拉公式形式一致，只是物理量的参照基准不同。
-
陀螺力矩
欧拉方程$I \dot{\omega}+\omega \times I\omega=\tau$中的$\omega \times I\omega$项在物理引擎中被称为陀螺力矩（gyroscopic torque），因为其有力矩的量纲。陀螺力矩的产生是由于角速度与角动量方向不一致，$\omega \times I\omega$叉乘结果不为零。
在某些游戏物理引擎中为了避免复杂计算影响实时性和稳定性，通常会忽略陀螺力矩，即直接使用公式：$$I\dot{\omega}=\tau$$
那么对于物理引擎中的刚体转动1中描述的那个问题，添加$\omega \times I\omega$项后准确的计算如下面Mathematica代码所示：
画出结果，可以发现与bullet2.83、ODE物理引擎中模拟的结果一致（仿真时间设为5s）：
bullet物理引擎从2.83版本开始才默认开启了陀螺力矩的计算。
// bullet 2.83 enum btRigidBodyFlags { BT_DISABLE_WORLD_GRAVITY = 1, ///BT_ENABLE_GYROPSCOPIC_FORCE flags is enabled by default in Bullet 2.83 and onwards. ///and it BT_ENABLE_GYROPSCOPIC_FORCE becomes equivalent to BT_ENABLE_GYROSCOPIC_FORCE_IMPLICIT_BODY ///See Demos/GyroscopicDemo and computeGyroscopicImpulseImplicit BT_ENABLE_GYROSCOPIC_FORCE_EXPLICIT = 2, BT_ENABLE_GYROSCOPIC_FORCE_IMPLICIT_WORLD=4, BT_ENABLE_GYROSCOPIC_FORCE_IMPLICIT_BODY=8, BT_ENABLE_GYROPSCOPIC_FORCE = BT_ENABLE_GYROSCOPIC_FORCE_IMPLICIT_BODY, }; // bullet 2.82 enum btRigidBodyFlags { BT_DISABLE_WORLD_GRAVITY = 1, ///The BT_ENABLE_GYROPSCOPIC_FORCE can easily introduce instability ///So generally it is best to not enable it. ///If really needed, run at a high frequency like 1000 Hertz: ///See Demos/GyroscopicDemo for an example use BT_ENABLE_GYROPSCOPIC_FORCE = 2 };