Langevin Equation - Stochastic Differential Equation. What are the subtleties?

Question

I am trying to find out the motion of a particle in 3D governed by the Langevin equation, numerically.

Anyway, the Langevin equation is given by

$$m \ddot{x} = -(6\pi a\nu) \dot{x} + F_b $$

where $F_b$ is due to random fluctuations.

From various sources I've read, $F_b$ is treated as stochastic term. According to wiki http://en.wikipedia.org/wiki/Langevin_dynamics it translates to

$$m \ddot{x} = -(6\pi a\nu) \dot{x} + \sqrt{2\gamma K_BTm}R(t)$$

where $\gamma$ is a friction term, and $R(t)$ is a delta-correlated stationary Gaussian process with zero-mean.

If my guess is correct, there isn't an explicit form of $R(t)$?

So I'm wondering if anyone can explain what $R(t)$ means and how I can go about trying to implement it into a simulation for example.

Answer 1

$R(t)$ is a function of time that represents complicated time-dependence of forces due to other molecules on the studied molecule.

Since only correlation function is assumed, there is no single unique function $R(t)$ assumed; although not all, many functions would be appropriate. You can generate many of them in computer using Cholesky decomposition of correlation matrix or discrete Fourier transform methods (faster).

Answer 2

The exact solution of your equation can be written as

$$x(t) = x(0) + \frac{m}{6 \pi a \nu} \dot{x}(0) - \frac{m}{6 \pi a \nu} \dot{x}(0) \, \text{e}^{-\frac{6\pi a\nu}{m}t} + \frac{1}{m}\int_0^t \text{d}\tau_1 \int_0^{\tau_1} \text{d}\tau_2 \text{e}^{\frac{6\pi a\nu}{m} \left(\tau_2-\tau_1\right)} F(\tau_2) \, . $$

$x(0)$ and $\dot{x}(0)$ are the initial conditions. Then whatever you want to compute, you can write it as an expression that depends on $F(t)$ and take its average over the fluctuations of $F(t)$.

For example

\begin{align*} \langle x(t) \rangle = & \langle x(0) \rangle + \frac{m}{6 \pi a \nu} \langle \dot{x}(0) \rangle- \frac{m}{6 \pi a \nu} \langle \dot{x}(0) \rangle \, \text{e}^{-\frac{6\pi a\nu}{m}t} \\ & + \frac{1}{m}\int_0^t \text{d}\tau_1 \int_0^{\tau_1} \text{d}\tau_2 \text{e}^{\frac{6\pi a\nu}{m} \left(\tau_2-\tau_1\right)} \langle F(\tau_2) \rangle \, .\end{align*}

If your particle starts at rest and at the origin,

$$ \langle x(0) \rangle = 0 \, \qquad \langle \dot{x}(0) \rangle = 0 \, ,$$

and if it experiences a constant force,

$$ \langle F(t) \rangle = F \, ,$$

then you find

\begin{align*}\langle x(t) \rangle & = \frac{F}{m}\int_0^t \text{d}\tau_1 \int_0^{\tau_1} \text{d}\tau_2 \text{e}^{\frac{6\pi a\nu}{m} \left(\tau_2-\tau_1\right)} \\ & = \frac{F}{m} \frac{m}{6 \pi a \mu} \left[ t + \frac{m}{6 \pi a \nu} \left( \text{e}^{-\frac{6 \pi a \nu}{m}t}-1\right) \right] \, .\end{align*}

You see that you can choose the statistics of $F(t)$ freely. Then if you know the moments of $F(t)$, you can compute the moments of $x(t)$. It's all about computing integrals. Typically one chooses Gaussian statistics with

$$ \langle F(t) \rangle = 0 \, , \qquad \langle F(t_1) F(t_2) \rangle = D \, \delta(t_1-t_2) \, .$$

If you insist on solving this problem numerically, you need to discretise time

$$ t \in \left[0,\infty\right[ \rightarrow t \in \left\{t_i\right\}_{i = 1,..,N} \, . $$

Then you can sample $F(t)$ according to your favourite probability distribution

$$ P\left[F(t_1),F(t_2),..,F(t_N)\right] \, .$$

For every sample you get a discretised function, $F(t_i)$ and you can switch to finite difference derivatives (for example) to solve your differential equation. Then you average at the end.

Answer 3

Solving this numerically is pretty much like Runge-Kutta methods except that $R(t)$ is not represented by a usual function but is usually a number generated pseudorandom generators supplied by your langage. You can use for instance Heun Midpoint Method. Assume you are at state $x_0,v_0$ at time $t_0$ and you want to find your state $x_1, v_1$ at time $t_1$. You first find a first order approximation: $$ \tilde{r} = r_0 + \Delta t v_0 $$ $$ \tilde{v} = v_0 - \frac{\Delta t}{m}( \xi v + \nabla U(r_0)) +R(t_0) $$ where $U$ is the potential energy of your system so $-\nabla U$ is the force due to internal interactions and $\xi$ is the friction coefficient. Now you use this approximation to find your real next step as

$$ {r}_1 = r_0 + \Delta t \frac{1}{2}(v_0+\tilde{v}) $$ $$ {v}_1 = v_0 - \frac{1}{2} \frac{\Delta t}{m}( \xi (v_0 +\tilde{v}) + (\nabla U(r_0)+\nabla U(\tilde{r}))) +R(t_0). $$ This is accurate up to second order in time I believe. There are first order ones as well such as Euler-Maruyana. As you can see it is not very different from usual deterministic integrators.