Lecture 22 Animation Cont
一、Single particle simulation
First study motion of a single particle
Later, generalize to a multitude of particles
To start, assume motion of particle determined by a velocity vector field that is a function of position and time: v(x, t)
假设把粒子放在某个速度场内，那么粒子应该会随着速度场而运动，那么就要模拟这个粒子速度随时间变化的过程。
Ordinary Differential Equation (ODE)
Computing position of particle over time requires solving a firstorder ordinary differential equation:
“First-order” refers to the first derivative being taken.
“Ordinary” means no “partial” derivatives, i.e. x is just a function of t.
We can solve the ODE, subject to a given initial particle position x0, by using forward numerical integration
这里用常微分方程去求速度。
（一）、Euler’s Method (a.k.a. Forward Euler, Explicit Euler)
Simple iterative methodCommonly usedVery inaccurateMost often goes unstable最简单的一种方法就是把时间分成小段，每一个小块计算时间△t。知道某个时刻t的位置和速度，计算下一时刻t+△t的位置和速度。
Euler’s Method - Errors
With numerical integration, errors accumulate Euler integration is particularly bad
Example:
如图，用不同大小的步长△t，得到的结果会不同，△t分的越小，得到的结果就会越准。
Instability of the Euler Method
对于这样一个环型速度场，无论取多小的步长△t，得到的轨迹都不可能沿着一个环走，因为速度的方向一直在改变，每一步总会有偏差，那么这个偏差是不断累积的，就会偏的越来越多最后飞出去。
Two key problems:
Inaccuracies increase as time step Δt increasesInstability 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 proceedsAccuracy may not be critical in graphics applications
Instability
Errors can compound, causing the simulation to diverge even when the underlying system does notLack of stability is a fundamental problem in simulation, and cannot be ignored
（二）、Combating 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
1、Midpoint Method
中点法
Compute Euler step (a)Compute derivative at midpoint of Euler step (b)Update position using midpoint derivative ©
中点法：
①、由于一开始出发点有位置和速度，可以直接用欧拉方法模拟得到点a
②、取中点b，考虑中点所在的速度
③、利用中点所在的速度，回到出发点，再算一次欧拉方法，得到点c
这样算出来的值就比直接用欧拉方法算出来的效果好很多。
Average velocity at start and end of stepBetter results
将式子整理写开可以发现，中点法比欧拉方法多了一个二次的项，相当于在模拟抛物线，因此更准确
2、Adaptive Step Size
自适应方法
Adaptive step size
Technique for choosing step size based on error estimateVery practical techniqueBut may need very small steps!
Repeat until error is below threshold:
Compute xT an Euler step, size TCompute xT/2 two Euler steps, size T/2Compute error || xT – xT/2 ||If (error > threshold) reduce step size and try again
原始点在用欧拉方法计算经过△t后运动到了xT点，显然不是很准，那么这时候把时间减半变成△t/2，用△t/2通过欧拉方法算两遍，第一遍算到了原始点和xT的中点，第二次算到了xT/2这个点。然后比较xT和xT/2的位置，如果差的很远，那么就应该把时间再分的更小去计算。直到两个点的位置差不多，说明时间已经分的足够小了，这时候得到的结果是比较准确的。
3、Implicit Euler Method
隐式欧拉方法
Informally called backward methodsUse derivatives in the future, for the current step
下一时刻位置是用下一时刻的速度和加速度计算的。
Solve nonlinear problem for xt+△tand x˙\dot{x}x˙t+△tUse root-finding algorithm, e.g. Newton’s methodOffers much better stability
但是这个式子并不是很好解。假设下一时刻加速度知道，那么可以用各种优化方法来解。隐式欧拉方法可以提供很好的稳定性。
How to determine / quantize “stability”?
We use the local truncation error (every step) / total accumulated error (overall)
Absolute values do not matter, but the orders w.r.t. step
Implicit Euler has order 1, which means that
--Local truncation error: O(h²) and
--Global truncation error: O(h) (h is the step, i.e. ∆t)
Understanding of O(h)
--If we halve h, we can expect the error to halve as well
如何定义这个方法是不是稳定的以及多么稳定？
通常对数值方法会定义两个概念：局部截断误差（每一步会产生的误差）和整体误差（每一步的局部截断误差的累积）
研究这两个数没有太大意义，但是研究他们的阶（误差与△t的关系）有意义。
隐式欧拉方法是一阶的，也就是局部误差是O(h2)，整体误差是O(h)。这里h指的是取的步长△t。
如何理解O(h)？如果把h（△t）减小一半，那么期望误差也会跟着减小1/2。
对于 O(h2)来讲，如果把h（△t）减小一半，那么误差会减小到1/4。
如果有O(h3)，那么如果把h（△t）减小一半，那么误差会减小到1/8。
因此阶数越高误差越小。
4、Runge-Kutta Families
A family of advanced methods for solving ODEs
Especially good at dealing with non-linearityIt’s order-four version is the most widely used, a.k.a. RK4
龙格库塔方法特别适合解微分方程，尤其是对于非线性的。
RK4：
5、Position-Based / Verlet Integration
Idea:
After modified Euler forward-step, constrain positions of particles to prevent divergent, unstable behaviorUse constrained positions to calculate velocityBoth of these ideas will dissipate energy, stabilize
Pros / cons
Fast and simpleNot physically based, dissipates energy (error)
不基于物理的方法，通过调整不同位置使得满足某种限制。但是由于不是基于物理的方法，因此不满足能量守恒。
二、Rigid Body Simulation
Simple case
Similar to simulating a particleJust consider a bit more properties
对于刚体的模拟类似于对单个粒子的模拟，但是会更多的考虑刚体本身的物理量。
三、Fluid Simulation
A Simple Position-Based Method
Key idea
Assuming water is composed of small rigid-body spheresAssuming 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 particlesYou need to know the gradient of the density anywhere w.r.t. each particle’s positionUpdate? Just gradient descent!
整个水体是由很多不可压缩的刚体小球组成。通过模拟这些小球的运动位置，来模拟整个水体的运动，最后再渲染出来。
这里假设整个水是不可压缩的。也就是在任何一个时刻任何一个地方的密度相同。
通过任何一个时刻小球的分布都可以知道某个小球周围的密度。如果有任何一个地方的密度和之前平静的水在这个位置的密度不一样，那么就需要通过移动小球把这个地方的密度修正回来。通过不断的根据小球分布计算不同位置的密度，不断的修正，就可以将水模拟出来。这种方式是不基于物理的。
Eulerian vs. Lagrangian
Two different views to simulating large collections of matters
在物理中模拟大规模物质的两种方法：
质点法（拉格朗日方法）：如刚刚对水体的模拟，认为水是由各个圆形的小水滴组成的，通过模拟各个小水滴的运动来模拟整个的水体运动。
网格法（欧拉方法）：把空间划分成网格，计算不同网格中不同时间密度的变化来模拟。
Material Point Method (MPM)
Hybrid, combining Eulerian and Lagrangian views
Lagrangian: consider particles carrying material propertiesEulerian: use a grid to do numerical integrationInteraction: particles transfer properties to the grid, grid performs update, then interpolate back to particles
结合拉格朗日方法和欧拉法。
首先认为不同的粒子具有某些材质属性（如上图这是一只会融化的兔子，粒子间有粘性、粒子有质量等），然后融化过程通过网格去计算，算出来之后再把信息写回每个粒子中去。
