【发布时间】:2022-03-02 00:31:25
【问题描述】:
我正在尝试使用 Python 研究以下延迟微分方程的行为:
y''(t) = -y(t)/τ^2 - 2y'(t)/τ - Nd*f(y(t-T))/τ^2,
其中f是一个截止函数,当其参数的绝对值在1到10之间,否则等于0(见图1),Nd,@ 987654328@ 和 T 是常量。
为此,我使用了 JiTCDDE 包。这为上述方程提供了一个合理的解。然而,当我尝试在等式的右侧添加噪声时,我得到的解在几次振荡后稳定到非零常数。这不是方程的数学解(唯一可能的常数解等于零)。我不明白为什么会出现这个问题以及是否可以解决它。
我在下面重现我的代码。在这里,为了简单起见,我用高频余弦代替了噪声,该余弦被引入方程系统作为虚拟变量的初始条件(余弦可以直接引入系统,但对于一般噪音,这似乎是不可能的)。为了进一步简化问题,我还删除了涉及f 函数的术语,因为没有它也会出现问题。图2显示了代码给出的函数图。
from jitcdde import jitcdde, y, t
import numpy as np
from matplotlib import pyplot as plt
import math
from chspy import CubicHermiteSpline
# Definition of function f:
def functionf(x):
return x/4*(1+symengine.erf(x**2-Bmin**2))*(1-symengine.erf(x**2-Bmax**2))
#parameters:
τ = 42.9
T = 35.33
Nd = 8.32
# Definition of the initial conditions:
dt = .01 # Time step.
totT = 10000. # Total time.
Nmax = int(totT / dt) # Number of time steps.
Vt = np.linspace(0., totT, Nmax) # Vector of times.
# Definition of the "noise"
X = np.zeros(Nmax)
for i in range(Nmax):
X[i]=math.cos(Vt[i])
past=CubicHermiteSpline(n=3)
for time, datum in zip(Vt,X):
regular_past = [10.,0.]
past.append((
time-totT,
np.hstack((regular_past,datum)),
np.zeros(3)
))
noise= lambda t: y(2,t-totT)
# Integration of the DDE
g = [
y(1),
-y(0)/τ**2-2*y(1)/τ+0.008*noise(t)
]
g.append(0)
DDE = jitcdde(g)
DDE.add_past_points(past)
DDE.adjust_diff()
data = []
for time in np.arange(DDE.t, DDE.t+totT, 1):
data.append( DDE.integrate(time)[0] )
plt.plot(data)
plt.show()
顺便说一句,我注意到即使没有噪音,解在零点似乎也是不连续的(对于负时间,y 设置为零),我不明白为什么。
【问题讨论】:
-
您的示例代码没有按原样运行。当用一些有根据的猜测来解决这个问题时,我得到了与你相似的结果,但收敛到零。如果我用简单的余弦替换输入法,并且如果我将相应的 ODE 仅与
scipy.integrate.ode集成,则此结果仍然存在,这强烈表明我观察到的是系统的实际动态。但是,收敛到另一个值正是您的问题所在。 -
另外,我最近添加了an interface to add inputs 和callbacks。两者都允许您处理输入噪声,而无需摆弄虚拟动态变量以及随之而来的所有错误。
-
嗨@Wrzlprmft,谢谢你的cmets。是的,代码中有错字,对此感到抱歉,现在应该修复它。是的,在这个简化版本中,问题似乎是收敛到错误的值。实际上,我为原始方程(我的帖子开头的那个)获得了一个非常相似的图,而预计会出现振荡(这可以在 Mathematica 的帮助下进行验证)。我只是认为修复更简单的情况将有助于解决最初的问题。我会尝试使用新界面直接添加噪音,我会告诉你在这种情况下是否有效。
-
@Wrzlprmft 顺便问一下,您知道为什么解决方案似乎在 0 处不连续吗?
-
请注意,就您的问题而言,我的第一条评论完全适用:这似乎是系统的实际动态,正如 ODE 模拟所证明的那样,具有不同模块产生相同结果。
标签: python-3.x differential-equations numerical-stability jitcode-jitcdde-jitcsde