A collocation method of lines for two-sided space-fractional advection-diffusion equations with variable coefficients

TL;DR

This paper introduces a spectral collocation-based Method of Lines for solving space-fractional advection-diffusion equations with variable coefficients, demonstrating high accuracy and efficiency through numerical examples.

Contribution

It develops a new basis function set and applies a spectral collocation Method of Lines to effectively solve complex space-fractional PDEs with variable coefficients.

Findings

Achieves exponential accuracy in numerical solutions.

Successfully handles equations with discontinuous coefficients.

Demonstrates efficiency through benchmark problems.

Abstract

We present the Method Of Lines (MOL), which is based on the spectral collocation method, to solve space-fractional advection-diffusion equations (SFADEs) on a finite domain with variable coefficients. We focus on the cases in which the SFADEs consist of both left- and right-sided fractional derivatives. To do so, we begin by introducing a new set of basis functions with some interesting features. The MOL, together with the spectral collocation method based on the new basis functions, are successfully applied to the SFADEs. Finally, four numerical examples, including benchmark problems and a problem with discontinuous advection and diffusion coefficients, are provided to illustrate the efficiency and exponentially accuracy of the proposed method.

Figures

Tables

Table 1. Table 1: The l 2 subscript 𝑙 2 l_{2} and l ∞ subscript 𝑙 l_{\infty} errors for Example 2, obtained by the proposed method with n = 1 , 2 , 3 𝑛 1 2 3 n=1,2,3 and some values of α 𝛼 \alpha and β 𝛽 \beta with γ = 2 𝛾 2 \gamma=2 .

				$\beta =1.2$				$\beta =1.4$				$\beta =1.6$				$\beta =1.8$
$\alpha$		$n$		$E_2\left(n\right)$		$E_{\infty }\left(n\right)$		$E_2\left(n\right)$		$E_{\infty }\left(n\right)$		$E_2\left(n\right)$		$E_{\infty }\left(n\right)$		$E_2\left(n\right)$		$E_{\infty }\left(n\right)$
0.20		1		5.5e-03		4.6e-01		5.2e-03		4.6e-01		5.0e-03		4.6e-01		4.6e-03		4.5e-01
0.20		2		8.5e-14		1.4e-11		4.4e-14		7.9e-12		5.6e-14		1.1e-11		1.4e-14		3.0e-12
0.20		3		8.3e-14		1.4e-11		4.3e-14		7.8e-12		5.6e-14		1.1e-11		1.4e-14		3.0e-12
0.40		1		6.4e-03		5.6e-01		6.2e-03		5.6e-01		5.9e-03		5.6e-01		5.5e-03		5.4e-01
0.40		2		8.2e-14		1.4e-11		4.3e-14		8.0e-12		5.5e-14		1.1e-11		1.3e-14		3.0e-12
0.40		3		8.1e-14		1.4e-11		4.2e-14		7.8e-12		5.5e-14		1.1e-11		1.3e-14		3.0e-12
0.60		1		7.6e-03		6.8e-01		7.3e-03		6.9e-01		7.0e-03		6.8e-01		6.6e-03		6.7e-01
0.60		2		7.8e-14		1.3e-11		4.1e-14		7.6e-12		5.4e-14		1.1e-11		1.3e-14		3.0e-12
0.60		3		7.9e-14		1.4e-11		4.1e-14		7.9e-12		5.3e-14		1.1e-11		1.3e-14		2.9e-12
0.80		1		9.0e-03		8.4e-01		8.7e-03		8.4e-01		8.3e-03		8.4e-01		7.8e-03		8.2e-01
0.80		2		7.4e-14		1.3e-11		4.0e-14		7.7e-12		5.2e-14		1.1e-11		1.3e-14		2.9e-12
0.80		3		7.4e-14		1.3e-11		4.0e-14		7.7e-12		5.2e-14		1.1e-11		1.3e-14		2.9e-12

Table 2. Table 2: Different values of parameters α 𝛼 \alpha , β 𝛽 \beta , K α subscript 𝐾 𝛼 K_{\alpha} and K β subscript 𝐾 𝛽 K_{\beta} considered in Example 3 for problem with s ​ ( x , t ) = 0 𝑠 𝑥 𝑡 0 s(x,t)=0 (see Fig.  4 ).

	Case 1	Case 2	Case 3	Case 4	Case 5	Case 6	Case 7	Case 8	Case 9
$K_{\alpha }$	0.0	0.0	0.0	0.0	0.0	0.0	0.2	0.3	0.4
$\alpha$	-	-	-	-	-	-	0.1	0.5	0.9
$K_{\beta }$	0.1	0.2	0.3	0.1	0.2	0.3	0.01	0.01	0.01
$\beta$	2.0	2.0	2.0	1.5	1.5	1.5	2.0	2.0	2.0

Equations

$$\displaystyle\frac{{\partial u(x,t)}}{{\partial t}}+\Big{(}{c_{\alpha}^{+}}(x,t){{}_{0}^{R}D_{x}^{\alpha}}+{c_{\alpha}^{-}}(x,t){{}_{x}^{R}D_{\ell}^{\alpha}}\Big{)}u(x,t)=\left({c_{\beta}^{+}}(x,t){{}_{0}^{R}D_{x}^{\beta}}+{c_{\beta}^{-}}(x,t){{}_{x}^{R}D_{\ell}^{\beta}}\right)u(x,t)+s(x,t),$$

$$u(x,0)=f(x),$$

$$u(0,t)=u(\ell ,t)=0.$$

$${}_a^R{D}_x^{\sigma }g(x)$$

$${}_x^R{D}_b^{\sigma }g(x)$$

$$P_0^{(a,b)}(x):=1,P_1^{(a,b)}(x):=\frac{1}{2}[(a+b+2)x+a-b],$$

$$P_{j+1}^{(a,b)}(x):=(A_j^{a,b}x-B_j^{a,b})P_j^{(a,b)}(x)-C_j^{a,b}{P}_{j-1}^{(a,b)}(x),j=1,2,\dots ,m,$$

$$A_j^{a,b}:=\frac{(a+b+2j+1)(a+b+2j+2)}{2(j+1)(a+b+j+1)},$$

$$B_j^{a,b}:=\frac{(b^2-a^2)(a+b+2j+1)}{2(j+1)(a+b+j+1)(a+b+2j)},$$

$$C_j^{a,b}:=\frac{(a+j)(b+j)(a+b+2j+2)}{(j+1)(a+b+j+1)(a+b+2j)}.$$

$${\int }_{-1}^1{P}_j^{(a,b)}(x)P_i^{(a,b)}(x)w^{(a,b)}(x)dx={\gamma }_j^{a,b}{\delta }_{ij},$$

$${\gamma }_j^{a,b}:=\frac{2^{a+b+1}\Gamma (j+a+1)\Gamma (j+b+1)}{(2j+a+b+1)j!\Gamma (j+a+b+1)}.$$

$$P_j^{(a,b)}(-x)=(-1{)}^j{P}_j^{(b,a)}(x),$$

$$\frac{d^m}{d{x}^m}{P}_j^{(a,b)}(x)=d_{j,m}^{a,b}{P}_{j-m}^{(a+m,b+m)}(x),j\ge m\in \mathbb{N},$$

$$d_{j,m}^{a,b}:=\frac{\Gamma (j+m+a+b+1)}{2^m\Gamma (j+a+b+2)}.$$

$$P_k^{(a,b-1)}\left(\frac{2x}{\ell }-1\right)=\frac{k+a+b}{2k+a+b}{P}_k^{(a,b)}\left(\frac{2x}{\ell }-1\right)+\frac{k+a}{2k+a+b}{P}_{k-1}^{(a,b)}\left(\frac{2x}{\ell }-1\right).$$

$${\int }_{-1}^1(1-x{)}^a(1+x{)}^{b}f(x)\text{d}x\simeq \sum _{i=0}^m{\omega }_i^{(a,b)}f({\xi }_i^{(a,b)}),$$

$$w_i^{(a,b)}:=\frac{2^{a+b+1}\Gamma (m+a+2)\Gamma (m+b+2)}{(m+1)!\Gamma (m+a+b+2)(1-x_i^2)[(P_{m+1}^{(a,b)}{)}^{\prime }(x_i)]}.$$

$${}_0^R{D}_x^{\alpha }\left[x^b{P}_k^{(a,b)}\left(\frac{2x}{\ell }-1\right)\right]=\frac{\Gamma (k+b+1)}{\Gamma (k+b-\alpha +1)}{x}^{b-\alpha }{P}_k^{(a+\alpha ,b-\alpha )}\left(\frac{2x}{\ell }-1\right).$$

$${\mathcal{J}}_k^{(a,b,\rho ,\theta )}(x):=x^{\rho }(\ell -x{)}^{\theta }{P}_k^{(a,b)}\left(\frac{2x}{\ell }-1\right),k=0,1,\dots .$$

$$\frac{d^{k}}{dx^{k}}\mathcal{J}^{(a,b,\rho,\theta)}_{k}(x)\Big{|}_{x=0}=\frac{d^{k}}{dx^{k}}\mathcal{J}^{(a,b,\rho,\theta)}_{k}(x)\Big{|}_{x=\ell}=0.$$

$${\mathcal{J}}_0^{(a,b,\rho ,\theta )}(x)=x^{\rho }(\ell -x{)}^{\theta },{\mathcal{J}}_1^{(a,b,\rho ,\theta )}(x)=\frac{1}{2}\left[(a+b+2)\left(\frac{2x}{\ell }-1\right)+a-b\right]x^{\rho }(\ell -x{)}^{\theta },$$

$${\mathcal{J}}_{j+1}^{(a,b,\rho ,\theta )}(x)=\left(A_j^{a,b}\left(\frac{2x}{\ell }-1\right)-B_j^{a,b}\right){\mathcal{J}}_j^{(a,b,\rho ,\theta )}(x)-C_j^{a,b}{\mathcal{J}}_{j-1}^{(a,b,\rho ,\theta )}(x),j=1,2,\dots ,m,$$

$${\int }_0^{\ell }{\mathcal{J}}_i^{(a,b,\rho ,\theta )}(x){\mathcal{J}}_j^{(a,b,\rho ,\theta )}(x)w^{(a,b,\rho ,\theta )}(x)dx={\left(\frac{\ell }{2}\right)}^{a+b+1}{\gamma }_j^{a,b}{\delta }_{ij},$$

$$u(x,t)\simeq {\tilde{u}}_n(x,t)=\sum _{k=0}^n{a}_k(t){\phi }_k(x),$$

$${\phi }_k(x):={\lambda }_k{\mathcal{J}}_k^{(1,1,1,1)}(x)={\lambda }_{k}x(\ell -x)P_k^{(1,1)}\left(\frac{2x}{\ell }-1\right)$$

$${\lambda }_k:=\frac{1}{{\ell }^3}\frac{\left(k+2\right)\left(2k+3\right)}{k+1}.$$

$$\frac{\partial }{\partial t}u(x,t)\simeq \frac{\partial }{\partial t}{\tilde{u}}_n(x,t)=\sum _{k=0}^n{\dot{a}}_k(t){\phi }_k(x),$$

$${}_0^R{D}_x^{\alpha }u(x,t)\simeq {}_0^R{D}_x^{\alpha }{\tilde{u}}_n(x,t)=\sum _{k=0}^n{a}_k(t){\psi }_k^{\alpha +}(x),{}_0^R{D}_x^{\beta }u(x,t)\simeq {}_0^R{D}_x^{\beta }{\tilde{u}}_n(x,t)=\sum _{k=0}^n{a}_k(t){\psi }_k^{\beta +}(x),$$

$${}_x^R{D}_{\ell }^{\alpha }u(x,t)\simeq {}_x^R{D}_{\ell }^{\alpha }{\tilde{u}}_n(x,t)=\sum _{k=0}^n{a}_k(t){\psi }_k^{\alpha -}(x),{}_x^R{D}_{\ell }^{\beta }u(x,t)\simeq {}_x^R{D}_{\ell }^{\beta }{\tilde{u}}_n(x,t)=\sum _{k=0}^n{a}_k(t){\psi }_k^{\beta -}(x),$$

$${\psi }_k^{\alpha +}(x):={}_0^R{D}_x^{\alpha }{\phi }_k(x),{\psi }_k^{\beta +}(x):={}_0^R{D}_x^{\beta }{\phi }_k(x),$$

$${\psi }_k^{\alpha -}(x):={}_x^R{D}_{\ell }^{\alpha }{\phi }_k(x),{\psi }_k^{\beta -}(x):={}_x^R{D}_{\ell }^{\beta }{\phi }_k(x).$$

$$\begin{gathered} \sum _{k=0}^n{\dot{a}}_k(t){\phi }_k(x)+\left(c_{\alpha }^{+}(x,t)\sum _{k=0}^n{a}_k(t){\psi }_k^{\alpha +}(x)+c_{\alpha }^{-}(x,t)\sum _{k=0}^n(-1{)}^k{a}_k(t){\psi }_k^{\alpha -}(x)\right) \\ =\left(c_{\beta }^{+}(x,t)\sum _{k=0}^n{a}_k(t){\psi }_k^{\beta +}(x)+c_{\beta }^{-}(x,t)\sum _{k=0}^n(-1{)}^k{a}_k(t){\psi }_k^{\beta -}(x)\right)+s(x,t). \end{gathered}$$

$$\begin{gathered} \sum _{k=0}^n{\dot{a}}_k(t){\phi }_k({\hat{\xi }}_i^{(1,1)})+\left(c_{\alpha }^{+}({\hat{\xi }}_i^{(1,1)},t)\sum _{k=0}^n{a}_k(t){\psi }_k^{\alpha +}({\hat{\xi }}_i^{(1,1)})+c_{\alpha }^{-}({\hat{\xi }}_i^{(1,1)},t)\sum _{k=0}^n(-1{)}^k{a}_k(t){\psi }_k^{\alpha -}({\hat{\xi }}_i^{(1,1)})\right) \\ =\left(c_{\beta }^{+}({\hat{\xi }}_i^{(1,1)},t)\sum _{k=0}^n{a}_k(t){\psi }_k^{\beta +}({\hat{\xi }}_i^{(1,1)})+c_{\beta }^{-}({\hat{\xi }}_i^{(1,1)},t)\sum _{k=0}^n(-1{)}^k{a}_k(t){\psi }_k^{\beta -}({\hat{\xi }}_i^{(1,1)})\right)+s({\hat{\xi }}_i^{(1,1)},t). \end{gathered}$$

$$\mathbf{M}\,\dot{\mathbf{a}}(t)+\mathbf{C}^{+}_{\alpha}(t)\big{[}{\mathfrak{D}}^{+}_{\alpha}\ \mathbf{a}(t)\big{]}+\mathbf{C}^{-}_{\alpha}(t)\big{[}{\mathfrak{D}}^{-}_{\alpha}\ \mathbf{a}(t)\big{]}=\mathbf{C}^{+}_{\beta}(t)\big{[}{\mathfrak{D}}^{+}_{\beta}\ \mathbf{a}(t)\big{]}+\mathbf{C}^{-}_{\beta}(t)\big{[}{\mathfrak{D}}^{-}_{\beta}\ \mathbf{a}(t)\big{]}+\mathbf{s}(t),$$

$$\mathbf{M}:=\left[\begin{array}{cccc}{\phi }_0({\hat{\xi }}_0^{(1,1)})& {\phi }_0({\hat{\xi }}_1^{(1,1)})& \dots & {\phi }_0({\hat{\xi }}_n^{(1,1)})\\ {\phi }_1({\hat{\xi }}_0^{(1,1)})& {\phi }_1({\hat{\xi }}_1^{(1,1)})& \dots & {\phi }_1({\hat{\xi }}_n^{(1,1)})\\ ⋮& ⋮& \ddots & ⋮\\ {\phi }_n({\hat{\xi }}_0^{(1,1)})& {\phi }_n({\hat{\xi }}_1^{(1,1)})& \dots & {\phi }_n({\hat{\xi }}_n^{(1,1)})\end{array}\right],\mathbf{a}(t):=\left[\begin{array}{c}{a}_0(t)\\ a_1(t)\\ ⋮\\ a_n(t)\end{array}\right],\mathbf{s}(t):=\left[\begin{array}{c}s({\hat{\xi }}_0^{(1,1)},t)\\ s({\hat{\xi }}_1^{(1,1)},t)\\ ⋮\\ s({\hat{\xi }}_n^{(1,1)},t)\end{array}\right],$$

Full text

\jnlcitation\cname

, , , and , (\cyear2019), \ctitleA collocation method of lines for two-sided space-fractional advection-diffusion equations with variable coefficients, \cjournalMathematical Methods in the Applied Sciences.

††articletype: Research Article

1 Introduction

Fractional Calculus (FC) started with some speculations of Leibniz and Euler, and it has been developed progressively up until now. The initial fundamental works on the subject of FC are gathered in the books of Miller and Ross 1, Kilbas et al2, Hilfer 3, Podlubny 4 and Kiryakova 5. Without any doubt, in today’s world, the concept of FC arises in many application’s fields. In fact, the recent developments of FC have shown that a large number of complex systems can be accurately modelled by ordinary and partial differential equations with fractional derivatives. These facts have created a new revolution of development in various fields of physics and engineering. 6, 7, 8

In the last decade, fractional Advection-Diffusion Equations (ADEs) have received considerable attention in physics, mathematics and applications. A fractional ADE contains derivatives of fractional order in time, space or time-space. These three types of fractional ADEs have been used increasingly in the modelling of physical phenomena. Fractional ADEs are more proper for describing some phenomena such as anomalous diffusions, which arise in complex dynamics 9, transport of passive tracers carried by fluid flow in porous mediums 10, 11, velocity of particle motion in blood flow or underground water 12, etc. Various methods have been developed for solving fractional ADEs. We refer to some of them. Zheng et al13 used the finite element method to solve a space-fractional ADE with non-homogeneous initial-boundary conditions. Yang and Liu14 considered a numerical method for ADEs with Riesz space fractional derivatives on a bounded domain. Ding and Zhang15 presented a numerical technique for solving ADEs with Riesz space fractional derivatives. A homotopy analysis method were used to solve ADM with Riesz space fractional derivatives by Ray and Sahoo16. The Crank–Nicolson–Galerkin and backward Euler approximate schemes for two-dimensional space-fractional ADEs were investigated by Yanmin et al17. To solve ADE with two-sided space-time fractional derivative, a high-accuracy Legendre tau method is developed and a convergence analysis provided by Bhrawy and Zaky 18. A meshless method for solving time-fractional ADEs is considered by Tayebi et al19. Zhu et al20 presented a differential quadrature method for solving time- and space- fractional ADEs. In the very recent paper 21, Li et al study fractional convection operators from physical and mathematical points of view and present a numerical scheme for space-fractional ADEs. As other numerical methods for fractional ADEs, we can refer to Zhai et al22, Bhrawy et al23 and Li and Deng24. For some numerical methods for other equations, which are relevant to fractional ADEs, we can refer the interested readers to Meerschaert and Tadjeran25, Dehghan and Abbaszadeh26, Feng et al27, 28, Esmaeili29, Zaky30 and references therein.

Here, the following two-sided space-fractional ADE with variable coefficients is considered:

[TABLE]

In Eq. (1a), the advection coefficients ${c_{\alpha}^{+}}(x,t)\geq 0$ and ${c_{\alpha}^{-}}(x,t)\geq 0$, and the diffusion coefficients ${c_{\beta}^{+}}(x,t)\geq 0$ and ${c_{\beta}^{-}}(x,t)\geq 0$, are some given functions, and ${{}_{0}^{R}D_{x}^{\sigma}}$ and ${{}_{0}^{R}D_{x}^{\sigma}}$ are the left- and right-sided Riemann–Liouville fractional derivatives of order $\sigma$, respectively. We remind that

[TABLE]

where $m$ is the integer number such that $m-1\leq\sigma<m$. It should be noted that the problem is considered on the domain $0<x<\ell$ and $0\leq t\leq T$. We consider the case $0<\alpha\leq 1$ and $1<\beta\leq 2$, where the parameters $\alpha$ and $\beta$ are the fractional orders (fractors) of the spatial variable. Moreover, the known function $s(x,t)$ is a source term.

We use the Method Of Lines (MOL) to solve problem (1). The MOL is one of the best methods for solving time-dependent partial differential equations numerically. It proceeds in two steps. In the first step, the spatial variable(s) are discretized by using one from various possible techniques, such as spectral, finite element, finite difference or meshless methods. In fact, the main idea of this method is to convert the partial differential equation to a system of initial values problems (IVPs). In the second step, a time integration method, such as a Runge–Kutta method, is utilized to solve the resulting system of IVPs. The MOL has been successfully applied to solve various types of problems in engineering and science.31, 32, 33, 34, 35

In this paper, a suitable set of basis functions for spatial discretization is constructed by using Jacobi polynomials. Moreover, the right and left space-fractional derivatives of the basis functions are derived in closed form, which help us to simplify the spatial discretization of the problem. Then, the unknown solution is expanded in terms of the basis functions, and by using a collocation technique, the problem is reduced to an explicit system of IVPs. By solving this system of IVPs, an approximation for the solution of the problem is obtained.

In Section 2, we briefly review some of the standard definitions of Jacobi polynomials. The main aim of this paper is presented in Section 3, where we construct and develop our MOL for solving problem (1). In Section 4, some numerical examples are provided to show the efficiency of the proposed method. Some concluding remarks are provided at the end of the paper, in Section 5 of conclusion.

2 Preliminaries on Jacobi polynomials

One of the most practical classes of orthogonal polynomials on a closed interval $[-1,1]$ is the Jacobi polynomials, denoted by $P_{j}^{(a,b)}(x)$, $j=0,1,\ldots$. It is worthy to note that the Chebyshev, Legendre and Gegenbauer polynomials, are some important special cases of the Jacobi polynomials. The Jacobi polynomials have great flexibility to develop efficient numerical methods for solving a wide range of fractional differential models. Consequently, in the last decade, the Jacobi polynomials have been widely used to solve fractional problems. 36, 37, 38, 39, 40, 41, 42, 43 Our method uses the Jacobi polynomials too. Thus, in this section we briefly review the Jacobi polynomials, Jacobi quadrature rules, and a relevant theorem on the fractional derivatives of Jacobi polynomials.

An efficient formula for evaluating the Jacobi polynomials $P_{j}^{(a,b)}(x)$, $a,b>-1$, $x\in[-1,1]$, is given by the following three-term recurrence equation:

[TABLE]

where

[TABLE]

The Jacobi polynomials, with respect to the weight function $w^{(a,b)}(x):=(1-x)^{a}(1+x)^{b}$, are orthogonal:

[TABLE]

where

[TABLE]

We remind the following two properties of the Jacobi polynomials, which are useful for our purposes:

[TABLE]

where

[TABLE]

The Jacobi polynomials satisfy the following important relation:

[TABLE]

Since the Jacobi polynomials are orthogonal in $[-1,1]$, all zeros of $P_{m+1}^{(a,b)}(\tau)$ are simple and belong to the interval $(-1,1)$ 44. These zeros are called the Jacobi–Gauss nodes with parameters $a$ and $b$, which we denote by $\{\xi^{(a,b)}_{i}\}_{i=0}^{m}$. It is worth mentioning that there is no closed-form formula to obtain the Jacobi–Gauss nodes. However, We can use stable and accurate methods to determine them 44. The Jacobi–Gauss quadrature rule with parameters $a$ and $b$ is based on the Jacobi–Gauss nodes $\{\xi^{(a,b)}_{i}\}_{i=0}^{m}$ and can be used for approximating the integral of a function over the interval $[-1,1]$ with weight $(1-x)^{a}(1+x)^{b}$ as

[TABLE]

where $\omega^{(a,b)}_{i}$, $i=0,\ldots,m$, are the Jacobi–Gauss quadrature weights, which are computed by

[TABLE]

The $(n+1)$-points Jacobi–Gauss quadrature rule is exact whenever $f(x)$ is a polynomial of degree less or equal than $2n+1$, i.e., in the Jacobi–Gauss quadrature rule, the degree of exactness is $2n+1$.

The following theorem plays a key role in establishing our method.

Theorem 2.1** (Zayernouri and Karniadakis45).**

Let $-1<b-\alpha<b$ and $b>0$. Then, for $x\in[0,\ell]$, we have

[TABLE]

3 The proposed method of lines

The first step of our method consists to introduce appropriate basis polynomials (or functions) and expand the unknown solution in terms of them.

3.1 Modified Jacobi basis functions

A good choice of the basis functions is crucial from the numerical point of view. To do so, we begin by presenting a new set of basis functions and then some useful properties for them are given.

Definition 3.1**.**

Let $a,b>-1$. The Modified Jacobi Functions (MJFs) on $[0,\ell]$ are denoted by $\mathcal{J}^{(a,b,\rho,\theta)}_{k}(x)$, $k=0,1,\ldots$, and are defined as

[TABLE]

Remark 3.2**.**

It is worth to point out that our new class of MJFs is a generalization of the following four classes of well-known basis functions:

• •

if $\rho=\theta=0$ in Definition 3.1, then the classical Jacobi polynomials 46 on $[0,\ell]$ are obtained;

• •

MJFs for $\theta=0$, $\rho=b=\alpha$, and $a=\alpha-q>-1$ are reduced to the basis functions presented in Esmaeili and Shamsi36;

• •

MJFs for $\rho=b$, $\theta=0$, and $\theta=a$, $\rho=0$, are reduced to the Jacobi-polyfractonomials (or generalized Jacobi functions) 45, 47;

• •

MJFs for $\theta=a$ and $\rho=b$ are reduced to the basis functions defined in Mao and Karniadakis48.

In the following, we give some simple but important properties of MJFs.

{property}

[Integer derivative]

If $k$ is a positive integer number and $\rho,\theta>k$, $k=0,1,\ldots$, then

[TABLE]

For stable and accurate evaluation of MJFs, the following three-term recurrence relation is advised.

{property}

[Three-term recurrence relation]

The MJFs satisfy the following recurrence relation:

[TABLE]

where $A_{j}^{a,b}$, $B_{j}^{a,b}$ and $C_{j}^{a,b}$ are defined in (2).

It is easy to verify that the MJFs are orthogonal with respect to the weight function.

{property}

[Orthogonality] The MJFs are orthogonal with respect to the weight function $w^{(a,b,\rho,\theta)}(x)$ as follows:

[TABLE]

where $w^{(a,b,\rho,\theta)}(x):=x^{b-2\rho}(\ell-x)^{a-2\theta}$ and $\gamma_{j}^{a,b}$ is defined in (4).

3.2 Spectral collocation method for spatial discretization

Now, we are in the position to discretize the spatial variable. For this purpose, we start to approximate the unknown function $u(x,t)$ in the problem (1) by $\tilde{u}_{n}(x,t)$ as follows:

[TABLE]

where

[TABLE]

and

[TABLE]

It is interesting to point out that, in our collocation method, we choose $\varphi_{k}(x)$, $k=0,\ldots,n$, as the basis functions. We should mention that our motivation to use these basis functions is that they have some good features. First, they satisfy the boundary conditions (1c), thus $\tilde{u}_{n}(x,t)$ automatically satisfies the boundary conditions (1c). Second, as we will show in the next Subsection 3.3, the left and also right fractional derivatives of these basis functions can be obtained in a closed form, which simplifies the discretization stage.

Remark 3.3**.**

We note that, in order to improve the efficiency of the method, it is desirable that MJFs with $\rho,\,\theta\in(0,1)$ are considered as the basis function. However, by this consideration, we cannot derive both left and right derivatives in a closed form based on Jacobi polynomials. Accordingly, we develop our method with $\rho=\theta=1$.

Now, from (8), we can obtain the following approximations for $\tfrac{\partial}{\partial t}u$, ${}_{0}^{R}D_{x}^{\alpha}u$ and ${}_{0}^{R}D_{x}^{\beta}u$:

[TABLE]

where

[TABLE]

Note that the computation of $\psi_{k}^{\alpha+}$, $\psi_{k}^{\beta+}$, $\psi_{k}^{\alpha-}$ and $\psi_{k}^{\beta-}$ are not straightforward. Thus, in the next Subsection 3.3, we present a novel method for evaluating them.

By substituting the approximations (11), (12) and (13) into Eq. (1a), we get

[TABLE]

Let $\hat{\xi}_{i}^{(1,1)}$, $i=0,\ldots,n$, be the Jacobi–Gauss nodes associated with the interval $[0,\ell]$, which are defined as $\hat{\xi}_{i}^{(1,1)}:=\frac{\ell}{2}(\xi_{i}^{(1,1)}+1)$. By collocating (16) at $x=\hat{\xi}_{i}^{(1,1)}$, $i=0,\dots,n$, we arrive at

[TABLE]

The above equations can be expressed in the following matrix form:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

By collocating the initial condition (1b) at $x=\hat{\xi}_{i}^{(1,1)}$, $i=0,\ldots,n$, we obtain that

[TABLE]

The above equations can be expressed in the following matrix form:

[TABLE]

where $\mathbf{F}:=[f(\hat{\xi}_{0}^{(1,1)}),f(\hat{\xi}_{1}^{(1,1)}),\ldots,f(\hat{\xi}_{n}^{(1,1)})]^{T}$ and $\mathbf{a}(0):=[a_{0}(0),a_{1}(0),a_{2}(0),\ldots,a_{n}(0)]^{T}$.

In summary, the main problem (1) is reduced to the following system of ordinary differential equations with initial condition:

[TABLE]

If we compute the inverse of the mass matrix $\mathbf{M}$, then we can derive an explicit form of the above IVP. In the next theorem, we present the explicit form of the inverse of the mass matrix $\mathbf{M}$.

Theorem 3.4**.**

The mass matrix $\mathbf{M}$ is nonsingular and the $(i,j)$th entry of it can be obtained explicitly as

[TABLE]

where $w_{i}^{(1,1)},\ i=0,1,2,\ldots,n$, are the Gauss–Jacobi weights.

Proof 3.5**.**

To prove the theorem, it is sufficient to show that the inner product of the $i$th raw of $\mathbf{M}$ and the $j$th column of $\mathbf{M}^{-1}$ is equal to $\delta_{ij}$, i.e.,

[TABLE]

To prove the above equations, we begin by using the orthogonality property (3), which allows us to write that

[TABLE]

By noting that $P_{i}^{(1,1)}(x)P_{j}^{(1,1)}(x)$ is a polynomial of degree at most $2n$, we conclude that the integral in the above equation can be computed exactly by using the Jacobi–Gauss quadrature rule based on the nodes $\{\xi_{j}\}_{j=0}^{n}$. In this way, from the above equation, we get

[TABLE]

The above equations can be reformulated as

[TABLE]

Using Eqs. (9)–(10) and by noting that $\xi_{k}^{(1,1)}=\frac{2\hat{\xi}_{k}^{(1,1)}}{\ell}-1$ and $\hat{\xi}_{k}^{(1,1)}(\ell-\hat{\xi}_{k}^{(1,1)})=\tfrac{\ell^{2}}{2^{2}}(1+\xi_{k}^{(1,1)})(1-\xi_{k}^{(1,1)})$, we arrive at Eq. (18) and thus the proof of the theorem is completed.

By using $\mathbf{M}^{-1}$, we can convert the implicit IVP (17) into the following explicit IVP:

[TABLE]

By solving the above system of IVPs, with a well-developed time-marching method 49, the unknown vector function $\mathbf{a}(t)$ is obtained numerically. Then, by using (8), an approximation of $u(x,t)$ is obtained.

3.3 On the computation of the left and right fractional derivatives of the basis functions

As we have seen, the fractional differentiation matrices ${{\mathbf{{\mathfrak{D}}}}}_{+}^{\alpha+}$, ${{\mathbf{{\mathfrak{D}}}}}_{+}^{\alpha-}$, ${{\mathbf{{\mathfrak{D}}}}}_{+}^{\beta+}$, and ${{\mathbf{{\mathfrak{D}}}}}_{+}^{\beta-}$ have a key role in our method. These matrices simplify the spatial discretization process by replacing the left and right fractional differentiations with matrix-vector products. To derive these matrices, we can use Eqs. (14) and (15). However, before it, we need to obtain the functions $\psi_{k}^{\alpha+}$, $\psi_{k}^{\alpha-}$, $\psi_{k}^{\beta+}$ and $\psi_{k}^{\beta-}$. At first glance, it seems that these functions can be easily obtained, by taking the analytical left/right fractional derivative of $\varphi_{k}(x)$. However, taking the analytical fractional derivative, especially for large $n$, is sophisticated. In this respect, in what follows, we present some theorems, which play key roles in obtaining the mentioned functions.

Theorem 3.6**.**

Let $\varphi_{k}(x)$, $k=0,1,\ldots$, be defined in (9) and $\psi_{k}^{\sigma+}(x):={}_{0}^{R}D_{x}^{\sigma}\varphi_{k}(x)$. Then, for the cases $0<\sigma<1$ or $1<\sigma<2$, we have

[TABLE]

and, for $k\geq 1$,

[TABLE]

When $\sigma=2$, we get

[TABLE]

[TABLE]

where $\hat{P}^{(a,b)}_{k}\left(x\right):=\displaystyle P^{(a,b)}_{k}\left(\frac{2x}{\ell}-1\right)$.

Proof 3.7**.**

The evaluation of $\psi_{0}(x)$ is trivial. For $k\geq 1$, we have

[TABLE]

and thus

[TABLE]

Using Theorem 2.1, one has

[TABLE]

Moreover, using the property (7), with $a=1$ and $b=2$, we get

[TABLE]

Multiplying both sides of (22) by $\lambda_{k}x^{2}$, and applying the left-side fractional derivative, we arrive at

[TABLE]

Using Theorem 2.1 once again, we conclude that

[TABLE]

The proof of this theorem is concluded by plugging Eqs. (21) and (23) into relation (20). When $\sigma=2$, the use of property (6) completes the proof.

The next theorem is an important result, which has a key role for us to obtain a closed form for the right Riemann–Liouville fractional derivatives of the considered basis functions.

Theorem 3.8**.**

Let $\sigma>0$ and $m-1\leq\sigma<m$. Then

[TABLE]

Proof 3.9**.**

We have

[TABLE]

By applying the change of variable $t\rightarrow\ell-\tau$, we get

[TABLE]

and, by noting that $\varphi_{k}(\ell-\tau)=\varphi_{k}(\tau)$, we conclude that

[TABLE]

The proof is complete.

Theorem 3.8 states that to evaluate the right-sided fractional derivative of order $\sigma$ of $\varphi_{k}(x)$, we can use the left-sided fractional derivative of order $\sigma$ of $\varphi_{k}(x)$, then substitute $x$ by $\ell-x$ and, finally, multiply it by $(-1)^{m}$. This property is related with the concept of duality as introduced by Caputo and Torres 50.

In summary, by the explicit formulas in Theorems 3.6 and 3.8, for functions $\psi_{k}^{\sigma+}$ and $\psi_{k}^{\sigma-}$, we can simply compute the elements of the matrices ${{\mathbf{{\mathfrak{D}}}}}_{+}^{\alpha+}$, ${{\mathbf{{\mathfrak{D}}}}}_{+}^{\alpha-}$, ${{\mathbf{{\mathfrak{D}}}}}_{+}^{\beta+}$, and ${{\mathbf{{\mathfrak{D}}}}}_{+}^{\beta-}$.

4 Numerical experiments

In this section we illustrate the proposed method of lines by using concrete numerical experiments. We have implemented our method using Matlab on a 3.5 GHz Core i7 PC with 8 GB of RAM. Moreover, to solve IVP (19), the solver ode45 was used. This solver is based on an explicit Dormand–Prince Runge–Kutta (4, 5) formula 51, 52. By using a suitable error estimation, the solver generates an adaptive mesh $t_{0},\dots,t_{k}$ and obtains the values of the solution of IVP at this mesh. This formula is a single-step solver, i.e., in computing $\mathbf{a}(t_{k})$, it needs only the solution at the immediately preceding time point, $\mathbf{a}(t_{k-1})$. In ode45, we can adjust the accuracy of the solution through the input parameters AbsTol and RelTol. In our numerical experiments, we set AbsTol=$10^{-14}$ and RelTol=$10^{-12}$.

We consider four different examples. For the problems with a known exact solution, we assess the accuracy of our method by reporting the following $l_{2}$ and $l_{\infty}$ errors:

[TABLE]

Example 1

For the first example, we consider problem (1) with the following data borrowed from Meerschaert and Tadjeran53:

[TABLE]

In this problem, ehe exact solution is

[TABLE]

By applying the proposed MOL with $n=4$, the obtained approximated solution and the error function are plotted in Fig. 1. Moreover, the errors $E_{2}(n)$ and $E_{\infty}(n)$ for various values of $n$, together with their CPU times, are shown in Fig. 2. As plotted in Fig. 2 (left), the errors $E_{2}(n)$ and $E_{\infty}(n)$ decrease very rapidly until, in $n=4$, they take over about $1e-13$. It is noted that, for $n>4$, the round-off errors on computer prevent any further improvement. However, it is worthwhile to note that, for $n>4$, the errors remain about $1e-13$, which shows that the presented MOL is numerically stable.

Example 2

For the second example, we consider problem (1) with the following data:

[TABLE]

The exact solution for this problem is

[TABLE]

This example is taken from Li et al21. In Table 1, we report the obtained results by applying the proposed method with $\gamma=2$ and some values of $\alpha$ and $\beta$. As it is seen, accurate results are obtained even by using a small $n$.

Example 3

In this example, we consider problem (1) with the data

[TABLE]

This example, with $s(x,t):=0$ and different values of parameters $K_{\alpha}$ and $K_{\beta}$, is treated in Le et al21 and Yang et al14. When $s(x,t)=0$, the exact solution to this problem is, however, not known. For $\alpha=0.5$, $\beta=1.5$ and

[TABLE]

the exact solution of this problem is given by

[TABLE]

We point out that FresnelC and FresnelS are the Fresnel cosine and Fresnel sine integral functions, respectively, which are defined as

[TABLE]

In contrast with the last two examples, here the exact solution is not a polynomial type function and thus this example is more proper for assessing the accuracy and convergence rate of our method. With this purpose, the errors $E_{2}(n)$ and $E_{\infty}(n)$ for various values of $n$ with $K_{\alpha}=2$ and $K_{\beta}=0.1$ are plotted in Fig. 3 (left). It can be observed from this figure that the errors decrease rapidly and the spectral convergence rate is gained by the method. Moreover, the obtained solution with $n=30$ and its corresponding error are plotted in Fig. 3 (right).

In addition, in order to show the effects of parameters $\alpha$, $\beta$, $K_{\alpha}$ and $K_{\beta}$, the obtained solutions for Example 3 with $s(x,t)=0$ are plotted in Fig. 4 for the parameter values of Table 2.

Example 4

As a last example, we consider problem (1) with discontinuous data. Let

[TABLE]

and consider the diffusion coefficients in the following two cases:

[TABLE]

By applying the proposed method with $n=200$ to the problem in Case I and Case II, the solutions are obtained after 3.1 and 3.2 seconds, respectively. These solutions are plotted in Fig. 5. As we can see, the obtained solutions are confirmed from a physical point of view. We conclude that the presented method can be applied with success to challenging problems with discontinuous data.

5 Conclusion

We proposed an efficient and accurate numerical method for solving two-sided space-fractional advection-diffusion equations. We started by introducing a new set of basis functions, called the Modified Jacobi Functions (MJFs). The method of lines, together with the spectral collocation method based on the introduced basis functions, were successfully applied to the considered problems. In order to increase the efficiency of the proposed method, from the numerical point of view, the left- and right-sided Riemann–Liouville fractional differentiation matrices were derived. Several numerical examples were provided, showing good efficiency and high accuracy (exponential accuracy). The obtained results confirm that our method not only works well for problems with smooth solutions (see Examples 1, 2 and 3) but also can be easily applied for problems with discontinuities on the initial condition, advection or diffusion coefficients (see Example 4). Obtaining some theoretical estimates for the approximation errors and convergence analysis would be desirable. This work is currently in progress. A further interesting topic for research is using MJFs, with $\rho,\,\theta\in(0,1)$, as the basis functions to solve the considered problem, by a Galerkin scheme.

Acknowledgments

Torres was partially supported by FCT through the R&D Unit CIDMA, project UID/MAT/04106/2019. The authors are very grateful to anonymous referees for carefully reading their manuscript and for several comments and suggestions, which helped them to improve the paper.