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.

  • 3
    $\begingroup$ Why would you do this numerically? It can be solved analytically. $\endgroup$ Commented Sep 1, 2015 at 8:50
  • $\begingroup$ What kind of probability/statistics background do you have? $\endgroup$
    – Kyle Kanos
    Commented Sep 1, 2015 at 10:12
  • $\begingroup$ @StevenMathey I'm trying to study the mean squared displacement of many particles $\endgroup$
    – Candy Man
    Commented Sep 2, 2015 at 2:01
  • $\begingroup$ @KyleKanos Not much - I took a course in Statistical Mechanics before, but thats really just about it. $\endgroup$
    – Candy Man
    Commented Sep 2, 2015 at 2:03

3 Answers 3


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


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.

  • $\begingroup$ A question, in the solution to the equation what is the reason for appearance of $e^{\frac{6\pi a \nu}{m}t}$ in the integrand? Wouldn`t it also be the solution if we simply integrated $F(t)$ since it is, I am guessing usually an integrable function of the variable $t$? $\endgroup$
    – Sina
    Commented Feb 20, 2016 at 14:42
  • $\begingroup$ @Sina Sorry, but I don't really understand your question. The first equation of my answer holds for an arbitrary function $F(t)$. Since we know almost nothing about $F(t)$, we can not perform the integration over $\tau_2$ or $\tau_1$ (since it is the boundary of the $\tau_2$-integration). Does this help? $\endgroup$ Commented Feb 20, 2016 at 18:31
  • $\begingroup$ Sorry I understood now. Though isnt $F(t)$ which contains the random force atleast supposed to be integrable? Or not? $\endgroup$
    – Sina
    Commented Feb 21, 2016 at 0:33
  • $\begingroup$ @Sina $F(t)$ is white noise in time. That means that $F(t)$ and $F(t+dt)$ are completely uncorrelated even if $dt$ is very small. $F(t)$ is indeed integrable but is highly irregular. The Riemann definition of integrals must be generalised to apply it here. You can read up on Ito calculus for more details. $\endgroup$ Commented Feb 21, 2016 at 8:45

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.