On numerical solution of two-dimensional variable-order fractional diffusion equation arising in transport phenomena

Nomenclature

VO: Variable-order.

CO: Constant-order.

C-N: Crank-Nicolson.

EDG: Explicit decoupled group.

CPU: Central processing unit.

VO-TFDE: Variable-order time fractional diffusion equation.

Symbols

$M_1$: Number of spatial steps in the $x$-direction.

$M_2$: Number of spatial steps in the $y$-direction.

$N$: Number of temporal steps.

$\Delta x$: Length of spatial step along $x$-direction.

$\Delta y$: Length of spatial step along $y$-direction.

$\tau$: Temporal step size.

$T$: Final time.

$K_1$: Diffusion coefficient.

$K_2$: Diffusion coefficient.

$\alpha$: Fractional order.

1.   Introduction

A plethora of problems in engineering applications can be described by integer-order differential equations ^[1]. On the other hand, fractional differential equations have numerous applications in many sciences, which makes them one of the most prominent tools of fractional calculus that have attracted the interest of scholars and researchers in recent years. For instance, Lozynskyy et al. ^[2] applied the Caputo-Fabrizio operator to describe the elastic moment in a long shaft and produced an enhanced mathematical model in the form of a fractional ordinary differential equation for a two-mass system with concentrated parameters. Asjad et al. ^[3] successfully described the convection heat transfer in clay nanofluid by utilizing a Maxwell model. The authors used the newly developed hybrid fractional derivative in ^[4] to generalize the classical Maxwell model to the novel fractional one. Liu et al. ^[5] established a new Green-Ampt model containing a Caputo fractional derivative to describe an important phenomenon in geotechnical engineering; that is, the infiltration process in slopes. The researchers argued that the outcomes of their model were better than those of Chen and Young's model and very close to the actual experimental results. Kumar et al. ^[6] provided an opportunity for a better understanding of climate change by introducing a fractional-order model that describes the plankton-oxygen dynamics based on the Liouville-Caputo fractional derivative. Mustafa Inc et al. ^[7] considered a novel Covid-19 model based on the Caputo fractional operator to study the dynamics of the disease spread. The results of the numerical simulations were in good match with the real data and showed the advantages of the fractional derivative. Shen et al. ^[8] got more insight into the case where a patient has two infections at the same time by considering a fractional co-infection model under the Caputo fractional derivative.

Differential equations with constant-order (CO) fractional derivatives have been considered as powerful instruments for characterizing the long memory of dynamic processes and real-world phenomena. This means modeling the case in which the current state of the system depends on all its historical states rather than only some recent ones. More theoretical and comprehensive discussions of fractional calculus can be found in ^{[9,10,11,12,13,14]}. In recent years, many researchers have found systems exhibiting variable memory that cannot be described adequately by CO fractional operators. This fact led to the generalization of the CO fractional calculus to the so-called variable-order (VO) fractional calculus, where the derivatives of fractional operators are functions of time and/or space variables. Sun et al. ^[15] introduced an important comparative study between constant and VO operators in describing the memory property of systems. It was Samko and Ross ^[16] who gave the definitions of VO derivatives for the first time in 1993. Later, Lorenzo and Hartley ^[17,18] and Coimbra ^[19] constructed the mathematical framework of VO operators and discussed some of their potential applications in mechanics. These works paved the way for the development of new definitions, which served as a springboard for the field of VO fractional calculus. The list of the definitions in the literature includes: Riemann-Liouville, Caputo, Grunwald-Letnikov, Hadamard, Caputo-Fabrizio, Atangana-Baleanu and more. These definitions are reported to describe complex dynamic problems effectively. Consequently, with the extension of CO fractional models to their counterparts, VO models have grown in a remarkable manner. Examples are the diffusion equation ^[20], reaction-diffusion equation ^[21], Burgers equation cite{r21}, telegraph equation ^[23], mobile/immobile equation ^[24], Schrödinger equation ^[25] and more. For more details about the VO fractional calculus, the interested reader is referred to ^[26,27,28].

In the last few years, VO fractional differential equations have achieved tremendous success in modeling physical problems of evolutionary nature-type. The reason for this is the ability of VO operators to adjust the system's order based on its present or even its previous response effectively. For example, Sweilam et al. ^[29] proposed a VO fractional model to simulate the healing process of the cancer disease in biological systems. The established model involved VO Caputo derivatives with two control variables to reduce the number of cancer cells. Xiang and others ^[30] made use of the VO Caputo derivative and came up with a novel fractional model for describing the stress responses of glassy polymers. The authors noticed a relationship between the varying order and the behavior of microstructures, which opens up a new avenue for manufacturing new polymers with enhanced properties. Liu et al. ^[31] suggested a fractional damage creep model to depict the mechanical behavior of creep, i.e., the time-dependent behavior under constant stress. Their model relied on the Caputo derivative of the VO function, which is related to the relaxation time of the rheological model. In another study, Fei et al. ^[32] developed a Maxwell model with VO fractional derivative to describe the creep phenomenon of a salt rock. A fractional model with a VO time derivative for an economic system was proposed and analyzed in ^[33]. Fan et al. ^[34] constructed a VO fractional model and showed its superiority over the Anand model in reporting the mechanical properties of sintered nano-silver particles. To get more insight into the applications of VO fractional calculus, the interested reader is referred to the valuable works ^[35,36].

Fractional diffusion equations are among the most important models of fractional calculus that are applied to describe the anomalous diffusion phenomena in complex systems. Until date, a plethora of numerical methods have been developed to handle the diffusion models with CO fractional derivatives ^{[37,38,39,40,41,42,43,44,45,46]}. However, real data has revealed that the mean square displacement scales as a function of time and/or space variables, which cannot be characterized by CO fractional diffusion models ^[35]. As a result, fractional diffusion equations with variable orders were introduced to deal with the aforementioned issue. It is worth mentioning that solving variable-order fractional diffusion equations analytically is highly complex and cannot be achieved in many cases due to the variable exponent. This stimulated numerical researchers to design efficient approximation methods for solving these equations, namely, finite difference methods ^[47,48,49], finite element methods ^[20,50], wavelet methods ^[51], collocation methods ^[52,53], reproducing kernel methods ^[54] and others.

Due to the nonlocal property of VO operators, standard numerical schemes lead to costly simulations, particularly when a high-dimensional problem or mesh refinement is considered. Explicit group difference methods are combinations of finite difference schemes and grouping strategies on standard and skewed grids. In fact, explicit group difference methods can efficiently reduce the spectral radius of the iteration matrix and increase the convergence rates of the approximation methods. Besides, these methods can be implemented on parallel computers, which makes them very advantageous in practice. In the past few years, several researchers have successfully extended the explicit group difference methods from solving integer-order differential equations to handling various types of CO fractional differential equations ^{[55,56,57,58,59,60]}. For more details, the interested reader can refer to the cited studies. However, to the best of our knowledge, there have been no trials to utilize explicit group difference methods to determine numerical solutions of differential equations in the frame of VO fractional derivatives. In this line of reasoning, this article aims to provide a grouping-based numerical method to deal with the two-dimensional VO time fractional diffusion equation (VO-TFDE) of the following form:

where $\Omega = [0, L]\times[0, L]$ is a closed subset in $\mathbb{R}^2$ representing the solution domain, $\partial\Omega$ is its boundary, $K_1$ and $K_2$ are positive constants denoting the diffusion coefficients, $f(x, y, t)$, $h(x, y, t)$ and $g(x, y)$ are known smooth functions and $u(x, y, t)$ is the unknown function. $0 < \alpha(x, y, t) < 1$, ${}^C_0D_t^{\alpha(x, y, t)} u(x, y, t)$ is the VO Caputo derivative defined as

The above definition of the Caputo derivative has received a great deal of interest in the field of applications since it allows the traditional initial and boundary conditions that have solid physical interpretations to be imposed in the fractional models. The aim of this study is to construct a new grouping-based numerical method, namely, the explicit decoupled group (EDG) method for the two-dimensional VO-TFDE (1.1). We start by introducing a Crank-Nicolson (C-N) finite difference method, which serves as a reference scheme to examine the efficiency of the proposed method. Thereafter, the EDG method is developed based on a new difference scheme that is derived by rotating the standard grid and applying the Taylor series expansion. The resulting method is computationally advantageous since it considers only half of the mesh points in the iteration process at each time level, which means the numerical solutions can be computed in a short time. The stability, convergence and numerical implementations of the proposed methods are also presented in this work. As far as we know, explicit group difference methods with stability and convergence analyses for VO fractional models have not appeared in the literature.

The rest of this paper is as follows: The next section presents a C-N finite difference scheme to solve the two-dimensional VO-TFDE (1.1). In section three, we propose the EDG method for the mentioned equation. The stability and convergence analyses are included in sections four and five, respectively. Some numerical experiments are reported in section six to validate the efficiency and accuracy of the presented methods. Finally, a brief conclusion is expressed in section seven.

2.   The C-N difference scheme

In this section, we present a C-N numerical scheme for solving the VO-TFDE (1.1). Such a discretization scheme will reduce the considered problem into a system of simultaneous algebraic equations, which can be solved by a proper numerical technique. To this end, we let $x_i = i\Delta x$, $0\le i\le M_1$, $y_j = j\Delta y$, $0\le j\le M_2$, $t_k = k\tau$, $0\le k\le N$, where $M_1$, $M_2$ and $N$ are some positive integers. Here, $\Delta x = L/M_1$, $\Delta y = L/M_2$ and $\tau = T/N$ represent the spatial and temporal step sizes, respectively. For simplicity, define $u_{i, j}^k$ as the numerical approximation at the point $(x_i, y_j, t_k)$ and $f(x_i, y_j, t_k) = f_{i, j}^k$. Then, the discretization of the second order space derivatives is expressed as:

For the approximation of the Caputo time VO derivative, we utilize the following discretization scheme ^[61]:

where $0 < \alpha(x, y, t) < 1$, $\alpha^{i, j, k+1/2} = \alpha(x_i, y_j, t_{k+1/2})$ and

The local error term has the following bound,

Now, setting Eqs (2.1), (2.2) and (2.3) into (1.1), we get the following,

By omitting the local error terms and replacing $u_{i, j}^k$ with its approximation $\mathfrak{U}_{i, j}^k$, the desired result of the C-N difference scheme is obtained as

Let

The above difference scheme can be represented in matrix form as,

Here, $A$ and $B$ are pentadiagonal matrices given by

where

From the structure of the above matrices, it's evident that the matrix $A$ is a strictly diagonally dominant matrix. Hence, we deduce that the C-N difference scheme (2.5) is uniquely solvable. In the next section, the EDG method is proposed.

3.   The fractional EDG method

The EDG method was first introduced by Abdullah ^[62] as an efficient and reliable Poisson solver. The success of the EDG technique in handling the Poisson equation has initiated interest in its application to solving various types of differential equations under the frame of integer-order and CO fractional derivatives. Here, we extend the formulation of the EDG method for solving the two-dimensional VO-TFDE (1.1). In this endeavor, another discretization scheme derived by rorating the standard grid 45 clockwise is necessitated to approximate Eq (1.1). Thus, the rotated C-N finite difference approximation can be written as follows,

By dropping the local error terms and replacing $u_{i, j}^k$ with its approximation $\mathfrak{U}_{i, j}^k$, the following expression is obtained,

where $\mathfrak{C_1} = K_1/(4(\Delta x)^2)$ and $\mathfrak{C_2} = K_2/(4(\Delta y)^2)$.

Consider the location points ($i, j$), ($i+1, j+1$), ($i+1, j$) and ($i, j+1$). The application of Eq (3.2) to any group of these four points will lead to the ($4\times4$) system of equation represented in the following matrix form (see ^[62] for extra details),

where

and

The last ($4\times4$) system of equations can be decoupled into two ($2\times2$) systems whose explicit forms are as follows,

and

where

Figure 1 sheds the light on the grid shape of the EDG method. From the figure, the grid points of the discretized solution domain are branched into two distinct types, namely, black filled diamonds and white filled diamonds. One can easily verify that the implementation of Eq (3.4) involves only points of type $\blacklozenge$, while the implementation of Eq (3.5) requires only points of type $\lozenge$. As a result, the computations of Eqs (3.4) and (3.5) can be performed independently.

By keeping these observations into consideration, the EDG method can be completed in two steps. First, the solution values on one type of points are evaluated iteratively by utilizing either Eq (3.4) or Eq (3.5). Second, after convergence is achieved in the first step, the solution values on the residual grid points are computed directly once using Eq (2.5). In this way, the computational cost is reduced effectively as only half of the grid points are involved in the iterative process, which ultimately reduces the execution time of the whole computation process. The application of the EDG method for solving the VO-TFDE will result in a linear system of algebraic equations. Following a similar manner in section two, it can be easily verified that the coefficient matrix is a strictly diagonally dominant matrix and, hence, the resulting system is uniquely solvable.

4.   Stability analysis

In this section, we analyze the stability of the proposed finite difference approximations via the von-Neumann method. The stability of numerical schemes is one of the most important issues in the realm of numerical analysis that must be addressed. We start by introducing the following lemma that will be utilized in the stability analysis,

Lemma 1. The coefficients $\mathcal{W}_s^{i, j, k}$, $0\le i\le M_1-1$, $0\le j\le M_2-1$ defined in Eqs (2.5) and (3.2) satisfy the following properties ^[61],

(i) $\mathcal{W}_s^{i, j, k}\ge\mathcal{W}_{s+1}^{i, j, k}, \ 1\le s\le k = 1, 2, \dots, N-1.$

(ii) $\sum_{s = 1}^{k-1}\left(\mathcal{W}_{k-s}^{i, j, k}-\mathcal{W}_{k-s+1}^{i, j, k}\right) = \mathcal{W}_1^{i, j, k}-\mathcal{W}_k^{i, j, k}$.

Let $\hat{\mathfrak{U}}_{i, j}^k$ and $\bar{\mathfrak{U}}_{i, j}^k$ represent the numerical solutions of Eqs (2.5) and (3.2), respectively, then the errors can be defined as:

Next, by setting (4.1) into (2.5), we obtain

Similarly, replacing (4.2) into (3.2), we have

Then, $\xi^k(x, y)$ and $\varepsilon^k(x, y)$ can be expanded in the Fourier series as

where $I = \sqrt{-1}$ and $\eta^k$ and $\mu^k$ are given by

Based on the definition of $l^2$ norm and the Pareval equality, we obtain

Suppose that the solutions of Eqs (4.3) and (4.4) can be written as

in which $\gamma_1 = 2\pi Z_1/L$ and $\gamma_2 = 2\pi Z_2/L$.

Substituting $\xi_{i, j}^k = \eta^k e^{I(\gamma_1i\Delta x+\gamma_2 j\Delta y)}$ into Eq (16), we get

where

Next, we prove the following result.

Lemma 2. If $\eta^k\ (1\le k\le N)$ satisfies Eq (4.10) and $2\ge 3^{1-\alpha^{i, j, k}}$, then, $\lvert\eta^k\lvert\le\lvert\eta^0\lvert$.

Proof. We utilize mathematical induction to complete the proof. First, let $k = 0$, and we get

Since $\mathfrak{A}^{i, j, k}$, $\rho_1$ and $\rho_2$ are nonnegative numbers, then,

Now, suppose that $\left| {{\eta^{m+1}}} \right|\le\left| {{\eta^0}} \right|$, $m = 0, 1, 2, \dots, k-1$.We show it holds for $m = k$. Using Eq (4.10) and Lemma 1, we get

If $\mathfrak{A}^{i, j, k}-\rho_1-\rho_2-\sigma^{i, j, k}\mathcal{W}_1^{i, j, k} > 0$, then

If $\mathfrak{A}^{i, j, k}-\rho_1-\rho_2-\sigma^{i, j, k}\mathcal{W}_1^{i, j, k} < 0$, then

Here,

□

Now, we prove the next theorem.

Theorem 1. If $2\ge 3^{1-\alpha^{i, j, k}}$, then the finite difference scheme defined by Eq (2.5) is stable.

Proof. From Lemma 2 and the Parseval equality, we have

from which we get,

□

Setting $\varepsilon_{i, j}^k = \mu^k e^{I(\gamma_1i\Delta x+\gamma_2 j\Delta y)}$ into Eq (4.4) gives

in which

Lemma 3. If $\mu^k\ (1\le k\le N)$ satisfies Eq (4.11) and $2\ge 3^{1-\alpha^{i, j, k}}$, then, $\lvert\mu^k\lvert\le\lvert\mu^0\lvert$.

Proof. First, putting $k = 0$ in Eq (4.11), we obtain

Next, assume that $\left| {{\mu^{m+1}}} \right|\le\left| {{\mu^0}} \right|$, $m = 0, 1, 2, \dots, k-1$. We prove it is true for $m = k$. With the help of Eq (4.11) and Lemma 1, we obtain

If $\mathfrak{A}^{i, j, k}-\kappa_1-\kappa_2-\sigma^{i, j, k}\mathcal{W}_1^{i, j, k} > 0$, then

If $\mathfrak{A}^{i, j, k}-\kappa_1-\kappa_2-\sigma^{i, j, k}\mathcal{W}_1^{i, j, k} < 0$, then

Here,

□

Now, Lemma 2 is employed to prove the next theorem.

Theorem 2. If $2\ge 3^{1-\alpha^{i, j, k}}$, then, the finite difference scheme defined by Eq (3.2) is stable.

Proof. From Lemma 3 and the Parseval equality, we have

which gives the desired result

□

5.   Convergence analysis

This section is devoted to analyze the convergence of the proposed methods. Let $\zeta_{i, j}^k$ be the difference between the exact and numerical solutions at the point $(x_i, y_j, t_k)$, then the error can be written as,

By setting Eq (5.1) into Eq (2.5), the error equation follows immediately as

where $R_{i, j}^{k+1/2}$ denotes the truncation error at the point $(x_i, y_j, t_{k+1/2})$. From now on, we suppose that $C$ is a constant that may take different values at different locations. From Eq (2.4), we have

where

Next, $\zeta^k(x, y)$ and $R^{k+1/2}(x, y)$ can be expanded in the Fourier series as

such that

Considering the definition of $l^2$ norm and the Parseval equality, we get

Suppose the solutions of Eq (5.2) are as follows:

Substitute Eq (5.6) into Eq (5.2) and simplify to emerge with the following expression,

where $\rho_1$ and $\rho_2$ are as defined in Eq (4.10).

Lemma 4. If $\psi^k \; (1\le k\le N)$ satisfies Eq (5.7), then, we have

Proof. From Eq (5.5), there is a positive constant $C$ such that

We complete the proof by mathematical induction. First, setting $k = 0$ in Eq (5.7) and noticing that $\psi^0 = 0$, we obtain

Next, assume that

Putting $m = k$ in Eq (5.7) and using Lemma 1, we get

If $\mathfrak{A}^{i, j, k}-\rho_1-\rho_2-\sigma^{i, j, k}\mathcal{W}_1^{i, j, k} > 0$, then

If $\mathfrak{A}^{i, j, k}-\rho_1-\rho_2-\sigma^{i, j, k}\mathcal{W}_1^{i, j, k} < 0$, then

Here,

Therefore, under the above condition we have

□

We utilize Lemma 4 to prove the next result.

Theorem 3. If $2\ge 3^{1-\alpha^{i, j, k}}$, then the finite difference scheme defined by Eq (2.5) is $l^2$ convergent with convergence order $O(\tau+(\Delta x)^2+(\Delta y)^2)$.

Proof. With the help of Lemma 4 and Eqs (5.4) and (5.5), we obtain

Taking Eq (5.3) into consideration leads to

where $C^* = CT$ as $(k+1)\tau\le T$. □

Theorem 4. If $2\ge 3^{1-\alpha^{i, j, k}}$, then the finite difference scheme defined by Eq (3.2) is $l^2$ convergent with convergence order $O(\tau+(\Delta x)^2+(\Delta y)^2)$.

Proof. The proof can be established in analogous way to Theorem 3. □

Theorems 5.1 and 5.2 indicate the convergence of the numerical schemes (2.5) and (3.2), respectively. This means that as the spatial and temporal step sizes approach zero, the numerical solutions approach to the exact solution of the considered problem. These considerations will be verified in the next section.

6.   Numerical experiments and discussion of results

6.1.   Numerical experiments

After the construction and analysis of the C-N and EDG solution methods in the previous sections, we turn our attention to simulating the two-dimensional VO-TFDE (1.1) using the aforementioned numerical schemes. In this endeavor, three numerical examples are chosen to illustrate the capability, accuracy and efficiency of the presented methods. The exact analytical solution of each example is known and is utilized for comparison with the evaluated numerical solutions. Unless otherwise specified, the solution domain is set to be $[0, 1]\times[0, 1]\times(0, 1]$. Furthermore, we assume that $\Delta x = \Delta y = h$ throughout numerical experiments. All computations are carried out using the Julia Programming Language and run on a PC with the configuration: Intel(R) Core(TM) i7-8550U CPU 2.00 GHz 8 GB RAM and the Windows 10 (64-bit) operating system. Matlab is also employed for generating some figures. The computational cost and the accuracy of the proposed methods will be reported based on the criteria of CPU time (in seconds) and the maximum absolute error ($MAE$), respectively. The numerical outcomes of these criteria are attained for different values of mesh sizes and various VO functions. Finally, we include some graphical representations to facilitate the discussion of our results.

The $MAE$ between the exact solution and the numerical solution is given by

To compute the convergence order ($CO$), we use

Example 1. Consider the following two-dimensional VO-TFDE:

with the initial and boundary conditions

and the exact analytical solution is given by

Tables 1 gives the CPU times and the numerical errors of the proposed methods for a variety of mesh sizes (1/6, 1/12, 1/18, 1/24) with fixed $\tau = 0.01$ and various values of $\alpha(x, y, t)$. In view of this table, it can be seen that the C-N and EDG methods produce accurate results in solving Example 1. Furthermore, as the mesh size increases, the observed maximum errors for all chosen VO derivatives $\alpha(x, y, t)$ decrease. Figure 2 sketches the CPU time of the C-N and EDG methods according to different sizes of space steps with $\alpha(x, y, t) = \frac{10-(xyt)^3}{80}$ and $\alpha(x, y, t) = \frac{9-x^5+t^3}{120}$. We can see that the EDG method significantly outperforms the C-N method in terms of CPU timing. The time reduction percentage of the EDG method compared to the C-N method is (52.74-79.44)%. For instance, Table 1 shows that we need to take $h = 1/24$ in order to obtain the maximum error of 2.4740E-4 using the C-N method with $\alpha(x, y, t) = \frac{10-(xyt)^3}{80}$, which costs 44.8959 seconds. On the other hand, the EDG method with $h = 1/24$ achieves a maximum error of 3.4746E-4 and the corresponding computing time is only 11.4655 seconds. These results and comparisons justify our effort to establish the EDG method. Figure 3 shows a comparison of the numerical solutions of the proposed methods and the exact solution of Example 1 for $h = 1/24$, $\tau = 0.01$, $x = 0.5$, $\alpha(x, y, t) = \frac{11-\cos(xt)^2}{11}$ and $T = 0.25, 0.5, 0.75, 1$. One can observe from this figure that the behavior of the numerical solutions is the same as that of the exact solution. Figure 4 displays 3D surface plots of the numerical solutions and the maximum errors for the C-N and EDG methods at $h = 1/24$, $\tau = 0.01$ and $\alpha(x, y, t) = \frac{15+\sin(yt)^5}{16}$. It follows that the numerical solutions of the C-N and EDG methods are in good agreement with the exact solutions, which illustrate the accuracy and effectiveness of the proposed methods.

Example 2. In this example, we consider the following two-dimensional VO-TFDE:

subject to the initial and boundary conditions

with the exact analytical solution written as

We now utilize the C-N and EDG methods to solve the given problem in Example 2. The CPU times and the maximum errors for several mesh sizes (1/4, 1/10, 1/16, 1/22) with fixed $\tau = 0.01$ and various values of $\alpha(x, y, t)$ are shown in Table 2. Figure 5 depicts the CPU times of the proposed methods against different values of spatial step sizes for $\alpha(x, y, t) = \frac{e^{xyt}+\cos(xyt)}{30}$ and $\alpha(x, y, t) = \frac{5+(xy)^3-(xt)^4}{50}$. From the data in these tables and figures, it's seen that the C-N and EDG methods generate comparable results in terms of accuracy, but the latter method results in better simulations in terms of CPU timing. From Table 2, the improvement percentage in CPU time of the EDG method compared to the C-N method is (40.36-83.56)%. This example clearly lends itself to a comparison similar to that of Example 1. The numerical solutions of the proposed methods along with the exact solution of Example 2 is portrayed in Figure 6 for $h = 1/22$, $\tau = 0.01$, $x = 0.5$, $\alpha(x, y, t) = \frac{16-e^{xyt}}{17}$ and $T = 0.25, 0.5, 0.75, 1$. Figure 7 highlights the 3D surface plots of the numerical solutions and the maximum errors of the C-N and EDG methods at $h = 1/22$, $\tau = 0.01$ and $\alpha(x, y, t) = \frac{3+(xy)^2-(xt)^3}{30}$. These figures demonstrate the effectiveness of the proposed methods as the numerical solutions are found to match well with the exact analytical solution.

Example 3. Here, the following two-dimensional VO-TFDE is considered:

with the initial and boundary conditions

and the exact analytical solution of the following form:

Table 3 highlights the numerical results of solving Example 3 using CN and EDG methods with $T = 1$, $\tau = 0.01$, $h^{-1} = 12, 24$ and different VO functions. From the data in this table, it can be seen that the numerical solutions converge to the exact solution where maximum errors decrease as mesh sizes increase. A reduction of CPU timing and, hence, computational burden is also noted in using the EDG method over the CN method for solving the considered problem. This can be beneficial for capturing the long-term behavior of physical phenomena that is governed by VO fractional models. From the previous discussion, it follows that the proposed methods are accurate and the EDG method is superior to the C-N method in terms of CPU timing. Finally, Table 4 highlights the computational orders of the proposed methods for solving Examples 1–3. From this table, the second order convergence of the proposed methods can be anticipated, which is in line with the theoretical considerations.

6.2.   Discussion of results

In the previous subsection, three numerical examples were solved using the proposed methods, and the results were illustrated in Tables 1–3 along with Figures 2–7. In view of the said tables, it can be seen that the numerical solutions of the proposed methods approach the exact solution as the spatial step size decreases. This is in line with our theoretical analysis. In addition, the EDG method was observed to be computationally efficient compared to the C-N method in solving the test problems, where the former required less CPU time compared to the latter. Such reduction in computational burden is attributed to the smaller number of mesh points comprised in the iterative process as discussed earlier. The experimental results are in line with the considerations of previous sections and demonstrate the feasibility, accuracy and efficiency of the proposed methods.

7.   Conclusions

The VO fractional diffusion equation and its numerical solution are of practical importance for describing complex diffusion processes in which the system memory changes as a function of time and/or space variables. In this paper, the two-dimensional VO-TFDE (1.1) has been considered to account for its numerical solution. A C-N difference scheme along with an EDG method constructed based on a skewed finite difference approximation on the rotated grid have been developed for solving the VO fractional problem. The stability and convergence analyses have been investigated in detail. Numerical test problems were provided to illustrate the capability and performance of the proposed methods.

Based on the numerical simulations and the discussion of experimental results, some conclusions are drawn.

First, both C-N and EDG methods perform well for the VO-TFDE as their numerical solutions are found to be in good agreement with the exact analytical solutions. See Figures 3, 4, 6 and 7.

Second, the proposed methods are stable and convergent, and the maximum absolute errors are observed to decrease as the mesh sizes increase. See Tables 1 and 2.

Third, the EDG method reduces the computational burden of solving the VO-TFDE and results in faster simulations compared to the C-N difference method. See Figures 2 and 5.

Finally, the current work has successfully developed the C-N and EDG numerical methods to solve the fractional diffusion equation involving the VO time derivative defined in section one. The extension of the proposed methods to other classes of linear and nonlinear fractional differential equations with time-invariant type VO derivative ^[63] is an interesting line of future research. The construction of the proposed methods based on higher-order difference approximations for enhancing the accuracy is another direction for further study.

Use of AI tools declaration

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

Acknowledgments

The authors would like to acknowledge the support provided by the Deanship of Research Oversight and Coordination (DROC) at King Fahd University of Petroleum & Minerals (KFUPM) for funding this work through project (No. EC221011). We are thankful for the anonymous reviewers for their valuable comments and suggestions.

Conflict of interest

The authors declare no conflicts of interest.