An efficient modified hybrid explicit group iterative method for the time-fractional diffusion equation in two space dimensions

1.   Introduction

In recent years, fractional calculus has gained a great deal of attention from the research community due to its widespread applications in physics, chemistry, fluid mechanics, viscoelasticity, finance, control systems and other areas of science and engineering. The phenomena in the aforesaid fields can be modelled very successfully by means of equations containing fractional derivatives and fractional integrals; and therefore, the investigation of fractional differential equations has become a hot topic for many researchers. Fractional kinetic equations including Fokker-Planck equation, fractional diffusion equation, fractional cable equation and fractional advection-diffusion equation are proved to be powerful instruments for modelling transport dynamics in complex systems such as electrodiffusion of ions in spiny dendrites, transport of proteins molecules, movement of a solution in an aquifer, pollution of underground water, movement under the effect of optical tweezers and more ^[1,2,3,4,5]. To get more insight into the applications of fractional calculus, the interested reader can refer to the classical books in ^[6,7], besides the recent books in ^[8,9].

In this work, we study the following two-dimensional time-fractional diffusion equation with a non-homogeneous source term:

where $f(x, y, t)$, $p(x, y, t)$ and $g(x, y)$ are known functions, $\Omega = (0, L)\times(0, L)$ is a bounded domain in $\mathbb{R}^2$, $\partial\Omega$ is the boundary, $a_x > 0$ and $a_y > 0$ are the diffusion coefficients, and $\alpha \; (0 < \alpha\le1)$ is the anomalous diffusion exponent. The term ${}^C_0D_t^\alpha u(x, y, t)$ stands for the Caputo fractional derivative of the function $u(x, y, t)$ that is defined as,

From the above definition, it can be observed that the fractional diffusion problems (1.1)–(1.3) corresponds to the classical diffusion model as $\alpha = 1$. Fractional diffusion equation is one of the most fundamental equations in the literature, and its capability to model many problems in science and engineering is a well established fact. Generally speaking, it has been widely used for describing random walks and modelling phenomena that are governed by anomalous diffusion in various fields such as porous systems, nuclear magnetic resonance and transport in fractal geometries ^[10,11,12]. Since most fractional differential equations are difficult to handle analytically, many researchers have resorted to numerical techniques, especially the finite difference method for solving fractional diffusion problems ^{[4,10,11,13,14,15,16,17,18,19,20,21,22,23,24,25]}. The complexity of fractional differential equations stems from the presence of non-local fractional derivatives that have the property of global dependence on time or space. This forms the principal obstacle to the development of efficient simulation algorithms in terms of CPU time and memory usage, especially for multi-dimensional problems ^[26,27,28]. To surmount this issue, techniques including fast Poisson solver ^[14], parallel implementation ^[15,29], preconditioning ^[15] and multigrid method ^[30] were suggested in the literature. Therefore, the question of how we can remove the fractional derivative from our computations to reduce the complexity and establish efficient solution algorithms is considered a big question arising in the numerical simulations of fractional differential equations. In this regard, Ren et al. ^[31] have introduced a relatively new approach for solving Caputo-type fractional differential equations. In this approach, a Laplace transform technique is proposed to approximate the fractional differential equation by its corresponding integer-order differential equation, which can be solved with less effort by using a suitable numerical or analytical method (see ^[32,33]). Recently in ^[34], a new hybrid standard point (HSP) iterative method based on a combination of the Laplace transform technique and implicit difference scheme has been proposed to solve the time-fractional diffusion problems (1.1)–(1.3). The authors showed that their method produces accurate numerical solutions with an optimal computational complexity of $O(N)$ and optimal memory requirement of $O(M_s)$, where $N$ and $M_s$ are the total number of time levels and spatial unknowns, respectively.

It is well known that the explicit group iterative methods ^{[35,36,37,38,39]} based on finite difference approximations are proposed for the numerical solutions of integer-order differential equations. From one side, explicit group methods are unconditionally stable as the conventional implicit schemes. On the other side, they reduce the number of spatial unknowns taken in the iterative process and thus lead to computationally efficient algorithms. Although the grouping methods have been successfully employed to abroad spectrum of classical differential models, very little work has been reported to deal with fractional differential equations using grouping techniques (see ^{[40,41,42,43]}). In addition, the handling of nonlinear and variable order fractional differential equations by explicit group methods, besides their parallel implementation, is another topic that is still at its infancy. Therefore, this subject needs a major development.

The main purpose of this paper is to combine the Laplace transform technique with an explicit group scheme to solve the two-dimensional time-fractional diffusion Eqs (1.1)–(1.3). The unconditional stability, convergence and solvability of the resulting method, namely the modified hybrid explicit group (MHEG) method are rigorously proved. To illustrate the efficiency and feasibility of the proposed method, the fast HSP iterative method based on the recent work in ^[34] is also presented.

The rest of this paper is organized as follows. In Section 2, we briefly review the existing HSP iterative method ^[34] for solving the problems (1.1)–(1.3). In Section 3, the MHEG iterative method is proposed, and its unconditional stability, convergence and solvability are rigorously proved in Section 4. In Section 5, numerical simulations are performed to indicate the efficiency of the proposed method and support our theoretical analysis. Finally, this work ends with a brief summary in Section 6.

2.   Overview of the existing HSP iterative method

Since the Caputo fractional derivative in (1.1) is non-local and has the character of history dependence, conventional difference schemes (such as an implicit or explicit scheme with a particular discretiztion formula for the fractional derivative) for the problems (1.1)–(1.3) result in costly simulations in terms of CPU time and memory consumption. In order to surmount the computational challenge, and utilizing the Laplace transform technique, the Caputo's time fractional derivative is approximated as follows ^[34]:

By substituting (2.1) into (1.1), the original time-fractional diffusion Eq (1.1) is reduced to its corresponding integer-order partial differential equation (PDE) with the following initial and boundary conditions ^[34]:

where $A_x = a_x/\alpha$, $A_y = a_y/\alpha$ and $r = 1/\alpha$ are positive constants. Next, we let

where $h = L/M$ and $\tau = T/N$ are the uniform space and time step sizes, respectively. We also introduce the grid functions given by

Here, the notations $\delta_t U_{i, j}^{k+1}$, $\delta_x^2 U_{i, j}^{k+1}$ and $\delta_y^2 U_{i, j}^{k+1}$ are defined as follows:

Now, applying the difference operators in (2.5) to the approximating PDEs (2.2)–(2.4) and disregarding the truncation errors, we obtain the implicit difference scheme in the following form:

where $q_1 = A_x\tau/h^2$, $q_2 = A_y\tau/h^2$ and $u_{i, j}^k$ is the numerical approximation of the function $U_{i, j}^k$ after neglecting the truncation errors.

Considering the fact that the spatial and temporal discretizations usually lead to large and sparse systems of linear equations to be solved, iterative solvers, such as Gauss-Seidel method, are more practical and economic than the direct solvers. Here, and at each time level, the HSP method functions through the iterative evaluation of solutions at all mesh points $\blacklozenge$ shown in Figure 1 by using Eq (2.7) until certain convergence is achieved. Once the converged solution values $u^{k+1} = (u_{11}^{k+1}, u_{1, 2}^{k+1}, \ldots, u_{1, M-1}^{k+1}, \ldots, u_{M-1, 1}^{k+1}, u_{M-1, 2}^{k+1}, \ldots, u_{M-1, M-1}^{k+1})^T$ are attained, they are utilized as an initial guess for the next time level $u^{k+2} = (u_{11}^{k+2}, u_{1, 2}^{k+2}, \ldots, u_{1, M-1}^{k+2}, \ldots, u_{M-1, 1}^{k+2}, u_{M-1, 2}^{k+2}, \ldots, u_{M-1, M-1}^{k+2})^T$. Such process is repeated until it reaches the targeted time level. The HSP iterative method has been shown to reduce the memory requirement, computational complexity as well as the CPU time greatly when compared to the fractional standard point (FSP) iterative methods derived based on the conventional implicit difference schemes. For extra details, please refer to ^[34]. Seeking for more efficient solution algorithm, we propose our method in the next section.

3.   Formulation of the MHEG iterative method

Consider the standard mesh of points with the spatial step size $2h = 2L/M$ illustrated in Figure 2. We begin this section by discretizing the derivatives of the space variables in (2.2) around these $2h$ spaced points. Recalling the $\delta_t$ definition, the spatial difference operators $\delta_x^2$ and $\delta_y^2$ are defined as follows:

Applying the difference operators in (3.1) to the approximating PDE (2.2), we obtain

After simplification and imposing the initial and boundary conditions, the following implicit scheme is attained:

where $u_{i, j}^k$ is the approximate numerical solution of $U_{i, j}^k$ after omitting the truncation error terms. The MHEG method can now be constructed by applying Eq (3.3) to a group of four mesh points of the solution domain. This will result in the following ($4\times4$) system of equations:

where

The coefficients matrix in (3.6) can be inverted to get the MHEG equation given as below:

where

Back to Figure 2, it can be seen that the discretized solution domain comprises three different types of mesh points ($\blacklozenge$, $\square$, ○). Here, the MHEG method functions through the iterative evaluation of solutions at $\blacklozenge$ mesh points by utilizing Eq (3.7) until a certain convergence criterion is met. Afterwards, the solution values at the remaining mesh points of types $\square$ and ○ are computed directly once in a particular sequence. The computational cost of the iterative method still depends on the number of unknowns involved in the iteration process. For the MHEG method, it is obvious that the number of mesh points taken in the iteration process is lesser than that of the prespecified HSP method, which is expected to speed up the convergence rate of the proposed method. For convenience, the MHEG method in combination with the Gauss-Seidel iterative solver is elaborated in Algorithm 1.

4.   Stability and convergence

In this section, we turn to consider the stability and convergence of the MHEG scheme defined by (3.7) by using the matrix analysis method. For the convenience of our analysis, three remarks are given at first.

Remark 4.1. Let $z\in \mathbb{N}$. The symbol $\|A_{z\times z}\|$ denotes the infinity norm of matrix $A = [a_{i, j}]$ which is defined as,

Remark 4.2. $A_{z\times z}$ is called a stricatly diagonally dominant (SDD) matrix if $|a_{i, i}| > r_i(A)$, where $r_i(A)$ is the $i$-th deleted absolute row sum (i.e., $r_i (A) = \sum_{j\neq i, j = 1}^{z}\left| {{a_{i, j}}} \right|$, $1\le i\le z$).

Remark 4.3. Given that $A_{z\times z}$ is a SDD matrix, then $A_{z\times z}$ is non-singular and the following norm's bound does hold ^[44]:

4.1.   Stability analysis

Here, we analyze the stability of the MHEG scheme (3.7). For simplicity and without loss of generality, we assume that $q_1 = q_2 = q = \tau/h^2$. Indeed, the MHEG scheme (3.7) can be represented in the following matrix form:

where $u^k$is an $\frac{(M-2)^2}{4}$-dimensional vector expressed in the block form

and

Taking note of the above block matrices, the matrix representation in (4.1) can be rewritten in the more elaborated general form

Next, we prove the following result:

Theorem 4.1. The MHEG scheme defined by (4.2) is unconditionally stable.

Proof. Suppose $u^k$ and $\tilde{u}^k$ are respectively the exact and approximate solutions of (4.2). Then the error at time level $k$ is expressed as $e^k = u^k-\tilde{u}^k$. From Remarks 4.2 and 4.3, it follows that $A$ is invertible and therefore Eq (4.2) can be rewritten as

The error function $e^k$ satisfies (4.3) and leads to the following round-off error equation:

where

and $\psi_{i, j}^{k+1} = u_{i, j}^{k+1}-\tilde{u}^{k+1}.$

In order to demonstrate the stability, we shall prove that $\|e^{k+1}\|\le\|e^0\|$ for all $0\le k\le N-1$. To this end, the mathematical induction will be used. For $k = 0$, we obtain

Using the fact that the matrix and the vector infinity norms are consistent, yields

Define $r_i(A)$ as in Remark 4.2. Then using Remark $4.3$, we get

Next, suppose that

We will show the above inequality is true for $s = k$. Noticing (4.4) and (4.5), we get

Thus, the proof is completed.

4.2.   Convergence analysis

We now analyze the convergence of the MHEG scheme (3.7) by following a similar approach to that in the previous subsection. As a start, we suppose that the truncation errors on each four-point group of mesh points at any time level are represented by the block vector of the following form:

Noticing Eq (3.2) and since $i$, $j$ and $k$ are finite, there is a positive constant $C^*$ such that

Subtracting (4.2) from the following equation:

will result in error equation given by

where

and $\phi_{i, j}^{k+1} = U_{i, j}^{k+1}-u_{i, j}^{k+1}$.

The convergence property is given in the next theorem.

Theorem 4.2. Assume problem (1.1) has a smooth solution $u(x, y, t)$. Then, the MHEG scheme defined by (3.7) is convergent and the convergence order is $O(\tau+h^2)$.

Proof. We will use mathematical induction for the proof. For $k = 0$ and using that $E^0 = 0$, we obtain

In view of of Remark 4.3 and utilizing (4.6), yields

where $C_0 = C^*/(1+(r-1)\tau)$.

Now, assume that

We prove the above result is true for $s = k$. Noticing (4.7) and (4.8), we get

where $C_k = C_{k-1}+C^*$ as $\lim_{n \to \infty}\tau = 0$.

Hence, the proof is completed by induction.

4.3.   Solvability

Theorem 4.3 The MHEG scheme defined by (3.7) is uniquely solvable.

Proof. In view of (4.1), The coefficients matrix $A$ of the MHEG scheme (3.7) is an SDD matrix. Consequently, and based on Remark 4.3, it follows that $A$ is a non-singular matrix. This proves the existence and uniqueness of the solution of the MHEG scheme.

5.   Numerical simulations

In this section, we report on numerical simulations for (1.1)–(1.3). In this endeavor, two test problems with known exact solutions are given to demonstrate the accuracy and efficiency of the newly developed MHEG method. For comparison purposes, the fast HSP method proposed recently in ^[34] is applied to solve the two test problems. All simulations using the both methods are performed in Mathematica 11.3 software and run on a computer with an Intel (R) Core (TM) i7-8550U CPU and 8.00 GB of RAM. In these simulations, the Gauss-Seidel iterative solver with a stopping criterion selected to be $10^{-5}$ are employed to generate the corresponding numerical results. To get more insight into these results, and based on the total number of arithmetic operations to be implemented before and after the convergence process, an analysis on the computational complexity of the MHEG and HSP iterative methods is established and presented in Table 1.

Test problem 5.1. We consider the following diffusion equation of fractional order:

subject to the boundary and initial conditions

in the spatial domain $(0, L)\times(0, L) = (0, 1)^2, \ T = 1$, and the exact analytic solution is $u(x, y, t) = t^2\sin(x)\sin(y).$

In this test, numerical outputs including maximum absolute error $E_{max} = \max_{1\le i, j\le M-1}|u(x_i, y_j, t_N)-u_{i, j}^N|$, CPU time (in seconds), number of iterations ($Ite$) as well as total number of arithmetic operations are computed for different values of $\tau, \ h$ and $\alpha$ and listed in Table 2. It can be observed that the MHEG iterative method results in significantly faster simulations without compromising too much of the accuracy when compared to the HSP iterative method. To illustrate this further, Figures 3 and 4 show respectively the graphs of CPU time and total operations of both methods against various mesh sizes. Form these figures, it is apparent that the shape of the computing time plots is consistent with the shape of the total computing effort plots. This indicates that the experimental results are in good agreement with our theoretical analysis. The comparison of the exact and the numerical solutions besides the graphical error representation when $\tau = h = 1/18$ and $\alpha = 0.55$ are depicted in Figures 5 and 6, respectively. A comparison of the numerical errors of the test problems for different fractional orders shows that the accuracy of the proposed method is greatest at $\alpha = 0.55$, followed by $\alpha = 0.75$ and $\alpha = 0.95$. It is clear that the numerical solution matches well with the exact solution.

Test problem 5.2. We consider another diffusion equation of fractional order as follows:

with the boundary and initial conditions

in the spatial domain $\Omega = (0, 1)^2, \ T = 1$, and the exact analytic solution $u(x, y, t) = t^2e^{(x+y)}.$

The numerical results of solving this problem for various values of $\tau$, $h$ and $\alpha$ are presented in Table 3. Figure 7 highlights the CPU time plots versus several mesh sizes. In view of Table 3 and Figure 7, it can be seen that the MHEG iterative method reduces the CPU time greatly without jeopardizing the accuracy of numerical solutions compared to the HSP iterative method. This can be attributed to the reduction in the number of executed arithmetic operations that leads to a lower computational complexity as shown in Figure 8. Figure 9 introduces the comparison between the exact and the numerical solutions, whereas Figure 10 presents the 3D plot of the maximum errors when $\tau = h = 1/30$ and $\alpha = 0.75$. Again, these demonstrate the effectiveness of the proposed MHEG method. Table 4 illustrates the experimental convergence order in the temporal direction using the formula $\rho(\tau, h) = \log_2\left(E_{max}(2\tau, h)/E_{max}(\tau, h)\right)$. It can be observed that the performance is in good agreement with the theoretical analysis.

6.   Conclusions

In this paper, a new MHEG iterative method for solving the two-dimensional diffusion equation with time fractional derivative has been developed. The unique solvability, unconditional stability and convergence are proved by the matrix analysis method. Moreover, the feasibility and efficiency of the proposed numerical method are confirmed through the computational outputs drawn from the conducted numerical simulations. The advantages of being uncomplicated, easy to implement and diminishing the amount of computational complexity indicate that the proposed method and numerical analysis in this article can be extended to solve other types of non-linear and variable order fractional differential equations. In addition, the parallel implementation of the proposed MHEG method is an interesting line to study which can be considered as a part of future work.

Acknowledgments

This research was supported by Ministry of Higher Education Malaysia, Fundamental Research Grant Scheme with Project Code: FRGS/1/2019/STG06/USM/01/1. We would like to express our gratitude to the anonymous reviewers for their valuable comments.

Conflict of interest

The authors have no conflict of interest.