python numerical-methods differential-equations stochastic

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

I have an stochastic differential equation (SDE) that I am trying to solve using Milsteins method but am getting results that disagree with experiment.

The SDE is

$\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$

which I have broken up into 2 first order equations:

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$

Then I have used the Ito form:

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

So that for eq1:

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

and for 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}$$$

My python code used to attempt to solve this is like so:

# 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

Where in this case p is velocity and q is position.

I then get the following plots of q and p:

I expected the resulting plot of position to look something like the following, which I get from experimental data (from which the constants used in the model are determined):

Have I implemented Milstein's method correctly?

If I have, what else might be wrong my process of solving the SDE that'd causing this disagreement with the experiment?

Solution

You missed a term in the drift coefficient, note that to the right of dp there are two dt terms. Thus

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

which is actually the part that makes the oscillator into an oscillator. With that corrected the solution looks like

And no, you did not implement the Milstein method as there are no derivatives of b_p which are what distinguishes Milstein from Euler-Maruyama, the missing term is +0.5*b'(X)*b(X)*(dW**2-dt).

There is also a derivative-free version of Milsteins method as a two-stage kind-of Runge-Kutta method, documented in wikipedia or the original in arxiv.org (PDF).

The step there is (vector based, duplicate into X=[p,q], K1=[k1_p,k1_q] etc. to be close to your conventions)

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)