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)
print(yInt)
Outputs:
[[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!
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)