Python - Using A Kronecker Delta With ODEINT
I'm trying to plot the output from an ODE using a Kronecker delta function which should only become 'active' at a specific time = t1. This should give a sawtooth like response where the initial value decays down exponentially until t=t1 where it rises again instantly before decaying down once again. However, when I plot this it looks like the solver is seeing the Kronecker delta function as zero for all time t. Is there anyway to do this in Python?
from scipy import KroneckerDelta import scipy.integrate as sp import matplotlib.pyplot as plt import numpy as np def dy_dt(y,t): dy_dt = 500*KroneckerDelta(t,t1) - 2y return dy_dt t1 = 4 y0 = 500 t = np.arrange(0,10,0.1) y = sp.odeint(dy_dt,y0,t) plt.plot(t,y)
I think the problem could be internal rounding errors, because 0.1 cannot be represented exactly as a python
float. I would try
import math def dy_dt(y,t): if math.isclose(t, t1): return 500 - 2*y else: return -2y
Also the documentation of
odeint suggests using the
args parameter instead of global variables to give your derivative function access to additional arguments and replacing
import scipy.integrate as sp import matplotlib.pyplot as plt import numpy as np import math def dy_dt(y, t, t1): if math.isclose(t, t1): return 500 - 2*y else: return -2*y t1 = 4 y0 = 500 t = np.linspace(0, 10, num=101) y = sp.odeint(dy_dt, y0, t, args=(t1,)) plt.plot(t, y)
I did not test the code so tell me if there is anything wrong with it.
When testing my code I took a look at the
t values for which
dy_dt is evaluated. I noticed that
odeint does not only use the
t values that where specified, but alters them slightly:
... 3.6636447422787928 3.743098503914526 3.822552265550259 3.902006027185992 3.991829287543431 4.08165254790087 4.171475808258308 ...
Now using my method, we get
math.isclose(3.991829287543431, 4) # False
because the default tolerance is set to a relative error of at most 10^(-9), so the
odeint function "misses" the bump of the derivative at 4. Luckily, we can fix that by specifying a higher error threshold:
def dy_dt(y, t, t1): if math.isclose(t, t1, abs_tol=0.01): return 500 - 2*y else: return -2*y
dy_dt is very high for all values between 3.99 and 4.01. It is possible to make this range smaller if the
num argument of
linspace is increased.
Your problem is not a problem of python but a problem of numerically solving an differential equation: You need to alter your derivative for an interval of sufficient length, otherwise the solver will likely miss the interesting spot. A kronecker delta does not work with numeric approaches to solving ODEs.
- → What are the pluses/minuses of different ways to configure GPIOs on the Beaglebone Black?
- → Django, code inside <script> tag doesn't work in a template
- → React - Django webpack config with dynamic 'output'
- → GAE Python app - Does URL matter for SEO?
- → Put a Rendered Django Template in Json along with some other items
- → session disappears when request is sent from fetch
- → Python Shopify API output formatted datetime string in django template
- → Shopify app: adding a new shipping address via webhook
- → Shopify + Python library: how to create new shipping address
- → shopify python api: how do add new assets to published theme?
- → Access 'HTTP_X_SHOPIFY_SHOP_API_CALL_LIMIT' with Python Shopify Module