I am using odeint in scipy to integrate a function. The function is actually a vector (it's the returned value from the L_total function below). The problem I have is that for some entries of this vector the integration is very easy, meaning that it does not require many internal substeps, however for some other values the integration should be more detailed. I assume that odeint takes one number of substeps and that is the one that is applied to every entry of the vector. If this assumption is correct is there a way to assign a maximum number of internal substeps in scipy in an efficient way depending on the entry in the function vector? I know that by using an IF statement I might be able of changing the number of maximum number of steps allowed but then I'm afraid that it will take more computational time and that is what I am aiming to decrease. My code looks like this:
def L_total(self, r, d, u):
tem = self.T(u)
condlist = [numpy.log10(tem) <= 2, numpy.log10(tem) >= 2]
choicelist = [self.L1(r, d, tem), self.L2(r, tem)]
return numpy.select(condlist, choicelist) #This is the returned value for the scipy integration, i.e. the vector I want to integrate
def integrate_function(self, i_0, time_final, r, d):
def f(i, t):
return - self.L_total(r, d, i)
time = numpy.linspace(0., time_final, 2)
result = odeint(f, i_0, time)
return result[-1]
The last function works fine for a while until the values of the variable tem increase, then I get a message saying 'Excess work done on this call' and the suggestion is to increase the number of timesteps allowed:
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.3319309749049D+03
Excess work done on this call (perhaps wrong Dfun type)
Increasing the number of the argument mxstep sounds like the solution except for the fact that as long as log(tem)< 2 the integration is easy and the number of internal timesteps can be small, I don't want to waste time on this easy integration, the idea is to only increase the mxstep value for certain values of the vector. A way to do it is:
def integrate_function(self, i_0, time_final, r, d, tem):
def f(i, t):
return - self.L_total(r, d, i)
time = numpy.linspace(0., time_final, 2)
if numpy.log10(tem)<=2:
result = odeint(f, i_0, time, mxstep=10)
else:
result = odeint(f, i_0, time, mxstep=1000)
return result[-1]
however the length of the vector function is very large so it will consume lots of computational time since I have to do this loop inside of another loop. What I need is efficiency since my code already takes a lot of time.
Subtimesteps chosen automatically by odeint are the same for all equations in the system. This is because normally, when solving a system of ODEs, one cannot work with the components separately: all of them must be known at the same point of time, in order to find the derivative at that time.
For an uncoupled system, where equations have nothing to do with each other, the above poses a possibly unnecessary constraint on the time steps. In this case solving individual equations in a loop may actually be faster than working with them as a system.