使用 JiTCDDE 的意外解决方案
Unexpected solution using JiTCDDE
我正在尝试使用 Python 研究以下时滞微分方程的行为:
y''(t) = -y(t)/τ^2 - 2y'(t)/τ - Nd*f(y(t-T))/τ^2,
其中 f
是一个截止函数,当其参数的绝对值在 1 和 10 之间时,它基本上等于恒等式,否则等于 0(见图 1),并且 Nd
、τ
和 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 设置为等于零),我不明白为什么。
随着评论的揭晓,你的问题最终归结为:
step_on_discontinuities
假定相对于积分时间的延迟很小,并执行放置在延迟分量指向积分开始的那些时间的步骤(在您的情况下为 0
)。这样 initial discontinuities 就处理好了。
但是,使用延迟的虚拟变量实现输入会给系统带来很大的延迟,在您的情况下 totT
。
step_on_discontinuities
的相应步骤将在 totT
本身,即在所需的积分时间之后。
因此,当您在代码中达到 for time in np.arange(DDE.t, DDE.t+totT, 1):
时,DDE.t
就是 totT
。
因此,在你真正开始整合和观察之前,你已经迈出了一大步,这看起来像是一个不连续的过程,会导致奇怪的结果,特别是你看不到你输入的效果,因为它已经在这一点上“结束”了。
为避免这种情况,请使用 adjust_diff
或 integrate_blindly
而不是 step_on_discontinuities
.
我正在尝试使用 Python 研究以下时滞微分方程的行为:
y''(t) = -y(t)/τ^2 - 2y'(t)/τ - Nd*f(y(t-T))/τ^2,
其中 f
是一个截止函数,当其参数的绝对值在 1 和 10 之间时,它基本上等于恒等式,否则等于 0(见图 1),并且 Nd
、τ
和 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 设置为等于零),我不明白为什么。
随着评论的揭晓,你的问题最终归结为:
step_on_discontinuities
假定相对于积分时间的延迟很小,并执行放置在延迟分量指向积分开始的那些时间的步骤(在您的情况下为 0
)。这样 initial discontinuities 就处理好了。
但是,使用延迟的虚拟变量实现输入会给系统带来很大的延迟,在您的情况下 totT
。
step_on_discontinuities
的相应步骤将在 totT
本身,即在所需的积分时间之后。
因此,当您在代码中达到 for time in np.arange(DDE.t, DDE.t+totT, 1):
时,DDE.t
就是 totT
。
因此,在你真正开始整合和观察之前,你已经迈出了一大步,这看起来像是一个不连续的过程,会导致奇怪的结果,特别是你看不到你输入的效果,因为它已经在这一点上“结束”了。
为避免这种情况,请使用 adjust_diff
或 integrate_blindly
而不是 step_on_discontinuities
.