Immersed Interface Method for Fokker-Planck Equation with Discontinuous Drift

Show more

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

${v}^{\prime}(t)=-\sigma \left(v(t)\right)+\xi (t).$ (1.1)

Here $\sigma \left(v(t)\right)$ denotes the sign of the velocity $v(t)$ , representing the dry friction force. $\xi (t)$ is the Gaussian white noise with zero mean and delta correlation $\langle \xi (t)\xi ({t}^{\prime})\rangle =2D\delta (t-{t}^{\prime})$ with $D>0$ . The notation $\langle \xb7\rangle $ stands for the average overall possible realizations of the noise, and $\delta $ is the Dirac delta function. The transition probability distribution $p(v,t|{v}_{0},0)$ of (1.1) satisfies the following Fokker-Planck equation [4] [5],

$\frac{\partial p(v,t|{v}_{0},0)}{\partial t}=\frac{\partial \left(\sigma (v)p(v,t|{v}_{0},0)\right)}{\partial v}+D\frac{{\partial}^{2}p(v,t|{v}_{0},0)}{\partial {v}^{2}}.$ (1.2)

The corresponding initial condition is $p(v,0|{v}_{0},0)=\delta (v-{v}_{0})$ if $v(0)={v}_{0}$ for (1.1).

Since Equation (1.2) has a discontinuous drift $\sigma (v)$ , 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 $D=1$ for convenience. At the discontinuous point $v=0$ , we have the matching condition for the solution,

${p}^{+}={p}^{-}$ (2.1)

where $+$ and $-$ stand for the limiting values from the right- and left-hand sides of $v=0$ . Integrating (1.2) across the discontinuity, we find

${p}_{v}^{+}={p}_{v}^{-}-2{p}^{-}$ (2.2)

and then

${p}_{v}^{-}={p}_{v}^{+}+2{p}^{+}$ (2.3)

by replacing ${p}^{-}$ with ${p}^{+}$ in (2.2).

It follows from (2.1) that ${p}_{t}^{+}={p}_{t}^{-}$ , that is

${\left(\sigma (v)p\right)}_{v}^{+}+{p}_{vv}^{+}={\left(\sigma (v)p\right)}_{v}^{-}+{p}_{vv}^{-}$

according to Equation (1.2). Then using the relations (2.1)-(2.3) we have

${p}_{vv}^{+}={p}_{vv}^{-}-2{p}_{v}^{-}+2{p}^{-}$ (2.4)

and

${p}_{vv}^{-}={p}_{vv}^{+}+2{p}_{v}^{+}+2{p}^{+}.$ (2.5)

For the numerical scheme, we have first to truncate the computational domain $(-\infty ,+\infty )$ to a finite domain. Without loss of generality, let us assume the finite domain to be $[-{L}_{v},{L}_{v}]$ , where ${L}_{v}$ is a positive constant. Then we assume the probability vanishes at the boundary, i.e.,

$p(\pm {L}_{v},t|{v}_{0},0)=0$ . (2.6)

A uniform grid with step $h=2{L}_{v}/M$ is chosen here, where $M$ is a positive constant.

Therefore, the grid points can be expressed as ${v}_{i}=-{L}_{v}+(i-1)h$ , $i=1,2,\cdot \cdot \cdot ,M+1$ with the discontinuous point being between ${v}_{j}$ and ${v}_{j+1}$ , ${v}_{j}\le 0<{v}_{j+1}$ .

We hope to develop finite difference scheme of the form

$\frac{{p}_{i}^{n+1}-{p}_{i}^{n}}{\tau}={\gamma}_{i,1}{p}_{i-1}^{n}+{\gamma}_{i,2}{p}_{i}^{n}+{\gamma}_{i,3}{p}_{i+1}^{n}+{C}_{i}^{n}$ , $i=2,3,\cdot \cdot \cdot ,M$ , (2.7)

where $\tau $ is the time-step size. This means that we need to determine the coefficients $\gamma $ and the correction term ${C}_{i}^{n}$ so that

${\gamma}_{i,1}{p}_{i-1}^{n}+{\gamma}_{i,2}{p}_{i}^{n}+{\gamma}_{i,3}{p}_{i+1}^{n}+{C}_{i}^{n}\approx {\left\{\frac{\partial \left(\sigma (v)p\right)}{\partial v}+\frac{{\partial}^{2}p}{\partial {v}^{2}}\right\}}_{v={v}_{i}}.$ (2.8)

At a regular grid point ${v}_{i}$ , $i\ne j,\text{\hspace{0.17em}}j+1$ , the coefficients $\gamma $ in the explicit difference scheme (2.7) are obtained by the standard approximation as follows

${\gamma}_{i,1}=\frac{1}{{h}^{2}}-\frac{\sigma ({v}_{i-1/2})}{2h},$ ${\gamma}_{i,2}=\frac{\sigma ({v}_{i+1/2})-\sigma ({v}_{i-1/2})}{2h}-\frac{2}{{h}^{2}},$ ${\gamma}_{i,3}=\frac{1}{{h}^{2}}+\frac{\sigma ({v}_{i+1/2})}{2h},$ (2.9)

and the correction term ${C}_{i}^{n}=0.$

At the irregular grid point ${v}_{j}$ , we expand ${p}_{j-1}^{n}$ , ${p}_{j}^{n}$ and ${p}_{j+1}^{n}$ in Taylor series at the discontinuous point $v=0$ to obtain

${p}_{j-1}^{n}={p}^{-}+{v}_{j-1}{p}_{v}^{-}+\frac{1}{2}{v}_{j-1}^{2}{p}_{vv}^{-}+O({h}^{3}),$ (2.10)

${p}_{j}^{n}={p}^{-}+{v}_{j}{p}_{v}^{-}+\frac{1}{2}{v}_{j}^{2}{p}_{vv}^{-}+O({h}^{3}),$ (2.11)

${p}_{j+1}^{n}={p}^{+}+{v}_{j+1}{p}_{v}^{+}+\frac{1}{2}{v}_{j+1}^{2}{p}_{vv}^{+}+O({h}^{3}).$ (2.12)

For Equation (2.12), using (2.1), (2.2) and (2.4), we have

${p}_{j+1}^{n}=(1-2{v}_{j+1}+{v}_{j+1}^{2}){p}^{-}+({v}_{j+1}-{v}_{j+1}^{2}){p}_{v}^{-}+\frac{1}{2}{v}_{j+1}^{2}{p}_{vv}^{-}+O({h}^{3}).$ (2.13)

Furthermore, substituting (2.10), (2.11) and (2.13) into (2.8) we have

$\begin{array}{l}{\gamma}_{j,1}\left({p}^{-}+{v}_{j-1}{p}_{v}^{-}+\frac{1}{2}{v}_{j-1}^{2}{p}_{vv}^{-}\right)+{\gamma}_{j,2}\left({p}^{-}+{v}_{j}{p}_{v}^{-}+\frac{1}{2}{v}_{j}^{2}{p}_{vv}^{-}\right)\\ +{\gamma}_{j,3}\left((1-2{v}_{j+1}+{v}_{j+1}^{2}){p}^{-}+({v}_{j+1}-{v}_{j+1}^{2}){p}_{v}^{-}+\frac{1}{2}{v}_{j+1}^{2}{p}_{vv}^{-}\right)+{C}_{j}^{n}+O(h)\\ ={\left\{\frac{\partial \left(\sigma (v)p\right)}{\partial v}+\frac{{\partial}^{2}p}{\partial {v}^{2}}\right\}}_{v={0}^{-}}+O(h).\end{array}$ (2.14)

Then by arranging terms we obtain

$\begin{array}{l}{p}^{-}\left({\gamma}_{j,1}+{\gamma}_{j,2}+(1-2{v}_{j+1}+{v}_{j+1}^{2}){\gamma}_{j,3}\right)+{p}_{v}^{-}\left({v}_{j-1}{\gamma}_{j,1}+{v}_{j}{\gamma}_{j,2}+({v}_{j+1}-{v}_{j+1}^{2}){\gamma}_{j,3}\right)\\ +{p}_{vv}^{-}\left(\frac{1}{2}{v}_{j-1}^{2}{\gamma}_{j,1}+\frac{1}{2}{v}_{j}^{2}{\gamma}_{j,2}+\frac{1}{2}{v}_{j+1}^{2}{\gamma}_{j,3}\right)+{C}_{j}^{n}+O(h)\\ =-{p}_{v}^{-}+{p}_{vv}^{-}+O(h).\end{array}$ (2.15)

Comparing both sides of (2.15), one obtains three equations for ${\gamma}_{j,1}$ , ${\gamma}_{j,2}$ and ${\gamma}_{j,3}$ as follows

${\gamma}_{j,1}+{\gamma}_{j,2}+(1-2{v}_{j+1}+{v}_{j+1}^{2}){\gamma}_{j,3}=0,$ (2.16)

${v}_{j-1}{\gamma}_{j,1}+{v}_{j}{\gamma}_{j,2}+({v}_{j+1}-{v}_{j+1}^{2}){\gamma}_{j,3}=-1,$ (2.17)

${v}_{j-1}^{2}{\gamma}_{j,1}+{v}_{j}^{2}{\gamma}_{j,2}+{v}_{j+1}^{2}{\gamma}_{j,3}=2,$ (2.18)

and the correction term ${C}_{j}^{n}=0.$

Therefore, one can solve (2.16) - (2.18) to attain the coefficients of Equation (2.7) for $i=j.$

In a similar way, we can compute the coefficients at the irregular grid point ${v}_{j+1}$ from the equations

$(1+2{v}_{j}+{v}_{j}^{2}){\gamma}_{j+1,1}+{\gamma}_{j+1,2}+{\gamma}_{j+1,3}=0,$ (2.19)

$({v}_{j}+{v}_{j}^{2}){\gamma}_{j+1,1}+{v}_{j+1}{\gamma}_{j+1,2}+{v}_{j+2}{\gamma}_{j+1,3}=1,$ (2.20)

${v}_{j}^{2}{\gamma}_{j+1,1}+{v}_{j+1}^{2}{\gamma}_{j+1,2}+{v}_{j+2}^{2}{\gamma}_{j+1,3}=2,$ (2.21)

and the correction term ${C}_{j+1}^{n}=0.$

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]:

$p(v,t|{v}_{0},0)=\frac{1}{D}\stackrel{^}{p}\left(\frac{1}{D}v,\frac{1}{D}t|\frac{1}{D}{v}_{0},0\right),$ (3.1)

where

$\stackrel{^}{p}(x,\tau |{x}^{\prime},0)=\frac{{e}^{-\tau /4}}{2\sqrt{\pi \tau}}{e}^{-(|x|-|{x}^{\prime}|)/2}{e}^{-{(x-{x}^{\prime})}^{2}/(4\tau )}+\frac{{e}^{-\left|x\right|}}{4}\left[1+\mathrm{erf}\left(\frac{\tau -\left(\left|x\right|+\left|{x}^{\prime}\right|\right)}{2\sqrt{\tau}}\right)\right]$ (3.2)

is the transition probability distribution in non-dimensional units and

$\mathrm{erf}(z)=\frac{2}{\sqrt{\pi}}{\displaystyle {\int}_{0}^{z}{e}^{-{t}^{2}}dt}$ (3.3)

is the error function. In addition, when $t\to \infty $ the Fokker-Planck Equation (1.2) admits a steady stationary state

$p(v)=\frac{1}{2D}{e}^{-\left|v\right|/D}.$ (3.4)

Let $D=1$ , ${v}_{0}=2$ and the computing interval be $[-10,10]$ . We choose the space-step $h=0.01$ and the time-step $\tau =5\times {10}^{-5}$ . For simplicity, we take the analytic distribution (3.1) at time $t=0.01$ as the initial condition for computing. Figure 1 shows the comparison of numerical and analytical results of the probability distribution $p(v,t|{v}_{0},0)$ 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 ${L}_{2}$ and ${L}_{\infty}$ errors between the numerical solutions and the exact solutions defined by

${\Vert Error\Vert}_{{L}_{2}}=\sqrt{h{\displaystyle \underset{i=1}{\overset{M}{\sum}}{({p}_{i}-{\stackrel{\xaf}{p}}_{i})}^{2}}},$ (3.5)

${\Vert Error\Vert}_{{L}_{\infty}}=\underset{1\le i\le M}{\mathrm{max}}|{p}_{i}-{\stackrel{\xaf}{p}}_{i}|,$ (3.6)

where ${p}_{i}$ is the numerical solution and ${\stackrel{\xaf}{p}}_{i}$ is the exact solution. Then we calculate the order of accuracy. A small time-step $\tau =1\times {10}^{-5}$ and ${v}_{0}=2$ are chosen and the problem is recalculated from time $t=0.01$ to $t=1$ . As illustrated in Table 1, the scheme is approximated second order in the velocity direction.

Figure 1. Transition probability distribution $p(v,t|{v}_{0},0)$ 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 $t=1$ and $\tau =1\times {10}^{-5}$ .

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).

References

[1] Touchette, H., Erik, V.D.S. and Just, W. (2010) Brownian Motion with Dry Friction: Fokker-Planck Approach. Journal of Physics A Mathematical & Theoretical, 43, 3045-3067. https://doi.org/10.1088/1751-8113/43/44/445002

[2] Browne, S., Whitt, W. and Dshalalow, J.H. (1995) Piecewise-Linear Diffusion Processes. Advances in Queuing: Theory, methods, and Open Problems, 4, 463-480.

[3] Gennes, P.G.D. (2005) Brownian Motion with Dry Friction. Journal of Statistical Physics, 119, 953-962. https://doi.org/10.1007/s10955-005-4650-4

[4] Risken, H. (1984) The Fokker-Planck Equation: Methods of Solution and Applica-tions. Springer-Verlag, Berlin. https://doi.org/10.1007/978-3-642-96807-5

[5] Gardiner, C.W. (1985) Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences. 2nd Edition, Springer-Verlag, Berlin.
https://doi.org/10.1007/978-3-662-02452-2

[6] Leveque, R.J. and Li, Z. (1994) The Immersed Interface Method for Elliptic Equa-tions with Discontinuous Coefficients and Singular Sources. SIAM Journal on Numerical Analysis, 31, 1019-1044. https://doi.org/10.1137/0731054

[7] Li, Z., McTigue, D.F. and Heine, J.T. (1997) Short Communication: A Numerical Method for Diffusive Transport with Moving Boundaries and Discontinuous Material Properties. International Journal for Numerical & Analytical Methods in Geomechanics, 21, 653-662.
https://doi.org/10.1002/(SICI)1096-9853(199709)21:9%3C653::AID-NAG894%3E3.0.CO;2-5

[8] Li, Z. (1997) Immersed Interface Methods for Moving Interface Problems. Numerical Algorithms, 14, 269-293. https://doi.org/10.1023/A:1019173215885

[9] Li, Z., Wang, D. and Zou, J. (1998) Theoretical and Numerical Analysis on a Thermo-Elastic System with Discontinuities. Journal of Com-putational and Applied Mathematics, 92, 37-58. https://doi.org/10.1016/S0377-0427(98)00044-2

[10] Wiegmann, A. and Bube, K.P. (1998) The Immersed Interface Method for Nonlinear Differential Equations with Discontinuous Coefficients and Singular Sources. SIAM Journal on Numerical Analysis, 35, 177-200.
https://doi.org/10.1137/S003614299529378X

[11] Li, Z., Jaiman, R.K. and Khoo, B.C. (2016) An Immersed Interface Method for Flow Past Circular Cylinder in the Vicinity of a Plane Moving Wall. International Journal for Numerical Methods in Fluids, 81, 611-639. https://doi.org/10.1002/fld.4198

[12] Li, Z. (2003) An Overview of the Immersed Interface Method and its Applications. Taiwanese Journal of Mathematics, 7, 1-49.
https://doi.org/10.11650/twjm/1500407515

[13] Li, Z. and Ito, K. (2006) The Immersed Interface Method: Numerical Solutions of PDEs Involving Interfaces and Irregular Domains. Society for Industrial and Applied Mathematics, Philadelphia. https://doi.org/10.1137/1.9780898717464

[14] Chen, Y., Baule, A., Touchette, H. and Just W. (2013) Weak-Noise Limit of a Piece-wise-Smooth Stochastic Differential Equation. Physical Review E, 88, 052103.
https://doi.org/10.1103/PhysRevE.88.052103