Numerical methods for the two-dimensional Fokker-Planck equation governing the probability density function of the tempered fractional Brownian motion

Abstract

In this paper, we study the numerical schemes for the two-dimensional Fokker-Planck equation governing the probability density function of the tempered fractional Brownian motion. The main challenges of the numerical schemes come from the singularity in the time direction. When 0 < H < 0.5, a change of variables $\partial \left (t^{2H}\right )=2Ht^{2H-1}\partial t$ avoids the singularity of numerical computation at t = 0, which naturally results in nonuniform time discretization and greatly improves the computational efficiency. For H > 0.5, the time span dependent numerical scheme and nonuniform time discretization are introduced to ensure the effectiveness of the calculation and the computational efficiency. The stability and convergence of the numerical schemes are demonstrated by using Fourier method. By numerically solving the corresponding Fokker-Planck equation, we obtain the mean squared displacement of stochastic processes, which conforms to the characteristics of the tempered fractional Brownian motion.

Appendices

Appendix A: Numerical stability

For 0 < H < 0.5, the Fourier series of u^k(x,y) is

$$ u^{k}(x,y)=\sum\limits_{p_{1}=-\infty}^{+\infty}\sum\limits_{p_{2}=-\infty}^{+\infty} \hat{u}_{p_{1},p_{2}}^{k}\exp\left( i\frac{2p_{1}\pi x}{L}+i\frac{2p_{2}\pi y}{L^{\prime}}\right), $$

(A.1)

where

$$ \hat{u}_{p_{1},p_{2}}^{k}=\frac{1}{LL^{\prime}}{{\int}_{0}^{L}}{\int}_{0}^{L^{\prime} } u^{k}(x,y)\exp\left( -i\frac{2p_{1}\pi x}{L}-i\frac{2p_{2}\pi y}{L^{\prime}}\right)\mathrm{d} x \mathrm{d} y\quad p_{1},p_{2}=0,\pm1,{}\cdots. $$

There exists Parseval equation

$$ \left\|u^k(x,y)\right\|_{L^2}^2=LL^{\prime}\sum\limits_{p_1=-\infty}^{+\infty} \sum\limits_{p_2=-\infty}^{+\infty}\left|\hat{u}_{p_1,p_2}^k\right|^2. $$

From (2.4), we get

$$ \begin{array}{l} u^{k+1}\left( x+x_{m},y+y_{n}\right)-u^{k}\left( x+x_{m},y+y_{n}\right)\\ =\frac{r}{h^{2}}{\delta_{x}^{2}}u^{k+1}\left( x+x_{m},y+y_{n}\right)^{k+1}+ \frac{r}{l^{2}}{\delta_{y}^{2}}u^{k+1}\left( x+x_{m},y+y_{n}\right). \end{array} $$

(A.2)

Substituting (A.1) into (A.2) leads to

$$ \begin{array}{l} \sum\limits_{p_{1}=-\infty}^{+\infty}\sum\limits_{p_{2}=-\infty}^{+\infty} \hat{u}_{p_{1},p_{2}}^{k}Q\left( p_{1},p_{2}\right)\\ =\sum\limits_{p_{1}=-\infty}^{+\infty}\sum\limits_{p_{2}=-\infty}^{+\infty} \hat{u}_{p_{1},p_{2}}^{k+1}Q\left( p_{1},p_{2}\right) \left\{\left( 1+\frac{2r}{h^{2}}+\frac{2r}{l^{2}}\right)\right.\\ ~~~~\left.-\frac{r}{h^{2}} \left[\exp\left( i\frac{2p_{1}\pi h}{L}\right)+\exp\left( -i\frac{2p_{1}\pi h}{L}\right)\right] - \frac{r}{l^{2}}\left[\exp\left( i\frac{2p_{2}\pi l}{L^{\prime}}\right)+\exp\left( -i\frac{2p_{2}\pi l}{L^{\prime}}\right)\right]\right\}, \end{array} $$

(A.3)

where

$$ Q\left( p_1,p_2\right)=\exp\left( i\frac{2p_1\pi x}{L}+i\frac{2p_2\pi y}{L^{\prime}}\right)\exp\left( i\frac{2p_1\pi m h}{L}+i\frac{2p_2\pi n l}{L^{\prime}}\right). $$

Since the two sides of (A.3) are the Fourier series, we have

$$ \hat{u}_{p_{1},p_{2}}^{k+1}=G_{1}\left( p_{1}h,p_{2}l\right)\hat{u}_{p_{1},p_{2}}^{k}, $$

(A.4)

where

$$ G_1\left( p_1h,p_2l\right)=\frac{1}{1+\frac{2r}{h^2}\left( 1-\cos\frac{2p_1\pi h}{L}\right)+\frac{2r}{l^2}\left( 1-\cos\frac{2p_2\pi l}{L^{\prime}}\right)}. $$

This implies that

$$ 0\le G_1\left( p_1h,p_2l\right)\le1. $$

Combining Parseval equation and (A.4) results in

$$ \begin{array}{@{}rcl@{}} \left\|u^{k}(x,y)\right\|_{L^{2}}^{2} &=&LL^{\prime}\sum\limits_{p_{1}=-\infty}^{+\infty}\sum\limits_{p_{2}=-\infty}^{+\infty} \left|\hat{u}_{p_{1},p_{2}}^{k}\right|^{2}\\ &<&\left\|u^{0}(x,y)\right\|_{L^{2}}^{2}. \end{array} $$

As H > 0.5, for $t\ge t_{\max \nolimits }=t_{k_{1}}$, using the same process, we have

$$ \hat{u}_{p_{1},p_{2}}^{k+1}=G_{2}\left( p_{1}h,p_{2}l\right)\hat{u}_{p_{1},p_{2}}^{k}, $$

(A.5)

where

$$ G_2\left( p_1h,p_2l\right)=\frac{1}{1+\frac{2r_1}{h^2}\left( 1-\cos\frac{2p_1\pi h}{L}\right)+\frac{2r_1}{l^2}\left( 1-\cos\frac{2p_2\pi h}{L^{\prime}}\right)} $$

and

$$ r_1=\frac{{\varGamma}(H+1/2)}{H\sqrt{\pi}(2\lambda)^H} \left[\lambda t_{k+1}K_{H-1}(\lambda t_{k+1})\tau\right]. $$

For k ≤ k₁, with the proof being completely the same as the case that 0 < H < 0.5, there exists

$$ \left\|u^{k}(x,y)\right\|_{L^{2}}^{2}<\left\|u^{0}(x,y)\right\|_{L^{2}}^{2}. $$

(A.6)

For k > k₁, combining (A.5) and (A.6) leads to

$$ \begin{array}{@{}rcl@{}} \left\|u^{k}(x,y)\right\|_{L^{2}}^{2} &=&LL^{\prime}\sum\limits_{p_{1}=-\infty}^{+\infty}\sum\limits_{p_{2}=-\infty}^{+\infty} \left[G_{2}\left( p_{1}h,p_{2}l\right)\right]^{2k-2k_{1}}\left|\hat{u}_{p_{1},p_{2}}^{k_{1}}\right|^{2}\\ &<&\left\|u^{0}(x,y)\right\|_{L^{2}}^{2}. \end{array} $$

Appendix B: Convergence

We use notations

$$ \begin{array}{l} Lu(x,y,t)=\frac{\partial u(x,y,t)}{\partial \left( t^{2H}\right)}-\frac{{\varGamma}(H+1/2)\lambda t^{1-H}K_{H-1}(\lambda t)}{2H\sqrt{\pi}(2\lambda)^{H}} \left[\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}\right]u(x,y,t),\\ L^{(1)}u_{m,n}^{k}=\frac{u_{m,n}^{k+1}-u_{m,n}^{k}}{{\varDelta} \left( t_{k}^{2H}\right)}-\frac{{\varGamma}(H+1/2)\lambda t_{k+1}^{1-H}K_{H-1}(\lambda t_{k+1})}{2H\sqrt{\pi}(2\lambda)^{H}}\left[\frac{{\delta_{x}^{2}}}{h^{2}}+\frac{{\delta_{y}^{2}}}{l^{2}}\right]u_{m,n}^{k+1}. \end{array} $$

As 0 < H < 0.5, performing the Taylor expansion at $t_{k}^{2H}$, there exist

$$ \begin{array}{l} \frac{{\varDelta} t_{k}}{{\varDelta}\left( t_{k}^{2H}\right)} =\frac{t_{k+1}-t_{k}}{\tau}\\ =\frac{t_{k}^{1-2H}}{2H}-\frac{(1-2H)t_{k}^{1-4H}}{8H^{2}}\tau+O\left( \tau^{2}\right) \end{array} $$

(B.1)

and

$$ \frac{\left( {\varDelta} t_{k}\right)^{2}}{{\varDelta} \left( t_{k}^{2H}\right)}=\frac{t_{k}^{2-4H}}{4H^{2}}\tau+O\left( \tau^{2}\right). $$

(B.2)

Letting $R_{m,n}^{k}=L^{(1)}u_{m,n}^{k}-[Lu(x,y,t)]_{m,n}^{k},$ and using (B.1) and (B.2) lead to

$$ \begin{array}{l} R_{m,n}^{k}=-\frac{(1-2H)t_{k}^{1-4H}}{8H^{2}}\tau\left( \frac{\partial u(x,y,t)}{\partial t}\right)_{m,n}^{k}-\frac{t_{k}^{2-4H}}{8H^{2}}\tau\left( \frac{\partial^{2} u(x,y,t)}{\partial t^{2}}\right)_{m,n}^{k}+O\left( \tau^{2}+h^{2}+l^{2}\right)\\ =O\left( \tau+h^{2}+l^{2}\right) \end{array}. $$

(B.3)

For $e_{m,n}^{k}=u_{m,n}^{k}-u(x_{m},y_{n},t_{k})$, from (1.3), (2.4), and (B.3), we have

$$ \begin{array}{l} e^{k+1}\left( x+x_{m},y+y_{n}\right)-e^{k}\left( x+x_{m},y+y_{n}\right)\\ =r\left[\frac{{\delta_{x}^{2}} }{h^{2}}+\frac{{\delta_{y}^{2}} }{l^{2}}\right]e^{k+1}\left( x+x_{m},y+y_{n}\right)+\tau R^{k}\left( x+x_{m},y+y_{n}\right). \end{array} $$

Following the proof process of numerical stability and using the expansion similar to (A.3), there exists

$$ \left\|e^{k+1}\right\|_{L^2}^2<\left\|e^k+\tau R^k\right\|_{L^2}^2, $$

leading to

$$ \begin{array}{@{}rcl@{}} \left\|e^{k}\right\|_{L^{2}}&<&\left\|e^{k-1}\right\|_{L^{2}}+\tau \left\|R^{k-1}\right\|_{L^{2}}\\ &\le&\left\|e^{0}\right\|_{L^{2}}+k\tau \max\limits_{0\le i\le k}\left\|R^{i}\right\|_{L^{2}}\\ &\le&t_{k}^{2H}\max\limits_{0\le i\le k}\left\|R^{i}\right\|_{L^{2}}\\&=&O\left( \tau+h^{2}+l^{2}\right). \end{array} $$

For H > 0.5, when $t>t_{\max \nolimits }$, by Taylor expansion at ${t_{k}^{H}}$, we have

$$ \frac{{\varDelta} t_{k}}{{\varDelta} \left( {t_{k}^{H}}\right)} =\frac{t_{k}^{1-H}}{H}-\frac{(1-H)t_{k}^{1-2H}}{2H^{2}}\tau+O\left( \tau^{2}\right) $$

(B.4)

and

$$ \frac{\left( {\varDelta} t_{k}\right)^{2}}{{\varDelta} \left( {t_{k}^{H}}\right)}=\frac{t_{k}^{2-2H}}{H^{2}}\tau+O\left( \tau^{2}\right), $$

(B.5)

which implies that

$$ \begin{array}{l} R_{m,n}^{k}=-\frac{(1-H)t_{k}^{1-2H}}{2H^{2}}\tau\left( \frac{\partial u(x,y,t)}{\partial t}\right)_{m,n}^{k}-\frac{t_{k}^{2-2H}}{2H^{2}}\tau\left( \frac{\partial^{2} u(x,y,t)}{\partial t^{2}}\right)_{m,n}^{k}+O\left( \tau^{2}+h^{2}+l^{2}\right)\\ =O\left( \tau+h^{2}+l^{2}\right). \end{array} $$

For k ≤ k₁,

$$ \left\|e^{k}\right\|_{L^{2}}<\left\|e^{0}\right\|_{L^{2}}+t_{k}^{2H}\max\limits_{0\le i\le k}\left\|R^{i}\right\|_{L^{2}}=O\left( \tau+h^{2}+l^{2}\right). $$

(B.6)

When k ≥ k₁, combining (B.6) leads to

$$ \begin{array}{@{}rcl@{}} \left\|e^{k}\right\|_{L^{2}}&<&\left\|e^{k_{1}}\right\|_{L^{2}}+\tau (k-k_{1}) \left\|R^{k-1}\right\|_{L^{2}}\\ &\le&\left\|e^{0}\right\|_{L^{2}}+t_{k_{1}}^{2H} \max_{0\le i\le k_{1}}\left\|R^{i}\right\|_{L^{2}}+{t_{k}^{H}} \max_{k_{1}\le i\le k}\left\|R^{i}\right\|_{L^{2}}\\ &\le&\left( {t_{k}^{H}}+t_{\max}^{2H}\right)\max_{0\le i\le k}\left\|R^{i}\right\|_{L^{2}}\\&=&O\left( \tau+h^{2}+l^{2}\right). \end{array} $$