[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