I was messing around with some numerical integration functions. I wrote an arbitrary differential equation to test my understanding, the code is as follows:

import numpy as np

from scipy.integrate import odeint

intT = np.array([0,1,2,3,4,6,8,10,13,16,19,20,21,22,27,40,80])

def dydt(y,t):
    dydt = np.random.choice([0,1], 1)[0]
    return dydt

yInt = odeint(dydt, 0, intT)


[[0.00000000e+000] [0.00000000e+000] [0.00000000e+000] [0.00000000e+000] [0.00000000e+000] [0.00000000e+000] [0.00000000e+000] [2.21317648e-001] [3.22131765e+000] [3.62733355e+000] [0.00000000e+000] [2.42092166e-322] [2.70777855e-316] [2.81469831e-316] [0.00000000e+000] [0.00000000e+000] [7.11454530e-322]]

/usr/lib/python3/dist-packages/scipy/integrate/odepack.py:236: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information. warnings.warn(warning_msg, ODEintWarning)

This gave me very weird results, and that was not what I was expecting. If I change return value of dydt into 1, then this simple integration surely just works properly. What I was actually trying to do is to apply a predetermined window function to my differential equation, this code above is a simplified version I would imagine. But I just don't get what is wrong here. I wonder if you have any thoughts.


Thanks to @Lutz Lehmann. I was trying to integrate some breathing function, who spikes every 100 time intervals. I had another shot as follows, this doesn't work either. And I notice this seems to be very tricky, if I change my window function and differential equation, sometimes it works, sometimes it doesn't. I kind of understand the smoothness issue, but I just wonder if there is a way to handle this problem. I initially had a numerical integration script by myself, but it works way slower, that's why I turned to the scipy integrate library.

intT = (np.arange(1, 1000) + 0.1)

def arbWindow(t):
    w = np.exp(1/t - t/10)
    return w
def dydt(y,t):
    dydt = arbWindow(t%100) - 0.01*y**2
    return dydt

yInt = odeint(dydt, 0, intT)

This use of the numerical solver is completely wrong. The numerical ODE solvers are for problems that have a smooth right side. As long as the existence of an exact solution is ensured, they can also be applied to problems where the right side is piecewise smooth, but that will tend to slow down the integration with very small step sizes at the discontinuities.

But here the right side function is not even a function, as the value changes randomly even if the same arguments are given. From inside the algorithm this will look as if the function is highly oscillatory, or discontinuous on all scales, requiring smaller and smaller step sizes without finding a scale where the problem looks smooth at this zoom level. This then leads to the Excess work done on this call error message.

  • $\begingroup$ Thank you, that explains it. $\endgroup$ – contain100pctrecycledfibre Jan 6 at 10:57
  • $\begingroup$ I tried another more physically meaningful function, I wonder if you can take a look. I amended it to my question. Thanks $\endgroup$ – contain100pctrecycledfibre Jan 6 at 13:19
  • $\begingroup$ This new system should work. For such long time intervals, you might need to set the error tolerances. Perhaps you should shift the shift 0.1 directly to the division, w=np.exp(1/(t+0.1)-t/10), as t%100 can produce values very close to zero. $\endgroup$ – Lutz Lehmann Jan 6 at 13:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.