[SciPy-User] odeint problems with a simple ode

[Laura Matrajt]

Hi there,
 I have a simple ODE that I was trying to integrate with odeint and ran
into weird problems:
the ode is basically this:

ds/dt = 0 if t <b
            constant if b<= t <= b+a
            0 t>b+a

this integrates to a stair-type function, with the "stair" part being a
line.
I know this can be integrated analytically but this is just part of a more
complicated ode.

If one tries to integrate this using odeint, the solver will work fine as
long as b<a but will fail for b>a.
In fact, I checked the values that t was taking, and it is going all over
the place.
I know that this function is non differentiable in two points, but it is so
easy that I was hoping the integrator would be able to handle this. In
fact, other solvers (my own RK4, matlab, etc) can deal with this without
any problem.

Any ideas of how/why this is happening would be very much appreciated.
thanks,


Below is a minimal (not)working example:

import matplotlib

import numpy as np

import pylab as plt

from scipy.integrate import odeint


def frac(t, a, fraction, initS0, b):

    if (t>0) & (t<= b):

        frac = 0

    elif (t> b) & (t<= b + a):

        frac = (1.0*initS0*fraction/a)

    else:

        frac = 0

    return frac


def linearModel(y,t,params):

    [a, fraction, initS0,b] = params

    out = [-frac(t, a, fraction, initS0, b),

           frac(t,a, fraction, initS0, b)

           ]


    return out


if __name__ == '__main__':

    a = 15.0

    b = 10.0

    time = np.linspace(0,100, 101)

    initLM = 1000.0

    initCondLM = [initLM, 0]

    paramsLM = [a,0.1,initLM, b]

    LM = odeint(linearModel, initCondLM, time, args=(paramsLM,))

    plt.plot(time,LM[:,0],'b')

    plt.plot(time,LM[:,1],'r')
    plt.show()

--
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20140916/39d31a93/attachment.html>