How to stop ODE simulations after impact using solve_ivp event detection

by

Summary

The simulation fails to terminate upon impact because it uses a fixed time-stepping approach via scipy.integrate.odeint. Instead of detecting the collision event, the integrator blindly follows the mathematical differential equations, pushing the object into the negative altitude space (inside the planet) once the orbital path crosses the surface radius.

Root Cause

The issue stems from two fundamental architectural flaws in the simulation logic:

  • Lack of Event Detection: The integrator is told to solve the equations for a fixed time array (tout). It has no instruction to monitor the state variables (specifically the radius $r$) to stop when a boundary condition is met.
  • Mathematical Continuity vs. Physical Reality: The physics model defines gravity and drag as continuous functions. Without a “stop” mechanism, the differential equations remain valid mathematically even when $r < R_{planet}$, causing the object to “tunnel” through the surface and continue simulating as if it were traveling through the planet’s core.

Why This Happens in Real Systems

In high-fidelity engineering simulations, this is known as a Boundary Violation. It happens because:

  • Fixed Step Integration: If the time step is too large, the object can “leap” from a positive altitude to a negative altitude in a single step, completely missing the exact moment of impact.
  • Discontinuous Physics: Real-world collisions involve instantaneous changes in velocity (impulse) or a total change in state (destruction). Differential equations (ODEs) assume smoothness; they cannot naturally handle the discontinuity of a collision without explicit event handling.

Real-World Impact

  • Inaccurate Data: Post-impact trajectories are physically impossible, leading to incorrect calculations of total energy dissipation or terminal velocity.
  • Numerical Instability: As the object enters the planet, the density and gravity models may produce extreme values, leading to floating-point overflows or “NaN” errors that crash the entire simulation pipeline.
  • Resource Waste: In large-scale Monte Carlo simulations, simulating “ghost” trajectories after a collision wastes significant computational budget.

Example or Code (if necessary and relevant)

To fix this, we must transition from scipy.integrate.odeint to scipy.integrate.solve_ivp, which supports Event Detection.

import numpy as np
import scipy.integrate as sci

def collision_event(t, state, Rplanet):
    x, y, z = state[0], state[1], state[2]
    r = np.sqrt(x**2 + y**2 + z**2)
    return r - Rplanet

collision_event.terminal = True
collision_event.direction = -1

def derivatives(t, state, Rplanet, mplanet, G, mass, aeroModel, S, CD):
    x, y, z, vx, vy, vz = state
    r = np.sqrt(x**2 + y**2 + z**2)

    # Gravity
    accel = (G * mplanet / r**3) * np.array([x, y, z])

    # Aero
    altitude = r - Rplanet
    rho = aeroModel.getDensity(altitude)
    v_vec = np.array([vx, vy, vz])
    v_mag = np.linalg.norm(v_vec)
    qinf = -0.5 * rho * v_mag * S
    aero_f = qinf * CD * v_vec

    total_accel = accel + (aero_f / mass)
    return [vx, vy, vz, total_accel[0], total_accel[1], total_accel[2]]

# Implementation usage:
# sol = sci.solve_ivp(derivatives, [0, period], stateInitial, 
#                     events=collision_event, args=(Rplanet, mplanet, G, mass, aeroModel, S, CD))

How Senior Engineers Fix It

A senior engineer solves this by implementing Root-Finding Event Handlers:

  1. Switch to Event-Aware Solvers: Move away from legacy solvers (odeint) to modern solvers (solve_ivp) that explicitly allow defining events.
  2. Define Termination Criteria: Create a function that returns $0$ when the object hits the surface ($r – R_{planet} = 0$) and set the terminal attribute to True.
  3. Zero-Crossing Detection: Use a solver that employs a root-finding algorithm (like Brent’s method) to find the exact timestamp where the altitude transitions from positive to negative.
  4. State Cleanup: Ensure the simulation outputs only the valid trajectory up to the event, preventing the downstream analysis of invalid “inside-planet” data.

Why Juniors Miss It

  • The “Black Box” Trap: Juniors often treat integrators as magic boxes that “just work,” without realizing that the solver is purely a mathematical tool that doesn’t know what a “planet” is.
  • Focus on Physics, Not Logic: They spend time refining the aerodynamics or gravity formulas but forget the control logic that governs the simulation lifecycle.
  • Missing the Discontinuity: They assume that if the math is correct, the simulation must be correct, overlooking the fact that the physical boundary condition is a logic constraint, not just a mathematical one.

Leave a Comment