5

Python を使用して、次の遅延微分方程式の動作を調査しようとしています。

y''(t) = -y(t)/τ^2 - 2y'(t)/τ - Nd*f(y(t-T))/τ^2,

ここfで、 はカットオフ関数であり、その引数の絶対値が 1 から 10 の間であり、それ以外の場合は 0 に等しい場合 (図 1 を参照)、およびNdτおよびTは定数です。

図 1: 関数 f のプロット

このために、JiTCDDE パッケージを使用しています。これにより、上記の方程式の妥当な解が得られます。それにもかかわらず、方程式の右側にノイズを追加しようとすると、数回の振動の後、ゼロ以外の定数に安定する解が得られます。これは方程式の数学的な解ではありません (ゼロに等しい唯一の可能な定数解)。この問題が発生する理由と、解決できるかどうかはわかりません。

以下にコードを再現します。ここでは、簡単にするために、ノイズを高周波コサインで置き換えました。これは、ダミー変数の初期条件として連立方程式に導入されます (コサインはシステムに直接導入することもできますが、一般的なノイズでは、これは不可能と思われます)。問題をさらに単純化するために、f関数がなくても問題が発生するため、関数に関する用語も削除しました。図 2 は、コードによって与えられた関数のプロットを示しています。

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

4

1 に答える 1