[SciPy-User] scipy.integrate.odeint returns zeros when changing function
Jose Guzman
sjm.guzman@googlemail....
Fri Jun 22 03:16:19 CDT 2012
On 21/06/12 17:50, Jose Guzman wrote:
> 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
>
>
Obviously the way I defined the multiple pulse was erroneous and had
nothing to do with the integration method. You can find the solution here:
http://aleph.sagemath.org/?q=7d522c13-3579-40c1-a2ff-4f9a3941f2d6
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20120622/3294038a/attachment.html
More information about the SciPy-User
mailing list