Python

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

我认为问题可能是内部舍入错误，因为 0.1 不能精确表示为 python float。我会尝试

import math

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

odeint 的文档还建议使用 args 参数而不是全局变量来让您的导数函数访问其他参数，并将 np.arange 替换为 np.linspace：

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)

我没有测试代码，如果有什么问题请告诉我。

编辑：

在测试我的代码时，我查看了评估 dy_dt 的 t 值。我注意到 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 函数 "misses" 导数的凸起为 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 不适用于求解 ODE 的数值方法。

Answer 2

在使用时间的简单 Kronecker delta 的情况下，您可以运行像这样分段的颂歌：

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 功能。

Python - 将 Kronecker Delta 与 ODEINT 结合使用

Python - Using a Kronecker Delta with ODEINT

delta

odeint