[SciPy-User] scipy.integrate.odeint returns zeros when changing function

Warren Weckesser warren.weckesser@enthought....
Thu Jun 21 10:35:23 CDT 2012


On Thu, Jun 21, 2012 at 10:20 AM, Jose Guzman <sjm.guzman@googlemail.com> wrote:
> Hello everybody!
>
> I am trying to solve a basic 2-state kinetic scheme of the form:
>
> # Two-state kinetic model
> #
> #        alpha*[cc]
> # (1-n) ----------> n
> #       <----------
> #          beta
>
>
> where [cc] is a function step function. To solve it numerically, I used
> scipy.integrate.odeint. Everything works fine if I define a step function
> that takes the value 1 between 2 and 3 .
>
>
> from numpy import linspace
> from scipy.integrate import odeint
>
> # define a square pulse with the concentration of Glu
>
> def pulse(t):
>     """ t in ms """
>     if t<=2 or t>=3:
>         return 0.0
>     else:
>         return 1.0 # in mM
>
> # independent variable
> time = linspace(0,10, 1000) # 10 ms
>
>
> # define differential equation
> def diff(n, t, alpha, beta):
>     """
>     alpha in 1/(mM*ms)
>     beta in 1/ms
>     dn/dt = alpha*cc*(1-n) - beta*n
>     """
>     cc = pulse(t) # square pulse
>     return alpha*cc*(1-n)-beta*n
>
> # and solve it for alpha = 1 and beta = 1
> y = odeint(func = diff, y0=0, t= time, args=(1,1))
>
>
> However, if i change the pulse function to return 1 between, say 5 and 6 the
> solution to the differential equation returns 0 for all the interval. I
> cannot figure out why.
>
> # new function
> def pulse(t):
>     """ t in ms """
>     if t<=5 or t>=6:
>         return 0.0
>     else:
>         return 1.0 # in mM
> y = odeint(func = diff, y0=0, t= time, args=(1,1)) # returns all zeros
>
> I would greatly appreciate some help here
>


Short answer: try adding the keyword argument hmax=0.1 in your call to odeint.

Longer answer:  If the pulse were not there, the exact solution to
your problem (with initial condition 0) is y(t) = 0.  The solver is
adaptive, and with a solution this nice, it will quickly adapt to very
large step sizes.  It is likely that the solver has simply skipped
right past your pulse by taking a step greater than 1 time unit.
Setting hmax=0.1 tells the solver to never take a step larger than
0.1.

(I'll refrain from my usual advice against using discontinuous
functions in a solver that is based on certain smoothness
assumptions...)

Warren


> Best
>
> Jose
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>


More information about the SciPy-User mailing list