Simulating motion

Why it matters

Newton's second law turns almost every mechanics problem into a differential equation: a = F(x) / m. Only a few cases (a free fall, a spring, a single planet) have neat solutions. For everything else (three bodies, a molecule of a thousand atoms, a galaxy) you step the equations forward in time on a computer. The choice of stepping method decides whether your simulated planet keeps orbiting or spirals away.

The setup: an orbit in convenient units

Use astronomical units (AU) for distance and years for time. Then the Sun's gravitational parameter is GM = 4 pi^2 AU^3 / yr^2, and the Earth's orbit is a circle of radius 1 with speed 2 pi AU/yr. The acceleration of a planet at position (x, y) is

r = sqrt(x^2 + y^2)
a = -GM * (x, y) / r^3

Euler's method and why it fails

The obvious approach updates position and velocity with the current values:

x_new = x + v * dt
v_new = v + a(x) * dt

Its error per step is of order dt^2, and worse, the errors always push the same way for orbits: energy rises every revolution. With dt = 0.001 yr for one year, the orbit's energy grows by about 7% and the radius drifts from 1 to about 1.08 AU. Over a thousand years the planet escapes.

Velocity Verlet

Velocity VerletA time-stepping method that updates velocity with the average of the accelerations at both ends of the step; second order and symplectic. uses the acceleration at both ends of the step:

x_new = x + v * dt + 0.5 * a * dt^2
a_new = a(x_new)
v_new = v + 0.5 * (a + a_new) * dt

It costs one force evaluation per step, the same as Euler, but its error is second order and it is symplecticDescribes an integrator that respects the structure of Newtonian mechanics, so energy wobbles but does not drift.: it respects the structure of Newtonian mechanics, so the energy wobbles slightly but does not drift. For the same one-year Earth orbit with dt = 0.001, the planet returns to within about 1e-4 AU of its start and the relative energy error stays tiny.

import math

GM = 4 * math.pi ** 2

def accel(x, y):
    r3 = (x * x + y * y) ** 1.5
    return -GM * x / r3, -GM * y / r3

x, y, vx, vy, dt = 1.0, 0.0, 0.0, 2 * math.pi, 0.001
ax, ay = accel(x, y)
for _ in range(1000):                 # one year
    x += vx * dt + 0.5 * ax * dt * dt
    y += vy * dt + 0.5 * ay * dt * dt
    nax, nay = accel(x, y)
    vx += 0.5 * (ax + nax) * dt
    vy += 0.5 * (ay + nay) * dt
    ax, ay = nax, nay
print(x, y)                           # close to (1, 0)

Checking a simulation

Never trust an integrator just because it runs. Check quantities that physics says are conserved:

• Energy per unit mass: E = v^2 / 2 - GM / r. For a circular Earth orbit it is -2 pi^2, about -19.74.
• Angular momentum per unit mass: L = x * vy - y * vx, here 2 pi.
• ConvergenceThe check that halving dt changes the answer by the expected factor (about 4 for a second-order method).: halve dt and rerun. For a second-order method the error should drop by about 4.

From two bodies to N bodies

With many bodies, each feels the pull of all the others, so a direct force calculation costs about N^2 operations per step. Codes add a small softening lengthA small length added in N-body codes to avoid infinite forces when two bodies pass very close. to avoid infinite forces when two bodies pass very close, and large simulations use tree or mesh methods to approach N log N. Molecular dynamics uses exactly the same Verlet step with interatomic forces instead of gravity. Libraries such as SciPy's solve_ivp offer adaptive Runge-Kutta methods; they are great for short, accurate runs, but for long conservative simulations a symplectic method like Verlet usually behaves better.

Choosing the time step

The step must be small compared with the fastest motion in the system. A rule of thumb for orbits is at least a few hundred steps per orbit; for close encounters you need much smaller steps or an adaptive scheme. If halving dt changes your answer noticeably, it was too large.

Common mistakes

• Using the new position with the old acceleration twice, which turns Verlet back into a first-order method.
• Forgetting the 0.5 factors in the position and velocity updates.
• Wrong sign on gravity. Acceleration points toward the Sun: -GM x / r^3.
• Judging by a plot alone. A drifting energy shows a problem long before the orbit looks wrong.

Recap

• Euler is simple but drifts; velocity Verlet is just as cheap and conserves energy far better.
• Work in natural units (AU, years, GM = 4 pi^2) to keep numbers near 1.
• Validate with conserved quantities and a time-step convergence test.

Try it

Open the Do tab to implement velocity Verlet for an orbit and check that energy is conserved.