⇦Chinese Physics Letters, 2020, Vol. 37, No. 6, Article code 060201 Solution to the Fokker–Planck Equation with Piecewise-Constant Drift * Bin Cheng (程彬)1, Ya-Ming Chen (陈亚铭)2**, Xiao-Gang Deng (邓小刚)2,3 Affiliations 1College of Computer, National University of Defense Technology, Changsha 410073, China 2College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China 3Chinese Academy of Military Science, Beijing 100071, China Received 4 February 2020, online 26 May 2020 *Supported by the National Natural Science Foundation of China (Grant Nos. 11972370 and 61772542).
**Corresponding author. Email: chenym-08@163.com Citation Text: Cheng B, Chen Y M and Deng X G 2020 Chin. Phys. Lett. 37 060201    Abstract We study the solution to the Fokker–Planck equation with piecewise-constant drift, taking the case with two jumps in the drift as an example. The solution in Laplace space can be expressed in closed analytic form, and its inverse can be obtained conveniently using some numerical inversion methods. The results obtained by numerical inversion can be regarded as exact solutions, enabling us to demonstrate the validity of some numerical methods for solving the Fokker–Planck equation. In particular, we use the solved problem as a benchmark example for demonstrating the fifth-order convergence rate of the finite difference scheme proposed previously [Chen Y and Deng X Phys. Rev. E 100 (2019) 053303]. DOI:10.1088/0256-307X/37/6/060201 PACS:02.60.Cb, 52.65.Ff, 02.30.Uu © 2020 Chinese Physics Society Article Text Piecewise-smooth stochastic systems are used as models in a wide range of applications, from macroscopic earthquake motions[1,2] to microscopic Brownian motions.[3–5] Many of them can be modelled by stochastic differential equations (SDEs) with piecewise-smooth drifts. The interplay between noise and discontinuities may result in some distinctive features that are different from the smooth cases.[6–13] Particularly we study the model described by the Langevin equation[14] $$ \dot{v}(t)=\varPhi(v)+\sqrt{2D}\xi(t),~~ \tag {1} $$ where the dot stands for the time derivative, $\varPhi(v)$ is the drift that is piecewise-smooth, and $D>0$ denotes the strength of the Gaussian white noise $\xi(t)$ characterized by the zero mean $\langle\xi(t)\rangle=0$ and the correlation $\langle\xi(t)\xi(t^\prime)\rangle=\delta(t-t^\prime)$. Here $\delta$ represents the Dirac delta function. For the initial value $v(0)=v_0$, let us denote the propagator (or transition probability distribution) of the process $v(t)$ as $p(v,t|v_0,0)$, which obeys the following Fokker–Planck equation[15] $$ \partial_t p(v,t|v_0,0) =- \partial_v [{\varPhi}(v) p] +D \partial_v^2 p,~~ \tag {2} $$ where the initial value condition is $$ p(v,0|v_0,0)=\delta(v-v_0).~~ \tag {3} $$ One of the simplest piecewise-smooth cases is the model of Brownian motion with pure dry (also called solid or Coulomb) friction,[4] where the drift is a sign function of $v$. In that case, the propagator can be obtained in closed analytic form using spectral decomposition methods[16] or Laplace transforms[17,18] to solve Eq. (2). Some related analytic results are also available for other simple cases, including the cases with two-value drift[19,20] or piecewise-linear drift.[16,21,22] However, to the best of our knowledge, it is not possible to find in the literature an analytic propagator of the case with multiple jumps in the drift. Therefore, it is meaningful theoretically to derive a corresponding analytic result. On the other hand, an analytic result is necessary for us to check the validity of some numerical methods that can be applied to solving the case with multiple jumps, including the traditional Euler–Maruyama scheme,[23,24] the so-called exact simulation sampling method[25,26] and finite different schemes.[27,28] We consider the case with two jumps studied in Ref. [26]: $$ \varPhi(v)= \left\{ \begin{array}{lll} 0, && v < 0, \\ 1, && 0 < v < 1,\\ 0, && v>1. \end{array} \right.~~ \tag {4} $$ As we will apply matching conditions at each position of the jumps [see, e.g., Eqs. (8) and (9)] to derive a unique solution to Eq. (2), the values of $\varPhi(v)$ at $v=0,1$ are unnecessary here for deriving the analytic solution. However, for convenience of choosing an initial condition for the numerical scheme considered later [see Eq. (17)], we will set the value at each jump as the average of the left and right limits, for example, $\varPhi(0)=[ {\varPhi}(0^-)+\varPhi(0^+) ]/2$. Like the case with a single jump in the drift, we may apply a spectral decomposition method to solve the problem (2)-(4). However, the computation would be too cumbersome. Here we alternatively solve the considered problem using the Laplace transform method. Generally speaking, Eq. (2) with piecewise linear drift can be solved easier in Laplace space than in time space. Moreover, the correlation function of the associated stochastic process (1) can also be studied more conveniently in Laplace space.[18] For the case with a single jump, it can be traced back to Ref. [17] to find a result of using the method of Laplace transform to obtain the propagator of Eq. (2). Related results can also be found in Refs. [18,29,30]. For the case with two jumps (4) considered here, it is noted that the method can be applied similarly. In the following, we present the solution to the problem (2)-(4) following the procedure presented in Ref. [17]. Solutions in Laplace space.—Denote the Laplace transform of the propagator $p(v,t|v_0,0)$ with respect to time $t$ as $$ \widetilde{p}(v,s)=\int_0^\infty p e^{\rm -s t} dt.~~ \tag {5} $$ Then using the differentiation property of the Laplace transform and the initial value condition (3) we obtain from Eq. (2) that $$ s\widetilde{p}(v,s)-\delta(v-v_0) =-\frac{d}{d v}[{\varPhi}(v)\widetilde{p}]+D \frac{d^2}{d v^2} \widetilde{p}.~~ \tag {6} $$ Since the above equation contains the Dirac delta function and the piecewise-constant drift (4), we must pay special attention to the corresponding singular point $v=v_0$ and discontinuous points $v=0$ and $v=1$. To solve Eq. (6), it is necessary to compare the value of $v_0$ with the two discontinuous points. Thus, we consider the three cases: (i) $v_0 < 0$, (ii) $0 < v_0 < 1$, (iii) $v_0>1$. In each case, we study separately the four subdomains divided by the three special points. Then we impose some matching conditions at the interfaces of the subdomains to determine the solution uniquely. It is noted that the Laplace transform $\tilde{p}(v,s)$ (5) is continuous with respect to $v$, since it inherits this property from the propagator $p(v,t|v_0,0)$ (2). Thus, we have the following continuous conditions at the interfaces of the subdomains: $$\begin{align} \tilde{p}(v_0^-,s)&=\tilde{p}(v_0^+,s),~~ \tag {7} \end{align} $$ $$\begin{align} \tilde{p}(0^-,s)&=\tilde{p}(0^+,s),~~ \tag {8} \end{align} $$ $$\begin{align} \tilde{p}(1^-,s)&=\tilde{p}(1^+,s),~~ \tag {9} \end{align} $$ where the superscripts $-$ and $+$ indicate the left and right limits, respectively. For the similar reason, at the positions of the jumps we have the continuous conditions of the current: $$ \tilde{f}(v,s)=-{\varPhi}(v)\tilde{p}(v,s)+D\frac{d}{dv}\tilde{p}(v,s),~~ \tag {10} $$ namely, $$\begin{align} \tilde{f}(0^-,s)&=\tilde{f}(0^+,s),~~ \tag {11} \end{align} $$ $$\begin{align} \tilde{f}(1^-,s)&=\tilde{f}(1^+,s).~~ \tag {12} \end{align} $$ In addition, at the singular point $v=v_0$ we have $$ \tilde{f}(v_0^+,s)-\tilde{f}(v_0^-,s)=-1,~~ \tag {13} $$ which can be obtained by integrating Eq. (6) over the interval $(v_0 -\varepsilon,v_0 +\varepsilon)$ and considering the limit $\varepsilon\rightarrow 0$. We first consider the case with $v_0 < 0$. In this case, we solve Eq. (6) separately for the four intervals $(-\infty,v_0)$, $(v_0,0)$, $(0,1)$ and $(1,\infty)$. By using the vanishing condition $\tilde{p}(\pm \infty, s)=0$ for $s$ with sufficiently large real part, it is not difficult to find the solution to Eq. (6) as follows: $$ \tilde{p}= \begin{cases} A_1e^{\sqrt{2\,s}v}, ~~~~v\in (-\infty,v_0),\\ A_2e^{\sqrt{2\,s}v} +B_2e^{-\sqrt{2\,s}v}, ~~~~v\in (v_0,0), \\ A_3e^{(1+\sqrt{1+2\,s})v} \!+\!B_3e^{-(1+\sqrt{1+2\,s})v}, ~~\!v\in (0,1), \\ B_4e^{-\sqrt{2\,s}v}, ~~~~v\in (1,\infty), \end{cases}~~ \tag {14} $$ where the coefficients $A_i$ ($i=1,2,3$) and $B_i$ ($i=2,3,4$) are determined by the matching conditions (7)-(9) and (11)-(13), resulting in six linear equations that can be solved easily. Here the detailed expressions of the coefficients are not presented due to the length constraint. In fact, it is unnecessary to write them down here since the expressions can be easily obtained using some symbolic computing software. We only need to know that these coefficients depend on $s$ while they are independent of $v$. This statement applies for the following two cases as well. For the second case with $0 < v_0 < 1$, we shall consider the four intervals $(-\infty,0)$, $(0,v_0)$, $(v_0,1)$ and $(1,\infty)$, respectively. The corresponding solution to Eq. (6) can be expressed as $$ \tilde{p}\!=\!\! \begin{cases} \!\!\!\!A_1e^{\sqrt{2\,s}v}, ~~ v\in (-\infty,0),\\ \!\!\!\!A_2e^{(1+\sqrt{1+2\,s})v} \!+\!B_2e^{-(1+\sqrt{1+2\,s})v}, ~~v\!\!\in\!\! (0,\!v_0), \\ \!\!\!\!A_3e^{(1+\sqrt{1+2\,s})v} \!+\!B_3e^{-(1+\sqrt{1+2\,s})v}, ~~v\!\!\in\!\! (v_0,\!1), \\ \!\!\!\!B_4e^{-\sqrt{2\,s}v}, ~~ v\in (1,\infty). \end{cases}~~ \tag {15} $$ Similarly to the previous two cases, we consider the four intervals $(-\infty,0)$, $(0,1)$, $(1,v_0)$ and $(v_0,\infty)$ separately for the case with $v_0>1$. Then we solve Eq. (6) to obtain $$ \tilde{p}= \begin{cases} \!\!\!A_1e^{\sqrt{2\,s}v}, ~~~~ v\in (-\infty,0),\\ \!\!\!A_2e^{(1+\sqrt{1+2\,s})v} \!+\!B_2e^{-(1+\sqrt{1+2\,s})v}, ~~v\!\in\! (0,1), \\ \!\!\!A_3e^{\sqrt{2\,s}v} +B_3e^{-\sqrt{2\,s}v}, ~~~~ v\in (1,v_0), \\ \!\!\!B_4e^{-\sqrt{2\,s}v}, ~~~~ v\in (v_0,\infty). \end{cases}~~ \tag {16} $$ To summarize, depending on the values of $v_0$, the propagator of the problem (2) and (3) in Laplace space is given by Eqs. (14)-(16). Solutions in time space.—Due to the complexity of the expressions (14)-(16) and their corresponding coefficients, we are unable to perform an inverse Laplace transform analytically to obtain the solution in time space. Thus we resort to some numerical inversion algorithm. Here, as was applied successfully in Ref. [18], we use the so-called Talbot method.[31–33] Although it may take time to understand this method, it is actually very convenient to realize this method by using a corresponding Mathematica implementation (available at https://library.wolfram.com/infocenter/MathSource/5026/). Corresponding to the three above-mentioned cases, we choose the initial value to be $v_0=-1$, $v_0=0.5$ and $v_0=2$, respectively, and produce the profiles of the propagator as shown in Fig. 1. As is expected, for each case there are two kinks appearing at the positions of the jumps in the drift. However, the propagators are smooth at $v=v_0$ although the expressions in Laplace space look different; see Eqs. (14)-(16). It should be explained that in Fig. 1 we just chose the times suitable to illustrate the kinks in the propagators. One could certainly produce the propagators at other times conveniently.
cpl-37-6-060201-fig1.png
Fig. 1. Propagators of the problem (2)-(4) for different values of $v_0$. The "exact" solutions denote the numerical inversion of Eqs. (14)-(16), while the "numeric" results are obtained using the finite different scheme presented in Ref. [28]. Here the numeric results are produced using the numbers of solution points $N_1=150$, $N_2=10$ and $N_3=140$ in the three subdomains divided by $v=0$ and $v=1$. (a) The case with $v_0=-1$ at time $t=2$; (b) the case with $v_0=0.5$ at time $t=1$; (c) the case with $v_0=2$ at time $t=4$.
It is worth noting that the numerical results obtained from Eqs. (14)-(16) can be regarded as exact solutions since the used numerical inversion is very accurate. Therefore, the obtained results can serve as reference solutions for testing the validity of some numerical methods. Here we are interested in studying the fifth-order finite difference scheme proposed in Ref. [28], which is suitable for solving the Fokker–Planck equation (2) with a single jump or multiple jumps in the drift. For the case with a single jump, the fifth-order convergence rate of the scheme has been demonstrated using the analytic propagator of Brownian motion with pure dry friction. However, the validity of the scheme for the case with two jumps is only checked by showing its convergence numerically. On the one hand, it is not very clear whether the numerical solutions converge to the right one since it may happen that a convergent scheme would result in a wrong solution; see, e.g., Refs. [34,35]. On the other hand, it is not certain that the scheme achieves the design fifth-order convergence rate for the case with two jumps. Therefore, it is necessary to use the obtained solutions from Eqs. (14)-(16) to check the validity of the scheme. To develop the finite difference scheme presented in Ref. [28] we have to truncate the computational domain of $v$ to be finite, and then divide it into three subdomains by using the discontinuous points $v=0$ and $v=1$. For the three subdomains, we denote the numbers of solution points as $N_1$, $N_2$ and $N_3$ subsequently. We do not attempt to introduce the details of the scheme here as one can consult the mentioned reference to find more information. To implement the scheme, we still set the initial condition to be Gaussian as Ref. [28], i.e., $$ p\left(v, \tau_{0} | v_{0}, 0\right)=\frac{1}{\sqrt{4 \pi D \tau_{0}}} e^{\left[v-v_{0}-{\varPhi}\left(v_{0}\right) \tau_{0}\right]^{2} /\left(4 D \tau_{0}\right)}.~~ \tag {17} $$ Here $D=0.5$ and $\tau_{0}=0.01$ are chosen. Actually, we would obtain essentially the same numerical result if we use the numerical Laplace inversion of Eqs. (14)-(16) as an initial condition. In addition, we truncate the computational domain to be $[-15,15]$ and set the current $$ f(v,t)=-\varPhi(v) p+D \partial_v p~~ \tag {18} $$ of Eq. (2) to be zero at the boundaries. As that in Ref. [28], the third-order Runge–Kutta scheme with the same time step is used.
Table 1. Accuracy test of the finite difference scheme presented in Ref. [28] for different values of $v_0$. Here $N_i$ ($i=1,2,3$) are the numbers of solution points in the subdomains divided by $v=0$ and $v=1$, and the total solution points $N=N_1+N_2+N_3-2$ since the solution points of different subdomains overlap at the interfaces. The $L^{\infty}$ error is defined as the same as that in Ref. [27] [see Eq. (33) therein]. The convergence rate is computed by $-\ln(E_M/E_K)/\ln(M/K)$ with $E_M$ and $E_K$ denoting the corresponding errors of the cases with $M$ and $K$ total solution points, respectively.
$v_0=-1$, $t=2$ $v_0=0.5$, $t=1$ $v_0=2$, $t=4$
$N_1$ $N_2$ $N_3$ $N$ $L^{\infty}$ error Rate $L^{\infty}$ error Rate $L^{\infty}$ error Rate
75 5 70 148 8.95$\times 10^{-3}$ 3.82$\times 10^{-1}$ 6.73$\times 10^{-3}$
150 10 140 298 8.49$\times 10^{-6}$ 9.95 4.29$\times 10^{-3}$ 6.41 1.86$\times 10^{-5}$ 8.42
300 20 280 598 3.32$\times 10^{-7}$ 4.66 2.50$\times 10^{-5}$ 7.39 7.49$\times 10^{-7}$ 4.61
600 40 560 1198 1.12$\times 10^{-8}$ 4.87 1.63$\times 10^{-7}$ 7.25 2.56$\times 10^{-8}$ 4.86
1200 80 1120 2398 3.64$\times 10^{-10}$ 4.94 2.89$\times 10^{-9}$ 5.81 8.30$\times 10^{-10}$ 4.94
As we can see from Fig. 1, the results of the finite difference scheme match well with those of the corresponding Laplace inversion. By comparing the results of these two methods, we also obtain the test of accuracy as illustrated in Table 1, showing that the design fifth-order convergence rate is achieved for different values of $v_0$, as expected. This indicates that the finite difference scheme does produce the right solution at a fifth-order convergence rate. Therefore, the validity of the scheme is verified. In summary, we have provided the exact solution in Laplace space of the Fokker–Planck equation with piecewise-constant drift that admits two jumps. To solve the considered problem, the domain is divided into four subdomains by the initial value $v_0$ and the positions of the two jumps. Then the expressions of the propagator are obtained separately in different subdomains. Finally, at the interfaces of the subdomains six matching conditions are provided to determine the solution in Laplace space uniquely. To recover the propagator in time space, the very efficient Talbot method is recommended. From the results, we observe that at the position of each jump there is a kink in the propagator. Regarding the results obtained from Laplace inversion as exact solutions, we are able to demonstrate the validity of some numerical methods that can handle the case with two jumps in the drift. In particular, we show that the finite difference scheme proposed in Ref. [28] converges to the right solution at a fifth-order convergence rate, indicating the validity of the scheme. Furthermore, we are confident that the considered scheme is suitable for use in investigating the solutions of the Fokker–Planck equations with more complex non-smooth drifts. References Accumulated Slip of a Friction-Controlled Mass Excited by Earthquake MotionsBiaxial slip of a mass on a foundation subjected to earthquake motionsBrownian motors: noisy transport far from equilibriumBrownian Motion with Dry FrictionGranular Brownian motion with dry frictionStick–slip motion of solids with dry friction subject to random vibrations and an external fieldEffect of Coulombic friction on spatial displacement statisticsSingular features in noise-induced transport with dry frictionRectification of asymmetric surface vibrations with dry friction: An exactly solvable modelFirst-passage time of Brownian motion with dry frictionLarge-deviation properties of Brownian motion with dry frictionNonequilibrium dynamics of a pure dry friction model subjected to colored noiseSome new advance on the research of stochastic non-smooth systemsLangevin equation with Coulomb frictionBrownian motion with dry friction: Fokker–Planck approachAnalysis of a Nonlinear First‐Order System with a White Noise InputExact power spectra of Brownian motion with solid frictionTrivariate Density of Brownian Motion, Its Local and Occupation Times, with Application to Stochastic ControlA path integral approach to random motion with nonlinear frictionWeak-noise limit of a piecewise-smooth stochastic differential equationA numerical method for SDEs with discontinuous driftStrong convergence for the Euler–Maruyama approximation of stochastic differential equations with discontinuous coefficientsExact sampling of diffusions with a discontinuity in the driftExact Simulation of Brownian Diffusions with Drift Admitting JumpsNumerical solutions of Fokker-Planck equations with drift-admitting jumpsFifth-order finite-difference scheme for Fokker-Planck equations with drift-admitting jumpsSpectral density of piecewise linear first order systems excited by white noiseFirst order piecewise linear systems with random parametric excitationThe Accurate Numerical Inversion of Laplace TransformsMulti-precision Laplace transform inversionA Unified Framework for Numerically Inverting Laplace TransformsAN ANALYSIS OF THREE DIFFERENT FORMULATIONS OF THE DISCONTINUOUS GALERKIN METHOD FOR DIFFUSION EQUATIONSAnalysis and Comparison of Different Approximations to Nonlocal Diffusion and Linear Peridynamic Equations
[1] Crandall S H, Lee S S and Williams Jr J H 1974 J. Appl. Mech. 41 1094
[2] Crandall S H and Lee S S 1976 Ingenieur-Arch. 45 361
[3] Reimann P 2002 Phys. Rep. 361 57
[4] De Gennes P G 2005 J. Stat. Phys. 119 953
[5] Gnoli A, Puglisi A and Touchette H 2013 Europhys. Lett. 102 14002
[6] Baule A, Touchette H and Cohen E G D 2011 Nonlinearity 24 351
[7] Menzel A M and Goldenfeld N 2011 Phys. Rev. E 84 011122
[8] Baule A and Sollich P 2012 Europhys. Lett. 97 20001
[9] Baule A and Sollich P 2013 Phys. Rev. E 87 032112
[10] Chen Y and Just W 2014 Phys. Rev. E 89 022103
[11] Chen Y and Just W 2014 Phys. Rev. E 90 042102
[12] Geffert P M and Just W 2017 Phys. Rev. E 95 062111
[13] Xu W, Wang L, Feng J, Qiao Y and Han P 2018 Chin. Phys. B 27 110503
[14] Hayakawa H 2005 Physica D 205 48
[15]Risken H 1989 The Fokker–Planck Equation: Methods of Solution and Applications (Berlin: Springer)
[16] Touchette H, Der Straeten E V and Just W 2010 J. Phys. A 43 445002
[17] Caughey T K and Dienes J K 1961 J. Appl. Phys. 32 2476
[18] Touchette H, Prellberg T and Just W 2012 J. Phys. A 45 395002
[19] Karatzas I and Shreve S E 1984 Ann. Probab. 12 819
[20]Simpson D J W and Kuske R 2014 Disc. Contin. Dyn. Syst. B 19 2889
[21] Baule A, Cohen E G D and Touchette H 2010 J. Phys. A 43 025003
[22] Chen Y, Baule A, Touchette H and Just W 2013 Phys. Rev. E 88 052103
[23] Leobacher G and Szölgyenyi M 2016 BIT Numer. Math. 56 151
[24] Ngo H L and Taguchi D 2017 Stat. Probab. Lett. 125 55
[25] Papaspiliopoulos O, Roberts G O and Taylor K B 2016 Adv. Appl. Probab. 48 249
[26] Dereudre D, Mazzonetto S and Roelly S 2017 SIAM J. Sci. Comput. 39 A711
[27] Chen Y and Deng X 2018 Phys. Rev. E 98 033302
[28] Chen Y and Deng X 2019 Phys. Rev. E 100 053303
[29] Atkinson J D and Caughey T K 1968 Int. J. Non-Linear Mech. 3 137
[30] Atkinson J D and Caughey T K 1968 Int. J. Non-Linear Mech. 3 399
[31] Talbot A 1979 IMA J. Appl. Math. 23 97
[32] Abate J and Valkó P P 2004 Int. J. Numer. Methods Eng. 60 979
[33] Abate J and Whitt W 2006 INFORMS J. Comput. 18 408
[34] Zhang M and Shu C W 2003 Math. Mod. Meth. Appl. Sci. 13 395
[35] Tian X and Du Q 2013 SIAM J. Numer. Anal. 51 3458