Explicit Euler

Deep dive

How it works

In: f, t₀, y₀, step h, step count · Out: y at t₀ + N·h · State: none between steps (y updated in place, one scratch vector) · Cost: 1 f-evaluation + n multiply-adds per step; global error O(h)

1. Evaluate k = f(t_n, y_n).
2. y_{n+1} = y_n + h·k; t_{n+1} = t_n + h.
3. Error: Taylor gives y(t + h) = y + h·y' + h²/2·y''(ξ), so the local error is h²/2·|y''|; over T/h steps it accumulates to O(h).
4. Stability: applied to y' = λy, the method multiplies y by R(z) = 1 + z, z = hλ, each step. Decay requires |R(z)| ≤ 1: a disc of radius 1 centred at −1. For a system, every eigenvalue of the Jacobian, times h, must lie in that disc.

 Im z
  |      ___
  |    /     \        stable region of explicit Euler:
--+---|---*---|---    disc |1 + z| ≤ 1, centre -1
 -2    \ -1  /  0     the imaginary axis is outside (except 0):
  |      ‾‾‾          oscillators always gain energy

Code

/// One step: y <- y + h f(t, y). `k` is scratch of the same length as y.
pub fn euler_step<F: FnMut(f64, &[f64], &mut [f64])>(f: &mut F, t: f64, y: &mut [f64], h: f64, k: &mut [f64]) {
    f(t, y, k);
    for (yi, ki) in y.iter_mut().zip(k.iter()) {
        *yi += h * ki;
    }
}

/// Integrate from t0 to t1 in n equal steps.
pub fn integrate<F: FnMut(f64, &[f64], &mut [f64])>(mut f: F, t0: f64, y0: &[f64], t1: f64, n: usize) -> Vec<f64> {
    let h = (t1 - t0) / n as f64;
    let mut y = y0.to_vec();
    let mut k = vec![0.0; y.len()];
    for i in 0..n {
        euler_step(&mut f, t0 + i as f64 * h, &mut y, h, &mut k); // t from i, not accumulated: no drift
    }
    y
}

/// Estimate the Jacobian's spectral radius at (t, y) by power iteration on
/// finite-difference products J v ≈ (f(y + εv) − f(y)) / ε. Then h < 2 / rho for real negative spectra.
pub fn spectral_radius<F: FnMut(f64, &[f64], &mut [f64])>(f: &mut F, t: f64, y: &[f64], iters: usize) -> f64 {
    let n = y.len();
    let (mut f0, mut f1) = (vec![0.0; n], vec![0.0; n]);
    f(t, y, &mut f0);
    let mut v: Vec<f64> = (0..n).map(|i| 1.0 + 0.1 * i as f64).collect(); // generic start vector
    let mut rho = 0.0;
    let ynorm = y.iter().map(|x| x * x).sum::<f64>().sqrt();
    for _ in 0..iters {
        let vnorm = v.iter().map(|x| x * x).sum::<f64>().sqrt();
        let eps = f64::EPSILON.sqrt() * (1.0 + ynorm) / vnorm;
        let yp: Vec<f64> = y.iter().zip(&v).map(|(a, b)| a + eps * b).collect();
        f(t, &yp, &mut f1);
        for i in 0..n { v[i] = (f1[i] - f0[i]) / eps; } // v <- J v
        let norm = v.iter().map(|x| x * x).sum::<f64>().sqrt();
        rho = norm / vnorm;
        if norm == 0.0 { break; }
    }
    rho
}

/// Euler-Maruyama for dX = a(X) dt + b(X) dW (scalar noise per component).
/// `normal` must return independent N(0, 1) samples. Strong order 0.5, weak order 1.
pub fn euler_maruyama_step(
    a: impl Fn(&[f64], &mut [f64]), b: impl Fn(&[f64], &mut [f64]),
    x: &mut [f64], dt: f64, normal: &mut impl FnMut() -> f64, scratch: &mut [f64],
) {
    let n = x.len();
    let (da, db) = scratch.split_at_mut(n);
    a(x, da);
    b(x, db);
    let sq = dt.sqrt();
    for i in 0..n { x[i] += da[i] * dt + db[i] * sq * normal(); }
}

Complexity and constants

• Cost per step: one f-evaluation plus n fused multiply-adds. For cheap f (a few flops per component) the update is memory-bound: order of 16-24 bytes moved per component per step.
• Error example: y' = −y on [0, 1]: n = 100 gives error ≈ 1.8·10⁻³, n = 1000 gives ≈ 1.8·10⁻⁴. RK4 with n = 10 already reaches order of 10⁻⁷.
• Rounding: each step adds order of ε·|y| (ε = 2.2·10⁻¹⁶), so with N = T/h steps total error ≈ C·h + N·ε·|y|. Below h of order 10⁻⁸·T the error grows again.
• Stability limits: λ = −1000 forces h ≤ 0.002 regardless of accuracy needs (stiffness). Diffusion on a grid: h ≤ Δx²/(2α) (1D), Δx²/(4α) (2D, square grid).

Efficiency

• Cost: one f-evaluation plus n multiply-adds per step; for cheap f the update is memory-bound, order of 16-24 bytes moved per component per step. Accuracy is the expensive part: y' = −y on [0, 1] needs n = 1000 steps for an error of 1.8·10⁻⁴, while RK4 reaches order of 10⁻⁷ with n = 10 (numbers above). Bound: memory bandwidth, then the stability limit h ≤ 2/|λ|.
• Levers:
  • higher-order method instead of a smaller h — 4 evaluations per step but orders of magnitude lower error (example above) — smooth, non-stiff problems
  • h = c/ρ with c < 2 from a power-iteration estimate of the spectral radius — largest stable step, no blow-ups — stability-limited problems
  • implicit Euler or BDF — step set by accuracy instead of stability (Hairer et al. 1993) — stiff systems
  • in-place update with a caller-provided scratch buffer — no allocation per step — small systems, tight loops
  • batch many independent paths or particles in SIMD/GPU lanes — throughput scales with bandwidth — Monte Carlo and large n
• Naive → standard → state of the art:
  • forward Euler with h shrunk until the error is small · n = 1000 for 1.8·10⁻⁴ (above)
  • RK4 at fixed h · n = 10 for order of 10⁻⁷ (above)
  • adaptive embedded RK, or implicit BDF if stiff · error-controlled step size, work follows the solution

Implementation notes

• Compute t as t0 + i·h, not by summing h, to avoid drift in long runs.
• Write f to fill a caller-provided buffer (&mut [f64]); allocation per step dominates for small systems.
• For a sanity check of a new integrator, verify the order: halve h and check the error ratio (2 for Euler, 16 for RK4).
• When stability (not accuracy) limits h, estimate the spectral radius (above) once in a while and set h = c/ρ with c < 2, or switch to an implicit method.
• Gradient descent viewpoint: learning rate η is h, the Hessian's largest eigenvalue L plays |λ|, so η < 2/L is the same stability bound.

Pitfalls

• Using it for springs, orbits or cloth: the energy grows exponentially; the simulation "explodes" for any step size given enough time.
• Variable h chosen without error control: no guarantee on accuracy.
• Euler-Maruyama with dt instead of √dt on the noise term: wrong diffusion scaling.
• Mistaking instability for a bug in f: if the blow-up rate changes with h, it is the integrator.

Variants and alternatives

In production

Finance: Euler-Maruyama in Monte Carlo option pricing (Glasserman). Machine learning: forward-Euler and Euler-Maruyama samplers for diffusion models (score-based SDE sampling), the fixed-grid euler solver in torchdiffeq, and SGD itself as Euler on the gradient flow. Simulation: many simple particle systems (often a mistake for anything spring-like).

Implementations

Sources

• Ernst Hairer, Syvert P. Nørsett, Gerhard Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed., Springer, 1993. book
• Randall J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM, 2007. book
• Peter E. Kloeden, Eckhard Platen, Numerical Solution of Stochastic Differential Equations, Springer, 1992. book
• Paul Glasserman, Monte Carlo Methods in Financial Engineering, Springer, 2004. book
• Yang Song et al., Score-Based Generative Modeling through Stochastic Differential Equations, ICLR 2021. paper