Kepler integrators, Adams-type multistep methods and global trajectory optimisation, in Python.
A comparison of numerical integrators applied to the Kepler problem, rewritten in Python from the MATLAB programs developed for my M.Sc. thesis in Mathematics on symplectic integrators, and extended with the Adams-type multistep methods of a numerical analysis course.
The MATLAB versions are in numerical-methods-matlab.
Implemented methods:
- Explicit Euler (order 1)
- Symplectic Euler (order 1)
- Implicit midpoint (order 2), solved with
scipy.optimize.fsolve - Störmer–Verlet (order 2)
- Adams–Bashforth, 1 to 4 steps (order 1 to 4)
- Adams–Moulton, 1 to 4 steps (order 2 to 5), implicit equation solved by fixed-point iteration
- ABM4 predictor–corrector in PECE form (order 5)
Kepler problem in normalised units (mu = 1, semi-major axis a = 1) with eccentricity e = 0.7, starting from pericentre. The state is y = (x, y, vx, vy). The conserved quantities are the energy E = -0.5 and the angular momentum L = sqrt(1 - e^2). The reference is the exact solution of Kepler's equation, solved with Newton's method (kepler_position in convergence.py). DOP853 (tolerances 1e-13) provides the starting values of the multistep methods.
kepler.py: force, initial state, energy, angular momentum, the four integrators,integrateand the reference solution.convergence.py: convergence study. Computes the global position error at t = 3 against the exact solution for several step sizes, estimates the observed order with a log-log fit and savesfigures/convergence.png.orbits.py: orbit comparison for e = 0.7; savesfigures/orbits.png.multistep.py: Adams–Bashforth, Adams–Moulton and ABM4 (PECE) multistep methods.multistep_convergence.py: observed orders of the multistep methods on the Kepler problem.multistep_orbits.py: Kepler orbits computed with the multistep methods.gauss_jackson.py: Gauss–Jackson (Störmer–Cowell) integrator for second-order problems.local_optimisation.py: Newton and gradient methods (steepest descent, conjugate gradient, least squares) for unconstrained minimisation, tested on the Rosenbrock functiontrajectory_optimisation.py: global optimisation of an Earth–Mars or Earth–Venus transfer (Lambert solver and differential evolution), with a circular coplanar model and a model with real ephemerides (astropy); saves pork-chop plots infigures/.compare_local_global.py: comparison of a local method (L-BFGS-B from random starting points) with differential evolution on the Earth–Mars costparallel_speedup.py: strong-scaling test of the search with several worker processes; savesfigures/speedup.png.tests/test_kepler.py,tests/test_multistep.py,tests/test_gauss_jackson.py,tests/test_trajectory_optimisation.py,tests/test_local_optimisation.py,tests/test_compare_local_global.py,tests/test_parallel_speedup.py: automated tests with pytest.requirements.txt: required libraries.
python3 -m venv .venv
source .venv/bin/activate
pip install -r requirements.txt
python3 convergence.py
python3 orbits.py
python3 trajectory_optimisation.py mars
python3 trajectory_optimisation.py venus --start 2026-10-01 --end 2028-12-01
python3 trajectory_optimisation.py mars --model circular
pytest -v
Global position error at t = 3 for the Kepler problem (e = 0.7, pericentre start), measured against the exact solution of Kepler's equation. All four methods use the same step counts (n = 250 … 4000).
| Method | Observed order |
|---|---|
| Explicit Euler | 0.93 |
| Symplectic Euler | 0.95 |
| Implicit midpoint | 2.00 |
| Störmer–Verlet | 2.00 |
The final time is t = 3 rather than one full period: at t = 2π the symmetry of the orbit hides the first-order error of symplectic Euler and the fitted order would read about 2.
Cross-validation. The errors agree with the MATLAB implementation (numerical-methods-matlab) to the four significant digits compared (20 values: 4 methods × 5 step counts).
Orbits obtained with e = 0.7, 5 orbits and 200 steps per orbit, compared with the exact ellipse (dashed). All start from pericentre.
With this relatively coarse step:
- Explicit Euler diverges: the energy grows at every step and the orbit opens up until it leaves the ellipse.
- Symplectic Euler stays bounded, but the numerical orbit is distorted with respect to the exact ellipse.
- Implicit midpoint keeps the shape and size of the ellipse, but the major axis rotates at every revolution (a phase error, absent in the exact problem).
- Störmer–Verlet is the closest to the exact orbit, with a very small precession.
Symplectic and symmetric methods do not show the systematic energy drift of explicit Euler: the dominant error is a phase error, and the energy stays bounded.
multistep.py adds Adams-type linear multistep methods: Adams–Bashforth (1 to 4 steps, order 1 to 4), Adams–Moulton (1 to 4 steps, order 2 to 5, implicit equation solved by fixed-point iteration) and the predictor–corrector scheme ABM4 in PECE form (AB4 predictor, AM4 corrector). Starting values are computed with DOP853.
The observed orders on the test problem y' = -y, y(0) = 1 match the theoretical ones (checked by tests/test_multistep.py). On the Kepler problem (e = 0.7, position error at t = 3 against the exact solution of Kepler's equation, n = 250 … 4000 steps, multistep_convergence.py) the observed orders are:
| Method | AB1 | AB2 | AB3 | AB4 | AM1 | AM2 | AM3 | AM4 | ABM4 (PECE) |
|---|---|---|---|---|---|---|---|---|---|
| Observed order | 0.93 | 1.99 | 2.88 | 3.95 | 2.00 | 2.91 | 3.97 | 5.31 | 5.00 |
AB1 coincides with explicit Euler. AM4 reads slightly above 5 because its finest errors (about 4e-12) are close to the accuracy of the starting values. With 1,000 steps, ABM4 (PECE) reaches an error of about 5e-8, against about 1.3e-5 for AB4. The errors agree with the MATLAB implementation (numerical-methods-matlab): 41 of the 45 values (9 methods × 5 step counts) match to four significant digits; the other four are at the finest steps, and AM4 at n = 4000 differs by about 3%.
Orbits computed with a deliberately coarse step (2 periods, 1,000 steps, multistep_orbits.py) show the differences in accuracy. The position error at the final time, against the exact solution, ranges from 3.4 for AB1 and 1.1 for AB2, whose orbits lose their shape, to 4.9e-3 for AM3, 1.8e-3 for ABM4 (PECE) and 1.2e-4 for AM4, which follow the exact orbit closely. At equal order, the implicit Adams–Moulton methods are more accurate than the explicit Adams–Bashforth ones.
These methods were first implemented in MATLAB for the Scientific Computing course (2023/24) and are rewritten here from the formulas.
gauss_jackson.py is a fixed-step multistep method for second-order problems y'' = f(t, y), in backward-difference form: Störmer predictor, Cowell corrector (one PECE correction, default order 8), with velocities from Adams–Bashforth–Moulton. The coefficients are generated in exact rational arithmetic from their generating functions, and the starting values come from DOP853 (rtol 1e-13).
from gauss_jackson import gauss_jackson
t, Y, V = gauss_jackson(f, t0, y0, v0, h, n_steps, order=8)Checked by tests/test_gauss_jackson.py:
- the coefficients equal the tabulated Störmer, Cowell, Adams–Bashforth and Adams–Moulton values;
- polynomial solutions are reproduced exactly;
- Kepler orbit, e = 0.7, against the exact solution (800 → 1600 steps): observed orders 6.2 for order 4 and 7.7 for order 6, i.e. about order + 2 for the position error; for order 8 the error reaches about 1e-11, the floor set by the starting values and round-off, so only the first refinement is meaningful (observed 8.5);
- 10 periods at 400 steps per period: maximum energy error 9e-7, against 1.4e-5 for classical RK4;
- J2 perturbation: the nodal precession rate matches -1.5 J2 n (R/p)^2 cos i within 2%.
This is the difference form, not the summed-ordinate form of Berry and Healy. Reference: D. A. Vallado, Fundamentals of Astrodynamics and Applications, Sec. 7.5.
trajectory_optimisation.py searches the departure date and time of flight of an
Earth–Mars or Earth–Venus transfer with scipy.optimize.differential_evolution.
Each candidate transfer is a Lambert problem solved with universal variables.
The cost is the sum of the departure and arrival excess speeds. Two models are available:
- Circular model: circular coplanar planets, two-body dynamics. The minimum found (Mars 5.594 km/s, Venus 5.202 km/s) matches the Hohmann transfer, which validates the Lambert solver and the search.
- Ephemeris model (default): real heliocentric positions and velocities from the astropy built-in ephemeris (cubic Hermite interpolation), 3D Lambert problem. In the departure window 2026-10-01 to 2028-12-01 the search finds, for Mars, departure 2026-10-31, 310 days, arrival 2027-09-06, 5.61 km/s; for Venus, departure 2028-01-10, 195 days, arrival 2028-07-24, 8.11 km/s.
The destination and the model are chosen from the command line (see Usage). The model uses patched conics (no gravity assists, no parking orbit) and the results are not compared with mission data. The astropy built-in ephemeris is less accurate than the JPL ones.
The parallel and local-versus-global experiments below use the circular Earth–Mars cost.
parallel_speedup.py measures the strong scaling of the search with an increasing
number of worker processes. To mimic a more expensive model, the cost function is
evaluated 20 times per call (REPS = 20, artificial load). On a MacBook Pro with
10 cores the speedup is 1.89, 3.67 and 5.93 with 2, 4 and 8 processes (efficiency
0.94, 0.92 and 0.74). The scaling is consistent with Amdahl's law with an effective
parallel fraction of 0.95, which includes the cost of starting the processes and of
moving data.
References: Curtis, Orbital Mechanics for Engineering Students (Lambert's problem); Storn and Price, J. Global Optimization 11 (1997), 341-359; Virtanen et al., Nature Methods 17 (2020), 261-272.
compare_local_global.py runs a local method (scipy.optimize.minimize,
L-BFGS-B with finite-difference gradients) from 100 random starting points on
the Earth-Mars cost, and compares it with differential evolution.
| Method | Reaches the 5.594 km/s minimum | Evaluations |
|---|---|---|
| L-BFGS-B, 100 random starts | 57 of 100 | about 80 per run |
| Differential evolution | always (4 seeds) | 6000-12000 per run |
The local method is much cheaper but depends on the starting point; the worst
run ends at 10.47 km/s. Differential evolution costs more evaluations and finds
the minimum from any seed. The local method here is SciPy's, not the one from
local_optimisation.py, because the cost function has no analytic gradient or Hessian.
local_optimisation.py is a Python version of the methods from the
Scientific Computing labs 7-9 (originally written in MATLAB):
Newton, Newton with Armijo backtracking, steepest descent with Armijo
backtracking, steepest descent with exact line search, conjugate gradient
(the last two for symmetric positive definite linear systems) and conjugate
gradient for least squares problems.
On the Rosenbrock function from (-1.2, 1), tolerance 1e-6 on the gradient norm:
| Method | Iterations |
|---|---|
| Newton | 6 |
| Newton + Armijo, step never reset | 174 |
| Newton + Armijo, step reset to 1 at every iteration | 21 |
| Steepest descent + Armijo | 32847 |
A local method needs a good starting point; the global search above (differential evolution) does not.
The repository has automated tests (pytest): orders of convergence of the integrators, multistep methods, the Gauss–Jackson integrator, the Lambert problem and global optimisation (circular and real-ephemeris models), the parallel speedup experiment, the local optimisation methods and the local versus global comparison. They run on every push with GitHub Actions.
- Comparison with the MATLAB implementations beyond the convergence studies (long-time energy and angular momentum).
- Study of other eccentricities and longer integration times.
- Planar circular restricted three-body problem.