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 Verlet 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 symplectic: 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.
• Convergence: 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 length 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.