Help!!!!!! having problems with ODEINT
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 * * -- Laura
On Thu, Jul 21, 2011 at 6:22 PM, Laura Matrajt <matrajt@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@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
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@enthought.com> wrote:
On Thu, Jul 21, 2011 at 6:22 PM, Laura Matrajt <matrajt@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@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
-- Laura
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@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@enthought.com> wrote:
On Thu, Jul 21, 2011 at 6:22 PM, Laura Matrajt <matrajt@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@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
-- Laura
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
Anne Archibald <aarchiba <at> physics.mcgill.ca> writes:
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
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,
On Fri, Jul 22, 2011 at 17:49, Laura Matrajt <matrajt@gmail.com> wrote:
Anne Archibald <aarchiba <at> physics.mcgill.ca> writes:
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
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,
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-mult... Hope that helps, Kevin
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-mult...
Hope that helps, Kevin
Hi Kevin, the example was good. I now have implemented my Jacobian. Thank you very much!!!!!
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-mult...
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!
On Mon, Jul 25, 2011 at 17:53, Laura Matrajt <matrajt@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-mult...
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
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!
Hi Laura, On Mon, Jul 25, 2011 at 4:53 PM, Laura Matrajt <matrajt@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-mult...
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
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!
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
Hi Laura,
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
participants (6)
-
Anne Archibald -
bhanukiran perabathini -
Kevin Dunn -
Laura Matrajt -
Rob Clewley -
Warren Weckesser