Численное решение дифференциального уравнения, содержащего дельта-функцию Дира c

Question

Я пытаюсь использовать scipy для численного решения следующего дифференциального уравнения

x''+x=\sum_{k=1}^{20}\delta(t-k\pi), y(0)=y'(0)=0.

Вот код

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()

Однако то, что я получил от python, является нулевым решением, что отличается от того, что я получил от Mathematica. Вот математический код и график

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}]

введите описание изображения здесь

Мне кажется, что scipy игнорирует дельта-функцию Дира c. Где я не прав? Любая помощь приветствуется.

Answer 1

odeint:

1. не обрабатывает симпозиционные символы c объекты

2. маловероятно, что он когда-либо сможет обработать Дира c Дельта-условия.

Лучше всего, вероятно, превратить дельты диры c в граничные условия: предположим, что функция является непрерывной в месте расположения дельты диры c, но скачки первой производной. Интегрирование по бесконечно малому интервалу вокруг местоположения дельта-функции дает вам граничное условие для производной слева и справа от дельты.