# Python - Using A Kronecker Delta With ODEINT

## 10 October 2019 - 1 answer

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 `np.arange` by `np.linspace`:

``````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.

EDIT:

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
``````

Now `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.

TL;DR

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.