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

[Jonathan Stickel]

On 3/9/11 15:28 , scipy-user-request at scipy.org wrote:
> Date: Wed, 9 Mar 2011 14:20:12 -0800 (PST)
> From: Randy Williams<skapunxter at yahoo.com>
> Subject: [SciPy-User] Possible to integrate an ODE just until the
> 	solution	reaches a certain value?
> To:scipy-user at scipy.org
> Message-ID:<539506.47663.qm at web39409.mail.mud.yahoo.com>
> Content-Type: text/plain; charset=us-ascii
>
> Greetings,
>
>
> 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.  I am trying to solve for the position of the projectile as well as the
> other components in my mechanism.  At some point in time, the projectile
> separates from the mechanism, and becomes airborne.  The equations governing the
>
> system change at that point in time, but because it's a function of position
> (which i'm solving for), I don't know up front what timespan to integrate over.
>
> 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.
> 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.
>
> Thanks,
> Randy

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).

HTH,
Jonathan


import numpy as np
from scipy.integrate import ode
from matplotlib.pyplot import *

def dfdt(t, f, a,b,c,d):
     x = f[0]
     y = f[1]
     dxdt = np.sin(a*x + b*y)
     dydt = np.cos(c*x + d*y)
     return [dxdt, dydt]

f0 = [1., 0.]
a = 1.0
b = -2.0
c = -1.0
d = 1.0

t = [0.0]
dt = 0.1
f = [f0]
result = f[0][0]
solver = ode(dfdt)
solver.set_initial_value(f0,t[0])
solver.set_f_params(a,b,c,d)
while solver.successful() and result < 4.0 and t[-1]<100.0:
     t.append(t[-1]+dt)
     solver.integrate(t[-1])
     f.append(solver.y)
     result = f[-1][0]