Numerical Methods

Overview

High-altitude balloon trajectory simulation presents unique numerical challenges. The system evolves across multiple time scales: rapid thermodynamic responses occur in seconds, while full flights span hours. Altitude ranges from sea level to 40+ km introduce dramatic changes in atmospheric properties. This section details the advanced numerical methods that ensure accurate, stable, and efficient simulation.

Dormand-Prince 8(7) Method

The simulator employs the Dormand-Prince 8(7) method, a member of the Runge-Kutta family specifically designed for solving non-stiff ordinary differential equations with high accuracy. This 8th-order method with 7th-order error estimation provides exceptional accuracy for smooth balloon dynamics.

Mathematical Foundation

For the system of ODEs:

The Dormand-Prince method computes 13 stages:

The 8th-order solution and 7th-order error estimate are:

Butcher Tableau

The method is defined by its Butcher tableau. Key coefficients include:

Adaptive Timestep Control

Efficient simulation requires dynamic timestep adjustment based on solution behavior. The simulator implements sophisticated PI (Proportional-Integral) control for smooth timestep adaptation.

Error-Based Step Size Control

Given error estimate $\text{err}$ and tolerance $\text{tol}$, the optimal timestep is:

The PI controller smooths timestep changes:

Multi-Component Error Scaling

Different state variables have different scales and tolerances:

def compute_error_scale(state, atol, rtol):
    """Compute error scaling for each component"""
    scales = {}

    # Position: meters
    scales['position'] = atol['position'] + rtol * np.abs(state.position)

    # Velocity: m/s
    scales['velocity'] = atol['velocity'] + rtol * np.abs(state.velocity)

    # Temperature: Kelvin
    scales['temperature'] = atol['temperature'] + rtol * np.abs(state.temperature)

    # Mass: kg (use relative tolerance only for small masses)
    scales['mass'] = atol['mass'] + rtol * np.abs(state.mass)

    return scales

def compute_weighted_error(error_vec, scales):
    """Compute RMS norm of scaled error"""
    scaled_errors = []

    for component, error in error_vec.items():
        scaled_error = error / scales[component]
        scaled_errors.extend(scaled_error.flatten())

    return np.sqrt(np.mean(np.square(scaled_errors)))

Timestep Constraints

Physical constraints impose additional timestep limits:

Multi-Scale Integration

Balloon dynamics span multiple time scales, requiring specialized handling:

Time Scale Separation

Stiffness Detection and Handling

The simulator monitors for stiffness using the eigenvalue ratio:

def detect_stiffness(jacobian):
    """Detect system stiffness via eigenvalue analysis"""
    eigenvalues = np.linalg.eigvals(jacobian)
    lambda_max = np.max(np.abs(eigenvalues))
    lambda_min = np.min(np.abs(eigenvalues[eigenvalues != 0]))

    stiffness_ratio = lambda_max / lambda_min

    if stiffness_ratio > 1000:
        return True, stiffness_ratio
    return False, stiffness_ratio

Implementation Architecture

Core Integration Loop

class DormandPrince87Integrator:
    """High-order adaptive integrator for balloon dynamics"""

    def __init__(self, atol=1e-10, rtol=1e-10):
        self.atol = atol
        self.rtol = rtol
        self.setup_coefficients()

    def integrate(self, dynamics, state, t_start, t_end):
        """Integrate from t_start to t_end"""
        t = t_start
        h = self.estimate_initial_timestep(dynamics, state, t)

        trajectory = [state.copy()]
        times = [t]

        while t < t_end:
            # Adaptive step
            state_new, error, h_new = self.step(dynamics, state, t, h)

            if error <= 1.0:  # Accept step
                state = state_new
                t += h
                trajectory.append(state.copy())
                times.append(t)

            # Update timestep
            h = min(h_new, t_end - t)
            h = self.apply_timestep_limits(h, state)

        return times, trajectory

    def step(self, dynamics, state, t, h):
        """Single Dormand-Prince step with error estimation"""
        # Compute stages
        k = np.zeros((13, len(state)))

        k[0] = dynamics.compute_derivatives(t, state)

        for i in range(1, 13):
            t_i = t + self.c[i] * h
            state_i = state + h * np.sum(self.a[i,:i] * k[:i].T, axis=1)
            k[i] = dynamics.compute_derivatives(t_i, state_i)

        # 8th order solution
        state_8 = state + h * np.sum(self.b * k.T, axis=1)

        # 7th order solution for error
        state_7 = state + h * np.sum(self.b_hat * k.T, axis=1)

        # Error estimation
        error_vec = state_8 - state_7
        error = self.compute_error_norm(error_vec, state)

        # PI controller for timestep
        h_new = self.pi_timestep_control(h, error)

        return state_8, error, h_new

State Vector Management

class SimulationState:
    """Complete state vector for balloon simulation"""

    def __init__(self):
        # Primary state variables
        self.position = np.zeros(3)      # [x, y, z] in ECEF
        self.velocity = np.zeros(3)      # [vx, vy, vz] in ECEF
        self.temperature_film = 273.15   # Kelvin
        self.temperature_gas = 273.15    # Kelvin
        self.mass_gas = 0.0             # kg

        # Derived quantities (not integrated)
        self.pressure_gas = 101325       # Pa
        self.volume = 0.0               # m³
        self.radius = 0.0               # m
        self.altitude = 0.0             # m MSL

    def to_array(self):
        """Convert to array for integration"""
        return np.concatenate([
            self.position,
            self.velocity,
            [self.temperature_film],
            [self.temperature_gas],
            [self.mass_gas]
        ])

    def from_array(self, array):
        """Update from integration array"""
        self.position = array[0:3]
        self.velocity = array[3:6]
        self.temperature_film = array[6]
        self.temperature_gas = array[7]
        self.mass_gas = array[8]

Dense Output Interpolation

For smooth trajectory output and event detection, the integrator provides dense output:

def dense_output(self, t_eval, h, y_n, k_stages):
    """Evaluate solution at arbitrary points within step"""
    theta = (t_eval - self.t_n) / h

    # 7th order interpolation polynomials
    b_theta = np.zeros(13)
    for i in range(13):
        b_theta[i] = sum(self.d[i,j] * theta**j for j in range(8))

    return y_n + h * np.sum(b_theta * k_stages.T, axis=1)

Error Analysis & Stability

Local Truncation Error

The local truncation error for DP8(7) is:

For typical balloon dynamics with $|y^{(9)}| \sim 10^{-3}$:

Stability Region

The stability region for explicit Runge-Kutta methods is bounded. For DP8(7):

Global Error Propagation

Over a complete trajectory, errors accumulate as:

def estimate_global_error(self, trajectory, dynamics):
    """Estimate accumulated global error"""

    # Estimate Lipschitz constant
    L = self.estimate_lipschitz_constant(dynamics, trajectory)

    # Time span
    t_span = trajectory[-1].time - trajectory[0].time

    # Global error bound
    global_error = (np.exp(L * t_span) - 1) / L * self.atol

    return {
        'position_error': global_error * 1.0,  # m
        'velocity_error': global_error * 0.1,  # m/s
        'temperature_error': global_error * 0.01,  # K
        'relative_error': global_error / np.mean([s.altitude for s in trajectory])
    }

Performance Optimization

Computational Efficiency

The DP8(7) method requires 13 function evaluations per step. Performance optimizations include:

# Vectorized stage computation
k_stages = np.zeros((13, state_dim))
for i in range(13):
    # Compute all components simultaneously
    k_stages[i] = dynamics.derivatives_vectorized(t + c[i]*h,
                                                  state + h*a[i]@k_stages)

# Reuse Jacobian for multiple steps
if steps_since_jacobian > 5 or error > 0.1:
    jacobian = dynamics.compute_jacobian(state)
    steps_since_jacobian = 0

Memory Management

class MemoryEfficientIntegrator:
    """Minimize memory allocation during integration"""

    def __init__(self, state_dim):
        # Pre-allocate work arrays
        self.k_stages = np.zeros((13, state_dim))
        self.state_temp = np.zeros(state_dim)
        self.error_weights = np.zeros(state_dim)

    def step(self, dynamics, state, t, h):
        """Step with minimal allocation"""
        # Reuse pre-allocated arrays
        np.copyto(self.state_temp, state)

        # In-place operations
        for i in range(13):
            self.compute_stage_inplace(i, dynamics, t, h)

        # Update state in-place
        self.apply_step_inplace(state, h)

        return self.compute_error()

Parallel Integration

For ensemble simulations or sensitivity analysis:

def parallel_integrate(dynamics_list, states, t_span, n_cores=4):
    """Integrate multiple trajectories in parallel"""

    with multiprocessing.Pool(n_cores) as pool:
        # Prepare work packages
        work = [(dynamics, state, t_span) for dynamics, state
                in zip(dynamics_list, states)]

        # Parallel execution
        trajectories = pool.starmap(integrate_single, work)

    return trajectories

Special Case Handling

Burst Detection

Accurate burst detection requires special handling:

class BurstDetector:
    """Detect balloon burst using root finding"""

    def __init__(self, burst_diameter):
        self.burst_diameter = burst_diameter
        self.interpolator = None

    def check_burst(self, state_prev, state_curr, h):
        """Check for burst between states"""
        if state_curr.diameter >= self.burst_diameter:
            # Set up interpolation
            self.setup_interpolation(state_prev, state_curr, h)

            # Find exact burst time via root finding
            t_burst = self.find_burst_time()
            state_burst = self.interpolate_state(t_burst)

            return True, t_burst, state_burst
        return False, None, None

    def find_burst_time(self):
        """Binary search for burst time"""
        def burst_condition(theta):
            state = self.interpolate_state(theta)
            return state.diameter - self.burst_diameter

        # Bracket and solve
        theta_burst = scipy.optimize.brentq(burst_condition, 0, 1)
        return self.t_prev + theta_burst * self.h

Ground Impact

def detect_ground_impact(state, terrain_model):
    """Detect and handle ground contact"""
    ground_altitude = terrain_model.get_elevation(state.lat, state.lon)

    if state.altitude <= ground_altitude:
        # Compute impact conditions
        impact_velocity = np.linalg.norm(state.velocity)
        impact_angle = np.arccos(-state.velocity[2] / impact_velocity)

        return {
            'impact': True,
            'location': (state.lat, state.lon),
            'velocity': impact_velocity,
            'angle': np.degrees(impact_angle)
        }
    return {'impact': False}

Numerical Instability Recovery

def recover_from_instability(self, state, h):
    """Recover from numerical difficulties"""

    # Reduce timestep aggressively
    h_recovery = h * 0.1

    # Check for NaN or infinity
    if not np.all(np.isfinite(state.to_array())):
        raise SimulationError("Non-finite values detected")

    # Check physical constraints
    if state.temperature_gas < 0 or state.mass_gas < 0:
        raise SimulationError("Unphysical state detected")

    # Limit extreme velocities
    v_mag = np.linalg.norm(state.velocity)
    if v_mag > 1000:  # m/s - unrealistic for balloon
        state.velocity *= 100 / v_mag
        warnings.warn("Velocity limited due to instability")

    return h_recovery