integration basics 번역

https://gafferongames.com/post/integration_basics/

Introduction

Hi, I’m Glenn Fiedler and welcome to Game Physics.

If you have ever wondered how the physics simulation in a computer game works then this series of articles will explain it for you. I assume you are proficient with C++ and have a basic grasp of physics and mathematics. Nothing else will be required if you pay attention and study the example source code.

A physics simulation works by making many small predictions based on the laws of physics. These predictions are actually quite simple, and basically boil down to something like “the object is here, and is traveling this fast in that direction, so in a short amount of time it should be over there”. We perform these predictions using a mathematical technique called integration.

Exactly how to implement this integration is the subject of this article.

물리가 컴퓨터게임에서 어떻게 시뮬레이션되는지 궁금했다면 다음 일련의 게시글들이 도움이 될 겁니다. 여러분이 C++에 능숙하며 물리와 수학에 대한 기초가 있다고 가정합니다. 집중해서 예제 소스코드를 공부하는 것 외에 특별히 더 필요한 것은 없습니다.

물리 시뮬레이션은 물리법칙에 근거해서 아주 작은 예측을 함으로써 이루어집니다. 이 예측이란 꽤나 간단해서, 결국은 “이 객체는 지금 여기있고, 저 방향으로 이 만큼 빠른 속도로 움직이고 있으니까 잠시 후엔 저기에 있을 거야”와 같다고 할 수 있습니다. 우린 이러한 예측을 수학적인 기법인 적분을 통해 수행합니다.

이 게시글의 주제는 바로 이 적분을 어떻게 구현할 것인가에 관한 내용입니다.

 

Integrating the Equations of Motion

You may remember from high school or university physics that force equals mass times acceleration.
  F = ma
We can switch this around to see that acceleration is force divided by mass. This makes intuitive sense because heavier objects are harder to throw.
  a = F/m
Acceleration is the rate of change in velocity over time:
  dv/dt = a = F/m
Similarly, velocity is the rate of change in position over time:
  dx/dt = v
This means if we know the current position and velocity of an object, and the forces that will be applied to it, we can integrate to find its position and velocity at some point in the future.

운동 방정식 미분하기

고등학교나 대학 물리에서 배운 힘은 질량 곱하기 가속도라는 걸 기억하실지 모르겠습니다.
F = ma
이 공식을 변형하면 가속도가 힘 나누기 질량이라는 걸 알 수 있습니다. 무거운 객체가 던지기 어렵다는 걸 생각하면 직관적으로 이해할 수 있습니다.
a = F/m
가속도는 시간에 따른 속도의 변화량입니다.
dv/dt = a = F/m
마찬가지로 속도는 시간에 따른 위치의 변화량입니다.
dx/dt = v
이렇게 객체의 현재 위치와 속도, 적용된 힘을 알고있다면 이를 적분해서 미래의 특정 지점에서의 위치와 속도를 알 수 있습니다.

Numerical Integration 수치 적분

For those who have not formally studied differential equations at university, take heart for you are in almost as good a position as those who have. This is because we’re not going to analytically solve the differential equations as you would do in first year mathematics. Instead, we are going to numerically integrate to find the solution.

대학교에서 미분 방정식을 배우지 않은 분들도 걱정할 필요가 없습니다. 우리는 대학교 첫 해에 하듯 분석적으로 미분 방정식을 풀 것이 아니라 수치적으로 해결책을 찾을 것이기 때문입니다.

Here is how numerical integration works. First, start at an initial position and velocity, then take a small step forward to find the velocity and position at a future time. Then repeat this, moving forward in small time steps, using the result of the previous calculation as the starting point for the next.

수치 적분법은 다음과 같이 이루어집니다. 먼저 초기 위치와 속도로 시작하여, 작은 스텝을 밟아 미래 시간에서의 속도와 위치를 찾습니다. 이를 반복하여 여러 스텝들을 밟아 나가며 이전 계산의 결과를 다음 계산의 시작 지점으로 삼습니다.

 

But how do we find the change in velocity and position at each step?

The answer lies in the equations of motion.

Let’s call our current time t, and the time step dt or ‘delta time’.

We can now put the equations of motion in a form that anyone can understand:

그런데 각 스텝에서의 속도와 위치의 변화량은 어떻게 알 수 있을까요?

답은 운동 방정식에 들어있습니다.

현재 시간을 t, 타임 스텝을 dt 또는 delta time이라고 합시다.

이제 운동 방정식을 누구나 이해할 수 있는 식으로 나타낼 수 있습니다.

   acceleration = force / mass
   change in position = velocity * dt
   change in velocity = acceleration * dt

This makes intuitive sense because if you’re in a car traveling 60 kilometers per-hour, in one hour you’ll be 60 kilometers further down the road. Similarly, a car accelerating 10 kilometers per-hour-per-second will be moving 100 kilometers per-hour faster after 10 seconds.

이는 직관적으로 이해할 수 있습니다. 만약 시속 60킬로미터로 차가 움직인다면, 한 시간안에 60킬로미터의 길을 지나게 될 것입니다. 마찬가지로 초당 시속 10킬로미터로 차가 가속한다면 10초후에는 시속 100킬로미터 더 빨리 움직일 것입니다.

Of course this logic only holds when acceleration and velocity are constant. But even when they’re not, it’s still a pretty decent approximation to start with.

물론 이 로직은 가속도와 속도가 상수일 때만 가능합니다. 근데 상수가 아니더라도 꽤나 근사한 값을 제공합니다.

Let’s put this into code. Starting with a stationary object at the origin weighing one kilogram, we apply a constant force of 10 newtons and step forward with time steps of one second:

이제 이걸 코드로 표현해봅시다. 무게 1 킬로그램의 멈춰있는 객체에 10N의 상수 힘을 적용해서 일 초의 타임 스텝을 두고 나아간다고 해봅시다.

   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;
   }

Here is the result:

   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

As you can see, at at each step we know both the position and velocity of the object. This is numerical integration.

이처럼 각 스텝에서의 객체의 위치와 속도를 모두 알 수 있습니다. 이를 수치 적분이라고 합니다.

Explicit Euler

What we just did is a type of integration called explicit euler.

우리가 방금 수행한 적분법을 explicit euler(explicit 오일러 방법) 라고 합니다.

To save you future embarrassment, I must point out now that Euler is pronounced “Oiler” not “yew-ler” as it is the last name of the Swiss mathematician Leonhard Euler who first discovered this technique.

미래에 있을 망신을 피하기 위해 여러분에게 Euler는 “오일러” 라고 발음하지 “율러”가 아니라는 걸 알려드리고자 합니다. 이는 이 기법을 발견한 스위스 수학자 레온하르트 오일러의 이름입니다.

Euler integration is the most basic numerical integration technique. It is only 100% accurate when the rate of change is constant over the timestep.

오일러 적분은 제일 기본적인 수치 적분법입니다. 이는 타임 스텝에 따른 변화량이 상수일 때에만 100% 정확합니다.

Since acceleration is constant in the example above, the integration of velocity is without error. However, we are also integrating velocity to get position, and velocity is increasing due to acceleration. This means there is error in the integrated position.

위의 예제에서 가속도가 상수이므로, 속도에 대한 적분은 오류가 없습니다. 하지만 우리는 속도를 적분하여 위치를 구하고 있고, 속도는 가속도에 의해 증가합니다. 이는 적분된 위치에 오류가 있다는 걸 의미합니다.

Just how large is this error? Let’s find out!

이 오류가 얼마나 큰 걸까요? 알아내 봅시다.

There is a closed form solution for how an object moves under constant acceleration. We can use this to compare our numerically integrated position with the exact result:

객체가 등가속도 운동을 통해 어떻게 움직이는지에 대한 이미 알려진 공식을 사용해서 정확한 값과 우리가 수치 적분으로 얻은 값을 비교해보도록 하겠습니다.

   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

After 10 seconds, the object should have moved 500 meters, but explicit euler gives a result of 450 meters. That’s 50 meters off after just 10 seconds!

This sounds really, really bad, but it’s not common for games to step physics forward with such large time steps. In fact, physics usually steps forward at something closer to the display framerate.

Stepping forward with dt = 1⁄100 yields a much better result:

10초 후에, 객체는 500미터를 움직였어야 하지만 explicit euler의 결과는 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

With a small enough timestep explicit euler gives decent results for constant acceleration, but what about the case where acceleration isn’t constant?

A good example of non-constant acceleration is a spring damper system.

In this system a mass is attached to a spring and its motion is damped by some kind of friction. There is a force proportional to the distance of the object that pulls it towards the origin, and a force proportional to the velocity of the object, but in the opposite direction, which slows it down.

Now the acceleration is definitely not constant throughout the timestep, but is a continously changing function that is a combination of the position and velocity, which are themselves changing continuously over the timestep.

This is an example of a damped harmonic oscillator. It’s a well studied problem and there’s a closed form solution that we can use to check our numerically integrated result.

Let’s start with an underdamped system where the mass oscillates about the origin while slowing down.

Here are the input parameters to the mass spring system:

충분히 작은 타입스텝일 때 explicit euler는 상수 가속도에 대한 괜찮은 결과를 내고 있습니다만, 상수 가속도가 아닌 경우는 어떨까요?

비상수 가속도에 대한 좋은 예로는 spring damper system이 있습니다.

이 시스템에서 질량은 스프링에 부착되어 일종의 마찰력에 의해 운동이 감쇠(damped)되어집니다.  객체가 원점으로부터 당겨진 거리에 비례하는 힘, 객체의 속도에 비례하는 힘이 있고, 이를 느려지게 하는 힘이 반대방향으로 작용합니다.

이제 타임 스텝을 밟는 동안 가속도는 확실히 상수가 아니며, 타임 스텝동안 계속해서 변화하는 위치와 속도의 조합에 의해 계속해서 변하는 함수가 됩니다.

예를 들면 damped harmonic oscillator가 있습니다. 이는 잘 연구된 문제로 우리의 수치 적분 결과 확인에 사용할만한 공식이 존재합니다.

먼저 질량이 느려지면서 진원에 대하여 진동하는 저감쇠 시스템으로 시작해봅시다.

다음은 질량 스프링 시스템에 대한 입력 인자들입니다 :

• Mass: 1 kilogram

• Initial position: 1000 meters from origin

• Hooke’s law spring coefficient: k = 15

• Hooke’s law damping coefficient: b = 0.1

And here is a graph of the exact solution:

정확한 결과에 대한 그래프는 다음과 같습니다.

When we apply explicit euler to integrate this system, we get the following result, which has been scaled down vertically to fit:

explicit 오일러로 이 시스템을 적분하면 다음과 같은 결과가 나오는데, 틀에 맞추기 위해 수직적으로 축소된 상태입니다.

Instead of damping and converging on the origin, it gains energy over time!

감쇠를 하며 근원으로 수렴하기보다, 시간에 따라 에너지를 얻고있습니다!

이 시스템은 explict 오일러와 dt=1/100로 적분하면 불안정해집니다.

안타깝지만 우리가 이미 작은 타임 스텝으로 적분하고 있기 때문에 정확성을 올릴만한 실용적이 옵션은 많지 않습니다. 여러분이 타임 스텝을 더 줄인다고 하더라도 위와같은 스프링 tightness 계수 k를 보게될 겁니다.

Semi-implicit Euler 반 암시적 오일러

Another integrator to consider is semi-implicit euler.

Most commercial game physics engines use this integrator.

Switching from explicit to semi-implicit euler is as simple as changing:

또다른 적분법으로는 반 암시적 오일러가 있습니다.

대부분의 상용 게임 물리 엔진은 이 적분법을 사용합니다.

explicit에서 semi-implicit 오일러로 바꾸는 건 아주 간단합니다.

   position += velocity * dt;
   velocity += acceleration * dt;
to:

   velocity += acceleration * dt;
   position += velocity * dt;

 

Even though semi-implicit euler has the same order of accuracy as explicit euler (order 1), we get a much better result when integrating the equations of motion because it is symplectic.

비록 semi-implicit 오일러가 explicit 오일러와 같은  order of accuracy 를 가지지만, 사교(symplectic)하기에 운동 방정식을 적분할 때 좀 더 좋은 결과를 냅니다.

 

Many different integration methods exist

And now for something completely different.

Implicit euler is an integration technique that is well suited for simulating stiff equations that become unstable with other methods. The drawback is that it requires solving a system of equations per-timestep.

Implicit 오일러는 다른 방법으로는 불안정한 stiff 방정식을 시뮬레이팅 할 때 적합한 적분법입니다. 단점은 타임스텝마다 방정식을 풀어야 한다는 것입니다.

Verlet integration provides greater accuracy than implicit euler and less memory usage when simulating a large number of particles is. This is a second order integrator which is also symplectic.

벨벳 적분법은 많은 수의 파티클들을 시뮬레이팅할 때 Implicit 오일러보다 훨씬 정확하고 더 적은 메모리를 사용합니다. 이는 symplectic한 second order의 적분법입니다.

There is a whole family of integrators called the Runge-Kutta methods. In fact, explicit euler is considered part of this family, but it also includes higher order integrators, the most classic of these being the Runge Kutta order 4 or simply RK4.

또한 룽게-쿠타 방법이라고 불리는 군의 적분기들이 있습니다. explicit 오일러도 이 군에 속하지만, 좀더 높은 order의 적분기이자 제일 클래식한 적분기인 룽게 쿠타 4 또는 RK4도 포함하고 있습니다.

This family of integrators is named for the German physicists who discovered them: Carl Runge and Martin Kutta. This means the ‘g’ is hard and the ‘u’ is a short ‘oo’ sound. I am sorry to inform but this means we are talking about the ‘roon-geh koo-ta’ methods and not a ‘runge cutter’, whatever that is :)

(룽게 쿠타 이름의 유래)

The RK4 is a fourth order integrator, which means its accumulated error is on the order of the fourth derivative. This makes it very accurate. Much more accurate than explicit and implicit euler which are only first order.

RK4는 4차 적분기로, 누적된 오차가 4차 도함수의 order에 놓인다는 뜻입니다. (???) 덕분에 매우 정확해집니다. 오직 1차일 뿐인 explicit/implicit 오일러보다 더 정확합니다.

But although it’s more accurate, that’s not to say RK4 is automatically “the best” integrator, or even that it is better than semi-implicit euler. It’s much more complicated than this. Regardless, it’s an interesting integrator and is well worth studying.

그러나 더 정확하다고 해서 RK4가 자동적으로 “최고’의 적분기가 된다거나, semi-implicit 오일러보다 더 낫다고 할 수는 없습니다. 좀 더 복잡하게 따져봐야 합니다. 그렇지만 상당히 흥미로운 적분기이며 공부할 가치가 충분합니다.

Implementing RK4

There are many great explanations of the mathematics behind RK4 already. For example: here, here and here. I highly encourage you to follow the derivation and understand how and why it works at a mathematical level. But, seeing as the target audience for this article are programmers, not mathematicians, we’re all about implementation here. So let’s get started.

Before we go any further let’s define the state of an object as a struct in C++ so we have both position and velocity stored conveniently in one place:

  struct State
   {
       float x;      // position
       float v;      // velocity
   };

We also need a struct to store the derivatives of the state values:

  struct Derivative
   {
       float dx;      // dx/dt = velocity
       float dv;      // dv/dt = acceleration
   };

Next we need a function to advance the physics state ahead from t to t+dt using one set of derivatives, and once there, recalculate the derivatives at this new state:

  Derivative evaluate( const State & initial,
                        double t,
                        float dt,
                        const Derivative & d )
   {
       State state;
       state.x = initial.x + d.dx*dt;
       state.v = initial.v + d.dv*dt;

       Derivative output;
       output.dx = state.v;
       output.dv = acceleration( state, t+dt );
       return output;
   }

The acceleration function is what drives the entire simulation. Let’s set it to the spring damper system and return the acceleration assuming unit mass:

  float acceleration( const State & state, double t )
   {
       const float k = 15.0f;
       const float b = 0.1f;
       return -k * state.x - b * state.v;
   }

What you write here is of course simulation dependent, but you must structure your simulation so you can calculate the acceleration inside this method given the current state and time, otherwise it won’t work with the RK4 integrator.

Finally we get to the integration routine itself:

  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;
   }

The RK4 integrator samples the derivative at four points to detect curvative. Notice how derivative a is used when calculating b, b is used when calculating c, and c into d. This feedback of the current derivative into the calculation of the next is what gives the RK4 integrator its accuracy.

Importantly, each of these derivatives a,b,c and d will be different when the rate of change in these quantities is a function of time or a function of the state itself. For example, in our spring damper system acceleration is a function of the current position and velocity which change throughout the timestep.

Once the four derivatives have been evaluated, the best overall derivative is calculated as a weighted sum derived from the taylor series expansion. This combined derivative is used to advance the position and velocity forward, just as we did with the explicit euler integrator.

Semi-implicit euler vs. RK4

Let’s put the RK4 integrator to the test.

Obviously, since it is a higher order integrator (4th order vs. 1st order) it will be visibly more accurate than semi-implicit euler, right?

Wrong. Both integrators are so close to the exact result that it’s impossible to make out any difference at this scale. Both integrators are stable and track the exact solution very well with dt=1⁄100.

Zooming in confirms that RK4 is more accurate than semi-implicit euler, but is it really worth the complexity and extra runtime cost of RK4? It’s hard to say.

Let’s push a bit harder and see if we can find a significant difference between the two integrators. Unfortunately, we can’t look at this system for long periods of time because it quickly damps down to zero, so let’s switch to a simple harmonic oscillator which oscillates forever without any damping.

Here’s the exact result we’re aiming for:

To make it harder on the integrators, let’s increase delta time to 0.1 seconds.

Next, we let the integrators run for 90 seconds and zoom in:

After 90 seconds the semi-implicit euler solution (orange) has drifted out of phase with the exact solution because it has a slightly different frequency, while the green line of RK4 matches the frequency, but is losing energy!

We can see this more clearly by increasing the time step to 0.25 seconds.

RK4 maintains the correct frequency but loses energy:

While semi-implicit euler does a better job at conserving energy, on average:

But drifts out of phase. What an interesting result! As you can see it’s not simply the case that RK4 has a higher order of accuracy and is “better”. It’s much more nuanced than this.

Conclusion

We’ve implemented three different integrators and compared their results.

1. Explicit euler

2. Semi-implicit euler

3. Runge Kutta order 4 (RK4)

So which integrator should you use in your game?

My recommendation is semi-implicit euler. It’s cheap and easy to implement, it’s much more stable than explicit euler, and it tends to preserve energy on average even when pushed near its limit.

If you really do need more accuracy than semi-implicit euler, I recommend you look into higher order symplectic integrators designed for hamiltonian systems. This way you’ll discover more modern higher order integration techniques that are better suited to your simulation than RK4.

And finally, if you are still doing this in your game:

  position += velocity * dt;
   velocity += acceleration * dt;

Please take a moment to change it to this:

  velocity += acceleration * dt;
   position += velocity * dt;

You’ll be glad you did :)