Numerical integration of Fokker-Planck equation allowing for negative drift?

Question

The Fokker-Planck equation (a.k.a Kolmogorov forward equation or Smoluchowski equation) describes the evolution of a probability density function and numerical integration of the FPE should conserve probability and ensure probabilites remain non-negative. The Chang-Cooper method is a popular finite-difference scheme which satisfies these constraints, however it requires that in the FPE of the form

$$\frac{\partial u}{\partial t} = \frac{1}{A(x)}\frac{\partial}{\partial x} \big[B(x,t)u + C(x,t)\frac{\partial u}{\partial x}\big]$$

the functions $A(x), B(x,t)$ and $C(x,t)$ are positive. For some problems, however, such as the Fokker-Planck equation describing Brownian motion of a particle inside a potential, $B(x,t)$ may not be a positive function.

Is there a method that exists for integrating the FPE that allows for a negative $B(x,t)$ term?

Answer 1

I am familiar with the equation in a different application, it is a so called convection-diffusion equation. There is no problem to solve it for negative values of B. The scheme you mention (Chang-Cooper method) is known in my field as the exponential upwind scheme and it is valid for all nonzero values of B. Due to a lack of time I put here screen shots of my paper with a slightly different notations:

Frolkovic, P. "Maximum principle and local mass balance for numerical solutions of transport equation coupled with variable density flow." Acta Mathematica Universitatis Comenianae 1.68 (1998): 137-157.

that you can find easily on internet. I do not know which paper you use to define the Chang-Cooper method, but I hope you can understand it using my explanation below.

Here is the screenshot with the stationary form of the equation where the unknown function is $c$. If you take a single representative constant values of $q$ (your $B$) and $d>0$ (your $C>0$), you obtain an exact solution in a form of exponential interpolation, see (40), where two values of numerical solution $c_i$ and $c_j$ in two grid points are interpolated (you can think in 1D that $j=i+1$):

An important role is played by so called Peclet number $P$, see above, that can be positive or negative.

If you use this interpolation you can derive the exponential upwind scheme, see the next screenshot. The idea can be explained by using a standard linear interpolation for the convection term and standard finite difference for the diffusion term where the diffusion coefficient is multiplied by a value of $K \ge 1$, it means that the effective numerical diffusion is enlarged, see (42) below (and you should recognize the formula (43) in your definition of the Chang-Cooper method)

Finally I add screenshots of an output from Mathematica where you see that the exponential interpolation is working for positive and negative values of $q$ (your $B$) for $P=-100,-10,-1,1,10,100$ and where you see a plot of different values of $K$ for different values of Peclet number $P$.