Numerical approximation of a variable-order time fractional advection-reaction-diffusion model via shifted Gegenbauer polynomials

1.   Introduction

With the development of fractional derivative, operators without singular kernel function has become a research topic for the fact that it can better describe nonlocal dynamics systems. Some new classes of them are available including the Caputo-Fabrizio derivative ^[1], the Atangana-Baleanu derivative ^[2], etc. See ^{[3,4,5,6,7,8,9]} and references therein for further details. On the other hand, it is worth noting that the fractional calculus has extended to study the variable-order (VO) models ^[10,11], which are generalizations of fixed-order fractional derivatives. For example, some VO fractional differential diffusion processes can better simulate the temperature change than integer order classical models.

Many numerical treatments have recently been proposed to approximate VO fractional differential equations. VO fractional derivatives are global operators in essence, it is more convenient to use global technology such as spectral methods to deal with VO operators^[12]. Spectral methods can obtain the desired solutions with a small degree of freedom, which makes improved accuracy with a significant reduction in the computational cost^[13]. Here, we focus on the polynomial-based spectral collocation methods. Lagrange polynomials were applied to simulate two chaotic practice problems modeled by VO fractional differential equations ^[14]. A collocation technique based on a reproducing kernel function was developed for VO fractional initial value problems and terminal value problems ^[15]. Shifted Jacobi polynomials were used to approximate VO time fractional Emder-Fowler equation ^[16], two-dimensional VO time fractional cable equations^[17]. Gegenbauer polynomials are a special of Jacobi polynomials. As a significant kind of orthogonal polynomials, Gegenbauer polynomials are the generalization of Legendre polynomials and Chebyshev polynomials, which have been widely used in mathematical physics, engineering technology, scientific computing and other fields. Shifted Gegenbauer operational matrices jointly with the Tau method were developed for solving fractional differential equations ^[18]. Shifted Gegenbauer polynomials were applied to solve time-fractional delay differential equations ^[19]. Fractional optimal control problems were solved via Gegenbauer cardinal functions ^[20]. A numerical scheme in ^[21] was proposed to solve the nonlinear VO fractional Schr$\ddot{o}$dinger equation by using the shifted Legendre polynomials.

Diffusion process is one of the most famous processes in nature, such as the numerical simulation of urban air quality ^[22], transport process of pollutants in groundwater ^[23], numerical determination the pollution sources in a river ^[24], transport process of bimolecular reaction in porous media ^[25]. The fractional form of diffusion equation can model the random collision of solute molecules with fluid molecules causes diffusion and produces fluxes from high concentrated to low concentration regions. Numerical solutions of space fractional diffusion equation defined by Caputo derivative were obtained via shifted Gegenbauer polynomials ^[26]. A compact implicit difference scheme was proposed to approximate time-fractional diffusion-wave equation^[27]. An operational matrix was derived for multidimensional VO time fractional anomalous diffusion equations based on the shifted Chebyshev collocation methods^[28]. Gegenbauer spectral method was used to approximate convection-diffusion equation with time fractional derivative ^[29]. By using the shifted Chebyshev polynomial of the second kind as the cardinal function, an operational matrix scheme was given in ^[30] for approximating VO time fractional nonlinear reaction-diffusion equation. A homotopy perturbation method was applied for time fractional advection-reaction-diffusion equation with Liouville-Caputo derivative ^[31]. A numerical method was proposed to approximate time fractional advection-reaction-diffusion equation by using shifted Legendre polynomials ^[32].

Given the following VO time-fractional derivative model for advection-reaction-diffusion equation

Here, $\sigma$ and $\tau$ represent spatial and time nodes respectively, $u(\sigma, \tau)$ is concentration of the solute in fluid at finite distance, parameters $\kappa$ and $\mu$ are the constant coefficients on $\mathbb{R}$. Variable order $\theta(\sigma, \tau)\in C(0, 1)$, $\delta (\sigma, \tau, u(\sigma, \tau))$ and $f(\sigma, \tau)$ are continuous nonlinear source terms.

This paper aims to investigate numerically problem (1.1) with Atangana-Baleanu-Caputo defined time-fractional derivative without singular kernel ^[2]. For solving Eq (1.1), the Gegenbauer spectral method will be used. Firstly, the considered model is approximated by the shifted Gegenbauer polynomials (SGPs) with undetermined coefficients. A nonlinear algebraic equations are derived via the SGPs operational matrices by choosing suitable collocation points. Finally, the desired problem is transformed into solving the obtained nonlinear system.

This work proceeds as follows. Preliminaries on variable-order fractional derivatives with Mittage-Leffler kernel are provided in Section 2. Review of some useful properties of the SGPs is presented in Section 3. The operational matrices of cardinal functions for the SGPs are obtained in Section 4. Based on the Gegenbauer spectral method, the proposed numerical scheme is formulated in Section 5. A few illustrative examples with different initial conditions and boundary conditions are shown in Section 6. The last section is devoted to a brief conclusion.

2.   Variable order fractional calculus with Mittage-Leffler kernel

This section devotes to some useful preliminaries of Atangana-Baleanu-Caputo (ABC) derivatives with fractional order.

To begin with we introduce the Mittag-Leffler function ${\bf E}_{\alpha_1, \alpha_2}(\tau)$ that will commonly be encountered in fractional calculus as follows ^[33].

Notation $\boldsymbol{E}_{\alpha_1}(\tau)$ is used for $\alpha_2 = 1$. Where $\Gamma(\cdot)$ is the Gamma function and $\Gamma(x) = (x-1)!$ if $x\in\mathbb{Z}^{+}$. The relationship between the Gamma function and the Beta function is

The definition of the ABC fractional derivative for the order ${\bar \theta }$ is given by

where ${\bf C}(\bar \theta)$ is a normalization function taking the form

The definition of the ABC fractional derivative in the sense of variable order $\theta(\sigma, \tau)$ is given by

where ${{\bf C}}(\theta (\sigma, \tau)) = 1 - \theta (\sigma, \tau) + \frac{{\theta (\sigma, \tau)}}{{\Gamma (\theta (\sigma, \tau))}}$ is a normalization function, and $u(\sigma, \tau)\in C^{1}([0, X] \times [0, T])$ is a real function.

Corollary 2.1. (^[34]) Let $k \in \mathbb{N} \cup \{0\}$, there is

3.   The shifted Gegenbauer polynomials and their properties

Given a closed interval $[0, T]$, we use the shifted Gegenbauer polynomials (SGPs) of order $m$ as the cardinal function ^[35], which is denoted by $G_{T, m}^\lambda (\tau)$. It is defined by

with the orthogonal property

In Eq (3.2), the weight function is $w_T^\lambda (\tau) = {(T - \tau)^{\lambda - 0.5}}{\tau ^{\lambda - 0.5}}$ and

If $u(x)\in L^2[0, T]$, then it can be approximated by a polynomial $u_p$ of order $p$ in terms of the SGPs

with the coefficients $c_m$ are determined by the orthogonality condition

Let

then Eq (3.4) is expressed in matrix form as

Likewise, the function $u(\sigma, \tau) \in {L^2}([0, T]\times [0, T])$ with two variables can be approximated by the double SGPs of degrees $q$ and $p$ as

where $U = [u_{mn}]$ is a $(q + 1)\times (p + 1)$ unknown coefficients matrix with

for $m = 1, 2, \ldots, q+1, \, n = 1, 2, \ldots, p+1$.

4.   Operational matrices

In this section, some novel operator matrices about the SGPs are derived. Let

Lemma 4.1. The basis vector $\Psi_{T, p}(\tau)$ in (3.6) satisfies

where the transform matrix ${\bf \Theta}_{T, p} = [\xi_{T, mn}]$ is an $(p+1)$-order square matrix with

Proof. By expressing the element $\varphi_{p, \hat m}(\tau)(\hat m = 0, 1, \ldots, p)$ of ${\Phi_p}(\tau)$ in terms of the SGPs, we have

where

Eq (3.5) yields

By recalling that $w_T^\lambda (\tau) = {(T - \tau)^{\lambda - 0.5}}{\tau ^{\lambda - 0.5}}$ and using Eqs (3.1), (3.3) and (4.1), we get

where ${I_\tau } = \int_0^T {{{(T - \tau)}^{\lambda - 0.5}}{\tau ^{\lambda + {\hat m} + k - 0.5}}} d\tau$.

Let $\tau = Tx$ and from (2.2), it is

Substituting Eq (4.6) in (4.5) gives

Finally, the desired result can be obtained by changing the indices $m = \hat m + 1, n = \hat n + 1$ and replacing $\hat \Theta _{T, m-1}^T$ and ${{\hat \xi }_{T, (m - 1)(n - 1)}}$ by $\Theta _{T, m}^T$ and $\xi_{T, mn}$, respectively.

We give an example to illustrate for $\lambda = 1$, $T = 1$ and $m = 4$.

Lemma 4.2. The derivatives of the monomial function vector ${\Phi_p}(\tau)$ in (4.1) satisfy

and

in which $\mathbb{D}_p^{(r)}$ represents the $r^{\mathit{\text{th}}}$ power of the strictly lower triangular matrix $\mathbb{D}_p^{(1)}$ with

Proof. The result of the lemma can be implied by the definitions of vectors ${\Phi _p}(\tau)$ and $\mathbb{D}_p^{(1)}$.

Lemma 4.3. ${\Phi_p}(\tau)$ denotes the monomial function vector in (4.1), then

with

Proof. It is easy to complete the proof by using Corollary 2.1.

Theorem 4.1. The derivative of SGPs vector ${\Psi _{T, p}}(\tau)$ in (3.6) satisfies

where ${{\bf D}}_{T, p}^{(1)}$ represents the $(p+1)$-order operational matrix of the SGPs given by

Generally, the $r$-order derivative of ${\Psi _{T, p}}(\tau)$ satisfies

in which ${{\bf D}}_{T, p}^{(r)}$ represents the $r^{\mathit{\text{th}}}$ power of the matrix ${{\bf D}}_{T, p}^{(1)}$.

Proof. It can be proved simply from Eq (4.1) and Lemma 4.2.

Theorem 4.2. The SGPs vector ${\Psi _{T, p}}(\tau)$ in (3.6) satisfies

with

Here the $(p+1)$-order matrix ${{\bf Z}}_{T, p}^{(\theta (\sigma, \tau))}$ is the operator matrix of SGPs in the sense of variable order ${\theta (\sigma, \tau)}$.

Proof. The proof is completed by using Lemma 4.2 and Eq (4.10).

5.   Numerical scheme

The unknown solution $u(\sigma, \tau)$ of Eq (1.1) is approximated by the SGPs as follows

where ${\bf{U}} = \left[{{u_{mn}}} \right]_{(q + 1) \times (p + 1)}$ is an unknown matrix, ${\Psi _{X, q}}(\xi)$ and ${\Psi _{T, p}}(\tau)$ represent the vectors mentioned in Eq (3.8).

Using Eqs (4.11), (4.12) and (5.1), one has

and

Moreover, Theorem 4.2 results in

Combining Eqs (5.1)–(5.4) and (1.1), we have

From the initial and boundary conditions in Eq (1.1) and (5.1), we have

and

Finally, Eqs (5.5)–(5.7) are combined to the following system

For solving the unknown matrix ${\bf{U}}$ in the above $(q + 1) \times (p + 1)$ nonlinear algebraic equations, we can choose Gaussian nodes

and

Then by substituting ${\bf{U}}$ into Eq (5.1), an approximate solution of Eq (1.1) can be obtained.

6.   Numerical examples and results

The effectiveness of the proposed scheme is shown through some numerical simulations. All computations are carried out in MATLAB software. For $i = 1, 2$, let ${\varepsilon _i} = \max \left({\left| {u_E^i - u_N^i} \right|} \right)$ be the $i$-th maximum absolute error (MAE), with ${u_E^i}$ and ${u_N^i}$ are the analytical solution and numerical solution, respectively. ${\eta _i} = \left({{q_i} + 1} \right)\left({{p_i} + 1} \right)$ represents the number of the SGPs used in the $i$-th approximations. Then the order of the convergence (CO) corresponding to the approximate solutions is defined by

Example 1. If $\delta (\sigma, \tau, u(\sigma, \tau)) + f(\sigma, \tau) = 0$, then the model (1.1) reduces to ^[36]

If $\theta (\sigma, \tau) = 1$, the solution of Eq (6.1) reads

Set $\kappa = 0.1, \mu = 0.25, \lambda = 1$ and $q = 7, p = 7$. Figure 1 shows the numerical results for $u(\sigma, 0.4)$, $u(0.4, \tau)$ for different constant functions $\theta (\sigma, \tau)$. Numerical simulations for $u(\sigma, 0.6)$, $u(0.6, \tau)$ with some values of $\theta (\sigma, \tau)$ are displayed in Figure 2. Comparison of the approximation solution with $\theta (\sigma, \tau) = {\rm{0}}{\rm{.75 + 0}}{\rm{.15sin(2}}\pi \sigma \tau {\rm{)}}$ and the exact solution with $\theta (\sigma, \tau) = 1$ are depicted in Figure 3. It can be seen that the numerical solutions are convergent. Table 1 compares the absolute errors of the method proposed in this paper (abbreviated as SGPs method) with those obtained in ^[36] for various choices of $\theta (\sigma, \tau)$ and $q = p = 4$ at $t = 0.8$. It is observed that we achieved an excellent approximation for the exact solution. In Table 2, the $L_2$ errors are presented for $q = p = 6$ and various values $\theta (\sigma, \tau)$ at $x = 0.6$.

Example 2. If $\kappa = 1, \mu = 0$ and $f(\sigma, \tau) = 0$, then the model (1.1) becomes ^[30]

One finds that the solution of Eq (6.2) is $u(\sigma, \tau) = \left({1 + {{\rm{e}}^{\frac{{ - 1}}{{\sqrt 2 }}\left({\sigma + \frac{{(1 - 2\gamma)}}{{\sqrt 2 }}\tau } \right)}}} \right)^{ - 1}$ in the case $\theta (\sigma, \tau) = 1$. For $\lambda = 0.5, q = 9, p = 12$, the obtained numerical results of $u(\sigma, 0.2)$ and $u(3, \tau)$ with different constant functions $\theta (\sigma, \tau)$ are shown in Figure 4. Results of some functions $\theta (\sigma, \tau)$ with $\lambda = 1, q = 9, p = 9$ are shown in Figure 5. Comparison of the numerical solution with $\lambda = 1, q = 9, p = 9$, $\theta (\sigma, \tau) = 0.75+0.2\sin(\sigma \tau)$ and analytical solution with $\lambda = 1, q = 9, p = 9, \theta (\sigma, \tau) = 1$ are depicted in Figure 6. It can be seen that the numerical solutions are convergent.

In order to compare the results with the works in ^[30], we discuss the following case

The analytical solution of this case is given by $u(\sigma, \tau) = {\tau ^3}\sin (\sigma)$. This example has been solved for two different functions ${\theta (\sigma, \tau)}$ by utilizing the presented method (abbreviated as SGPs method) with $p = 3$ and some values of $q$. In Tables 3 and 4, the MAE errors are listed for the present method (abbreviated as SGPs method) and the method proposed in ^[30]. It can be seen that the SGPs method is able to achieve smaller errors.

Example 3. If $\kappa = \mu = 1$ and $f(\sigma, \tau){\rm{ = 0}}$, then the model (1.1) is

Let $\delta (\sigma, \tau, u(\sigma, \tau)) = u(\sigma, \tau)\frac{{{\partial ^2}u(\sigma, \tau)}}{{\partial {\sigma ^2}}} - u{(\sigma, \tau)^2} + u(\sigma, \tau)$ and $\theta (\sigma, \tau) = 1$, then the solution of Eq (6.4) reads $u(\sigma, \tau) = {e^{\sigma + \tau }}$. When $\lambda = 0.5, X = 0.5, T = 0.5, q = 10, p = 8$, numerical results of $u(\sigma, 0.2)$ and $u(0.2, \tau)$ with constant function $\theta (\sigma, \tau)$ are shown in Figure 7. Comparison of the numerical solution with $\lambda = 1, X = 0.5, T = 1, q = 9, p = 9$, $\theta (\sigma, \tau) = 0.75 + 0.15\sin (2\pi \sigma \tau)$ and the analytical solution with $\theta (\sigma, \tau) = 1$ are depicted in Figure 8. It can be seen that the numerical solutions are convergent.

Example 4. If $\kappa = \mu = 1$ and $\delta (\sigma, \tau, u(\sigma, \tau)) = 2u(\sigma, \tau)$, then the model (1.1) reads like

Given

the solution of Eq (6.5) is $u(\sigma, \tau) = {\rm{sin(}}\sigma{\rm{)}}{{\rm{e}}^{{\rm{(- \mathsf{ τ})}}}}$. Comparison of the numerical solution when $X = 0.5, T = 1, \lambda = 1, q = 6, p = 6, \theta (\sigma, \tau) = {\rm{0}}{\rm{.7 + 0}}{\rm{.2sin(}}\sigma \tau {\rm{)}}$ and the exact solution are depicted in Figure 9. Table 5 displays the results of the MAEs and CO for different parameters in Example 4. It is clear that the SGPs method has a good convergence rate.

Example 5. If $\kappa = \mu = 1$ and $\delta (\sigma, \tau, u(\sigma, \tau)) = 0$, then the model (1.1) is

For force function

the solution of Eq (6.6) reads $u(\sigma, \tau) = {\tau ^5}\sin \pi \sigma$. Comparison of the numerical solution when $\lambda = 1.2,$ $q = 14, \; p = 9$, $\theta (\sigma, \tau) = 0.8 - 0.2{e^{ - 2\sigma \tau }}$ and the exact solution are displayed in Figure 10. It can be seen that the numerical solution is convergent.

Finally, we try to use the proposed computational method to solve another type problem.

Example 6. Consider the VO fractional burgers equation^[37]

It is hard to find the analytical solution of Eq (6.7) with ABC fractional derivative, but we are able to predict by the proposed method simultaneously the solutions with ordinary and variable fractional derivative with some initial parameters. Choose $X = 1, \; T = 5$, numerical results of some spatial and time nodes with different $\theta (\sigma, \tau)$ are shown in Figures 11 and 12 respectively. Choose $\lambda = 0.5$, the 3D diagrams of the solutions predicted by the proposed algorithm when $\theta (\sigma, \tau) = 0.4$ and $\theta (\sigma, \tau) = {\text{0}}{\text{.75 + 0}}{\text{.2sin(2}}\sigma \tau)$ respectively are shown in Figure 13. It is clear that the numerical solution is convergent.

7.   Conclusions

Shifted orthogonal polynomials combined with spectral methods (such as Galerkin method, collocation method and tau method) can transform complex variable fractional differential equations into algebraic equations, so as to reduce the complexity and improve the accuracy. We propose a numerical approach to solve a VO time fractional model for advection-reaction-diffusion equation with Atangana-Baleanu-Caputo derivative via the shifted Gegenbauer polynomial. By using a collocation approach and obtaining the operational matrix, all that remains is to solve a nonlinear equations. Some examples are given to illustrate the effectiveness of the proposed numerical algorithm.

Acknowledgments

This work was supported by National Natural Science Foundation of China (No. 11971094).

Conflict of interest

The authors declare no conflicts of interest.