本文是一份用于游戏物理引擎的分析力学快速入门文章。
主要是为了编写物理引擎而记录的相关知识。和物理引擎关联不大的知识不做记录（混沌系统，相对论，哈密顿力学等）。
分析力学的教程推荐看B站任延宇的《理论力学》1。书籍则是Goldstain的《Classical Mechanics》2。但GoldStain的书写的一拖四，下定义的时候从来不说定义，并且刚体动力学那一章说的一坨。建议先看视频搞清楚整体脉络再看书。
另外推荐一本Taylor的《Classical Mechanics》3。是一本很老的书，补充了很多牛顿力学中的知识。就是内容太啰嗦，适合完全没学过物理的人看。可以没事看看用来查漏补缺。
质点运动学与动力学#
单质点运动学与动力学#
定义空间中的原点$O$。这样，$O$到空间中任意一点的位置构成一个矢量，称为位矢(radius vector)$\pmb{r}$，其表明了物体的位置。
而位置对时间的导数是速度：
$$ \pmb{v} = \frac{d\pmb{r}}{dt} = \dot{\pmb{r}} $$速度对时间的导数是加速度：
$$ \pmb{a} = \frac{d\pmb{v}}{dt} = \dot{\pmb{v}} = \ddot{\pmb{r}} $$牛顿第二定律表示：
$$ \pmb{F} = m\pmb{a} $$其实有一个等价的，使用动量描述的公式：
$$ \begin{aligned} \pmb{p} &= m\pmb{v} \\ \pmb{F} &= \dot{\pmb{p}} \end{aligned} $$定义角动量为：
$$ \pmb{L} = \pmb{r} \times \pmb{p} $$这也称为角动量矩。因为矩的定义是$\pmb{r} \times \cdot$
所以力矩就是$\pmb{N} = \pmb{r} \times \pmb{F}$
而力矩与角动量的关系则是：
$$ \pmb{N} = \dot{\pmb{L}} $$注意：单质点下的角动量，是指质点绕着$O$点移动的运动度量。单个质点本身显然没有旋转。
功与能量#
功是能量变化的度量。做功会传递能量。
功的计算式为力沿着路径的第二类曲线积分：
$$ W = \int_{L} \pmb{F} d\pmb{r} $$注意：以下情况会导致能量变化：
- 做功
- 热传递
- 质量迁移：物质带着能量进出系统
- 辐射：电磁波将能量带入/出
做功只是其中一种方式。只是在经典力学中，只研究功对能量的改变。
动能：通过功的定义式直接得出：
$$ \begin{aligned} W = \int_1^2 \pmb{F} d\pmb{r} = m\int_1^2 \frac{d\pmb{v}}{dt}\cdot \pmb{v}dt = \frac{1}{2}m(v_2^2 - v_1^2) \end{aligned} $$这里$T = \frac{1}{2}mv^2$就是动能(KE - Kinetic Energy)。而这条公式就是动能定理，即：
功是动能的改变量
接下来是 保守力(Conservation Force) 的定义：
沿环路做功为0的力为保守力，即：
$$ \oint_L \pmb{F_{cons}}d\pmb{r} = 0 $$
而 势能(PE - Potential Energy) 则是由保守力定义的：
势能是只和位置有关的函数，用$V(\pmb{r})$表示。保守力是势能的梯度：
$$ \pmb{F}_{cons} = -\nabla{V(\pmb{r})} $$
研究势能时，必须先选取势能零点(zero level of PE)，这是势能为0的点。
势能零点可以任意选取。比如重力势能的势能零点一般就选取在地球表面。
由势能定义我们知道：
$$ \begin{aligned} & \pmb{F}_{cons} = -\frac{\partial V}{\partial d\pmb{s}} \\ & \Rightarrow \\ & \pmb{F}_{cons} d\pmb{r} = -dV \\ W &= \int \pmb{F}_{cons}d\pmb{r} = V_1 - V_2 \end{aligned} $$所以保守力做正功，势能反而会减少。
只有保守力做功的系统称为保守系统(Conservation System)。
在保守系统中，保守力做功同时也会改变动能。这意味着动能的改变量和势能的改变量是一样的，于是有：
$$ \begin{aligned} & T_2 - T_1 = V_1 - V_2 \\ & \Rightarrow \\ & T_1 + V_1 = T_2 + V_2 \end{aligned} $$这意味着任意时刻，动能和势能的总和是不变的。这就是机械能守恒定律。
机械能就是动能加势能。
也就是说：只有保守力做功时，机械能守恒。
注意此时系统内能量依然有流动（如果保守力做正功，能量从势能流到动能），只是总能量不变。
随时间改变的势能#
有些势能是随时间改变的（$V(\pmb{r}, t)$）。比如两个电荷，其中一个电荷的电量随时间流逝，那么另一个电荷受到的电势能也会随之改变。
在经典力学中，只考虑不随时间改变的势能。
守恒律#
除了能量守恒还有如下守恒率：
合力为0($\sum \pmb{F}_i = 0$)时，动量守恒
合力矩为0($\sum \pmb{N}_i = 0$)时，角动量守恒
守恒律的意义#
物理学中，绝大部分时间都在求解微分方程。而守恒率可以帮助我们降低方程的阶数，从而更方便地求解。
守恒率比起牛顿定律适用范围更广。比如牛顿第三定律在电场中就不正确：
假设有垂直向上的电场，电场中有两个电荷沿着互相垂直的方向移动。那么电场给两个电荷的力就不满足牛顿第三定律。但是动量守恒依旧满足。
不过经典力学也不研究场。
多质点动力学#
多质点下，主要研究 N 个质点在内力 + 外力作用下的整体运动，核心是把复杂系统"简化 + 找守恒量“。
定义质心为：
$M = \frac{\sum m_i \pmb{r}_i}{\sum m_i}$
由于牛顿第三定律，系统内部的作用力可以两两抵消，动量守恒变成：
合外力为0时，动量守恒
而角动量则变成质心本身的角动量+相对质心的质点的角动量，变成：
合力矩为0时，角动量守恒
动能则变成：
$T = \frac{1}{2}Mv_{cm}^2 + \frac{1}{2}\sum m_i v_i^{‘2}$
其中$v_i^{’}$是质点相对质心的速度。即动能变成质心本身的动能+其他质点相对质心的动能。
而势能则相对复杂，一般只是对所有质点自己的势能求和得到系统势能。
因为势能分为：
- 内势：系统内各个质点之间的势能
- 外势能：整个系统的势能。如果在均匀场且耦合正比于质量，可以视为质心的势能，否则不行
所以没有一个通用计算方式，只能逐点势能计算后求和。
而机械能守恒依旧是：
系统只有保守力做功时机械能守恒
分析力学#
指用广义坐标、能量等标量函数和变分原理来研究力学的一整套方法，与"用矢量和力"的牛顿力学相对。
使用数学分析的方式重构牛顿力学。本质上还是牛顿力学，但提供了更好的数学工具，这些工具有更好的普适性（在非牛顿力学，比如量子力学，相对论等中都有作用）。
现在的教程中一般会按照这个路径去
广义坐标和约束#
广义坐标#
能唯一确定系统位置的一组任意参数，这些参数彼此独立
比如在笛卡尔坐标系中自由运动的小球，其广义坐标就是笛卡尔坐标中的$x,y,z$。
但在有约束的情况下（比如一个滑块从斜面上滑下来），其不需要使用整个笛卡尔坐标系表示。因为滑块总是紧贴斜面，所以完全可以用“滑块到斜面顶端的距离”这一个量就能表示滑块的位置（即使用图中绿色坐标轴建系）。
或者另一种情况：一个硬杆的一段被固定在天花板上。此时只需要杆到天花板的角度$\theta$即可描述。
也就是说，在很多时候，描述系统的坐标系可能不是固定的（可能是球坐标系，笛卡尔坐标系或者甚至是一维的），甚至量也不是同一的（比如上面例子中，一个是距离一个是角度）。这种就叫做广义坐标。
广义坐标使用$q$表示。一组广义坐标就是$(q_1,q_2, \cdots, q_i)$。注意：$q_i$之间不需要在一个坐标系下（比如$q_1$是角度，$q_2$是长度，$q_3$又是其他什么量）。只要这些量组合起来能唯一确定系统的位置就行。
那么广义坐标对时间的导数就是广义速度，记为$(\dot{q_1}, \dot{q_1}, \cdots, \dot{q_i})$。同理还有广义加速度。
而“彼此独立”是指一个广义坐标不能被其他广义坐标表示。其实就是不能有冗余坐标。
约束#
对系统运动的限制。它规定哪些位置/速度是可能的，哪些是不可能的
比如上面例子中“物体总是紧贴斜面”就是一种约束。而“连杆的一头固定在天花板”又是一种约束。
约束方程：
即约束的方程，其固定了一些广义坐标的取值范围
约束有很多分类方法：
- 完整(holonomic)/非完整(nonholonomic)约束：只和位置有关，可能和时间有关的约束$f(q_1,q_2,\cdots,q_i, t)$。这也是我们最常见，最需要研究的约束。 本质上来说，完整约束就是此约束对应的微分方程可解。非完整约束就是不可解。
- 定常/时变约束：是否和时间相关
- 几何约束/微分约束：
- 几何约束就是只和位置有关，可能和时间有关（其实就是完整约束）
- 微分约束就是还和速度有关。此时是微分方程，如果此方程不可解就是非完整约束，否则其实可以变成几何约束，所以依旧是完整约束
其实只需要记住完整/非完整约束即可。
拉格朗日方程#
历史上有两个拉格朗日方程：
- 力学的拉格朗日方程：直接从虚功原理和达郎贝尔原理出发推导出来，没有变分法
- 欧拉-拉格朗日方程：最早是欧拉和拉格朗日通过变分法发展出的一个数学方程，没什么物理意义。但后来哈密顿发展他的力学时发现，从最小作用量原理可以直接得到这个方程，赋予了他物理意义。
只是后来发现1中的方程就是2中的特例，所以后面统一叫欧拉-拉格朗日方程。
并且现在的教程路径都是先说哈密顿原理，然后直接从哈密顿原理，经过变分法推导出欧拉-拉格朗日方程。
虚功原理#
是处理静力学的一种方法。
首先定义虚位移：
一个假象的，很小的位移。用$\delta \pmb{r}$表示
注意：虚位移并没有真的移动物体。而是假设移动了一点点。
那么，在虚位移上做的功就是虚功。整个系统的虚功和就是：
$$ \delta W = \sum{(\pmb{F}_a + \pmb{F}_c}) \cdot \delta \pmb{r} $$这里$\pmb{F}_c$是约束力，$\pmb{F}_a$是主动力（即除了约束力之外的力）。
接下来定义理想约束：
“约束力在系统允许的任意虚位移上做的总功为零"的约束。即：
$$ \sum \pmb{F}_c \cdot \delta \pmb{r} = 0 $$
这就是说，要么约束力垂直于允许的位移，要么总约束力大小刚好抵消。
接下来定义理想系统：
只有理想约束的系统称为理想系统
那么就可以定义虚功原理：
理想约束下的静止系统保持平衡的充要条件为主动力在虚位移上所做虚功总和为零，即
$$ \sum \pmb{F}_a \cdot \delta \pmb{r} = 0 $$
即无论你往哪个方向挪动一点点物体，物体的主动力都不作功，那么此时系统肯定是平衡的（静止或做匀速直线运动，不过一般只关注静止情况）。
接下来了解广义坐标下的虚功原理。对于任意广义坐标$\pmb{r} = \pmb{r}(q_1, q_2, \cdots, q_n)$，其虚位移为：
$$ \delta \pmb{r} = \sum \frac{\partial \pmb{r}}{\partial q_i} \delta q_i $$求虚位移其实只是对$\pmb{r}$进行全微分。
那么虚功原理变为：
$$ \sum (\pmb{F}_a \cdot \frac{\partial \pmb{r}}{\partial q_i}) \delta q_i = 0 $$即
$$ \pmb{F}_a \cdot \frac{\partial \pmb{r}}{\partial q_i} = 0 $$D’Alembert原理（达朗伯原理）#
其实只是数学上的一个变换。由牛顿第二定律$\pmb{F} = \dot{\pmb{p}}$，我们将等式右边挪到左边就有：
$$ \pmb{F} - \dot{\pmb{p}} = 0 $$那么就将一个动力学问题变成了一个静力学问题。那么利用虚功原理，就有：
$$ \begin{aligned} & \sum{(\pmb{F} - \dot{\pmb{p}})} \cdot \delta \pmb{r} = 0 \\ & \Rightarrow \\ & \sum{(\pmb{F}_a + \pmb{F}_c - \dot{\pmb{p}})} \cdot \delta \pmb{r} = 0 \\ & \Rightarrow \\ & \sum{(\pmb{F}_a - \dot{\pmb{p}})} \cdot \delta \pmb{r} = 0 \\ \end{aligned} $$拉格朗日方程#
从达郎贝尔原理可以推导出拉格朗日方程。主要的公式是：
$$ \begin{cases} & \sum{(\pmb{F}_a - \dot{\pmb{p}})} \cdot \delta \pmb{r} = 0 \\ & \delta \pmb{r} = \sum \frac{\partial \pmb{r}}{\partial q_i} \delta q_i \\ & \pmb{v} = \sum \frac{\partial \pmb{r}}{\partial q_i}\dot{q_i} + \frac{\partial \pmb{r}}{\partial t} \end{cases} $$利用这些公式联立，最后可以得到拉格朗日方程：
$$ \begin{aligned} & \frac{d}{dt}(\frac{\partial T}{\partial \dot{q_i}}) - \frac{\partial T}{\partial q_i} = Q_i \\ & Q_i = \sum \pmb{F}_i\cdot \frac{\partial \pmb{r}}{\partial q_i} \end{aligned} $$其中$T$是动能。
而如果系统只有保守力，那么就有$\pmb{F}_i = -\nabla_i V(\pmb{r})$。那么就有：
$$ Q_i = -\sum \nabla_i V(\pmb{r}) \cdot \frac{\partial \pmb{r}}{\partial q_i} = -\frac{\partial V}{\partial q_i} $$那么移到拉格朗日方程左边，就有：
$$ \frac{d}{dt}(\frac{\partial T}{\partial \dot{q_i}}) - \frac{\partial (T - V)}{\partial q_i} = 0 $$而我们假设势能不随时间改变，所以左边的$\frac{d}{dt}\frac{\partial T}{\partial \dot{q}_i}$也可以加一个$V$得到$\frac{d}{dt}\frac{\partial (T-V)}{\partial \dot{q}_i}$。那么整个式子就是：
$$ \frac{d}{dt}(\frac{\partial (T - V)}{\partial \dot{q_i}}) - \frac{\partial (T - V)}{\partial q_i} = 0 $$令$L = T -V$，最后的拉格朗日方程就是：
$$ \frac{d}{dt}(\frac{\partial L}{\partial \dot{q_i}}) - \frac{\partial L}{\partial q_i} = 0 $$
其中$L$称为“拉格朗日量”。
注意：拉格朗日方程有两个形式，一个是带广义力$Q_i$的，一个是不带广义力的简洁版。
另外，利用“循环坐标”可以从拉格朗日方程推导出动量守恒。而当$L$不显含时间$t$的话，还能推导出能量守恒。这里不再赘述。
哈密顿原理与欧拉-拉格朗日方程#
哈密顿原理#
在固定时间区间内，系统从初始位形$q(t_1)$到最终位形$q(t_2)$的运动路径，是使下面作用量泛函取驻值(stationary value)的结果：
$$ S = \int_{t_1}^{t_2} L(q, \dot{q}, t)dt $$
所谓泛函就是：输入是函数，输出是数值的函数。
其中函数$S$就是作用量。
此定理也叫最小作用量原理。
哈密顿原理是一种假设。是从最速降线问题进行一般化推广而来。
此原理不仅作用于经典力学。对于任意的领域，只要能找到对应的函数$L$，就能直接套这个方法解决。
变分法和欧拉-拉格朗日方程#
变分法是一种找到哈密顿原理中作用量$S$的驻点的一种方法。本质上是数学上的一种手段。
驻点，从实函数来说，就是指函数的导数/偏导数为0的点。或者更一般地，就是让函数$f(x_1, x_2, \cdots, x_n)$的微分$df=0$的点。
那么，对于泛函也是一样的操作：我们需要找到某个函数，让作用量$S$的微分$\delta S = 0$。
那么对于泛函$S[y(x)]$，我们可以找一个距离他很接近的函数$\bar{y(x)}$，然后做差来得到他的微分。不过由于$S$是泛函，我们需要换一个称呼：“变分”，用$\delta$表示：
“变分"主要是指”函数（整条曲线/路径）的无穷小改变"——即 $\delta y$。它和"微分"是平行的两件事：微分管"自变量的改变”，变分管"函数的改变"。
而$\delta S$就是“作用量的变分“。
$$ \begin{aligned} \delta S &= S[\bar{y(x)}] - S[(y(x))] \\ & = \int_{x_1}^{x_2}L(x, \bar{y},\bar{y}^{'})dx - \int_{x_1}^{x_2}L(x, y, y^{'})dx = 0 \end{aligned} $$
注：研究变分的过程中，$x$没有发生改变（是固定值），即$\delta x = 0$。
同时，显然在$x_1$和$x_2$处，两函数值是相同的（因为起始状态和最终状态肯定是确定的），所以有：
$$ \delta y(x_1) = \delta y(x_2) = 0 $$接下来，对$S[\bar{y(x)}]$做多元函数泰勒展开，展开到一阶（其余高阶项可以忽略）：
$$ S[\bar{y(x)}] = \int_{x_1}^{x_2}(L(x, y, y^{'}) + \frac{\partial L}{\partial y}\delta y + \frac{\partial L}{\partial y^{'}}\delta y^{'})dx $$这样，$S[\bar{y(x)}] - S[y(x)]$就变成：
$$ \delta S = \int_{x_1}^{x_2} (\frac{\partial L}{\partial y}\delta y + \frac{\partial L}{\partial y^{'}}\delta y^{'})dx = 0 $$数学上可以证明，变分$\delta y^{’}=\frac{d}{dx}(\delta y)$，所以有：
$$ \delta S = \int_{x_1}^{x_2} (\frac{\partial L}{\partial y}\delta y + \frac{\partial L}{\partial y^{'}}\frac{d}{dx}(\delta y))dx = 0 $$然后将加号后面那一半变成：
$$ \frac{d}{dx}(\frac{\partial L}{\partial y^{'}}\delta y) - \delta y \frac{d}{dx}(\frac{\partial L}{\partial y^{'}}) $$这样$\delta S$就是：
$$ \begin{aligned} \delta S &= \int_{x_1}^{x_2} (\frac{\partial L}{\partial y}\delta y + \frac{d}{dx}(\frac{\partial L}{\partial y^{'}}\delta y) - \delta y \frac{d}{dx}(\frac{\partial L}{\partial y^{'}}))dx \\ &= \int_{x_1}^{x_2} \frac{d}{dx}(\frac{\partial L}{\partial y^{'}}\delta y) dx + \int_{x_1}^{x_2}(\frac{\partial L}{\partial y}\delta y - \delta y \frac{d}{dx}(\frac{\partial L}{\partial y^{'}}))dx \\ &= \left. (\frac{\partial L}{\partial y^{'}}\delta y)\right|_{x_1}^{x_2} + \int_{x_1}^{x_2} \delta y(\frac{\partial L}{\partial y} - \frac{d}{dx}(\frac{\partial L}{\partial y^{'}}))dx \\ \end{aligned} $$而由于$y(x)$在$x_1, x_2$处的变分为0。所以$\left. (\frac{\partial L}{\partial y^{’}}\delta y)\right|_{x_1}^{x_2} = 0$。
而$\delta y$是两函数变分，这个东西不一定恒取0。所以为了让$\delta S = 0$，必须是：
$$ \frac{\partial L}{\partial y} - \frac{d}{dx}(\frac{\partial L}{\partial y^{'}}) = 0 $$
这个公式就是欧拉-拉格朗日公式。
其与前面说的拉格朗日公式一模一样。只需要取$L = T-V$就变成之前的那个公式。但欧拉-拉格朗日公式并不规定$L$的形式，这个公式有更好的普适性。
哈密顿力学#
TODO：未学习，感觉和物理引擎关联不大，等到需要时再补充。
刚体动力学#
物体内任意两点之间距离永远不变的物体称为刚体
刚体自由度#
刚体动力学主要就是在研究刚体的旋转方程。
首先必须明确，对于一个刚体，需要多少自由度去定义（默认笛卡尔坐标系下）：
使用刚体中任意三个不共线的质点$P_1, P_2, P_3$即可确定刚体位置。其约束方程为：
$$ \begin{aligned} & (p_1 - p_2) ^ 2 = d_1^2 \\ & (p_2 - p_3) ^ 2 = d_2^2 \\ & (p_3 - p_1) ^ 2 = d_3^2 \\ \end{aligned} $$产生了3个约束方程。
而其余$N-3$个点满足方程：
$$ \begin{aligned} & (p_i - p_1)^2 = a_i^2 \\ & (p_i - p_2)^2 = b_i^2 \\ & (p_i - p_3)^2 = c_i^2 \\ \end{aligned} $$产生了$3(N-3)$个约束方程。那么系统总约束方程个数就是 $3N - 6$。
而系统中一共有$N$个质点，其需要$3N$个自由度描述，所以显然刚体自由度就是
$$ 3N - (3N - 6) = 6 $$而我们知道，刚体的位置可以通过刚体中任意一点的位置来表示，占据3个自由度。那么剩下3个自由度就是刚体旋转的自由度。
但究竟怎么描述旋转，其实有很多方法（Proper Eular angles, Axis-Angle, Tait–Bryan angles）。后面会一一介绍。
接下来只研究刚体转动的情况。假设刚体无位移。
一般转动#
我们假设刚体以他的质心为原点进行转动。那么转动前的坐标记为$\pmb{i}, \pmb{j}, \pmb{k}$，而转动后的记为$\pmb{i^{’}}, \pmb{j^{’}}, \pmb{k^{’}}$。那么，我们可以对转动后的任意一个轴，使用其相对于转动前的轴的向量余弦来表示：
$$ \begin{aligned} & \cos(\pmb{i^{'}},\pmb{i}) = \pmb{i^{'}}\cdot \pmb{i} \\ & \cos(\pmb{i^{'}},\pmb{j}) = \pmb{i^{'}}\cdot \pmb{j} \\ & \cos(\pmb{i^{'}},\pmb{k}) = \pmb{i^{'}}\cdot \pmb{k} \\ \end{aligned} $$左右两边同时乘上$\pmb{i}, \pmb{j}, \pmb{k}$，就可以得到$\pmb{i^{’}}$相对于$\pmb{i},\pmb{j},\pmb{k}$的关系：
$$ \begin{aligned} & \cos(\pmb{i^{'}},\pmb{i}) \pmb{i} = \pmb{i^{'}} \\ & \cos(\pmb{i^{'}},\pmb{j}) \pmb{j} = \pmb{i^{'}} \\ & \cos(\pmb{i^{'}},\pmb{k}) \pmb{k} = \pmb{i^{'}} \\ \end{aligned} $$再将三式相加，就可以用旧坐标表示新坐标：
$$ 3\pmb{i^{'}} = \cos(\pmb{i^{'}},\pmb{i}) \pmb{i} + \cos(\pmb{i^{'}},\pmb{j}) \pmb{j} + \cos(\pmb{i^{'}},\pmb{k}) \pmb{k} $$$\pmb{j^{’}},\pmb{k^{’}}$同理。这样就得到新坐标和旧坐标的关系式（用$\alpha_{ij}$来表示系数）：
$$ \begin{aligned} & \pmb{i^{'}} = \alpha_{11} \pmb{i} + \alpha_{12}\pmb{j} + \alpha_{13} \pmb{k} \\ & \pmb{j^{'}} = \alpha_{21} \pmb{i} + \alpha_{22}\pmb{j} + \alpha_{23} \pmb{k} \\ & \pmb{k^{'}} = \alpha_{31} \pmb{i} + \alpha_{32}\pmb{j} + \alpha_{33} \pmb{k} \\ \end{aligned} $$那么 $\alpha_{ij}$构成一个矩阵。这个矩阵就是表示旋转的矩阵，记为$A$。
那么对于刚体中任意点$\pmb{r}$的旋转，都可以使用此线性变换得到：
$$ \pmb{r}^{'} = A\pmb{r} $$那么至于要怎么表示这个旋转，历史上有多个方式：
欧拉角（Proper Eular Angles）#
注意：此欧拉角是欧拉最开始发明的表示角的方式。在其他领域中，欧拉角的含义各不相同。
欧拉角是使用三个角度表示的角：
- 进动角$\phi$(precession)：刚体绕自己的转轴（即$\pmb{k}$）旋转的角度
- 章动角$\theta$(nutation)：在进动之后，刚体绕旋转过的$\pmb{i}$轴转动的角度
- 自转角$\psi$(intrinsic)：在章动之后，刚体绕旋转过的$\pmb{k}$轴转动的角度
也就是说，每次欧拉角会产生三个坐标轴，我们分别记为$x,y,z$, $x^{’},y^{’},z^{’}$和$x^{’’},y^{’’},z^{’’}$。
而欧拉角的整个绕轴顺序是 $z,x^{’},z^{’’}$。
这也是能够造成“万向节死锁”的一种角度表示方法。其矩阵为：
$$ \begin{aligned} & D = \begin{bmatrix} \cos \phi & \sin \phi & 0 \\ -\sin \phi & \cos \phi & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} \\ & C = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos \theta & \sin \theta \\ 0 & -\sin \theta & \cos \theta \\ \end{bmatrix} \\ & B = \begin{bmatrix} \cos \psi & \sin \psi & 0 \\ -\sin \psi & \cos \psi & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} \end{aligned} $$泰特–布莱恩角（Tait–Bryan angles）#
即一般计算机图形学中的角。绕着全局的x, y, z轴依次旋转（旋转顺序可变，比如y-x-z, z-x-y等）。也是我们最熟悉的旋转表示法。在此就不赘述了。
四元数表示法#
四元数也可以表示角，同时也有矩阵形式。老生常谈的事情，这里就不赘述了。
轴-角 表示法#
用一根旋转轴$\pmb{n}$和绕此轴旋转的角度$\theta$来表示。除了上面的四元数表示，还有罗德里格斯公式（ Rodrigues’ rotation formula）。其公式很易推导，这里直接给结果：
对于任意点$\pmb{r}$：
$$ \pmb{r^{'}} = \pmb{r} \cos \phi + (\pmb{n} \times \pmb{r}) \sin \phi + \pmb{n}(\pmb{n} \cdot \pmb{r})(1 -\cos \phi) $$定点转动的欧拉定理#
当物体绕着某点转动的时候，存在如下欧拉定理：
任意次的绕点转动，都可以表示为一次绕轴转动。此轴过此转动点
此轴就是旋转矩阵$A$的特征向量。
这里证明很好证：
首先，$A$的行列式只能为$\pm1$。因为旋转矩阵一定是正交矩阵，所以有$A^T = A^{-1}$，那么有$AA^T = I$ ，等价于$det(AA^T) = det(A)det(A^T) = 1$。而$\det(A) = \det(A^T)$，所以有$det(A) = \pm1$。
进一步可以断言，$det(A) = 1$。因为当不发生转动的时候，$A = I$为单位阵，其$det(A) = 1$。而当发生一点点转动时，由于转动的连续性（连续函数，不存在间断点），这一点点转动不可能使得$det(A)$突变到$-1$。所以$det(A)$始终为1（这里应该是有严格数学证明的，不过我没找到）。
接下来证明$A$有一个特征值1。根据特征值公式得求$det(A - E) = 0$。这很好证明：
$$ \begin{aligned} & AA^T = I \\ \Rightarrow & AA^T - A = I - A \\ \Rightarrow & A(A^T - I) = I - A \\ \Rightarrow & det(A)det(A^T - I) = det(I - A) \\ \Rightarrow & det(A^T - I) = det(I - A) \\ & \text{由于A是3x3矩阵} \\ \Rightarrow & det(A^T - I) = (-1)^3det(A - I) \\ \Rightarrow & det(A - I) = (-1)^3det(A - I) \\ \Rightarrow & det(A - I) = 0 \\ \end{aligned} $$那么此特征值对应一个特征向量$\pmb{n}$。而在此特征向量上的所有点，经过$A$变换后都不变。所以$\pmb{n}$就是转轴。
无限小转动分析#
我们无法直接将牛顿定律内容搬入刚体动力学，因为：
- 由于转动的表示$\pmb{r^{’}} = A \pmb{r}$。我们知道，当做多次转动时，这些转动的顺序是不能互换的。这样我们就无法定义角速度（因为角速度是矢量，满足交换性）。所以我们通过引入无限小转动来让极限情况下，转动操作满足交换性。
- 牛顿力学只在惯性系下成立，而惯性系不含旋转。
无限小转动分析就是在处理这些问题。
定义一个无限小转动矩阵$\epsilon$为$t \rightarrow t+dt$时间内的转动，即$\epsilon \rightarrow 0$。
那么一次无限小转动的矩阵就可以记为：
$$ A = I + \epsilon $$使用此转动表示法，就可以让转动满足交换律。假设有两个转动 $A_1 = I + \epsilon_1$和$A_2 = I + \epsilon_2$，那么有：
$$ \begin{aligned} & A_1A_2 = (I + \epsilon_1)(I + \epsilon_2) = I + \epsilon_1 + \epsilon_2 + \epsilon_1\epsilon_2 \end{aligned} $$这里$\epsilon_1\epsilon_2$是高阶无穷小，可以忽略。最后结果就是$A_1A_2 = A_2A_1$。
而$A$的逆$A^{-1}$其实是$I - \epsilon$。这一点可以通过算$AA^{-1}$得到。
而$A$的转置$A^T$其实是$I - \epsilon^T$，而旋转操作一定表示一个正交矩阵（因为旋转后三轴互相垂直），而又由正交矩阵性质可知$A^T = A^{-1}$，那么就有：
$$ I + \epsilon = I - \epsilon^T \Rightarrow \epsilon + \epsilon^T = 0 $$那么显然，$\epsilon$的对角线元素必须是0，而上三角区域和下三角区域的各个元素必是相反数：
$$ \epsilon = \begin{bmatrix} 0 & -d\Omega_3 & d\Omega_2 \\ d\Omega_3 & 0 & -d\Omega_1 \\ -d\Omega_2 & d\Omega_1 & 0 \\ \end{bmatrix} $$这里的$d\Omega_i$只是一种符号，表示极小的值。
于是，旋转矩阵就可以表示为：
$$ A = I + \epsilon = \begin{bmatrix} 1 & -d\Omega_3 & d\Omega_2 \\ d\Omega_3 & 1 & -d\Omega_1 \\ -d\Omega_2 & d\Omega_1 & 1 \\ \end{bmatrix} $$那么，对于任意点$\pmb{r}$的无限小转动之差就是：
$$ \pmb{r^{'}} - \pmb{r} = d\pmb{r} = \epsilon \pmb{r} $$可以将这个矩阵乘开得到：
$$ \begin{aligned} & d\pmb{r}_1 = -\pmb{r}_1 d\Omega_3 + \pmb{r}_3 d\Omega_2 \\ & d\pmb{r}_2 = \pmb{r}_3 d\Omega_1 - \pmb{r}_1 d\Omega_3 \\ & d\pmb{r}_3 = -\pmb{r}_1 d\Omega_2 + \pmb{r}_2 d\Omega_1 \\ \end{aligned} $$你可以发现这变成了一个叉乘：
$$ \begin{aligned} d\pmb{r} =d\pmb{\Omega} \times \pmb{r} & & (1) \end{aligned} $$其中：
$$ d\pmb{\Omega} = \begin{bmatrix} d\Omega_1 & d\Omega_2 & d\Omega_3 \end{bmatrix} ^T $$而根据罗德里格斯公式：
$$ \pmb{r^{'}} = \pmb{r} \cos \phi + (\pmb{n} \times \pmb{r}) \sin \phi + \pmb{n}(\pmb{n} \cdot \pmb{r})(1 -\cos \phi) $$当$\theta \rightarrow 0$时，通过对$\sin\theta, \cos\theta$进行泰勒展开并略去二阶项（忽略高阶无穷小）后，可以得到：
$$ \begin{aligned} & \sin\theta = \theta \\ & \cos\theta = 1 \\ \end{aligned} $$那么原公式变为：
$$ \begin{aligned} & \pmb{r^{'}} = \pmb{r} + (\pmb{n} \times \pmb{r}) d\phi \\ & d\pmb{r} = \pmb{r^{'}} - \pmb{r} = (\pmb{n} \times \pmb{r}) d\phi \\ \end{aligned} $$这样，我们就得到了，在无限小转动下质点的位移变化量，此表示只和角度$d\phi$有关。
和$(1)$式做比较，还可以发现：
$$ \begin{cases} & d\pmb{r} = d\pmb{\Omega} \times \pmb{r} \\ & d\pmb{r} = \pmb{n}d\phi \times \pmb{r} \\ \end{cases} $$那么有
$$ d\pmb{\Omega} = \pmb{n}d\phi $$对两边除以$dt$，就可以得到角速度定义：
$$ \pmb{\omega} = \pmb{n}\frac{d\phi}{dt} $$这也就解释了为什么角速度是矢量。
而同时也可以得到角速度和速度的转换公式：
$$ \frac{d\pmb{r}}{dt} = \pmb{\omega} \times \pmb{r} $$额外扩展：角速度其实是轴矢量，在宇称变换下不改变方向4。
刚体运动的分解#
对于刚体上任意的向量$\pmb{G}$（可以是刚体上一点的位矢，或刚体的角动量等）（假定刚体不移动），其总能分解成两个矢量和：
$$ d\pmb{G}_{space} = d\pmb{G}_{body} + d\pmb{G}_{diff} $$其中$d\pmb{G}{space}$是在全局坐标下看此向量，$d\pmb{G}{body}$是在物体自己的局部坐标系下看此向量（一般以质心为原点，坐标轴随刚体一起旋转），而$d\pmb{G}_{diff}$则是这两个$d\pmb{G}$的差。即在两个坐标系下观察出结果的差：
比如我要观察$B$点，上图中$\vec{FB}$就是$\pmb{G}{space}$，$\vec{EB}$就是$\pmb{G}{body}$。而$\vec{FE}$就是$\pmb{G}_{diff}$。
那么显然，$d\pmb{G}{diff}$只是刚体的转动量$d\pmb{G}{rot}$，且$d\pmb{G}_{body} = 0$，于是式子可以写成：
$$ d\pmb{G}_{space} = d\pmb{G}_{body} + d\pmb{G}_{rot} = 0 + d\pmb{\Omega} \times \pmb{G} $$两边除$dt$就有：
$$ \frac{d\pmb{G}_{space}}{dt} = \frac{d\pmb{\Omega}}{dt} \times \pmb{G} = \pmb{\omega} \times \pmb{G} $$或者更一般地，可以定义算子：
$$ \begin{aligned} (\frac{d \cdot}{dt})_{space} = (\frac{d \cdot}{dt})_{body} + \pmb{\omega} \times \cdot & & (2) \end{aligned} $$正如图上所画，当选取质心为局部坐标系的原点时，刚体绕任意点的旋转可以被分解为：
- 刚体质心的转动，可以应用质点动力学/运动学
- 刚体在自己局部坐标系下的转动
这也是物理引擎中处理刚体问题的方法。
欧拉角表示无限小转动#
经过上面步骤，我们可以用欧拉角来表示无限小转动：
$$ \pmb{\omega} = \dot{\phi} \pmb{z} + \dot{\theta} \pmb{x^{'}} + \dot{\psi} \pmb{z^{''}} $$即将绕不同轴的小转动组合起来。
但这里有个问题，$\pmb{\omega}$中包含三个不相同的坐标系中的轴。这不利于计算。我们可以将其统一到最后的坐标系$x^{’’},y^{’’},z^{’’}$中。
显然，我们通过旋转方程可以得到各个轴之间的变换矩阵，所以可以利用这些矩阵将其他轴都变换成目标轴。矩阵很多也很简单，这里就不再赘述了（可见《理论力学》3.3.3节1）。总之最后的结果为：
$$ \begin{aligned} & \omega_{x^{''}} = \dot{\phi}\sin{\theta}\sin{\psi} + \dot{\theta}\cos{\psi} \\ & \omega_{y^{''}} = \dot{\phi}\sin{\theta}\cos{\psi} - \dot{\theta}\sin{\psi} \\ & \omega_{z^{''}} = \dot{\phi}\cos{\theta} + \dot{\psi} \\ \end{aligned} $$刚体动力学方程#
首先有：
刚体的角速度和刚体的局部坐标系选取无关
即，你可以在刚体上随便选原点和坐标系朝向，只要这个坐标系是"绑死在刚体上、跟着它一起转"的，那么你算出来的角速度 $\pmb{\omega}$ 永远是一样的。
惯性张量(Inertia Tensor)#
当刚体绕着固定点旋转时，刚体上某点的角动量是：
$$ \pmb{L} = m_i(\pmb{r}_i \times \pmb{v}_i) $$其中$\pmb{r}_i$是刚体上此质点相对于此旋转点的位矢。
而由上一节可知，$\pmb{v_i} = \pmb{\omega} \times \pmb{r}_i$，并且利用向量三重积公式$a \times (b \times c) = b(a \cdot c) - c(a \cdot b)$，我们就有：
$$ \begin{aligned} \pmb{L} & = m_i[\pmb{r}_i \times (\pmb{\omega} \times \pmb{r}_i)] \\ & = m_i[\pmb{\omega}\pmb{r}_i^2 - \pmb{r}_i(\pmb{r}_i \cdot \pmb{\omega})] \\ \end{aligned} $$然后把所有点乘写开，就可以得到$\pmb{L}$的各个分量表示：
$$ \begin{aligned} L_x & = \omega_x m_i(r_i^2 - x_i^2) - \omega_y m_i x_i y_i - \omega_z m_i x_i z_i \\ & = I_{xx} \omega_x + I_{xy} \omega_y + I_{xz} \omega_z \end{aligned} $$其他两个如法炮制：
$$ \begin{aligned} L_y & = I_{yx} \omega_x + I_{yy} \omega_y + I_{yz} \omega_z \\ L_z & = I_{zx} \omega_x + I_{zy} \omega_y + I_{yz} \omega_z \\ \end{aligned} $$其中：
$$ \begin{aligned} I_{xx} & = m_i({r_i^2 - x_i}) \\ I_{xy} & = -m_i x_i y_i \\ \end{aligned} $$可以看到$I_{ij}$组成了一个3x3矩阵。这就是惯性张量(Inertia Tensor)矩阵。
刚才我们算的是刚体上一个点的惯性张量。而显然整个刚体的惯性张量就是对所有点的惯性张量求和：
$$ \begin{aligned} & I_{xx} =\int_V \rho(\pmb{r})(r^2 - x^2)dV \\ & I_{ij} =\int_V \rho(\pmb{r})(r^2 \delta_{ij} - x_jx_k)dV\\ \end{aligned} $$其中，$\delta_{ij}$是一种符号：
$$ \delta_{ij} = \begin{cases} & 1, & i = j \\ & -1, & i \ne j \\ \end{cases} $$这个符号只是经常用来简化公式，没什么特别的含义。
完整的惯性张量为：
$$ \begin{bmatrix} \int_V \rho(\pmb{r})(r^2 - x^2)dV & -\int_V \rho(\pmb{r})xydV & -\int_V \rho(\pmb{r})xzdV \\ -\int_V \rho(\pmb{r}) yx dV & \int_V \rho(\pmb{r})(r^2 - y^2)dV & -\int_V \rho(\pmb{r})yz dV \\ -\int_V \rho(\pmb{r})zx dV & -\int_V \rho(\pmb{r})zy dV & \int_V \rho(\pmb{r})(r^2 - z^2)dV \\ \end{bmatrix} $$或者写成对质量积分的形式（这个形式更常见，推导见《理论力学》3.5.2节1）：
$$ \begin{bmatrix} \int (y^2 + z^2)dm & -\int xydm & -\int xzdm \\ -\int yx dm & \int (x^2 + z^2)dm & -\int yz dm \\ -\int zx dm & -\int zy dm & \int (x^2 + y^2)dm \\ \end{bmatrix} $$可以看到他是个对称矩阵。
所以最后我们就得到刚体的角动量为：
$$ \pmb{L} = \pmb{I}\pmb{\omega} $$或者用分量形式表示（只是将矩阵乘开）：
$$ \begin{aligned} \pmb{L} &= (I_{11}+I_{12}+I_{13})\omega_1 + (I_{21}+I_{22}+I_{23})\omega_2 (I_{31}+I_{32}+I_{33})\omega_3 \end{aligned} $$计算$\pmb{I}$的公式了解一下即可。在物理引擎中刚体是很规整的，其惯性张量基本上都有固定公式。
注意：惯性张量是“张量”，这是一种数学概念。其一个重要定律是张量是坐标无关的。也就是说一个刚体的惯性张量的值是固定的。但是他的分量是随着坐标系改变的（类比向量和向量在不同坐标系下的分量）。
对张量的研究是数学中的一个独立分支，在此不再赘述。想要直观理解张量可以看这个，这个和这个视频。
转动惯量(Moment of inertia)#
通过变换动能公式可以得到旋转时的动能：
$$ \begin{aligned} T & = \frac{1}{2}m_i\pmb{v}_i^2 \\ & = \frac{1}{2}m_i \pmb{v}_i\cdot(\pmb{\omega}\times \pmb{r}_i) \\ & \text{利用混合积轮换公式} a\cdot(b\times c) = b \cdot(c\times a) = c\cdot (a\times b) \\ & = \frac{1}{2}\pmb{\omega} \cdot m_i(\pmb{r}_i \times \pmb{b}_i) \\ & = \frac{1}{2}\pmb{\omega} \cdot \pmb{L} \\ & = \frac{1}{2}\pmb{\omega} \cdot \pmb{I} \cdot \pmb{\omega} \\ \end{aligned} $$令$\pmb{n}$为$\pmb{\omega}$的方向($\pmb{\omega} = \omega \pmb{n}$)，最后的公式就变为：
$$ T = \frac{1}{2}\omega \pmb{n}\cdot \pmb{I} \cdot \pmb{n} \omega = \frac{1}{2} \omega^2 I $$这里$I$是标量，为：
$$ \begin{aligned} I & = \pmb{n} \cdot \pmb{I} \cdot \pmb{n} \\ & = m_i[r_i^2 - (\pmb{r}_i \cdot \pmb{n})^2] \end{aligned} $$就叫做绕这个轴的转动惯量(moment of inertia)。
其他物理书中对转动惯量的定义都是“物体上质点到轴的距离为$r$，转动惯量为$\int r^2 dm$"。这里其实也一样。因为$r_i^2 - (\pmb{r}_i \cdot \pmb{n})^2$就是$r^2$。
注意转动惯量和惯性张量不同，他依赖于轴的位置和朝向。
转动惯量有一个“平行轴定理”：
刚体绕任意一根轴的转动惯量，等于它绕过质心且与之平行的那根轴的转动惯量，加上"把总质量集中到质心"后对这根轴的转动惯量：
$$ I = I_{cm} + M d^2 $$其中 $I_{cm}$是以过质心的那根轴的转动惯量。$M$是刚体的总质量，$d$是两根轴之间的距离。
主轴定理#
因为转动惯量$\pmb{I}$是一个实对称矩阵，所以显然有如下性质：
- 其可以被正交相似对角化成 $TIT^{-1}$，其中$I = diag(I_1,I_2, I_3)$
- $I_1,I_2,I_3$是$\pmb{I}$的特征值（也叫做主转动惯量(principal moments)），而对于的特征向量$\pmb{e_1},\pmb{e_2},\pmb{e_3}$就是$\pmb{I}$的惯性主轴(principal axes of inertia)（也就是说主轴有三个，且两两互相垂直）
惯性主轴有如下性质：
- 一般刚体转动时角动量和角速度不一定同向，但绕主轴转动时一定同向。这意味着：
- 在无外力矩时，绕主轴的自由旋转是稳定的（角动量守恒，且刚体方向不晃动）
- 如果绕非主轴转，刚体会摇晃，要维持固定转轴的话得不停地施加力矩
- 在刚体的局部坐标下(Body Frame)，主轴总是常向量，不随时间改变
如果将旋转中心视为刚体重心的话，主轴就一定过质心，此时称为质心主轴。一般说主轴都是说质心主轴。此时质心和主轴构成刚体的一个局部坐标系。
这个坐标系带来一个非常非常好的性质：
- 在此坐标系下惯性张量是对角阵，且其对角线元素总是常量。这意味着对于任意的物体，我们可以预计算其惯性张量，并且只需要存三个分量
在物理引擎中，大部分的基本几何体，只要选取其几何中心为质心，其主轴就可以自然地被选取：
- 球体：选取任意正交轴
- 盒体：选取平行于棱的三轴
- 胶囊体/圆柱体：选取中轴线上的轴，以及垂直于中轴线的平面上的任意两轴
而对于凸多边形，则需要以其几何中心为原点先算一次惯性张量，然后做正交对角化得到真正的主轴和惯性张量。
这也意味着，如果在引擎中改变了物体的质心，必须重算惯性张量和主轴。一般会使用惯性张量形式的平行轴定理5，从旧$I$得到新的，然后再算惯性主轴。
沙勒定理(Chasles Theorem)#
比起 “刚体运动的分解” 中说的分解方法，沙勒定理则更严格：
刚体在空间中的任意一个有限位移，都可以分解为：绕某一根轴的转动，加上沿这一根轴的平移
注意：平移是沿着轴的朝向平移。
这种运动又称为螺旋运动，就和螺丝转进螺母的运动一样。
欧拉动力学方程#
接下来进入力矩和转动的关系。也十分简单：
根据单质点动力学我们知道：
$$ \pmb{N} = \frac{d\pmb{L}}{dt} $$那么为了方便，选取质心主轴作为坐标系（即body frame），根据无限小转动分析中的公式$(2)$，有：
$$ \begin{aligned} \pmb{N} & = \frac{d\pmb{L}}{dt} = (\frac{d\pmb{L}}{dt})_{body} + \pmb{\omega} \times \pmb{L} \\ & = (\frac{d\pmb{I}\pmb{\omega}}{dt})_{body} + \pmb{\omega} \times \pmb{L} \\ & \text{由于是质心主轴坐标系，I不改变，所以} \frac{d\pmb{I}\pmb{\omega}}{dt} = \pmb{I}\frac{d\pmb{\omega}}{dt} \\ & = \pmb{I}\frac{d\pmb{\omega}}{dt} + \pmb{\omega}\times \pmb{L} \end{aligned} $$这就是欧拉力学方程。或者使用展开后的形式：
$$ \begin{aligned} N_1 &= I_1 \dot{\omega_1} - \omega_2\omega_3(I_2 - I_3) \\ N_2 &= I_2 \dot{\omega_2} - \omega_3\omega_1(I_3 - I_1) \\ N_3 &= I_3 \dot{\omega_3} - \omega_1\omega_2(I_1 - I_2) \\ \end{aligned} $$注意：如果不选取质心主轴作为body frame，那么公式中的导数数量会变得很多，这会导致微分方程求解十分困难。所以一般就只用此body frame的形式。
刚体能量#
前提：选取质心主轴作为body frame。
动能（克尼希定理）：
$$ T = \frac{1}{2}M\pmb{v}_{cm}^2 + \frac{1}{2}\pmb{\omega}\pmb{I}_{cm}\pmb{\omega} $$势能则需要不同情况不同分析。
对于在均匀重力场中的重力势能来说，就是重力作用在质心的势能$V = M\pmb{g}$。
刚体内部势能为常数。因为内部各点相对距离固定不变。在拉格朗日量里一般直接设为0或并入最后的常数。
如果是在其他势能，一般的做法就只是对刚体上所有点的能量求和。
刚体静力学#
根据网课1补充的一些静力学知识。
刚体平衡#
刚体平衡条件：
- 刚体初始时处于平衡状态
- 维持平衡条件
- 主矢为0：$\sum \pmb{F}_i = 0$
- 主力矩为0：$\sum{\pmb{N}_i} = 0$
其中，主矢为所有矢量相加，而主力矩就是所有力矩相加。
虽然力矩的值随着选取的中心点$O$改变，但在平衡问题中可以任意选取中心点：
在刚体平衡问题中，简化中心可以任意选取
这里$O$点的力矩为$\pmb{N}_O = \sum \pmb{r}_i \times \pmb{F}_i$
而$O^{’}$点则是$N_{O^{’}} = \sum (\pmb{r}_i + \pmb{r}_0) \times \pmb{F}_i = \pmb{N} + \pmb{r}_0 \times \sum \pmb{F}_i$
而平衡的条件就是$\sum \pmb{F}_i = 0$，所以可得$\pmb{N}_O = \pmb{N}_O^{’}$。所以任意点的主力矩都是一样的。
力系的简化#
力偶：
大小相等，方向相反，但不共线的一对平行力
力偶只会对刚体施加转动效果。
力的可传递性原理：
沿着力的作用线移动力，力的作用效果不变
力的平移定理：
力从$O$点平移到$O^{’}$点时，要额外加上$O$到$O^{’}$点的力偶矩
力偶矩为 $(\pmb{r}_1 - \pmb{r_2}) \times \pmb{F}_1$。