[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