Convergence of a numerical method for solving discontinuous Fokker-Planck equations

In studies of molecular motors, the stochastic motion is modeled using the Langevin equation. If we consider an ensemble of motors, the probability density is governed by the corresponding Fokker-Planck equation. Average quantities, such as average velocity, effective diffusion coefficient, and randomness parameter, can be calculated from the probability density. A numerical method was previously developed to solve Fokker-Planck equations [ H. Wang, C. Peskin, and T. Elston, J. Theoret. Biol., 221 ( 2003), pp. 491-511]. It preserves detailed balance, which ensures that if the system is forced to an equilibrium, the numerical solution will be the same as the Boltzmann distribution. Here we study the convergence of this numerical method when the potential has a finite number of discontinuities at half numerical grid points. We prove that this numerical method is stable and is consistent with the differential equation. Based on the consistency analysis, we propose a modified version of this numerical method to eliminate the first order error term caused by the discontinuity. We also show that in the presence of discontinuities, detailed balance is a necessary condition for converging to the correct solution. This explains why the central difference method converges to a wrong solution.