Simulate Stochastic Differential equation in Python

Question

I am trying to solve SDE for Brownian particle and Langevein Dynamics. At first I tried to simulate 2D brownian motion with normal random number generator, The code is:

import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline
dt = .001  # Time step.
T = 2.  # Total time.
n = int(T / dt)  # Number of time steps.
t = np.linspace(0., T, n)  # Vector of times.
sqrtdt = np.sqrt(dt)
y = np.zeros(n)
x = np.zeros(n)


for i in range(n-1):
  x[i + 1] = x[i] +  np.random.normal(0.0,1.0)
  y[i + 1] = y[i] +  np.random.normal(0.0,1.0)


fig, axs = plt.subplots(1, 1, figsize=(12, 12))
plt.plot(y, x, label ='Position')
plt.title("Simulation of Brownian motion") 
plt.show()

Now when I am trying to simulate the same process with the help of forward Euler Method, the governing equation is

m^dv/_dt=η

using the following code,

import numpy as np
import matplotlib.pyplot as plt


%matplotlib inline
dt = .001  # Time step.
T = 2.  # Total time.
n = int(T / dt)  # Number of time steps.
t = np.linspace(0., T, n)  # Vector of times.
sqrtdt = np.sqrt(dt)
v_x = np.zeros(n)
v_y = np.zeros(n)

y = np.zeros(n)
x = np.zeros(n)
for i in range(n-1):
  v_x[i + 1] = v_x[i] +  sqrtdt * np.random.normal(0.0,1.0)
  v_y[i + 1] = v_y[i] +  sqrtdt * np.random.normal(0.0,1.0)
  x[i+1] = x[i] + (v_x[i]*dt)
  y[i+1] = y[i] + (v_y[i]*dt)

fig, axs = plt.subplots(1, 1, figsize=(12, 8))
plt.plot(y, x, label ='Position')
plt.title("Simulation of Brownian motion") 
plt.show()

The result is this,

I want to figure out my mistake. Please help

Answer 1

Well, that's not really a programming question. These lines

for i in range(n-1):
  v_x[i + 1] = v_x[i] +  sqrtdt * np.random.normal(0.0,1.0)
  v_y[i + 1] = v_y[i] +  sqrtdt * np.random.normal(0.0,1.0)
  x[i+1] = x[i] + (v_x[i]*dt)
  y[i+1] = y[i] + (v_y[i]*dt)

are just simply not true, because it's a SDE.

The general form of the equation is dx = a(t, x)dt + b(t, x)dW, where a(t, x) is deterministic, b(t, x) is stochastic in nature (Wiener process). Making it numeric it becomes

x[n+1] = x[n] + dx = x[n] + a(t, x[n])dt + b(t, x[n]) sqrt(dt) ξ, where ξ is normally distributed with mean 0 and variance 1. The sqrt(dt) comes from the properties of the Wiener process.

Instead of using Euler method you should go for Euler-Maruyama. These are the right equations:

for i in range(n - 1):
    x[i + 1] = x[i] + b_x(t, x) * sqrtdt * np.random.normal(0.0, 1.0)
    y[i + 1] = y[i] + b_y(t, y) * sqrtdt * np.random.normal(0.0, 1.0)

and in your case b_x(t, x) = b_y(t, y) = 1

Answer 2

Try the package sdeint. I think, it will do exactly what you want to do.