Notebooks · posts/euler-loses-the-planet

Forward Euler loses the planet

Adam Zettel · September 21, 2026 · written with an agent

Ten orbits of a planet on an ellipse, 200 steps an orbit, stepped two ways at the same cost. Forward Euler gains energy until the planet is unbound. Leapfrog's energy error stays under one percent, and it keeps the angular momentum to rounding.

r ¨ = - \frac{r}{∣ r ∣ ^{3}}, E = \frac{1}{2} ∣ r ˙ ∣^{2} - \frac{1}{∣ r ∣}

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

The orbit, in units that keep the numbers small

A planet goes around a sun that sits still at the origin. The only force on it is the sun's gravity, which pulls toward the origin, so its acceleration is

r ¨ = - \frac{GM r}{∣ r ∣ ^{3}},

where $r$ is the planet's position, $G$ is the gravitational constant and $M$ is the sun's mass. The program works in units where $GM = 1$ and the orbit's semi-major axis $a$ (half the long diameter of the ellipse) is 1. Kepler's third law then says one orbit takes $2 π a^{3} / GM = 2 π$ time units.

The eccentricity ECC says how stretched the ellipse is: 0 is a circle, and at 0.6 the closest approach to the sun (the perihelion) is $1 - 0.6 = 0.4$ and the farthest point is $1 + 0.6 = 1.6$ . The program follows the planet for ten orbits with 200 time steps per orbit, and records how far the energy has wandered at SAMPLES evenly spaced moments.

[1]lines 31–35

pi = 3.141592653589793
ECC = 0.6
ORBITS = 10
PER_ORBIT = 200
SAMPLES = 200

The number that should not change

With only gravity acting, the planet's energy per unit mass, kinetic plus potential,

E = \frac{1}{2} ∣ v ∣^{2} - \frac{1}{∣ r ∣},

has the same value at every moment of the true motion. That makes it a check on any numerical method: whatever change in $E$ the computer reports is error the method made.

The planet starts at perihelion, at $(0.4, 0)$ , moving straight up. Its speed there is $(1 + e) / (1 - e) = 2$ (from the vis-viva equation $v^{2} = 2/ r - 1/ a$ , with $r = 0.4$ ). Before running: what is $E$ ? Any closed orbit has $E = - 1/ (2 a)$ , so $- \frac{1}{2}$ here. It is negative because the planet is bound: it lacks the energy to climb out of the sun's pull. If a method pushed $E$ above zero, the simulated planet would leave and never come back.

[2]lines 54–62

#| Energy per unit mass of the state (x, y, vx, vy): kinetic plus potential.
def energy(x, y, vx, vy) {
    return 0.5 * (vx * vx + vy * vy) - 1.0 / sqrt(x * x + y * y)
}

# the state at perihelion: position (x, y), then velocity (vx, vy)
START = [1.0 - ECC, 0.0, 0.0, sqrt((1.0 + ECC) / (1.0 - ECC))]
start_energy = energy(START[0], START[1], START[2], START[3])
show start_energy

The printed start_energy is $- 0.5$ , as predicted. So a method has to add $0.5$ to the energy to set the planet loose, and every bit it adds is error.

Forward Euler lets the start of the step stand in for the whole step

A computer cannot follow the motion continuously. It jumps ahead by a time step $Δ t$ , again and again. Forward Euler is the simplest way to make each jump: take the velocity and the acceleration at the start of the step and pretend they stay fixed for the whole step,

r_{n + 1} = r_{n} + Δ t v_{n}, v_{n + 1} = v_{n} + Δ t a (r_{n}) .

The state lives in one array, s = [x, y, vx, vy], and a step overwrites it.

[3]lines 79–93

#| One forward Euler step of length dt on the state s = [x, y, vx, vy]. The
#| acceleration is taken at the old position and the position moves with the
#| old velocity, so everything the step uses is from its start.
def euler_step(s: Arr, dt) {
    x = s[0]
    y = s[1]
    r3 = pow(x * x + y * y, 1.5)
    ax = -x / r3
    ay = -y / r3
    s[0] = x + dt * s[2]
    s[1] = y + dt * s[3]
    s[2] = s[2] + dt * ax
    s[3] = s[3] + dt * ay
    return 0.0
}

Before running one step: does it raise the energy or lower it? At perihelion the velocity points straight up and the pull points straight left, toward the sun. The position update moves the planet along a straight line, up and off the curve, which leaves it a little farther from the sun, so the potential energy rises. The velocity update adds a sideways kick at right angles to the velocity, which makes the velocity longer, so the kinetic energy rises too. Both parts go up.

[4]lines 103–108

DT = 2.0 * pi / PER_ORBIT
after_one = 1.0 * START
euler_step(after_one, DT)
euler_one_step_gain = energy(after_one[0], after_one[1], after_one[2], after_one[3]) - start_energy
show DT
show euler_one_step_gain

The step is $Δ t = 2 π /200 \approx 0.0314$ , and one Euler step from perihelion adds $0.0496$ to the energy (euler_one_step_gain), a tenth of the $0.5$ that would set the planet loose. Perihelion is where the planet moves fastest and the pull is strongest, so it is where a step's error is largest; there are 2000 steps to come.

Leapfrog kicks, drifts and kicks

Leapfrog does the same work per step, one evaluation of the force, in a different order. It gives the velocity half a step's worth of acceleration (a kick), then moves the position a full step with that updated velocity (a drift), then gives the velocity the other half kick using the force at the new position:

v_{n + \frac{1}{2}} = v_{n} + \frac{Δ t}{2} a (r_{n}), r_{n + 1} = r_{n} + Δ t v_{n + \frac{1}{2}}, v_{n + 1} = v_{n + \frac{1}{2}} + \frac{Δ t}{2} a (r_{n + 1}) .

The second kick's force is the next step's first kick's force, so a program that keeps it pays for one force per step, the same as Euler. This one recomputes it, to keep each step on its own.

[5]lines 131–160

#| One leapfrog step of length dt on the state s = [x, y, vx, vy]: half a kick
#| to the velocity, a full drift of the position, and the other half kick with
#| the force at the new position.
def leapfrog_step(s: Arr, dt) {
    r3 = pow(s[0] * s[0] + s[1] * s[1], 1.5)
    s[2] = s[2] - 0.5 * dt * s[0] / r3
    s[3] = s[3] - 0.5 * dt * s[1] / r3
    s[0] = s[0] + dt * s[2]
    s[1] = s[1] + dt * s[3]
    r3 = pow(s[0] * s[0] + s[1] * s[1], 1.5)
    s[2] = s[2] - 0.5 * dt * s[0] / r3
    s[3] = s[3] - 0.5 * dt * s[1] / r3
    return 0.0
}

after_one_leap = 1.0 * START
leapfrog_step(after_one_leap, DT)
leap_one_step_change = energy(after_one_leap[0], after_one_leap[1], after_one_leap[2], after_one_leap[3]) - start_energy
show leap_one_step_change
# what ten orbits would add if every step erred the way this one did
leap_if_errors_added = ORBITS * PER_ORBIT * leap_one_step_change
show leap_if_errors_added
# how many times Euler's first step changes the energy more than leapfrog's
step_ratio = euler_one_step_gain / leap_one_step_change
show step_ratio

test "one Euler step changes the energy 398 times more than one leapfrog step" {
    # the prose states three digits
    assert_near(step_ratio, 398.0, 0.5)
}

#| One leapfrog step of length dt on the state s = [x, y, vx, vy]: half a kick
#| to the velocity, a full drift of the position, and the other half kick with
#| the force at the new position.
def leapfrog_step(s: Arr, dt) {
    r3 = pow(s[0] * s[0] + s[1] * s[1], 1.5)
    s[2] = s[2] - 0.5 * dt * s[0] / r3
    s[3] = s[3] - 0.5 * dt * s[1] / r3
    s[0] = s[0] + dt * s[2]
    s[1] = s[1] + dt * s[3]
    r3 = pow(s[0] * s[0] + s[1] * s[1], 1.5)
    s[2] = s[2] - 0.5 * dt * s[0] / r3
    s[3] = s[3] - 0.5 * dt * s[1] / r3
    return 0.0
}

after_one_leap = 1.0 * START
leapfrog_step(after_one_leap, DT)
leap_one_step_change = energy(after_one_leap[0], after_one_leap[1], after_one_leap[2], after_one_leap[3]) - start_energy
show leap_one_step_change
# what ten orbits would add if every step erred the way this one did
leap_if_errors_added = ORBITS * PER_ORBIT * leap_one_step_change
show leap_if_errors_added
# how many times Euler's first step changes the energy more than leapfrog's
step_ratio = euler_one_step_gain / leap_one_step_change
show step_ratio

test "one Euler step changes the energy 398 times more than one leapfrog step" {
    # the prose states three digits
    assert_near(step_ratio, 398.0, 0.5)
}

From the same starting point, one leapfrog step changes the energy by $0.000125$ (leap_one_step_change), 398 times less than Euler's step (step_ratio). That is not yet the difference that matters. If each of the 2000 steps added that much, ten orbits would add $0.25$ (leap_if_errors_added). What matters is whether the errors of many steps pile up or cancel.

Running ten orbits

The next function runs either method for ten orbits and records the energy error $∣ E - E_{0} ∣$ , where $E_{0} = - \frac{1}{2}$ is the starting energy, at SAMPLES evenly spaced times. It records the base-10 logarithm, because the errors range over many powers of ten and a logarithm puts them on one readable scale; the $1 0^{- 16}$ added inside keeps the logarithm finite when the error is exactly zero. The planet's final state is left in FINAL for the checks further down.

[6]lines 178–203

#| Integrate from perihelion and return log10 |E - E0| at SAMPLES times.
#| method 0 is forward Euler, 1 is leapfrog. FINAL gets the state at the end.
FINAL = zeros(len(START))
def energy_error(method, per_orbit) {
    dt = 2.0 * pi / per_orbit
    steps = ORBITS * per_orbit
    every = steps / SAMPLES
    s = 1.0 * START
    e0 = energy(s[0], s[1], s[2], s[3])
    out = zeros(SAMPLES)
    for n to steps {
        if method == 0.0 {
            euler_step(s, dt)
        }
        else {
            leapfrog_step(s, dt)
        }
        if mod(n + 1.0, every) == 0.0 {
            out[(n + 1.0) / every - 1.0] = log(abs(energy(s[0], s[1], s[2], s[3]) - e0) + 1.0e-16) / log(10.0)
        }
    }
    for k to len(s) {
        FINAL[k] = s[k]
    }
    return out
}

Forward Euler, ten orbits

Before running: 2000 steps, each adding a little energy. Does the total reach the $0.5$ that unbinds the planet?

[7]lines 212–222

YEARS = linspace(ORBITS / SAMPLES, ORBITS, SAMPLES)
euler = energy_error(0.0, PER_ORBIT)
euler_final = pow(10.0, euler[SAMPLES - 1.0])
euler_end = 1.0 * FINAL
euler_end_energy = energy(euler_end[0], euler_end[1], euler_end[2], euler_end[3])
show euler_final
show euler_end_energy

test "forward Euler has gained more than the binding energy 1/2: the planet is unbound" {
    assert_gt(euler_final, 0.5)
}

It does. After ten orbits Euler has added $0.546$ (euler_final), so the planet's energy is $+ 0.046$ (euler_end_energy): positive, which means the simulated planet is no longer in orbit at all. It is flying away. The real planet would go around forever.

Leapfrog, ten orbits, same step

Before running: if leapfrog's errors simply added up, the $0.25$ from the cell above would be the answer, half of what it takes to unbind the planet.

[8]lines 235–249

leap = energy_error(1.0, PER_ORBIT)
leap_worst = pow(10.0, max(leap))
leap_end = 1.0 * FINAL
show leap_worst
# the worst error as a share of the orbit's binding energy, 1/2
leap_worst_share = leap_worst / 0.5
show leap_worst_share

test "leapfrog's worst energy error is under one percent of the orbit's energy" {
    assert_lt(leap_worst, 0.01 * 0.5)
}
test "leapfrog's worst energy error is 0.74% of the orbit's energy" {
    # to the two digits the prose states
    assert_near(leap_worst_share, 0.0074, 0.00005)
}

The largest energy error leapfrog makes in ten orbits is $0.0037$ (leap_worst), which is $0.74%$ of the orbit's $0.5$ (leap_worst_share) and almost seventy times less than $0.25$ . Its errors do not add up.

The figure shows both runs on a logarithmic scale. Before looking: Euler's error is largest at perihelion, and the planet reaches perihelion once an orbit. Does Euler's error grow steadily, or in jumps?

[9]lines 260–303

figure("Energy error over ten orbits, equal step size", "orbits", "$\\log_{10}$ of the energy error")
plot(YEARS, euler, "forward Euler")
plot(YEARS, leap, "leapfrog")

# Euler's error after the first orbit, and how long the orbit it put the planet
# on takes: a bound orbit with energy E has a = 1 / (2 |E|) and period a^(3/2)
euler_error_after_1_orbit = pow(10.0, euler[SAMPLES / ORBITS - 1.0])
euler_orbit_period = pow(1.0 / (2.0 * abs(start_energy + euler_error_after_1_orbit)), 1.5)
show euler_error_after_1_orbit
show euler_orbit_period
# leapfrog's worst error in the first orbit and in the tenth
leap_first_orbit = -20.0
leap_tenth_orbit = -20.0
for k to SAMPLES / ORBITS {
    leap_first_orbit = max(leap_first_orbit, leap[k])
    leap_tenth_orbit = max(leap_tenth_orbit, leap[SAMPLES - SAMPLES / ORBITS + k])
}
leap_first_orbit = pow(10.0, leap_first_orbit)
leap_tenth_orbit = pow(10.0, leap_tenth_orbit)
show leap_first_orbit
show leap_tenth_orbit
# leapfrog's typical error, the median of the samples, and the share of the
# samples that round to the same four decimals
leap_sorted = zeros(SAMPLES)
for k to SAMPLES {
    leap_sorted[k] = pow(10.0, leap[k])
}
sort(leap_sorted)
leap_median = 0.5 * (leap_sorted[SAMPLES / 2 - 1] + leap_sorted[SAMPLES / 2])
leap_near_median = 0.0
for k to SAMPLES {
    if abs(leap_sorted[k] - leap_median) < 0.00005 {
        leap_near_median = leap_near_median + 1.0
    }
}
leap_near_median = leap_near_median / SAMPLES
show leap_median
show leap_near_median

test "leapfrog's error sits near 0.0036 for most of each orbit" {
    assert_near(leap_median, 0.0036, 0.00005)
    # most: more than half of the samples are within 0.00005 of it
    assert_gt(leap_near_median, 0.5)
}

figure("Energy error over ten orbits, equal step size", "orbits", "$\\log_{10}$ of the energy error")
plot(YEARS, euler, "forward Euler")
plot(YEARS, leap, "leapfrog")

# Euler's error after the first orbit, and how long the orbit it put the planet
# on takes: a bound orbit with energy E has a = 1 / (2 |E|) and period a^(3/2)
euler_error_after_1_orbit = pow(10.0, euler[SAMPLES / ORBITS - 1.0])
euler_orbit_period = pow(1.0 / (2.0 * abs(start_energy + euler_error_after_1_orbit)), 1.5)
show euler_error_after_1_orbit
show euler_orbit_period
# leapfrog's worst error in the first orbit and in the tenth
leap_first_orbit = -20.0
leap_tenth_orbit = -20.0
for k to SAMPLES / ORBITS {
    leap_first_orbit = max(leap_first_orbit, leap[k])
    leap_tenth_orbit = max(leap_tenth_orbit, leap[SAMPLES - SAMPLES / ORBITS + k])
}
leap_first_orbit = pow(10.0, leap_first_orbit)
leap_tenth_orbit = pow(10.0, leap_tenth_orbit)
show leap_first_orbit
show leap_tenth_orbit
# leapfrog's typical error, the median of the samples, and the share of the
# samples that round to the same four decimals
leap_sorted = zeros(SAMPLES)
for k to SAMPLES {
    leap_sorted[k] = pow(10.0, leap[k])
}
sort(leap_sorted)
leap_median = 0.5 * (leap_sorted[SAMPLES / 2 - 1] + leap_sorted[SAMPLES / 2])
leap_near_median = 0.0
for k to SAMPLES {
    if abs(leap_sorted[k] - leap_median) < 0.00005 {
        leap_near_median = leap_near_median + 1.0
    }
}
leap_near_median = leap_near_median / SAMPLES
show leap_median
show leap_near_median

test "leapfrog's error sits near 0.0036 for most of each orbit" {
    assert_near(leap_median, 0.0036, 0.00005)
    # most: more than half of the samples are within 0.00005 of it
    assert_gt(leap_near_median, 0.5)
}

In jumps. Euler's error is $0.31$ after the first orbit (euler_error_after_1_orbit), nearly all of it gained in the few steps around the starting perihelion. That leaves the planet on a larger orbit, one that takes $4.3$ of the original orbits to go around (euler_orbit_period), so the next close pass to the sun comes near orbit 4.3 on the plot. That pass adds the rest, the total goes past $0.5$ , and the planet leaves; there is no third pass.

Leapfrog's error sits near $0.0036$ for most of each orbit and dips as the planet swings past the sun. Its worst in the first orbit is $0.0037$ (leap_first_orbit) and its worst in the tenth is $0.0037$ (leap_tenth_orbit): the same level, nine orbits later.

Why leapfrog's errors do not pile up

Leapfrog is symplectic. Each leapfrog step is the exact motion of a slightly different system, one whose energy function differs from the true $E$ by an amount that shrinks with $Δ t^{2}$ . The steps hold that modified energy almost exactly constant, so the true $E$ can only wobble around it; it cannot drift away. Forward Euler has no such conserved quantity, so each step's error is free to add to the last.

A second conserved quantity makes the same point more directly. The angular momentum per unit mass, $L = x v_{y} - y v_{x}$ , is constant for any force that points at the origin. At the start it is $0.4 \times 2 = 0.8$ , which equals $1 - e^{2}$ for this orbit.

Before running: which method keeps $L$ at $0.8$ ? A leapfrog kick changes the velocity along $r$ , and $r \times r = 0$ , so a kick leaves $L$ alone. A drift changes the position along $v$ , and $v \times v = 0$ , so a drift leaves it alone too. Leapfrog should keep $L$ to rounding error. Euler moves the position and the velocity at the same time, from the old values; working out the product gives a change of $Δ t^{2} (v_{x} a_{y} - v_{y} a_{x})$ per step, which is not zero.

[10]lines 336–350

#| The angular momentum per unit mass of a state [x, y, vx, vy].
def angular_momentum(s: Arr) {
    return s[0] * s[3] - s[1] * s[2]
}

exact_momentum = sqrt(1.0 - ECC * ECC)
euler_momentum_drift = angular_momentum(euler_end) - exact_momentum
momentum_drift = abs(angular_momentum(leap_end) - exact_momentum)
show exact_momentum
show euler_momentum_drift
show momentum_drift

test "leapfrog keeps the angular momentum sqrt(1 - e^2) to rounding" {
    assert_lt(momentum_drift, 1.0e-12)
}

Euler's angular momentum has grown by $0.330$ (euler_momentum_drift), about two fifths of its starting value. Leapfrog's drift is at the level of rounding error, $2.9 \times 1 0^{- 15}$ (momentum_drift), and the test holds it under $1 0^{- 12}$ .

How the error shrinks with the step

A method is order p if its error shrinks like $Δ t^{p}$ : halve the step and the error drops by a factor of $2^{p}$ . Leapfrog is second order, so half the step should give a quarter of the error. The cell below runs leapfrog again with 400 steps per orbit and measures $p$ from the two worst errors:

p = \frac{lo g _{10} ( error at Δ t ) - lo g _{10} ( error at Δ t /2 )}{lo g _{10} 2} .

[11]lines 367–374

FINE_PER_ORBIT = 2.0 * PER_ORBIT
leap_fine = energy_error(1.0, FINE_PER_ORBIT)
order = (max(leap) - max(leap_fine)) / (log(2.0) / log(10.0))
show order

test "and it is second order: half the step, a quarter of the error" {
    assert_near(order, 2.0, 0.1)
}

The measured order is $2.006$ (order). Twice the work per orbit buys four times the accuracy, and the error stays bounded however long the run. Forward Euler is first order, and its error keeps growing as the run gets longer.

[12]lines 382–383

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors