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

[Rob Clewley]

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. An
example of using PyDSTool events to switch between sub-systems is
given in IF_squarespike_model.py at

http://pydstool.bzr.sourceforge.net/bzr/pydstool/revision/1#tests/IF_squarespike_model.py

which demonstrates an "integrate and fire" neuron model with a fixed
rectangular pulse for a spike. There are several other demos of hybrid
models provided in the package, or you can ask me.

-Rob