title: Have I implemented Milstein's method/Euler-Maruyama correctly?
body: <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>