python - 我是否正确实现了 Milstein 的方法/Euler-Maruyama？

Question

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

SDE是

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

eq1:

eq2:

然后我用了伊藤形式:

所以对于 eq1:

对于 eq2:

我用来尝试解决这个问题的 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 的图:

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

我是否正确实现了 Milstein 的方法？

如果我有，我解决 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 中。或 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)