explainer
Numerical integration: compare drift, phase, and step cost
Advance an oscillator with forward Euler, velocity-first symplectic Euler, and classical RK4. Compare each method with the analytic solution and separate energy drift, phase error, step cost, and stability.
What you will learn
- Apply three integration methods to the same position-velocity state.
- Explain why updating velocity first changes the Euler method's behavior.
- Compare state error, energy error, and phase error against a known solution.
- Relate step size to convergence, computational work, and oscillator stability.
Before you start
- ODEs, initial conditions, and forward Euler steps
- Derivatives and local rate predictions
- Sine, cosine, and radians
An undamped oscillator should keep the same energy. Forward Euler adds energy every time it takes a positive step. Updating velocity before position changes that behavior, while classical RK4 spends more work per step to follow the solution more closely.
Numerical integration of an ODE builds a sequence of approximate states from a rate law and an initial state. Here, a known exact solution lets us measure what each method gets wrong.
Start with an oscillator whose solution is known
Use a normalized position q and velocity v:
q′ = v, v′ = −q
(q(0), v(0)) = (1, 0)
This is equivalent to q″ = −q. Its solution is q(t) = cos t and v(t) = −sin t. Differentiating verifies both equations. In the (q, v) plane, the state travels clockwise around the unit circle.
The dimensionless energy is E = ½(q² + v²). Along the exact solution, E′ = qq′ + vv′ = qv − vq = 0, so E remains 0.5.
Normalization matters. For a physical spring with mass m, stiffness k, displacement scale L, and time τ, set ω = √(k/m), t = ωτ, q = x/L, and v = ẋ/(Lω). Then q and v are dimensionless, and E is physical energy divided by kL². A dimensionless step h corresponds to physical duration h/ω.
This model can describe a simplified compliant robot joint near equilibrium, before adding damping, forcing, or contact. Simulation software advances positions and velocities from computed accelerations; MuJoCo documents this integration step explicitly. In machine learning, neural ODEs use a learned rate function and an ODE solver to compute a hidden state's evolution. Our oscillator is a controlled example of solver error, with no learned parameters or measured robot data.
Advance both variables from the old state
With a positive step h, forward Euler uses the rates at the beginning of the interval:
qₙ₊₁ = qₙ + hvₙ
vₙ₊₁ = vₙ − hqₙ
Keep the old q until both expressions are evaluated. Reusing an already updated variable would implement another method.
Substitute these updates into the energy formula. The cross terms cancel, leaving Eₙ₊₁ = (1 + h²)Eₙ. For any nonzero state and any h greater than zero, energy increases. After n steps, Eₙ = (1 + h²)ⁿE₀.
Smaller steps reduce this error over a fixed time interval. They do not create a positive step size for which forward Euler keeps this oscillator bounded forever. Convergence as h approaches zero and behavior over arbitrarily long time are different questions.
Update velocity before position
The velocity-first semi-implicit Euler method changes the order of operations:
vₙ₊₁ = vₙ − hqₙ
qₙ₊₁ = qₙ + hvₙ₊₁
For this oscillator, it is a symplectic Euler method, also called Euler-Cromer. Position uses the new velocity. Both calculations are explicit here because acceleration depends only on position. More general symplectic Euler formulas can require an implicit solve. Hairer's lecture notes state this distinction and identify the method as first order.
Symplectic methods preserve a geometric structure of Hamiltonian motion. In this two-dimensional example, the update preserves area in the (q, v) plane. It does not preserve the exact energy E at each step.
For a fixed h, direct substitution shows that this particular update preserves ½(q² + v² − hqv), apart from floating-point error. That modified quadratic differs from physical energy. For 0 < h < 2, its level curves are ellipses that keep q and v bounded. The step restriction matters.
Sample four slopes with classical RK4
Let z = (q, v) and F(z) = (v, −q). Classical fourth-order Runge-Kutta, or RK4, builds four rate estimates:
k₁ = F(z)
k₂ = F(z + ½hk₁)
k₃ = F(z + ½hk₂)
k₄ = F(z + hk₃)
zₙ₊₁ = z + ⅙h(k₁ + 2k₂ + 2k₃ + k₄)
Here each k is a rate, so h appears in the final weighted sum. For a time-dependent F(t, z), evaluate the four rates at times t, t + h/2, t + h/2, and t + h. These intermediate states are predictions used inside one step.
Driscoll and Braun give the classical RK4 formula and its four-stage construction. Higher order improves small-step accuracy under suitable smoothness conditions. It does not make the method exact or remove every stability restriction.
Calculate the first half-unit step
Start at (1, 0), choose h = 0.5, and compare the first update.
Forward Euler gives q₁ = 1 + 0.5(0) = 1 and v₁ = 0 − 0.5(1) = −0.5. Symplectic Euler first gets the same v₁, then uses it to obtain q₁ = 1 + 0.5(−0.5) = 0.75.
For RK4, the four rates are:
| Rate | Value |
|---|---|
| k₁ | (0, −1) |
| k₂ | (−0.25, −1) |
| k₃ | (−0.25, −0.9375) |
| k₄ | (−0.46875, −0.875) |
The weighted update gives (q₁, v₁) = (337/384, −23/48), about (0.877604, −0.479167). The exact state is (cos 0.5, −sin 0.5), about (0.877583, −0.479426).
| Method | Energy after one step |
|---|---|
| Exact solution | 0.500000 |
| Forward Euler | 0.625000 |
| Symplectic Euler | 0.406250 |
| Classical RK4 | 0.499895 |
The symplectic result already demonstrates that preserving geometric structure does not mean exact energy preservation. RK4 has a smaller error here and uses four acceleration evaluations. Each Euler method uses one.
Compare all methods over the same interval
Every run below starts at (1, 0) and ends at t = 10. The selected method appears in both diagrams; the table compares all three methods at the same h.
The default forward Euler run uses h = 0.25 and 40 steps. Its final energy is about 5.651029, despite the unchanged exact value 0.5. Choose Small steps to reduce h to 0.1. Error decreases over the same interval, while the number of steps rises to 100.
Choose Compare at h = 0.5, then switch the plotted method. The symplectic path stays bounded at this step size but misses the exact endpoint. RK4 follows more closely. Inspect the final coordinates and comparison values because each plot rescales to include its selected trajectory.
The amber segments simply connect numerical samples. Coarse steps can skip much of the motion between samples, so a tidy-looking segment does not prove an accurate trajectory.
Separate energy, phase, and stability
Energy error describes a wrong distance from the origin in this normalized phase plane. Phase error describes being at the wrong point around the cycle. Two states can share exactly the same energy and still be far apart. For example, (1, 0) and (0, −1) both have energy 0.5 but differ by a quarter-cycle.
The readout uses the phase angle atan2(−v, q), then compares the numerical and exact final angles. It reports the shortest signed difference between −180° and 180°. A positive result means ahead of the exact phase. Wrapping discards full revolutions, so this value alone cannot detect an entire missed cycle.
For symplectic Euler, eliminating v from the updates gives a position recurrence with characteristic equation λ² − (2 − h²)λ + 1 = 0. When 0 < h < 2, its two distinct roots lie on the unit circle. At h = 2 they merge, and this initial state grows: the first position-velocity pairs are (−3, −2), (5, 4), and (−7, −6). The Stability boundary preset reaches (−11, −10) at t = 10. Steps above 2 also produce growth for this initial state.
Classical RK4 has its own restriction. Substitution into this oscillator gives the per-step energy multiplier 1 − h⁶/72 + h⁸/576. It is below 1 for 0 < h < 2√2, causing artificial damping, and exceeds 1 when h > 2√2. At the endpoint 2√2 it equals 1, but phase error remains. The experiment's largest h is 2.5, inside that RK4 interval and outside the symplectic interval.
These boundaries belong to these methods on this particular oscillator. The ODE stability lesson develops the distinction between a model's dynamics and its numerical approximation.
Compare accuracy at a stated computational cost
For smooth problems over a fixed finite interval, the two Euler methods have global error O(h), while classical RK4 has global error O(h⁴) as h approaches zero. Starting a step from the exact state, their one-step errors are O(h²) and O(h⁵), respectively. MIT's numerical ODE guide distinguishes local and accumulated error orders.
The Taylor expansion lesson explains how omitted higher-order terms determine local error.
In the small-step regime, halving h typically reduces the leading global error by about two for a first-order method and sixteen for a fourth-order method. Roundoff, vanishing leading error terms, or steps outside that regime can change the observed ratio. A single coarse comparison cannot establish an order.
Each Euler step here evaluates acceleration once; each RK4 step evaluates it four times. At h = 0.25, that means 40 versus 160 acceleration evaluations. The experiment keeps the same step size across methods, so its accuracy comparison uses different computational budgets.
These counts measure one operation. Runtime also includes other work and depends on the implementation.
A useful simulation check repeats a run with smaller steps and inspects the quantities that matter: timing, state error, conserved quantities, or task outcomes. Real robot models may include contact and fast damping that require different methods. Accuracy on this smooth oscillator establishes neither a universal solver choice nor an accurate physical model.
Reproduce the experiment in Python
This standard-library program implements all three updates directly. Each method starts afresh, and the integer loop count makes every final comparison occur at t = 10.
from math import cos, sin, hypot
def advance(method, state, h):
q, v = state
if method == "Euler":
return q + h * v, v - h * q
if method == "Symplectic":
new_v = v - h * q
return q + h * new_v, new_v
if method != "RK4":
raise ValueError("Unknown method")
k1 = (v, -q)
k2 = (v + h * k1[1] / 2, -(q + h * k1[0] / 2))
k3 = (v + h * k2[1] / 2, -(q + h * k2[0] / 2))
k4 = (v + h * k3[1], -(q + h * k3[0]))
return tuple(state[i] + h * (
k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]
) / 6 for i in range(2))
def energy(state):
q, v = state
return (q * q + v * v) / 2
for method in ("Euler", "Symplectic", "RK4"):
first = advance(method, (1.0, 0.0), 0.5)
print(f"{method} first=({first[0]:.6f}, {first[1]:.6f}), "
f"energy={energy(first):.6f}")
exact = (cos(10), -sin(10))
print(f"Exact at 10=({exact[0]:.6f}, {exact[1]:.6f})")
for method in ("Euler", "Symplectic", "RK4"):
state = (1.0, 0.0)
for _ in range(40):
state = advance(method, state, 0.25)
error = hypot(state[0] - exact[0], state[1] - exact[1])
print(f"{method} at 10: error={error:.6f}, "
f"energy={energy(state):.6f}")
Expected output:
Euler first=(1.000000, -0.500000), energy=0.625000
Symplectic first=(0.750000, -0.500000), energy=0.406250
RK4 first=(0.877604, -0.479167), energy=0.499895
Exact at 10=(-0.839072, 0.544021)
Euler at 10: error=2.390297, energy=5.651029
Symplectic at 10: error=0.089779, energy=0.446303
RK4 at 10: error=0.000325, energy=0.499933
Try it yourself
Exercise 1. Starting from (q, v) = (0, 1), take one step with h = 0.5 using forward Euler and velocity-first symplectic Euler. Calculate both energies. Does symplectic Euler always decrease the physical energy?
Show solution 1
Forward Euler gives (q₁, v₁) = (0.5, 1). Symplectic Euler first obtains v₁ = 1 − 0.5(0) = 1, then q₁ = 0 + 0.5(1) = 0.5. Both energies are ½(0.25 + 1) = 0.625, above the initial 0.5.
Symplectic Euler can increase or decrease physical energy on a step. Its preserved modified quadratic here is 0.625 − ½(0.5)(0.5)(1) = 0.5. The earlier worked example decreased physical energy because it started at another point in the cycle.
Exercise 2. For a fixed interval of length 10, forward Euler uses h = 0.1 and RK4 uses h = 0.4. How many acceleration evaluations does each use? If a numerical endpoint has the correct energy, is its state necessarily accurate?
Show solution 2
Forward Euler takes 100 steps at one evaluation each. RK4 takes 25 steps at four evaluations each. Both use 100 acceleration evaluations. Equal counts make this a useful comparison of this operation's cost; they do not guarantee equal wall-clock time or accuracy.
Correct energy fixes the radius in the normalized phase plane, leaving the phase unconstrained. A numerical result at (0, −1) when the exact state is (1, 0) has the correct energy 0.5 but state error √2. Check the full state and timing as well as energy.
Sources and further study
- Ernst Hairer, Geometric Numerical Integration, Lecture 2: Symplectic integrators, for symplectic Euler, its order, and separable Hamiltonian systems.
- Tobin Driscoll and Richard Braun, Fundamentals of Numerical Computation: Runge-Kutta methods, for the classical four-stage method and comparisons using function-evaluation cost.
- MIT, Numerical Solutions for Differential Equations, for local and accumulated error and method order.
- MuJoCo computation documentation: numerical integration, for integrating robot positions and velocities from accelerations.
- Chen, Rubanova, Bettencourt, and Duvenaud, Neural Ordinary Differential Equations, for a machine-learning model whose output requires solving an ODE.