Fokker-Planck equation

The **Fokker-Planck equation**, also known as **Kolmogorov forward equation**, is a hyperbolic partial differential equation, the evolution equation of the probability density of a certain class of stochastic differential equations.

$${\mathrm{dX}}_{t}=\phantom{\rule{thickmathspace}{0ex}}f({X}_{t},t)\mathrm{dt}+\phantom{\rule{thickmathspace}{0ex}}g({X}_{t},t)\phantom{\rule{thickmathspace}{0ex}}{\mathrm{dW}}_{t}$$

has an associated Fokker-Planck equation (as physicists say) aka Kolmogorov forward equation (as mathematicians say), that describes the time evolution of the probability distribution of the Langevin equation: From the physicist’s POV this probability distribution tells us what the probability is to find the particle at a certain time in a certain place.

(The **Smoluchowski Equation** is a special kind of Fokker-Planck equation.)

For a one dimensional Ito stochastic differential equation of the form

$${\mathrm{dX}}_{t}=\phantom{\rule{thickmathspace}{0ex}}f({X}_{t},t)\mathrm{dt}+\phantom{\rule{thickmathspace}{0ex}}g({X}_{t},t)\phantom{\rule{thickmathspace}{0ex}}{\mathrm{dW}}_{t}$$

the probability density is a solution of this one-dimensional Fokker-Planck equation:

$$\frac{\partial}{\partial t}P(x,t)=(-\frac{\partial}{\partial x}f(x,t)+\frac{1}{2}\frac{{\partial}^{2}}{\partial {x}^{2}}{g}^{2}(x,t))P(x,t)$$

The mathematical theorem making all of this precise is the Feynman-Kac formula.

The differential operator on the right side is sometimes called the **Fokker-Planck operator**.

While the Fokker-Planck equation describes the evolution of the probability distribution **forwards in time**, the **backward** Fokker-Planck equation, or **Kolmogorov backward equation**, describes the evolution of the probability distribution **backwards in time**. The equation is an evolution equation with the **formal adjoint** of the Fokker-Planck operator:

$$\frac{\partial}{\partial t}P(x,t)=(-f(x,t)\frac{\partial}{\partial x}+\frac{1}{2}{g}^{2}(x,t)\frac{{\partial}^{2}}{\partial {x}^{2}})P(x,t)$$

The Fokker-Planck equation describes the time evolution of a probability distribution, and therefore satisfies a conservation law that describes that no probability is lost, the probability of the trivial event is still 1. For the sake of simplicity we state it in one spatial dimension only, we have:

$$\frac{\partial}{\partial t}P(x,t)+\frac{\partial}{\partial x}J(x,t)=0$$

with the **probability current**

$$J(x,t):=f(x,t)\phantom{\rule{thickmathspace}{0ex}}P(x,t)-\frac{1}{2}\phantom{\rule{thickmathspace}{0ex}}\frac{\partial}{\partial x}g(x,t)\phantom{\rule{thickmathspace}{0ex}}P(x,t)$$

A **stationary solution** is a solution that does not depend on time and is therefore an eigenvector of the Fokker-Planck operator to the eigenvalue $0$. In one dimension without explicit time dependency of the Fokker-Planck operator, it can be calculated in closed form:

$${P}_{s}(x)=\frac{N}{{g}^{2}(x)}\mathrm{exp}(2{\int}_{0}^{x}\phantom{\rule{thickmathspace}{0ex}}\frac{f(\chi )}{{g}^{2}(\chi )}\phantom{\rule{thickmathspace}{0ex}}d\chi )$$

$N\in {\mathbb{R}}_{+}$ is a normalization constant to ensure that

$${\int}_{-\mathrm{\infty}}^{\mathrm{\infty}}{P}_{s}(x)\mathrm{dx}=1$$

This formula is obviously valid only if the function ${P}_{s}$ can be normalized in this way, of course.

A particle in one dimension in a potential $V(x)$ that is described by a Langevin equation of the form

$${\mathrm{dX}}_{t}=-V\prime (x)\mathrm{dt}+\sqrt{2D}\phantom{\rule{thickmathspace}{0ex}}{\mathrm{dW}}_{t}$$

has the following Fokker-Planck equation:

$$\frac{\partial}{\partial t}P(x,t)=(\frac{\partial}{\partial x}V\prime (x)+D\frac{{\partial}^{2}}{\partial {x}^{2}})P(x,t)$$

Therefore, the stationary solution is according to our general formula, this:

$${P}_{s}(x)=\frac{N}{D}\mathrm{exp}(-\frac{1}{D}V(x))$$

This is true for all potential functions that are differentiable and such that $\mathrm{exp}(-\frac{1}{D}V(x))\in {L}^{1}(\mathbb{R})$.

In our toy example for stochastic resonance, the stationary solution to the time independent system with the bistable potential

$$V(x):=\frac{1}{4}{x}^{4}-\frac{1}{2}{x}^{2}$$

is, according to the preceding paragraph, simply

$${P}_{s}(x)=\frac{N}{D}\mathrm{exp}(-\frac{1}{D}V(x))$$

The following picture shows both the potential and the stationary solution (not normalized though).

For applications it is often of interest, when a system or a particle leaves a certain area and crosses a threshold. This kind of problem is often called a First-hitting-time problem, and these can often be handled using techniques involving Fokker-Planck equations.

Here we will consider an example discussed here:

- Hopf bifurcation, Azimuth Project.

We have transformed the SDE system to polar coordinates as an example application of the Itô formula.

The noise matrix of the SDE system is

$$B=\left(\begin{array}{cc}\lambda \mathrm{cos}(\varphi )& \lambda \mathrm{sin}(\varphi )\\ -\frac{\lambda}{r}\mathrm{sin}(\varphi )& \frac{\lambda}{r}\mathrm{cos}(\varphi )\end{array}\right)$$

In the multidimensional case we have to calculate the diffusion matrix $D$ which is

$$D=B{B}^{T}=\left(\begin{array}{cc}{\lambda}^{2}& 0\\ 0& \frac{{\lambda}^{2}}{{r}^{2}}\end{array}\right)$$

The Fokker-Planck equation for the probability distribution $p(r,\varphi ,t)$ is therefore

$${\partial}_{t}p=-{\partial}_{r}(\beta \phantom{\rule{thickmathspace}{0ex}}r-{r}^{3}+\frac{{\lambda}^{2}}{2\phantom{\rule{thickmathspace}{0ex}}r})p+\frac{1}{2}{\lambda}^{2}{\partial}_{\mathrm{rr}}p+\frac{{\lambda}^{2}}{{r}^{2}}{\partial}_{\varphi \varphi}p$$

Note that the Fokker-Planck equation has no explicit dependency on the angle $\varphi $: We expected this, since the original system is rotationally invariant.

Some handwaving tells us that there should be a stationary equillibirium solution, and that this should be rotationally invariant, too. If we assume that $p$ actually depends on $r$ only, the Fokker-Planck equation simplifies to

$$0=-{\partial}_{r}(\beta \phantom{\rule{thickmathspace}{0ex}}r-{r}^{3}+\frac{{\lambda}^{2}}{2\phantom{\rule{thickmathspace}{0ex}}r})p+\frac{1}{2}{\lambda}^{2}{\partial}_{\mathrm{rr}}p$$

If we further impose the conditions that the function value and its first derivative should vanish in $0$, we get as solution

$$p(r)=Nr\mathrm{exp}(-\frac{1}{2{\lambda}^{2}}{r}^{4}+\frac{\beta}{{\lambda}^{2}}{r}^{2})$$

The normalization constant N has to be determined numerically.

We’ll consider two examples, $\beta =-0.25$ and $\lambda =0.1$ and $\lambda =0.3$. For the first case we get $N\approx 56.5$, for the second $N\approx 9.3$.

The following picture compares the analytic solution to a histogram that was obtained by running the SDE simulation used for This Weeks Finds 308, with taking a sample of $r$ every two time units:

(Green is $\lambda =0.3$ and red is $\lambda =0.1$.)

- Fokker-Planck equation, Wikipedia

The following book is a classic reference for the practicioner, Risken explains model building and approximate solution methods like linear response theory? using a minimum of mathematical background:

- Hannes Risken:
*The Fokker-Planck equation. Methods of solution and applications.*(ZMATH)

The mathematically inclined reader will find a nice exposition in the following classic textbook, chapter 4 “Brownian Motion and Partial Differential Equations”:

- Ioannis Karatzas and Steven E Shreve:
*Brownian motion and stochastic calculus.*(ZMATH)

category: mathematical methods