Clever Geek Handbook

Did I follow the Milstein / Euler-Maruyama method correctly?

I have a Stochastic Differential Equation (SDE) that I am trying to solve using the Milstein method, but I get results that are not consistent with experiment.

SDE

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% 7D & space;% 5Cdfrac% 7BdW (t)% 7D% 7Bdt% 7D & space; = & space; 0 & space; $$

which I split into 2 first order equations:

eq1: data-src =

eq2: 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

Then I used the Ito form:

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 $

So for eq1:

data-src =

and for eq2:

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

My python code used to try to solve this looks like this:

# 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 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 speed and q is position.

Then I get the following q and p plots:

p plotted over time

q built over time

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

experimental position over time

Did I follow Milstein's method correctly?

If I have what else might be wrong, my SDE decision process causing this disagreement with the experiment?

+3


source to share


1 answer


You missed a term in the drift coefficient, you will notice that dp

there are two terms to the right of this dt

. Thus,

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

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

One possible solution

And no, you have not implemented Milstein's method, since there are no derivatives of b_p

that distinguish Milstein from Euler-Maruyama, the missing term +0.5*b'(X)*b(X)*(dW**2-dt)

.



There is also a derivative version of the Milstein method as a two-step form of the Runge-Kutta method, documented in wikipedia or original at arxiv.org (PDF) .

This step is (vector based, duplicate on X=[p,q]

, K1=[k1_p,k1_q]

etc. to be closer 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)
+3


source






More articles:


All Articles