Comparison of Simulation Methods — Integrators, Fluid Models & Monte Carlo

1. Integrator Comparison Table

The integrator is the piece of code that advances position and velocity by a timestep dt. All of the differences below come down to three properties: order of accuracy (how fast the local error shrinks as dt shrinks), energy behaviour (does the simulation gain or lose energy over long runs), and cost (force evaluations per step).

"Symplectic" methods exactly conserve a nearby (shadow) Hamiltonian, so energy error stays bounded forever instead of drifting away — critical for long-running simulations like orbital mechanics or particle systems that must not visibly gain energy.

2. Explicit Euler

The simplest possible integrator: evaluate the derivative at the current state, then step forward using that single slope estimate.

// Explicit (forward) Euler — DO NOT use for anything oscillatory
function stepEuler(body, dt) {
  const a = body.force.scale(body.invMass);
  body.pos = body.pos.add(body.vel.scale(dt)); // uses OLD velocity
  body.vel = body.vel.add(a.scale(dt));
}

Because position uses the velocity from before the force was applied, explicit Euler systematically adds energy to any oscillating system (springs, orbits, pendulums). A frictionless spring simulated with explicit Euler will visibly spiral outward and eventually explode, no matter how small dt is made — the error only shrinks, it never disappears qualitatively. It is included here purely as the baseline every other method improves on.

3. Semi-Implicit (Symplectic) Euler

A one-line fix: update velocity first, then use the new velocity to update position. Same cost as explicit Euler, dramatically better energy behaviour.

function stepSemiImplicit(body, dt) {
  const a = body.force.scale(body.invMass);
  body.vel = body.vel.add(a.scale(dt));      // velocity first
  body.pos = body.pos.add(body.vel.scale(dt)); // then position, using new v
}

4. Verlet & Velocity Verlet

Position (Störmer) Verlet drops velocity entirely and derives it implicitly from two stored positions. It is second- order accurate for the same single force evaluation as Euler, which makes it the default choice for cloth, ropes, and particle chains.

function stepVerlet(p, dt) {
  const a = p.force.scale(p.invMass);
  const next = p.pos
    .scale(2)
    .sub(p.prevPos)
    .add(a.scale(dt * dt));
  p.prevPos = p.pos;
  p.pos = next;
}

Velocity Verlet keeps an explicit velocity, giving the same second-order accuracy but with a usable v(t) for damping, drag, or velocity-dependent forces:

5. Runge-Kutta 4 (RK4)

RK4 samples the derivative four times per step — at the start, twice at the midpoint, and once at the end — and combines them into a weighted average. It is fourth-order accurate: halving dt shrinks the local error by roughly 16×, versus 4× for Verlet and 2× for Euler.

function rk4Step(state, dt, derivative) {
  const k1 = derivative(state);
  const k2 = derivative(addScaled(state, k1, dt / 2));
  const k3 = derivative(addScaled(state, k2, dt / 2));
  const k4 = derivative(addScaled(state, k3, dt));
  return addScaled(
    state,
    combine(k1, k2, k3, k4), // k1 + 2k2 + 2k3 + k4
    dt / 6
  );
}

This site's Lorenz Attractor and Double Pendulum simulations use RK4 because chaotic systems amplify small numerical errors exponentially — a low-order integrator's error becomes visually wrong trajectories within seconds. RK4's downside: it is not symplectic, so orbital simulations run with RK4 slowly gain or lose energy over thousands of orbits (use Leapfrog or Verlet instead for N-body work that must run indefinitely).

6. Implicit Euler & Stiff Systems

Implicit (backward) Euler evaluates the force at the future state instead of the current one, which requires solving an equation (often a sparse linear system) each step rather than just plugging in numbers.

This extra work buys unconditional stability: even with a huge dt, an implicit integrator will not blow up — it just damps out high-frequency motion, which is usually a feature. That property matters most for stiff systems, where forces vary over very different time scales (e.g. a very taut cloth spring next to a loose one). Baraff and Witkin's 1998 "Large Steps in Cloth Simulation" paper popularised implicit integration for cloth precisely because explicit methods needed impractically tiny timesteps to stay stable on stiff springs.

7. Lagrangian vs. Eulerian Fluid Simulation

Beyond single-particle integrators, continuous media (fluids, smoke, terrain erosion) can be modelled from two opposite viewpoints.

Lagrangian: follow the material

Represent the fluid as a set of moving particles, each carrying mass, velocity, and density. Smoothed Particle Hydrodynamics (SPH), used in this site's SPH Fluid simulation, is the classic example: each particle's properties are computed by smoothing (kernel-weighting) over nearby particles.

SPH naturally conserves mass (particles ARE the mass), handles free surfaces and splashing without extra bookkeeping, and its cost scales with the number of particles rather than the size of empty space. Its downside is the O(N·k) neighbour search every step and noisier surfaces at low particle counts.

Eulerian: fix a grid, let material flow through it

Discretise space into a fixed grid and store velocity/density/ pressure at each cell, solving the Navier-Stokes equations with finite differences. Jos Stam's "Stable Fluids" method (semi-Lagrangian advection + pressure projection) is the standard real-time technique:

• Advection: for each cell, trace backward along the velocity field and sample the previous grid — unconditionally stable regardless of timestep, unlike naive forward advection.
• Pressure projection: solve a Poisson equation (Jacobi or conjugate-gradient iteration) to remove divergence and enforce incompressibility.

Grid methods give smooth, artifact-free fields and are ideal for smoke, gases, and reaction-diffusion patterns (see this site's Reaction-Diffusion simulation, which uses a fixed-grid ping-pong GPU texture). Their weakness is resolution-bound cost — a 256×256 grid spends equal compute on empty air and the interesting fluid region — and grid methods lose fine detail (small droplets, thin sheets) below one cell's size.

8. Monte Carlo Methods

Not every simulation advances a deterministic state through time. Monte Carlo methods use repeated random sampling to estimate a result that would otherwise require an intractable integral or exhaustive search — trading determinism for a statistical answer that converges as sample count grows.

That 1/√N convergence rate is both Monte Carlo's biggest strength and its biggest weakness: it is the same rate no matter how many dimensions the problem has, which is why Monte Carlo beats grid-based quadrature for high-dimensional integrals (path tracing light transport, financial option pricing), but it converges slowly compared to a good deterministic method when the dimensionality is low.

• Direct sampling: generate uniform random points and average — used by this site's Genetic Algorithm simulation for mutation and crossover selection.
• Markov Chain Monte Carlo (MCMC): sample from complex probability distributions by taking a random walk that preferentially visits high-probability regions (Metropolis- Hastings). Used for statistical mechanics and Bayesian inference.
• Monte Carlo path tracing: in rendering, randomly sample light paths bouncing through a scene to solve the rendering equation — the noisy-then-converges look of path-traced renders is exactly this 1/√N convergence made visible.

9. How to Choose

• Real-time game/particle physics: semi-implicit Euler, optionally with 2–4 substeps per frame for extra stability — cheapest option with acceptable energy behaviour.
• Cloth, ropes, soft bodies: position Verlet with constraint relaxation, or implicit/PBD solvers if the material is stiff and must not blow up at large timesteps.
• Orbits, N-body, anything that must run indefinitely without gaining energy: Leapfrog or velocity Verlet — symplectic integrators are the correct default here, not RK4.
• Chaotic systems where trajectory accuracy matters more than long-run energy (Lorenz, double pendulum): RK4, or adaptive-step RK45 if you need error control.
• Splashing liquids, free surfaces: SPH (Lagrangian).
• Smoke, gas, reaction-diffusion patterns: grid- based Eulerian solvers (Stable Fluids, ping-pong shaders).
• High-dimensional estimation, rendering, optimisation search: Monte Carlo / MCMC / genetic algorithms.