3
$\begingroup$

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.

$\endgroup$
4
  • 3
    $\begingroup$ Why would you do this numerically? It can be solved analytically. $\endgroup$ Sep 1, 2015 at 8:50
  • $\begingroup$ What kind of probability/statistics background do you have? $\endgroup$
    – Kyle Kanos
    Sep 1, 2015 at 10:12
  • $\begingroup$ @StevenMathey I'm trying to study the mean squared displacement of many particles $\endgroup$
    – Candy Man
    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
    Sep 2, 2015 at 2:03

3 Answers 3

3
$\begingroup$

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

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
4
  • $\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, 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$ 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
    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$ Feb 21, 2016 at 8:45
1
$\begingroup$

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.

$\endgroup$

Your Answer

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