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

Jose Guzman sjm.guzman@googlemail....
Thu Jun 21 10:50:05 CDT 2012


On 21/06/12 17:35, Warren Weckesser wrote:
> 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.
Wov, this was absolutely amazing! I was not aware that odeint uses an 
adaptive ODE solver.

Just one more question. Now I am applying a complex pulse:

def pulse(t):
     """ t in ms """
     if t<=1 or t>=2:
         return 0.0
     else:
         return 1.0 # in mM

# construct a 10 times pulse
step = np.linspace(0,10, 1000)
p = list()
# 10 pulses
for i in range(10):
     p += [pulse(t) for t in step]
time = np.linspace(0, 10*10, 1000*10)

However, the solver only shows the first peak:
y = odeint(func = diff_eq, y0=0, t= time, args=(4,1), hmax=0.01)

you can have a look here: 
http://aleph.sagemath.org/?q=3e9bc8e9-dbb0-4a1d-b495-dbd93dd6e594

Is there any other option in odeint that I should be awared of? Should I 
use another ODE solver?

Thanks a lot Warren! I really appreciate your help!

Jose




More information about the SciPy-User mailing list