【问题标题】:Unexpected solution using JiTCDDE使用 JiTCDDE 的意外解决方案
【发布时间】: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 inputscallbacks。两者都允许您处理输入噪声,而无需摆弄虚拟动态变量以及随之而来的所有错误。
  • 嗨@Wrzlprmft,谢谢你的cmets。是的,代码中有错字,对此感到抱歉,现在应该修复它。是的,在这个简化版本中,问题似乎是收敛到错误的值。实际上,我为原始方程(我的帖子开头的那个)获得了一个非常相似的图,而预计会出现振荡(这可以在 Mathematica 的帮助下进行验证)。我只是认为修复更简单的情况将有助于解决最初的问题。我会尝试使用新界面直接添加噪音,我会告诉你在这种情况下是否有效。
  • @Wrzlprmft 顺便问一下,您知道为什么解决方案似乎在 0 处不连续吗?
  • 请注意,就您的问题而言,我的第一条评论完全适用:这似乎是系统的实际动态,正如 ODE 模拟所证明的那样,具有不同模块产生相同结果。

标签: python-3.x differential-equations numerical-stability jitcode-jitcdde-jitcsde


【解决方案1】:

随着 cmets 的揭幕,您的问题最终归结为:

step_on_discontinuities 假设延迟相对于积分时间来说很小,并执行放置在延迟组件指向积分开始的那些时间的步骤(在您的情况下为 0)。这样initial discontinuities就被处理了。

但是,使用延迟的虚拟变量实现输入会在系统中引入很大的延迟,在您的情况下为 totTstep_on_discontinuities 的相应步骤将在 totT 本身,即在所需的集成时间之后。 因此,当您在代码中达到for time in np.arange(DDE.t, DDE.t+totT, 1): 时,DDE.t 就是totT。 因此,在您真正开始整合和观察之前,您已经迈出了一大步,这可能看起来像是一个不连续性并导致奇怪的结果,特别是您看不到输入的效果,因为此时它已经“结束”了。 为避免这种情况,请使用adjust_diffintegrate_blindly 而不是step_on_discontinuities

【讨论】:

  • 非常感谢您发现问题!我一定会接受这个答案。我可以问你解决它的最佳方法是什么吗?我应该在 step_on_discontinuities 中设置一个小的最大步长吗?还是我应该只使用新界面添加噪音,以避免大延迟totT?或者你会推荐使用adjust_diff?
  • 查看我的编辑。即使使用jitcdde_input,也不应使用step_on_discontinuities。事实上,这个问题提醒我JiTCDDE应该在这种情况下抛出错误,it now does