Skip to content

ODE Integrators

The numerical propagator integrates the equations of motion forward (or backward) in time using one of two families of solvers:

  1. Adaptive Runge-Kutta (default) — embedded error estimation, PID step-size controller, accept/reject logic. Butcher tableaux from Verner (2010) — the RKV98.IIa.Efficient, RKV87.IIa.Robust and RKV65.IIIXb.Efficient sets published on his delightful web page — plus Tsitouras (2011) for rkts54 and the RODAS4 Rosenbrock method of Hairer & Wanner (1996), §IV.7, for rodas4.
  2. Gauss-Jackson 8 — 8th-order fixed-step multistep predictor-corrector specialised for 2nd-order ODEs (Berry & Healy 2004), the integrator used for the U.S. Space Surveillance Network catalog. For smooth long-duration propagation it needs several times fewer force evaluations than rkv98 at comparable accuracy (satkit benchmarks show 3–10× on GEO and MEO arcs).

Integrator Choices

Select via the integrator parameter of propsettings:

Integrator Order Type Dense Output Notes
rkv98 9(8) adaptive RK, 21 stages (16 + 5 for dense output) 8th-degree Default. Best accuracy for precision work. With enable_interp=False it runs as rkv98_nointerp.
rkv98_nointerp 9(8) adaptive RK, 16 stages None Same stepping accuracy, 24% fewer force evaluations per step; what rkv98 becomes when interpolation is off.
rkv87 8(7) adaptive RK, 17 stages (13 + 4 for dense output) 7th-degree Good balance of speed and accuracy.
rkv65 6(5) adaptive RK, 10 stages 6th-degree Faster, moderate accuracy.
rkts54 5(4) adaptive RK, 7 stages (FSAL) 4th-degree Fastest. Good for quick propagations.
rodas4 4(3) Rosenbrock, 6 stages None L-stable (implicit). For stiff problems. No STM support.
gauss_jackson8 8 fixed-step multistep 5th-order Hermite High-efficiency for smooth long-duration propagation. No STM support.

Higher-order RK methods can take larger steps for the same accuracy, so despite more stages per step they often require fewer total function evaluations. The default rkv98 covers almost all situations. Reach for rodas4 on stiff problems (very low perigee, re-entry); reach for gauss_jackson8 on long smooth arcs.

import satkit as sk

# Faster adaptive RK
settings = sk.propsettings(integrator=sk.integrator.rkts54)

# 8(7) integrator with a higher-degree gravity model
settings = sk.propsettings(
    integrator=sk.integrator.rkv87,
    gravity_model=sk.gravmodel.egm96,
    gravity_degree=16,
)

# Implicit Rosenbrock for stiff low-perigee dynamics
settings = sk.propsettings(
    integrator=sk.integrator.rodas4,
    gravity_degree=8,
)

# Gauss-Jackson 8 for long-duration GEO with a 120 s fixed step
settings = sk.propsettings(
    integrator=sk.integrator.gauss_jackson8,
    gj_step_seconds=120.0,
)

Note

The rodas4 and gauss_jackson8 integrators do not support state transition matrix propagation (output_phi=True). Attempting to use output_phi=True with either will raise a RuntimeError.

Integrator step budget: max_steps

Every integrator stops with a max-steps error once it exceeds propsettings.max_steps total steps. This is a runaway-propagation safeguard, not a quality knob. The default of 1_000_000 comfortably covers the longest realistic arcs — roughly 700 days of gauss_jackson8 at a 60 s step, or millions of adaptive RK steps at typical tolerances. Lower it if you want to fail-fast on configurations that would take a very long time; raise it if you hit the limit on a genuine long-arc propagation.

Gauss-Jackson step-size selection

gauss_jackson8 uses a fixed step size (gj_step_seconds) which the user must choose based on the orbit regime. Typical values (satkit's guidance; Berry & Healy 2004 discuss the accuracy-vs-step trade-off):

  • LEO (400-800 km): 30-60 s
  • MEO: 60-300 s
  • GEO: 300-600 s
  • HEO / eccentric transfer: use rkv98 instead — GJ8's fixed step wastes accuracy at apogee and misses resolution at perigee.

GJ8 is also unsuitable for propagation across discontinuities such as eclipse boundaries or impulsive maneuvers — use an adaptive RK method for those cases. The integrator needs ≥ 9 steps of startup, so propagations shorter than ~9 × gj_step_seconds will fail; use an adaptive RK integrator for such short intervals.

Error Tolerances

The adaptive RK integrators (and rodas4) accept the same step every time both:

  • the absolute error estimate falls below abs_error (default 1e-8)
  • the relative error estimate falls below rel_error (default 1e-8)

For sub-meter precision over a day, tighten both to 1e-10 to 1e-13. For coarse mission planning, 1e-6 to 1e-8 is usually fine.

gauss_jackson8 ignores both — its accuracy is set by the fixed step gj_step_seconds.

Starting Step and Warm Start

An adaptive integrator has to guess its first step. The textbook heuristic (Hairer–Nørsett–Wanner) is deliberately conservative and scale sensitive: for an orbit in metres and seconds at 1e-9 tolerance it starts rkv98 at a fraction of a millisecond against a working step of a few hundred seconds, and the controller spends a dozen accepted steps growing into it — about half the force evaluations of a one-hour arc. satkit therefore does not use it. By default the first step is derived from the initial state, the tolerances and the integrator order:

h0 = 1.5 · |r| / |v| · tol^(1/(p+1)),    tol = rel_error + abs_error / |r|

The settled stride of an order-p method scales as tol^(1/(p+1)), and the constant is fit to LEO strides. Across rkts54 to rkv98 and tolerances from 1e-6 to 1e-12 this lands within a factor of about 2.5 of the settled stride (for rkv98 at 1e-9 in LEO: 170 s predicted, 270 s settled), which the controller closes within a step or two.

Two settings give you control over this:

  • propsettings.initial_step_secs overrides the default. It is a magnitude (backward propagation applies the sign) and is clamped to the arc length. None restores the state-derived default. Ignored by gauss_jackson8.
  • propresult.next_step_secs is the step the integrator would take next — its working stride at the end of the arc, not the final step that was shortened to land exactly on the end time. Feed it to the next arc to continue at full stride:
ps = sk.propsettings(abs_error=1e-9, rel_error=1e-9)
res = sk.propagate(state, t0, duration_secs=3600.0, propsettings=ps)
ps.initial_step_secs = res.next_step_secs          # warm start
res = sk.propagate(res.state, res.time_end, duration_secs=3600.0, propsettings=ps)

Chained arcs then cost the same as one continuous propagation. The step controller and the error tolerances are unchanged; only the starting point is.

See Also