包含 Dirac delta 函数的微分方程的数值解 [英] Numerical solution to a differential equation containing a Dirac delta function

问题描述

我正在尝试使用 scipy 对以下微分方程进行数值求解

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.

这是代码

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 得到的不同.这是数学代码和图表

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

在我看来，scipy 忽略了 Dirac delta 函数.我哪里错了?任何帮助表示赞赏.

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

推荐答案

Dirac delta 不是函数.将其写为积分中的密度仍然只是一种符号表示.作为数学对象，它是连续函数空间上的泛函.delta(t0,f)=f(t0)，不多不少.

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.

人们可以近似评估，或筛选"delta 运算符对连续函数的影响.通常的好近似具有 N*phi(N*t) 的形式，其中 N 是一个大数，而 phi 是一个非负函数，通常有一个有点紧凑的形状，有一个完整的.流行的例子是盒函数、帐篷函数、高斯钟形曲线……所以你可以采取

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)

给出

请注意，t 数组的密度也会影响准确性，而 tcrit 关键点没有太大影响.

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

另一种方法是记住delta是max(0,x)的二阶导数，因此可以构造一个函数，它是右侧的两倍原语，

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

所以现在等价于

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

这篇关于包含 Dirac delta 函数的微分方程的数值解的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！