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)