Consider the Langevin equation in the overdamped regime,
$$ 0 = -\gamma \dot{\mathbf{x}} -\nabla U(\mathbf{x}) +\boldsymbol{\eta}(t) \, $$
where $\boldsymbol{\eta}$ is the usual white-noise term, $U$ a potential for the force and $\gamma$ the damping coefficient (or a "damping matrix"). Thanks to the good reference provided in the accepted answer, I found how to derive the associated Fokker-Planck equation for this system.
Assuming that we have our Fokker-Planck equation for the particle distribution $P(\mathbf{x},t)$, I imagine (but I am not sure) that the average velocity of particles is given by
$$\langle \dot{\mathbf{x}}(t) \rangle = -\int d^Nx \, P(\mathbf{x},t) \gamma^{-1} \nabla U(\mathbf{x}). $$
Now here is my doubt: at $t=0$ we could choose a certain $P(\mathbf{x},0)$, find $P(\mathbf{x},t)$ with the Fokker Planck and calculate $\langle \dot{\mathbf{x}}(t) \rangle$ as above. Alternatively, we could sample $M$ different initial conditions $\mathbf{x}_i(0)$ from $P(\mathbf{x},0)$, evolve each $\mathbf{x}_i(t)$ for $i=1...M$ with the Langevin equation and obtain
$$\langle \dot{\mathbf{x}}(t) \rangle \approx M^{-1} \sum_i \dot{\mathbf{x}}_i(t) .$$
If this is correct, which method in general is more convenient from the numerical point of view? I see a big difference: simulating a single PDE (the Fokker-Planck) and performing an integral VS simulating an high number $M$ of ODEs (but performing a simple sum).
