Euler vs Runge–Kutta: Time Stepping in Simulations

A test problem with a known answer

The frictionless spring x″ = −x is written as a first-order system for position and velocity:

Starting from x = 1, v = 0 the exact solution is x(t) = cos t, and the energy E = (x² + v²)/2 stays exactly ½ forever. Any drift in E is error introduced by the numerical method.

Three stepping rules

• Euler (forward Euler): use the slope at the start of the step for the whole step. x ← x + h·v, v ← v − h·x (both from the old values). One derivative evaluation per step; error per unit time proportional to h.
• Semi-implicit (symplectic) Euler: update the velocity first, then move with the new velocity. Same cost as Euler, same first-order accuracy, but it preserves the geometric structure of mechanical systems, so energy oscillates within a narrow band instead of drifting.
• Classical Runge–Kutta (RK4): evaluate the slope four times — at the start, twice at the midpoint, at the end — and take a weighted average. Error per unit time proportional to h⁴.

const f = ([x, v]) => [v, -x];

const euler = (y, h) => { const k = f(y); return [y[0] + h * k[0], y[1] + h * k[1]]; };

const semiImplicit = ([x, v], h) => { const v2 = v - h * x; return [x + h * v2, v2]; };

function rk4(y, h) {
  const k1 = f(y);
  const k2 = f([y[0] + h / 2 * k1[0], y[1] + h / 2 * k1[1]]);
  const k3 = f([y[0] + h / 2 * k2[0], y[1] + h / 2 * k2[1]]);
  const k4 = f([y[0] + h * k3[0], y[1] + h * k3[1]]);
  return [y[0] + h / 6 * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]),
          y[1] + h / 6 * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1])];
}

def rk4(f, y, h):
    k1 = f(y)
    k2 = f([a + h / 2 * b for a, b in zip(y, k1)])
    k3 = f([a + h / 2 * b for a, b in zip(y, k2)])
    k4 = f([a + h * b for a, b in zip(y, k3)])
    return [a + h / 6 * (p + 2 * q + 2 * r + s)
            for a, p, q, r, s in zip(y, k1, k2, k3, k4)]

Results

Three things stand out.

• Euler pumps energy into the spring. Each step multiplies the energy by exactly 1 + h². Over 100 periods at h = 0.01 that is (1.0001)^62,832 ≈ e^6.28 ≈ 535. The simulated spring spirals outward — a qualitative error no amount of plotting will reveal as "slightly inaccurate".
• Semi-implicit Euler is no more accurate per step, but it is faithful. Its energy stays within a band of width proportional to h and never runs away. For long simulations of mechanical systems that property matters more than the order of accuracy.
• RK4 is dramatically more accurate. Reducing h from 0.1 to 0.05 cut its error from 3.9 × 10⁻⁶ to 2.7 × 10⁻⁷, a factor of about 15 — close to the 2⁴ = 16 predicted for a fourth-order method. It still loses a little energy over very long runs, because it is not symplectic.

Choosing a method for a simulation

• Smooth, non-stiff models integrated over moderate times — the SIR model, logistic growth, Lotka–Volterra, projectiles with drag — are well served by RK4 with a fixed step. The OpenSim.zip models use it and report a conserved quantity or exact solution as a check.
• Long-running mechanical systems where energy conservation matters more than phase accuracy — orbits, molecular dynamics — favor symplectic methods such as semi-implicit Euler or velocity Verlet.
• Stiff systems (rates that differ by orders of magnitude, such as chemical kinetics) need implicit methods; explicit methods including RK4 require impractically small steps to stay stable. The same stability limit shows up in the diffusion simulation, where the explicit scheme is only stable for DΔt/Δx² ≤ ½.
• For serious work, adaptive-step solvers (such as the Dormand–Prince pair behind scipy.integrate.solve_ivp’s default RK45) choose the step automatically from an error estimate.

Always test the step: halve it and rerun. If the answer changes in a way that matters, the step was too large.