[SciPy-User] Help!!!!!! having problems with ODEINT
Laura Matrajt
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
--
Laura
Warren Weckesser
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...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user
Laura Matrajt
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!
Laura
Anne Archibald
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
Laura Matrajt
Hi Anne,
thank you very very much for your response.I appreciate your thorough response.
I have one last question that is driving me crazy: why is this crashing my code
eventhough it is just a warning?
Also, regarding the Jacobian, I searched the documentation online and have one
little question:
suppose my system of odes looks like
dy/dt = f(y,t)
the Jacobian is only with respect to y?
I can compute it analytically (I believe) but I am unsure about the correct
syntax to pass it. Do you by any chance about a worked example?
I am really sorry I am probably doing dumb questions. I am pretty overwhelmed by
this right now and i searched in the internet for examples and I couldn't find
one with the Jacobian on it.
THanks again,
Kevin Dunn
Hi Laura,
There's an example in the SciPy documentation on using a Jacobian:
http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html
This is for the ``ode`` function, not the ``odeint`` function though.
I've just posted an example on using ``ode`` with coupled ODEs, no Jacobian
over here:
http://scipy-central.org/item/13/0/integrating-an-initial-value-problem-multiple-odes
Hope that helps,
Kevin
Laura Matrajt
> Hi Laura,
>
> There's an example in the SciPy documentation on using a Jacobian:
> http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html
>
> This is for the ``ode`` function, not the ``odeint`` function though.
>
> I've just posted an example on using ``ode`` with coupled ODEs, no Jacobian
> over here:
>
http://scipy-central.org/item/13/0/integrating-an-initial-value-problem-multiple-odes
>
> Hope that helps,
> Kevin
>
Hi Kevin,
the example was good. I now have implemented my Jacobian. Thank you very much!!!!!
Laura Matrajt
>
> Kevin Dunn <kgdunn <at> gmail.com> writes:
>
> > Hi Laura,
> >
> > There's an example in the SciPy documentation on using a Jacobian:
> > http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html
> >
> > This is for the ``ode`` function, not the ``odeint`` function though.
> >
> > I've just posted an example on using ``ode`` with coupled ODEs, no Jacobian
> > over here:
> >
>
http://scipy-central.org/item/13/0/integrating-an-initial-value-problem-multiple-odes
> >
> > Hope that helps,
> > Kevin
> >
>
> Hi Kevin,
> the example was good. I now have implemented my Jacobian. Thank you very
much!!!!!
>
Hi all:
sorry to bother you again! I implemented the Jacobian as Anne kindly suggested
to me with the help of Kevin's pointers to the correct webpage AND I increased
the maximum number of steps as Warren kindly said.
I am now getting a new message:
lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
In above, R1 = 0.1209062893646E+03 R2 = 0.9059171791494E-18
It is just a warning, and my code continues to run. But I am really worried
about this being a bug. The problem is that I am coupling a system of ODE's with
a stochastic process. Mainly, I simulate a day of an epidemic, stop the
integrator, change some of the initial conditions stochastically (not just
randomly, I do follow some rules) and I run the ODE again and so on.
I have run this millions of times (and I am not exagerating about the millions)
and it doesn't produce any warnings, but every now and then (~15 times) it does
it. I don't know if this is a bug or just that not all of my domain will be good
for the ODE's...
If anyone has any suggestion of how to think about this problem, I will really
appreciate it!!!!!!
thanks to all the people that replied to me previously, you helped me sooo much
already!
Kevin Dunn
> Laura Matrajt <matrajt <at> gmail.com> writes:
>
>>
>> Kevin Dunn <kgdunn <at> gmail.com> writes:
>>
>> > Hi Laura,
>> >
>> > There's an example in the SciPy documentation on using a Jacobian:
>> > http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html
>> >
>> > This is for the ``ode`` function, not the ``odeint`` function though.
>> >
>> > I've just posted an example on using ``ode`` with coupled ODEs, no Jacobian
>> > over here:
>> >
>>
> http://scipy-central.org/item/13/0/integrating-an-initial-value-problem-multiple-odes
>> >
>> > Hope that helps,
>> > Kevin
>> >
>>
>> Hi Kevin,
>> the example was good. I now have implemented my Jacobian. Thank you very
> much!!!!!
>>
>
>
> Hi all:
> sorry to bother you again! I implemented the Jacobian as Anne kindly suggested
> to me with the help of Kevin's pointers to the correct webpage AND I increased
> the maximum number of steps as Warren kindly said.
> I am now getting a new message:
>
> lsoda-- warning..internal t (=r1) and h (=r2) are
> such that in the machine, t + h = t on the next step
> (h = step size). solver will continue anyway
> In above, R1 = 0.1209062893646E+03 R2 = 0.9059171791494E-18
I would find the set(s) of initial conditions that give this warning,
and repeat the integration with a different integrator and/or
different tolerance settings on the integrator. Then compare the
trajectories. If there is negligible difference, then can probably
ignore the warning.
That being said, I've not used the lsoda integrator before, so I have
no idea how serious this warning might be. Which is why I recommend
you try a different integrator also.
Kevin
Warren Weckesser
On Mon, Jul 25, 2011 at 4:53 PM, Laura Matrajt <mat...@gmail.com> wrote:
> Laura Matrajt <matrajt <at> gmail.com> writes:
>
>>
>> Kevin Dunn <kgdunn <at> gmail.com> writes:
>>
>> > Hi Laura,
>> >
>> > There's an example in the SciPy documentation on using a Jacobian:
>> > http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html
>> >
>> > This is for the ``ode`` function, not the ``odeint`` function though.
>> >
>> > I've just posted an example on using ``ode`` with coupled ODEs, no Jacobian
>> > over here:
>> >
>>
> http://scipy-central.org/item/13/0/integrating-an-initial-value-problem-multiple-odes
>> >
>> > Hope that helps,
>> > Kevin
>> >
>>
>> Hi Kevin,
>> the example was good. I now have implemented my Jacobian. Thank you very
> much!!!!!
>>
>
>
> Hi all:
> sorry to bother you again! I implemented the Jacobian as Anne kindly suggested
> to me with the help of Kevin's pointers to the correct webpage AND I increased
> the maximum number of steps as Warren kindly said.
> I am now getting a new message:
>
> lsoda-- warning..internal t (=r1) and h (=r2) are
> such that in the machine, t + h = t on the next step
> (h = step size). solver will continue anyway
> In above, R1 = 0.1209062893646E+03 R2 = 0.9059171791494E-18
Something about your system is causing the solver to reduce its
internal step size down to about 1e-17 (and it can't go any smaller
than that). Do you actually have a discontinuity in your equations?
Is your system singularly perturbed, with shock-like or boundary layer
solutions? As Anne said: "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."
If you can isolate the initial conditions and parameters that lead to
this warning, you could plot the solution that it generates to see
what is going on.
Warren
bhanukiran perabathini
Rob Clewley
>> Hi all:
>> sorry to bother you again! I implemented the Jacobian as Anne kindly suggested
>> to me with the help of Kevin's pointers to the correct webpage AND I increased
>> the maximum number of steps as Warren kindly said.
>> I am now getting a new message:
>>
>> lsoda-- warning..internal t (=r1) and h (=r2) are
>> such that in the machine, t + h = t on the next step
>> (h = step size). solver will continue anyway
>> In above, R1 = 0.1209062893646E+03 R2 = 0.9059171791494E-18
>
The other recommendations are a good place to start fixing this,
except the one suggesting you unsubscribe :) But, short of finding an
implementation of a true stochastic DE solver (which I'm afraid I
can't help with as I'm not a big expert) you should find it easier to
introduce specific noise signals to your system if you use PyDSTool's
external input signals that can appear in the RHS function of your DE.
Alternatively, you could run your problem as a "hybrid" model where a
daily "event" in your code will cause the discrete state transition to
the new values drawn from the noise distribution.
I am not sure if this will fix your problem with the step size going
to zero, but that will depend on your parameters and exactly how
you've implemented the running of the system with the discontinuity in
state values (and how large they are). But since you are doing this
step millions of times, PyDSTool should speed up your runs
significantly, especially if you also use a C-based integrator that
will compile your RHS too.
If you decide to try this out there is an example interp_vode_test.py
in the package's tests directory that uses external input signals, and
others that demonstrate hybrid systems. But I'd be willing to help you
get it working for you if you send me your code, as using this to
simulate systems with noise is not yet a well-documented feature.
Best,
Rob