Immersed Interface Method for Fokker-Planck Equation with Discontinuous Drift ()
1. Introduction
In recent years, piecewise-smooth stochastic systems (governed by piecewise-smooth stochastic differential equations) are usually used to describe biological and physical systems. Although for some simple piecewise-linear stochastic differential equations, analytical solutions of the transition probability distribution can be obtained [1] [2], it is difficult to attain analytical expressions for many other cases. Hence, we need to develop some effective numerical methods to deal with the difficulty in order to know more dynamical behaviors of the systems.
In this paper, we attempt to solve numerically a Fokker-Planck equation with discontinuous drift, which results from a so-called Brownian motion with pure dry friction [3]. This dry friction model can be described as the following piecewise linear Langevin equation
(1.1)
Here
denotes the sign of the velocity
, representing the dry friction force.
is the Gaussian white noise with zero mean and delta correlation
with
. The notation
stands for the average overall possible realizations of the noise, and
is the Dirac delta function. The transition probability distribution
of (1.1) satisfies the following Fokker-Planck equation [4] [5],
(1.2)
The corresponding initial condition is
if
for (1.1).
Since Equation (1.2) has a discontinuous drift
, we must deal with it carefully. The IIM is a sharp interface method which can accurately capture discontinuities in the solution and the flux. This method has been used for many problems, such as elliptic interface problems [6], parabolic interface problems [7], moving interface problems [8] and many other applications [9] [10] [11] (see [12] [13] for excellent reviews). To the best of our knowledge, there is no literature about the IIM for solving Fokker-Planck equations with discontinuous drift so far. Hence, our goal is to solve it.
The rest of this paper is organized as follows. In Section 2, we derive the IIM for the Fokker-Planck Equation (1.2). The numerical results are compared with the analytical solutions in Section 3. In addition, the accuracy of the scheme is also obtained. Finally, conclusions are made in Section 4.
2. The Scheme
We set
for convenience. At the discontinuous point
, we have the matching condition for the solution,
(2.1)
where
and
stand for the limiting values from the right- and left-hand sides of
. Integrating (1.2) across the discontinuity, we find
(2.2)
and then
(2.3)
by replacing
with
in (2.2).
It follows from (2.1) that
, that is
according to Equation (1.2). Then using the relations (2.1)-(2.3) we have
(2.4)
and
(2.5)
For the numerical scheme, we have first to truncate the computational domain
to a finite domain. Without loss of generality, let us assume the finite domain to be
, where
is a positive constant. Then we assume the probability vanishes at the boundary, i.e.,
. (2.6)
A uniform grid with step
is chosen here, where
is a positive constant.
Therefore, the grid points can be expressed as
,
with the discontinuous point being between
and
,
.
We hope to develop finite difference scheme of the form
,
, (2.7)
where
is the time-step size. This means that we need to determine the coefficients
and the correction term
so that
(2.8)
At a regular grid point
,
, the coefficients
in the explicit difference scheme (2.7) are obtained by the standard approximation as follows
(2.9)
and the correction term
At the irregular grid point
, we expand
,
and
in Taylor series at the discontinuous point
to obtain
(2.10)
(2.11)
(2.12)
For Equation (2.12), using (2.1), (2.2) and (2.4), we have
(2.13)
Furthermore, substituting (2.10), (2.11) and (2.13) into (2.8) we have
(2.14)
Then by arranging terms we obtain
(2.15)
Comparing both sides of (2.15), one obtains three equations for
,
and
as follows
(2.16)
(2.17)
(2.18)
and the correction term
Therefore, one can solve (2.16) - (2.18) to attain the coefficients of Equation (2.7) for
In a similar way, we can compute the coefficients at the irregular grid point
from the equations
(2.19)
(2.20)
(2.21)
and the correction term
3. Numerical results
For the Fokker-Planck Equation (1.2), using spectral decomposition method, one can get the transition probability distribution in closed analytic form [1] [14]:
(3.1)
where
(3.2)
is the transition probability distribution in non-dimensional units and
(3.3)
is the error function. In addition, when
the Fokker-Planck Equation (1.2) admits a steady stationary state
(3.4)
Let
,
and the computing interval be
. We choose the space-step
and the time-step
. For simplicity, we take the analytic distribution (3.1) at time
as the initial condition for computing. Figure 1 shows the comparison of numerical and analytical results of the probability distribution
at different times. It can be seen that the numerical solutions (points) coincide with the exact solutions (solid lines), indicating the effectiveness of the Scheme (2.7).
To see the accuracy of the scheme numerically, we consider the
and
errors between the numerical solutions and the exact solutions defined by
(3.5)
(3.6)
where
is the numerical solution and
is the exact solution. Then we calculate the order of accuracy. A small time-step
and
are chosen and the problem is recalculated from time
to
. As illustrated in Table 1, the scheme is approximated second order in the velocity direction.
Figure 1. Transition probability distribution
of Fokker-Planck Equation (1.2) with solid lines corresponding to the exact solutions, points to the numerical solutions, and dashed line to the stationary solution.
Table 1. Accuracy test in the velocity direction for
and
.
4. Conclusions
We have used the IIM to solve a Fokker-Planck equation with discontinuous drift in this paper. The numerical results show that the developed scheme is effective and has second order of accuracy. Moreover, the scheme can be readily extended to other dry friction models and the numerical results obtained are important references to see whether the dry friction effect exists in engineering applications.
Acknowledgements
This work was supported by the National Natural Science Foundation of China (Grant Nos. 11571366 and 11601517).