Animation (cont.)

Single Particle Simulation

我们规定了一个物体任何一个时刻的速度，又知道一开始出现在某一个位置上，我们想求解的是：在某个时间，物体会出现在哪。

First study motion of a single particle

Later, generalize to a multitude of particles

例如我们想解某一个粒子在一个速度场中怎么运动。

📖

速度场（velocity vector field），和光场、磁场等是一个概念。也就是只要我们知道了位置，我们就知道了这个时刻的速度。这是一种理想化的情况。实际情况中很难知道，一般只能知道一个位置受多少力。

右图假设有粒子沿着类似水流的方向走。任何位置 x 在任何时间都有一个速度。

例如右图左下角，线的最左端是一个物体的初始位置，线的切线方向就是物体的速度，往前走一点就会改变一点速度，随着速度场运动下来。

To start, assume motion of particle determined by a velocity vector field that is a function of position and time:

Ordinary Differential Equation (ODE)

我们希望做出一个模拟：

对任何一个位置都知道其速度（一阶导数，表示加速度，Dot Notation for Derivatives 中提到过）。

速度的定义就是：

这种写法就是一阶常微分方程（Ordinary）。一阶指的是我知道这个量的微分是多少，希望知道这个量是多少。

“Ordinary” means no “partial” derivatives, i.e. x is just a function of t

📖

在初等数学中就有各种各样的方程，比如线性方程、二次方程、高次方程、指数方程、对数方程、三角方程和方程组等等。这些方程都是要把研究的问题中的已知数和未知数之间的关系找出来，列出包含一个未知数或几个未知数的一个或者多个方程式，然后取求方程的解。但是在实际工作中，常常出现一些特点和以上方程完全不同的问题。比如：物质在一定条件下的运动常微分方程常微分方程变化，要寻求它的运动、变化的规律；某个物体在重力作用下自由下落，要寻求下落距离随时间变化的规律；火箭在发动机推动下在空间飞行，要寻求它飞行的轨道，等等，要以现有数据求得出形式上的函数解析式，而不是以已知函数来计算特定的未知数。

常微分方程（数学概念）_百度百科

常微分方程，属数学概念。学过中学数学的人对于方程是比较熟悉的；在初等数学中就有各种各样的方程，比如线性方程、二次方程、高次方程、指数方程、对数方程、三角方程和方程组等等。这些方程都是要把研究的问题中的已知数和未知数之间的关系找出来，列出包含一个未知数或几个未知数的一个或者多个方程式，然后取求方程的解。但是在实际工作中，常常出现一些特点和以上方程完全不同的问题。 [1] ...

https://baike.baidu.com/item/%E5%B8%B8%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B/818194

常微分方程指不存在对其他变量的微分（导数），也就是单变量的微分方程。

Solving for Particle Position

We can solve the ODE, subject to a given initial particle position , by using forward numerical integration

现在我们知道速度，想知道任何时刻它的位置。

我们想知道任何时刻的位置，就要给它一个起始位置。给任意时间就能告诉我们它会在哪个位置。这样我们就能把点随着时间运动的轨迹给求出来。

📖

如何求解？

Euler’s Method

我们要求一定时间后，点的位置在哪。那么我们就将时间细分成很多个时间的小块，如果当前是时间，然后我们就对每一个小块计算后在哪，不断这么做，就能解出来了。

这里每一步的步长就是，表示从某一个时刻，到这个时刻加上后位置会怎么变化，就相当于对时间进行离散化的操作。

如何做？最简单的是欧拉方法。在上一个时刻有位置、速度、加速度，要算的位置、速度、加速度。

Euler’s Method (a.k.a. Forward Euler 前向欧拉, Explicit Euler 显式欧拉)

Simple iterative method

Commonly used

Very inaccurate

Most often goes unstable

会迅速地变得不稳定

下一帧的加速度和位置都可以由上一帧的位置、速度、加速度算出来。

Errors 误差

With numerical integration, errors accumulate

Euler integration is particularly bad

回到速度场上，不考虑加速度，只考虑速度和位置。那么在任何一个位置上都有对应的速度，我们要计算下一个后会出现在哪，从而递归地计算下去。

如果我们用不同的步长，越小就会越准确从而减小误差。

Instability of the Euler Method 不稳定

考虑到稳定性，如右图，不管我们取多大的步长，会发现无论如何得到的线都不可能沿着螺旋形走，从一个螺旋形的切线变到另一个螺旋形的切线，一直走，最后一定会离开螺旋形的速度场。但现实上来说，如果任何一个时刻都有和位置垂直的速度，粒子无论如何都应该是做圆周运动。

如由下图，它的场也会更奇怪。折线左下角会有一个速度把粒子打到上面去，之后又发现上面有个速度要往下，就往下打，之后速度又向上。正确地模拟的话应该要按照曲线慢慢走到水平，这样的模拟会出现越模拟越大的情况。

📖

信号与系统中，这种叫正反馈。出现了问题后，问题会无限的放大。

📖

这种和误差是两码事，不管步长取多小，最后一定都会出现和实际结果无限远的情况。

The Euler method (explicit / forward)

Two key problems:

Inaccuracies increase as time step increases

Instability is a common, serious problem that can cause simulation to diverge

Errors and Instability

Solving by numerical integration with finite differences leads to two problems:

一切用数值方法去解螺旋方程都会出现的问题 Errors

Errors at each time step accumulate. Accuracy decreases as simulation proceeds

由于每一步都有误差，因此累积起来肯定也有误差。步长较小可以降低误差。

Accuracy may not be critical in graphics applications

有时候误差可能问题不大，例如某些图形学应用中，模拟看起来效果还行，而不用关注真实的物理。

Instability

Errors can compound, causing the simulation to diverge even when the underlying system does not

diverge：有任何的某一种模拟方法，但不管如何模拟，最后的结果和实际的结果差的越来越远，这种就是基础的问题。因此大家都试图来解决不稳定的问题。

Lack of stability is a fundamental problem in simulation, and cannot be ignored

绝地求生中车的物理模拟就是一个例子，不稳定就会出现这样的问题。

💡

我们之前用欧拉方法来解微分方程，那么有没有其他方法来解决这种不稳定性？

Combating Instability

Some Methods to Combat Instability

Midpoint method / Modified Euler 中点法

Average velocities at start and endpoint

Adaptive step size 自适应改变步长

Compare one step and two half-steps, recursively, until error is acceptable

Implicit methods 隐式法

Use the velocity at the next time step (hard)

Position-based / Verlet integration 不基于物理的方法

Constrain positions and velocities of particles after time step

Midpoint Method

我们不希望欧拉方法在模拟的过程中结果离得越来越远。中点法一开始有个位置方向，可以直接用欧拉方法来模拟某个，例如下图中，点从原来的位置（最左侧中间）到达 a 点，但我们不用这个点。我们用原点和 a 点的中点 b 点，考虑中点的速度，然后再回到原始的出发点，用中点 b 点的速度重新算一遍欧拉方法。下图中，点走的不是弧线方向，而是下面的直线方向到达 c 点。

💡

中点法比欧拉方法要准确，只不过每次要用两次欧拉方法。

Midpoint method

Compute Euler step (a)

Compute derivative at midpoint of Euler step (b)

Update position using midpoint derivative (c)

Modified Euler

Average velocity at start and end of step

Better results

式子展开之后可以发现，中点法比欧拉方法多了二次的项。中点法模拟出了抛物线一样的运动轨迹。

原本欧拉方法可以看成是局部线性的模型。中点方法相当于算出了局部的二次的模型，因此比一次的模型准确。

💡

中点法有一个更好的应用：用中点法可以得到两个点：a 点和 c 点。我们能不能将两个方法来结合起来？我们可以用中点，但不是用中点方法的中点，而是应用的中点的思想。

Adaptive Step Size

应用一个时间，比如右图上一开始在最左边的点，用欧拉方法得到的是点，这个点不准，于是我们可以将时间减半分成两个，算两遍。算第一遍算到的是原点和点中间的点，第二遍用第一遍得到的点的速度算出点。如果和分得挺远，则意味着将分成两个这样做更准确，我们应该将时间再拆分得更短一点。

这会给我们一个标准，是否将划分成两部分，取决于差得远不远。如果两个点相差不远，说明用的步长已经足够小了，就不用继续分下去了。这是自适应的思路。

Adaptive step size

Technique for choosing step size based on error estimate

Very practical technique

But may need very small steps!

Repeat until error is below threshold:

Compute an Euler step, size

Compute two Euler steps, size

Compute error

If (error > threshold) reduce step size and try again

Implicit Euler Method 隐式欧拉方法、反向法

Implicit methods

Informally called backward methods

说成反向或后向更好理解，因为其核心思想就是我们用的 derivatives （梯度、导数）用的永远都是下一个导数。

Use derivatives in the future, for the current step

欧拉方法用的是上一个时间的位置、速度、加速度来更新这个时刻的速度和位置。

对于隐式欧拉方法，用的速度是下一帧的速度，用的加速度也是下一帧的加速度。

但是用的下一时刻的速度、加速度我们还不知道，我们解出来是一个方程组，方程组解出来一个位置、速度，因此这种情况比较不好解。我们对于一个普通的运动是很好解的，但如果速度、加速度不是简单的线性方程叠加的，就不好解了。

这个时刻我们知道位置，可以认为下一个时刻加速度知道，可以解下一时刻的位置、速度，用两个等式，两个未知数肯定有解。我们可以用各种优化方法或求根公式来解，例如用牛顿法来求函数的零点。

Solve nonlinear problem for and

Use root-finding algorithm, e.g. Newton’s method

Offers much better stability

如果要解出某个值，和数值的解法肯定要慢很多，但隐式欧拉方法能提供很好的稳定性。

💡

怎么定义一个方法是稳定的，有多么稳定？

How to determine / quantize “stability”?

We use the local truncation error (every step) / total accumulated error (overall)

一般数值方法的定义会定义两个概念：局部的 truncation error 截断误差（每一步会产生多少误差）、total accumulated error 整体最后算出来的误差（局部误差加起来）。

Absolute values do not matter, but the orders w.r.t. step

研究这两个数是没有意义的，应该研究他们的阶。阶指的是他们取的的关系。前面说过用更小的误差就能减小。但是误差是如何随着的减小而减小的，它们是什么关系。
下面就是结论，隐式方法是一阶的。

Implicit Euler has order 1, which means that

Local truncation error 局部误差: and
Global truncation error 整体误差: ( is the step 步长, i.e. )

越小，最后误差肯定越小。

Understanding of

If we halve h, we can expect the error to halve as well
指：如果（或）减小一半，那么误差也会减小到一半。
不同的方法有不同的阶数，阶数越高越好。

Runge-Kutta Families 龙格-库塔方法

📖

数值分析中，Runge-Kutta 法（Runge-Kutta methods）是用于非线性常微分方程的解的重要的一类隐式或显式叠代法。 Wikipedia 龍格－庫塔法

A family of advanced methods for solving ODEs

Especially good at dealing with non-linearity

很擅长解 ODE，特别对于非线性的情况。
例如场比较曲的时候，中点法就比前向欧拉法要好，因为前向欧拉法认为它是线性的模型，而中点法认为是个平方级的模型。

It’s order-four version is the most widely used, a.k.a. RK4

📖

在各种 Runge-Kutta 法当中有一个方法十分常用，以至于经常被称为「RK4」或者就是「Runge-Kutta法」。该方法主要是在已知方程导数和初始值信息，利用计算机仿真时应用，省去求解微分方程的复杂过程。 Wikipedia 龍格－庫塔法

RK4 解决的问题和各种欧拉方法解法的都是一样的

例如有某一个一阶导数，位置为，则（速度，位置的导数）和时间有关系，是一个函数，一开始会出现在某个位置上：。

知道了初始情况后，时间的计算也是一样的，下一帧时间在当前时间上加，只是在更新位置的时候和其他方法不一样。更新时用上一帧的量加上乘以某个平均值

平均值其中的几个量其实都是速度场不同的位置以及不同的时间的值是多少。这其实就是中点法的推广，只不过中间的量是经过精心设计的，最后的阶数也是可以严格证明的。

📖

数值分析课中可以深入了解，对图形学也是很有用的。

Position-Based / Verlet Integration 不基于物理的方法

最简单的方法，只通过调整各种位置，使得最后满足某种限制。

有时候不基于物理的方法还挺好用，例如可以认为某个弹簧是弹簧，只不过在拉开位置之后立刻会恢复原状，可以将其认为是个劲度系数无限大的弹簧。这个模拟很好用，可以理解为当弹簧拉开多长之后，会立刻想办法调整弹簧的两个端点，使得弹簧能立刻回到原始的长度。也就是说中间经过了一个很复杂的过程，通过非物理的简化方式，就能直接改变其位置。

后面讲流体会提到，怎么调整粒子不同的位置。

Idea:

After modified Euler forward-step, constrain positions of particles to prevent divergent, unstable behavior

Use constrained positions to calculate velocity

Both of these ideas will dissipate energy, stabilize

Pros / cons

Fast and simple

不需要经过物理，只需要把位置调好就对了。

Not physically based, dissipates energy (error)

因为不经过物理，有时候不太能保证一些性质，例如能量守恒。
如果我们在绳子模拟器上应用了 Verlet Integration，会立刻发现绳子停下来，有很大的能量损失。

讲完了粒子，来讲刚体的模拟。

Rigid Body Simulation

刚体不会发生形变，也就是让内部所有的点都按照同一种方式去运动。也就是说刚体的运动其实就是粒子的运动，能模拟一个粒子就一定能模拟一个刚体，只不过刚体的模拟中，人们会更多的考虑其他物理量，例如下图中位置、刚体朝向等。

Simple case

Similar to simulating a particle

Just consider a bit more properties

下图就是四个量分别对时间求导之后的结果，例如位置对时间求导之后就是速度，角度对时间求导就是角速度，速度求导之后变成加速度，角速度求导之后就是角加速度。

式子写出来，自然就能给定时间，用任何的数值方法例如欧拉方法，去求出任何时间之后，刚体对应的位置，以及它们是怎么旋转的。

Fluid Simulation

A Simple Position-Based Method

通过模拟形成整个水体积的小球的位置，来模拟整个水、浪花的运动。

Key idea

Assuming water is composed of small rigid-body spheres

首先认为整个水体是由不可压缩的刚体小球组成的。只要能精确的模拟整个小球的运动，剩下的就是渲染的问题了。

Assuming the water cannot be compressed (i.e. const. density)

这里有个重要的假设：水在任何地方都是不可压缩的。
在任何时刻，水在任何位置，其密度应该是一样的。

So, as long as the density changes somewhere, it should be “corrected” via changing the positions of particles

有了这样的假设，就有了 Position-Based 基本的解决思路：给我们任何一个时刻，对于这些小球的分布、位置，都可以知道其密度。如果说有任何一个地方的密度变成和一开始平静的水的密度不一样的情况，就需要把密度给修正过来。通过移动小球的位置，使得任何一个密度不正确的地方给修正回来。
修正的思路其实就是在模拟一个个水花（小球）怎么去运动。

You need to know the gradient of the density anywhere w.r.t. each particle’s position

为了做这个修正，我们需要知道任何一个点的密度对所有的小球位置的梯度（导数）是多少。
任何一个点，如果考虑周围一圈粒子都会影响它的密度，那么如果移动任何一个粒子，那么都会让它的密度发生改变。而上面说的梯度就是密度对于这一个粒子位置的导数，那么我们当然也能求这个密度对所有其他粒子位置变化的导数。
例如上图中，如果考虑左下角粒子对于右上角飞着的粒子，飞着的粒子的导数就是 0，因为这个粒子的位置变化不太可能影响到左下角的店。如果考虑左下角粒子的密度对于周围粒子的位置关系，如果这个粒子变化一下位置，就很有影响到其密度。
任何一个点的密度都是任何一个其他的粒子的位置的函数，当然能求出其导数。怎么做不用关心。

Update? Just gradient descent!

之后就是要让不正确的密度回到正确的密度，我们也知道如何调整小球的位置使得其密度往我们要的方向去变化，这是个梯度下降（机器学习）的过程。如果要让我们的目标与结果相似，如何调整输入让结果和目标相似，这个过程就是标准的梯度下降。

通过在任何的局部去调整不同的点的密度变化，自然就可以在任何一个时刻，保证水的密度始终在任何一个地方都是不变的。

这些过程中没有解物理的方程，也没有用到什么欧拉方程，这不是个基于物理的方法。

📖

这种修正会不会让粒子停不下来？会，会出现一种情况，让粒子不断的运动，水从左到右再到左到右。Position-Based 本身不是基于物理，它应该会有能量损失，最后粒子应该会停下来。当然也可以人为的加入一些能量损失，让其在更改的过程中产生能量的衰减。实际的过程中，人们会用一些复杂的方法让其停下来。

Eulerian vs. Lagrangian

这里介绍物理模拟中，在模拟大规模的物质的过程用到两个不同的思路。

模拟水时，我们认为水是由圆形的小水滴组成的，这些小水滴逐个模拟，得到的结果就是对的，这种方法就是拉格朗日方法，俗称质点法。

如果要模拟一群小鸟的移动，就盯着某一只小鸟是如何移动的，然后把每一只都正确模拟了就没问题了，当然还要考虑它们的相互作用。

另一种方法叫欧拉方法，这和我们前面说的解常微分方程的欧拉方法不是一回事。这里说的欧拉方法是指如何去看待模拟的一系列物体。我们将空间分成不同网格，不管网格里的东西是出去了还是进来了，只考虑网格随着不同时间应该如何变化。

如下图右下角，把空间切成条状的格子，那么在时间应该显示黑色的鸟，在时间应该显示蓝色的鸟（中应该还能看到蓝色的鸟飞到这个网格，但下图没有表示出来）。我们盯着看的是某一个固定空间，在不同时刻会出现不同的鸟。

如果要对流体进行模拟，就要把空间切成不同网格，然后知道网格中每一个格子的密度应该如何去变化。我们盯着格子来看而不是盯着物体来看。

Material Point Method (MPM) 物质点方法

最近有人把两种方法（欧拉方法、拉格朗日法）结合到一块。

Hybrid, combining Eulerian and Lagrangian views

Lagrangian: consider particles carrying material properties

认为不同的粒子都具有材质属性，下图兔子中有不同的粒子，粒子本身有粘性、质量。

Eulerian: use a grid to do numerical integration

然后去模拟如何融化，这些过程在格子里面做，信息都记录在格子上。

Interaction: particles transfer properties to the grid, grid performs update, then interpolate back to particles

然后把信息写回不同的粒子上去。

📖

模拟仿真非常火热的话题，如何尽快地模拟。