积分基础-如何积分运动方程

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;

Next:Fix Your Timestep!

原文链接