CODEWITHSUNDEEP
pythonscipy
xpaul
xpaul

Reputation: 87

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 odeint
import numpy as np
import matplotlib.pyplot as plt
from sympy import DiracDelta

def f(t):
    sum = 0
  for i in range(20):
     sum = sum + 1.0*DiracDelta(t-(i+1)*np.pi)
        
  return sum

def ode(X, t):
  x = X[0]
  y = X[1]
  dxdt = y
  dydt = -x + f(t)
return [dxdt, dydt]

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

sol = odeint(ode, X0, t)

x = sol[:, 0]
y = sol[:, 1]

plt.plot(t,x, t, y)
plt.xlabel('t')
plt.legend(('x', 'y'))

# phase portrait
plt.figure()
plt.plot(x,y)
plt.plot(x[0], y[0], 'ro')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

However what I got from python is zero solution, which is different from what I got from Mathematica. Here are the mathematica code and the graph

so=NDSolve[{x''(t)+x(t)=\sum _{i=1}^{20} DiraDelta (t-i \pi ),x(0)=0,x'(0)=0},x(t),{t,0,80}]

enter image description here

It seems to me that scipy ignores the Dirac delta function. Where am I wrong? Any help is appreciated.

Upvotes: 1

Views: 2098

Answers (2)

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

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

plot of the solution

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,

def u(t): return sum(np.maximum(0,t-(i+1)*np.pi) for i in range(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)

plot of the second variant

Upvotes: 2

ev-br
ev-br

Reputation: 26040

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.

Upvotes: 1

Related Questions