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

Question

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

SDE 是

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

方程1：

方程2：

然后我使用了 Ito 形式：

所以对于 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 的图：

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

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

如果我有，那么我解决 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 方法，记录在维基百科或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)