From derivative to difference
An equation of motion is an ordinary differential equation (ODE): given the state now, it tells you the rate of change now. Newton's second law is second order, m\ddot{x}=F(x,\dot{x},t); rewrite it as two first-order equations for position and velocity and you have the standard form every solver eats.
Any ODE system reduces to this first-order form, where \mathbf{y} stacks all positions and velocities. A solver's whole job is to march \mathbf{y} forward from the initial state \mathbf{y}(t_0).
To march forward we need the derivative as a difference. From the Taylor series come the finite-difference formulas — the crude forward difference and the more accurate central difference.
The central difference cancels the h^2 term of the Taylor expansion, so its error falls as h^2 rather than h — one extra order of accuracy for the same amount of work.
The simplest stepper: Euler's method
The Euler method is the finite-difference idea applied to time: pretend the slope is constant across one small step h, and take it.
Forward (explicit) Euler — the 'just step in the current direction' rule. Transparent, but crude.
Its local error per step is O(h^2), but accumulated over the \sim 1/h steps of a fixed time interval, the global error is only O(h) — first order. Worse, for oscillatory or orbital problems explicit Euler systematically pumps energy in: run a planet with it and the orbit spirals outward, a purely numerical catastrophe with no basis in physics.
The workhorse: Runge-Kutta
The Runge-Kutta family buys accuracy by sampling the slope several times within a step and averaging them cleverly. The classic RK4 evaluates \mathbf{f} four times per step.
RK4: four slope samples — one at the start, two at the midpoint, one at the end — averaged with weights 1:2:2:1. Its global error is O(h^4), making it the default general-purpose integrator.
RK4 is the sensible default for a generic ODE. But it is not special about physics: it does not know that a Hamiltonian system ought to conserve energy and phase-space volume. Over a million orbits, even RK4's energy slowly drifts. For long-term mechanics we want a method built to respect the geometry of the problem.
The physicist's favourite: Verlet and symplectic integrators
For a Newtonian system whose force depends only on position, Verlet integration is beautiful. Add the Taylor expansions for x(t+h) and x(t-h); the odd terms cancel, leaving a scheme that needs only the force to step.
Position-Verlet: time-reversible and third-order accurate per step, using just one force evaluation and no explicit velocity.
Verlet is symplectic: it exactly preserves phase-space volume and keeps the energy bounded — oscillating around the true value rather than drifting away — for astronomically long runs. That is why molecular-dynamics and solar-system codes use Verlet-family integrators, not RK4, even though RK4 is locally more accurate. The right method respects the physics, not merely the arithmetic.
Worked example: one step, two integrators
Take the mass-spring oscillator \ddot{x}=-\omega^2 x with \omega=1, starting from rest at x_0=1 (so v_0=0), and take a single step of size h=0.1. The exact answer is x(0.1)=\cos(0.1)=0.99500. Watch how the two methods fare on the same one force evaluation.
The setup: a pure oscillator, whose true motion conserves the 'energy' E \propto x^2 + v^2 = 1 for all time.
Euler. The acceleration now is a_0=-\omega^2 x_0=-1. Then x_1 = x_0 + h v_0 = 1 and v_1 = v_0 + h a_0 = -0.1. The position is unchanged at $1$ (error $+0.005$), and the energy \propto x^2+v^2 has already grown from $1 to 1.01$ — Euler is injecting energy, exactly the outward-spiral pathology.
Verlet. It needs a prior point; seed it with x_{-1}=x_0 - h v_0 + \tfrac{1}{2} a_0 h^2 = 1 - 0 - 0.005 = 0.995. Then x_1 = 2x_0 - x_{-1} + a_0 h^2 = 2 - 0.995 - 0.01 = 0.995. That matches \cos(0.1)=0.99500 to all five digits — with the same single force evaluation Euler used.