The Performance of Numerical Integrators with Varying Delta Times

I always like using a straightforward method for the main game loop and the ticking of tasks. Accordingly, the main game loop merely provides the delta-time (the time between the current and the previous iteration of the loop) to all the different tasks. The following snippet shows the basic outline of how I'm handling the primary game loop and ticking of tasks:

// Returns the current state of the high resolution clock in nanoseconds
uint64_t ns();
// Returns the minimum of the provided floats
float min(float a, float b);
// Stub task functions
void task0(float dt); void task1(float dt);

static void execute_tasks()
{
  static uint64_t lastUpdateNs = ns();
  const uint64_t currentNs = ns();

  // Calculate the delta time and clamp it to a maximum of 10 frames per second
  const float dt = min((currentNs - lastUpdateNs) * 1e-9f, 0.1f);
  // Update the cached timestamp for the next iteration
  lastUpdateNs = currentNs;

  // Execute tasks using the given delta-time
  {
    task0(dt);
    task1(dt);

    // ...
  }
}

static void game_loop()
{
  // Extremely simplified game loop
  while (true)
  {
    execute_tasks();
  }
}

Remember that it's always a good practice to utilize the highest resolution that the system clock provides. If you've got access to nanosecond precision, that's great, so let's make sure to use it.

Since the delta-time can and most certainly will differ from iteration to iteration, this is what you call a variable delta-time or timestep. Using a variable timestep works well in specific scenarios, like, e.g., incrementing timers or simple simulations that still provide a deterministic output if the delta-time fluctuates or is not necessarily small.

Sadly, it can be hard to determine if a specific piece of code will behave deterministically independent of the current delta-time. Accordingly, I always follow this rule of thumb: switch to fixed time-stepping for gameplay-critical tasks and whenever the visual results of non-critical systems differ between high and low framerates. Luckily it's straightforward to extend the proposed game loop to a hybrid that allows us to effectively decouple the simulation of various systems from the variable rendering frame rate.

To better understand the issues that arise from a varying delta-time, in the following section, we will take a closer look at some numerical integration schemes and how precisely they operate under varying conditions.

The Performance of Numerical Integration Schemes

This section will take a closer look at some typical numerical integration schemes that should be in everyone's toolbox. Then, using a simple example, we will highlight the ups and downs, making it easier to make the right choice when working on gameplay and simulation features.

Let's start with a straightforward one-dimensional particle simulation as an example first. Let $\vec{p} = 0$ be the position, $\vec{a} = -9.81$ the constant acceleration, and $\vec{v} = 0$ the initial velocity of a single particle. Using a numerical integrator, it's possible to incrementally approximate $\vec{p}(t)$ and $\vec{v}(t)$ at any point in time $t$ in the future.

The accuracy of the integration heavily depends on the chosen time step and the order of the integration scheme. The global error of a first-order integration scheme is directly proportional to the delta-time $O(\Delta t)$. In our context, the global error defines the deviation of the estimated position $\vec{p}(t)$ compared to the ground truth, summing the local errors of all the steps needed to reach this point. Accordingly, the global error of a second-order integration scheme is proportional to the square of the delta-time $O(\Delta t^2)$. So if we go ahead and halve the delta-time used with a first-order integrator, the error approximately decreases by a factor of two. Suppose we halve the delta-time used in conjunction with a second-order integrator the error roughly reduces by a factor of four.

The simplest form of numerical integration is called Euler's Method or Explicit Euler Integration. Its simplicity and efficiency probably make it the most used scheme today. But is it the best option available? Let's take a closer look!

Explicit Euler Integration (First Order)

Explicit Euler Integration is a first-order integration scheme that utilizes extrapolation based on the current frame's velocity to predict the position and velocity at the next point in time $t + \Delta t$:

A simple implementation in C++ can look like this:

template<typename T>
static void integrate_euler_explicit(T& p, T& v, T a, float dt)
{
  p = p + v * dt;
  v = v + a * dt;
}

The main downside of this scheme is that calculating the position is only accurate if the velocity has not changed from $t$ to $t + \Delta t$. If the velocity is changing, which is the case in our constant acceleration example, it strongly tends to overshoot the correct result.

Luckily, due to the fact that we use a constant acceleration, the velocity term is consistently accurate regardless of the chosen timestep $\Delta t$.

As already mentioned, with a first-order integrator, the error in position increases approximately linearly with the timestep size. Here are the results of our single-particle example at differently sized timesteps $\Delta t$ ranging from $\frac{1}{1000}\textrm{ s}$ to $\frac{1}{20}\textrm{ s}$:

It's evident that this integrator only works well in our example using small timesteps in the range of $\Delta t < \frac{1}{60}\textrm{ s}$. Moreover, for all larger timesteps, the tendency to overshoot is visible. Even worse, the error increases visibly over time. After a simulation of $0.5\textrm{ s}$ at $\Delta t = \frac{1}{20}\textrm{ s}$ we end up with an error of $\approx 0.11\textrm{ m}$ at the particle's position. Simulating for $10\textrm{ s}$, the error increases to a shocking $\approx 2.439\textrm{ m}$.

The performance is far below our expectations. Especially since we're only inspecting an oversimplified simulation, as soon as the complexity of the underlying simulation extends, the errors will be even more noticeable at reasonably sized timesteps. One might argue that positional errors don't matter concerning a particle simulation as long as the results are smooth. But given the fact that we might very well use timesteps ranging from $\frac{1}{144}\textrm{ s}$ to $\frac{1}{30}\textrm{ s}$ as well as features based on the particle's position, like, e.g., vector fields, the outcome of the simulations will vary drastically. In addition to all of this, the integrator also misses another vital property: stability. Physically-based simulations, like, e.g., a mass-spring system, tend to explode completely. We will discuss this further in one of the following sections.

It's pretty simple to advance this integrator to the so-called semi-implicit variant of Euler integration. Let's check if this allows us to improve the situation.

Semi-Implicit Euler Integration (First Order, Symplectic)

By using $\vec{v}(t + \Delta t)$ instead of $\vec{v}(t)$ to calculate the position $\vec{p}(t + \Delta t)$ it's possible to transform the explicit Euler integration to its semi-implicit variant:

Looking at the implementation, this merely means switching the two lines advancing $\vec{p}$ and $\vec{v}$ and thus simply updating $\vec{v}$ before $\vec{p}$:

template<typename T>
static void integrate_euler_semi_implicit(T& p, T& v, T a, float dt)
{
  v = v + a * dt;
  p = p + v * dt;
}

The term implicit in the context of numerical integration describes the fact that we're using results from future points in time. Since we're only using the future velocity to calculate the position, which stills uses the current frame's acceleration, it remains only partially implicit and hence the term semi-implicit.

Let's examine the results of the integrator for our simple particle example:

This time there is a tendency to undershoot and overshoot. While the average resulting positions describe the trajectory better, the end-position errors are in the same ballpark as for the Explicit Euler Method, yielding $\approx 0.11\textrm{ m}$ after simulating for $0.5\textrm{ s}$ and $\approx 2.439\textrm{ m}$ after $10\textrm{ s}$, once again using a relatively big timestep of $\frac{1}{20}\textrm{ s}$. Given that this is still a first-order integrator, the accuracy once again increases linearly with reducing the timestep size.

This form of integration is also known as Symplectic Euler and, as the name already suggests, this makes it the first symplectic integrator depicted in this article. Symplectic integrators enable us to solve Hamilton's equations numerically with certain guarantees about the outcome. In detail, they guarantee to conserve the 2-form $dp \wedge dq$ where $dp$ and $dq$ depict the infinitesimal change in position and momentum and $\wedge$ the exterior/wedge product. Respectively the results are confined to a parallelogram with sides $dp$ and $dq$. Given the nature of the parallelogram, using a symplectic integrator might perturb the results, but the area defined by the resulting positions and momentums will remain stable in phase space.

A harmonic oscillator can be used as a simple example from classical mechanics to show the strength of implicit integration schemes. For a given position $\vec{p}$ the harmonic oscillator introduces a restoring force relative to a rest position. Assuming the rest position at the origin, the force is calculated as follows:

$F = -k\vec{p}$

Assuming the constant is set to $k = 1$, a particle mass $m = 1$ and a initial position of $\vec{p} = -1$ we can directly use the force as acceleration in our integration scheme. The closed-form solution for this particular setup is

$f_{\vec{p}}(t) = \cos(t + \pi)$

$f_{\vec{v}}(t) = \sin(t + \pi).$

Plotting the results with $\Delta t = \frac{1}{20}\textrm{ s}$ this yields the following results for the semi-implicit and explicit Euler variants:

We're quickly starting to drift off the exact position using the explicit Euler scheme. Instead of conserving it, the system is gaining energy, which is an absolute no-go in the context of physics simulations. However, using the semi-implicit Euler scheme, the results are almost spot-on. To make the error even more apparent, we can go ahead and plot the results in phase space, plotting the velocity along the x-axis and the position along the y-axis:

While the results derived from the semi-implicit scheme closely follow the expected ellipse, the results from the explicit method, on the other hand, spiral away from the exact solution with an increasing error as time goes by.

Here's an even more extreme example utilizing a timestep $\Delta t = \frac{1}{4}\textrm{ s}$:

This example clearly shows off the area and thus energy-conserving properties of the semi-implicit scheme. The results of the explicit scheme, on the other hand, quickly diverge from the expected results and gain even more energy as in the previous example.

Modified Euler Integration / Midpoint Method (Second Order)

Let's have a look at our first second-order integrator, the Modified Euler Integration or Midpoint Method. This method provides relatively high accuracy at the cost of a bit more computational complexity.

Instead of extrapolating the position based on the current velocity, we're advancing $\vec{v}$ half a timestep based on the current acceleration $\vec{a}(t)$ first:

This midpoint velocity $\vec{v_\mathrm{mid}}(t + \Delta t)$ is then used for calculating the new position $\vec{p}(t + \Delta t)$ comparable to the previous Euler methods:

Afterward, we're doing something similar for integrating the velocity component: we're once again advancing it comparable to the semi-explicit/explicit Euler variants, but this time we're evaluating the acceleration at the time position $t + \frac{\Delta t}{2}$:

The implementation is once again straightforward:

// Returns the acceleration at the time position t
float acceleration_func(float t);

template<typename T>
static void integrate_euler_modified(T& p, T& v, T a, T t, float dt)
{
  const T vm = v + 0.5 * a * dt;
  p = p + vm * dt;
  v = v + acceleration_func(t + 0.5) * dt;
}

Given a more complex environment, where, e.g., the acceleration term might be based on a vector field and thus the position, it's important to directly use the up-to-date position $\vec{p}(t + \Delta t)$ for integrating the velocity, yielding:

Let's take a look at the results when applied to our particle example with constant acceleration:

These are the best results we've seen up to this point. Even for the largest timesteps, the position is close to spot-on at any given point in time. Checking the error, we end up with a position error of $\approx 4.41e^{-29}\textrm{ m}$ after simulating $0.5\textrm{ s}$ and $\approx 2.10e^{-22}\textrm{ m}$ after $10\textrm{ s}$ with our biggest timestep of $\frac{1}{20}\textrm{ s}$.

For this simple example, the error is almost negligible. But what about the harmonic oscillator? Let's check:

The results are almost comparable to applying the explicit Euler variant. But instead of gaining energy, we're losing energy this time. This drift is, once again, caused by the fact that this integrator is not symplectic.

In the following sections, we're going to look at two second-order integrators that are also symplectic.

Verlet Integration (Second Order, Symplectic)

The standard Verlet integration scheme does not keep track of the previous timestep's velocity. Instead, we're keeping a record of the last position and use this information to deduce the implied velocity and update the position accordingly:

Since it is necessary to have the velocity at hand in many situations, it's possible to derive the velocity from the integrated position and the previous one, of course:

Calculating velocity this way, we end up with an error order of $O(\Delta t)$. It's possible to enlargen the interval used for the calculation to advance to $O({\Delta t}^2)$, but in this case, the results for the velocity term will lag behind:

Once again, a simple implementation in C++ can look as follows:

template<typename T>
void integrate_verlet(T& p, T& pPrev, T& v, T a, float dt)
{
  const T pNew = pos * 2.0 - pPrev + a * dt * dt;
  v = (pNew - pos) / dt;
  // Swap with the above. More precise but lags
  // v = (pNew - pPrev) / (2.0 * dt);

  pPrev = p;
  p = pNew;
}

It's important to note that this scheme is not self-starting. This time we need two timesteps worth of positional information upfront: $\vec{p}(t - \Delta t)$ and $\vec{p}(t)$. Accordingly, it's necessary to utilize another integrator to make an educated guess for $\vec{p}(t)$ - if the information is not easy to deduce otherwise. Simply using the starting position twice will influence the quality of the simulation negatively.

For our constant acceleration example, we end up with results in the same error order in the position term as when we've applied the second-order Midpoint Method:

But this time, since the integrator is also symplectic, we're also achieving more reliable results when applying it to the harmonic oscillator:

While it's evident that the Verlet integrator is very lightweight and accurate, there are a couple of downsides:

• Velocity is not deducible for $t + \Delta t$ in the same error order as the position.
• In specific scenarios, like, e.g., in the context of collision reactions, it's necessary to patch the past positions to incorporate the new velocity.
• Since the integrator is not self-starting, it's necessary to utilize another integrator to get going.
• The integrator can't be applied when the delta-time is varying or fluctuating. There are options to correct this behavior, which are known as Time-Corrected Verlet Integration.

In the following section, we will look at another variant of the Verlet integration scheme that's addressing these downsides.

Velocity Verlet Integration (Second Order, Symplectic)

Compared to the basic Verlet integrator, the Velocity Verlet Integration incorporates velocity for deducing the position term. In addition, if we're not dealing with a simple constant acceleration, it's also necessary to keep track of the acceleration. This fact makes it the most memory-hungry integrator up to this point.

The velocity term can be simplified for the constant acceleration case $\vec{a}(t + \Delta t) = \vec{a}(t)$:

Implementation wise, this is a possible option in C++:

// Returns the acceleration at the time position t
float acceleration_func(float t);

template<typename T>
static void integrate_verlet_velocity(T& p, T& v, T& a, float t, float dt)
{
  p = p + v * dt + 0.5 * a * dt * dt;
  // Remember the "previous acceleration"
  const T a0 = a;
  a = acceleration_func(t + dt);
  // Calculate the new velocity based on the current and previous acceleration
  v = v + (a + a0) * 0.5 * dt;
  // Swap with the above in case the acceleration term is constant
  //v = v + a * dt;
}

Let's take another look at our simple particle example with constant acceleration:

As expected, we're once again in the same ballpark as with all our second-order integrators. Let's quickly check the harmonic oscillator:

Those are the most precise results up to this point for a timestep as big as $\frac{1}{4}\textrm{ s}$. The results clearly show that the precision of the velocity term, compared to the basic Verlet integration scheme, has improved, advancing the global error of the velocity to $O({\Delta t}^2)$.

For the price of storing the acceleration, in addition to the position and velocity, the Velocity Verlet Integration eliminates the major downsides of the basic Verlet Integration:

• We can calculate the velocity in time and the same global error order $O({\Delta t}^2)$ as the position.
• The integrator is self-starting.
• We can directly manipulate the velocity when it's necessary to do so.
• Works correctly when the delta-time is fluctuating.

Where Every Numerical Integrator Breaks

Up to this point, we've only looked at some elementary examples. Sadly, it's rarely possible to account for all the different situations when incorporating, e.g., a particle simulation into a game. As soon as the system gains in complexity, it's almost impossible to ensure that the math is spot-on in all the different configurations. It's even more complicated if the simulation reacts to the player's actions. This makes it practically impossible to ensure that a simulation always turns out the same if we're utilizing variable timesteps.

Let's assume we're once again dealing with our single-particle example. This time we're skipping the gravitational acceleration making the particle rest at its initial position by default. We then go ahead and apply an intense global wind force to the particle:

$f_{\vec{a}}(t) = \sin(x t 2\pi) 150$

Setting the frequency $x$ to $x = 4$, making it change rather slowly over time, and using the Velocity Verlet integrator, we still end up with relatively acceptable results at all the different timesteps:

But imagine the wind frequency is a property exposed in the context of the particle system as a user-modifiable parameter. So it might very well be that someone goes ahead and sets $x = 16$. Since the results look pleasing, the user goes ahead and commits the results. Taking another look at the current state of your game, you're baffled since the effect behaves almost completely broken on your setup. But why? Let's take a look at this concrete example:

It's visible that the simulation degenerates with an increasing timestep size reaching a complete drift at $\Delta t = \frac{1}{20}\textrm{ s}$. But why is that?

It's not possible to accurately simulate any given system with any timestep. In this concrete example, we're violating the sampling theorem, which states that the sampling frequency must be two times the highest frequency we're working with. Accordingly, we're undersampling our wind acceleration function, which shows as phase shifts or complete chaos. Depending on the integrator used, it might even be necessary to utilize an even higher sampling frequency.

For the $\Delta t = \frac{1}{20}\textrm{ s}$ case, we're almost always sampling the acceleration when it's accelerating downwards, explaining the extreme drift in the result.

Based on these results, one can imagine the issues that will arise if we incorporate vector fields, user-induced force fields, and various other modules into our simulation.

Conclusion

As can be seen, variable timesteps can quickly cause issues, especially if the framerate on the end user side is either low or inconsistent. This fact is crucial in cases where a deterministic result is desired, e.g., if particle effects should look or a gameplay system should behave the same, independent of how powerful the user's computer configuration is.

You can quickly achieve deterministic results by choosing a reasonably small but fixed timestep and sub-stepping critical systems in case the overall framerate does not suffice or skipping updates if the outer game loop is overperforming. The following code example demonstrates the general idea:

// Result of our sample task
struct TaskResult;
// Linearly interpolates between two results
TaskResult lerp(const TaskResult& x, const TaskResult& y, float a);
// Stub task function
TaskResult my_precious_task(float dt);

static TaskResult step_fixed(float dt)
{
  // 60 steps per second
  constexpr float timestep = 0.016f;
  // The "current" and "previous" result of the task
  TaskResult result;
  static TaskResult resultPrev;

  // Accumulator which is used to calculate the amount
  // of steps we'll have to take and the fraction
  // for interpolation
  static float acc = 0.0f;
  // Accumulate the time passed since the last call of this function
  acc += dt;
  
  // Skips or sub-steps the task depending on the variable
  // delta time of the outer game loop
  while (acc >= timestep)
  {
    // Execute a single step of this task
    resultPrev = result;
    result = my_precious_task(timestep);

    // Update the accumulator
    acc -= timestep;
  }

  // The results of the task need some interpolation here, especially
  // if they play a visual role, like so:
  return lerp(resultPrev, result, acc / timestep);
}

As you can see, there is a need for extra care to achieve perfectly smooth results, addressed in the example by introducing interpolation. We'll discuss this in more detail in some future blog posts.

Further Reading

• https://gafferongames.com/post/fix_your_timestep/
• https://box2d.org/files/ErinCatto_NumericalMethods_GDC2015.pdf