# Exercise session on the Langevin equation¶

In [1]:
from __future__ import division, print_function
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import math
import random
plt.rcParams['font.size'] = 16


## Impatient¶

https://try.jupyter.org/

## Langevin equation for the velocity¶

$$\dot v = - \gamma v + \xi$$
In [2]:
# Langevin equation for the velocity

v = 0
dt = 0.01
gamma = 0.1
v_factor = math.sqrt(2*dt)
F_data = []
v_data = []
for t in range(100000):
F = random.gauss(0,1)
v = v*(1-gamma*dt) + v_factor*F
F_data.append(F)
v_data.append(v)
F_data = np.array(F_data)
v_data = np.array(v_data)

In [3]:
plt.plot(v_data)

Out[3]:
[<matplotlib.lines.Line2D at 0x7ff1d061b1d0>]
In [4]:
tau = np.linspace(-len(v_data)*dt, len(v_data)*dt, 2*len(v_data)-1)
plt.plot(tau, np.correlate(v_data, v_data, mode='full')/len(v_data))
plt.plot(tau, np.exp(-gamma*np.abs(tau))/gamma)
plt.xlim(-20, 20)

Out[4]:
(-20, 20)
In [5]:
plt.plot(tau, np.correlate(v_data, F_data, mode='full')/len(v_data))
plt.xlim(-20, 20)

Out[5]:
(-20, 20)

## Langevin equation with colored noise¶

The noise is defined by its autocorrelation $$\langle F(\tau) F \rangle = 2 D m^2 \omega_c e^{-\omega_c |\tau|}$$ and the Langevin equation is $$\dot v = - \gamma v + F$$

In [6]:
# Langevin equation with colored noise

v = 0
F = 0
dt = 0.01 ; N_steps = 10000
gamma = 2
wc = 1
D = 1
v_factor = math.sqrt(2*D*wc)*dt
F_factor = math.sqrt(2*dt)
F_data = []
v_data = []
for t in range(N_steps):
F = F*(1-wc*dt) + F_factor*random.gauss(0,1)
v = v*(1-gamma*dt) + v_factor*F
F_data.append(F)
v_data.append(v)
time = np.arange(N_steps)*dt
F_data = np.array(F_data)
v_data = np.array(v_data)

In [7]:
plt.plot(time, F_data)
plt.plot(time, v_data)

Out[7]:
[<matplotlib.lines.Line2D at 0x7ff1d03762d0>]
In [8]:
n = len(F_data)
C_F = np.correlate(F_data, F_data, mode='full')/n
n_points = 101
tau = np.linspace(-n_points*dt, n_points*dt, 2*n_points)
plt.plot(tau, C_F[n-n_points:n+n_points])
plt.plot(tau, np.exp(-wc*np.abs(tau))/wc)

Out[8]:
[<matplotlib.lines.Line2D at 0x7ff1d061bc10>]
In [9]:
def vv(tau):
pref = 2*D*wc
t1 = 1/(2*gamma)*(1/(gamma+wc)-1/(gamma-wc))*np.exp(-gamma*np.abs(tau))
t2 = 1/(gamma-wc)/(gamma+wc)*np.exp(-wc*np.abs(tau))
return pref*(t1+t2)

In [10]:
n = len(v_data)
C_v = np.correlate(v_data, v_data, mode='full')/n
n_points = 401
tau = np.linspace(-n_points*dt, n_points*dt, 2*n_points)
plt.plot(tau, C_v[n-n_points:n+n_points])
plt.plot(tau, vv(np.abs(tau)))

Out[10]:
[<matplotlib.lines.Line2D at 0x7ff1cf961c90>]

## On white noise¶

To generate white noise, one use a random number generator with the normal (i.e. gaussian distributed) distribution. Below, I show its autocorrelation that should resemble a Dirac, have zero average and standard deviation 1.

In [11]:
data = np.random.normal(size=10000)

print('Average', np.mean(data))
print('Standard deviation', np.std(data))
plt.plot(np.correlate(data, data, mode='full') / len(data))
plt.title(r'$\langle \xi(\tau) \xi\rangle$')
plt.xlabel(r'$\tau$');

Average 0.00599770905559
Standard deviation 0.999392388272

In [ ]: