python - algebraic constraint to terminate ODE integration with scipy -

i'm using scipy 14.0 solve system of ordinary differential equations describing dynamics of gas bubble rising vertically (in z direction) in standing still fluid because of buoyancy forces. in particular, have equation expressing rising velocity u function of bubble radius r, i.e. u=dz/dt=f(r), , 1 expressing radius variation function of r , u, i.e. dr/dt=f(r,u). rest appearing in code below material properties. i'd implement account physical constraint on z which, obviously, limited liquid height h. consequently implemented sort of z<=h constraint in order stop integration in advance if needed: used set_solout in order so. situation code runs , gives results, set_solout not working @ (it seems z_constraint never called actually...). know why? there more clever idea, may in order interrupt when z=h (i.e. final value problem) ? right way/tool or should reformulate problem?

thanks in advance

emi

from scipy.integrate import ode   db0 = 0.001 # init bubble radius y0, t0 = [ db0/2 , 0. ], 0. #init conditions h = 1  def y_(t,y,g,p0,rho_g,mi_g,sig_g,h):     r = y[0]      z = y[1]     z_ = ( r**2 * g * rho_g ) / ( 3*mi_g )  #velocity     r_ = ( r/3 * g * rho_g * z_ ) / ( p0 + rho_g*g*(h-z) + 4/3*sig_g/r ) #r dynamics         return [r_, z_]  def z_constraint(t,y):     h = 1  #should rather variable..     z = y[1]      if z >= h:         flag = -1      else:         flag = 0     return flag  r = ode( y_ ) r.set_integrator('dopri5') r.set_initial_value(y0, t0) r.set_f_params(g, 5*1e5, 2000, 40, 0.31, h) r.set_solout(z_constraint)  t1 = 6 dt = 0.1  while r.successful() , r.t < t1:     r.integrate(r.t+dt) 

you're running this issue. set_solout work correctly, must called right after set_integrator, before set_initial_value. if introduce modification code (and set value g), integration terminate when z >= h, want.

to find exact time when bubble reached surface, can make change of variables after integration terminated solout , integrate respect z (rather t) z = h. paper describes technique m. henon, physica 5d, 412 (1982); may find this discussion helpful. here's simple example in time t such y(t) = 0.5 found, given dy/dt = -y:

import numpy np scipy.integrate import ode  def f(t, y):     """exponential decay: dy/dt = -y."""     return -y  def solout(t, y):     if y[0] < 0.5:         return -1     else:         return 0  y_initial = 1 t_initial = 0  r = ode(f).set_integrator('dopri5') r.set_solout(solout) r.set_initial_value(y_initial, t_initial)  # integrate until solout constraint violated r.integrate(2)  # new system t independent variable: see henon's paper details. def g(y, t):     return -1.0/y  r2 = ode(g).set_integrator('dopri5') r2.set_initial_value(r.t, r.y)  r2.integrate(0.5)  y_final = r2.t t_final = r2.y  # error: difference between found , analytical solution print t_final - np.log(2)