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

<pre><code>def a_p(t, p, q):
    return -(Gamma0 - Omega0*eta*q**2)*p - Omega0**2*q
</code></pre>

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

<p><a href="https://i.stack.imgur.com/4gCdX.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/4gCdX.png" alt="One possible solution"></a></p>

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

<hr>

<p>There is also a derivative-free version of Milsteins method as a two-stage kind-of Runge-Kutta method, documented in <a href="https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_method_%28SDE%29" rel="nofollow noreferrer">wikipedia</a> or the original in <a href="https://arxiv.org/pdf/1210.0933.pdf" rel="nofollow noreferrer">arxiv.org (PDF)</a>.</p>

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

<pre><code>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)
</code></pre>
Have I implemented Milstein's method/Euler-Maruyama correctly?
<p>I have an stochastic differential equation (SDE) that I am trying to solve using Milsteins method but am getting results that disagree with experiment. </p>

<p>The SDE is</p>

<p><img src="https://latex.codecogs.com/gif.latex?$$&space;%5Cdfrac%7Bd%5E2q(t)%7D%7Bdt%5E2%7D&space;&plus;&space;(%5CGamma_0&space;-&space;%5COmega_0&space;%5Ceta&space;q(t)%5E2)%5Cdfrac%7Bdq(t)%7D%7Bdt%7D&space;&plus;&space;%5COmega_0%5E2&space;q(t)&space;-&space;%5Csqrt%7B%5Cdfrac%7B2%5CGamma_0&space;k_B&space;T_0%7D%7Bm%7D%7D&space;%5Cdfrac%7BdW(t)%7D%7Bdt%7D&space;=&space;0&space;$$" title="$$ \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 $$" /></p>

<p>which I have broken up into 2 first order equations:</p>

<p>eq1:  <img src="https://latex.codecogs.com/gif.latex?dq&amp;=p%5C,dt%5C%5C" title="dq&amp;=p\,dt\\" /></p>

<p>eq2:  <img src="https://latex.codecogs.com/gif.latex?dp&space;=&space;-(%5CGamma_0&space;-&space;%5COmega_0&space;%5Ceta&space;q(t)%5E2)p~dt&space;-&space;%5COmega_0%5E2&space;q(t)~dt&space;&plus;&space;%5Csqrt%7B%5Cdfrac%7B2%5CGamma_0&space;k_B&space;T_0%7D%7Bm%7D%7D&space;dW" title="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" /></p>

<p>Then I have used the Ito form:</p>

<p><img src="https://latex.codecogs.com/gif.latex?$%7B%5Cmathrm&space;%7Bd%7D%7DX_%7Bt%7D=a(X_%7Bt%7D)%5C,%7B%5Cmathrm&space;%7Bd%7D%7Dt&plus;b(X_%7Bt%7D)%5C,%7B%5Cmathrm&space;%7Bd%7D%7DW_%7Bt%7D$" title="${\mathrm {d}}X_{t}=a(X_{t})\,{\mathrm {d}}t+b(X_{t})\,{\mathrm {d}}W_{t}$" /></p>

<p>So that for eq1: </p>

<p><img src="https://latex.codecogs.com/gif.latex?%5C%5C&space;$$&space;a(q_t)&space;=&space;p&space;$$&space;%5C%5C&space;$$&space;b(q_t)&space;=&space;0&space;$$&space;%5C%5C" title="\\ $$ a(q_t) = p $$ \\ $$ b(q_t) = 0 $$ \\" /></p>

<p>and for eq2:</p>

<p><img src="https://latex.codecogs.com/gif.latex?%5C%5C&space;$$a(p_t)&space;=&space;-(%5CGamma_0-%5COmega_0%5Ceta&space;q(t)%5E2)p$$&space;%5C%5C&space;$$b(p_t)&space;=&space;%5Csqrt%7B%5Cdfrac%7B2%5CGamma_0&space;k_B&space;T_0%7Dm%7D$$" title="\\ $$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}$$" /></p>

<p>My python code used to attempt to solve this is like so:</p>

<pre><code># 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 =&gt; 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
</code></pre>

<p>Where in this case p is velocity and q is position.</p>

<p>I then get the following plots of q and p:</p>

<p><a href="https://i.stack.imgur.com/OYQ5X.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/OYQ5X.png" alt="p plotted with time"></a></p>

<p><a href="https://i.stack.imgur.com/9VOdB.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/9VOdB.png" alt="q plotted with time"></a></p>

<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):</p>

<p><a href="https://i.stack.imgur.com/BDQRY.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/BDQRY.png" alt="experimental position with time"></a></p>

<p>Have I implemented Milstein's method correctly?</p>

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