Fokker–Planck equation

In one spatial dimension x, for an Itô process driven by the standard Wiener process ${\displaystyle W_{t}}$ and described by the stochastic differential equation (SDE)

$${\displaystyle dX_{t}=\mu (X_{t},t)\,dt+\sigma (X_{t},t)\,dW_{t}}$$

with drift ${\displaystyle \mu (X_{t},t)}$ and diffusion coefficient ${\displaystyle D(X_{t},t)=\sigma ^{2}(X_{t},t)/2}$, the Fokker–Planck equation for the probability density ${\displaystyle p(x,t)}$ of the random variable ${\displaystyle X_{t}}$ is [9]

Link between the Itô SDE and the Fokker–Planck equation

In the following, use ${\displaystyle \sigma ={\sqrt {2D}}}$.

Define the infinitesimal generator ${\displaystyle {\mathcal {L}}}$ (the following can be found in Ref.[10]):

$${\displaystyle {\mathcal {L}}p(X_{t})=\lim _{\Delta t\to 0}{\frac {1}{\Delta t}}\left(\mathbb {E} {\big [}p(X_{t+\Delta t})\mid X_{t}=x{\big ]}-p(x)\right).}$$

The transition probability ${\displaystyle \mathbb {P} _{t,t'}(x\mid x')}$, the probability of going from ${\displaystyle (t',x')}$ to ${\displaystyle (t,x)}$, is introduced here; the expectation can be written as

$${\displaystyle \mathbb {E} (p(X_{t+\Delta t})\mid X_{t}=x)=\int p(y)\,\mathbb {P} _{t+\Delta t,t}(y\mid x)\,dy.}$$

Now we replace in the definition of ${\displaystyle {\mathcal {L}}}$, multiply by ${\displaystyle \mathbb {P} _{t,t'}(x\mid x')}$ and integrate over ${\displaystyle dx}$. The limit is taken on

$${\displaystyle \int p(y)\int \mathbb {P} _{t+\Delta t,t}(y\mid x)\,\mathbb {P} _{t,t'}(x\mid x')\,dx\,dy-\int p(x)\,\mathbb {P} _{t,t'}(x\mid x')\,dx.}$$

Note now that

$${\displaystyle \int \mathbb {P} _{t+\Delta t,t}(y\mid x)\,\mathbb {P} _{t,t'}(x\mid x')\,dx=\mathbb {P} _{t+\Delta t,t'}(y\mid x'),}$$

which is the Chapman–Kolmogorov theorem. Changing the dummy variable ${\displaystyle y}$ to ${\displaystyle x}$, one gets

$${\displaystyle {\begin{aligned}\int p(x)\lim _{\Delta t\to 0}{\frac {1}{\Delta t}}\left(\mathbb {P} _{t+\Delta t,t'}(x\mid x')-\mathbb {P} _{t,t'}(x\mid x')\right)\,dx,\end{aligned}}}$$

which is a time derivative. Finally we arrive to

$${\displaystyle \int [{\mathcal {L}}p(x)]\mathbb {P} _{t,t'}(x\mid x')\,dx=\int p(x)\,\partial _{t}\mathbb {P} _{t,t'}(x\mid x')\,dx.}$$

From here, the Kolmogorov backward equation can be deduced. If we instead use the adjoint operator of ${\displaystyle {\mathcal {L}}}$, ${\displaystyle {\mathcal {L}}^{\dagger }}$, defined such that

$${\displaystyle \int [{\mathcal {L}}p(x)]\mathbb {P} _{t,t'}(x\mid x')\,dx=\int p(x)[{\mathcal {L}}^{\dagger }\mathbb {P} _{t,t'}(x\mid x')]\,dx,}$$

then we arrive to the Kolmogorov forward equation, or Fokker–Planck equation, which, simplifying the notation ${\displaystyle p(x,t)=\mathbb {P} _{t,t'}(x\mid x')}$, in its differential form reads

$${\displaystyle {\mathcal {L}}^{\dagger }p(x,t)=\partial _{t}p(x,t).}$$

Remains the issue of defining explicitly ${\displaystyle {\mathcal {L}}}$. This can be done taking the expectation from the integral form of the Itô's lemma:

$${\displaystyle \mathbb {E} {\big (}p(X_{t}){\big )}=p(X_{0})+\mathbb {E} \left(\int _{0}^{t}\left(\partial _{t}+\mu \partial _{x}+{\frac {\sigma ^{2}}{2}}\partial _{x}^{2}\right)p(X_{t'})\,dt'\right).}$$

The part that depends on ${\displaystyle dW_{t}}$ vanished because of the martingale property.

Then, for a particle subject to an Itô equation, using

$${\displaystyle {\mathcal {L}}=\mu \partial _{x}+{\frac {\sigma ^{2}}{2}}\partial _{x}^{2},}$$

it can be easily calculated, using integration by parts, that

$${\displaystyle {\mathcal {L}}^{\dagger }=-\partial _{x}(\mu \cdot )+{\frac {1}{2}}\partial _{x}^{2}(\sigma ^{2}\cdot ),}$$

which bring us to the Fokker–Planck equation:

$${\displaystyle \partial _{t}p(x,t)=-\partial _{x}{\big (}\mu (x,t)\cdot p(x,t){\big )}+\partial _{x}^{2}\left({\frac {\sigma (x,t)^{2}}{2}}\,p(x,t)\right).}$$

While the Fokker–Planck equation is used with problems where the initial distribution is known, if the problem is to know the distribution at previous times, the Feynman–Kac formula can be used, which is a consequence of the Kolmogorov backward equation.

The stochastic process defined above in the Itô sense can be rewritten within the Stratonovich convention as a Stratonovich SDE:

$${\displaystyle dX_{t}=\left[\mu (X_{t},t)-{\frac {1}{2}}{\frac {\partial }{\partial X_{t}}}D(X_{t},t)\right]\,dt+{\sqrt {2D(X_{t},t)}}\circ dW_{t}.}$$

It includes an added noise-induced drift term due to diffusion gradient effects if the noise is state-dependent. This convention is more often used in physical applications. Indeed, it is well known that any solution to the Stratonovich SDE is a solution to the Itô SDE.

The zero-drift equation with constant diffusion can be considered as a model of classical Brownian motion:

$${\displaystyle {\frac {\partial }{\partial t}}p(x,t)=D_{0}{\frac {\partial ^{2}}{\partial x^{2}}}\left[p(x,t)\right].}$$

This model has discrete spectrum of solutions if the condition of fixed boundaries is added for ${\displaystyle \{0\leq x\leq L\}}$:

$${\displaystyle p(0,t)=p(L,t)=0,}$$

$${\displaystyle p(x,0)=p_{0}(x).}$$

It has been shown[11] that in this case an analytical spectrum of solutions allows deriving a local uncertainty relation for the coordinate-velocity phase volume:

$${\displaystyle \Delta x\,\Delta v\geq D_{0}.}$$

Here ${\displaystyle D_{0}}$ is a minimal value of a corresponding diffusion spectrum ${\displaystyle D_{j}}$, while ${\displaystyle \Delta x}$ and ${\displaystyle \Delta v}$ represent the uncertainty of coordinate–velocity definition.

More generally, if

$${\displaystyle d\mathbf {X} _{t}={\boldsymbol {\mu }}(\mathbf {X} _{t},t)\,dt+{\boldsymbol {\sigma }}(\mathbf {X} _{t},t)\,d\mathbf {W} _{t},}$$

where ${\displaystyle \mathbf {X} _{t}}$ and ${\displaystyle {\boldsymbol {\mu }}(\mathbf {X} _{t},t)}$ are N-dimensional vectors, ${\displaystyle {\boldsymbol {\sigma }}(\mathbf {X} _{t},t)}$ is an ${\displaystyle N\times M}$ matrix and ${\displaystyle \mathbf {W} _{t}}$ is an M-dimensional standard Wiener process, the probability density ${\displaystyle p(\mathbf {x} ,t)}$ for ${\displaystyle \mathbf {X} _{t}}$ satisfies the Fokker–Planck equation

with drift vector ${\displaystyle {\boldsymbol {\mu }}=(\mu _{1},\ldots ,\mu _{N})}$ and diffusion tensor ${\textstyle \mathbf {D} ={\frac {1}{2}}{\boldsymbol {\sigma \sigma }}^{\mathsf {T}}}$, i.e.

$${\displaystyle D_{ij}(\mathbf {x} ,t)={\frac {1}{2}}\sum _{k=1}^{M}\sigma _{ik}(\mathbf {x} ,t)\sigma _{jk}(\mathbf {x} ,t).}$$

If instead of an Itô SDE, a Stratonovich SDE is considered,

$${\displaystyle d\mathbf {X} _{t}={\boldsymbol {\mu }}(\mathbf {X} _{t},t)\,dt+{\boldsymbol {\sigma }}(\mathbf {X} _{t},t)\circ d\mathbf {W} _{t},}$$

the Fokker–Planck equation will read:[10]^: 129

$${\displaystyle {\frac {\partial p(\mathbf {x} ,t)}{\partial t}}=-\sum _{i=1}^{N}{\frac {\partial }{\partial x_{i}}}\left[\mu _{i}(\mathbf {x} ,t)\,p(\mathbf {x} ,t)\right]+{\frac {1}{2}}\sum _{k=1}^{M}\sum _{i=1}^{N}{\frac {\partial }{\partial x_{i}}}\left\{\sigma _{ik}(\mathbf {x} ,t)\sum _{j=1}^{N}{\frac {\partial }{\partial x_{j}}}\left[\sigma _{jk}(\mathbf {x} ,t)\,p(\mathbf {x} ,t)\right]\right\}}$$

In general, the Fokker–Planck equations are a special case to the general Kolmogorov forward equation

$${\displaystyle \partial _{t}\rho ={\mathcal {A}}^{*}\rho }$$

where the linear operator ${\displaystyle {\mathcal {A}}^{*}}$ is the Hermitian adjoint to the infinitesimal generator for the Markov process.[12]

Consider an overdamped Brownian particle under external force ${\displaystyle F(r)}$:[14]

$${\displaystyle m{\ddot {r}}=-\gamma {\dot {r}}+F(r)+\sigma \xi (t)}$$

where the ${\displaystyle m{\ddot {r}}}$ term is negligible (the meaning of "overdamped"). Thus, it is just ${\displaystyle \gamma dr=F(r)dt+\sigma dW_{t}}$. The Fokker–Planck equation for this particle is the Smoluchowski diffusion equation:

$${\displaystyle \partial _{t}P(r,t|r_{0},t_{0})=\nabla \cdot [D(\nabla -\beta F(r))P(r,t|r_{0},t_{0})]}$$

Where ${\displaystyle D}$ is the diffusion constant and ${\displaystyle \beta ={\frac {1}{k_{\text{B}}T}}}$. The importance of this equation is it allows for both the inclusion of the effect of temperature on the system of particles and a spatially dependent diffusion constant.

Derivation of the Smoluchowski Equation from the Fokker–Planck Equation

Starting with the Langevin Equation of a Brownian particle in external field ${\displaystyle F(r)}$, where ${\displaystyle \gamma }$ is the friction term, ${\displaystyle \xi }$ is a fluctuating force on the particle, and ${\displaystyle \sigma }$ is the amplitude of the fluctuation.

$${\displaystyle m{\ddot {r}}=-\gamma {\dot {r}}+F(r)+\sigma \xi (t)}$$

At equilibrium the frictional force is much greater than the inertial force, ${\displaystyle \left\vert \gamma {\dot {r}}\right\vert \gg \left\vert m{\ddot {r}}\right\vert }$. Therefore, the Langevin equation becomes,

$${\displaystyle \gamma {\dot {r}}=F(r)+\sigma \xi (t)}$$

Which generates the following Fokker–Planck equation,

$${\displaystyle \partial _{t}P(r,t|r_{0},t_{0})=\left(\nabla ^{2}{\frac {\sigma ^{2}}{2\gamma ^{2}}}-\nabla \cdot {\frac {F(r)}{\gamma }}\right)P(r,t|r_{0},t_{0})}$$

Rearranging the Fokker–Planck equation,

$${\displaystyle \partial _{t}P(r,t|r_{0},t_{0})=\nabla \cdot \left(\nabla D-{\frac {F(r)}{\gamma }}\right)P(r,t|r_{0},t_{0})}$$

Where ${\displaystyle D={\frac {\sigma ^{2}}{2\gamma ^{2}}}}$. Note, the diffusion coefficient may not necessarily be spatially independent if ${\displaystyle \sigma }$ or ${\displaystyle \gamma }$ are spatially dependent.

Next, the total number of particles in any particular volume is given by,

$${\displaystyle N_{V}(t|r_{0},t_{0})=\int _{V}drP(r,t|r_{0},t_{0})}$$

Therefore, the flux of particles can be determined by taking the time derivative of the number of particles in a given volume, plugging in the Fokker–Planck equation, and then applying Gauss's Theorem.

$${\displaystyle \partial _{t}N_{V}(t|r_{0},t_{0})=\int _{V}dV\nabla \cdot \left(\nabla D-{\frac {F(r)}{\gamma }}\right)P(r,t|r_{0},t_{0})=\int _{\partial V}d\mathbf {a} \cdot j(r,t|r_{0},t_{0})}$$

$${\displaystyle j(r,t|r_{0},t_{0})=\left(\nabla D-{\frac {F(r)}{\gamma }}\right)P(r,t|r_{0},t_{0})}$$

In equilibrium, it is assumed that the flux goes to zero. Therefore, Boltzmann statistics can be applied for the probability of a particles location at equilibrium, where ${\displaystyle F(r)=-\nabla U(r)}$ is a conservative force and the probability of a particle being in a state ${\displaystyle r}$ is given as ${\displaystyle P(r,t|r_{0},t_{0})={\frac {e^{-\beta U(r)}}{Z}}}$.

$${\displaystyle j(r,t|r_{0},t_{0})=\left(\nabla D-{\frac {F(r)}{\gamma }}\right){\frac {e^{-\beta U(r)}}{Z}}=0}$$

$${\displaystyle \Rightarrow \nabla D=F(r)\left({\frac {1}{\gamma }}-D\beta \right)}$$

This relation is a realization of the fluctuation–dissipation theorem. Now applying ${\displaystyle \nabla \cdot \nabla }$ to ${\displaystyle DP(r,t|r_{0},t_{0})}$ and using the Fluctuation-dissipation theorem,

$${\displaystyle {\begin{aligned}\nabla \cdot \nabla DP(r,t|r_{0},t_{0})&=\nabla \cdot D\nabla P(r,t|r_{0},t_{0})+\nabla \cdot P(r,t|r_{0},t_{0})\nabla D\\&=\nabla \cdot D\nabla P(r,t|r_{0},t_{0})+\nabla \cdot P(r,t|r_{0},t_{0}){\frac {F(r)}{\gamma }}-\nabla \cdot P(r,t|r_{0},t_{0})D\beta F(r)\end{aligned}}}$$

Rearranging,

$${\displaystyle \Rightarrow \nabla \cdot \left(\nabla D-{\frac {F(r)}{\gamma }}\right)P(r,t|r_{0},t_{0})=\nabla \cdot D(\nabla -\beta F(r))P(r,t|r_{0},t_{0})}$$

Therefore, the Fokker–Planck equation becomes the Smoluchowski equation,

$${\displaystyle \partial _{t}P(r,t|r_{0},t_{0})=\nabla \cdot D(\nabla -\beta F(r))P(r,t|r_{0},t_{0})}$$

for an arbitrary force ${\displaystyle F(r)}$.

Brownian motion follows the Langevin equation, which can be solved for many different stochastic forcings with results being averaged (canonical ensemble in molecular dynamics). However, instead of this computationally intensive approach, one can use the Fokker–Planck equation and consider the probability ${\displaystyle p(\mathbf {v} ,t)\,d\mathbf {v} }$ of the particle having a velocity in the interval ${\displaystyle (\mathbf {v} ,\mathbf {v} +d\mathbf {v} )}$ when it starts its motion with ${\displaystyle \mathbf {v} _{0}}$ at time 0.

Brownian dynamics simulation for particles in 1-D linear potential compared with the solution of the Fokker–Planck equation

Theory

Starting with a linear potential of the form ${\displaystyle U(x)=cx}$ the corresponding Smoluchowski equation becomes,

$${\displaystyle \partial _{t}P(x,t|x_{0},t_{0})=\partial _{x}D(\partial _{x}+\beta c)P(x,t|x_{0},t_{0})}$$

Where the diffusion constant, ${\displaystyle D}$, is constant over space and time. The boundary conditions are such that the probability vanishes at ${\displaystyle x\rightarrow \pm \infty }$ with an initial condition of the ensemble of particles starting in the same place, ${\displaystyle P(x,t|x_{0},t_{0})=\delta (x-x_{0})}$.

Defining ${\displaystyle \tau =Dt}$ and ${\displaystyle b=\beta c}$ and applying the coordinate transformation,

$${\displaystyle y=x+\tau b,\ \ \ y_{0}=x_{0}+\tau _{0}b}$$

With ${\displaystyle P(x,t,|x_{0},t_{0})=q(y,\tau |y_{0},\tau _{0})}$ the Smoluchowki equation becomes,

$${\displaystyle \partial _{\tau }q(y,\tau |y_{0},\tau _{0})=\partial _{y}^{2}q(y,\tau |y_{0},\tau _{0})}$$

Which is the free diffusion equation with solution,

$${\displaystyle q(y,\tau |y_{0},\tau _{0})={\frac {1}{\sqrt {4\pi (\tau -\tau _{0})}}}e^{-{\frac {(y-y_{0})^{2}}{4(\tau -\tau _{0})}}}}$$

And after transforming back to the original coordinates,

$${\displaystyle P(x,t|x_{0},t_{0})={\frac {1}{\sqrt {4\pi D(t-t_{0})}}}\exp {\left[{-{\frac {(x-x_{0}+D\beta c(t-t_{0}))^{2}}{4D(t-t_{0})}}}\right]}}$$

Simulation

The simulation on the right was completed using a Brownian dynamics simulation.[16][17] Starting with a Langevin equation for the system,

$${\displaystyle m{\ddot {x}}=-\gamma {\dot {x}}-c+\sigma \xi (t)}$$

where ${\displaystyle \gamma }$ is the friction term, ${\displaystyle \xi }$ is a fluctuating force on the particle, and ${\displaystyle \sigma }$ is the amplitude of the fluctuation. At equilibrium the frictional force is much greater than the inertial force, ${\displaystyle \left|\gamma {\dot {x}}\right|\gg \left|m{\ddot {x}}\right|}$. Therefore, the Langevin equation becomes,

$${\displaystyle \gamma {\dot {x}}=-c+\sigma \xi (t)}$$

For the Brownian dynamic simulation the fluctuation force ${\displaystyle \xi (t)}$ is assumed to be Gaussian with the amplitude being dependent of the temperature of the system ${\textstyle \sigma ={\sqrt {2\gamma k_{\text{B}}T}}}$. Rewriting the Langevin equation,

$${\displaystyle {\frac {dx}{dt}}=-D\beta c+{\sqrt {2D}}\xi (t)}$$

where ${\textstyle D={\frac {k_{\text{B}}T}{\gamma }}}$ is the Einstein relation. The integration of this equation was done using the Euler–Maruyama method to numerically approximate the path of this Brownian particle.

Being a partial differential equation, the Fokker–Planck equation can be solved analytically only in special cases. A formal analogy of the Fokker–Planck equation with the Schrödinger equation allows the use of advanced operator techniques known from quantum mechanics for its solution in a number of cases. Furthermore, in the case of overdamped dynamics when the Fokker–Planck equation contains second partial derivatives with respect to all spatial variables, the equation can be written in the form of a master equation that can easily be solved numerically.[18] In many applications, one is only interested in the steady-state probability distribution ${\displaystyle p_{0}(x)}$, which can be found from ${\textstyle {\frac {\partial p(x,t)}{\partial t}}=0}$. The computation of mean first passage times and splitting probabilities can be reduced to the solution of an ordinary differential equation which is intimately related to the Fokker–Planck equation.

In mathematical finance for volatility smile modeling of options via local volatility, one has the problem of deriving a diffusion coefficient ${\displaystyle {\sigma }(\mathbf {X} _{t},t)}$ consistent with a probability density obtained from market option quotes. The problem is therefore an inversion of the Fokker–Planck equation: Given the density f(x,t) of the option underlying X deduced from the option market, one aims at finding the local volatility ${\displaystyle {\sigma }(\mathbf {X} _{t},t)}$ consistent with f. This is an inverse problem that has been solved in general by Dupire (1994, 1997) with a non-parametric solution.[19][20] Brigo and Mercurio (2002, 2003) propose a solution in parametric form via a particular local volatility ${\displaystyle {\sigma }(\mathbf {X} _{t},t)}$ consistent with a solution of the Fokker–Planck equation given by a mixture model.[21][22] More information is available also in Fengler (2008),[23] Gatheral (2008),[24] and Musiela and Rutkowski (2008).[25]

Every Fokker–Planck equation is equivalent to a path integral. The path integral formulation is an excellent starting point for the application of field theory methods.[26] This is used, for instance, in critical dynamics.

A derivation of the path integral is possible in a similar way as in quantum mechanics. The derivation for a Fokker–Planck equation with one variable ${\displaystyle x}$ is as follows. Start by inserting a delta function and then integrating by parts:

$${\displaystyle {\begin{aligned}{\frac {\partial }{\partial t}}p{\left(x',t\right)}&=-{\frac {\partial }{\partial x'}}\left[D_{1}(x',t)p(x',t)\right]+{\frac {\partial ^{2}}{\partial {x'}^{2}}}\left[D_{2}(x',t)p(x',t)\right]\\[5pt]&=\int _{-\infty }^{\infty }\mathrm {d} x\left(\left[D_{1}\left(x,t\right){\frac {\partial }{\partial x}}+D_{2}\left(x,t\right){\frac {\partial ^{2}}{\partial x^{2}}}\right]\delta \left(x'-x\right)\right)p\!\left(x,t\right).\end{aligned}}}$$

The ${\displaystyle x}$-derivatives here only act on the ${\displaystyle \delta }$-function, not on ${\displaystyle p(x,t)}$. Integrate over a time interval ${\displaystyle \varepsilon }$,

$${\displaystyle p(x',t+\varepsilon )=\int _{-\infty }^{\infty }\,\mathrm {d} x\left(\left(1+\varepsilon \left[D_{1}(x,t){\frac {\partial }{\partial x}}+D_{2}(x,t){\frac {\partial ^{2}}{\partial x^{2}}}\right]\right)\delta (x'-x)\right)p(x,t)+O(\varepsilon ^{2}).}$$

Insert the Fourier integral

$${\displaystyle \delta {\left(x'-x\right)}=\int _{-i\infty }^{i\infty }{\frac {\mathrm {d} {\tilde {x}}}{2\pi i}}e^{{\tilde {x}}{\left(x-x'\right)}}}$$

for the ${\displaystyle \delta }$-function,

$${\displaystyle {\begin{aligned}p(x',t+\varepsilon )&=\int _{-\infty }^{\infty }\mathrm {d} x\int _{-i\infty }^{i\infty }{\frac {\mathrm {d} {\tilde {x}}}{2\pi i}}\left(1+\varepsilon \left[{\tilde {x}}D_{1}(x,t)+{\tilde {x}}^{2}D_{2}(x,t)\right]\right)e^{{\tilde {x}}(x-x')}p(x,t)+O(\varepsilon ^{2})\\[5pt]&=\int _{-\infty }^{\infty }\mathrm {d} x\int _{-i\infty }^{i\infty }{\frac {\mathrm {d} {\tilde {x}}}{2\pi i}}\exp \left(\varepsilon \left[-{\tilde {x}}{\frac {(x'-x)}{\varepsilon }}+{\tilde {x}}D_{1}(x,t)+{\tilde {x}}^{2}D_{2}(x,t)\right]\right)p(x,t)+O(\varepsilon ^{2}).\end{aligned}}}$$

This equation expresses ${\displaystyle p(x',t+\varepsilon )}$ as functional of ${\displaystyle p(x,t)}$. Iterating ${\displaystyle (t'-t)/\varepsilon }$ times and performing the limit ${\displaystyle \varepsilon \rightarrow 0}$ gives a path integral with action

$${\displaystyle S=\int \mathrm {d} t\left[{\tilde {x}}D_{1}(x,t)+{\tilde {x}}^{2}D_{2}(x,t)-{\tilde {x}}{\frac {\partial x}{\partial t}}\right].}$$

The variables ${\displaystyle {\tilde {x}}}$ conjugate to ${\displaystyle x}$ are called "response variables".[27]

Although formally equivalent, different problems may be solved more easily in the Fokker–Planck equation or the path integral formulation. The equilibrium distribution for instance may be obtained more directly from the Fokker–Planck equation.

• Frank, Till Daniel (2005). Nonlinear Fokker–Planck Equations: Fundamentals and Applications. Springer Series in Synergetics. Springer. ISBN 3-540-21264-7.
• Gardiner, Crispin (2009). Stochastic Methods (4th ed.). Springer. ISBN 978-3-540-70712-7.
• Pavliotis, Grigorios A. (2014). Stochastic Processes and Applications: Diffusion Processes, the Fokker–Planck and Langevin Equations. Springer Texts in Applied Mathematics. Springer. ISBN 978-1-4939-1322-0.
• Risken, Hannes (1996). The Fokker–Planck Equation: Methods of Solutions and Applications. Springer Series in Synergetics (2nd ed.). Springer. ISBN 3-540-61530-X.