2

Scipy is moving away from odeint towards solve_ivp, which no longer supports passing in additional arguments for the dynamics function. Instead, lambdas are recommended. However, when I try the same for events, they do not work properly. Any thoughts?

MWE (modified from the doc page):

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# dynamics of a simple mass with ballistic flight and a bit of drag
def cannon(t, y, p): 
    return [y[2],y[3], -p['friction']*y[2], p['gravity']-p['friction']*y[3]]

# termination condition: cannonball hits the ground
# this event does not require parameters, but more complex events might
def hit_ground1(t, y, p): 
    return y[1]
hit_ground1.terminal = True
hit_ground1.direction = -1

def hit_ground2(t,y):
    return y[1]
hit_ground2.terminal = True
hit_ground2.direction = -1

p = {'gravity':-1,'friction':0} # paramters as a dict
y0 = [0, 0, 0, 10] # initial conditions
t_span = [0, 22] # integration time a bit over what is necessary

# we can handle dynamics with parameters by using lambdas
# but somehow the same doesn't seem to work with events
sol1 = solve_ivp(fun=lambda t,x:cannon(t,x,p), t_span=t_span, 
    y0=y0, events=hit_ground2, max_step=0.01)    
sol2 = solve_ivp(fun=lambda t,x:cannon(t,x,p), t_span=t_span, 
    y0=y0, events=lambda t,x:hit_ground1(t,x,p), max_step=0.01)

print(sol1.t[-1]) # terminates correctly
print(sol2.t[-1]) # continues integrating
plt.plot(sol1.t,sol1.y[1], linewidth=3)
plt.plot(sol2.t,sol2.y[1],'--',linewidth=3)
plt.show()

enter image description here

Steve Heim
  • 799
  • 12
  • 25
  • Shortly after your question, [`args` have been added again](https://github.com/scipy/scipy/commit/d9f2979ac59b6c0b3c23af08faca04cfeb1fecb9), at least for the RHS function. Cannot find anything for the `event` funtion though... Something I would need for my usecase. Maybe I will try the `lambda` approach. – winkmal Oct 05 '21 at 12:15

1 Answers1

2

The properties of the event, terminal and direction, do not get transferred to your lambda expression. You need to save your lambda into a variable and add the properties there instead of on the hit_ground1 function.

def hit_ground1(t, y, p):
    return y[1]

ground_event = lambda t,x:hit_ground1(t,x,p)
ground_event.terminal = True
ground_event.direction = -1

Use this event and it should work as expected.

 sol2 = solve_ivp(fun=lambda t,x:cannon(t,x,p), t_span=t_span,
    y0=y0, events=ground_event, max_step=0.01)
winkmal
  • 622
  • 1
  • 6
  • 16
Alex bGoode
  • 551
  • 4
  • 16
  • Steve, running print(sol1.t_events); print(sol2.t_events) the solver caught the events, t=20. appeared in both methods. Don't know why the sol2, hitground1 also caught t=0.. This doesn't happen in the solution by Alex – famaral42 Oct 06 '19 at 23:36