Solve_ivp Integration Gets Stuck If Initial Condition Triggers Event With Terminal = True

I work in spacecraft trajectory design and I am implementing an algorithm for finding periodic orbits that needs to detect whenever a trajectory crosses a specific plane. I am usin

Solution 1:

First check the system of differential equations. How do the linear terms make sense? The third component of the acceleration has its central force part proportional to the third coordinate Y[2]. It is quite possible that with this substantial correction your original approach using a non-terminal event will work.

If you decide to implement the ODE function in the vectorized version, then be fully compliant with the requirement. You can not suppose that the input vector contains only one state, make the output Ydot have the same format as the input Y.

defeq_motion(t,Y):

    Ydot = np.zeros(Y.shape)
    r = (Y[0]**2 + Y[1]**2 + Y[2]**2)**0.5

    Ydot[0] = Y[3]
    Ydot[1] = Y[4] 
    Ydot[2] = Y[5]
    Ydot[3] = - Y[0]/r**3 + 2*Y[4] + 3*Y[0]
    Ydot[4] = - Y[1]/r**3 - 2*Y[3]         
    Ydot[5] = - Y[2]/r**3 - Y[2]          

    return Ydot

It is better to also update t0 along with Y0, so set

t0, tfinal = 0, 10

Use event direction to prohibit the solver to detect the last event again, ideally you would from Y0 and eq_motion(t0,Y0) detect the current sign of the event in positive time direction, here positive, and set the direction opposite to it

event_xz_plane.direction = -1

n_plane_intersections = 4for i in np.arange(n_plane_intersections):

    integrate = solve_ivp(fun=eq_motion, t_span=[t0,tfinal], y0=Y0, method = 'LSODA', 
                            vectorized = True, rtol = 1e-13, atol = 1e-14,
                            events = (event_xz_plane),dense_output=True)

    print(integrate.message)
    print('Number of integrations before event: {}'.format(integrate.y.shape[1]-1))
    print(integrate.t_events[0],integrate.y_events[0])
    print('\n')

    # attempting to chain the different integration arcs
    Y0 = integrate.y[:,-1]
    t0 = integrate.t[-1]
    event_xz_plane.direction *= -1

This gives the output

A termination event occurred.
Number of integrations before event: 288
[0.96633212] [[ 3.54891748e-01 -6.67868538e-17 -8.77207053e-01  5.03973235e-02
  -7.78052855e-01 -1.76219528e-02]]


A termination event occurred.
Number of integrations before event: 263
[2.03187275] [[ 7.62458860e-03  1.66901228e-15  1.66823579e-01 -2.71527186e-01
   3.27865867e+00  1.07443749e-01]]


A termination event occurred.
Number of integrations before event: 300
[3.22281063] [[ 3.86773384e-01 -1.92554306e-16 -8.24376230e-01 -6.86117012e-02
  -8.96669320e-01  2.07695448e-01]]


A termination event occurred.
Number of integrations before event: 258
[4.10715034] [[-2.07546473e-02 -9.56266316e-16  1.42957911e-01  9.37459643e-02
   3.56344204e+00  7.24687825e-02]]

To compare with, the Dormand-Prince 853 method is much faster (and Radau much slower) and gives the same results

A termination event occurred.
Number of integrations before event: 67
[0.96633212] [[ 3.54891748e-01 -2.46764414e-16 -8.77207053e-01  5.03973235e-02
  -7.78052855e-01 -1.76219528e-02]]


A termination event occurred.
Number of integrations before event: 52
[2.03187275] [[ 7.62458861e-03 -5.10008702e-16  1.66823579e-01 -2.71527186e-01
   3.27865867e+00  1.07443749e-01]]


A termination event occurred.
Number of integrations before event: 58
[3.22281063] [[ 3.86773384e-01  1.21430643e-17 -8.24376230e-01 -6.86117012e-02
  -8.96669320e-01  2.07695448e-01]]


A termination event occurred.
Number of integrations before event: 52
[4.10715034] [[-2.07546473e-02  9.19403442e-17  1.42957911e-01  9.37459642e-02
   3.56344204e+00  7.24687825e-02]]