Python使用Kronecker Delta和ODEINT

Question

我试图用Kronecker delta函数来绘制一个ODE的输出，这个函数只应该在一个特定的时间t1变为“active”。 这应该会给出一个锯齿状的响应，其中初始值呈指数衰减，直到t=t1，在再次衰减之前，它会立即再次上升。 然而，当我绘制这个图时，看起来解算器看到Kronecker delta函数在所有时间t都为零。在Python中有这样做的方法吗？你知道吗

from scipy import KroneckerDelta
import scipy.integrate as sp
import matplotlib.pyplot as plt
import numpy as np

def dy_dt(y,t):

    dy_dt = 500*KroneckerDelta(t,t1) - 2y

return dy_dt

t1 = 4
y0 = 500
t = np.arrange(0,10,0.1)

y = sp.odeint(dy_dt,y0,t)

plt.plot(t,y)

Answer 1

在使用时间的简单Kronecker delta的情况下，您可以按如下方式运行ode：

from scipy.integrate import odeint
import matplotlib.pyplot as plt
import numpy as np

def dy_dt(y,t):
    return -2*y

t_delta = 4
tend = 10
y0 = [500]

t1 = np.linspace(0,t_delta,50)
y1 = odeint(dy_dt,y0,t1)

y0 = y1[-1] + 500  # execute Kronecker delta
t2 = np.linspace(t_delta,tend,50)
y2 = odeint(dy_dt,y0,t2)

t = np.append(t1, t2)
y = np.append(y1, y2)
plt.plot(t,y)

对于复杂情况的另一种选择是solve_ivp的events功能。你知道吗

Answer 2

我认为问题可能是内部舍入错误，因为0.1不能准确地表示为pythonfloat。我会努力的

import math

def dy_dt(y,t): 
    if math.isclose(t, t1):
        return 500 - 2*y
    else:
        return -2y

此外，odeint的文档还建议使用args参数而不是全局变量来为派生函数提供访问其他参数的权限，并用np.linspace替换np.arange：

import scipy.integrate as sp 
import matplotlib.pyplot as plt 
import numpy as np 
import math

def dy_dt(y, t, t1): 
    if math.isclose(t, t1):
        return 500 - 2*y
    else:
        return -2*y

t1 = 4 
y0 = 500 
t = np.linspace(0, 10, num=101) 
y = sp.odeint(dy_dt, y0, t, args=(t1,)) 
plt.plot(t, y)

我没有测试代码，所以告诉我是否有什么问题。你知道吗

---

编辑：

在测试我的代码时，我查看了t值，其中dy_dt是为这些值计算的。我注意到odeint不仅使用指定的t值，而且会稍微改变它们：

...
3.6636447422787928
3.743098503914526
3.822552265550259
3.902006027185992
3.991829287543431
4.08165254790087
4.171475808258308
...

现在用我的方法，我们得到

math.isclose(3.991829287543431, 4) # False

因为默认公差设置为最多10^（-9）的相对误差，所以odeint函数会“忽略”4处导数的凹凸。幸运的是，我们可以通过指定更高的错误阈值来解决这个问题：

def dy_dt(y, t, t1): 
    if math.isclose(t, t1, abs_tol=0.01):
        return 500 - 2*y
    else:
        return -2*y

现在dy_dt对于3.99和4.01之间的所有值都非常高。如果linspace的num参数增加，则可以使该范围变小。你知道吗

---

TL；DR

你的问题不是python的问题，而是一个数值求解微分方程的问题：你需要在足够长的区间内改变导数，否则解算器很可能会错过有趣的地方。kronecker delta不适用于求解ODEs的数值方法。你知道吗