[SciPy-User] Possible to integrate an ODE just until the solution reaches a certain value?
Jonathan Stickel
jjstickel@vcn....
Thu Mar 10 12:28:48 CST 2011
On 3/10/11 11:00 , scipy-user-request@scipy.org wrote:
> From: Rob Clewley<rob.clewley@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@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
More information about the SciPy-User
mailing list