Langevin Equation - Stochastic Differential Equation. What are the subtleties? 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.
 A: $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).
A: 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.
A: 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.
