The equations of motion contain singular terms that require careful handling.

The \(1/v\) Singularity¶

Problem¶

At liftoff, velocity is zero: \(v(0) = 0\).

The flight-path angle equation contains:

Both terms \(\to \infty\) as \(v \to 0\).

Solution: Regularized Inverse¶

Replace \(1/v\) with a smooth approximation:

Properties¶

Asymptotic behavior:

At zero:

At \(v = v_\epsilon\):

Choice of \(v_\epsilon\)¶

This is small enough to not affect the physics (\(v > 100\) m/s for most of flight) but large enough to provide stability at liftoff.

Implementation:

# Guard the 1/v terms near v -> 0 to avoid numerical singularity.
# For v >> v_eps this is exactly the standard equation.
# [Eq. 12] Regularized inverse velocity
inv_v = v / (v * v + v_eps * v_eps)

The \(1/r\) Singularity¶

Problem¶

If numerical error pushes \(r < R_E\) (rocket underground), terms like \(\mu/r^2\) explode.

Solution: Safe Radius¶

Implementation¶

# Guard r to avoid singularity at r=0 [Eq. 13]
r_safe = jnp.maximum(r, earth.r_e * 0.99)

Event Detection¶

Additionally, ground impact is detected as a termination event:

def _cond_ground(t, y, args, **kwargs):
    """Event: ground impact (failure)."""
    return y[_R] - args["earth"].r_e

Time Monotonicity¶

Problem¶

Adaptive ODE solvers with event detection can produce non-monotonic time samples:

• Step rejected and retried
• Event interpolation
• Floating-point accumulation

Solution: Filtering¶

def _strictly_increasing_mask(ts: Array) -> Array:
    """Create mask for strictly increasing time samples."""
    n = len(ts)
    mask = np.ones(n, dtype=bool)
    last_t = -np.inf

    for i in range(n):
        if ts[i] > last_t:
            last_t = ts[i]
        else:
            mask[i] = False

    return mask

Apply mask to all trajectory arrays.

Summary¶

References¶

1. Hairer, E. et al. (1993). Solving ODEs I. Springer. Section II.8 - "Numerical Difficulties".

2. Ascher, U.M. & Petzold, L.R. (1998). Computer Methods for ODEs and DAEs. SIAM. Chapter 3 - "Basic Methods".