Numerical Solution To A Differential Equation Containing A Dirac Delta Function

I am trying to use scipy to numerically solve the following differential equation x''+x=\sum_{k=1}^{20}\delta(t-k\pi), y(0)=y'(0)=0. Here is the code from scipy.integrate import o

Solution 1:

Dirac delta is not a function. Writing it as density in an integral is still only a symbolic representation. It is, as mathematical object, a functional on the space of continuous functions. delta(t0,f)=f(t0), not more, not less.

One can approximate the evaluation, or "sifting" effect of the delta operator by continuous functions. The usual good approximations have the form N*phi(N*t) where N is a large number and phi a non-negative function, usually with a somewhat compact shape, that has integral one. Popular examples are box functions, tent functions, the Gauß bell curve, ... So you could take

def tentfunc(t): return max(0,1-abs(t))
N = 10.0
def rhs(t): return sum( N*tentfunc(N*(t-(i+1)*np.pi)) for i in range(20))

X0 = [0, 0]
t = np.linspace(0, 80, 1000)

sol = odeint(lambda x,t: [ x[1], rhs(t)-x[0]], X0, t, tcrit=np.pi*np.arange(21), atol=1e-8, rtol=1e-10)

x,v = sol.T
plt.plot(t,x, t, v)

which gives

Note that the density of the t array also influences the accuracy, while the tcrit critical points did not do much.

---

Another way is to remember that delta is the second derivative of max(0,x), so one can construct a function that is the twice primitive of the right side,

defu(t): returnsum(np.maximum(0,t-(i+1)*np.pi) for i inrange(20))

so that now the equation is equivalent to

(x(t)-u(t))'' + x(t) = 0

set y = x-u then

y''(t) + y(t) = -u(t)

which now has a continuous right side.

X0 = [0, 0]
t = np.linspace(0, 80, 1000)

sol = odeint(lambda y,t: [ y[1], -u(t)-y[0]], X0, t, atol=1e-8, rtol=1e-10)

y,v = sol.T
x=y+u(t)
plt.plot(t,x)

Solution 2:

odeint :

1. does not handle sympy symbolic objects

2. it's unlikely it can ever handle Dirac Delta terms.

The best bet is probably to turn dirac deltas into boundary conditions: assume that the function is continuous at the location of the Dirac delta, but the first derivative jumps. Integrating over infinitesimal interval around the location of the delta function gives you the boundary condition for the derivative just left and just right from the delta.