计算机动画的本质上是解微分方程，因此本文主要介绍微分方程常用的一些解法。

particle system

首先从最简单的内容开始，也就是particle system。

上面是particle system的一个简单的案例，也就是一个弹簧+一个block，block在弹簧的作用下左右运动，其中弹簧不受力的情况下block所在的位置为原点，这里用q表示位置，那么上述状态可以用q=0来表示，定义向左为负，向右为正，相应的速度与加速度也可以基于下面的公式得到：

动能kinetic energy与势能potential energy则由如下的公式进行表示：

其中m与k分别表示block的质量与弹簧的弹性系数，而弹簧的弹力可通过下述公式计算得到：

根据前面的公式可以知道，弹力是conservative force的一种。

所谓的conservative force指的是将一个particle从一个点移动到另一个点做的功是跟移动的路径无关的力，在这个定义下，如果将一个particle沿着任意路径移动最后回到原点，那么力的整体做功（微分下可以用表示）为0，在conservative force的作用下，不论沿着哪条路径，只要保证首尾点是相同的，那么都会导致势能的变化幅度是相同的，常见的力中，重力跟弹力是conservative force，而摩擦力则不是conservative force，从这个角度来理解，conservative force是跟位置相关的力，而这种力作用的效果就是将位置的变化转换为势能的变化。

一个力场（force field）如果满足下面三个公式中的任意一个，那么就可以说明这个力场是conservative的：

1. 力的curl是0向量
   
   在2D空间中，上述公式可以简化为：
   

2. 闭环移动时，做功为0
   

3. 力可以通过对势能的梯度取反来表示：
   

前面两个公式计算比较复杂，通常我们会使用第三个公式进行conservative force的判断，根据前面弹簧粒子的公式，我们可以判断弹力就是conservative force。

再说回弹簧粒子系统，根据能量守恒定律，我们有：

其中C是常量，对这个公式的左右两边对时间进行求导，就可以得到：

上述这个公式有两个解：

其中第一个解的含义是速度不变，也就是物件静止，这种情况我们是不关心的，我们主要关心第二个解，这个解我们换个写法：

用中文的话翻译一下，就是物体运动的加速度大小与物体的受力成正比，跟物体的质量成反比，这就是牛顿第二定律（作为回顾：牛顿第一定律也称为惯性定律，可以表述为任何物体都要保持匀速直线运动或静止状态，直到外力迫使它改变运动状态为止；牛顿第三定律可以表述为相互作用的两个物体之间的作用力和反作用力总是大小相等，方向相反，作用在同一条直线上）。

下面再来看一个复杂一点的案例，block从一个变成两个：

动能公式就变成了：

用矩阵来表示就可以写成：

这里需要注意的是，我们前面说过，单个block的情况下，我们是以弹簧原始长度时block所在的位置为原点的，而实际上我们这里选择的坐标系是比较随性的，而不是有固定的一套标准，而使用不同的坐标系，得到的公式是不一样的，虽然他们输出的结果是一样的（能量不随坐标系的变化），而这给了我们一个启发，在实际使用中要谨慎选择坐标系，不同的坐标系计算公式不一样，导致计算效率会受影响，而这对于结果的输出并没有任何收益。

数值积分

重写一下前面的牛顿第二定律的公式：

这是一个以q为变量的二阶微分方程，在物理动画模拟中，我们最终需要的是计算各个时刻的位置，也就是解出每个time step下的q，而求解微分方程的方式就是数值积分。

为了解前面这个二阶微分方程，我们这里尝试将之改写成一个一阶的常微分方程（目的是降低求解的复杂度？），令

那么有

上面为了进行一般化处理，用替代了，那么怎么来理解这个y以及y的微分呢？

如下图所示，我们以q作为横坐标，作为纵坐标绘制一个坐标系（这个坐标系这里称之为phase space，当然，这里只是一个近似，真实的phase space的纵坐标是动量，不过不影响），当我们取点来看，其微分结果为，也就是一条水平方向（向右）的向量，如下面所示：

同样，如果我们取点来看，其微分结果为，也就是一条垂直方向（向下）的向量，如下面所示：

实际上，对单block的弹簧系统而言，对于每个y坐标，在phase space中得到的y的微分，称为y的速度向量，都可以得到这个向量是沿着顺时针方向的，如果我们做特殊处理，取k = m，那么将所有点连起来，我们就能得到一个正圆，如下图所示：

不同半径的正圆对应于弹簧和谐运动的振幅，通过这个正圆也可以看到，在q为0的时候，速度的绝对值最大，反之速度为0的时候，则是位置的绝对值最大的时候；另外，对于任意一个y，其速度都是这个正圆上的切线，也就是说下一刻的位置同样是在正圆上，如果没有摩擦力的存在，那么就意味着这个block将在这个圆上永无止尽的运动下去。

在进行数值积分方法介绍之前，我们先来看下不同的积分方法的评价标准，总的来说，对于一个数值积分方法好不好，我们通常有如下三个评判标准：

1. stability，稳定性，整个系统会随着时间的变迁变得稳定，不论是匀速运动，还是随着能量的衰减逐渐趋于静止，或者在某个有限的能量的范围内来回振荡
2. accuracy，精确度，指的是我们通过数值积分模拟出来的结果跟真实的精确解之间的误差
3. convergence，收敛度，指的是我们在进行数值积分模拟计算的时候，模拟的time step在不断收缩的情况下，其模拟的精度是否不断逼近真实解。

将上面的公式做一下近似，这里用explicit euler（显式欧拉）公式进行模拟，即中y取k：

从而得到

按照上面这个公式，由于是在y点的切线，也就是说计算得到的下一个坐标点将落在当前坐标点所规范的圆形外面，那么经过这个公式计算出来的运动轨迹将不能保证是圆形的，而是半径不断扩大的螺旋曲线，也就是说，这个方法是不稳定的，即系统的能量越来越大（半径对应着振幅）：

在显式欧拉方法下，我们看下整体方案的误差，根据泰勒展式：

可以看到，显式欧拉只使用了之前的部分，也就是说是属于一阶微分方程，其误差自然就可以算出来等于了。

如果前面不用显式欧拉，而是使用implicit euler（隐式欧拉），即中y取k+1：

在这种情况下，右边是一个方程，在弹簧系统中，这个方程是可以显式表达的，但是在其他的系统尤其是复杂系统中，这个方程是没有办法表达出来的，也就是说这个公式通过解析方法是没有办法计算得到的。

根据上面的公式，仿造前面显式欧拉的逻辑，我们可以看出来，y(k)是比y(k+1)半径要大的，也就是说，在这个系统中，半径在不断减小，能量在不断衰退，最终衰减到0，那么就处于静止，也就是说这种方式得到的是一个稳定的解（当然，这里需要说清楚的是，不是所有的系统使用隐式欧拉都是稳定的，只有那些能量守恒的系统才是稳定的），不过呢由于没有解析解，因此求解会很麻烦。

下面要介绍的是一个叫做symplectic euler（辛欧拉，或者半显式欧拉，半隐式欧拉方案）的数值积分方法，根据前面的公式，我们有：


将之转换为微分方程：

又回到前面的问题，右边的向量中的q是应该取k还是k+1时刻的值，前面我们说过，取k的值是显式欧拉的方案，好处是数值是已知的，这个方程可以解出来的，缺点则是系统不稳定，而取k+1的时候则是隐式欧拉的方案，好处是方案是稳定的，而缺点则是方程求解困难，基于这一点，我们这里考虑使用一个混合的方案，即上述方程中第二个分量使用显式欧拉方案，第一个分量则使用隐式欧拉方案，如下所示：

根据上面的公式，我们可以先根据第二个分量的公式，计算出，之后这个数值可以代入到第一个公式，求得，这里借用了显式欧拉的好处，只是需要确定的是，是否同时带入了显式欧拉的缺点——不稳定。

在弹簧系统中，上述公式可以表示为：

从而求得：

根据这公式通过不断迭代可以得到：

好像也不能说明啥，写成矩阵形式：

这个等式右边的2x2的矩阵，我们可以很容易计算出其行列式=1，而根据线性代数的理论，我们可以知道，矩阵的行列式恰好等于变换前后两个形状的面积之比，也就是说，某个区域的数据集合经过这个矩阵变化后得到的新的区域，两个区域的面积之比就恰好等于矩阵的行列式，而这里行列式为1就说明经过这个矩阵的变换，数据是基本维持稳定的，不会膨胀，也不会坍缩，而是能量守恒的（自然也是稳定的系统）。

而虽然上述结论是在弹簧系统中推导出来的，但实际上这个结论的应用面是可以推广到一般的（当然，不是所有的系统都可以使用），对于系统内仅包含动能和势能的系统而言，在经典力学范畴中，这个结论是成立的，即通过辛欧拉算法得到的解都是稳定的（甚至能量守恒的）。