ODE Solver Method Comparison — Kepler Orbit¶
This example compares ivprk (RK45), ivpbs (DOP853), ivpam (LSODA),
ode45, and ivp_dense on the Kepler two-body orbit problem.
The Kepler problem is a conservative Hamiltonian system with a known closed orbit over one period \(t \in [0, 2\pi]\), making it an ideal benchmark for comparing energy conservation and position accuracy across methods:
\[\mathbf{r}'' = -\frac{\mathbf{r}}{|\mathbf{r}|^3}\]
written as the first-order system \([x,\, y,\, v_x,\, v_y]' = [v_x,\, v_y,\, -x/r^3,\, -y/r^3]\) with eccentricity \(e = 0.1\).
Example Code¶
"""IMSL ODE method comparison example: Kepler orbit (two-body problem).
Demonstrates and compares ``ivprk`` (RK45), ``ivpbs`` (DOP853), ``ivpam``
(LSODA), ``ode45``, and ``ivp_dense`` on the Kepler two-body orbit problem:
r'' = -r / |r|^3 → [x, y, vx, vy]' = [vx, vy, -x/r^3, -y/r^3]
Initial conditions (eccentricity e=0.1):
x(0) = 1 - e = 0.9, y(0) = 0
vx(0) = 0, vy(0) = sqrt((1+e)/(1-e))
Integration over one full orbit: t ∈ [0, 2π].
Outputs:
- Table of final-position error and energy conservation per method, printed to
stdout.
- SVG plot saved to ``test_output/example_imsl_ode_methods.svg``.
"""
from __future__ import annotations
from pathlib import Path
from typing import Dict
import matplotlib.pyplot as plt
import numpy as np
from differentialequations import ivpam, ivpbs, ivp_dense, ivprk, ode45
# ---------------------------------------------------------------------------
# Kepler problem definition
# ---------------------------------------------------------------------------
ECCENTRICITY = 0.1
X0 = 1.0 - ECCENTRICITY
VY0 = np.sqrt((1.0 + ECCENTRICITY) / (1.0 - ECCENTRICITY))
Y0_KEPLER = [X0, 0.0, 0.0, VY0]
T_SPAN = (0.0, 2.0 * np.pi)
T_EVAL_DENSE = np.linspace(0.0, 2.0 * np.pi, 500)
def kepler(t: float, y: np.ndarray) -> list[float]:
"""Kepler two-body ODE right-hand side.
Args:
t (float): Time (unused — autonomous system).
y (np.ndarray): State ``[x, y, vx, vy]``.
Returns:
list[float]: Derivatives ``[vx, vy, ax, ay]``.
"""
x, yc, vx, vy = y
r3 = (x**2 + yc**2) ** 1.5
return [vx, vy, -x / r3, -yc / r3]
def kepler_energy(y: np.ndarray) -> np.ndarray:
"""Compute Kepler total mechanical energy per unit mass.
E = 0.5*(vx^2 + vy^2) - 1/r
Args:
y (np.ndarray): State array of shape ``(4, N)``.
Returns:
np.ndarray: Energy at each column, shape ``(N,)``.
"""
x, yc, vx, vy = y[0], y[1], y[2], y[3]
r = np.sqrt(x**2 + yc**2)
return 0.5 * (vx**2 + vy**2) - 1.0 / r
def run_demo_imsl_ode_methods() -> Dict[str, object]:
"""Run IMSL ODE method comparison on the Kepler orbit problem.
Solves one full orbit using ``ivprk``, ``ivpbs``, ``ivpam``, ``ode45``,
and ``ivp_dense``. Reports final-position error relative to the known
initial position (closed orbit) and energy drift.
Args:
None
Returns:
Dict[str, object]: Keys ``methods`` (list of names), ``errors``
(ndarray of final-position errors), ``energy_drift`` (ndarray),
and ``plot_path`` (str).
"""
solvers = {
"ivprk (RK45)": lambda: ivprk(kepler, T_SPAN, Y0_KEPLER),
"ivpbs (DOP853)": lambda: ivpbs(kepler, T_SPAN, Y0_KEPLER),
"ivpam (LSODA)": lambda: ivpam(kepler, T_SPAN, Y0_KEPLER),
"ode45": lambda: ode45(kepler, T_SPAN, Y0_KEPLER),
"ivp_dense": lambda: ivp_dense(kepler, T_SPAN, Y0_KEPLER, T_EVAL_DENSE),
}
E0 = kepler_energy(np.array(Y0_KEPLER).reshape(4, 1))[0]
results = {}
method_names = []
pos_errors = []
energy_drifts = []
print("\nIMSL ODE Method Comparison: Kepler Two-Body Orbit (e=0.1)")
print("=" * 72)
print(f"{'Method':<22} {'Steps':>6} {'NFev':>6} {'|r_final-r0|':>14} {'|ΔE/E0|':>12}")
print("-" * 72)
for name, solver_fn in solvers.items():
res = solver_fn()
results[name] = res
# Position error: orbit is closed after one period
xf, yf = res.y[0, -1], res.y[1, -1]
pos_err = np.sqrt((xf - X0) ** 2 + yf**2)
# Energy conservation
E_final = kepler_energy(res.y[:, -1].reshape(4, 1))[0]
energy_drift = abs((E_final - E0) / E0)
print(
f"{name:<22} {res.n_steps:>6d} {res.n_eval:>6d} "
f"{pos_err:>14.2e} {energy_drift:>12.2e}"
)
method_names.append(name)
pos_errors.append(pos_err)
energy_drifts.append(energy_drift)
print("=" * 72)
print(f"Reference energy E0 = {E0:.6f}")
# -----------------------------------------------------------------------
# Plot
# -----------------------------------------------------------------------
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "example_imsl_ode_methods.svg"
colors = ["#0369a1", "#dc2626", "#16a34a", "#d97706", "#7c3aed"]
styles = ["-", "--", "-.", ":", (0, (3, 1, 1, 1))]
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# Panel 1: orbital trajectories
ax = axes[0]
# Plot the dense result as the reference trajectory
dense_res = results["ivp_dense"]
ax.fill(
dense_res.y[0],
dense_res.y[1],
alpha=0.07,
color="#0369a1",
label=None,
)
for (name, res), color, ls in zip(results.items(), colors, styles):
ax.plot(res.y[0], res.y[1], color=color, linestyle=ls, linewidth=1.5, label=name)
ax.plot(0.0, 0.0, "k*", markersize=10, label="Focus (origin)")
ax.plot(X0, 0.0, "ko", markersize=5, label="Start")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("Kepler Orbit: All Methods Overlaid")
ax.set_aspect("equal")
ax.legend(fontsize=7, loc="upper right")
ax.grid(True, alpha=0.3)
# Panel 2: energy error over time (ivp_dense gives smooth curve)
ax = axes[1]
E_dense = kepler_energy(dense_res.y)
ax.plot(
dense_res.t,
np.abs((E_dense - E0) / E0),
color="#0369a1",
linewidth=1.5,
label="ivp_dense (RK45)",
)
ax.set_xlabel("t")
ax.set_ylabel(r"$|\Delta E / E_0|$")
ax.set_title("Energy Conservation (ivp_dense)")
ax.set_yscale("log")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
# Panel 3: final position error and energy drift bar chart
ax = axes[2]
x_pos = np.arange(len(method_names))
width = 0.35
bars1 = ax.bar(x_pos - width / 2, pos_errors, width, label="Position error", color="#0369a1", alpha=0.85)
bars2 = ax.bar(x_pos + width / 2, energy_drifts, width, label=r"$|\Delta E/E_0|$", color="#dc2626", alpha=0.85)
ax.set_yscale("log")
ax.set_xticks(x_pos)
ax.set_xticklabels([n.split(" ")[0] for n in method_names], rotation=20, ha="right", fontsize=8)
ax.set_ylabel("Error")
ax.set_title("Final Errors by Method")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3, axis="y")
fig.suptitle(
"IMSL ODE Methods: Kepler Two-Body Orbit (e=0.1, t∈[0,2π])",
fontsize=12,
fontweight="bold",
)
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
print(f"\nPlot saved to: {plot_path}")
return {
"methods": method_names,
"errors": np.array(pos_errors),
"energy_drift": np.array(energy_drifts),
"plot_path": str(plot_path),
}
if __name__ == "__main__":
run_demo_imsl_ode_methods()
Input (Console)¶
Run the script from the package root:
python examples/example_imsl_ode_methods.py
Console Output¶
IMSL ODE Method Comparison: Kepler Two-Body Orbit (e=0.1)
========================================================================
Method Steps NFev |r_final-r0| |ΔE/E0|
------------------------------------------------------------------------
ivprk (RK45) 33 224 1.57e-05 2.46e-06
ivpbs (DOP853) 13 158 5.32e-06 6.22e-07
ivpam (LSODA) 93 201 5.19e-05 7.05e-06
ode45 33 224 1.57e-05 2.46e-06
ivp_dense 500 74 1.57e-01 3.62e-02
========================================================================
Reference energy E0 = -0.500000
Plot saved to: test_output\example_imsl_ode_methods.svg
Plot Output¶
Generated SVG plot:
Left: orbital trajectories from all five methods overlaid on the Kepler ellipse (eccentricity 0.1). Centre: energy conservation error over time for the dense-output run. Right: bar chart comparing final position error and relative energy drift per method.¶
Description¶
ivprk and ode45 both invoke RK45 and produce identical results.
ivpbs (DOP853) achieves smaller error with fewer function evaluations due
to its higher-order scheme. ivpam (LSODA) uses an Adams
predictor-corrector and is efficient for mildly stiff regimes. ivp_dense
evaluates the solution at 500 uniformly-spaced points specified via
t_eval; note that its default (looser) tolerances produce a smoother
trajectory at the cost of slightly higher position and energy errors compared
with the tolerance-tightened adaptive solvers.