Classical mechanics

Around 1600, Johannes Kepler, working from Tycho Brahe's careful naked-eye observations of Mars, found three empirical laws: planets move on ellipses with the Sun at a focus; the line from the Sun to a planet sweeps out equal areas in equal times; and the square of the period is proportional to the cube of the orbit's size. Two generations later Isaac Newton showed that all three follow from a single differential equation — his second law with a force varying as the inverse square of distance. It was the first great triumph of mathematical physics, and its methods are still the ones engineers use to put a satellite into orbit over Nairobi or to design a building in Santiago that survives earthquakes.

This lesson develops classical mechanics as mathematics: Newton's laws as ordinary differential equations, conservation laws as the tool that makes them solvable, the Kepler problem solved completely, oscillators from the pendulum to resonance, and Lagrange's reformulation through the principle of stationary action. It assumes ordinary differential equations and multivariable calculus (gradients, the chain rule, line integrals); the orbits are the conic sections of school geometry.

The state of a particle at an instant is its position and its velocity: knowing where a ball is does not tell you where it will be, but knowing where it is and how it is moving does — once you know the forces. Newton's second law $m\mathbf r'' = \mathbf F$ is a second-order ODE, so by the Picard–Lindelöf theorem position and velocity at one instant determine the motion (for Lipschitz forces).

Solving such equations in closed form is usually impossible, but some combinations of position and velocity do not change along the motion — conserved quantities, such as energy and angular momentum. Each one confines the motion to a smaller set and lowers the effective dimension of the problem. For a ball rolling in a valley, conservation of energy alone tells you how high it will climb and how fast it moves at the bottom. For a planet, energy and angular momentum together are enough to find the orbit exactly.

A particle of mass $m > 0$ has position $\mathbf r(t) \in \R^3$, velocity $\mathbf v = \mathbf r'$, acceleration $\mathbf a = \mathbf r''$ and momentum $\mathbf p = m\mathbf v$.

The weight of a body is the gravitational force on it, $\mathbf W = m\mathbf g$ near the Earth's surface, where $|\mathbf g| \approx 9.8\ \mathrm{m/s^2}$. Mass is measured in kilograms; weight, a force, in newtons.

In one dimension, energy conservation reduces Newton's equation to a first-order one: $\tfrac12 m x'^2 + V(x) = E$ gives $x' = \pm\sqrt{2(E - V(x))/m}$. The motion is confined to $\{x : V(x) \le E\}$, reverses at turning points where $V(x) = E$, and between turning points $x_1 < x_2$ (with $V'(x_i) \ne 0$) it is periodic with period

For the pendulum $\theta'' = -(g/\ell)\sin\theta$ released from rest at angle $\theta_0$, this integral is an elliptic integral, $T = 4\sqrt{\ell/g}\,K\big(\sin\tfrac{\theta_0}{2}\big)$ — larger than the small-angle value $2\pi\sqrt{\ell/g}$, by about $18\%$ for $\theta_0 = 90^\circ$ (see the computation below).

The relative position $\mathbf r = \mathbf r_1 - \mathbf r_2$ satisfies $\mathbf r'' = \mathbf F_{12}/m_1 - \mathbf F_{21}/m_2 = \mathbf F_{12}\left(\frac1{m_1} + \frac1{m_2}\right)$, i.e. $\mu\mathbf r'' = \mathbf F_{12}$ with the reduced mass $\mu = \frac{m_1 m_2}{m_1 + m_2}$. So the two-body problem reduces to a one-body problem; for gravity, $\mathbf F_{12} = -Gm_1m_2\,\mathbf r/|\mathbf r|^3$ gives $\mathbf r'' = -G(m_1 + m_2)\,\mathbf r/|\mathbf r|^3$, exactly the problem of one body attracted to a fixed centre of "mass" $m_1 + m_2$.

Take the plane of motion with polar coordinates $(r, \theta)$, and unit vectors $\mathbf e_r = (\cos\theta, \sin\theta)$ and $\mathbf e_\theta = (-\sin\theta, \cos\theta)$. Differentiating $\mathbf r = r\,\mathbf e_r$ with $\mathbf e_r' = \theta'\mathbf e_\theta$ and $\mathbf e_\theta' = -\theta'\mathbf e_r$:

The angular momentum has magnitude $L = m r^2\theta'$.

For a central force $\mathbf F = f(r)\,\mathbf e_r$, the radial component of Newton's law is $m(r'' - r\theta'^2) = f(r)$. Using $\theta' = L/(mr^2)$ to eliminate $\theta$:

The radial motion is one-dimensional motion in the effective potential $V_{\text{eff}}$; the term $L^2/2mr^2$ is a "centrifugal barrier" that keeps a body with $L \ne 0$ away from the centre.

For two bodies (a planet that is not negligibly light), replace $m$ by the reduced mass $\mu$ and $GM$ by $G(M + m)$; Kepler's third law then holds up to the factor $1 + m/M$, which is $1.001$ for Jupiter. The value $a = -k/(2E)$, obtained by combining $a = p/(1 - e^2)$ with $e^2 - 1 = 2pE/k$, shows that the size of an orbit depends only on its energy.

Add friction proportional to velocity: $m x'' + c\,x' + kx = 0$, or $x'' + 2\gamma x' + \omega_0^2 x = 0$ with $2\gamma = c/m$ and $\omega_0^2 = k/m$. The energy $E = \tfrac12 m x'^2 + \tfrac12 kx^2$ now satisfies $E' = x'(mx'' + kx) = -c\,x'^2 \le 0$: it decreases. The characteristic roots $-\gamma \pm \sqrt{\gamma^2 - \omega_0^2}$ give three regimes (see the ODE lesson):

• underdamped ($\gamma < \omega_0$): $x(t) = Ce^{-\gamma t}\cos(\omega t - \varphi)$ with $\omega = \sqrt{\omega_0^2 - \gamma^2}$ — decaying oscillations;
• critically damped ($\gamma = \omega_0$): $x(t) = (A + Bt)e^{-\gamma t}$ — the fastest return without oscillating, roughly what a car's shock absorbers aim for;
• overdamped ($\gamma > \omega_0$): a sum of two decaying exponentials.

When the damping is small, the amplitude is largest for driving frequencies near the natural frequency $\omega_0$: at $\Omega = \omega_0$ it is $\frac{F_0}{2m\gamma\omega_0}$, which grows without bound as $\gamma \to 0$. This is resonance — useful in radio tuners and musical instruments, dangerous in machinery and structures.

Newton's law needs Cartesian coordinates and forces, including the forces of constraint that keep a bead on its wire or a pendulum bob on its circle. Lagrange (1788) found a formulation in any coordinates, with constraints built in. Its modern form is Hamilton's principle (1834): the motion between times $t_1$ and $t_2$ makes the action $S[q] = \int_{t_1}^{t_2}\mathcal L\,dt$ stationary among all paths with the same endpoints.

These are the simplest cases of Noether's theorem (Emmy Noether, 1918): every continuous symmetry of the action gives a conserved quantity — translations give momentum, rotations angular momentum, time translations energy.

Mathlib has no physics, but it has the calculus to state Newton's law exactly: a motion is a pair $\mathbf r, \mathbf v : \R \to E$ with HasDerivAt r (v t) t and HasDerivAt v (a t) t. Each conservation law is then proved as in the lesson: differentiate the conserved quantity, get $0$, and conclude it is constant.

First, uniform gravity — position, velocity and acceleration of a projectile:

/-- A projectile in uniform gravity: the height `y(t) = y₀ + v₀ t - g t²/2` has velocity
`v(t) = v₀ - g t`, whose derivative is the constant acceleration `-g`. -/
theorem projectile (y₀ v₀ g t : ℝ) :
    HasDerivAt (fun s => y₀ + v₀ * s - g * s ^ 2 / 2) (v₀ - g * t) t ∧
      HasDerivAt (fun s => v₀ - g * s) (-g) t := by
  constructor
  · have h := (((hasDerivAt_id' t).const_mul v₀).const_add y₀).sub
      ((hasDerivAt_pow 2 t).const_mul g |>.div_const 2)
    convert h using 1
    push_cast
    ring
  · have h := ((hasDerivAt_id' t).const_mul g).const_sub v₀
    convert h using 1
    ring

Conservation of energy, in any real inner product space (so in $\R^3$, but also for systems with many particles), for a force $-\nabla V$ given by Mathlib's HasGradientAt. The proof is the human one: $E' = m\langle \mathbf v, \mathbf v'\rangle + \langle\nabla V(\mathbf r), \mathbf r'\rangle = 0$.

/-- **Conservation of energy** for a conservative force, in any (real inner product) space.
Suppose the particle of mass `m ≠ 0` moves by Newton's law with force `F = -∇V`: `r' = v` and
`v' = -(1/m) ∇V(r)`. Then the total energy `E = ½ m ⟪v, v⟫ + V(r)` is constant. Proof: by the
product and chain rules, `E' = m ⟪v, v'⟫ + ⟪∇V(r), r'⟫ = -⟪v, ∇V(r)⟫ + ⟪∇V(r), v⟫ = 0`. -/
theorem energy_conservation {E : Type*} [NormedAddCommGroup E] [InnerProductSpace ℝ E]
    [CompleteSpace E] {m : ℝ} (hm : m ≠ 0) {V : E → ℝ} {gradV : E → E}
    (hV : ∀ x, HasGradientAt V (gradV x) x) {r v : ℝ → E}
    (hr : ∀ t, HasDerivAt r (v t) t) (hv : ∀ t, HasDerivAt v (-(1 / m) • gradV (r t)) t)
    (s t : ℝ) :
    m / 2 * inner ℝ (v s) (v s) + V (r s) = m / 2 * inner ℝ (v t) (v t) + V (r t) := by
  have hE : ∀ τ, HasDerivAt (fun τ => m / 2 * inner ℝ (v τ) (v τ) + V (r τ)) 0 τ := by
    intro τ
    have hkin := ((hv τ).inner ℝ (hv τ)).const_mul (m / 2)
    have hpot := (hasGradientAt_iff_hasFDerivAt.mp (hV (r τ))).comp_hasDerivAt τ (hr τ)
    have h := hkin.add hpot
    convert h using 1
    · rfl
    · simp only [InnerProductSpace.toDual_apply_apply, inner_smul_right, inner_smul_left,
        real_inner_comm (gradV (r τ)) (v τ), RCLike.conj_to_real]
      field_simp
      ring
  exact is_const_of_deriv_eq_zero (fun τ => (hE τ).differentiableAt) (fun τ => (hE τ).deriv) s t

Angular momentum. The library facts $\mathbf v \times \mathbf v = 0$ and $\mathbf r \cdot (\mathbf r \times \mathbf v) = 0$, then conservation of $\mathbf r \times \mathbf v$ under any force parallel to $\mathbf r$ (checked coordinate by coordinate), and its consequence that the motion is planar:

/-- The cross product of a vector with itself is zero (Mathlib's `cross_self`). It is why the
term `v × (m v)` drops out of `d/dt (r × m v) = v × m v + r × m v'`. -/
theorem cross_self_zero (v : Fin 3 → ℝ) : v ⨯₃ v = 0 :=
  cross_self v

/-- **Conservation of angular momentum** under a central force. If `r' = v` and `v' = k(t) r`
(acceleration always parallel to the position vector, as for any central force), then
`L = r × v` has derivative `v × v + r × (k r) = 0`, so it is constant. We check the derivative
coordinate by coordinate (`hasDerivAt_pi`). -/
theorem angular_momentum_conserved {r v : ℝ → Fin 3 → ℝ} {k : ℝ → ℝ}
    (hr : ∀ t, HasDerivAt r (v t) t) (hv : ∀ t, HasDerivAt v (k t • r t) t) (s t : ℝ) :
    r s ⨯₃ v s = r t ⨯₃ v t := by
  have hL : ∀ τ, HasDerivAt (fun τ => r τ ⨯₃ v τ) 0 τ := by
    intro τ
    rw [hasDerivAt_pi]
    have hri : ∀ i, HasDerivAt (fun τ => r τ i) (v τ i) τ := hasDerivAt_pi.mp (hr τ)
    have hvi : ∀ i, HasDerivAt (fun τ => v τ i) (k τ * r τ i) τ := by
      intro i
      simpa using hasDerivAt_pi.mp (hv τ) i
    intro i
    fin_cases i
    · have h := ((hri 1).mul (hvi 2)).sub ((hri 2).mul (hvi 1))
      convert h using 1
      · funext τ; simp [cross_apply]
      · simp; ring
    · have h := ((hri 2).mul (hvi 0)).sub ((hri 0).mul (hvi 2))
      convert h using 1
      · funext τ; simp [cross_apply]
      · simp; ring
    · have h := ((hri 0).mul (hvi 1)).sub ((hri 1).mul (hvi 0))
      convert h using 1
      · funext τ; simp [cross_apply]
      · simp; ring
  exact is_const_of_deriv_eq_zero (fun τ => (hL τ).differentiableAt) (fun τ => (hL τ).deriv) s t

/-- A vector is perpendicular to its cross product with any other vector:
`r · (r × v) = 0` (Mathlib's `dot_self_cross`). -/
theorem dot_self_cross_zero (r v : Fin 3 → ℝ) : r ⬝ᵥ (r ⨯₃ v) = 0 :=
  dot_self_cross r v

/-- **Central-force motion is planar.** Under the hypotheses of `angular_momentum_conserved`,
the position `r(t)` is perpendicular to the constant vector `L = r(0) × v(0)` at every time,
so (when `L ≠ 0`) the particle stays in the plane through the origin perpendicular to `L`. -/
theorem motion_in_plane {r v : ℝ → Fin 3 → ℝ} {k : ℝ → ℝ}
    (hr : ∀ t, HasDerivAt r (v t) t) (hv : ∀ t, HasDerivAt v (k t • r t) t) (t : ℝ) :
    r t ⬝ᵥ (r 0 ⨯₃ v 0) = 0 := by
  rw [angular_momentum_conserved hr hv 0 t]
  exact dot_self_cross_zero (r t) (v t)

Kepler's first law, in two formal pieces: the conic $u = (1 + e\cos\theta)/p$ solves Binet's equation $\ddot u + u = 1/p$, and a curve $r(1 + e\cos\theta) = p$ is the conic $(1 - e^2)x^2 + 2pex + y^2 = p^2$ in Cartesian coordinates:

/-- The orbit equation. With `u = 1/r` as a function of the polar angle `θ`, the inverse-square
law gives Binet's equation `u'' + u = 1/p`. Its solutions `u(θ) = (1 + e cos θ)/p` (for constants
`e`, `p` with `p ≠ 0`) indeed satisfy it: `u' = -(e/p) sin θ`, `u'' = -(e/p) cos θ`. -/
theorem binet_conic (e p θ : ℝ) :
    HasDerivAt (fun φ => (1 + e * cos φ) / p) (-(e / p) * sin θ) θ ∧
      HasDerivAt (fun φ => -(e / p) * sin φ) (1 / p - (1 + e * cos θ) / p) θ := by
  constructor
  · have h := (((hasDerivAt_cos θ).const_mul e).const_add 1).div_const p
    convert h using 1
    ring
  · have h := (hasDerivAt_sin θ).const_mul (-(e / p))
    convert h using 1
    ring

/-- **Kepler's first law, the geometric step.** If a point at distance `r ≥ 0` from the origin and
polar angle `θ` lies on the orbit `r (1 + e cos θ) = p`, then its Cartesian coordinates
`x = r cos θ`, `y = r sin θ` satisfy `(1 - e²) x² + 2 p e x + y² = p²`: a conic with a focus at
the origin — an ellipse if `0 ≤ e < 1` (a circle if `e = 0`), a parabola if `e = 1`, a hyperbola
branch if `e > 1`. Proof: `r = p - e x`, and `x² + y² = r²`. -/
theorem orbit_is_conic {r θ e p : ℝ} (horbit : r * (1 + e * cos θ) = p) :
    (1 - e ^ 2) * (r * cos θ) ^ 2 + 2 * p * e * (r * cos θ) + (r * sin θ) ^ 2 = p ^ 2 := by
  have hr : r = p - e * (r * cos θ) := by linarith
  have hpyth : (r * cos θ) ^ 2 + (r * sin θ) ^ 2 = r ^ 2 := by
    have := sin_sq_add_cos_sq θ
    calc (r * cos θ) ^ 2 + (r * sin θ) ^ 2 = r ^ 2 * (sin θ ^ 2 + cos θ ^ 2) := by ring
      _ = r ^ 2 := by rw [this, mul_one]
  calc (1 - e ^ 2) * (r * cos θ) ^ 2 + 2 * p * e * (r * cos θ) + (r * sin θ) ^ 2
      = ((r * cos θ) ^ 2 + (r * sin θ) ^ 2) - (e * (r * cos θ)) ^ 2
          + 2 * p * (e * (r * cos θ)) := by ring
    _ = r ^ 2 - (e * (r * cos θ)) ^ 2 + 2 * p * (e * (r * cos θ)) := by rw [hpyth]
    _ = (p - e * (r * cos θ)) ^ 2 - (e * (r * cos θ)) ^ 2 + 2 * p * (e * (r * cos θ)) := by
        rw [← hr]
    _ = p ^ 2 := by ring

The two algebraic computations behind the relation between energy and eccentricity, and behind the third law (the inputs — the perihelion distance, the area of the ellipse and the areal velocity — are taken from the human proofs):

/-- **Energy determines the shape.** At perihelion the radial velocity is `0` and
`r = p/(1 + e)`, so the energy is `E = L²/(2 m r²) - k/r`. With `p = L²/(m k)` this equals
`k (e² - 1)/(2p)`: hence `e < 1` (an ellipse) exactly when `E < 0`. -/
theorem eccentricity_energy {m k L e : ℝ} (hm : 0 < m) (hk : 0 < k) (hL : L ≠ 0)
    (he : 0 ≤ e) :
    let p := L ^ 2 / (m * k)
    let r := p / (1 + e)
    L ^ 2 / (2 * m * r ^ 2) - k / r = k * (e ^ 2 - 1) / (2 * p) := by
  intro p r
  have hp : p ≠ 0 := by positivity
  have h1e : 1 + e ≠ 0 := by positivity
  have hr : r ≠ 0 := div_ne_zero hp h1e
  have hpL : L ^ 2 = p * (m * k) := by
    simp only [p]; field_simp
  simp only [r]
  rw [hpL]
  field_simp
  ring

/-- **Kepler's third law, the algebraic step.** An ellipse with semi-latus rectum
`p = L²/(m k)` and semi-major axis `a` has semi-minor axis `b` with `b² = p a` (since
`a = p/(1 - e²)` and `b = a √(1 - e²)`). Its area `πab` is swept at the constant rate `L/(2m)`,
so the period is `T = πab / (L/(2m))`. Then `T² = 4π² m a³ / k`. -/
theorem kepler_third_law {m k L a b T : ℝ} (hm : 0 < m) (hk : 0 < k) (hL : L ≠ 0)
    (hb : b ^ 2 = L ^ 2 / (m * k) * a) (hT : T = π * a * b / (L / (2 * m))) :
    T ^ 2 = 4 * π ^ 2 * m * a ^ 3 / k := by
  have hm0 : m ≠ 0 := hm.ne'
  have hk0 : k ≠ 0 := hk.ne'
  rw [hT]
  calc (π * a * b / (L / (2 * m))) ^ 2 = 4 * π ^ 2 * m ^ 2 * a ^ 2 * b ^ 2 / L ^ 2 := by
        field_simp
        ring
    _ = 4 * π ^ 2 * m * a ^ 3 / k := by
        rw [hb]
        field_simp

The damped oscillator: $x = e^{-\gamma t}\cos\omega t$ solves $x'' + 2\gamma x' + (\gamma^2 + \omega^2)x = 0$:

/-- The underdamped oscillator `x'' + 2γ x' + ω₀² x = 0` with `ω₀² = γ² + ω²` is solved by
`x(t) = e^{-γt} cos(ωt)`: we compute its velocity `v` and check `v' = -2γ v - ω₀² x`. -/
theorem damped_oscillator (γ ω t : ℝ) :
    HasDerivAt (fun s => exp (-γ * s) * cos (ω * s))
        (exp (-γ * t) * (-γ * cos (ω * t) - ω * sin (ω * t))) t ∧
      HasDerivAt (fun s => exp (-γ * s) * (-γ * cos (ω * s) - ω * sin (ω * s)))
        (-2 * γ * (exp (-γ * t) * (-γ * cos (ω * t) - ω * sin (ω * t)))
          - (γ ^ 2 + ω ^ 2) * (exp (-γ * t) * cos (ω * t))) t := by
  have hl : ∀ c : ℝ, HasDerivAt (fun s => c * s) c t := fun c => by
    simpa using (hasDerivAt_id' t).const_mul c
  have he := (hl (-γ)).exp
  constructor
  · have h := he.mul (hl ω).cos
    convert h using 1
    ring
  · have h := he.mul ((((hl ω).cos).const_mul (-γ)).sub (((hl ω).sin).const_mul ω))
    convert h using 1
    simp only [Pi.sub_apply]
    ring

What is and is not formalised. Independently proved: the projectile, energy conservation for conservative forces, conservation of angular momentum and planarity under central forces, the Binet solution, the conic identity, the algebra of the energy–eccentricity relation and of the third law, and the damped oscillator. Cited from Mathlib: the two cross-product identities. Not formalised: the derivation of Binet's equation from Newton's law (the change of parameter from $t$ to $\theta$), the area of an ellipse and the areal-velocity argument of Kepler's second law, the reduction of the two-body problem, the resonance formula, and the Euler–Lagrange equations, whose proof needs differentiation under the integral sign and the fundamental lemma — all within reach of Mathlib, but not done here.

Integrating Newton's equations numerically is how real orbits are computed. But the choice of method matters: the velocity Verlet (leapfrog) method, used in molecular dynamics and celestial mechanics, respects the geometry of mechanics — it is symplectic and conserves angular momentum exactly for central forces — while the explicit Euler method spirals outwards and gains energy. The first experiment checks Kepler's predictions (eccentricity, perihelion, aphelion and period) against a simulation:

"""Kepler's laws from Newton's law, numerically - and why the integrator matters.

We integrate r'' = -r / |r|^3 (units with GM = 1) from r(0) = (1, 0),
v(0) = (0, 1.2). Theory: energy E = -0.28, angular momentum L = 1.2, an
ellipse with eccentricity e = sqrt(1 + 2 E L^2) = 0.44, semi-major axis
a = -1/(2E), perihelion a(1 - e) = 1 and period 2 pi a^(3/2). The velocity
Verlet (leapfrog) method conserves L exactly (up to rounding) and keeps E
within a small bounded band; explicit Euler spirals outwards, gaining energy.

Runs unchanged in the browser (Pyodide) using only the standard library.
"""

import math

def accel(x, y):
    """Inverse-square attraction towards the origin, GM = 1."""
    r3 = (x * x + y * y) ** 1.5
    return -x / r3, -y / r3

def energy(x, y, vx, vy):
    """Energy per unit mass: kinetic plus potential -1/r."""
    return 0.5 * (vx * vx + vy * vy) - 1 / math.hypot(x, y)

def ang_mom(x, y, vx, vy):
    """Angular momentum per unit mass, x vy - y vx."""
    return x * vy - y * vx

def verlet_step(state, h):
    """One velocity Verlet (leapfrog) step."""
    x, y, vx, vy = state
    ax, ay = accel(x, y)
    vx, vy = vx + h / 2 * ax, vy + h / 2 * ay
    x, y = x + h * vx, y + h * vy
    ax, ay = accel(x, y)
    return x, y, vx + h / 2 * ax, vy + h / 2 * ay

def euler_step(state, h):
    """One explicit Euler step."""
    x, y, vx, vy = state
    ax, ay = accel(x, y)
    return x + h * vx, y + h * vy, vx + h * ax, vy + h * ay

def run(step, state, h, n):
    """Apply n steps; return the list of states (including the initial one)."""
    states = [state]
    for _ in range(n):
        state = step(state, h)
        states.append(state)
    return states

def first_return_time(states, h):
    """Time of the first upward crossing of y = 0 with x > 0 after the start."""
    for i in range(1, len(states)):
        (x0, y0, _, _), (x1, y1, _, _) = states[i - 1], states[i]
        if i > 1 and y0 < 0 <= y1 and x1 > 0:
            return (i - 1 + (-y0) / (y1 - y0)) * h
    return None

def kepler_predictions(state):
    """Energy, angular momentum, eccentricity, semi-major axis and period."""
    E = energy(*state)
    L = ang_mom(*state)
    e = math.sqrt(1 + 2 * E * L * L)
    a = -1 / (2 * E)
    return E, L, e, a, 2 * math.pi * a**1.5

def main():
    start = (1.0, 0.0, 0.0, 1.2)
    E, L, e, a, T = kepler_predictions(start)
    print(f"Theory: E = {E:.4f}, L = {L:.4f}, e = {e:.4f}, a = {a:.4f}")
    print(f"        perihelion {a * (1 - e):.4f}, aphelion {a * (1 + e):.4f},"
          f" period {T:.4f}")
    h, n = 0.001, 60000  # about four orbits
    for name, step in [("Verlet", verlet_step), ("Euler", euler_step)]:
        states = run(step, start, h, n)
        radii = [math.hypot(s[0], s[1]) for s in states]
        energies = [energy(*s) for s in states]
        momenta = [ang_mom(*s) for s in states]
        print(f"{name} (h = {h}, {n} steps):")
        e_drift = max(abs(en - E) for en in energies)
        print(f"   max |E(t) - E(0)| = {e_drift:.2e}")
        drift = max(abs(m - L) for m in momenta)
        print(f"   max |L(t) - L(0)| = {drift:.2e}")
        print(f"   r ranges over [{min(radii):.4f}, {max(radii):.4f}]")
        ret = first_return_time(states, h)
        shown = f"{ret:.4f}" if ret else "none"
        print(f"   first return to the x-axis at t = {shown}")
    print()
    print("Epistemic status: numerical integration of one orbit. Verlet keeps L")
    print("constant to rounding error and E in a narrow band, matching Kepler's")
    print("predictions to several digits; Euler drifts outwards. The laws")
    print("themselves are theorems proved in the lesson, not consequences of")
    print("this simulation.")

if __name__ == "__main__":
    main()