Lecture 03: Kepler Orbits with REBOUND — Lecture Notes

Terms Developed in Context

The terms below connect the mathematical model, implementation choices, computed diagnostics, and limits of the scientific claim.

Symplectic integrator

A time-stepping map that preserves Hamiltonian phase-space geometry. In this lecture it is used to address What does a symplectic integrator preserve, and what does it not promise?

Model connection: Read the term against Kepler drift and the governing relation \(H=\sum_i\left(\frac{p_i^2}{2m_i}-\frac{GMm_i}{r_i}\right)+H_{\rm interaction}\). The relevant computational evidence is a controlled output, residual, invariant, or refinement sequence rather than visual smoothness alone.

Interpretive limit: An ideal two-body propagation is not an ephemeris. Earth-center state, Sun motion, Moon, other planets, and relativistic terms matter when comparing against DE441. This distinction separates numerical behavior from physical accuracy and model adequacy.

Kepler drift

Exact or specialized propagation under the dominant two-body Hamiltonian during a split step. In this lecture it is used to address What does a symplectic integrator preserve, and what does it not promise?

Model connection: Read the term against Kick step and the governing relation \(H=\sum_i\left(\frac{p_i^2}{2m_i}-\frac{GMm_i}{r_i}\right)+H_{\rm interaction}\). The relevant computational evidence is a controlled output, residual, invariant, or refinement sequence rather than visual smoothness alone.

Interpretive limit: An ideal two-body propagation is not an ephemeris. Earth-center state, Sun motion, Moon, other planets, and relativistic terms matter when comparing against DE441. This distinction separates numerical behavior from physical accuracy and model adequacy.

Kick step

Momentum update produced by the perturbing interaction part of a split Hamiltonian. In this lecture it is used to address What does a symplectic integrator preserve, and what does it not promise?

Model connection: Read the term against Timestep convergence and the governing relation \(H=\sum_i\left(\frac{p_i^2}{2m_i}-\frac{GMm_i}{r_i}\right)+H_{\rm interaction}\). The relevant computational evidence is a controlled output, residual, invariant, or refinement sequence rather than visual smoothness alone.

Interpretive limit: An ideal two-body propagation is not an ephemeris. Earth-center state, Sun motion, Moon, other planets, and relativistic terms matter when comparing against DE441. This distinction separates numerical behavior from physical accuracy and model adequacy.

Timestep convergence

Evidence obtained by comparing results from systematically refined step sizes. In this lecture it is used to address What does a symplectic integrator preserve, and what does it not promise?

Model connection: Read the term against Ephemeris and the governing relation \(H=\sum_i\left(\frac{p_i^2}{2m_i}-\frac{GMm_i}{r_i}\right)+H_{\rm interaction}\). The relevant computational evidence is a controlled output, residual, invariant, or refinement sequence rather than visual smoothness alone.

Interpretive limit: An ideal two-body propagation is not an ephemeris. Earth-center state, Sun motion, Moon, other planets, and relativistic terms matter when comparing against DE441. This distinction separates numerical behavior from physical accuracy and model adequacy.

Ephemeris

A time-tagged prediction of an object's position and velocity in a specified frame and dynamical model. In this lecture it is used to address What does a symplectic integrator preserve, and what does it not promise?

Model connection: Read the term against Symplectic integrator and the governing relation \(H=\sum_i\left(\frac{p_i^2}{2m_i}-\frac{GMm_i}{r_i}\right)+H_{\rm interaction}\). The relevant computational evidence is a controlled output, residual, invariant, or refinement sequence rather than visual smoothness alone.

Interpretive limit: An ideal two-body propagation is not an ephemeris. Earth-center state, Sun motion, Moon, other planets, and relativistic terms matter when comparing against DE441. This distinction separates numerical behavior from physical accuracy and model adequacy.

Problem and physical meaning

For a dominant central mass and weak perturbations, Wisdom-Holman splitting separates Keplerian drift from perturbing kicks. WHFast is fixed-step and symplectic; IAS15 is adaptive and high-order. Their errors have different structure and should be selected by the problem, not by a single accuracy slogan.

What does a symplectic integrator preserve, and what does it not promise? Begin from the scientific observable or state variable, and name what is directly specified versus inferred. The model should be simple enough to validate, but not quietly presented as the entire astrophysical system.

Model and derivation

The governing expression for this lecture is:

Define the symbols and the independent variable before numerical evaluation. Dimensional analysis checks that both sides have matching units; if the model is dimensionless, state the adopted scale factors and the transformation back to physical units.

For a dominant central mass and weak perturbations, Wisdom-Holman splitting separates Keplerian drift from perturbing kicks. WHFast is fixed-step and symplectic; IAS15 is adaptive and high-order. Their errors have different structure and should be selected by the problem, not by a single accuracy slogan.

Numerical implementation

Initialize units before particles; validate the two-body orbit; sweep fixed timesteps; record energy and angular momentum; compare a one-day Sun-Earth propagation with the JPL state in a consistent heliocentric frame.

Keep physical quantities and numerical arrays conceptually separate. Convert units at well-defined boundaries, preserve raw inputs, and record dtype, solver, tolerances/timestep, grid, and software versions. A numerical method is part of the model evidence, not an invisible implementation detail.

Worked, programmatically computed example

Using Astropy GM_sun, the circular 1-AU Kepler period is 365.25690 d. After one day, REBOUND and independent SciPy two-body runs differ from the DE441 Earth-center vector by 126.0 km and 126.1 km respectively; the residual is model/frame scope, not an accuracy claim.

These values are generated in src/astr474_computations.py, then shared with lecture/lab/assignment renderers. Re-run the generator after changing an input; do not hand-edit a displayed result without recomputing dependent artifacts.

At least one independent check is required: an analytic solution, a second numerical method, a conservation law, or a controlled convergence sequence. When these disagree, diagnose units, frame, discretization, source terms, and model scope before selecting a result.

Error budget and limitations

An ideal two-body propagation is not an ephemeris. Earth-center state, Sun motion, Moon, other planets, and relativistic terms matter when comparing against DE441.

Classify each uncertainty before combining terms: local/global solver error, timestep or spatial discretization, Monte Carlo sampling, parameter uncertainty, measurement uncertainty, and model discrepancy are not interchangeable. Conservation can test internal consistency but does not prove the model is physically complete.

Connection to evidence and research practice

Lab 02 uses the level-2 JPL vector record and distinguishes integrator error from model discrepancy.

For real parameters, the source record distinguishes the fetched product and epoch from derived values. For synthetic experiments, the initial conditions, generator, seed, and exact file are recorded. A simulation is not an observation, even when its parameters resemble a real object.

Misconceptions to avoid

• More displayed digits do not imply more physical accuracy.
• A successful solver exit or smooth plot does not establish convergence.
• Agreement between two solvers with the same wrong model is not independent physical validation.
• A formal covariance or Monte Carlo interval is conditional on the adopted model and noise process.

Study and review questions

1. What does a symplectic integrator preserve, and what does it not promise?
2. Which numerical parameter would you vary first, what observable would you track, and how would you know the sequence reached an asymptotic regime?
3. Design one independent check and state what failure would change your scientific interpretation.

Preparation for the next lecture: explain how this method's assumptions change when the problem becomes more coupled, stochastic, or data-driven.

Lecture-specific derivation and assumptions

For a circular orbit, equate centripetal acceleration v^2/r to GM/r^2, obtaining v^2=GM/r. With v=2*pi*r/P, rearrangement gives P^2=4*pi^2*r^3/(GM). For a two-body system use relative separation and G(M+m); the test-particle approximation drops m only when m/M is negligible. Real ephemeris vectors require matching target, center, epoch scale, reference plane, correction convention, and units before a residual has meaning.

Algorithm trace and computed-output interpretation

Convert Astropy GM_sun to AU^3/day^2, evaluate the circular 1-AU period, and compare with the 365.25-day calendar convention. For the supplied level-2 Horizons vectors, initialize the same state in REBOUND and an independent DOP853 integration, propagate one day, and subtract the second vector in the stated frame. Report both residuals and separate numerical convergence from missing Moon/planetary perturbations and frame/model mismatch.

Computed result: Using Astropy GM_sun, the circular 1-AU Kepler period is 365.25690 d. After one day, REBOUND and independent SciPy two-body runs differ from the DE441 Earth-center vector by 126.0 km and 126.1 km respectively; the residual is model/frame scope, not an accuracy claim.

The two-body period is an analytic benchmark, not a precise modern Earth-year prediction. The embedded vectors are traceable published ephemeris values but were not live-verified in this session; their request parameters and CSV checksum are recorded in the data provenance. Do not call a model-to-ephemeris residual an integrator error until the initial state and comparison frame are equivalent.

Two-meeting teaching arc

Each weekly module is taught in two 75-minute meetings. The sequence below includes a formative check in each meeting and is intended to use the lecture's worked computation, not merely present it.

1. Meeting 1 (75 min): 10 min state-vector context; 20 min Newtonian two-body reduction; 20 min Kepler-period derivation and units; 20 min Astropy constant check; 5 min frame-convention exit ticket.
2. Meeting 2 (75 min): 10 min recap; 20 min REBOUND WHFast/IAS15 assumptions; 25 min reproduce the one-day state comparison; 15 min classify residual sources; 5 min recommend an integrator for the next module.

Targeted study questions

1. Which exact center/frame transformation is required before comparing heliocentric and barycentric states?
2. How does using G(M+m) change the one-AU circular period for finite m?
3. Which residual trend under timestep refinement would distinguish solver error from model discrepancy?