Simulating engine burn with scipy ode solve

Question

I'm using ODE solver to calculate complex trajectories.

A simplified setup is like this:

def calc_dydt(t, y):
    julian_datetime = get_julian_datetime(t)

    bodies_r = {
        SUN: np.array([0.0, 0.0, 0.0]),
        EARTH: np.array(get_planet_xyz(EARTH, julian_datetime)),
        MARS: np.array(get_planet_xyz(MARS, julian_datetime)),
        JUPITER: np.array(get_planet_xyz(JUPITER, julian_datetime))
    }

    r = np.array(y[0:3])
    v = y[3:6]

    a = np.array([0.0, 0.0, 0.0])
    for body in [SUN, EARTH, MARS, JUPITER]:
        relative_r = r - bodies_r[body]
        relative_norm_r = np.linalg.norm(relative_r)
        a += -relative_r * MU[body] / relative_norm_r**3

    return np.concatenate([v, a])

And then the simulation loop

while solver.successful() and i < rounds:
    solver.integrate(solver.t + dt)
    results[i] = solver.y

Now I would like to simulate a burn at time $t_{burn}$ with some $\delta V$. I'm interested in both:

1. Instant increase of the velocity (simulating an engine with high specific impulse)

2. Gradual increase of the velocity over time (simulating low specific impulse)

How can I adjust my calc_dydt() method to add such planned burns?

For (1): Can I just add needed amount of $\delta V$ to the variable v without reflecting the engine thrust in the variable a in calc_dydt()?

Answer 1

For an instantaneous delta V, you definitely want to have the integrator stop exactly at the point in time where the change in velocity is to be applied. Dynamic step sized integrators stop where they want to stop. You'll need to force the issue and make the integrator stop at the desired point in time. You can specify a step size that makes a multistep integrator such as Adams-based techniques have a step end at the point in time where the instantaneous delta V is to be applied, but that too is undesirable. The instantaneous change in velocity invalidates the mathematics that underlie such techniques. Integrate up to the point where the delta V is to be applied, change the velocity, and integrate to the end (or to the next instantaneous delta V, if there is one).

What about finite burns? There's still a discontinuity here, but in the second derivative rather than first. Numerical ODE solvers can yield notoriously lousy results when the derivatives undergo a discontinuity. Discontinuities are particularly problematic with multistep integrators and with techniques that have a dynamic step size (e.g. Dormand-Prince). All of the scipy solvers are of one form or the other. It is best once again to force the integrator to stop and restart at such thruster boundaries.

One last comment:

def calc_dydt(t, y):
     julian_datetime = get_julian_datetime(t)
     bodies_r = {
        SUN: np.array([0.0, 0.0, 0.0]),
        EARTH: np.array(get_planet_xyz(EARTH, julian_datetime)),
        MARS: np.array(get_planet_xyz(MARS, julian_datetime)),
        JUPITER: np.array(get_planet_xyz(JUPITER, julian_datetime))
    }

   r = np.array(y[0:3])
    v = y[3:6]

   a = np.array([0.0, 0.0, 0.0])
    for body in [SUN, EARTH, MARS, JUPITER]:
        relative_r = r - bodies_r[body]
        relative_norm_r = np.linalg.norm(relative_r)
        a += -relative_r * MU[body] / relative_norm_r**3

   return np.concatenate([v, a])

Given that you are fixing the Sun at the origin means your calculation of acceleration is not quite valid. A body-centered frame such as this is not inertial. To be correct you'll either need to integrate in a solar system barycenter frame or use a heliocentric frame and account for the fact that the Sun is accelerating toward the planets. There are several questions and answers on this site that address such third body effects.

Answer 2

Here's an example using a "soft" normalized Gaussian bump for the impulse.

$$ \frac{1}{\sigma_1 \sqrt{2 \pi}} \exp\left(-\frac{1}{2}\left(\frac{t-t_0}{\sigma_1} \right)^2 \right) \mathbf{a_{bump}} $$

You can make it quite short, but even a short ramp up and down gives the integrator a chance to notice that things are changing and to reduce its internal step size accordingly. Remember that the time points you give it as input are usually interpolation points and the solution inside is on its own variable step-size grid.

When you turn on your ion engine you can also ramp up your thrust with

$$\frac{1}{2}\left(1 + \text{erf}\left( \frac{t-t_0} {\sigma_2} \right) \right) \hat{\mathbf{v}} $$

where in this case it's directed along the direction of motion.

One way to get a rough idea how well these behave is to run the same simulation with a wide range of abruptnesses sigma_one and sigma_two.

SciPy's odeint does a pretty good job, it switched between a non-stiff and stiff integrator internally. You can see some details of what's going on inside by examining the info dictionary it returns. However it failed for me on a very stiff problem as described in "Pythagorean Three Body Problem" - need some points from an accurate solution for comparison

You can see more about the accumulated error in need to understand better how rtol, atol work in scipy.integrate.odeint.

Solutions to next step in ODE solver testing for the “Pythagorean 3-Body Proxblem” are currently inconclusive; if you get it t work and can post a better answer I'll probably accept it!

If you want to learn how to get really accurate you can also think about reading answers to What does “symplectic” mean in reference to numerical integrators, and does SciPy's odeint use them?

Here is a Gaussian bump plane change after two periods and a retro-propulsive burn starting at four periods:

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf
from scipy.integrate import odeint as ODEint

def soft_impulse(t, t_zero, sigma):
    norm = 1. / (sigma * np.sqrt(2. * np.pi))
    return norm * np.exp(-0.5 * ((t - t_zero) / sigma)**2)

def deriv(X, t, t_zero, t_one, sigma_one, sigma_two, bump, retro):
    x, v = X.reshape(2, -1)
    vnorm = v / np.sqrt((v**2).sum())
    acc = -x * ((x**2).sum())**-1.5   # gravity
    acc += bump * soft_impulse(t, t_zero, sigma_one) # impulse
    acc += retro * vnorm * 0.5 * (1. + erf((t-t_one)/sigma_two))  # propulsion
    return np.hstack((v, acc))

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]

X0 = np.array([1, 0, 0] + [0, 1, 0], dtype=float)
times = np.linspace(0, 6*twopi, 200)

t_zero, t_one, sigma_one, sigma_two, retro = 2*twopi, 4*twopi, 0.2, 0.1, -0.05
bump = np.array([0, 0, 0.1])

answer, info = ODEint(deriv, X0, times, full_output=True, atol=1E-10, 
                      args=(t_zero, t_one, sigma_one, sigma_two, bump, retro))
x, v = answer.T.reshape(2, 3, -1)

plt.figure()
plt.subplot(2, 1, 1)
for thing in x:
    plt.plot(times/twopi, thing)
plt.title('x', fontsize=14)
plt.subplot(2, 1, 2)
for thing in v:
    plt.plot(times/twopi, thing)
plt.title('v', fontsize=14)
plt.xlabel('t / twopi', fontsize=14)
plt.show()