我是否正确实施了 Milstein 的 method/Euler-Maruyama？

Question

我有一个随机微分方程 (SDE)，我正在尝试使用 Milsteins 方法求解，但得到的结果与实验不符。

SDE 是

$\dfrac{d^2q(t)}{dt^2} + (\Gamma_0 - \Omega_0 \eta q(t)^2)\dfrac{dq(t)}{dt} + \Omega_0^2 q(t) - \sqrt{\dfrac{2\Gamma_0 k_B T_0}{m}} \dfrac{dW(t)}{dt} = 0$

我将其分解为 2 个一阶方程：

eq1: $dq&=p\,dt\$

eq2: $dp = -(\Gamma_0 - \Omega_0 \eta q(t)^2)p~dt - \Omega_0^2 q(t)~dt + \sqrt{\dfrac{2\Gamma_0 k_B T_0}{m}} dW$

那我就用伊藤形式了:

${\mathrm {d}}X_{t}=a(X_{t})\,{\mathrm {d}}t+b(X_{t})\,{\mathrm {d}}W_{t}$

因此对于 eq1：

$\ $$ a(q_t) = p $$ \ $$ b(q_t) = 0 $$ \$

对于 eq2：

$\ $$a(p_t) = -(\Gamma_0-\Omega_0\eta q(t)^2)p$$ \ $$b(p_t) = \sqrt{\dfrac{2\Gamma_0 k_B T_0}m}$$$

我的 python 用于尝试解决此问题的代码如下：

# set constants from real data
Gamma0 = 4000  # defines enviromental damping
Omega0 = 75e3*2*np.pi # defines the angular frequency of the motion
eta = 0 # set eta 0 => no effect from non-linear p*q**2 term
T_0 = 300 # temperature of enviroment
k_b = scipy.constants.Boltzmann 
m = 3.1e-19 # mass of oscillator

# set a and b functions for these 2 equations
def a_p(t, p, q):
    return -(Gamma0 - Omega0*eta*q**2)*p

def b_p(t, p, q):
    return np.sqrt(2*Gamma0*k_b*T_0/m)

def a_q(t, p, q):
    return p

# generate time data
dt = 10e-11
tArray = np.arange(0, 200e-6, dt)

# initialise q and p arrays and set initial conditions to 0, 0
q0 = 0
p0 = 0
q = np.zeros_like(tArray)
p = np.zeros_like(tArray)
q[0] = q0
p[0] = p0

# generate normally distributed random numbers
dwArray = np.random.normal(0, np.sqrt(dt), len(tArray)) # independent and identically distributed normal random variables with expected value 0 and variance dt

# iterate through implementing Milstein's method (technically Euler-Maruyama since b' = 0
for n, t in enumerate(tArray[:-1]):
    dw = dwArray[n]
    p[n+1] = p[n] + a_p(t, p[n], q[n])*dt + b_p(t, p[n], q[n])*dw + 0
    q[n+1] = q[n] + a_q(t, p[n], q[n])*dt + 0

在这种情况下，p 是速度，q 是位置。

然后我得到以下 q 和 p 的图：

我希望得到的位置图看起来像下面这样，这是我从实验数据中得到的（模型中使用的常数是从中确定的）：

我是否正确实施了米尔斯坦的方法？

如果我有，我解决 SDE 的过程可能还有什么问题导致与实验不一致？

Answer 1

您在漂移系数中遗漏了一项，请注意 dp 右侧有两项 dt。于是

def a_p(t, p, q):
    return -(Gamma0 - Omega0*eta*q**2)*p - Omega0**2*q

这其实就是把震荡器做成震荡器的部分。更正后的解决方案看起来像

不，您没有实施 Milstein 方法，因为没有 b_p 的导数，这正是 Milstein 与 Euler-Maruyama 的区别所在，缺少的项是 +0.5*b'(X)*b(X)*(dW**2-dt)。

还有一个 Milsteins 方法的无导数版本，作为一种两阶段的 Runge-Kutta 方法，记录在 wikipedia or the original in arxiv.org (PDF)。

那里的步骤（基于矢量，复制到 X=[p,q]、K1=[k1_p,k1_q] 等以接近您的约定）

S = random_choice_of ([-1,1])
K1 = a(X )*dt + b(X )*(dW - S*sqrt(dt))
Xh = X + K1
K2 = a(Xh)*dt + b(Xh)*(dW + S*sqrt(dt))

X = X + 0.5 * (K1+K2)

我是否正确实施了 Milstein 的 method/Euler-Maruyama？

Have I implemented Milstein's method/Euler-Maruyama correctly?

python

numerical-methods

stochastic

differential-equations