Investigation and analysis of the numerical approach to solve the multi-term time-fractional advection-diffusion model

1.   Introduction

Diffusion refers to the movement of substances in a physical system due to differences in the concentration of regions. The transfer of quantity due to the bulk movement of fluid is called advection ^[1,2]. In combination of the advection and diffusion equations, the transfer of substance or quantity (particle, energy or chemical substances) in a medium is explained, namely the advection-diffusion equation. The advection-diffusion equation is a partial differential equation used to describe physical phenomena, including two processes of advection and diffusion ^[3,4].

In recent years, researchers observed that the classical advection-diffusion equation is insufficient in explaining some physical phenomena. For example, in the classical equation of motion, the diffusion process is assumed to be linear, while there are many media in which the transport motion is not linear ^[5]. Therefore, the fractional advection-diffusion equation considers the non-linearity of the transport movement in the media and the feature of the memory effects (the effect of the previous transport processes on the current state of media) for the system to fix the failures of the classical equation. Because the fractional equation includes all-around movements, it is more accurate than the classical model ^[6,7,8,9]. The fractional advection-diffusion equation was used for the first time in 1990 to describe the physical phenomenon of the anomalous diffusion of particles in disordered systems. After that, other researchers used the fractional equation to describe various physical phenomena until 2000 ^[5].

The fractional advection-diffusion equation is a powerful tool for modeling transport phenomena in complex systems. For example, one of the applications of the advection-diffusion equation is modeling the transport of passive tracers in groundwater in porous media ^[10]. Apprehending the advection-diffusion instrument is important for forecasting and overseeing the transport procedures of substances in different environments. However, there are many transport procedures in complex systems controlled by non-Brownian diffusion, which is well described by fractional calculus, and fractional advection-dispersion equations recently play an important role in the modeling of transport processes in various physical ^{[11,12,13,14,15,16,17,18,19]} and natural ^{[20,21,22,23]} systems. However, solving such equations is not always possible, because finding an answer for the fractional advection-diffusion equation is more difficult than the classical advection-diffusion equation. In addition, the fractional advection-diffusion equation may not be appropriate for all media, therefore, it is crucial to choose appropriate equations to represent advection-diffusion processes in particular systems ^[24,25].

In general, various numerical methods can be employed to solve the fractional advection-diffusion equations, including finite difference, finite element, spectral, and Monte Carlo methods. However, the use of the Legendre-shifted polynomials ^[24] to approximate the fractional advection-diffusion equation is limited. The Legendre-shifted polynomials constitute a specific type of polynomials, which are suitable for approximating functions with special boundary conditions. The main feature of these polynomials is to calculate a more accurate solution using fewer terms than other methods ^[26,27]. A type of multi-time fractional advection-diffusion (MT-FADE) equation with boundary and initial conditions is defined as follows:

with initial $v(x, 0) = z(x), \; 0\leq x \leq 1$ and boundary conditions $v(0, \tau) = r(\tau), \; v(1, \tau) = s(\tau), \; 0\leq \tau \leq T$. In Eq (1.1), the fractional factor $\mathrm{P} (D_{\tau})$ is described as $\mathrm{P} (D_{\tau})v(x, \tau) = \sum\limits _{e = 0}^{\epsilon}b_{e} {}_0^{C}D^{\alpha _{e}}_{\tau} v(x, \tau), \; \epsilon \in \mathbb{N}$ with $0 < \alpha _{e} < \alpha _{e-1} < \ldots < \alpha _{0} \leq 1$ and $b_{e}\geq 0$. Terms $\frac{\partial ^{2\beta _{i}} v(x, \tau)}{\partial \vert x\vert ^{2\beta _{i}}}, \; i = 1, 2$ are the Riesz fractional derivative with respect to $x$ are described as

where ${}_xD^{2\beta _{i}}_{L} v(x, \tau)$ and ${}_xD^{2\beta _{i}}_{R} v(x, \tau)$ describe the left and right fractional Riemann-Liouville derivatives specified as below:

where $n-1 < 2\beta _{i}\leq n, \; n\in \mathbb{N}$.

It is noted that the well-known Eq (1.1) consists of some classic fractional instances. Using $e = 0$ in the fractional factor $\mathrm{P} (D_{\tau})$ results in traditional advection-diffusion equation ^[28,29]. Taking $K_{1} = 0, \; \beta _{2} = 1$, the multi-term time-fractional diffusion equations are obtained ^[30,31]. Taking $\alpha _{0} = 2, \; \alpha _{1} = 1$ and $K_{1} = 0$ in the case $1 < \alpha _{e} < \alpha _{e-1} < \ldots < \alpha _{0} \leq 2$ results in a spatially fractional Telegraph equation ^[32].

The MT-FAD model describes the changes in a quantity on multiple time scales. It is a generalization of the classical advection-diffusion equation. We often use these equations to model physical phenomena such as heat conduction, diffusion, and wave propagation ^[30,33]. There are two general analytical and numerical solutions to solve these equations. There are analytical ways for some MT-FAD model, although most of these answers are not accurate. In using numerical methods, we usually use finite difference methods ^[33,34]. In this way, we first discretize time and space and then solve the resulting system of equations using numerical techniques. In most cases, finite difference methods are the only possible way to solve MT-FAD model ^[30]. Many authors discussed the solution of MT-FAD model with numerical and analytical methods ^[35,36]. It is worth mentioning that Legendre-shifted polynomials are a specific type of polynomials, and they are suitable for approximating functions with special boundary conditions. The main feature of these polynomials is to calculate a solution more accurately using fewer terms than other methods. In this study, we proposed a novel approach to solve multi-term time-fractional advection-diffusion equation (1.1). We introduce Lagrange linear interpolation for temporal discretization and Legendre approximation for spatial discretization in Section 2. The highlights of the proposed numerical methods are unconditional stability and convergence, which are presented in Section 3. To demonstrate the effectiveness of the proposed numerical methods, we conducted some numerical experiments illustrated in Section 4. Eventually, the conclusion is summarized in Section 5.

2.   Discretization methods

This section consults the discretization of relation Eq (1.1). To do so, let $\delta h$ and $\delta\tau = \frac{T}{M}$ denote the sizes of the spatial and temporal discretization steps, respectively. We describe the partition in space as $\lbrace x_{i} \rbrace _{i = 0}^{M}$ that is the root of the basis polynomials and the partition in time as $\tau _{j} = j \delta\tau, j = 0, 1, \ldots, M$. We denote the discrete values of the variables $v(x, \tau)$ and $q(x, \tau)$ as

Considering the fact that $0 < \alpha _{e}\leq 1$, the Caputo fractional derivative can be utilized instead of the Riemann-Liouville.

2.1.   Time-discretization scheme

To start off on discretization, we use $S_{2}$ method in paper ^[37]. First, we must approximate unknown function for the case that $j = 0$ in the intervals $[\tau _{0}, \tau _{1}]$ by Lagrange linear interpolation, then, for the rest of the nodes $j \geq 0$ in the interval $[\tau _{j-1}, \tau _{j}]$ we apply Lagrange square interpolation. Practically, we have three nodes to conduct the approximation.

where for $k = 1$ and $k = 2$ the coefficients $\mathcal{S}^{\alpha _{e}}_{k, j}$ is as

and for $k\geq 3$, we have

in which

Applying the discretization of the relation (2.1) for the left side of Eq (1.1), we have

Letting $V_{i}^{k}$ as the approximate solution $v_{i}^{k}$ in Eq (2.2), then we have

2.2.   Spatial discretization scheme with Legendre approximation

Now, to get the full discrete of Eq (2.3), we employ the following series:

where $N$ and $M$ are the number point of spatial and temporal directions, respectively. Moreover, $\mathbb{L}_{\iota}^{\circledS}(x)$ is the shifted orthogonal polynomials in the domain $[0, 1]$ that is defined in paper ^[38] as

where

Then the unknown coefficients in Eq (2.4) are defined below:

Here is a formula for the approximation of the fractional derivative $\mathbb{L}_{i}^{\circledS}$ is carried out that we define it with the symbol $\mathbb{L}_{i}^{\gamma, \circledS}$. Suppose $\gamma > 0$, by using the Caputo linearity property, we get

in which $\lfloor \gamma \rfloor$ and $\lceil \gamma \rceil$ are the floor and ceiling of the fractional term $\gamma$ and

Notice that for $0\leq i < \lceil \gamma \rceil$ we have $\mathbb{L}_{i}^{\gamma, \circledS}(x) = 0$. Substituting the operators (2.5) and (2.6) in (2.3), we approximate $V^{k}(x)$ with respect to $x$ as below:

where

We use the collocation manner to obtain the coefficients $\sigma _{i}^k$ of relation (2.7). For this objective, we get the roots of the shifted Legendre polynomials, $\mathbb{L}_{i}^{\circledS}(x)$, as the collocation points and substitute them in Eq (2.7) to obtain linear equations at each time step $k$. $\mathbb{L}_{N-1}^{\circledS}(x)$ has $N-1$ roots, which needs two more conditions to obtain a system of linear equations with $N+1$ equations, which can be obtained using the following boundary conditions for $i = 0, 1, \ldots N$ and $k = 1, 2, \ldots M$.

To start the iterative method, we need the initial condition as

where

Below, we briefly describe the step-by-step numerical algorithm.

Step 1. Generate the basis functions, $\lbrace \mathbb{L}_{\iota}^{\circledS}(x) \rbrace _{i = 0}^{N}$, select the time step, $\delta \tau$, and produce a set of roots of $\mathbb{L}_{\iota}^{\circledS}(x)$, $\lbrace x_{i}\rbrace _{i = 0}^{N}$, as the collocation points.

Step 2. Define the approximate fractional time derivative (2.1) and construct the matrix $A_{e, i}(x)$, the vectors matrix $B_{e, i}(x)$ and $Q(x)$ in relation (2.7).

Step 3. Compute the coefficient matrices, $\sigma _{i}^k$, in relation (2.7) with considering the initial and boundary conditions.

Step 4. Finally, by substituting these coefficients in relation (2.4), one estimates the approximate solution.

3.   Study of the stability and convergence approach

This section investigates the order of convergence and the stability of the numerical method. Using the weak design (2.3) and defining the following function space in Hilbert space $L^{2}(\Omega)$ in $\Omega$ and standard norm $\Vert \theta(x)\Vert _{2}^2 = \langle \theta(x), \theta(x)\rangle$, we prove two theorems that describe the precision and efficiency of the numerical strategy explained in the earlier section.

where $D^\alpha$ is the fractional derivative.

Let ${V}^{k}(x)$ and $\overline{V}^{k}(x)$ be the exact and approximate solution of Eq (2.3), respectively. Multiplying $\varepsilon ^{k}(x)$ and integrating on $\Omega$ in the relation (2.3) and denoting $\varepsilon ^{k}(x) = \overline{V}^{k}(x)-{V}^{k}(x)$, we get the weak form of the relation (2.3) as

Before giving the stability and convergence of the weak scheme (3.1), we first give some lemmas.

Lemma 3.1. (See ^[39]) For $\forall f, g \in H^n_{\Omega}$ and $\forall x \in \mathbb{R}$, we have

Lemma 3.2. For $\forall \varepsilon \in H^n_{\Omega}$, we have

Proof. Using Lemma 3.1 and features of the interior product, we get

Theorem 3.1. When $0 < \alpha _{e}, \beta _{i} < 1, \; e = 0, 1, \ldots, \varepsilon, \; i = 1, 2$, the scheme (2.2) is unconditionally stable.

Proof. According to the Lemmas 3.1 and 3.2, the second and third terms of Eq (3.1) can be removed, then we gain

and using the Cauchy Schwarz inequality, one get

Since $\frac{ b_{e} { \delta \tau}^{-\alpha _{e}}}{\Gamma(2-\alpha _{e})} > 0$ and $1 < \mathcal{S}^{\alpha _{e}}_{k, k}\leq \frac{3}{2}$ from lemma of paper ^[40], we denote the above relation as form

where $\Gamma _{e} = \frac{ b_{e} { \delta \tau}^{-\alpha _{e}}}{\Gamma(2-\alpha _{e})}$. Applying the theory of induction applies to all $k$, we know that the above relation satisfies as below:

Using the properties of the expansion coefficients in the lemma of paper ^[40] that describe the characteristics of approximation coefficient $1 < \mathcal{D}_{k, k}\leq \frac{3}{2}$ and $\sum_{j = 0}^{k-1}\mathcal{S}_{k, j} = -\mathcal{D}_{k-1, k-1}$, it can be written as

where $c = \sum\limits _{e = 0}^{\epsilon}\Gamma _{e}\vert \mathcal{D}_{k, k}\vert$ is the positive constant. Then it shows that the discredited scheme (2.3) is unconditionally stable.

The convergence of the scheme (2.3) is given in the following theorem.

Theorem 3.2. Let $v^k(x) \in H^2_{\Omega}$ be the exact solution of (2.2), $V^k(x) \in H^2_{\Omega}$ be the solution of the semi-scheme (2.3). Then there exists a positive constant $C$ such that

Proof. Subtracting (2.2) from (2.3) and denoting $\xi ^{k}(x) = v^k(x)-V^k(x), k = 1, 2, \ldots, M$, we can obtain that

where $R^{k} = \mathcal{O}({ \delta \tau} ^ {3-\max {\alpha _{e}}})$. From Lemmas 3.1 and 3.2, we know that the second and third terms of the above relation are negative, thus it can be rewritten as

By the method employed to prove the earlier theorem, we know that there exists a positive fixed term $c$ such that

Since $\Vert \xi ^{0}(x)\Vert = 0$, then

This concludes the proof. □

4.   Numerical experiments

We have conducted some experiments to demonstrate the effectiveness of our proposed methods. The numerical calculation is carried out on a Dell Inspiron personal computer equipped with an Intel (R) Core i72630QM 2.00GHz processor and 8GB of memory, using Wolfram Mathematica version 11. Our unique approach helped to address various issues, and we evaluated the maximum norm error $L_{\infty}$-norm and $L_{2}$-norm error between the numerical solution and the exact solution at $T = 1$ as below:

respectively. In this case, $v_{i}^{M}$ is the exact solution and $\widehat{v}_{i}^{M}$ is the numerical solution with mesh step size $\delta h$ and grid point $(x_i, \tau _j)$, where $i = 0, 1, \ldots, N$ and $j = 0, 1, \ldots, M$. Moreover, the convergence rate is calculated as follows:

This formula measures how the error between the exact and numerical solutions decreases as the mesh step size $\delta h$ is decreased. A higher convergence rate indicates that the numerical solution is converging to the exact solution more rapidly.

The numerical experiments conducted to evaluate the precision of the new approach for discretizing equations. Equation (1.1) is presented in the tables. The convergence rate $Rate_\infty$ and maximal norm error $E_\infty(\delta h, \delta\tau)$ are demonstrated for various mesh step sizes $\delta h$. The findings suggest that the innovative approach can achieve temporal convergence accuracy of $O(3 - max(\alpha_ i))$, where $i = 1, 2$. Furthermore, the accuracy of the outcomes was compared to a previous study ^[33], and the proposed method demonstrated higher accuracy with lower errors than the method used in ^[33].

Example 4.1. Consider the following model with the diffusion coefficients are given $K_{1} = K_{2} = 2$ that the source term $q(x, \tau)$ is obtained from the exact solution $v(x, \tau) = \exp(-\tau)\sin(x)(x^2-x+1)$.

where the initial condition is $v(x, 0) = \sin(x)(x^2-x+1),$ and the boundary conditions are $v(0, \tau) = 0$ and $v(1, \tau) = \sin(1)$. The computing results of this model are shown in Table 1, wherein the maximal norm and $E_{2}$-norm errors have been provided; and the convergence rate has been found that it is equal to $\mathcal{O}(3-\max (\alpha _{i})), \; i = 1, 2$ for all temporal discretization levels. From this table, we can see that the convergence rate is $2.2$. Then the largest numbers of the parameters $\alpha _{i}$ are required to determine the convergence rate.

Example 4.2. We first consider the one-dimensional Eq (1.1) in the domain $(x, \tau)\in [0, 1]\times [0, 1]$, with $T = 1, K_{1} = K_{2} = 1$, and source term

To be specific, the initial condition is $v(x, 0) = 100(x^2-x^3),$ and the boundary conditions are $v(0, \tau) = v(1, \tau) = 0$. The precise solution for this issue can be expressed as $v(x, \tau) = 100(1 + \tau^2)(x^2 - x^3)$. We have carried out numerical experiments and presented the results in Tables 2 and 3. Table 2 shows the maximum norm and $E_2-norm$ errors for the new method, which uses parameters $\beta_1 = 0.2$ and $\beta_2 = 0.7$. Based on our analysis, we have determined that the convergence rate can be expressed as $O(3 - max(\alpha_ i))$.

Table 3 shows cases of a comparison between the error and convergence order of our novel method and the technique outlined in ^[33]. We used parameters $\alpha_0 = 0.5, \; \alpha_1 = 0.2, \; \beta_0 = 0.3$ and $\beta_1 = 0.8$. Our results demonstrate that our new method is more efficient than the method presented in ^[33]. When $\alpha_0 = 0.5$, the convergence order for our new method is $2.5$, while the convergence order for the method presented in ^[33] is only $1.5$. Additionally, the error of our new method is significantly lower than the error of the method presented in ^[33]. In conclusion, our new method outperforms the method presented in ^[33] in terms of efficiency and accuracy.

The stiffness matrix's condition numbers obtained for space partitions $N = 5$ and $N = 7$ with different values of $\beta_1$ and $\beta_2$ are shown in Figure 1. The two curves indicate that the new algorithm is robust in discretizing the equation to obtain the resulting systems. The condition numbers decrease as the number of discretization points increases. In Figure 2, the absolute error is shown with $N = 5, \; \alpha_1 = 0.5, \; \alpha_2 = \beta_1 = 0.2$ and $\beta_2 = 0.7$ parameters, while $M$ values are varied. As the number of discretization points increases, the absolute error decreases. In conclusion, the new algorithm effectively discretizes the equation and the resulting systems. When you increase the number of discretization points, it leads to a decrease in both the condition numbers and absolute error.

Example 4.3. The second example is similar to the first one, with the same accurate resolution of $v(x, \tau) = 100(1+\tau ^2)x^2(1-x)^{2}$ and diffusion coefficients $K_{1} = K_{2} = 1$. The source term for this example is

Tables 4 and 5 present the results of numerical experiments conducted on a new method. In Table 4, wherein the maximal norm and $E_{2}$-norm errors have been provided with parameters $\beta _{1} = 0.45, \; \beta _{2} = 0.7$, and the convergence rate has been found that it is equal to $\mathcal{O}(3-\max (\alpha _{i})), \; i = 1, 2$ for all temporal discretization levels. From this table, we can see that the convergence rate is $2.5$ when $\alpha _{1} = \alpha _{2} = 0.5$ but it is $2.35$ for $\alpha _{1} = 0.2, \; \alpha _{2} = 0.65$. Then, the largest numbers of the parameters $\alpha _{i}$ are required to determine the convergence rate. Comparing the error and convergence order of the method ^[33] with the new method is represented in Table 5 with the parameters $\alpha _{0} = 0.7, \; \alpha _{1} = 0.4, \; \beta _{0} = 0.3$ and $\beta _{1} = 0.85$. From this table, we can see that the proposed method is much more efficient than method of ^[33]. Moreover, the error of the new method is much less compared to the paper of ^[33].

Figure 3 illustrates the absolute error at different values of M, with parameters $N = 5, \; \alpha_1 = 0.3, \; \alpha_2 = 0.6, \; \beta_1 = 0.4$ and $\beta_2 = 0.8$. The figure shows that increasing the number of discretization points decreases the absolute error. In conclusion, the new technique is more efficient and accurate than the method from ^[33].

5.   Conclusions

In this study, we have introduced a novel approach for solving multi-term time-fractional advection-diffusion equations that utilizes a square interpolation method to discretize the temporal fractional derivative and a spectral method and employs Legendre polynomials to discretize the spatial derivatives. We then construct a matrix-approximate preconditioner to accelerate the solution of the discretized system. Our approach has been proven to be stable without any conditions and has a convergence rate of $O(3- max(\alpha_i))$. To the best of our knowledge, the theoretical proof obtained from this study represents the first results for these numerical schemes. We also presented experimental results that demonstrate the effectiveness of our proposed method in comparison with the methods shown in ^[33]. However, the proposed method is inherently a limitation, wherein the coefficient matrix of the linear system obtained from the discretization methods is not a sparse matrix, and this may increase a significant amount of memory and slow down the processing of that data for a large matrix of three-dimensional and large scale problems. In the future, we aim to expand our proposed approach to address variable coefficients in higher-order fractional differential equations and also plan to develop some fast numerical methods for these equations.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors would like to thank the support from the National Research Foundation of Korea under grant number NRF-2020K1A3A1A05101625 and from the research support for foreign professors of the College of Engineering at Seoul National University.

Conflict of interest

The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.