位置约束方程
直接限制物体的空间位置关系
$$ 对于有s个自由度的系统,其位置约束方程为：C(q_1,q_2,...,q_s,t)=0 $$
我们可以直接使用这个方程对物体位置进行更新，但是其计算成本高，适合离线或低实时性场景，并不适合物理引擎中的实时低性能场景。
例如将一个物体的质心约束在一个固定的点p上：
$$ C=\left | \vec{c_{a}}-\vec{p} \right | =0 $$
表示$c_a$到点p的距离为0
速度约束方程
速度约束方程限制物体的瞬时速度关系，一般用大写字母 C (Constraint equation)来表示。约束方程的一般形式：
$$ 对于有s个自由度的系统，对于其广义坐标q_1,q_2,...,q_s,\overset{.}{q_1}\overset{.}{q_2},...,\overset{.}{q_s},t \\ 其约束方程为：C(q_1,q_2,...,q_s,\overset{.}{q_1}\overset{.}{q_2},...,\overset{.}{q_s},t )=0 $$
速度约束通常更稳定、更易构建且更便于求解，接下来介绍一下在二维空间构建一组成对由速度约束推导出的冲量公式。
在二维空间中一个物体只有3个自由度，x，y轴平移和旋转。所以我们的位置约束方程为
$$ C(x(t),y(t),\theta(t))=C(\vec{p}(t))=0 $$
位置约束方程对时间微分即可得到速度约束方程
$$ \frac{d}{dt}C(\vec{p}(t))=J\frac{d}{dt}\vec{p}(t)=J\vec{v}=0 $$
其中$J$是雅可比矩阵(Jacobian Matrix)
先介绍一下几个概念
质量矩阵
在刚体动力学中，研究某个点的动力学状态时（如力的作用点），施加一个力不仅会引起平动，还可能引起转动（如果点不是质心）。要描述这种耦合效应，单一标量物理量（如质量或转动惯量）是不够的。
所以我们定义一个质量矩阵用于综合描述点的平动惯性和转动惯性。
$$ M=\begin{bmatrix} m & & & \\ & m & & \\ & & I & \\ \end{bmatrix} $$
其中m是物体的质量，I是物体的转动惯量。
我们可以写成分块矩阵的形式
$$ M=\begin{bmatrix} M_a & \\ & I_a \end{bmatrix} $$
动量矩阵
动量矩阵就可以表示为
$$ P=MV=\begin{bmatrix} m & & & \\ & m & & \\ & & & I\end{bmatrix} \begin{bmatrix} v_x \\ v_y \\ \omega \\ \end{bmatrix} $$
或者
$$ P_a=MV=\begin{bmatrix} M_a & & \\ & &I_a \end{bmatrix} \begin{bmatrix} \vec{v_a}\\ \omega_a \\ \end{bmatrix} $$
牛顿第二定律
牛顿第二定理就可以写作
$$ F=MA=\begin{bmatrix} m & & \\ & m & \\ & & I \end{bmatrix}\begin{bmatrix} a_x \\ a_y \\ \alpha \end{bmatrix}=\begin{bmatrix} M_a & \\ & I\end{bmatrix}\begin{bmatrix} \vec{a} \\ \alpha \end{bmatrix} $$
好了我们接着继续
我们需要更新速度使得满足约束方程
使用常微分方程数值方法中学到的前向欧拉法一阶数值积分（其实就是我们高中学的速度-加速度关系的离散形式）
$$ V_t= V_0+A \Delta t $$
Delta t就是我们的物理帧，那么我们只需求得A就可以更新速度了！
使用牛顿第二定律
$$ A=M^{-1}F $$
使用在拉格朗日乘数法可得（经典拉格朗日力学中约束力的数学表达）
$$ F=J^T\vec{\lambda_f }=J^T\begin{bmatrix} \lambda_x \\ \lambda_y \\ \lambda\omega \end{bmatrix} $$
其中λ_f 为拉格朗日乘子，表示力的大小，其分量可更新不同方向上的速度和角速度；J 为雅可比矩阵的转置
最后依次代入求解未知量为λ_f的一次线性方程方程组即可。
代入有
$$ J(V_0+M^{-1}J^T\vec{\lambda _f}\Delta t)=0 \\ $$
$$ J(V_0+M^{-1}J^T\vec{\lambda }_{imp})=0 \\ $$
$$ JM^{-1}J^T\vec{\lambda }_{imp}=-JV_0 $$
这正是我们需要的一次线性方程方程组
$$ AX=B\\其中A=JM^{-1}J^T \\ $$
$$ X=\vec{\lambda }_{imp} \\ $$
$$ B=-JV_0 $$
我们只需要求解该其次方程组的解即可得到更新所需的冲量(物理引擎一般都是使用冲量更新物体的动力学状态的)
