html5和ria网站设计/常州百度推广公司

1.前言

        上一篇介绍了动画与模拟的很多方法、模拟各种运动、基本知识。定义一个物体的速度和加速度，算出物体在任何时刻的位置，但是没有介绍具体怎么实现。这篇文章就是从上一篇的概念出发，介绍怎么把一个物体或多个物体运动的位置、不同时间出现在哪个位置怎么解出来。

2.单个粒子的模拟

2.1 定义

        我们规定了一个粒子任何一个时刻的速度，我们又知道这个粒子一开始的位置，那么我们如何求解某个时间之后它出现在哪里呢？这个很简单，高中的物理公式。假设一个物体做匀速直线运动 v，在 t0 时刻出现在某一个位置 p 上，那么在 t1 时刻出现的位置是：p + v*(t1-t0)。

        如果这个问题没有这么简单呢。如下图所示，要表示某一个粒子在一个速度场当中会做怎么样的运动呢？ 黑色箭头表示速度场，黑色曲线表示粒子运动轨迹。

        速度场与光场、磁场安全一样的概念，也就是说任何一个位置，只要知道位置就能知道这个时刻的速度，这只是一种理想化的情况。那么任何一个位置 x 和任何一个时间 t 都有一个速度 v。

2.2 常微分方程

        根据上面所说的，我们希望模拟一个如下图所示的单个粒子轨迹。

        那么，我们知道粒子在任何一个 x 位置的速度，如下图所示，前面讲过 x 上加一点表示一阶导数，也就是速度，两个点就是二阶导数，就是加速度。我们定义的如下公式写法就叫做一个常微分方程(ODE)，而是一阶的常微分方程。

        一阶常微分方程表示我们知道一个量的微分是多少，然后希望推出这个量是多少。

        而常微分方程的意思是不存在对其他变量的微分或者导数，就是只涉及一个微分的方程。带入到我们这里就是，我们知道速度，想求位置。

        想要粒子的任何时刻的位置，先要定义粒子的起始位置 x0，给任意时间 t ，都能知道在哪个位置，那么就可以随着时间运动，知道运动轨迹了。那么应该怎么解呢？

2.3 欧拉方法

        为了解出例子的运动轨迹，最简单的想法就是粒子在一定时间后位置在哪。那么就可以把时间细分成很多时间小块，每一个小块假设当前时间 t 加上某一个Δt 之后位置在哪，不断这么去做就可以。每一步考虑的Δt 就是步长，从某一个时刻到这个时刻 t 加上 Δt 的下一个时刻轨迹怎么变换。相当于对时间上进行了离散化的操作。

        具体使用欧拉(Euler)方法去实现，欧拉方法最简单的思路就是在上一个 t 时刻都有一个位置、速度和加速度，然后要计算 t+Δt 时刻的位置、速度和加速度。实现如下公式所示。

        根据上面的公式，我们每个下一个步长的位置和速度都可以根据上一个步长的速度和加速度求得。可以发现这里始终用上一个时刻的量去估计下一个时刻的量，这就叫做前向欧拉方法(Forward Euler Method)，或者叫显式欧拉方法(Explicit Euler Method)。

        但是欧拉方法非常不准确，并且有严重的不稳定问题。为什么呢？

        如下图所示，在任何一个位置上都有一个速度，通过欧拉方法可以得到下一个位置，不断模拟下去。可想而知，步长取的越小越密效果越好。

        至于稳定性问题，如下图所示的螺旋速度场，粒子应该是做一个圆周运动的轨迹，但是最终结果得到的是粒子从速度场中偏离了出去。因为步长再小也是离散的，所以无论取多小的步长，最终这个粒子一定会偏离出去。

        再看如下图所示的速度场，如果正确模拟，我们最终的粒子应该是沿着水平方向向右飞去，但是我们看到的结果是上下来回的穿梭，这同样是欧拉方法的不稳定性造成的问题。

        综上两种情况，都会产生不稳定问题，无论怎么取都与实际情况不符，这与误差是两码事。在图形学中误差不是关键问题，因为可以通过增加精度降低误差，而不稳定问题是比较关键，会偏离实际情况，所以都在想办法解决。

2.4 解决不稳定性

2.4.1 中点法

        中点法(Midpoint method)与欧拉方法类似，只不过做了一步改进。如下图所示，先根据起始位置和方向，用欧拉方法模拟某个Δt 时间到达 a 位置。但是我们不用 a，我们取原点与 a 的中点 b，考虑中点 b 所在速度。然后使用起始位置与中点 b 的速度用欧拉方法重新算一遍，算出Δt 时间后到达 c 位置 。

        中点法的核心仍然用的欧拉方法，只不过用了两次，第一次为了得到中点的速度，第二次用中点速度再算一次欧拉方法。因为速度场的变化在两点之间可能特别快，而中点的速度某种意义上来说比初始点更具有代表性。整理公式如下。

        从公式中可以看出多了一个2次方，路线会是一种抛物线的轨迹，就是局部的二次模型，当然比一次模型更加准确。

2.4.2 自适应步长

        如下图所示，对于一个步长Δt，我们用欧拉方法算出 XT 的位置。然后重新从起始位置用Δt/2算出一个位置，再用这个位置的速度用欧拉方法经过Δt/2算出 XT/2 位置。如果 XT 和 XT/2 位置差的很远，就说明我们取的步长过大了，所以要缩小步长，直到最终计算的两个结果相差在我们规定的范围内停下来，这就是自适应步长法(Adaptive Step Size)。

2.4.3 隐式欧拉方法

        隐式欧拉方法(Implicit Euler Method)又被称为后向(Backward)方法。显式欧拉方法是根据上一个时间的位置、速度与加速度来更新下一个时间的位置和速度，而隐式欧拉方法用的则是下一个时刻的速度和加速度，如下图所示。

        而下一刻的速度和加速度是未知的，我们要求解的是一个非线性的常微分方程。这种情况下就不好解了，所以求解通常需要用到牛顿迭代法，这里不详细介绍，有兴趣的朋友自行了解。我们只需要知道，隐式的欧拉方法更稳定。

2.4.4 稳定性的定义

        如何定义一个方法是不是稳定的，以及有多么稳定？

        通常会定义两个概念，一个叫做局部误差(local truncation error)，也就是每一个步长产生的误差，每一步的误差累积起来叫做全局误差(total accumulated error)，用这两个概念来衡量稳定性。

        正常情况下，人们认为研究这两个数的大小是没有意义的，我们需要研究的是这两个数和Δt的关系，也就是阶数，通常用O(h)表示。用更小的Δt误差自然就减小，那么误差是如何随着更小的Δt而减小的，它们之间又有什么关系，这就是研究的问题。

        这里直接给出结论，隐式欧拉方法是1阶的，也就是说局部误差是O(h^2)，全局误差是O(h)。h 就是步长Δt，Δt越小，误差越小。对O(h)的理解，如果我们把 h 减半，我们可以期望得到的误差也减半。也就是阶数越高越好，因为如果有个平方的阶数，h减半，那么误差可以减少到4分之一。

        所以人们去研究这些数值方法来计算对应的稳定性，就是用阶数表示稳定性。

2.4.5 龙格库塔法

        在数值的计算领域，有一类非常有名的方法叫做龙格库塔法(Runge-Kutta Families)。

        龙格库塔法不是一个方法而是一类方法，它们擅长解常微分方程(ODEs)，特别适合求解非线性的常微分方程，其中有一种龙格库塔的方法叫RK4用的最广泛，4就是4阶的意思，可想而知误差很小。RK4解决的问题跟各种欧拉方法一样，过程如下所示。

        一开始在 t0 时刻出现 y0 初始位置上，然后从某一个时间到下一个时间，就是加一个 h (就是Δt)，与之前方法类似，就是更新位置的时候不一样，更新位置用了很多 k 系数，这些都是精心设计的。这是一门数值分析的课程，这里不再讲述。

2.4.6 Position-Based / Verlet Integration

        Position-Based或者叫Verlet Integration不是一种基于物理的方法。之所以被称为“非物理”的方法，是因为它们主要关注在给定时间步长下尽可能准确地模拟物体的运动，而不是严格遵循物理规律。

        实际应用中，这些简化和近似手段通常可以带来较好的性能和视觉效果。其核心思想是通过直接操作物体上的粒子的位置信息，从而达到约束物体运动的目的(调整位置从而满足某种限制)。当然因为不基于物理，它有时很难保证它符合能量守恒定律。

3.刚体模拟

        刚体模拟(Rigid Body Simulation)非常简单，因为它不会发生弹性形变，说明它内部所有的点运动形式都是相同的。这说明刚体模拟其实就相当于一个粒子模拟，只不过刚体模拟要考虑更多的物理量，这些量对时间的导数是多少，比如位置，朝向，速度，角速度，如下图所示。

        考虑刚体的模拟无非就是在原本的位置和速度基础上加了这么一些量。比如位置 x 对时间的求导就是速度，角度对时间的求导就是角速度，速度对时间的求导就是加速度 F/M。我们同样可以用之前说的欧拉方法以及一系列方法去求解。

4.流体模拟

4.1 模拟水原理

        流体模拟(Fluid Simulation)就是我们之前所说的非物理的Position-Based方法的例子。通过模拟这些形成整个水体积的小球的位置，来模拟整个水的运动。如下图所示。

       首先我们认为水体是由许多不可压缩的刚体小球组成的。只要精确的模拟这些小球的运动，接下来浪花，水面就是渲染的事情了，而模拟和渲染是分开的过程。

        其次，假设水在任何一个地方都是不可压缩的，2L水不可能变成1L水，也就说任何位置密度都相同。所以，只要密度在某个地方发生变化(假设我们一开始平静的水面是各处密度相等的)，它应该通过改变小球的位置来“纠正”。

        为了做这个纠正，我们需要知道每个小球的密度对所有的小球的位置这一个梯度是多少，，或者导数是多少。

        梯度可以这么理解，任何一个小球考虑周围有一圈小球都会影响它的密度，如果移动任何一个小球都让它的密度发生一定的改变。这就是这个密度对于这一个小球位置的导数，那当然可以求这个密度对所有其他小球位置变换的导数。这样就可以得出结论：任何一个点处的密度都是任何一个其他小球位置的一个函数，那当然可以求它的导数，怎么做不用关心。

        接下来，任何一个位置上让它原本不正确的密度回到正确的密度上面去，并且又知道如何调整各个不同小球的位置使得它的密度往希望的方向去变换，这就是一个梯度下降(Gradient Descent)的过程。

        通过在任何一个局部去调整各个不同点的密度变化，那自然就可以在任意一个时刻都可以保证水的密度在任何一个地方都是不变的。整个过程没有任何方程，都是非物理的。任何时候都知道小球应该如何运动，那自然可以把这个过程随着时间重复进行模拟出效果。

4.2 拉格朗日方法和欧拉方法

        在物理模拟中，模拟大规模的物质的两个基本思路：拉格朗日(Lagrangian)方法和欧拉(Eulerian)方法。

        拉格朗日方法，俗称质点法就是我们上面提到的对于一大片水我们把它当成许多小水滴，然后我们找出每个小水滴的位置，最终整个水的位置就是正确的。如下图所示，我们模拟一群小鸟，把每个小鸟的位置都找出，我们就得出了整个鸟群的样子，当然这其中也会考虑它们之间的约束关系。 

         欧拉方法，或者叫网格法(这和我们之前讲的解常微分方程不是一回事，这里是指如何去看待模拟一系列大规模物体)。我们把空间分成许多网格，我们考虑一个网格内随着时间如何变化。如下图所示，我们看一个网格，t-1 的时间里网格里是蓝色的鸟，t 时间是黑色鸟，t+1 时间是橙色的鸟。

        两种方法的视角不同，拉格朗日视角是盯着物体看，欧拉视角是盯着空间看。

4.3 物质点法

        当然，也有两种方法结合的方法，叫做物质点法(MPM，Material Point Method)，或者叫材质点方法，既考虑欧拉方法，又考虑拉格朗日方法。我们先假设物体的各个粒子都有各自的材质属性，比如密度等等，然后模拟它的运动过程，但是这些过程我们在格子中考虑，最终格子中记录的信息我们再写回粒子中。