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 は負の時間でゼロに等しく設定されている) ことに気付きましたが、その理由はわかりません。