积分基础-如何积分运动方程
Integration Basics
How to integrate the equations of motion
Posted by Glenn Fiedler on Thursday,June 1,2004
Introduction
Hi, I’m Glenn Fiedler and welcome to Game Physics.
如果你曾经想知道pc游戏的物理模拟是如何工作的，那么这一系列的文章会告诉你。我假设你熟练使用C++，并且有一定的物理和数学基础。你学习源码将不会有额外的困难。
一个物理模拟，由多个基于物理定理的小的预测组成的。这些预测是非常简单的，并且大体上归结为类似“对象应该在这里，并且他向那个方向快速移动，因此过一段时间后，他应该在那”的东西。我们用叫做积分的数学工具来表现这些预测。
确切来说，如何实现这些积分，就是这篇文章的主题。
Integrating the Equations of Motion
你可能记得，高中或大学的物理课中，力等于质量乘加速度。
F=ma
我们可以看到，加速度等于力除质量。这样更符合直觉，因为较重的物理更难被扔出去。
a = F/m
加速度就是速度随时间的变化率：
dv/dt = a = F/m
相同的，速度就是位置随时间的变化率：
dx/dt = v
这意味着如果我们知道对象当前的位置和速度，和将要对其作用的力，我们可以求的它在将来某点的位置和速度。
Numerical Integration
对于那些没有在大学正式学习过微分方程的人，请放心，因为那些在大学里学习过微分方程的人并没有比你更有优势。这是因为我们将不会去分析解决微分方程。我们将进行数值化积分来解决问题。
数值化积分是这么工作的。首先，设定一个初始的位置和速度，然后向前一小步来找到这个未来点的速度和位置。然后重复这个过程，向前移动一个小时间段，使用计算结果作为下一步的起点。
但是我们如何算出每一步的速度和位置变化。
答案就在运动方程里（equations of motion）。
我们设当前时间为t,时间段为dt或者“delta time”。
现在我们可以把它放进运动方程，来帮助我们理解：
acceleration = force / mass
change in position = velocity * dt
change in velocity = acceleration * dt
这很符合直觉，因为如果你在一个每小时60公里的车上，1个小时你会沿着道路移动60公里。相同的，一辆车每秒加速10千米/小时，那么10秒后，它的速度是100千米/小时。
当然，这个逻辑只在速度和加速度恒定时成立。但是即使他们不是恒定的，作为一个开始也是一个好的近似结果。
让我们把它放进代码里。一个固定的物体，初始质量是1kg，我们提供10n的恒定力，向前推进的时间段设定为1s：
double t = 0.0;
float dt = 1.0f;
float velocity = 0.0f;
float position = 0.0f;
float force = 10.0f;
float mass = 1.0f;
while( t <= 10.0 ){
position = position + velocity * dt;
velocity = velocity + ( force / mass ) * dt;
t += dt;
}
下面是结果：
t=0: position = 0 velocity = 0
t=1: position = 0 velocity = 10
t=2: position = 10 velocity = 20
t=3: position = 30 velocity = 30
t=4: position = 60 velocity = 40
t=5: position = 100 velocity = 50
t=6: position = 150 velocity = 60
t=7: position = 210 velocity = 70
t=8: position = 280 velocity = 80
t=9: position = 360 velocity = 90
t=10: position = 450 velocity = 100
你可以看到，每一段我们可以知道对象的位置和速度。这就是数值化积分。
Explicit Euler
我们刚才做的是一种类型的数值化积分，叫做显式欧拉（explicit euler）。
为了让你们在将来不会尴尬，我要声明Euler的发音是“Oiler” 不是“yew-ler”这是第一个发明这个方法的瑞典数学家Leonhard Euler的姓。
欧拉积分是最基础的数值化积分技术。它只在变化率恒定时，百分百准确。
因为在上述例子中，加速度是恒定的，所以速度的结果是正确的。然而，我们使用了速度的积分来求出位置，速度是随时间增加的因为有加速度。这意味着位置的积分结果是错误的。
那么有多大的偏差呢？让我们来求出它！
对于物体在恒定加速度下的位移，又一个封闭形式的解。我们可以用这个方程的结果，来比较我们得出的结果：
s = ut + 0.5at^2
s = 0.0*t + 0.5at^2
s = 0.5(10)(10^2)
s = 0.5(10)(100)
s = 500 meters
在10秒之后，这个对象应该移动了500m，但是显式欧来给出的结果是450米。只用了10秒就产生了50米的偏移。
这听起来很糟糕，但是在常规游戏中，前向步进物理计算，不会使用这么长的步长。实际上，物理计算使用的步长接近于显式频率。
前向步进步长dt = 1⁄100 会得到一个更好的结果：
t=9.90: position = 489.552155 velocity = 98.999062
t=9.91: position = 490.542145 velocity = 99.099060
t=9.92: position = 491.533142 velocity = 99.199059
t=9.93: position = 492.525146 velocity = 99.299057
t=9.94: position = 493.518127 velocity = 99.399055
t=9.95: position = 494.512115 velocity = 99.499054
t=9.96: position = 495.507111 velocity = 99.599052
t=9.97: position = 496.503113 velocity = 99.699051
t=9.98: position = 497.500092 velocity = 99.799049
t=9.99: position = 498.498077 velocity = 99.899048
t=10.00: position = 499.497070 velocity = 99.999046
如你所见，这是一个很好的运算结果。对游戏来说足够了。
Why explicit euler is not (always) so great
在足够小的步长下，显式欧拉给出的结果很接近加速度恒定下实际的结果，但是如果加速度也不恒定呢？
一个很好的非恒定加速度的例子就是弹簧减震系统（spring damper system）。
在这个系统中，一个负载链接在弹簧上，它的动能通过一个参数来减小。有一个力与物体的距离成比例，将物体拉向原点，还有一个力与物体的速度成比例，但方向相反，使物体减速。
现在加速度一定不是随时间恒定的了，但是是一个不断变化的函数，它是位置和速度的组合，它们本身随着时间的推移不断变化。
这是一个阻尼谐振子的例子（[damped harmonic oscillator][]）。这是一个很好的研究对象，并且他有一个封闭的解，来检查我们的计算结果。
让我们从一个欠阻尼系统开始，它的振子将会减速。
这是弹簧阻尼系统的输入参数：
-质量：1kg
-初始位置：距离原点1000m
-虎克定律弹簧系数：k = 15
-虎克定律弹簧阻尼系数： b = 0.1
这是结果图像：
当我们在这个系统中使用显式欧拉，我们将得到如下结果，已经改变了比例：
不是随时间减速和接近原点，而是随时间获得能量！
这个系统在显式欧拉下并且dt=1/100时，是不稳定的。
不幸的是，因为我们已经使用了小的时间步长，我们没有很多实际的选择来提高准确度。即使你减小时间步长，这之上总有一个弹簧系数，你会看见它的影响。
Semi-implicit Euler
另一个可以考虑的方法是半隐式欧拉（semi-implicit euler）
很多商业的游戏物理引擎，使用这个方法。
将显式转化为半隐式，是很简单的：
position += velocity * dt
velocity += acceleration * dt
to:
velocity += acceleration * dt
position += velocity * dt
为半隐式欧拉积分器提供dt=1/100到弹簧阻尼系统里，结果就很接近实际结果了：
尽管半隐式欧拉有着和显式一个等级的精度（order 1），我们获得了一个更好的运动方程结果，因为它是symplectic
Many different integration methods exist
Implicit euler是一种积分技术，非常适合模拟在其他方法下变得不稳定的刚性方程。缺点是，它需要每一个时间步解一个方程组。
Verlet integration提供了比显式欧拉更高的精度，和更少的内存占用，在计算大量的粒子时。它是一个二阶积分器，并且也是symplectic。
积分器家族有个名称叫做Runge-Kutta methods。隐式欧拉是这个家族的一部分，但是它也包含了更高阶的积分器。最经典的是 Runge Kutta order 4 或者缩写为 RK4。
Runge Kutta积分器家族被发现他们的德国物理学家 Carl Runge 和 Martin Kutta 命名。
RK4是一个四阶积分器，也就是说它的累积误差是四阶导数。这让它很精确。比显式和隐式欧拉好很多，他们只有一阶。
但即使更精确，这也不意味RK4是默认“最好的”积分器，或者它比半隐式欧拉好。这个问题很复杂。
不管怎样，他说一个有趣的积分器并且它值得研究。
Implementing RK4
RK4早就有很多很棒的数学描述了。比如：here,here and here.我强烈建议你遵循推导法，理解它是如何和为什么在数学水平上工作的。但是，基于这篇文章主要面向程序员而不是数学家，我们只介绍实现。
让我们在C++里定义一个对象的状态，因此我们要把位置和速度储存在一起：
struct State{
float x; //position
float v; //velocity
};
我们也需要一个结构体来储存导数状态：
struct Derivative{
float dx // dx/dt = velocity
float dv; // dv/dt = acceleration
};
接下来，我们需要一个函数，用一组导数将物理状态从t提前到t+dt，然后在那里重新计算这个新状态下的导数。
Derivative evaluate(const State & initial,
double t,
float dt,
const Derivative & d){
State state;
state.x = inital.x + d.dx*dt;
state.v = inital.v + d.dv*dt;
Derivative output;
output.dx = state.v;
output.dv = acceleration( state, t+dt );
retrun output;
}
加速度方程式是用来驱动整个模拟的。让我们把它设定在弹簧阻尼系统里，并设定质量为单位质量：
float acceleration( const State & state, double t ){
const float k = 15.0f;
const float b = 0.1f;
return -k * state.x - b * state.v;
}
现在让我们实现RK4积分器：
void integrate( State & state,
double t,
float dt)
{
Derivative a,b,c,d;
a = evaluate( state, t, 0.0f, Derivative() );
b = evaluate( state, t, dt*0.5f, a );
c = evaluate( state, t, dt*0.5f, b );
d = evaluate( state, t, dt, c);
float dxdt = 1.0f / 6.0f *
( a.dx + 2.0f * ( b.dx + c.dx ) + d.dx );
float dvdt = 1.0f / 6.0f *
( a.dv + 2.0f * ( b.dv + c.dv ) + d.dv );
state.x = state.x + dxdt * dt;
state.v = state.v + dvdt * dt;
}
RK4 积分器在4个点采样，来检测曲线。注意在计算b时，是如何使用倒数的，b在计算c时被使用，c被d使用。这种导数反馈，就是RK4精确的原因。
重要的是，当变化率在等式中是时间的方程或者是状态的方程时，a,b,c,d都会不同。例如，在我们的弹簧阻尼系统中，加速度是相对于当前位置和速度变化的，而位置速度又随时间变化。
一旦对四个导数进行了评估，那么最好的整体导数是作为泰勒级数（taylor series）展开的加权和进行计算。这种组合导数被用来向前推进位置和速度，就像我们对显式欧拉积分所做的那样。
Semi-implicit euler vs. RK4
现在让我们把 RK4 积分器加入测试。
因为它是高阶积分器（4阶 vs 1阶）它会比半隐式欧拉更精确，对么？
警告。两个积分器都非常接近精确的结果，所以在这个尺度上是不可能有任何区别的。两个积分器都是稳定的而已符合实际结果在dt=1/100。
放大来确认 RK4 比半隐式欧拉更精确，但是使用更复杂更耗时的RK4值得么？这很难说。
让我们加大难度，来看看我们是否能找到两个积分器重要的差别。不幸的是，我们不能在这个系统观察很长时间，因为它很快就会衰减到0，因此让我么改为简谐振器（simple harmonic oscillator）,其会永远震荡不会衰减。
这是目标实际结果：
为了让积分器更难，让我们增加 delta time 到 0.1s。
然后，我们让积分器运行90s并方法结果：
在90s后，半隐式欧拉的结果（橙色）由于频率稍有不同，所以与精确解偏离了相位，绿色的 RK4 符合频率，但是在丢失能量！
当我们把dt增加到0.25s，就更明显了。
RK4 保留了正确的频率，但是丢失了能量：
平均来说，半隐式欧拉在保存能量上做的更好：
但是偏离了相位。多么有趣的结果！如你所见这不是简单的，RK4有更高的精度，那么它就“更好”。其中有很多细微的差别。
Conclusion
你应该使用在你游戏里使用哪个积分器？
我的推荐是 semi-implicit euler。它廉价并且易于实现，它比显式欧拉更稳定，即使在接近极限时，它也倾向于平均保存能量。
如果你真的需要比半隐式欧拉更精确的积分器，我推荐你看被设计用于 hamiltonian systems 的更高阶的 symplectic integrators。这样，您将发现比RK4更适合模拟的现代高阶积分技术。
最后，如果你在游戏里仍然这么写：
position += velocity * dt;
velocity += acceleration * dt;
请花点时间，改成如下：
velocity += acceleration * dt;
position += velocity * dt;
