Time integration

Time integration methods:

Simulating the motion of a physical system means solving its equations of motion, a system of ordinary differential equations (ODEs) that relate positions, velocities and accelerations. Since these equations can generally not be solved analytically, we approximate the solution step by step: given the state (position and velocity) at time $t$, a time integration method computes an approximation of the state at time $t + \Delta t$ using the currently acting forces. Repeating this update for a sequence of small time steps $\Delta t$ advances the simulation over time.

This example shows the motion of a single particle that has an initial velocity and is subject to a constant gravitational acceleration. Because the acceleration is constant here, this special case is one of the rare ones that admits a closed-form solution: $$\begin{align*} \mathbf x(t + \Delta t) &= \mathbf x(t) + \Delta t \mathbf v(t) + \frac12 \Delta t^2 \mathbf a^\text{grav} \\ \mathbf v(t + \Delta t) &= \mathbf v(t) + \Delta t \mathbf a^\text{grav} . \end{align*}$$ This exact trajectory is plotted in green and serves as ground truth against which the numerical solutions (blue), computed with the methods introduced below, can be compared. As soon as forces depend on the state itself in a more complex way — e.g. spring, contact or SPH pressure forces — no analytic solution is available anymore and numerical time integration becomes essential.

Three properties are commonly used to characterize and compare time integration methods:

• Order of accuracy: A method of order $p$ has a local error (introduced in a single step) of $O(\Delta t^{p+1})$ and, consequently, a global error (accumulated over a fixed time interval) of $O(\Delta t^{p})$. Halving $\Delta t$ for a method of order $p$ roughly reduces the global error by a factor of $2^p$.
• Computational cost: given by the number of evaluations of the acceleration/force function required per time step. This is important in practice since evaluating the forces (e.g. for cloth, deformable solids or SPH fluids) is usually the most expensive part of a simulation step.
• Stability and energy behavior: describes how the numerical solution behaves over long simulation times, in particular whether it stays close to the true trajectory or whether energy is spuriously gained, lost, or oscillates around the correct value.

Try the different methods in the simulation with a large time step (e.g. $\Delta t = 0.1$) and compare how far the blue (numerical) trajectory deviates from the green (exact) one.

Explicit Euler

The explicit Euler method (also called forward Euler) is the simplest time integration scheme for ODEs. It uses only quantities known at the current time $t$ to extrapolate the state to $t + \Delta t$: $$\begin{align*} \mathbf x(t + \Delta t) &= \mathbf x(t) + \Delta t \mathbf v(t) \\ \mathbf v(t + \Delta t) &= \mathbf v(t) + \Delta t \mathbf a(t). \end{align*}$$

Properties:

• First-order accurate ($O(\Delta t)$ global error); requires only a single force evaluation per step, so it is the cheapest of the three methods.
• Not symplectic: it systematically injects energy into the system because the position update uses the old velocity $\mathbf v(t)$ instead of the newly computed one. For an oscillating system (e.g. a mass on a spring) this causes the amplitude to grow over time, even though the exact solution conserves energy.
• Only conditionally stable: for stiff forces (e.g. stiff springs) $\Delta t$ must be chosen small enough, otherwise the simulation diverges numerically. For this reason explicit Euler is rarely used in practice and mainly serves as a baseline for comparison.

Symplectic Euler

The symplectic Euler method is also known as semi-implicit Euler: it first updates the velocity using the current acceleration, and then uses this already updated velocity to advance the position: $$\begin{align*} \mathbf v(t + \Delta t) &= \mathbf v(t) + \Delta t \mathbf a(t) \\ \mathbf x(t + \Delta t) &= \mathbf x(t) + \Delta t \mathbf v(t + \Delta t). \end{align*}$$

Properties:

• Still only first-order accurate, just like explicit Euler, and has the same cost of a single force evaluation per step — switching from explicit to symplectic Euler is essentially "free" in terms of performance.
• Symplectic: the method (approximately) preserves phase-space volume, which implies that the energy error stays bounded and oscillates around the true value instead of drifting away over time. This yields qualitatively much better long-term behavior than explicit Euler, even though the two methods have the same order and the same error at any single instant.
• Because it is cheap, simple to implement and behaves well over long simulations, symplectic Euler is the de-facto standard in real-time physics: it is used for particle systems, cloth, rigid bodies, SPH fluids and position-based dynamics alike.

Runge-Kutta 2

Higher-order integration methods are usually formulated for a general first-order ODE $\dot{\mathbf s}(t) = \mathbf f(t, \mathbf s(t))$, where $\mathbf s$ is the state of the system at time $t$. In our example the state consists of the particle's position and velocity, and the ODE is defined by: $$\begin{equation*} \mathbf f(t, \mathbf s(t)) = \begin{pmatrix} \dot{\mathbf x} \\ \dot{\mathbf v} \end{pmatrix} = \begin{pmatrix} \mathbf v(t) \\ \mathbf a(t) \end{pmatrix}, \quad\quad \mathbf s(t) = \begin{pmatrix} \mathbf x(t) \\ \mathbf v(t) \end{pmatrix}. \end{equation*}$$ Using this notation, the second-order Runge-Kutta method (also known as the midpoint method) is defined as: $$\begin{align*} \mathbf k_1 &= \Delta t \, \mathbf f(t, \mathbf s(t)) \\ \mathbf k_2 &= \Delta t \, \mathbf f\!\left(t + \tfrac12 \Delta t,\; \mathbf s(t) + \tfrac12 \mathbf k_1\right) \\ \mathbf s(t + \Delta t) &= \mathbf s(t) + \mathbf k_2. \end{align*}$$ Intuitively, $\mathbf k_1$ is a trial (explicit Euler) step used only to estimate the state at the midpoint $t + \frac12 \Delta t$; the actual step $\mathbf k_2$ is then taken using the derivative evaluated at that midpoint, which cancels out the leading-order error term of the Euler methods.

Properties:

• Second-order accurate ($O(\Delta t^2)$ global error): halving $\Delta t$ roughly quarters the error, so for smooth, non-stiff problems RK2 reaches a given accuracy with a much larger time step than either Euler method.
• Has two stages, i.e. the function $\mathbf f$ must be evaluated twice per time step — roughly twice the cost of an Euler step. In general, an order-$p$ Runge-Kutta method typically needs at least $p$ stages, so higher accuracy is traded against more force evaluations per step.
• Not symplectic: despite its higher accuracy, RK2 can still exhibit long-term energy drift on oscillatory problems; the higher order just makes the drift smaller and less visible over short simulation times.

Summary

The following table summarizes the properties discussed above:

Method	Order	Force evals / step	Symplectic	Typical use
Explicit Euler	1	1	No	Baseline / teaching example
Symplectic Euler	1	1	Yes	Default for real-time simulation (games, cloth, rigid bodies, SPH, PBD)
Runge-Kutta 2	2	2	No	Higher short-term accuracy for smooth, non-stiff problems

There is no single "best" method: the right choice depends on whether short-term accuracy, long-term (energy) stability, or per-step cost matters most for the application at hand. Higher-order methods such as Runge-Kutta 4, or implicit methods for stiff problems, extend these same trade-offs further and are covered in later examples.

References

• [HLW06] Ernst Hairer, Christian Lubich, Gerhard Wanner. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, 2nd edition. Springer, 2006.
• [BW98] David Baraff, Andrew Witkin. Large Steps in Cloth Simulation. In Proceedings of SIGGRAPH, 1998.