Tutorial 05: Comparison: Cowell vs Encke

This tutorial compares absolute numerical integration (Cowell’s Method) and perturbative numerical integration (Encke’s Method) within OpenGNC.

Theory Prerequisites

1. Cowell’s Method

Integrates the absolute equations of motion directly. The state integrated is the absolute position \(\mathbf{r}\) and velocity \(\mathbf{v}\):

\[ \ddot{\mathbf{r}} = -\frac{\mu}{r^3}\mathbf{r} + \mathbf{a}_{pert} \]

2. Encke’s Method

Useful for orbits where perturbations are small. It integrates the deviation \(\delta \mathbf{r}\) from a reference Keplerian orbit \(\mathbf{r}_{ref}\):

\[ \mathbf{r} = \mathbf{r}_{ref} + \delta \mathbf{r} \]

\[ \ddot{\delta \mathbf{r}} = \frac{\mu}{r_{ref}^3} \left( f(q)\mathbf{r}_{ref} - \delta\mathbf{r} \right) + \mathbf{a}_{pert} \]

where \(q = \frac{\delta\mathbf{r} \cdot (\delta\mathbf{r} - 2\mathbf{r}_{ref})}{r_{ref}^2}\). Encke’s method reduces numerical roundoff errors for bound orbits with small disturbances.

import numpy as np
import time
from opengnc.propagators.cowell import CowellPropagator
from opengnc.propagators.encke import EnckePropagator
from opengnc.disturbances.gravity import J2Gravity, TwoBodyGravity
from opengnc.integrators.rk4 import RK4

print("Propagator modules loaded.")

Propagator modules loaded.

1. Setup Simulation Environment

We define a LEO orbit and introduce a \(J_2\) disturbance to trigger the perturbative integrator loops on both solvers.

# Initial Orbit (m, m/s)
r_0 = np.array([7000000.0, 0.0, 0.0])
v_0 = np.array([0.0, 7546.0 * np.cos(np.radians(28)), 7546.0 * np.sin(np.radians(28))]) 

duration = 10000.0 # [s] - Approx 2 orbits
dt_step = 10.0 # Steps size [s]

def j2_perturbation(t, r, v):
    return J2Gravity().get_acceleration(r) - TwoBodyGravity().get_acceleration(r)

2. Cowell Propagation Benchmarking

cowell = CowellPropagator()

t0 = time.perf_counter()
r_f_cowell, v_f_cowell = cowell.propagate(r_0, v_0, duration, perturbation_acc_fn=j2_perturbation, dt_step=dt_step)
t_cowell = time.perf_counter() - t0

print(f"Cowell Final Pos:  {r_f_cowell} m")
print(f"Cowell Execution Time: {t_cowell:.6f} seconds")

Cowell Final Pos:  [-1310807.66896629 -6060113.86490409 -3226579.61689067] m
Cowell Execution Time: 0.096331 seconds

3. Encke Propagation Benchmarking

encke = EnckePropagator(rect_tol=1e-6) # Triggers rectification when deviation gets large

t0 = time.perf_counter()
r_f_encke, v_f_encke = encke.propagate(r_0, v_0, duration, perturbation_acc_fn=j2_perturbation, dt_step=dt_step)
t_encke = time.perf_counter() - t0

print(f"Encke Final Pos:  {r_f_encke} m")
print(f"Encke Execution Time: {t_encke:.6f} seconds")

Encke Final Pos:  [-1311079.24910403 -6060081.37923322 -3226564.18300095] m
Encke Execution Time: 0.337163 seconds

4. Comparison Summary

Let’s compare the magnitude difference between both methods and speed ratio.

pos_diff = np.linalg.norm(r_f_cowell - r_f_encke)
speedup = t_cowell / t_encke if t_encke > 0 else 1.0

print(f"Final Position Difference: {pos_diff:.6f} m")
print(f"Speed Advantage (Cowell/Encke ratio): {speedup:.2f}x")
print("\nNote: Encke is often faster/more accurate over long continuous orbits because it integrates smaller variables, reducing numerical noise limits on large loops.")

Final Position Difference: 273.951264 m
Speed Advantage (Cowell/Encke ratio): 0.29x

Note: Encke is often faster/more accurate over long continuous orbits because it integrates smaller variables, reducing numerical noise limits on large loops.