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$ Sep 1 '15 at 8:50
  • $\begingroup$ What kind of probability/statistics background do you have? $\endgroup$
    – Kyle Kanos
    Sep 1 '15 at 10:12
  • $\begingroup$ @StevenMathey I'm trying to study the mean squared displacement of many particles $\endgroup$
    – Candy Man
    Sep 2 '15 at 2:01
  • $\begingroup$ @KyleKanos Not much - I took a course in Statistical Mechanics before, but thats really just about it. $\endgroup$
    – Candy Man
    Sep 2 '15 at 2:03

$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
    Feb 20 '16 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$ Feb 20 '16 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
    Feb 21 '16 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$ Feb 21 '16 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, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.