On Numerical Simulations of Variable-Order Fractional Cable Equation Arising in Neuronal Dynamics

Abstract

In recent years, various complex systems and real-world phenomena have been shown to include memory and hereditary properties that change with respect to time, space, or other variables. Consequently, fractional partial differential equations containing variable-order fractional operators have been extensively resorted for modeling such phenomena accurately. In this paper, we consider the two-dimensional fractional cable equation with the Caputo variable-order fractional derivative in the time direction, which is preferable for describing neuronal dynamics in biological systems. A point-wise scheme, namely, the Crank–Nicolson finite difference method, along with a group-wise scheme referred to as the explicit decoupled group method are proposed to solve the problem under consideration. The stability and convergence analyses of the numerical schemes are provided with complete details. To demonstrate the validity of the proposed methods, numerical simulations with results represented in tabular and graphical forms are given. A quantitative analysis based on the CPU timing, iteration counting, and maximum absolute error indicates that the explicit decoupled group method is more efficient than the Crank–Nicolson finite difference scheme for solving the variable-order fractional equation.

1. Introduction

The concepts of integral and derivative of the integer order form the backbone of so-called classical calculus. Many quantities such as displacement, velocity, acceleration, jerk, jounce, area and volume can be described precisely based on classical calculus. However, it was discovered that systems of many real-world phenomena cannot be well modeled using classical calculus tools, particularly systems that depict memory effects. This means the case in which the current state of the system depends on all previous ones. For instance, the classical diffusion model is an efficient tool for describing transport processes where the mean square displacement is linear in time $<X^2>\sim Kt$, where K is the diffusion constant. However, in complex backgrounds such as porous media and biological systems, it was found that the mean square displacement is no longer linear in time $<X^2>\sim K{t}^{\gamma }$, where $\gamma$ is the diffusion exponent. This leads to a new phenomenon known as the anomalous transport process, where the diffusion becomes either slower or faster than the classical model. Obviously, such a process cannot be described adequately by the classical diffusion model. Fractional calculus is a generalization of classical calculus, where the differential and integral operators are allowed to take non-integer values [1]. Due to their non-local property, fractional differential operators serve as powerful tools for modeling complex phenomena that cannot be characterized by classical calculus. In line with that, fractional differential equations have been utilized extensively for modeling seemingly diverse practical problems in widespread fields of science, finance, mechanics, medicine and engineering [2,3]. There are numerous types of fractional-order derivatives that can be found in the literature. The definitions of Riemann–Lioville, Caputo, Grünwald–Letnikov, Riesz, and Hadamard are among the most popular derivatives of fractional calculus. These traditional derivatives are also referred to as power-law type fractional operators. On the other hand, several other kinds of fractional operators that involve Mittag–Leffler and exponential kernel types have been developed recently to overcome the initial singularities of the problems associated with the traditional fractional derivatives. The interested reader may refer to [4] for a full description of the fractional derivatives. Fractional calculus is one of the hottest topics of applied mathematics, where a diverse range of real-world problems can be effectively modeled in the form of fractional differential equations. In 2017, Matlob and Jamali [5] argued that there are no certain rules for selecting the type of fractional derivative suitable for modeling. However, in 2020, Tarasov and Tarasova [6] established a correspondence between the properties of the kernel of the fractional operator and the physical phenomenon under consideration. In the following discussion, some recent applications of fractional-order derivatives are briefly reviewed. In [7], Din et al. investigated the behavior of the climate change phenomenon under the frame of a fractional-order model. The authors employed the Caputo-type fractional derivative and reported that the fractional model generates better results than the integer-order mathematical model. Reyaz et al. [8] utilized the Caputo–Fabrizio fractional derivative to analyze the concentration, temperature, and velocity of a fluid flow with chemical reaction effects and thermal radiation. The authors noticed that the transition between the steady state and the unsteady state of the fluid can be controlled by the fractional-order gamma. In [9], Yang et al. researched the dynamics of the hepatitis B virus (HBV) via a fractional-order biliogical model defined in terms of a Caputo-type fractional derivative. Another fractional epidemic model under the Caputo fractional derivative was presented by Sharif Ullah et al. [10] to study the behavior of the COVID-19 outbreak. Rashid et al. [11] explored childhood diseases and their complications via an SIR model involving the Atangana–Baleanu fractional derivative. In [12], a study on the dynamics of fungicide application was proposed via a fractional mathematical model with the Caputo–Fabrizio operator. In [13], Lauria et al. established a novel intra-hour photovoltaic power forecasting method based on the Caputo-type fractional derivative. In [14], Anwar et al. considered two mathematical fractional models to account for the flow patterns and thermal behavior of a hybrid nanofluid. A comprehensive treatment of the theory, applications, and simulations of fractional derivative-based models can be found in [15].

Roughly speaking, the analytical treatment of fractional derivative-based models that are expressed in terms of fractional partial differential equations (FPDEs) is not an easy task [16]. Therefore, numerical methods have attracted much attention for dealing with FPDEs. The list of numerical methods that are available in the literature includes the finite difference method, finite element method [17,18], finite volume method, spectral method, collocation method, reproducing kernel method, mesh-free method, domain decomposition method, etc.

However, it has been verified that FPDEs with fixed fractional orders are not adequate for describing complex phenomena that exhibit variable memory with respect to time and/or space variables [19]. As a result, a new field of variable-order (VO) fractional calculus has emerged. It was Samko and Ross [20] who paved the way for this interesting field by generalizing the Riemann–Liouville and Marchaud fractional derivatives to their VO sense in 1993. Thereafter, other types of VO fractional differential operators were suggested by Lorenzo, Hartley and Coimbra [21,22,23]. Today, various definitions of VO fractional derivatives with specific meanings are available to handle real-world problems that depict systems with varying memory. The merits of using VO fractional derivatives instead of their corresponding constant-order (CO) fractional derivatives are illustrated in [24]. Another study revealed that it is much easier to describe the physical meaning of VO fractional operators [25]. Consequently, VO FPDEs have attracted the attention of numerous researchers and scholars as accurate models for describing a large variety of phenomena in various branches of science and engineering, such as mechanics, viscoelasticity, anomalous diffusion, wave propagation, control theory, ecology, and many others. Similar to the CO FPDEs, it is difficult to solve VO FPDEs analytically, and numerical techniques are very often resorted to, for instance, see [26,27,28,29]. A profound discussion of the definitions, applications, and numerical simulations of the VO fractional operators can be seen in the insightful review papers [30,31].

In this work, we consider the numerical solution of the non-homogeneous initial-boundary value problem of the two-dimensional variable-order fractional cable equation (VO FCE) in the form,

$$\begin{array}{cc}\hfill {}_0{}^{C}D_t^{\gamma (x,y,t)}u(x,y,t)& ={\mathcal{A}}_1\frac{{\mathsf{\partial }}^{2}u(x,y,t)}{\mathsf{\partial }{x}^2}+{\mathcal{A}}_2\frac{{\mathsf{\partial }}^{2}u(x,y,t)}{\mathsf{\partial }{y}^2}-\mu u(x,y,t)\hfill \\ \hfill & +f(x,y,t),(x,y,t)\in \Omega \times [0,T],\hfill \end{array}$$

(1)

with the initial-boundary conditions,

$$u(x,y,0)={\theta }_1(x,y),(x,y)\in \Omega$$

(2)

$$\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}u(0,y,t)={\theta }_2(y,t),u(L_1,y,t)={\theta }_3(y,t),(x,y,t)\in [0,T]\times \mathsf{\partial }\Omega ,$$

$$\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}u(x,0,t)={\theta }_4(x,t),u(x,L_2,t)={\theta }_5(x,t).$$

(3)

Here $\Omega =[0,L_1]\times [0,L_2]$ is the domain, and $\mathsf{\partial }\Omega$ is its boundary, ${\mathcal{A}}_1$, ${\mathcal{A}}_2$ and $\mu$ are positive constants. ${\theta }_i,i=1,2,\cdots ,5$ are known functions, and $u(x,y,t)$ is the unknown function. $0\le \gamma (x,y,t)\le 1$ and ${}_0{}^{C}D_t^{\gamma (x,y,t)}u(x,y,t)$ is the VO Caputo fractional derivative of u defined as,

$${}_0{}^{C}D_t^{\gamma (x,y,t)}u(x,y,t)=\left\{\begin{array}{c}\frac{1}{\Gamma (1-\gamma (x,y,t\left)\right)}{\int }_0^t{(t-\psi )}^{-\gamma (x,y,t)}\frac{\mathsf{\partial }u(x,y,\psi )}{\mathsf{\partial }\psi }d\psi ,0<\gamma (x,y,t)<1,\hfill \\ \frac{\mathsf{\partial }u(x,y,t)}{\mathsf{\partial }t},\gamma (x,y,t)=1.\hfill \end{array}\right.$$

The FCE can be viewed as a fundamental biological model that accounts for the voltage difference between the cell membrane and neurons. It was shown that the electrodiffusion of ions in spiny neuronal dendrites follows an anomalous pattern that cannot be captured by the classical cable equation [32]. In 2008, Henry and Langlands [33] introduced the FCE for the first time to study the anomalous diffusion phenomenon in the nerve cell. Since then, FCE has been employed to report on the behavior of different phenomena occurring in several fields, such as neuronal dynamics, control theory, the heterogeneous nature of neuronal tissue, and viscoelastic materials (see [34] and references cited therein). Consequently, the solution of FCE is of practical importance in different fields of application.

In the past few years, several numerical researchers have contributed significantly to solving FCEs since exact analytical solutions are mostly difficult to obtain. Here, we present the recent numerical treatments of the FCE. In [35], Bhrawy and Zaky proposed a spectral collocation method for solving the one-dimensional and two-dimensional VO FCE. In [36], Irandoust-Pakchin et al. generalized the application of the Chebyshev cardinal functions method for solving the one-dimensional VO FCE. In [37], Nagy and Sweilam researched a Crank–Nicolson numerical scheme with stability analysis to solve the one-dimensional VO FCE. In [38], Liu et al. presented finite difference/element approximations, which are formulated utilizing some second-order difference schemes in time and the Galerkin finite element method in space for solving one-dimensional and two-dimensional CO FCE. Stability and convergence analyses were also discussed. In [39], Sweilam and AL-Mekhlafi developed a non-standard compact finite difference scheme to solve the two-dimensional CO FCE. The fractional derivative is defined in the Atangana–Baleanu–Caputo sense, and the truncation error of their method is analyzed. An extended cubic B-Spline collocation method with stability analysis was introduced in [40] to solve the one-dimensional CO FCE. Salama and Ali [41] proposed a fast hybrid method based on a combination of the Laplace transform technique and finite difference approximations to account for the numerical solution of the two-dimensional CO FCE. A meshless numerical scheme with rigorous error analysis for the two-dimensional CO FCE was reported in [42]. Oruç [43] suggested a local hybrid kernel meshless method to deal with the two-dimensional CO FCE. Zheng et al. [44] utilized a finite difference scheme in the temporal direction and a Legendre spectral method in the spatial direction to solve the two-dimensional distributed-order FCE. The theoretical stability and convergence were established in their work. Kumar and Baleanu [45] derived an operational matrix-based numerical method for the solution of the two-dimensional CO FCE. Li and Rui [46] scrutinized an unconditionally stable block-centered finite difference method for solving the two-dimensional CO FCE. Mohebbi and Saffarian [47] applied second-order difference schemes for time discretization and a meshless approach for space discretization in solving the two-dimensional VO FCE. The theoretical stability and convergence have not been investigated. Yin et al. [48] applied the FBT-$\theta$ and FBN-$\theta$ methods proposed in [49] to solve the one-dimensional and two-dimensional CO FCE. The authors gave a detailed numerical analysis of stability and convergence. Sweilam et al. [50] analyzed the two-dimensional FCE and fractional reaction sub-diffusion equation by the weighted average finite difference method. The Sinc–Bernoulli collocation method was considered by Moshtaghi and Saadatmandi [51] to study the one-dimensional CO FCE. Mittal et al. [34] introduced a time–space Jacobi pseudospectral method for the numerical solution of the two-dimensional FCE. Stability, uniqueness, and error analysis were also given in their work. Li and Li [52] solved the two-dimensional CO FCE using a meshless finite point method. Other research work concerned with the numerical analysis of the CO FCE can be seen in [53,54,55,56].

Group-wise iterative methods, usually called explicit group methods, can be thought of as efficient alternatives to their corresponding point-wise iterative methods for solving various types of initial-boundary value problems. In summary, explicit group methods have some salient advantages, such as stability, being easy to implement, forming sparse algebraic systems, accelerating the rate of convergence, reducing arithmetic computations per iteration, and extension to multi-dimensional problems. In [57,58,59,60,61,62], different types of explicit group methods were developed based on standard and skewed difference schemes for solving a variety of CO FPDEs.

Based on the previous discussion, much research has been reported on the CO FCE, whereas numerical works on FCE involving VO fractional derivatives are very few and still far from adequate. In addition, to the best of our knowledge, VO FCE solved by utilizing explicit group methods has not emerged in the literature. Motivated by this, the aim of the current paper is to present a group-wise iterative method, namely the explicit decoupled group method (EDGM) for solving the two-dimensional VO FCE (1). The said method inherits the advantages of the family of explicit group methods for being accurate, stable, efficient, and extendable for multi-dimensional problems. The EDGM is established using Taylor series expansion on a skewed grid obtained by rotating the x-y axes 45° clockwise. Furthermore, we derive a point-wise iterative method, namely, the Crank–Nicolson finite difference method (CN–FDM), as a reference solution algorithm for the VO FCE. The aspects of stability, convergence, and computational complexity in terms of computing time and iteration count are given. The tabulated and graphical results extracted from numerical simulations show that the presented methods have comparable accuracy, while the EDGM results in faster simulations with lower computational complexity compared to the CN–FDM.

The plan of our paper is as follows. We construct the CN–FDM in Section 2 followed by the EDGM in Section 3. Then, Section 4 and Section 5 are devoted to the stability and convergence analyses of the proposed methods, respectively. Some numerical simulations to validate the accuracy and efficiency of the CN–FDM and EDGM are given in Section 6. Finally, the current work ends with a brief conclusion in Section 7.

2. Formulation of the Crank–Nicolson Finite Difference Method (CN–FDM)

In this section, each of the VO time fractional derivatives together with the spatial differential operators in Equation (1) are discretized to construct the CN–FDM. For this purpose, the following notations are defined as,

$$\begin{array}{c}\hfill t_k=k\tau ,0\le k\le N,x_i=i{h}_x,0\le i\le M_x,y_j=j{h}_y,0\le j\le M_y.\end{array}$$

Here, N, $M_x$ and $M_y$ are positive integers and denote the number of equidistant partitions in time and space directions, respectively. $h_x=L_1/M_x$, $h_y=L_2/M_y$ and $\tau =T/N$ are the spatial and temporal step sizes, respectively. In addition, define

$$\begin{array}{c}\hfill u_{i,j}^k=u(x_i,y_j,t_k),0\le i\le M_x,0\le j\le M_y,0\le k\le N,\end{array}$$

as the solution value located at the i-th and j-th coordinates and k-th time level. Obviously, $u^0={\theta }_1(x,y)$ represents the given initial condition, while $u^k$ refers to the unknown solution values for $1\le k\le N$.

The approximate formula for the Caputo VO time fractional derivative ($0\le \gamma (x,y,t)\le 1$) at the time node $t_{k+1/2}$ can be obtained by utilizing the following discretization scheme [63]:

$$\begin{array}{cc}\hfill & {}_0{}^{C}D_t^{i,j,k+1/2}u(x_i,y_j,t_{k+1/2})\hfill \\ \hfill & ={\sigma }_1\left[{\mathcal{H}}_1^{i,j,k}{u}_{i,j}^k-\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right)u_{i,j}^s-{\mathcal{H}}_k^{i,j,k}{u}_{i,j}^0+\frac{u_{i,j}^{k+1}-u_{i,j}^k}{2^{1-{\gamma }_{i,j,k+1/2}}}\right]\hfill \\ \hfill & +{\delta }^{k+1/2},\hfill \end{array}$$

(4)

where

$$\begin{array}{cc}\hfill {\sigma }_1& =\frac{1}{\Gamma (2-\gamma (x_i,y_j,t_{k+1/2})){\tau }^{\gamma (x_i,y_j,t_{k+1/2})}},\hfill \\ \hfill {\mathcal{H}}_s^{i,j,k}& ={(s+1/2)}^{1-\gamma (x_i,y_j,t_{k+1/2})}-{(s-1/2)}^{1-\gamma (x_i,y_j,t_{k+1/2})},\hfill \end{array}$$

and $r^{k+1/2}$ is the local truncation error bounded by,

$$\left|{\delta }^{k+1/2}\right|\le C\tau .$$

The right-hand side of Equation (1) can be discretized by applying central difference approximations at the time levels k and $k+1$, and taking the average as follows:

$$\begin{array}{cc}\hfill & {\mathcal{A}}_1\frac{{\mathsf{\partial }}^{2}u(x,y,t)}{\mathsf{\partial }{x}^2}+{\mathcal{A}}_2\frac{{\mathsf{\partial }}^{2}u(x,y,t)}{\mathsf{\partial }{y}^2}-\mu u(x,y,t)+f(x,y,t)\hfill \\ \hfill & =\frac{{\mathcal{A}}_1}{2}\left[\frac{u_{i+1,j}^{k+1}-2u_{i,j}^{k+1}+u_{i-1,j}^{k+1}}{h_x^2}+\frac{u_{i+1,j}^k-2u_{i,j}^k+u_{i-1,j}^k}{h_x^2}\right]\hfill \\ \hfill & +\frac{{\mathcal{A}}_2}{2}\left[\frac{u_{i,j+1}^{k+1}-2u_{i,j}^{k+1}+u_{i,j-1}^{k+1}}{h_y^2}+\frac{u_{i,j+1}^k-2u_{i,j}^k+u_{i,j-1}^k}{h_y^2}\right]\hfill \\ \hfill & -\mu \left(\frac{u_{i,j}^{k+1}+u_{i,j}^k}{2}\right)+f_{i,j}^{k+1/2}+O({\tau }^2+h_x^2+h_y^2).\hfill \end{array}$$

(5)

By substituting Equations (4) and (5) into Equation (1), the following expression is obtained:

$$\begin{array}{cc}\hfill & {\sigma }_1\left[{\mathcal{H}}_1^{i,j,k}{u}_{i,j}^k-\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right)u_{i,j}^s-{\mathcal{H}}_k^{i,j,k}{u}_{i,j}^0+\frac{u_{i,j}^{k+1}-u_{i,j}^k}{2^{1-{\gamma }_{i,j,k+1/2}}}\right]\hfill \\ \hfill & =\frac{{\mathcal{A}}_1}{2}\left[\frac{u_{i+1,j}^{k+1}-2u_{i,j}^{k+1}+u_{i-1,j}^{k+1}}{h_x^2}+\frac{u_{i+1,j}^k-2u_{i,j}^k+u_{i-1,j}^k}{h_x^2}\right]\hfill \\ \hfill & +\frac{{\mathcal{A}}_2}{2}\left[\frac{u_{i,j+1}^{k+1}-2u_{i,j}^{k+1}+u_{i,j-1}^{k+1}}{h_y^2}+\frac{u_{i,j+1}^k-2u_{i,j}^k+u_{i,j-1}^k}{h_y^2}\right]\hfill \\ \hfill & -\mu \left(\frac{u_{i,j}^{k+1}+u_{i,j}^k}{2}\right)+f_{i,j}^{k+1/2}+O(\tau +h_x^2+h_y^2).\hfill \end{array}$$

(6)

Omitting the error terms and replacing $u_{i,j}^k$ with its approximation $U_{i,j}^k$, and upon simplification, the CN–FDM for solving the Dirichlet-type boundary value problem (1)–(3) is obtained as follows:

$$\left\{\begin{array}{cc}\hfill & ({\mathcal{W}}_1+0.5\mu +2{\mathcal{P}}_1+2{\mathcal{P}}_2)U_{i,j}^{k+1}={\mathcal{P}}_1(U_{i+1,j}^{k+1}+U_{i-1,j}^{k+1}+U_{i+1,j}^k+U_{i-1,j}^k)\hfill \\ \hfill & +{\mathcal{P}}_2(U_{i,j+1}^{k+1}+U_{i,j-1}^{k+1}+U_{i,j+1}^k+U_{i,j-1}^k)+({\mathcal{W}}_1-0.5\mu -{\sigma }_1{\mathcal{H}}_1^{i,j,k}-2{\mathcal{P}}_1-2{\mathcal{P}}_2)U_{i,j}^k\hfill \\ \hfill & +{\sigma }_1\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right)U_{i,j}^s+{\sigma }_1{\mathcal{H}}_k^{i,j,k}{U}_{i,j}^0+f_{i,j}^{k+1/2},\hfill \\ \hfill & 1\le i\le M_x-1,1\le j\le M_y-1,0\le k\le N-1,\hfill \\ \hfill & {\mathcal{W}}_1=\frac{{\sigma }_1}{2^{1-{\gamma }_{i,j,k+1/2}}},{\mathcal{P}}_1=\frac{{\mathcal{A}}_1}{2h_x^2},{\mathcal{P}}_2=\frac{{\mathcal{A}}_2}{2h_y^2}.\hfill \\ \hfill & \mathrm{The}\mathrm{initial}\mathrm{and}\mathrm{boundary}\mathrm{conditions}\mathrm{are}\hfill \\ \hfill & U_{i,j}^0={\theta }_1(x_i,y_j),0\le i\le M_x,0\le j\le M_y,\hfill \\ \hfill & U_{i,j}^k{|}_{\mathsf{\partial }\Omega }=\theta (x_i,y_j,t_k),0\le k\le N.\hfill \end{array}\right.$$

(7)

Let

$$\begin{array}{cc}\hfill & U^k={[U_{1,1},U_{1,2},\cdots ,U_{1,M_y-1},U_{2,1},U_{2,2},\cdots ,U_{2,M_y-1},\cdots ,U_{M_x-1,1},U_{M_x-1,2},\cdots ,U_{M_x-1,M_y-1}]}^T,\hfill \\ \hfill & f^k={[f_{1,1},f_{1,2},\cdots ,f_{1,M_y-1},f_{2,1},f_{2,2},\cdots ,f_{2,M_y-1},\cdots ,f_{M_x-1,1},f_{M_x-1,2},\cdots ,f_{M_x-1,M_y-1}]}^T.\hfill \end{array}$$

The fully discrete scheme can be expressed in matrix form as,

$$\left\{\begin{array}{c}A{U}^1=B{U}^0+f^{1/2},k=0,\hfill \\ A{U}^{k+1}=B{U}^k+{\sigma }_1\sum _{m=1}^{k-1}\left({\mathcal{H}}_{k-m}{\mathcal{H}}_{k-m+1}\right)U^m+\sigma {\mathcal{H}}_k{U}^0+f^k,k\ge 1,\hfill \end{array}\right.$$

where A and B are pentadiagonal matrices defined as,

$$A=\left[\begin{array}{cccccccc}{\mathcal{S}}^{i,j,k}& \cdots & -{\mathcal{P}}_1& \cdots & -{\mathcal{P}}_2& \cdots & & 0\\ ⋮& {\mathcal{S}}^{i,j,k}& \cdots & -{\mathcal{P}}_1& \cdots & -{\mathcal{P}}_2& & ⋮\\ -{\mathcal{P}}_1& \ddots & {\mathcal{S}}^{i,j,k}& \ddots & \ddots & & \ddots & \\ ⋮& \ddots & \ddots & \ddots & & & \ddots & -{\mathcal{P}}_2\\ -{\mathcal{P}}_2& & -{\mathcal{P}}_1& \ddots & \ddots & & \ddots & ⋮\\ ⋮& -{\mathcal{P}}_2& \ddots & \ddots & & \ddots & & -{\mathcal{P}}_1\\ & & \ddots & & \ddots & & \ddots & ⋮\\ 0& \cdots & & -{\mathcal{P}}_2& \cdots & -{\mathcal{P}}_1& \cdots & {\mathcal{S}}^{i,j,k}\end{array}\right]$$

,

$$B=\left[\begin{array}{cccccccc}{\mathcal{T}}^{i,j,k}& \cdots & {\mathcal{P}}_1& \cdots & {\mathcal{P}}_2& \cdots & & 0\\ ⋮& {\mathcal{T}}^{i,j,k}& \cdots & {\mathcal{P}}_1& \cdots & {\mathcal{P}}_2& & ⋮\\ {\mathcal{P}}_1& \ddots & {\mathcal{T}}^{i,j,k}& \ddots & \ddots & & \ddots & \\ ⋮& \ddots & \ddots & \ddots & & & \ddots & {\mathcal{P}}_2\\ {\mathcal{P}}_2& & {\mathcal{P}}_1& \ddots & \ddots & & \ddots & ⋮\\ ⋮& {\mathcal{P}}_2& \ddots & \ddots & & \ddots & & {\mathcal{P}}_1\\ & & \ddots & & \ddots & & \ddots & ⋮\\ 0& \cdots & & {\mathcal{P}}_2& \cdots & {\mathcal{P}}_1& \cdots & {\mathcal{T}}^{i,j,k}\end{array}\right],$$

where

$${\mathcal{S}}^{i,j,k}={\mathcal{W}}_1+0.5\mu +2{\mathcal{P}}_1+2{\mathcal{P}}_2, {\mathcal{T}}^{i,j,k}={\mathcal{W}}_1-0.5\mu -{\sigma }_1{\mathcal{H}}_{i,j,k}^1-2{\mathcal{P}}_1-2{\mathcal{P}}_2.$$

In view of the structure of the above matrices, it can be seen that A is a strictly diagonally dominant matrix. This means that the matrix A is non-singular and the fully discrete scheme defined in Equation (7) has a unique solution. Figure 1 depicts the computational molecule of the CN–FDM with $M_x=M_y=10$. In practice, the CN–FDM is combined with the Gauss–Seidel iterative scheme to generate iterations on all ⧫ points at each time level until convergence is achieved. This point-wise iterative method terminates when the targeted time level N is reached.

In order to accelerate the rate of convergence, a new group-wise iterative method, namely the EDGM, is introduced in the next next section.

3. Formulation of the Explicit Decoupled Group Method (EDGM)

The idea of group-wise iterative methods is to generate iterations on fixed-size groups of points rather than on a single point in point-oriented iterative methods. It has been proven that the iteration matrices of group-wise iterative methods have better spectral properties than those of point-wise iterative methods, which makes the former methods computationally superior to the latter ones. In this section, we present the construction of the EDGM for solving the VO FCE (1). To this end, a new approximation scheme based on the skewed grid for the mentioned equation is derived. This will result in the following skewed CN–FDM of the form,

$$\begin{array}{cc}\hfill & {\sigma }_1\left[{\mathcal{H}}_1^{i,j,k}{u}_{i,j}^k-\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right)u_{i,j}^s-{\mathcal{H}}_k^{i,j,k}{u}_{i,j}^0+\frac{u_{i,j}^{k+1}-u_{i,j}^k}{2^{1-{\gamma }_{i,j,k+1/2}}}\right]\hfill \\ \hfill & =\frac{{\mathcal{A}}_1}{2}\left[\frac{u_{i+1,j-1}^{k+1}-2u_{i,j}^{k+1}+u_{i-1,j+1}^{k+1}}{2h_x^2}+\frac{u_{i+1,j-1}^k-2u_{i,j}^k+u_{i-1,j+1}^k}{2h_x^2}\right]\hfill \\ \hfill & +\frac{{\mathcal{A}}_2}{2}\left[\frac{u_{i+1,j+1}^{k+1}-2u_{i,j}^{k+1}+u_{i-1,j-1}^{k+1}}{2h_y^2}+\frac{u_{i+1,j+1}^k-2u_{i,j}^k+u_{i-1,j-1}^k}{2h_y^2}\right]\hfill \\ \hfill & -\mu \left(\frac{u_{i,j}^{k+1}+u_{i,j}^k}{2}\right)+f_{i,j}^{k+1/2}+O(\tau +h_x^2+h_y^2).\hfill \end{array}$$

(8)

After simplification and disregarding the error terms and utilizing $U_{i,j}^k$ as an approximation to $u_{i,j}^k$, the following expression is obtained:

$$\begin{array}{cc}\hfill & ({\mathcal{W}}_1+0.5\mu +2{\mathcal{Q}}_1+2{\mathcal{Q}}_2)U_{i,j}^{k+1}={\mathcal{Q}}_1(U_{i+1,j-1}^{k+1}+U_{i-1,j+1}^{k+1}+U_{i+1,j-1}^k\hfill \\ \hfill & +U_{i-1,j+1}^k)+{\mathcal{Q}}_2(U_{i+1,j+1}^{k+1}+U_{i-1,j-1}^{k+1}+U_{i+1,j+1}^k+U_{i-1,j-1}^k)+({\mathcal{W}}_1\hfill \\ \hfill & -0.5\mu -{\sigma }_1{\mathcal{H}}_1^{i,j,k}-2{\mathcal{Q}}_1-2{\mathcal{Q}}_2)U_{i,j}^k+{\sigma }_1\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right)U_{i,j}^s\hfill \\ \hfill & +{\sigma }_1{\mathcal{H}}_k^{i,j,k}{U}_{i,j}^0+f_{i,j}^{k+1/2},\hfill \end{array}$$

(9)

where ${\mathcal{Q}}_1=\frac{{\mathcal{A}}_1}{4h_x^2}$, ${\mathcal{Q}}_2=\frac{{\mathcal{A}}_2}{4h_y^2}$.

Next, consider the grid points of the discretized solution domain that are located at the coordinates of $(i,j)$, $(i+1,j+1)$, $(i+1,j)$, $(i,j+1)$. By applying the skewed difference formula (9) to any group of four points located at the said spatial coordinates, the following $(4\times 4)$ system of equations is obtained:

$$\left[\begin{array}{cccc}{\mathcal{I}}_1& -{\mathcal{Q}}_2& 0& 0\\ -{\mathcal{Q}}_2& {\mathcal{I}}_2& 0& 0\\ 0& 0& {\mathcal{I}}_3& -{\mathcal{Q}}_1\\ 0& 0& -{\mathcal{Q}}_1& {\mathcal{I}}_4\end{array}\right]\left[\begin{array}{c}{U}_{i,j}^{k+1}\\ U_{i+1,j+1}^{k+1}\\ U_{i+1,j}^{k+1}\\ U_{i,j+1}^{k+1}\end{array}\right]=\left[\begin{array}{c}rh{s}_{i,j}\\ rh{s}_{i+1,j+1}\\ rh{s}_{i+1,j}\\ rh{s}_{i,j+1}\end{array}\right],$$

(10)

where

$$\begin{array}{cc}\hfill {\sigma }_1& =\frac{1}{\Gamma (2-\gamma (x_i,y_j,t_{k+1/2})){\tau }^{\gamma (x_i,y_j,t_{k+1/2})}},{\mathcal{W}}_1=\frac{{\sigma }_1}{2^{1-\gamma (x_i,y_j,t_{k+1/2})}},\hfill \\ \hfill {\sigma }_2& =\frac{1}{\Gamma (2-\gamma (x_{i+1},y_{j+1},t_{k+1/2})){\tau }^{\gamma (x_{i+1},y_{j+1},t_{k+1/2})}},{\mathcal{W}}_2=\frac{{\sigma }_2}{2^{1-\gamma (x_{i+1},y_{j+1},t_{k+1/2})}},\hfill \\ \hfill {\sigma }_3& =\frac{1}{\Gamma (2-\gamma (x_{i+1},y_j,t_{k+1/2})){\tau }^{\gamma (x_{i+1},y_j,t_{k+1/2})}},{\mathcal{W}}_3=\frac{{\sigma }_3}{2^{1-\gamma (x_{i+1},y_j,t_{k+1/2})}},\hfill \\ \hfill {\sigma }_4& =\frac{1}{\Gamma (2-\gamma (x_i,y_{j+1},t_{k+1/2})){\tau }^{\gamma (x_i,y_{j+1},t_{k+1/2})}},{\mathcal{W}}_4=\frac{{\sigma }_4}{2^{1-\gamma (x_i,y_{j+1},t_{k+1/2})}},\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill {\mathcal{I}}_1& ={\mathcal{W}}_1+0.5\mu +2{\mathcal{Q}}_1+2{\mathcal{Q}}_2,{\mathcal{I}}_2={\mathcal{W}}_2+0.5\mu +2{\mathcal{Q}}_1+2{\mathcal{Q}}_2,\hfill \\ \hfill {\mathcal{I}}_3& ={\mathcal{W}}_3+0.5\mu +2{\mathcal{Q}}_1+2{\mathcal{Q}}_2,{\mathcal{I}}_4={\mathcal{W}}_4+0.5\mu +2{\mathcal{Q}}_1+2{\mathcal{Q}}_2,\hfill \end{array}$$

$$\begin{array}{cc}\hfill rh{s}_{i,j}& ={\mathcal{Q}}_1(U_{i+1,j-1}^{k+1}+U_{i-1,j+1}^{k+1}+U_{i+1,j-1}^k+U_{i-1,j+1}^k)+{\mathcal{Q}}_2(U_{i-1,j-1}^{k+1}\hfill \\ \hfill & +U_{i+1,j+1}^k+U_{i-1,j-1}^k)+({\mathcal{W}}_1-0.5\mu -{\sigma }_1{\mathcal{H}}_1^{i,j,k}-2{\mathcal{Q}}_1-2{\mathcal{Q}}_2)U_{i,j}^k\hfill \\ \hfill & +{\sigma }_1\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right)U_{i,j}^s+{\sigma }_1{\mathcal{H}}_k^{i,j,k}{U}_{i,j}^0+f_{i,j}^{k+1/2},\hfill \\ \hfill rh{s}_{i+1,j+1}& ={\mathcal{Q}}_1(U_{i+2,j}^{k+1}+U_{i,j+2}^{k+1}+U_{i+2,j}^k+U_{i,j+2}^k)+{\mathcal{Q}}_2(U_{i+2,j+2}^{k+1}\hfill \\ \hfill & +U_{i+2,j+2}^k+U_{i,j}^k)+({\mathcal{W}}_2-0.5\mu -{\sigma }_2{\mathcal{H}}_1^{i+1,j+1,k}-2{\mathcal{Q}}_1-2{\mathcal{Q}}_2)U_{i+1,j+1}^k\hfill \\ \hfill & +{\sigma }_2\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i+1,j+1,k}-{\mathcal{H}}_{k-s+1}^{i+1,j+1,k}\right)U_{i+1,j+1}^s+{\sigma }_2{\mathcal{H}}_k^{i+1,j+1,k}{U}_{i+1,j+1}^0+f_{i+1,j+1}^{k+1/2},\hfill \\ \hfill rh{s}_{i+1,j}& ={\mathcal{Q}}_1(U_{i+2,j-1}^{k+1}+U_{i+2,j}^k+U_{i,j+1}^k)+{\mathcal{Q}}_2(U_{i+2,j+1}^{k+1}\hfill \\ \hfill & +U_{i,j-1}^{k+1}+U_{i+2,j+1}^k+U_{i,j-1}^k)+({\mathcal{W}}_3-0.5\mu -{\sigma }_3{\mathcal{H}}_1^{i+1,j,k}-2{\mathcal{Q}}_1-2{\mathcal{Q}}_2)U_{i+1,j}^k\hfill \\ \hfill & +{\sigma }_3\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i+1,j,k}-{\mathcal{H}}_{k-s+1}^{i+1,j,k}\right)U_{i+1,j}^s+{\sigma }_3{\mathcal{H}}_k^{i+1,j,k}{U}_{i+1,j}^0+f_{i+1,j}^{k+1/2},\hfill \\ \hfill rh{s}_{i,j+1}& ={\mathcal{Q}}_1(U_{i-1,j+2}^{k+1}+U_{i+1,j}^k+U_{i-1,j+2}^k)+{\mathcal{Q}}_2(U_{i+1,j+2}^{k+1}\hfill \\ \hfill & +U_{i-1,j}^{k+1}+U_{i+1,j+2}^k+U_{i-1,j}^k)+({\mathcal{W}}_4-0.5\mu -{\sigma }_4{\mathcal{H}}_1^{i,j+1,k}-2{\mathcal{Q}}_1-2{\mathcal{Q}}_2)U_{i,j+1}^k\hfill \\ \hfill & +{\sigma }_4\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j+1,k}-{\mathcal{H}}_{k-s+1}^{i,j+1,k}\right)U_{i,j+1}^s+{\sigma }_4{\mathcal{H}}_k^{i,j+1,k}{U}_{i,j+1}^0+f_{i,j+1}^{k+1/2}.\hfill \end{array}$$

The system defined in (10) can be decoupled as follows:

$$\left[\begin{array}{c}{U}_{i,j}^{k+1}\\ U_{i+1,j+1}^{k+1}\end{array}\right]={\mathcal{K}}_1\left[\begin{array}{cc}-{\mathcal{I}}_2& -{\mathcal{Q}}_2\\ -{\mathcal{Q}}_2& -{\mathcal{I}}_1\end{array}\right]\left[\begin{array}{c}rh{s}_{i,j}\\ rh{s}_{i+1,j+1}\end{array}\right],$$

(11)

and

$$\left[\begin{array}{c}{U}_{i+1,j}^{k+1}\\ U_{i,j+1}^{k+1}\end{array}\right]={\mathcal{K}}_2\left[\begin{array}{cc}-{\mathcal{I}}_4& -{\mathcal{Q}}_1\\ -{\mathcal{Q}}_1& -{\mathcal{I}}_3\end{array}\right]\left[\begin{array}{c}rh{s}_{i+1,j}\\ rh{s}_{i,j+1}\end{array}\right],$$

(12)

where

$$\begin{array}{cc}\hfill {\mathcal{K}}_1^{-1}=& 3{\mathcal{Q}}_2^2+{\mathcal{W}}_1{\mathcal{W}}_2+0.5\mu {\mathcal{W}}_1+2{\mathcal{W}}_1{\mathcal{Q}}_1+2{\mathcal{W}}_1{\mathcal{Q}}_2+0.5\mu {\mathcal{W}}_2+0.25{\mu }^2\hfill \\ \hfill & +2\mu {\mathcal{Q}}_1+2\mu {\mathcal{Q}}_2+2{\mathcal{W}}_2{\mathcal{Q}}_1+4{\mathcal{Q}}_1^2+8{\mathcal{Q}}_1{\mathcal{Q}}_2+2{\mathcal{W}}_2{\mathcal{Q}}_2,\hfill \\ \hfill {\mathcal{K}}_2^{-1}=& 3{\mathcal{Q}}_1^2+{\mathcal{W}}_3{\mathcal{W}}_4+0.5\mu {\mathcal{W}}_3+2{\mathcal{W}}_3{\mathcal{Q}}_1+2{\mathcal{W}}_3{\mathcal{Q}}_2+0.5\mu {\mathcal{W}}_4+0.25{\mu }^2\hfill \\ \hfill & +2\mu {\mathcal{Q}}_1+2\mu {\mathcal{Q}}_2+2{\mathcal{W}}_4{\mathcal{Q}}_1+4{\mathcal{Q}}_2^2+8{\mathcal{Q}}_1{\mathcal{Q}}_2+2{\mathcal{W}}_4{\mathcal{Q}}_2.\hfill \end{array}$$

With reference to Figure 2, the computational molecule of the EDGM is obtained by arranging the nodal points of the solution domain into groups of four points. One can easily verify that the execution of Equation (11) corresponds to the solution values located at the ⧫ points, whereas the implementation of Equation (12) corresponds to the solution outcomes placed at the ◊ points. Therefore, iterative computations can be performed independently on either one type of point until convergence is attained. Regardless of the selected point type that takes part in the iterative process, the solution values located at the remaining grid points will be evaluated directly once using the CN–FDM given in Section 2.

In this work, the EDGM is combined with the Gauss–Seidel iterative solver to generate iterations at the ⧫ points. Once convergence is achieved, the solution outcomes at the residual ◊ points will be calculated directly once using Equation (7). The prescribed process is terminated when the targeted time level is reached. Compared to the point-wise iterative method presented in Section 2, the EDGM is expected to result in faster simulations since only half of the nodal points are involved in the iteration process, which reduces the computational complexity effectively.

4. Stability Analysis

In this section, the stability of the proposed discrete numerical schemes is investigated via the technique of von Neumann analysis. As a start, the following lemma is introduced.

Lemma 1

([63]). Let ${\mathcal{H}}_s^{i,j,k}$, $0\le i\le M_x-1$, $0\le j\le M_y-1$ be the coefficients defined in Equation (7), then the following hold:

(i)

${\mathcal{H}}_s^{i,j,k}\ge {\mathcal{H}}_{s+1}^{i,j,k},1\le s\le k=1,2,\cdots ,N-1.$

(ii)

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

Suppose that ${\widehat{U}}_{i,j}^k$ and ${\overline{U}}_{i,j}^k$ are the approximate solutions of Equations (7) and (9), respectively. Then, the errors are given by

$$\begin{array}{c}\hfill {\psi }_{i,j}^k=U_{i,j}^k-{\widehat{U}}_{i,j}^k,0\le i\le M_x,0\le j\le M_y,1\le k\le N,\end{array}$$

(13)

$${\varphi }_{i,j}^k=U_{i,j}^k-{\overline{U}}_{i,j}^k,0\le i\le M_x,0\le j\le M_y,1\le k\le N.$$

(14)

Using Equation (13), the round-off error equation for the discrete scheme (7) is as follows:

$$\begin{array}{cc}\hfill & ({\mathcal{W}}_1+0.5\mu +2{\mathcal{P}}_1+2{\mathcal{P}}_2){\psi }_{i,j}^{k+1}-{\mathcal{P}}_1({\psi }_{i+1,j}^{k+1}+{\psi }_{i-1,j}^{k+1})-{\mathcal{P}}_2({\psi }_{i,j+1}^{k+1}+{\psi }_{i,j-1}^{k+1})=\hfill \\ \hfill & {\mathcal{P}}_1({\psi }_{i+1,j}^k+{\psi }_{i-1,j}^k)+{\mathcal{P}}_2({\psi }_{i,j+1}^k+{\psi }_{i,j-1}^k)+({\mathcal{W}}_1-0.5\mu -{\sigma }_1{\mathcal{H}}_1^{i,j,k}-2{\mathcal{P}}_1\hfill \\ \hfill & -2{\mathcal{P}}_2){\psi }_{i,j}^k+{\sigma }_1\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right){\psi }_{i,j}^s+{\sigma }_1{\mathcal{H}}_k^{i,j,k}{\psi }_{i,j}^0.\hfill \end{array}$$

(15)

Similarly, utilizing Equation (14), the round-off error equation for the numerical scheme (9) is in the form

$$\begin{array}{cc}\hfill & ({\mathcal{W}}_1+0.5\mu +2{\mathcal{Q}}_1+2{\mathcal{Q}}_2){\varphi }_{i,j}^{k+1}-{\mathcal{Q}}_1({\varphi }_{i+1,j-1}^{k+1}+{\varphi }_{i-1,j+1}^{k+1})-{\mathcal{Q}}_2({\varphi }_{i+1,j+1}^{k+1}\hfill \\ \hfill & +{\varphi }_{i-1,j-1}^{k+1}){\mathcal{Q}}_1({\varphi }_{i+1,j-1}^k+{\varphi }_{i-1,j+1}^k)+{\mathcal{Q}}_2({\varphi }_{i+1,j+1}^k+{\varphi }_{i-1,j-1}^k)+({\mathcal{W}}_1\hfill \\ \hfill & -0.5\mu -{\sigma }_1{\mathcal{H}}_1^{i,j,k}-2{\mathcal{Q}}_1-2{\mathcal{Q}}_2){\psi }_{i,j}^k+{\sigma }_1\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right){\varphi }_{i,j}^s\hfill \\ \hfill & +{\sigma }_1{\mathcal{H}}_k^{i,j,k}{\varphi }_{i,j}^0.\hfill \end{array}$$

(16)

Then, the Fourier series expansions of ${\psi }^k(x,y)$ and ${\varphi }^k(x,y)$ are defined as

$$\begin{array}{cc}\hfill {\psi }^k(x,y)& =\sum _{Z_1=-\infty }^{\infty }\sum _{Z_2=-\infty }^{\infty }{\lambda }^k(Z_1,Z_2)e^{2\pi I(Z_{1}x/L+Z_{2}y/L)},\hfill \\ \hfill {\varphi }^k(x,y)& =\sum _{Z_1=-\infty }^{\infty }\sum _{Z_2=-\infty }^{\infty }{\rho }^k(Z_1,Z_2)e^{2\pi I(Z_{1}x/L+Z_{2}y/L)},\hfill \end{array}$$

(17)

where $I=\sqrt{-1}$ and ${\lambda }^k$ and ${\rho }^k$ have the following form:

$${\lambda }^k(Z_1,Z_2)=\frac{1}{L^2}{\int }_0^L{\int }_0^L{\psi }^k(x,y)e^{-2\pi I(Z_{1}x/L+Z_{2}y/L)}dxdy,$$

(18)

$${\rho }^k(Z_1,Z_2)=\frac{1}{L^2}{\int }_0^L{\int }_0^L{\varphi }^k(x,y)e^{-2\pi I(Z_{1}x/L+Z_{2}y/L)}dxdy.$$

(19)

With the help of the $l^2$ norm and Parseval’s equality, we obtain

$$\begin{array}{cc}\hfill \parallel {\psi }^k{\parallel }_2& ={\left(\sum _{Z_2=-\infty }^{\infty }\sum _{Z_1=-\infty }^{\infty }{\left|{\lambda }^k(Z_1,Z_2)\right|}^2\right)}^{1/2},\hfill \\ \hfill \parallel {\varphi }^k{\parallel }_2& ={\left(\sum _{Z_2=-\infty }^{\infty }\sum _{Z_1=-\infty }^{\infty }{\left|{\rho }^k(Z_1,Z_2)\right|}^2\right)}^{1/2}.\hfill \end{array}$$

(20)

Assume that the solutions of Equations (15) and (16) are expressed as

$${\psi }_{i,j}^k={\lambda }^k{e}^{I({\nu }_{1}i{h}_x+{\nu }_{2}j{h}_y)},{\varphi }_{i,j}^k={\rho }^k{e}^{I({\nu }_{1}i{h}_x+{\nu }_{2}j{h}_y)},$$

(21)

where ${\nu }_1=2\pi Z_1/L$ and ${\nu }_2=2\pi Z_2/L$.

Setting ${\psi }_{i,j}^k={\lambda }^k{e}^{I({\nu }_{1}i{h}_x+{\nu }_{2}j{h}_y)}$ into Equation (15) and simplifying the result, we obtain

$$\begin{array}{cc}\hfill {\lambda }^{k+1}=& \frac{{\mathcal{W}}_1-0.5\mu -{\chi }_1-{\chi }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}{\lambda }^k\hfill \\ \hfill & +\frac{1}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}\left[{\sigma }_1\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right){\lambda }^s+{\sigma }_1{\mathcal{H}}_k^{i,j,k}{\lambda }^0\right],\hfill \end{array}$$

(22)

where

$${\chi }_1=4{\mathcal{P}}_1{sin}^2\left(\frac{{\nu }_1{h}_x}{2}\right),{\chi }_2=4{\mathcal{P}}_2{sin}^2\left(\frac{{\nu }_2{h}_y}{2}\right).$$

Lemma 2.

Suppose that ${\lambda }^k$$(0\le k\le N-1)$ are the solutions of Equation (22) and $2\ge 3^{1-\gamma (x_i,y_j,t_{k+1/2})}$, then we have

$$|{\lambda }^{k+1}|\le |{\lambda }^0|,k=0,1,\cdots ,N-1.$$

(23)

Proof. 

To complete the proof, mathematical induction is utilized. First, choose $k=0$, and we obtain

$$|{\lambda }^1|=|\frac{{\mathcal{W}}_1-0.5\mu -{\chi }_1-{\chi }_2}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}\left|\right|{\lambda }^0|.$$

Since ${\mathcal{W}}_1$, ${\chi }_1$ and ${\chi }_2$ are non-negative constants, then

$$|{\lambda }^1|\le |{\lambda }^0|.$$

Now, assume that $|{\lambda }^{m+1}|\le |{\lambda }^0|$, $m=0,1,2,\cdots ,k-1$. We prove it is true for $m=k$. With the help of Equation (22) and Lemma 1, we obtain

$$\begin{array}{cc}\hfill |{\lambda }^{k+1}|\le & |\frac{{\mathcal{W}}_1-0.5\mu -{\chi }_1-{\chi }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}\left|\right|{\lambda }^k|\hfill \\ \hfill & +|\frac{1}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}|\left[{\sigma }_1\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right)|{\lambda }^s|+{\sigma }_1{\mathcal{H}}_k^{i,j,k}\left|{\lambda }^0\right|\right],\hfill \\ \hfill & \le |\frac{{\mathcal{W}}_1-0.5\mu -{\chi }_1-{\chi }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}\left|\right|{\lambda }^0|\hfill \\ \hfill & +\frac{1}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}\left[{\sigma }_1\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right)|{\lambda }^0|+{\sigma }_1{\mathcal{H}}_k^{i,j,k}\left|{\lambda }^0\right|\right],\hfill \\ \hfill & =\frac{|{\mathcal{W}}_1-0.5\mu -{\chi }_1-{\chi }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}|+{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu {\chi }_1+{\chi }_2}\left|{\lambda }^0\right|.\hfill \end{array}$$

If ${\mathcal{W}}_1-0.5\mu -{\chi }_1-{\chi }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}>0$, then

$$|{\lambda }^{k+1}|\le \frac{{\mathcal{W}}_1-0.5\mu -{\chi }_1-{\chi }_2}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}|{\lambda }^0|<|{\lambda }^0|.$$

If ${\mathcal{W}}_1-0.5\mu -{\chi }_1-{\chi }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}<0$, then

$$|{\lambda }^{k+1}|\le \frac{-{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2+2{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}\left|{\lambda }^0\right|.$$

Here,

$$\begin{array}{cc}\hfill |{\lambda }^{k+1}|& \le |{\lambda }^0|\hfill \\ \hfill & \iff \frac{-{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2+2{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}\le 1\hfill \\ \hfill & \iff -{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2+2{\sigma }_1{\mathcal{H}}_1^{i,j,k}\le {\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2\hfill \\ \hfill & \iff 2\ge 3^{1-\gamma (x_i,y_j,t_{k+1/2})}.\hfill \end{array}$$

□

Theorem 1.

Given that $2\ge 3^{1-\gamma (x_i,y_j,t_{k+1/2})}$, the fully discrete numerical scheme (7) is stable.

Proof. 

In view of lemma 2 and the Parseval equality, it follows that,

$$\parallel {\psi }^k{\parallel }_2^2=\sum _{Z_1=-\infty }^{\infty }\sum _{Z_2=-\infty }^{\infty }|{\lambda }^k(Z_1,Z_2){|}^2\le \sum _{Z_1=-\infty }^{\infty }\sum _{Z_2=-\infty }^{\infty }|{\lambda }^0(Z_1,Z_2){|}^2={\parallel {\psi }^0\parallel }_2^2,$$

from which we obtain

$$\parallel {\psi }^k\parallel \le \parallel {\psi }^0\parallel ,0\le k\le N.$$

□

Substituting ${\varphi }_{i,j}^k={\rho }^k{e}^{I({\nu }_{1}i{h}_x+{\nu }_{2}j{h}_y)}$ into Equation (16) yields

$$\begin{array}{cc}\hfill {\rho }^{k+1}=& \frac{{\mathcal{W}}_1-0.5\mu -{\eta }_1-{\eta }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\eta }_1+{\eta }_2}{\rho }^k\hfill \\ \hfill & +\frac{1}{{\mathcal{W}}_1+0.5\mu +{\eta }_1+{\eta }_2}\left[{\sigma }_1\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right){\rho }^s+{\sigma }_1{\mathcal{H}}_k^{i,j,k}{\rho }^0\right],\hfill \end{array}$$

(24)

where

$${\eta }_1=4{\mathcal{Q}}_1{sin}^2\left(\frac{{\nu }_1{h}_x-{\nu }_2{h}_y}{2}\right),{\eta }_2=4{\mathcal{Q}}_2{sin}^2\left(\frac{{\nu }_1{h}_x+{\nu }_2{h}_y}{2}\right).$$

Lemma 3.

Suppose that ${\rho }^k$$(0\le k\le N-1)$ are the solutions of Equation (24) and $2\ge 3^{1-\gamma (x_i,y_j,t_{k+1/2})}$, then we have

$$|{\rho }^{k+1}|\le |{\rho }^0|,k=0,1,\cdots ,N-1.$$

(25)

Proof. 

First, putting $k=0$ in Equation (24), we obtain

$$|{\rho }^1|=|\frac{{\mathcal{W}}_1-0.5\mu -{\eta }_1-{\eta }_2}{{\mathcal{W}}_1+0.5\mu +{\eta }_1+{\eta }_2}\left|\right|{\rho }^0|<|{\rho }^0|$$

Next, assume that $|{\rho }^{m+1}|\le |{\rho }^0|$, $m=0,1,2,\cdots ,k-1$. We prove it is true for $m=k$. With the help of Equation (24) and Lemma 1, we obtain

$$\begin{array}{cc}\hfill |{\rho }^{k+1}|\le & |\frac{{\mathcal{W}}_1-0.5\mu -{\eta }_1-{\eta }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\eta }_1+{\eta }_2}\left|\right|{\rho }^k|\hfill \\ \hfill & +|\frac{1}{{\mathcal{W}}_1+0.5\mu +{\eta }_1+{\eta }_2}|\left[{\sigma }_1\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right)|{\rho }^s|+{\sigma }_1{\mathcal{H}}_k^{i,j,k}\left|{\rho }^0\right|\right]\hfill \\ \hfill & \le |\frac{{\mathcal{W}}_1-0.5\mu -{\eta }_1-{\eta }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\eta }_1+{\eta }_2}\left|\right|{\rho }^0|\hfill \\ \hfill & +\frac{1}{{\mathcal{W}}_1+0.5\mu +{\eta }_1+{\eta }_2}\left[{\sigma }_1\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right)|{\rho }^0|+{\sigma }_1{\mathcal{H}}_k^{i,j,k}\left|{\rho }^0\right|\right]\hfill \\ \hfill & =\frac{|{\mathcal{W}}_1-0.5\mu -{\eta }_1-{\eta }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}|+{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\eta }_1+{\eta }_2}\left|{\rho }^0\right|.\hfill \end{array}$$

If ${\mathcal{W}}_1-0.5\mu -{\eta }_1-{\eta }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}>0$, then

$$|{\rho }^{k+1}|\le \frac{{\mathcal{W}}_1-0.5\mu -{\eta }_1-{\eta }_2}{{\mathcal{W}}_1+0.5\mu +{\eta }_1+{\eta }_2}|{\rho }^0|<|{\rho }^0|.$$

If ${\mathcal{W}}_1-0.5\mu -{\eta }_1-{\eta }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}<0$, then

$$|{\rho }^{k+1}|\le \frac{-{\mathcal{W}}_1+0.5\mu +{\eta }_1+{\eta }_2+2{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\eta }_1+{\eta }_2}\left|{\rho }^0\right|.$$

Here,

$$\begin{array}{cc}\hfill |{\rho }^{k+1}|& \le |{\rho }^0|\hfill \\ \hfill & \iff \frac{-{\mathcal{W}}_1+0.5\mu +{\eta }_1+{\eta }_2+2{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\eta }_1+{\eta }_2}\le 1\hfill \\ \hfill & \iff -{\mathcal{W}}_1+0.5\mu +{\eta }_1+{\eta }_2+2{\sigma }_1{\mathcal{H}}_1^{i,j,k}\le {\mathcal{W}}_1+0.5\mu +{\eta }_1+{\eta }_2\hfill \\ \hfill & \iff 2\ge 3^{1-\gamma (x_i,y_j,t_{k+1/2})}.\hfill \end{array}$$

□

Theorem 2.

Given that $2\ge 3^{1-\gamma (x_i,y_j,t_{k+1/2})}$, the fully discrete numerical scheme (9) is stable.

Proof. 

Utilizing Lemma 3 and Parseval equality, we obtain

$$\parallel {\varphi }^k{\parallel }_2^2=\sum _{Z_1=-\infty }^{\infty }\sum _{Z_2=-\infty }^{\infty }|{\rho }^k(Z_1,Z_2){|}^2\le \sum _{Z_1=-\infty }^{\infty }\sum _{Z_2=-\infty }^{\infty }|{\rho }^0(Z_1,Z_2){|}^2={\parallel {\varphi }^0\parallel }_2^2.$$

From the above result, it immediately follows,

$$\parallel {\varphi }^k\parallel \le \parallel {\varphi }^0\parallel ,0\le k\le N.$$

□

5. Convergence Analysis

In this section, we intend to analyze the convergence of the proposed numerical schemes.

Subtracting Equation (7) from Equation (6), the error equation can be obtained as

$$\begin{array}{cc}\hfill & ({\mathcal{W}}_1+0.5\mu +2{\mathcal{P}}_1+2{\mathcal{P}}_2){\Psi }_{i,j}^{k+1}-{\mathcal{P}}_1({\Psi }_{i+1,j}^{k+1}+{\Psi }_{i-1,j}^{k+1})-{\mathcal{P}}_2({\Psi }_{i,j+1}^{k+1}+{\Psi }_{i,j-1}^{k+1})=\hfill \\ \hfill & {\mathcal{P}}_1({\Psi }_{i+1,j}^k+{\Psi }_{i-1,j}^k)+{\mathcal{P}}_2({\Psi }_{i,j+1}^k+{\Psi }_{i,j-1}^k)+({\mathcal{W}}_1-{\sigma }_1{\mathcal{H}}_1^{i,j,k}-2{\mathcal{P}}_1-2{\mathcal{P}}_2){\Psi }_{i,j}^k\hfill \\ \hfill & +{\sigma }_1\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right){\Psi }_{i,j}^s+{\sigma }_1{\mathcal{H}}_k^{i,j,k}{\Psi }_{i,j}^0+R_{i,j}^{k+1/2},\hfill \end{array}$$

(26)

in which

$${\Psi }_{i,j}^k=u_{i,j}^k-U_{i,j}^k,1\le i\le M_x-1,1\le j\le M_y-1,1\le k\le N,$$

(27)

and $R_{i,j}^{k+1/2}$ is the truncation error at the nodal point $(x_i,y_j,t_{k+1/2})$. Henceforward, C is a constant that may take various values at different positions. Based on Equation (6), we have

$$|R_{i,j}^{k+1/2}|\le C(\tau +h_x^2+h_y^2),1\le i\le M_x-1,1\le j\le M_y-1,0\le k\le N,$$

(28)

where

$$C=\underset{1\le i\le M_x-1,1\le j\le M_y-1,0\le k\le N}{max}\{C_{i,j}^k\}.$$

The Fourier series expansions of ${\Psi }^k(x,y)$ and $R^{k+1/2}(x,y)$ can be written as

$${\Psi }^k(x,y)=\sum _{Z_2=-\infty }^{\infty }\sum _{Z_1=-\infty }^{\infty }{\xi }^k(Z_1,Z_2)e^{2\pi I(Z_{1}x/L+Z_{2}y/L)},$$

$$R^{k+1/2}(x,y)=\sum _{Z_2=-\infty }^{\infty }\sum _{Z_1=-\infty }^{\infty }{\phi }^k(Z_1,Z_2)e^{2\pi I(Z_{1}x/L+Z_{2}y/L)},$$

in which

$${\xi }^k(Z_1,Z_2)=\frac{1}{L^2}{\int }_0^L{\int }_0^L{\Psi }^k(x,y)e^{-2\pi I(Z_{1}x/L+Z_{2}y/L)}dxdy,$$

$${\phi }^k(Z_1,Z_2)=\frac{1}{L^2}{\int }_0^L{\int }_0^L{R}^{k+1/2}(x,y)e^{-2\pi I(Z_{1}x/L+Z_{2}y/L)}dxdy.$$

Based on the $l^2$ norm and Parseval’s equality, we have

$$\parallel {\Psi }^k{\parallel }_2={\left(\sum _{j=1}^{M_y-1}\sum _{i=1}^{M_x-1}{h}_y{h}_x{\left|{\Psi }_{i,j}^k\right|}^2\right)}^{1/2}={\left(\sum _{Z_1=-\infty }^{\infty }\sum _{Z_2=-\infty }^{\infty }{\left|{\xi }^k(Z_1,Z_2)\right|}^2\right)}^{1/2},$$

(29)

$$\parallel R^{k+1/2}{\parallel }_2={\left(\sum _{j=1}^{M_y-1}\sum _{i=1}^{M_x-1}{h}_y{h}_x{\left|R_{i,j}^{k+1/2}\right|}^2\right)}^{1/2}={\left(\sum _{Z_1=-\infty }^{\infty }\sum _{Z_2=-\infty }^{\infty }{\left|{\phi }^k(Z_1,Z_2)\right|}^2\right)}^{1/2}.$$

(30)

Next, assume that the solutions of Equation (26) are as follows:

$${\Psi }_{i,j}^k={\xi }^k{e}^{I({\nu }_{1}i{h}_x+{\nu }_{2}j{h}_y)},R_{i,j}^{k+1/2}={\phi }^{k+1/2}{e}^{I({\nu }_{1}i{h}_x+{\nu }_{2}j{h}_y)}.$$

(31)

Setting (31) into (26) and simplifying the result, we obtain

$$\begin{array}{cc}\hfill {\xi }^{k+1}=& \frac{{\mathcal{W}}_1-0.5\mu -{\chi }_1-{\chi }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}{\xi }^k\hfill \\ \hfill & +\frac{1}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}\left[{\sigma }_1\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right){\xi }^s+{\sigma }_1{\mathcal{H}}_k^{i,j,k}{\xi }^0+{\phi }^{k+1/2}\right],\hfill \end{array}$$

(32)

where ${\chi }_1$ and ${\chi }_2$ are as given in Equation (22).

Lemma 4.

Suppose that ${\xi }^k$$(1\le k\le N)$ are the solutions of Equation (32) and $2\ge 3^{1-\gamma (x_i,y_j,t_{k+1/2})}$, then we have

$$|{\xi }^{k+1}|\le C(k+1)\tau |{\phi }^{1/2}|,k=0,1,\cdots ,N-1.$$

(33)

Proof. 

From Equation (30), there exists a positive constant C such that

$$|{\phi }^{k+1/2}|\le C\tau |{\phi }^{1/2}|,0\le k\le N-1.$$

The proof is completed by mathematical induction. First, letting $k=0$ in Equation (31) and observing that ${\xi }^0=0$, we obtain

$$|{\xi }^1|=\frac{1}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}|{\phi }^{1/2}|\le |{\phi }^{1/2}|\le C\tau |{\phi }^{1/2}|.$$

Next, suppose that

$$|{\xi }^{m+1}|\le C(m+1)\tau |{\phi }^{1/2}|,m=0,1,2,\cdots ,k-1.$$

Setting $m=k$ in Equation (31) and utilizing Lemma 1, we obtain

$$\begin{array}{cc}\hfill |{\xi }^{k+1}|\le & |\frac{{\mathcal{W}}_1-0.5\mu -{\chi }_1-{\chi }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}\left|\right|{\xi }^k|\hfill \\ \hfill & +|\frac{1}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}|\left[{\sigma }_1\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right)|{\xi }^s|+{\sigma }_1{\mathcal{H}}_k^{i,j,k}|{\xi }^0|+|{\phi }^{k+1/2}|\right]\hfill \\ \hfill & \le \frac{|{\mathcal{W}}_1-0.5\mu -{\chi }_1-{\chi }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}|}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}Ck\tau \left|{\phi }^{1/2}\right|\hfill \\ \hfill & +\frac{1}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}\left[{\sigma }_{1}Ck\tau |{\phi }^{1/2}|\left(\sum _{s=1}^{k-1}\left({\mathcal{H}}_{k-s}^{i,j,k}-{\mathcal{H}}_{k-s+1}^{i,j,k}\right)+{\mathcal{H}}_k^{i,j,k}\right)+C\tau \left|{\phi }^{1/2}\right|\right]\hfill \\ \hfill & =\frac{|{\mathcal{W}}_1-0.5\mu -{\chi }_1-{\chi }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}|+{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}Ck\tau \left|{\phi }^{1/2}\right|\hfill \\ \hfill & +\frac{1}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}C\tau \left|{\phi }^{1/2}\right|\hfill \\ \hfill & =\left[\frac{|{\mathcal{W}}_1-0.5\mu -{\chi }_1-{\chi }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}|+{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}k+\frac{1}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}\right]C\tau \left|{\phi }^{1/2}\right|.\hfill \end{array}$$

If ${\mathcal{W}}_1-0.5\mu -{\chi }_1-{\chi }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}>0$, then

$$\begin{array}{cc}\hfill |{\xi }^{k+1}|& \le \left[\frac{{\mathcal{W}}_1-0.5\mu -{\chi }_1-{\chi }_2}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}k+\frac{1}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}\right]C\tau \left|{\phi }^{1/2}\right|\hfill \\ \hfill & \le C(k+1)\tau |{\phi }^{1/2}|.\hfill \end{array}$$

If ${\mathcal{W}}_1-0.5\mu -{\chi }_1-{\chi }_2-{\sigma }_1{\mathcal{H}}_1^{i,j,k}<0$, then

$$\begin{array}{cc}\hfill |{\xi }^{k+1}|& \le \left[\frac{-{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2+2{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}k+\frac{1}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}\right]C\tau \left|{\phi }^{1/2}\right|.\hfill \end{array}$$

Here,

$$\begin{array}{cc}\hfill |{\xi }^{k+1}|& \le C(k+1)\tau |{\phi }^{1/2}|\hfill \\ \hfill & \iff \frac{-{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2+2{\sigma }_1{\mathcal{H}}_1^{i,j,k}}{{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2}\le 1\hfill \\ \hfill & \iff -{\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2+2{\sigma }_1{\mathcal{H}}_1^{i,j,k}\le {\mathcal{W}}_1+0.5\mu +{\chi }_1+{\chi }_2\hfill \\ \hfill & \iff 2\ge 3^{1-\gamma (x_i,y_j,t_{k+1/2})}.\hfill \end{array}$$

Therefore, we have

$$|{\xi }^{k+1}|\le C(k+1)\tau |{\phi }^{1/2}|.$$

□

Theorem 3.

Given that $2\ge 3^{1-\gamma (x_i,y_j,t_{k+1/2})}$, the fully discrete numerical scheme (7) is $l^2$ convergent, and the convergence order is $O(\tau +h_x^2+h_y^2)$.

Proof. 

By employing Lemma 4 and Equations (29) and (30), we obtain

$$\begin{array}{cc}\hfill \parallel {\Psi }^k{\parallel }_2^2=\sum _{Z_1=-\infty }^{\infty }\sum _{Z_2=-\infty }^{\infty }{\left|{\xi }^k(Z_1,Z_2)\right|}^2& \le \sum _{Z_1=-\infty }^{\infty }\sum _{Z_2=-\infty }^{\infty }{C}^2{(k+1)}^2{\tau }^2{\left|{\phi }^{1/2}(Z_1,Z_2)\right|}^2\hfill \\ \hfill & =C^2{(k+1)}^2{\tau }^2{\parallel R^{1/2}\parallel }_2^2.\hfill \end{array}$$

Using Equation (28), we obtain

$$\begin{array}{cc}\hfill \parallel {\Psi }^k{\parallel }_2& \le C(k+1)\tau \parallel R^{1/2}{\parallel }_2\le C(k+1)\tau (\tau +h_x^2+h_y^2)\hfill \\ \hfill & \le C^{\ast }(\tau +h_x^2+h_y^2),\hfill \end{array}$$

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

Theorem 4.

The fully discrete numerical scheme (9) is $l^2$ convergent and the convergence order is $O(\tau +h_x^2+h_y^2)$.

Proof. 

The proof can be established in a similar way to Theorem 3. □

6. Numerical Results and Discussion

In this part, we investigate the viability, correctness and efficiency of the proposed group-wise EDGM in comparison to the point-wise CN–FDM presented in Section 2. To this end, two numerical simulations accompanied by their computational outputs in terms of CPU time $\left(CPU\right)$, number of iterations $\left(Iter\right)$ and maximum absolute error $\left(MAE\right)$ are provided in tabular and graphical forms. The Gauss–Seidel iterative scheme with stopping criteria ${10}^{-5}$ and the $l^{\infty }$ norm is utilized to solve the linear systems. All numerical experiments are conducted using Julia programming language and run on a PC with the configuration: Intel(R) Core(TM) i7-8550U and 8 GB of RAM. In the runs, a uniform grid spacing $h=h_x=h_y$ is used in the horizontal and vertical spatial directions, while the $MAE$ between the exact and numerical solutions is computed as

$$\begin{array}{c}\hfill MAE=\underset{1\le i\le M_x-1,1\le j\le M_y-1,1\le k\le N}{max}u(x_i,y_j,t_k)-U_{i,j}^k.\end{array}$$

Example 1.

In the first example, we consider the two-dimensional variable-order fractional cable equation:

$$\begin{array}{cc}\hfill {}_0{}^{C}D_t^{\gamma (x,y,t)}u(x,y,t)& =\frac{{\mathsf{\partial }}^{2}u(x,y,t)}{\mathsf{\partial }{x}^2}+\frac{{\mathsf{\partial }}^{2}u(x,y,t)}{\mathsf{\partial }{y}^2}-u(x,y,t)\hfill \\ \hfill & +f(x,y,t),(x,y,t)\in \Omega \times [0,1],\hfill \end{array}$$

with the initial-boundary conditions

$$\begin{array}{cc}\hfill u(x,y,0)& =0,(x,y)\in \Omega \hfill \\ \hfill u(0,y,t)& =t^{2}sin\left(\pi y\right),u(1,y,t)=t^{2}sin\left(\pi y\right),(x,y,t)\in [0,1]\times \mathsf{\partial }\Omega ,\hfill \\ \hfill u(x,0,t)& =t^{2}sin\left(\pi x\right),u(x,1,t)=t^{2}sin\left(\pi x\right),\hfill \end{array}$$

where $\Omega =(0,1)\times (0,1)$, and the source term is given by

$$\begin{array}{cc}\hfill f(x,y,t)& =\frac{2t^{2-\gamma (x,y,t)}}{\Gamma (3-\gamma (x,y,t\left)\right)}(sin\left(\pi x\right)+sin\left(\pi y\right))+({\pi }^2+1)t^2(sin\left(\pi x\right)+sin\left(\pi y\right)).\hfill \end{array}$$

The exact analytical solution of the stated problem is of the form

$$u(x,y,t)=t^2(sin\left(\pi x\right)+sin\left(\pi y\right)).$$

The considered problem is solved for different choices of $\gamma (x,y,t)$ and varying spatial and temporal step sizes. The CPU times, iteration numbers, and numerical errors in solving Example 1 are recorded in

Table 1. Comparison of numerical results of the proposed methods for Example 1 at $\tau =0.01$, $T=1$ and various choices of $\gamma (x,y,t)$.

		CN–FDM			EDGM
$\mathit{\gamma }(\mathit{x},\mathit{y},\mathit{t})$	${\mathit{h}}^{-\mathbf{1}}$	$\mathit{CPU}$	$\mathit{Ite}$	$\mathit{MAE}$	$\mathit{CPU}$	$\mathit{Ite}$	$\mathit{MAE}$
$\frac{1-xyt+sin\left(xyt\right)}{10}$	6	3.2214	24	2.4354 × 10−2	1.5446	11	2.5830 × 10−2
	12	57.0982	75	6.0448 × 10−3	14.1272	34	6.1778 × 10−3
	18	250.0227	145	2.5440 × 10−3	51.5105	66	2.6716 × 10−3
	24	623.4341	230	1.2371 × 10−3	148.0147	106	1.4167 × 10−3
$\frac{5+{\left(xy\right)}^3-{\left(yt\right)}^4}{50}$	6	4.3391	24	2.4449 × 10−2	1.6299	12	2.5933 × 10−2
	12	63.6731	78	6.0859 × 10−3	17.5187	35	6.2206 × 10−3
	18	296.8552	150	2.5441 × 10−3	72.8325	68	2.6855 × 10−3
	24	765.8251	237	1.2493 × 10−3	212.8129	109	1.4199 × 10−3
$\frac{3+{\left(xy\right)}^2-{\left(yt\right)}^3}{30}$	6	3.2722	24	2.4449 × 10−2	1.3767	12	2.5934 × 10−2
	12	58.3578	78	6.0862 × 10−3	16.7624	35	6.2144 × 10−3
	18	250.1333	150	2.5414 × 10−3	67.1791	68	2.6880 × 10−3
	24	723.0534	237	1.2615 × 10−3	169.6976	109	1.4235 × 10−3
$\frac{16-e^{xyt}}{17}$	6	1.3651	9	2.3432 × 10−2	0.8389	5	2.4777 × 10−2
	12	14.9717	21	5.7997 × 10−3	4.4889	11	5.9192 × 10−3
	18	64.7063	39	2.4295 × 10−3	17.991	18	2.5349 × 10−3
	24	226.331	61	1.1482 × 10−3	63.1113	28	1.3300 × 10−3
$\frac{e^{xyt}+sin\left(xyt\right)}{30}$	6	3.2977	25	2.4525 × 10−2	1.3007	12	2.6012 × 10−2
	12	44.7645	79	6.1073 × 10−3	12.388	36	6.2598 × 10−3
	18	215.872	152	2.5741 × 10−3	54.5599	69	2.7378 × 10−3
	24	671.5466	240	1.3141 × 10−3	157.9549	111	1.4666 × 10−3
$\frac{1-{\left(xyt\right)}^3+sin{\left(xyt\right)}^2}{10}$	6	3.7105	24	2.4446 × 10−2	2.1392	12	2.5927 × 10−2
	12	48.9569	77	6.0763 × 10−3	14.094	35	6.2183 × 10−3
	18	247.4537	150	2.5423 × 10−3	51.0625	68	2.6972 × 10−3
	24	658.7823	236	1.2699 × 10−3	153.9475	109	1.4156 × 10−3

and

Table 2. Comparison of numerical results of the proposed methods for Example 1 at $h^{-1}=30$, $T=1$ and different choices of $\gamma (x,y,t)$.

		CN–FDM			EDGM
$\mathit{\gamma }(\mathit{x},\mathit{y},\mathit{t})$	$\mathit{\tau }$	$\mathit{CPU}$	$\mathit{Ite}$	$\mathit{MAE}$	$\mathit{CPU}$	$\mathit{Ite}$	$\mathit{MAE}$
$(16-e^{xyt})/17$	4	7.6451	493	3.3805 × 10−2	2.3588	211	3.3457 × 10−2
	8	22.1864	388	7.9937 × 10−3	5.0926	167	7.6691 × 10−3
	16	49.5047	283	2.0338 × 10−3	12.9866	123	1.7656 × 10−3
	32	123.0862	191	8.4239 × 10−4	31.4933	84	5.8782 × 10−4
$5^{yt-3}$	4	13.3397	599	4.1630 × 10−2	3.6625	257	4.1215 × 10−2
	8	33.5994	546	1.1022 × 10−2	7.8271	236	1.0530 × 10−2
	16	101.0359	491	3.3731 × 10−3	24.7500	216	2.8573 × 10−3
	32	307.0350	435	1.5878 × 10−3	77.1075	194	9.8333 × 10−4
$(11-{cos}^2\left(xt\right))/11$	4	9.2405	483	3.2694 × 10−2	2.4086	206	3.2342 × 10−2
	8	20.9804	371	7.5679 × 10−3	5.6095	160	7.2467 × 10−3
	16	48.5232	260	1.9174 × 10−3	12.4909	113	1.6750 × 10−3
	32	111.5133	168	8.2150 × 10−4	28.1021	74	5.6681 × 10−4

. From these tables, it is observed that the EDGM and the CN–FDM can approximate the analytical solution accurately. The results also show that the EDGM outperforms the CN–FDM in terms of iteration counting and CPU timing as well. Figure 3a,b show 2D plots of CPU time when $\gamma (x,y,t)=\frac{1-{\left(xyt\right)}^3+sin{\left(xyt\right)}^2}{10}$ and $\frac{e^{xyt}+sin\left(xyt\right)}{30}$, respectively, whereas Figure 4a,b show 2D plots of iteration numbers when $\gamma (x,y,t)=\frac{1-{\left(xyt\right)}^3+sin{\left(xyt\right)}^2}{10}$ and $\frac{e^{xyt}+sin\left(xyt\right)}{30}$. From these figures, it is evident that the EDGM requires fewer iterations and consumes less computing time in handling the fractional problem compared to the CN–FDM. These improvements in time consumption and iteration count indicate that the EDGM is computationally superior to the CN–FDM. Figure 5 presents a graphical representation of the exact solution and the numerical solutions of the proposed methods at y = 0.5, $h=1/24$, $N=100$, $\gamma (x,y,t)=\frac{16-e^{xyt}}{17}$ and $T=0.25$, 0.5, 0.75 and 1. As seen in this figure, the numerical solutions match well with the exact solution of the problem under consideration. The behavior of the numerical errors using the CN–FDM and the EDGM is highlighted in Figure 6 at $T=1$, $h=1/24$, $N=100$ and $\gamma (x,y,t)=\frac{1-{\left(xyt\right)}^3+sin{\left(xyt\right)}^2}{10}$, from which it is indicated that the numerical solutions are in good agreement with the exact solution.

Example 2.

In this example, the following two-dimensional variable-order fractional cable equation is considered:

$$\begin{array}{cc}\hfill {}_0{}^{C}D_t^{\gamma (x,y,t)}u(x,y,t)& =\frac{{\mathsf{\partial }}^{2}u(x,y,t)}{\mathsf{\partial }{x}^2}+\frac{{\mathsf{\partial }}^{2}u(x,y,t)}{\mathsf{\partial }{y}^2}-u(x,y,t)\hfill \\ \hfill & +f(x,y,t),(x,y,t)\in \Omega \times [0,1],\hfill \end{array}$$

subject to the initial-boundary conditions,

$$\begin{array}{cc}\hfill u(x,y,0)& =0,(x,y)\in \Omega \hfill \\ \hfill u(0,y,t)& =t^3{(1-y^2)}^2,u(1,y,t)=0,(x,y,t)\in [0,1]\times \mathsf{\partial }\Omega ,\hfill \\ \hfill u(x,0,t)& =t^3{(1-x^2)}^2,u(x,1,t)=0,\hfill \end{array}$$

where $\Omega =(0,1)\times (0,1)$ with the source term defined as,

$$\begin{array}{cc}\hfill f(x,y,t)& =\frac{6t^{3-\gamma (x,y,t)}}{\Gamma (4-\gamma (x,y,t\left)\right)}{(1-x^2)}^2{(1-y^2)}^2-4t^3({(1-y^2)}^2(3x^2-1)\hfill \\ \hfill & +{(1-x^2)}^2(3y^2-1))+t^3{(1-x^2)}^2{(1-y^2)}^2.\hfill \end{array}$$

The exact analytical solution of the current problem is

$$u(x,y,t)=t^3{(1-x^2)}^2{(1-y^2)}^2.$$

Table 3. Comparison of numerical results of the proposed methods for Example 2 at $\tau =0.01$, $T=1$ and various choices of $\gamma (x,y,t)$.

		CN–FDM			EDGM
$\mathit{\gamma }(\mathit{x},\mathit{y},\mathit{t})$	${\mathit{h}}^{-\mathit{1}}$	$\mathit{CPU}$	$\mathit{Ite}$	$\mathit{MAE}$	$\mathit{CPU}$	$\mathit{Ite}$	$\mathit{MAE}$
$\frac{15-sin{\left(xyt\right)}^4}{16}$	6	1.1107	7	3.6260 × 10−3	0.8549	4	6.3340 × 10−3
	12	9.8215	14	9.1115 × 10−4	3.5592	8	1.4611 × 10−3
	18	41.471	26	3.1795 × 10−4	11.7159	13	6.2110 × 10−4
	24	109.8017	40	1.7590 × 10−4	32.605	19	2.9598 × 10−4
$\frac{9-y^5+t^3}{120}$	6	2.8476	19	3.9242 × 10−3	1.2287	10	6.7795 × 10−3
	12	34.0943	58	1.0029 × 10−3	9.621	28	1.5797 × 10−3
	18	148.4631	109	3.1527 × 10−4	41.7843	53	6.7007 × 10−4
	24	398.0526	168	2.3551 × 10−4	110.8827	83	3.2448 × 10−4
$\frac{15+sin\left(xt\right)}{16}$	6	1.0247	6	3.6211 × 10−3	0.7978	4	6.3249 × 10−3
	12	8.1283	14	9.2666 × 10−4	2.9754	7	1.4593 × 10−3
	18	32.5322	24	3.2975 × 10−4	9.9754	12	6.2310 × 10−4
	24	90.0535	38	1.7051 × 10−4	24.3636	18	3.0412 × 10−4
$\frac{11-cos{\left(yt\right)}^2}{11}$	6	1.0972	7	3.6316 × 10−3	0.7846	4	6.3415 × 10−3
	12	8.8651	15	9.1747 × 10−4	3.3896	8	1.4587 × 10−3
	18	38.1758	27	3.2260 × 10−4	11.0987	13	6.1325 × 10−4
	24	107.0003	42	1.8170 × 10−4	28.3911	20	2.9483 × 10−4
$5^{xt-3}$	6	2.9888	19	3.9502 × 10−3	1.4219	10	6.8113 × 10−3
	12	41.7307	59	1.0089 × 10−3	10.3176	29	1.5737 × 10−3
	18	180.8514	111	3.2972 × 10−4	51.3228	54	6.6835 × 10−4
	24	496.4411	170	2.3834 × 10−4	134.1644	85	3.4322 × 10−4
$\frac{10-{\left(xyt\right)}^3}{80}$	6	2.3473	19	3.9200 × 10−3	1.0168	9	6.7680 × 10−3
	12	37.1117	57	9.6831 × 10−4	10.5125	28	1.5614 × 10−3
	18	166.6839	108	3.2475 × 10−4	43.7882	52	6.7173 × 10−4
	24	462.7543	165	2.4535 × 10−4	130.107	82	3.1735 × 10−4

lists the numerical results of solving Example 2 using the EDGM and the CN–FDM for various choices of $\gamma (x,y,t)$, h and $\tau =0.01$. The corresponding results show that both EDGM and CN–FDM produce precise numerical solutions, which makes them feasible in terms of accuracy. It is well known that increasing the iteration’s number of a numerical algorithm will result in larger computational complexity, which ultimately leads to slower convergence. Figure 7 and Figure 8 establish a comparison between the EDGM and the CN–FDM in terms of CPU timing and iteration counting, respectively. One can see that the graphs of CPU timing are consistent with those of iteration counting. Obviously, the EDGM provides much efficient simulations in terms of CPU time as well as iteration numbers compared to the CN–FDM. Figure 9 demonstrates a comparison between the exact solution and numerical solutions of the EDGM and CN–FDM at $y=0.5$, $h=1/24$, $N=100$, $\gamma (x,y,t)=\frac{15+sin\left(xt\right)}{16}$ and $T=0.25$, 0.5, 0.75 and 1. Figure 10 sketches the maximum errors for the EDGM and the CN–FDM at $T=1$, $h=1/24$, $N=100$ and $\gamma (x,y,t)=\frac{10-{\left(xyt\right)}^3}{80}$. As these figures show, the proposed methods generate numerical solutions that are in good agreement with the exact solution of the given variable-order fractional problem.

In order to quantify the computational efficiency of the EDGM, the improvement percentages in aspects of CPU timing and iteration counting for solving Examples 1 and 2 are reported in Table 3 and

Table 4. The improvement percentages of the EDGM compared to the CN–FDM for various choices of $\gamma (x,y,t)$ in Example 1.

$\mathit{\gamma }(\mathit{x},\mathit{y},\mathit{t})$	CPU Timing	Iteration Counting
$\frac{1-xyt+sin\left(xyt\right)}{10}$	52.01–79.40%	3.91–54.67%
$\frac{5+{\left(xy\right)}^3-{\left(yt\right)}^4}{50}$	62.44–75.47%	50.00–55.13%
$\frac{3+{\left(xy\right)}^2-{\left(yt\right)}^3}{30}$	57.93–76.53%	50.00–55.13%
$\frac{16-e^{xyt}}{17}$	38.55–72.20%	44.44–54.10%
$\frac{e^{xyt}+sin\left(xyt\right)}{30}$	60.56–76.48%	52.00–54.61%
$\frac{1-{\left(xyt\right)}^3+sin{\left(xyt\right)}^2}{10}$	42.35–79.36%	50.00–54.67%

, respectively. For instance, the improvement percentages of CPU timing and iteration counting for Example 1 at $\gamma (x,y,t)=\frac{1-xyt+sin\left(xyt\right)}{10}$ indicate that the EDGM has reduced the CPU time approximately by 52.01–79.40% and the number of iterations by 53.91–54.67% compared to the CN–FDM. The improvement percentages of the EDGM for solving Examples 1 and 2 at other values of $\gamma (x,y,t)$ can be drawn similarly from Table 3 and Table 4.

7. Conclusions

In this paper, a point-wise CN–FDM and a group-wise EDGM are designed to deal with the numerical solution of the two-dimensional variable-order fractional cable equation. The CN–FDM and the EDGM are developed based on finite difference approximations derived on the standard and skewed grids, respectively. The stability and convergence of the proposed methods are proved via the technique of von Neumann analysis. The proposed numerical schemes are implemented, and the computational outputs are presented in tabular and graphical forms, from which we observe that the numerical solutions are in good agreement with the exact analytical solutions (see Table 1, Table 2 and Table 3 and Figure 5, Figure 6, Figure 9 and Figure 10). Furthermore, as can be clearly seen from Table 4 and

Table 5. The improvement percentages of the EDGM compared to the CN–FDM for various choices of $\gamma (x,y,t)$ in Example 2.

$\mathit{\gamma }(\mathit{x},\mathit{y},\mathit{t})$	CPU Timing	Iteration Counting
$\frac{15-sin{\left(xyt\right)}^4}{16}$	23.03–71.75%	42.86–52.50%
$\frac{9-y^5+t^3}{120}$	56.85–72.14%	47.37–51.72%
$\frac{15+sin\left(xt\right)}{16}$	22.14–72.95%	33.33–52.63%
$\frac{11-cos{\left(yt\right)}^2}{11}$	28.49–73.47%	42.85–52.38%
$5^{xt-3}$	49.41–75.28%	47.36–51.35%
$\frac{10-{\left(xyt\right)}^3}{80}$	56.68–73.73%	50.30–52.63%

and Figure 3, Figure 4, Figure 7 and Figure 8, the EDGM is shown to require cheap computational complexity in terms of CPU timing and iteration counting compared to the CN–FDM. In fact, the obtained numerical results in Table 1 and Table 3 reveal that the EDGM achieves improved percentages 22.14–79.40% in CPU timing and 42.85–55.13% in iteration counting as compared to the CN–FDM. As future works, the EDGM can be constructed based on higher-order discrete schemes to achieve better accuracy. In addition, the EDGM can be implemented on parallel computers and extended to solve other types of high-dimensional variable-order fractional problems.