Accurate Boundary Treatment for Riesz Space Fractional Diffusion Equations

Abstract

We consider numerical boundary treatment for solving the Cauchy problems of the Riesz space fractional diffusion equation with compact initial data in one and two space dimension(s). First, the Riesz space fractional equation is semi-discretized into a lattice system. Then we derive an equivalent decoupled form for its dynamics using kernel functions. Series expansions and path integration are devised to numerically evaluate the kernel functions with high accuracy. For the first time, this allows an accurate numerical boundary treatment for the Riesz space fractional diffusion equation. Numerical results demonstrate the effectiveness of the method. The methodology may be extended to treat other fractional partial differential equations.

Appendices

Appendix A: Series Expansions for $g_n(t,\alpha )$ and ${\dot{g}}_0(t,\alpha )$

Now we give the proof of Lemma 1.

Proof

First, we recall that [35]

$$\begin{aligned} \int _{0}^{\pi /2}(\cos \eta )^{\alpha }\cos (\beta \eta )d\eta = \displaystyle \frac{\pi }{2^{\alpha +1}}\frac{\Gamma (1+\alpha )}{\Gamma \left( 1+\displaystyle \frac{\alpha +\beta }{2}\right) \Gamma \left( 1+\displaystyle \frac{\alpha -\beta }{2}\right) }. \end{aligned}$$

(50)

We expand $g_{n}(t,\alpha )$ with respect to t.

$$\begin{aligned} \begin{array}{ll} g_{n}(t,\alpha )&{}=\displaystyle \frac{1}{2\pi }\int _{0}^{2\pi }e^{-\displaystyle \Big (2\sin \frac{\xi }{2}\Big )^{\alpha }t}e^{-in\xi }d\xi \\ &{}=\displaystyle \frac{1}{\pi }\int _{0}^{\pi }e^{-\displaystyle \Big (2\sin \frac{\xi }{2}\Big )^{\alpha }t}\cos (n\xi )d\xi \\ &{}=\displaystyle \frac{1}{\pi }\sum _{k=0}^{+\infty }\frac{(-2^{\alpha }t)^k}{k!}\int _{0}^{\pi }\Big (\sin (\xi /2)\Big )^{k\alpha }\cos (n\xi )d\xi \quad \quad \quad \pi -\xi \rightarrow \xi _1\\ &{}=\displaystyle \frac{(-1)^{n+1}}{\pi }\sum _{k=0}^{+\infty }\frac{(-2^{\alpha }t)^k}{k!}\int _{\pi }^{0}\Big (\cos (\xi _1/2)\Big )^{k\alpha }\cos (n\xi _1)d\xi _1 \quad \quad \quad \frac{\xi _1}{2} \rightarrow \eta \\ &{}=\displaystyle \frac{2(-1)^{n}}{\pi }\sum _{k=0}^{+\infty }\frac{(-2^{\alpha }t)^k}{k!}\int _{0}^{\pi /2}\Big (\cos \eta \Big )^{k\alpha }\cos (2n\eta )d\eta \\ &{}=\displaystyle \frac{2(-1)^{n}}{\pi }\sum _{k=0}^{+\infty }\frac{(-2^{\alpha }t)^k}{k!}\frac{\pi }{2^{k\alpha +1}} \displaystyle \frac{\Gamma (1+k\alpha )}{\Gamma (1+(k\alpha +2n)/2)\Gamma (1+(k\alpha -2n)/2)}\\ &{}=(-1)^{n}\displaystyle \sum _{k=0}^{+\infty }\frac{(-t)^k}{k!}\displaystyle \frac{\Gamma (1+k\alpha )}{\Gamma (1+(k\alpha +2n)/2)\Gamma (1+(k\alpha -2n)/2)},\nonumber \end{array} \end{aligned}$$

(51)

which leads to (20).

For the second expansion, we make change of variable $\eta =2\sin \displaystyle \frac{\xi }{2}$ ($\xi =2\arcsin (\eta /2)$) and then $\theta =2^\alpha - \eta ^\alpha$ to compute

$$\begin{aligned} g_{n}(t,\alpha )&=\displaystyle \frac{1}{\pi }\int _{0}^{\pi }e^{-\displaystyle \Big (2\sin \frac{\xi }{2}\Big )^{\alpha }t}\cos (n\xi )d\xi \nonumber \\&=\displaystyle \frac{1}{\pi }\int _{0}^{2 } e^{-\eta ^{\alpha }t} \displaystyle \frac{\cos (2n\arcsin \eta /2)}{\sqrt{1-(\eta /2)^2}}d\eta \nonumber \\&=\displaystyle \frac{1}{\pi }\Big (\int _{0}^{(2^\alpha -1)^{1/\alpha }}+\int _{ (2^\alpha -1)^{1/\alpha }}^{2}\Big ) e^{-\eta ^{\alpha }t} \displaystyle \frac{\cos (2n\arcsin \eta /2)}{\sqrt{1-(\eta /2)^2}}d\eta \nonumber \\&= \displaystyle \frac{1}{\pi }\int _{0}^{(2^{\alpha }-1)^{1/\alpha }}e^{-\eta ^{\alpha }t}\cos (2n\arcsin \eta /2) \displaystyle \frac{1}{\sqrt{1-(\eta /2)^2}}d\eta \nonumber \\&\quad +\frac{1}{\pi } \int _{0}^{1} e^{-2^{\alpha }t+\theta t} \cos \left( 2n\arcsin \displaystyle \frac{(2^\alpha -\theta )^{1/\alpha }}{2}\right) \nonumber \\&\quad \displaystyle \frac{1}{\sqrt{1-\left( \displaystyle \frac{(2^\alpha -\theta )^{1/\alpha }}{2}\right) ^2}} \displaystyle \frac{1}{\alpha }(2^\alpha -\theta )^{1/\alpha -1} d\theta \nonumber \\&= \displaystyle \frac{1}{\alpha \pi } \sum _{k=0}^{+\infty }a_{k}^{n}\frac{\gamma _{1}( (2^{\alpha }-1)t,(k+1)/\alpha )}{t^{(k+1)/\alpha }} + \displaystyle \frac{1}{\alpha \pi } \displaystyle \sum _{k=0}^{+\infty }b_{k}^{n}\displaystyle \frac{e^{-2^{\alpha }t}\gamma _{2}(t,k-1/2)}{t^{k+1/2}}. \end{aligned}$$

(52)

Here the Taylor expansion for $\displaystyle \frac{\cos (2n\arcsin \eta /2)}{\sqrt{1-(\eta /2)^2}}$ is denoted as $\displaystyle \sum _{k=0}^{+\infty }a_k^n \eta ^{k}$, and the Taylor expansion for $\displaystyle \frac{\sqrt{\theta }\cos \left( 2n\arcsin \displaystyle \frac{(2^\alpha -\theta )^{1/\alpha }}{2}\right) }{\sqrt{1-\left( \displaystyle \frac{(2^\alpha -\theta )^{1/\alpha }}{2}\right) ^2}} (2^\alpha -\theta )^{1/\alpha -1}$ is denoted as $\displaystyle \sum _{k=0}^{+\infty }b_k^n \eta ^{k}$, where a term $\sqrt{\theta }$ is inserted to eliminate the singularity at $\theta =0$. Thus we prove (21). $\blacksquare$

Next,we demonstrate the applicability of (22) in Remark 1 in the case where both n and t are big. In particular, we consider n large enough.

We have

$$\begin{aligned} \begin{array}{ll} g_{n}(t,\alpha )&{}=\displaystyle \frac{1}{2\pi }\int _{0}^{2\pi }e^ {-\displaystyle \Big (2\sin \frac{\xi }{2}\Big )^{\alpha }t}e^{-in\xi }d\xi \\ &{}=\displaystyle \frac{1}{\pi }Re\left[ \int _{0}^{\pi }e^{-\displaystyle \Big (2\sin \frac{\xi }{2}\Big )^{\alpha }t}e^{-in\xi }d\xi \right] . \end{array} \end{aligned}$$

The function $e^{-\displaystyle \Big (2\sin \frac{z}{2}\Big )^{\alpha }t}e^{-inz}$ is analytic, hence we may calculate the above integral over a new path from 0 to $\pi$ in the complex plane. In particular, we consider three subsequent straight line segments, namely, from 0 to $20(1-i)/n$, then from $20(1-i)/n$ to $\pi -20(1+i)/n$, and the last one from $\pi -20(1+i)/n$ to $\pi$. When t is big, the integral on the second and the third segments can be omitted. In fact, because n large enough, on the second segment we have $Re\Big (2\sin \frac{\xi }{2}\Big )^{\alpha }>0$ and

$$\begin{aligned} \Big |e^{-\displaystyle \Big (2\sin \frac{\xi }{2}\Big )^{\alpha }t}e^{-in\xi }\Big | \le e^{-\displaystyle Re\Big (2\sin \frac{\xi }{2}\Big )^{\alpha }t}e^{-20}. \end{aligned}$$

Similarly, integral on the third segment can be omitted as well. So it amounts to

$$\begin{aligned} \begin{array}{ll} \displaystyle \int _{0}^{20(1-i)/n}e^{-\displaystyle \Big (2\sin \frac{\xi }{2}\Big )^{\alpha }t}e^{-in\xi }d\xi &{}=\displaystyle \int _{0}^{2\sin (10(1-i)/n) } \frac{e^{-\eta ^{\alpha }t} e^{-2ni\arcsin \eta /2}}{\sqrt{1-(\eta /2)^2}}d\eta \\ &{}=\displaystyle \int _{0}^{2\sin (10(1-i)/n) } e^{-\eta ^{\alpha }t}\sum _{k=0} ^{+\infty }c^{n}_{k}\eta ^{k}d\eta .\\ \end{array} \end{aligned}$$

Here $c_k^n$ stand for the Taylor expansion coefficients of $\displaystyle \frac{e^{-2ni\arcsin \eta /2}}{\sqrt{1-(\eta /2)^2}}$, which can be computed through symbolic computation.

Accordingly, we can accurately compute the kernel function from

$$\begin{aligned} \begin{aligned}&g_n(t,\alpha )\approx \displaystyle \frac{1}{\alpha \pi } Re\left[ \sum _{k=0}^{+\infty } c_{k}^{n}\frac{1}{t^{(k+1)/\alpha }}\Big (\int _{0}^{(2\sin 10(1-i)/n)^{\alpha }t} e^{-s} s^{(k+1)/\alpha -1}ds\Big ) \right] \\&=\displaystyle \frac{1}{\alpha \pi } Re\left[ \sum _{k=0}^{+\infty }c_{k}^{n}\frac{1}{t^{(k+1)/\alpha }}\gamma _{1}\Big ( \Big (2\sin 10(1-i)/n\Big )^{\alpha }t, (k+1)/\alpha \Big )\right] .\\ \end{aligned} \end{aligned}$$

(53)

The function ${\dot{g}}_{0}(t,\alpha )$ is treated in the same way. We give the following proof of Lemma 2. $\square$

Proof

As a matter of fact, we consider $n=0$ in (20) and (21), and take time derivative. The (20) one gives

$$\begin{aligned} \begin{array}{ll} {\dot{g}}_{0}(t,\alpha )&{}=-\displaystyle \sum _{k=0}^{+\infty }\frac{(-t)^{k-1}k}{k!} \displaystyle \frac{\Gamma (1+k\alpha )}{\Gamma (1+(k\alpha +2n)/2)\Gamma (1+(k\alpha -2n)/2)}\\ &{}=-\displaystyle \sum _{k=0}^{+\infty }\frac{(-t)^{k}}{k!}\displaystyle \frac{\Gamma (1+k\alpha +\alpha )}{\Gamma (1+(k\alpha +\alpha +2n)/2)\Gamma (1+(k\alpha +\alpha -2n)/2)}, \end{array} \end{aligned}$$

(54)

which leads to (28). The (29) can be obtained in the same way. $\blacksquare$

Appendix B: Decoupled Governing System in Two Space Dimensions

Now we give the proof of Theorem 2.

Proof

To prove (42), we first show that $\{f_{m,n}(t,\vec {\kappa },\vec {\alpha })|m,n\in {\mathbb {Z}}\}$ form a fundamental solution to (41), i.e., it holds that

$$\begin{aligned} \begin{array}{lll} {\dot{f}}_{m,n}(t,\vec {\kappa },\vec {\alpha })&{}=\kappa _{1} \displaystyle \sum _{p=-\infty }^{+\infty }c_{p}(\alpha _{1}) f_{m-p,n} (t,\vec {\kappa },\vec {\alpha })+\kappa _{2}\sum _{q=-\infty }^{+\infty }c_{q}(\alpha _{2}) f_{m,n-q}(t,\vec {\kappa },\vec {\alpha })\\ &{}\quad +\displaystyle \frac{1}{4\pi ^2}\int _{0}^{2\pi }\int _{0}^{2\pi }e^{-im\xi -in\eta }\Big (\kappa _{1} (2\sin (\xi /2))^{\alpha _{1}}+\kappa _{2}(2\sin (\eta /2))^{\alpha _{2}}\Big )\\ &{} d\xi d\eta \delta (t), (m,n) \ne (0,0);\\ f_{m,n}(0,\vec {\kappa },\vec {\alpha })&{}=-\displaystyle \frac{1}{4\pi ^2}\int _{0}^{2\pi } \int _{0}^{2\pi }e^{-im\xi -in\eta }\Big (\kappa _{1}(2\sin (\xi /2))^{\alpha _{1}} +\kappa _{2}(2\sin (\eta /2))^{\alpha _{2}}\Big ) d\xi d\eta ;&{}\\ f_{0,0}(t,\vec {\kappa },\vec {\alpha })&{}=\delta (t). &{} \end{array} \end{aligned}$$

(55)

In fact, one can easily see for $(m,n)\ne (0,0)$,

$$\begin{aligned}&s{\tilde{f}}_{m,n}(s,\vec {\kappa },\vec {\alpha })-\kappa _{1}\displaystyle \sum _{p=-\infty }^{+\infty }c_{p}(\alpha _{1}) {\tilde{f}}_{m-p,n}(s,\vec {\kappa },\vec {\alpha })-\kappa _{2}\displaystyle \sum _{q=-\infty }^{+\infty }c_{q}(\alpha _{1}) {\tilde{f}}_{m,n-q}(s,\vec {\kappa },\vec {\alpha })\nonumber \\&\quad =\displaystyle \frac{1}{4\pi ^2{\tilde{g}}_{0,0}(s,\vec {\kappa },\vec {\alpha })} \int _{0}^{2\pi }\int _{0}^{2\pi }\nonumber \\&\quad \frac{se^{-im\xi -in\eta }- \kappa _{1}\displaystyle \sum _{p=-\infty }^{+\infty }c_{p}(\alpha _{1})e^{-i(m-p)\xi -in\eta } -\kappa _{2}\displaystyle \sum _{q=-\infty }^{+\infty }c_{q}(\alpha _{2})e^{-im\xi -i(n-q)\eta }}{s+\kappa _{1}(2\sin (\xi /2))^{\alpha _{1}} +\kappa _{2}(2\sin (\eta /2))^{\alpha _{2}}}d\xi d\eta \nonumber \\&\quad =\displaystyle \frac{1}{4\pi ^2{\tilde{g}}_{0,0}(s,\vec {\kappa },\vec {\alpha })}\int _{0}^{2\pi }\int _{0}^{2\pi } \frac{e^{-im\xi -in\eta }\left( s-\kappa _{1}\displaystyle \sum _{p=-\infty }^{+\infty }c_{p}(\alpha _{1})e^{ip\xi }-\kappa _{2} \displaystyle \sum _{q=-\infty }^{+\infty }c_{q}(\alpha _{2})e^{i{q}\eta }\right) }{s+\kappa _{1}(2\sin (\xi /2))^{\alpha _{1}} +\kappa _{2}(2\sin (\eta /2))^{\alpha _{2}}}d\xi d\eta \nonumber \\&\quad =\displaystyle \frac{1}{4\pi ^2{\tilde{g}}_{0,0}(s,\vec {\kappa },\vec {\alpha })}\int _{0}^{2\pi }\int _{0}^{2\pi } \frac{e^{-im\xi -in\eta }\Big (s+\kappa _{1}(2\sin (\xi /2))^{\alpha _{1}}+\kappa _{2}(2\sin (\eta /2))^{\alpha _{2}}\Big )}{s+\kappa _{1}(2\sin (\xi /2))^{\alpha _{1}}+\kappa _{2}(2\sin (\eta /2))^{\alpha _{2}}}d\xi d\eta \nonumber \\&\quad =\displaystyle \frac{1}{4\pi ^2{\tilde{g}}_{0,0}(s,\vec {\kappa },\vec {\alpha })}\int _{0}^{2\pi }\int _{0}^{2\pi }e^{-im\xi -in\eta }d\xi d\eta \nonumber \\&\quad =0. \end{aligned}$$

(56)

In the mean time, by definition we have ${\tilde{f}}_{0,0}(s,\vec {\kappa },\vec {\alpha })=1$, hence $f_{0,0}(t,\vec {\kappa },\vec {\alpha })=\delta (t)$.

Secondly, from (45)(46), direct calculations show that

$$\begin{aligned} \begin{array}{ll} &{}\kappa _{1}\displaystyle \sum _{p=-\infty }^{+\infty }c_{p}(\alpha _{1}){\tilde{g}}_{-p,0}(s,\vec {\kappa },\vec {\alpha })+\kappa _{2}\displaystyle \sum _{q=-\infty }^{+\infty }c_{q}(\alpha _{2}){\tilde{g}}_{0,-q}(s,\vec {\kappa },\vec {\alpha })\\ &{}\quad =\displaystyle \frac{1}{4\pi ^2}\int _{0}^{2\pi }\int _{0}^{2\pi }\displaystyle \frac{\kappa _{1}\displaystyle \sum _{p=-\infty }^{+\infty }c_{p}(\alpha _{1})e^{ip\xi }+\kappa _{2}\displaystyle \sum _{q=-\infty }^{+\infty }c_{q}(\alpha _{2})e^{iq\eta }}{s+\kappa _{1}(2\sin \xi /2)^{\alpha _{1}}+\kappa _{2}(2\sin \eta /2)^{\alpha _{2}}}d\xi d\eta \\ &{}\quad =\displaystyle \frac{1}{4\pi ^2}\int _{0}^{2\pi }\int _{0}^{2\pi }\displaystyle \frac{-\kappa _{1}(2\sin (\xi /2))^{\alpha _{1}}-\kappa _{2}(2\sin (\eta /2))^{\alpha _{2}}}{ s+\kappa _{1}(2\sin \xi /2)^{\alpha _{1}}+\kappa _{2}(2\sin \eta /2)^{\alpha _{2}}}d\xi d\eta \\ &{}\quad =s{\tilde{g}}_{0,0}(s,\vec {\kappa },\vec {\alpha })-1. \end{array} \end{aligned}$$

(57)

This relates the kernel functions defined in (43) by dividing ${\tilde{g}}_{0,0}(s,\vec {\kappa },\vec {\alpha })$ on the both side of (57),

$$\begin{aligned} {\tilde{F}}_{2}(s,\vec {\kappa },\vec {\alpha })=\kappa _{1} \displaystyle \sum _{p=-\infty }^{+\infty }c_{p}(\alpha _{1}){\tilde{f}}_{-p,0}(s,\vec {\kappa }, \vec {\alpha })+\kappa _{2}\displaystyle \sum _{q=-\infty }^{+\infty }c_{q}(\alpha _{2}) {\tilde{f}}_{0,-q}(s,\vec {\kappa },\vec {\alpha }). \end{aligned}$$

(58)

Hence from (56) and (57) we get a uniform expression for ${\tilde{f}}_{m,n}(s,\vec {\kappa },\vec {\alpha })$

$$\begin{aligned}&s{\tilde{f}}_{m,n}(s,{\tilde{\kappa }},{\tilde{\alpha }})=\kappa _{1} \displaystyle \sum _{p=-\infty }^{+\infty }c_{p}(\alpha _{1}){\tilde{f}}_{m-p,n} (s,\vec {\kappa },\vec {\alpha })\nonumber \\&\quad +\kappa _{2}\displaystyle \sum _{q=-\infty }^{+\infty }c_{q} (\alpha _{2}){\tilde{f}}_{m,n-q}(s,\vec {\kappa },\vec {\alpha }) +\delta _{0,0}^{m,n}(s-{\tilde{F}}_{2}(s,\vec {\kappa },\vec {\alpha })), \end{aligned}$$

(59)

with $\delta ^{m,n}_{0,0}$ the Kronecker delta.

Noticing ${\tilde{f}}_{0,0}(s,\vec {\kappa },\vec {\alpha })=1$, we are now ready to verify that the solution to (42), or equivalently,

$$\begin{aligned} s{\tilde{\phi }}_{m,n}(s)={\tilde{F}}_{2}(s,\vec {\kappa },\vec {\alpha }) {\tilde{\phi }}_{m,n}(s)+\sum _{ (p,q)}\phi _{p,q}(0){\tilde{f}}_{m-p,n-q} (s,{\tilde{\kappa }},{\tilde{\alpha }}) \end{aligned}$$

(60)

solves the original lattice system (41), namely,

$$\begin{aligned} s{\tilde{\phi }}_{m,n}(s)- \phi _{m,n}(0) =\kappa _{1}\sum _{p=-\infty }^{+\infty } c_p(\alpha _{1}) {\tilde{\phi }}_{m-p,n}(s)+\kappa _{2}\sum _{q=-\infty }^{+\infty } c_q(\alpha _{2}) {\tilde{\phi }}_{m,n-q}(s). \end{aligned}$$

(61)

As a matter of fact, we compute

$$\begin{aligned}&s{\tilde{\phi }}_{m,n}(s)- \phi _{m,n}(0) -\kappa _{1}\displaystyle \sum _{k=-\infty }^{+\infty } c_k(\alpha _{1}) {\tilde{\phi }}_{m-k,n}(s)-\kappa _{2}\displaystyle \sum _{l=-\infty }^{+\infty } c_l(\alpha _{2}) {\tilde{\phi }}_{m,n-l}(s) \\&\quad = \displaystyle \frac{s}{s-{\tilde{F}}_{2}(s,\vec {\kappa },\vec {\alpha })} \sum _{p,q} \phi _{p,q}(0){\tilde{f}}_{m-p,n-q}(s,{\tilde{\kappa }},{\tilde{\alpha }}) -\phi _{m,n}(0) \\&\qquad -\displaystyle \frac{\kappa _{1}}{s-{\tilde{F}}_{2}(s,\vec {\kappa },\vec {\alpha })} \displaystyle \sum _{k=-\infty }^{+\infty } c_k(\alpha _{1}) \displaystyle \sum _{p,q} \phi _{p,q}(0){\tilde{f}}_{m-k-p,n-q}(s,\vec {\kappa },\vec {\alpha }) \\&\qquad -\displaystyle \frac{\kappa _{2}}{s-{\tilde{F}}_{2}(s,\vec {\kappa },\vec {\alpha })} \displaystyle \sum _{l=-\infty }^{+\infty } c_l(\alpha _{2}) \displaystyle \sum _{p,q} \phi _{p,q}(0){\tilde{f}}_{m-p,n-l-q}(s,\vec {\kappa },\vec {\alpha }) \\&\quad = \displaystyle \frac{1}{s-{\tilde{F}}_{2}(s,\vec {\kappa },\vec {\alpha })} \sum _{p,q} \phi _{p,q}(0)\Big [s{\tilde{f}}_{m-p,n-q}(s,{\tilde{\kappa }},{\tilde{\alpha }})- \kappa _{1}\sum _{k=-\infty }^{+\infty } c_k(\alpha _{1}){\tilde{f}}_{m-k-p,n-q}(s,\vec {\kappa },\vec {\alpha })\\&\qquad - \kappa _{2}\displaystyle \sum _{l=-\infty }^{+\infty } c_l(\alpha _{2}){\tilde{f}}_{m-p,n-l-q}(s,\vec {\kappa },\vec {\alpha })\Big ] - \phi _{m,n}(0)\\&\quad = \displaystyle \frac{1}{s-{\tilde{F}}_{2}(s,\vec {\kappa },\vec {\alpha })} \sum _{p,q} \phi _{p,q}(0)\delta ^{m-p,n-q}_{0,0} (s- {\tilde{F}}_{2}(s,\vec {\kappa },\vec {\alpha }))- \phi _{m,n}(0)\\&\quad =0. \end{aligned}$$

$\blacksquare$