[SciPy-User] Possible to integrate an ODE just until the solution reaches a certain value?

[Jonathan Stickel]

On 3/10/11 11:00 , scipy-user-request at scipy.org wrote:
> From: Rob Clewley<rob.clewley at gmail.com>
> Subject: Re: [SciPy-User] SciPy-User] Possible to integrate an ODE
> 	just until the solution reaches a certain value?
> To: SciPy Users List<scipy-user at scipy.org>
>
> Hi,
>
>> >  I'm trying to model the dynamics of a catapult-like mechanism used to launch a
>> >  projectile, and have a system of ODEs which I need to numerically integrate over
>> >  time.
> Does your project require numerical precision or just "looks right" accuracy?
>
>>> >>
>>> >>  I would like the ODE solver to stop integrating once the the solution reaches
>>> >>  this certain value, and I will use the states at that point to compute the
>>> >>  initial conditions to another ODE describing the motion from that time onward.
> This is often referred to as a hybrid system.
>
>>> >>  Is there an ODE solver in Python/SciPy which will integrate from the initial t
>>> >>  until the solution reaches a certain value, or until a specific condition is
>>> >>  met? ?The ODE solvers in Matlab have "events" which will do this, but I'm trying
>>> >>
>>> >>  my best to stick with Python.
> PyDSTool is a pure python implementation of event-based hybrid systems
> of ODEs (or discrete mappings), but in your case it may only be
> worthwhile to set up if you need accurate calculations and/or possibly
> more complex hybrid models. (There's some syntax overhead in setting
> up hybrid models.)
>
>> >
>> >  If I understand what you are asking, you can do it with the ode class
>> >  integrator (scipy.integrate.ode). ?Below is a short toy example. ?The
>> >  key is how you setup your loop (while loop with solution criteria vs.
>> >  for loop over time).
>> >
> Just FYI, the example given using the scipy solver is only fine if you
> just want a "quick and dirty" demonstration. If you care about
> accuracy then this will not work: the "result<  4.0" condition does
> not guarantee that you will stop*at*  the point, typically you will
> stop somewhere close but before the point you wish to switch ODEs. You
> would have to (inefficiently) set dt to be very small to resolve the
> event accurately.
>
> An efficient and accurate way to do this is in the PyDSTool
> integrators or in Sundials, but the latter is not pure python.

You could add some more to my scipy example for a rudimentary dynamic 
time step to get a somewhat more precise answer:

...
dt0 = 0.1
dtlow = 1e-5
target = 4.0
while solver.successful() and result < target and t[-1]<100.0:
     t.append(t[-1]+dt)
     solver.integrate(t[-1])
     f.append(solver.y)
     prevres = result
     result = f[-1][0]

     r = (result - prevres)/dt
     dt1 = (target - result)/r
     if dt1 < dt0 and dt1 > 0:
         if dt1 > dtlow:
             dt = dt1
         else:
             dt = dtlow
     else:
         dt = dt0
...

But I am sure PyDSTool and Sundials are much better tools for this 
(haven't tried them yet).

Jonathan