[SciPy-User] Help!!!!!! having problems with ODEINT

Fri Jul 22 01:39:39 EDT 2011

Hi Laura,

odeint uses an adaptive method (I believe an "embedded 4-5
Runge-Kutta") which evaluates a small number of points and estimates
the error produced by using those points. If the error
is small - and odeint has parameters that let you set how big "small"
is - it'll return this estimate, otherwise it'll subdivide the region
and repeat the process on each subdivision. To avoid infinite loops,
this bails out if too many steps are taken, with a warning. That's
what you're seeing; increasing mxstep will make it go further before
giving up. It won't make a difference in cases where all those steps
aren't necessary.

It's worth thinking about why the integrator is taking all those
steps. It generally happens because there's some sort of kink or
complex behaviour in the solution there; this can be a genuine feature
of the solution, or (more often in my experience) it can be because
your RHS function has some discontinuity (perhaps due to a bug). So
it's worth checking out those segments. I believe that with
full_output you can check how many steps were needed and have your
program bail out with a set of parameters that trigger the problem.

Asking odeint for more intermediate points you don't actually care
about will be a less-efficient way of increasing maxstep - with one
exception. If you know there are problem points in your domain you can
put them in your list to make sure the integrator notices them.
Otherwise, just rely on it to figure out which parts of the
integration need extra work.

Passing the Jacobian in will definitely help reduce the number of
steps needed, if you can compute it analytically. (Otherwise let the
integrator build up a numerical Jacobian rather than trying to do it
yourself.)

Finally, I should say that apart from the stream of ugly messages,
when you hit mxstep what you get is still an estimate of the solution,
just not a very good one (that is, the integrator thinks there's too
much error). If you're okay with crummy solutions in a few corner
cases, you can just live with the warnings.

Anne

On 22 July 2011 01:11, Laura Matrajt <matrajt at gmail.com> wrote:
> Hi Warren,
>  thanks for your fast reply.
> As I said in my email, 90% of the time this runs perfectly well.
> Is 'mxstep' a number that will be computed at each iteration or is it just
> there in case of more steps are needed? I am worried that by increasing it I
> will make my code slower...
> Also, would increasing my tspan solve the issue?
> If the answer is yes, which of the two solutions would be better?
> Also, I am not passing the Jacobian to my function. Do you think that it
> will make a difference, both in terms of speed and in terms of this error to
> pass the Jacobian?
>
> THank you very very much! I was really worried about this!
>
>
> On Thu, Jul 21, 2011 at 6:08 PM, Warren Weckesser
> <warren.weckesser at enthought.com> wrote:
>>
>>
>> On Thu, Jul 21, 2011 at 6:22 PM, Laura Matrajt <matrajt at gmail.com> wrote:
>>>
>>> Hi all,
>>>  I am working with a system of 16 differential equations that simulates
>>> an epidemic in a city.  Because there are many cities interacting with each
>>> other, I need to run my ode's for a single day, stop them, modify the
>>> initial conditions and run them again.  Because the ode is running only for
>>> a day, I defined my tspan to have only two points, like this:
>>>
>>> tspan = tspan = linspace(day, day+1, 2)
>>>
>>> I wrote my equations in Python and I am using scipy.odeint to solve it.
>>>
>>> Here is my code:
>>> def advanceODEoneDay(self,day):
>>>
>>>         #rename variables for convenience
>>>         N0, N1 = self.children, self.adults
>>>         S0,S1,A0,A1,I0,I1,RA0,RA1,RI0,RI1 =
>>> self.S0,self.S1,self.A0,self.A1,self.I0,self.I1,self.RA0,self.RA1,self.RI0,self.RI1
>>>
>>>
>>>
>>>         #create a vector of times for integration:
>>>         tspan = linspace(day, day+1, 2)
>>>
>>>
>>>         #set initial conditions. To do this, I need to look in the array
>>> in the day
>>>         initCond = [S0[day,0], S0[day,1],
>>>                     S1[day,0], S1[day,1],
>>>                     A0[day,0], A0[day,1], A1[day,0], A1[day,1],
>>>                     I0[day,0], I0[day,1], I1[day,0], I1[day,1],
>>>                     RA0[day,0], RA1[day,0], RI0[day,0], RI1[day,0] ]
>>>
>>>
>>>         #run the ode:
>>>         sir_sol = odeint(sir2groups,initCond, tspan,
>>> args=(self.p,self.C,self.m,self.VEs,self.VEi,self.VEp,N0,N1))#,self.gamma,self.rho))
>>>
>>>
>>> Most of the time, it works just fine. However, there are some times where
>>> the following message appears:
>>>
>>> Excess work done on this call (perhaps wrong Dfun type).
>>> Run with full_output = 1 to get quantitative information.
>>>  lsoda--  at current t (=r1), mxstep (=i1) steps
>>>        taken on this call before reaching tout
>>>       In above message,  I1 =       500
>>>       In above message,  R1 =  0.1170754095027E+03
>>>
>>> I should mention that it is NOT only a warning. This is repeated over and
>>> over (thousands of times) and then it will break the rest of my code
>>> Ok, after searching in this mailing list, someone else posted a similar
>>> warning message and it was suggested to him that "
>>>
>>> In your case, you might simply be computing for to coarse a mesh in t,
>>>  so "too much work" has to be done for each step. "
>>>
>>>
>>> Is this what is happening to me? the problem is that I don't get this
>>> error every single time, so I don't even know how to run it with full_output
>>> = 1 to get the info...
>>>
>>>
>>>
>>>
>>> I really don't know what to do. Any help will be very very very very
>>> appreciated!
>>> thank you very
>>
>> The function odeint has the keyword argument 'mxstep' that determines the
>> maximum number of internal steps allowed between requested time values.  The
>> default is 500.  It is not unusual to have to increase this, especially in a
>> case like yours where you simply want the value at the end of a long time
>> interval.  Try increasing it to, say, mxstep=5000.
>>
>> Warren
>>
>>>
>>>
>>>
>>>
>>>
>>>
>>> --
>>> Laura
>>>
>>>
>>>
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>
>>
>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>
>
>
> --
> Laura
>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>