0

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.

xpaul
  • 77
  • 7

2 Answers2

2

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

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
1

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.

ev-br
  • 24,968
  • 9
  • 65
  • 78
  • Thanks. I tried to redefine the Dirac delta function by this: if |t|<0.02, return 1, else return 0. It still doesn't work. – xpaul Aug 03 '20 at 20:12
  • I don't think it's good direction. odeint necessarily has some assumptions about coefficients being smooth, and delta surely isn't. – ev-br Aug 03 '20 at 20:20
  • Do you have any idea to fix the problem? What else can I do? – xpaul Aug 04 '20 at 13:21
  • Like I said, I'd turn it into the linear algebra problem. Search for Kronig-Penney model. – ev-br Aug 04 '20 at 22:51
  • I just wanted to use python to numerical solve the IVP with impulses. – xpaul Aug 05 '20 at 20:47