Numerical Finite Difference Scheme for a Two-Dimensional Time-Fractional Semilinear Diffusion Equation

Fractional calculus, a branch of mathematics concerned with the properties of derivatives and integrals of non-integer orders, has been an area of interest since the inception of classical calculus. This mathematical field specializes in the resolution of time-dependent fractional differential equations, which entail fractional derivatives.

Regarded as a pivotal tool for the exploration of dynamical systems, fractional calculus is valued for its non-local operators that encapsulate the historical progression of dynamics. The usage of fractional calculus and fractional processes has become a favored approach when grappling with the unique properties of long-term memory effects found in multiple domains of applied sciences. These fields range from finance and economic dynamics, biological systems and bioinformatics, nonlinear waves and acoustics, to image and signal processing, transportation systems, geosciences, astronomy, and cosmology.

The development of fractional calculus has been incrementally built upon by an array of researchers including, but not limited to, Heaviside, Lagrange Riemann, Liouville, Gr¨unwald, Euler, Fourier, and Abel [1]. The contemporary utilization of fractional calculus owes its popularity to the versatility of the differintegral operator, which amalgamates both integer-order derivatives and integrals as special cases. The fractional integral, for instance, may be employed to describe the cumulation of a quantity when the order of integration is undetermined, and can be inferred as a parameter of a regression model, as elucidated by Podlubny [2] and Kisela [3]. The fractional derivative, on the other hand, is frequently used to depict damping. Other applications span control theory of dynamical systems, optics and signal processing, fluid flow, diffusive transport, probability and statistics, viscoelasticity, electrical networks, dynamical processes in self-similar and porous structures, electrochemistry of corrosion, and rheology.

Over the past decades, several analytical and numerical methods for solving fractional differential equations have been brought to bear. Analytical methods encompass the Fourier transformation method, the Laplace transformation method, and the green function method [4-6]. However, most fractional differential equations resist analytical solutions. Hence, the development of numerical schemes for these equations is imperative. Effective methods such as the finite difference method, the spectral method, and the finite element method have been utilized to solve time (space) fractional differential equations [7-9].

The numerical solutions of linear types of time (space)-fractional equations have been the focal point of many authors [9-13], while semilinear or nonlinear types have received less attention [14-17]. Nevertheless, research on these latter types is still in its nascent stages.

Since, most of the mathematical models, that describe real would phenomena, involving non-linear fractional partial differential equations, and most of previous studies have focused on linear types, the aim of this work is to suggest a high order efficient numerical technique for a two-dimensional semi-linear time fractional equation with Dirichlet boundary conditions. Namely, we propose the linearly implicit Crank-Nicolson finite difference scheme, and we show that the proposed scheme is consistent, conditionally stable and convergent with the help of mathematical induction. Moreover, we prove that the order of accuracy takes the form: $O\left(k^{2-\alpha}+h^2\right)$, where $\alpha$ is the fractional order of the time-derivative. In addition, the obtained numerical results of the studied test examples show that the proposed scheme can be used efficiently to compute the numerical solutions of the governed problem.

The rest of this paper is organized as follows: In section two, we state the mathematical formulation of the time fractional diffusion equation and present some literature review regarding the studied problem. In section three the derivation, consistency, stability and convergent of the finite difference scheme are considered in section four. Finally, some conclusions are given in the last section.

In this section, we propose a numerical approximation to problem (2). Namely, the linearly Crank-Nicolson method (CNM). Firstly, we derive its formula, then we investigate its consistency, stability and convergence.

In this segment, the grid dimensions in relation to space and time for the positive integers I and N are respectively represented by $h=\frac{1}{I}$ and $k=\frac{T}{N}$. The grid point in the space interval [0, 1] is denoted $x_i=i h, y_j=j h, i, j=0,1,2, \ldots, I$ and the grid points for time are designated $t_n=n k, n=0,1,2 \ldots . . N$. In addition, we denote $u_i^n=u\left(x_i, y_j, t_n\right)$.

3.1 Crank Nicolson formula

The goal of this section is to derive Crank-Nicolson scheme for problem (2), which is one of the most popular methods in practice. Moreover, it has high order of convergence in both space and time.

In order to approximate Eq. (2) and Eq. (4) at the mush point $\left(x_i, y_j, t_{n+\frac{1}{2}}\right)$, we use the following approximations [10, 13]:

$\left.\frac{\partial^\alpha u}{\partial t^\alpha}\right|_i ^{n+\frac{1}{2}}=\left[\widetilde{w}_1 u_{i, j}^n+\sum_{s=1}^{n-1}\left(\widetilde{w}_{n-s+1}-\widetilde{w}_{n-s}\right) u_{i, j}^s-\widetilde{w}_n u_{i, j}^0+\sigma\left(\frac{u_{i, j}^{n+1}+u_{i, j}^n}{2^{1-\alpha}}\right)\right]+O\left(k^{(2-\alpha)}\right)$                     (5)

$\left.\frac{\partial^2 u}{\partial x^2}\right|_i ^{n+\frac{1}{2}}=\left[\left(\frac{u_{i+1, j}^{n+1}\,\,+\,\,-\,\,2 u_{i, j}^{n+1}\,\,+u_{i-1, j}^{n+1}}{2 h^2}\quad \right)\right]+\left[\left(\frac{u_{i+1, j}^n\,\,+\,\,-\,\,2 u_{i, j}^n\,\,+u_{i-1, j}^n}{2 h^2}\quad\right)\right]+O\left(h^2\right)$                     (6)

$\begin{gathered}\left.\frac{\partial^2 u}{\partial y^2}\right|_i ^{n+\frac{1}{2}}=\left[\left(\frac{u_{i, j+1}^{n+1}+-2 u_{i, j}^{n+1}+u_{i, j,-1}^{n+1}}{2 h^2}\right)\right]+\left[\left(\frac{\left.u_{i, j+1}^n+-2 u_{i, j}^n+u_{i, j-1}^n\right)}{2 h^2}\right)\right]+O\left(h^2\right) \\ \sigma=\frac{1}{k^\alpha \Gamma(2-\alpha)}, \widetilde{w}_s=\sigma\left[\left(s+\frac{1}{2}\right)^{(1-\alpha)}-\left(s-\frac{1}{2}\right)^{(1-\alpha)}\right]\end{gathered}$                        (7)

The nonlinear term can be approximated using Taylor expansion [16]:

$f\left(x_i, y_j, t_{n+\frac{1}{2}}, u\left(x_i, y_j, t_{n+\frac{1}{2}}\right)\right)=f\left(x_i, y_j, t_{n+\frac{1}{2}}, \frac{3}{2} u\left(x_i, y_j, t_n\right)-\frac{1}{2} u\left(x_i, y_j, t_{n-1}\right)\right)+O\left(k^2\right)$

Thus

$f\left(x_i, y_j, t_{n+\frac{1}{2}}, u\left(x_i, y_j, t_{n+\frac{1}{2}}\right)\right)=f\left(x_i, y_j, t_{n+\frac{1}{2}}, \frac{3}{2} u_{i, j}^n-\frac{1}{2} u_{i, j}^{n-1}\right)+O\left(k^2\right)$                     (8)

By substituting the above forms (5)-(8) in problem (2), we can get the Crank-Nicolson approximated formula as follows:

$\begin{gathered}{\left[\widetilde{w}_1 u_{i, j}^n+\sum_{s=1}^{n-1}\left(\widetilde{w}_{n-s+1}-\widetilde{w}_{n-s}\right) u_{i, j}^s-\widetilde{w}_n u_{i, j}^0+\sigma\left(\frac{u_{i, j}^{n+1}-u_{i, j}^{n}}{2^{1-\alpha}}\right)\right]} \\ =\frac{1}{2 h^2}\left[\left(u_{i+1, j}^{n+1}+-2 u_{i, j}^{n+1}+u_{i-1, j}^{n+1}\right)+\left(u_{i+1, j}^n+-2 u_{i, j}^n+u_{i-1, j}^n\right)\right] \\ +\frac{1}{2 h^2}\left[\left(u_{i, j+1}^{n+1}+-2 u_{i, j}^{n+1}+u_{i, j-1}^{n+1}\right)+\left(u_{i, j+1}^n+-2 u_{i, j}^n+u_{i, j-1}^n\right)\right]+f\left(x_i, y_j, t_{n+\frac{1}{2}, }\frac{3}{2} u_{i, j}^n-\frac{1}{2} u_{i, j}^{n-1}\right) .\end{gathered}$

Multiplying both sides of the last equation by $\delta=k^\alpha \Gamma(2-\alpha) 2^{1-\alpha}$ , yields that:

$\begin{gathered}{\left[2^{1-\alpha} w_1 u^n+2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s+1}-w_{n-s}\right) u_{i, j}^s-2^{1-\alpha} w_n u_{i, j}^0+u_{i, j}^{n+1}-u_{i, j}^n\right]} \\ =\frac{\delta}{2 h^2}\left[\left(u_{i+1, j}^{n+1}+-2 u_{i, j}^{n+1}+u_{i-1, j}^{n+1}\right)+\left(u_{i+1, j}^n+-2 u_{i, j}^n+u_{i-1, j}^n\right)\right] \\ +\frac{\delta}{2 h^2}\left[\left(u_{i, j+1}^{n+1}+-2 u_{i, j}^{n+1}+u_{i, j-1}^{n+1}\right)+\left(u_{i, j+1}^n+-2 u_{i, j}^n+u_{i, j-1}^n\right)\right]+\delta f\left(x_i, y_j, t_{n+\frac{1}{2}}, \frac{3}{2} u_{i, j}^n-\frac{1}{2} u_{i, j}^{n-1}\right) .\end{gathered}$

So,

$\begin{gathered}u_{i, j}^{n+1}-\frac{\delta}{2 h^2}\left(u_{i+1, j}^{n+1}++u_{i-1, j}^{n+1}+u_{i, j+1}^{n+1}+u_{i, j-1}^{n+1}\right)+\frac{2 \delta}{h^2} u_{i, j}^{n+1}=\frac{\delta}{2 h^2}\left(u_{i+1, j}^n+u_{i-1, j}^n+u_{i, j+1}^n+u_{i, j-1}^n\right) \\ -\frac{2 \delta}{h^2} u_{i, j}^n+u_{i, j}^n-2^{1-\alpha} w_1 u_{i, j}^n+2^{1-\alpha} \sum_{s=1}^n\left(w_{n-s}-w_{n+1-s}\right) u_{i, j}^s+2^{1-\alpha} w_n u_{i, j}^0+\delta f\left(x_i, y_j, t_{n+\frac{1}{2}, }\frac{3}{2} u_{i, j}^n-\frac{1}{2} u_{i, j}^{n-1}\right) .\end{gathered}$

It follows

$\begin{gathered}(1+2 r) u_{i, j}^{n+1}-\frac{r}{2}\left(u_{i+1, j}^{n+1}++u_{i-1, j}^{n+1}+u_{i, j+1}^{n+1}+u_{i, j-1}^{n+1}\right)=\frac{r}{2}\left(u_{i+1, j}^n+u_{i-1, j}^n+u_{i, j+1}^n+u_{i, j-1}^n\right) \\ +\left(1-2 r-2^{1-\alpha} w_1\right) u_{i, j}^n+2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n+1-s}\right) u_{i, j}^s+2^{1-\alpha} w_n u_{i, j}^0+\delta f\left(x_i, y_j, t_{n+\frac{1}{2}} ,\frac{3}{2} u_{i, j}^n-\frac{1}{2} u_{i, j}^{n-1}\right) .\end{gathered}$

We obtain

$\begin{gathered}(1+2 r) u_{i, j}^{n+1}-\frac{r}{2}\left[u_{i+1, j}^{n+1}+u_{i-1, j}^{n+1}+u_{i, j+1}^{n+1}+u_{i, j-1}^{n+1}\right]=\left(1-2^{1-\alpha} w_1-2 r\right) u_{i, j}^n+\frac{r}{2}\left[u_{i+1, j}^n+u_{i-1, j}^n+u_{i, j+1}^n+u_{i, j-1}^n\right] \\ +2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right) u_{i, j}^s+2^{1-\alpha} w_n u_{i, j}^0+\delta f\left(x_i, y_j, t_{n+\frac{1}{2}}, \frac{3}{2} u_{i, j}^n-\frac{1}{2} u_{i j}^{n-1}\right)\end{gathered}$                         (9)

where, $\delta=k^\alpha \Gamma(2-\alpha) 2^{1-\alpha}, r=\frac{\delta}{h^2}, w_s=\left[\left(s+\frac{1}{2}\right)^{(1-\alpha)}-\left(s-\frac{1}{2}\right)^{(1-\alpha)}\right]$

For n=0, CN formula (9) becomes as follows:

$\begin{gathered}(1+2 r) u_{i, j}^1-\frac{r}{2}\left[u_{i+1, j}^1+u_{i-1, j}^1+u_{i, j+1}^1+u_{i, j-1}^1\right] \\ =(1-2 r) u_{i, j}^0+\frac{r}{2}\left[u_{i+1, j}^0+u_{i-1, j}^0+u_{i, j+1}^0+u_{i, j-1}^0\right]+\delta f\left(x_i, y_j, t_{n+\frac{1}{2}, }\frac{3}{2} u_{i, j}^0\right) .\end{gathered}$

Lemma 3.1 [13]

For s = 0,1, 2…, n

$\cdot \quad w_s>0$

$\cdot \quad w_{s+1}<w_s, w_{s+1}^{-1}>w_s^{-1}$

$\cdot \quad \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right)+w_n=w_1$.

The Crank – Nicolson difference formula (5) can be written in matrix form as follows:

$\left\{\begin{array}{l} A U^1=B U^0+\delta F^0+Z^1+Z^0 \\ A U^{n+1}=B U^n+2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right) U^s +2^{1-\alpha} w_n U^0+\delta F^n+Z^{n+1}+Z^n\end{array}\right\}$                    (10)

$A=\left[\begin{array}{ccccc}A^* & -\frac{r}{2} I & 0 & \cdots & 0 \\ -\frac{r}{2} I & A^* & -\frac{r}{2} I & \cdots & 0 \\ & & \ddots & & \\ & & &\ddots & \\ 0 & \cdots & 0 & -\frac{r}{2} I & A^*\end{array}\right]_{(I-1)^2 \times(I-1)^2}$

$B=\left[\begin{array}{ccccc}B^* & \frac{r}{2} I & 0 & \cdots & 0 \\ \frac{r}{2} I & B^* & \frac{r}{2} I & \cdots & 0 \\ & & \ddots & & \\ & & &\ddots & \\ 0 & \cdots & 0 & \frac{r}{2} I & B^*\end{array}\right]_{(I-1)^2 \times(I-1)^2}$

$A^*=\left[\begin{array}{ccccc} (1+2r) & -\frac{r}{2} & 0 & \cdots & 0 \\ -\frac{r}{2} & (1+2r) & -\frac{r}{2} & \cdots & 0 \\ & & \ddots & & \\ & & &\ddots & \\ 0 & \cdots & 0 & -\frac{r}{2} & (1+2r)\end{array}\right]_{(I-1) \times(I-1)}$

$B^*=\left[\begin{array}{ccccc} (1-2r-2^{1-\alpha}w_1) & \frac{r}{2} & 0 & \cdots & 0 \\ \frac{r}{2} & (1-2r-2^{1-\alpha}w_1) & \frac{r}{2} & \cdots & 0 \\ & & \ddots & & \\ & & &\ddots & \\ 0 & \cdots & 0 & \frac{r}{2} & (1-2r-2^{1-\alpha}w_1)\end{array}\right]_{(I-1) \times(I-1)}$

$U^n=\left[\begin{array}{c}u_1^n \\ u_2^n \\ \vdots \\ u_{I-1}^n\end{array}\right], F^n=\left[\begin{array}{c}f_1^n \\ f_2^n \\ \vdots \\ f_{I-1}^n\end{array}\right]$,

$\begin{aligned} & u_i^n=\left[\begin{array}{llll}u_{i, 1}^n & u_{i, 2}^n & \cdots & u_{i, I-1}^n\end{array}\right]^T, \\& f_i^n=\left[\begin{array}{llll}f_{i, 1}^n & f_{i, 2}^n & \cdots & f_{i, I-1}^n\end{array}\right]^T \\ & f_{i, j}^n=f\left(x_i, y_j, t_{n \frac{1}{2}}, \frac{3}{2} u_{i, j}^n-\frac{1}{2} u_{i j}^{n-1}\right) . \\ & \end{aligned}$

$Z^n=\frac{r}{2}\left[u_{0,1}^n+u_{1,0}^n, u_{0,2}^n, \ldots, u_{0, I-2}^n \,\,\,, u_{0, I-1}^n+u_{1, I}^n ; u_{2,0}^n, 0, \ldots, \right.$

$\left. 0,u_{2, I}^n; \ldots; u_{I, 1}^n+u_{I-1,0}^n \,\,\,, u_{I, 2}^n, \ldots ., u_{I, I-2}^n \,\,\,, u_{I, I-1}^n+u_{I-1, I}^n\right]^T$

Theorem 3.1 At each time level, (n+1), the linear system (10) is uniquely salable.

Proof: Since A is diagonally dominant with positive real diagonal entries, then A is positive definite and nonsingular [26]. Hence, the linear system (10) is uniquely solvable.

3.2 Stability analysis for C.N. method

Suppose that $\tilde{u}_i^n$ is the approximate solution of Eq. (2). As in studies [18, 27], we set $e_{i, j}^n=\tilde{u}_{i, j}^n-u_{i, j}^n,$$i, j=0,1,2, \ldots, I-1,$ $n=0,1,2, \ldots, N$.

Define $\|E^n\|=\operatorname{Max}_{\substack{1 \leq i \leq I-1 \\ 1 \leq j \leq I-1}} \,\,\left|e_{i,j}^n\right|$

Clearly, $e_{i,j}^n$  satisfies (9):

For n=0,

$\begin{gathered}(1+2 r) e_{i, j}^1-\frac{r}{2}\left[e_{i+1, j}^1+e_{i-1, j}^1+e_{i, j+1}^1+e_{i, j-1}^1\right]=(1-2 r) e_{i, j}^0- \\ \frac{r}{2}\left[e_{i+1, j}^0+e_{i-1, j}^0+e_{i, j+1}^0+e_{i, j-1}^0\right]+\delta\left[f\left(x_i, y_j, t_{\frac{1}{2}}, \frac{3}{2} \tilde{u}_{i, j}^0\right)-f\left(x_i, y_j, t_{\frac{1}{2}}, \frac{3}{2} u_{i, j}^0\right)\right]\end{gathered}$                          (11)

For n>0,

$\begin{gathered}(1+2 r) e_{i, j}^{n+1}-\frac{r}{2}\left[e_{i+1, j}^{n+1}+e_{i-1, j}^{n+1}+e_{i, j+1}^{n+1}+e_{i, j-1}^{n+1}\right] \\ =\left(1-2^{1-\alpha} w_1-2 r\right) e_{i, j}^n+\frac{r}{2}\left[e_{i+1, j}^n+e_{i-1, j}^n+e_{i, j+1}^n+e_{i, j-1}^n\right] \\ +2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right) e_{i, j}^n+2^{1-\alpha} w_n e_{i, j}^0 \\ +\delta\left[f\left(x_i, y_j, t_{n+\frac{1}{2}}, \frac{3}{2} \tilde{u}_{i, j}^n-\frac{1}{2} \tilde{u}_{i, j}^{n-1}\right)-f\left(x_i, y_j, t_{n+\frac{1}{2}}, \frac{3}{2} u_{i, j}^n-\frac{1}{2} u_{i, j}^{n-1}\right)\right] .\end{gathered}$                        (12)

Based on (10), the above error formula can be written in matrix form as follows:

$\left\{\begin{array}{l}A E^1=B E^0+\delta G^0 \\ A E^{n+1}=B E^n+2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right) E^s +2^{1-\alpha} w_n E^0+\delta G^n\end{array}\right\}$

E⁰ is given,

$E^n=\left[\begin{array}{c}E_1^n \\ E_2^n \\ \vdots \\ E_{I-1}^n\end{array}\right], G^n=\left[\begin{array}{c}G_1^n \\ G_2^n \\ \vdots \\ G_{I-1}^n\end{array}\right]$,

$E_i^n=\left[\begin{array}{llll}e_{i, 1}^n & e_{i, 2}^n & \cdots & e_{i, I-1}^n\end{array}\right]^T, \quad G_i^n=\left[\begin{array}{llll}g_{i, 1}^n & g_{i, 2}^n & \cdots & g_{i, I-1}^n\end{array}\right]^T$

$g_{i, j}^n=f\left(x_i, y_j, t_{n \frac{1}{2}}, \frac{3}{2} \tilde{u}_{i, j}^n-\frac{1}{2} \tilde{u}_{i j}^{n-1}\right)-f\left(x_i, y_j, t_{n \frac{1}{2}}, \frac{3}{2} u_{i, j}^n-\frac{1}{2} u_{i j}^{n-1}\right)$.

$i, j=0,1,2, \ldots, I-1, n=0,1,2, \ldots, N$

Definition 3.1 [14] For any arbitrary initial rounding error $E^0$, the difference approximation $A U^{n+1}=B U^n+b^n$ is stable, if there exists C>0, independent on h, k such that $\left\|E^n\right\| \leq C\left\|E^0\right\|$, or $\left\|\left(A^{-1} B\right)^n\right\| \leq C, \forall n$.

Theorem 3.2 The C.N. finite difference approximation is stable, if $\left(1-2^{1-\alpha} w_1-2 r\right) \geq 0$.

Proof: To prove this theorem, we use the Mathematical induction.

For n=0, set $\left\|E^1\right\|_{\infty}=\operatorname{Max}_{\substack{1 \leq i \leq I-1 \\ 1 \leq j \leq I-1}}\,\,\,\left|e_{i, j}^1\right|=\left|e_{p, q}^1\right|$, we have

$\begin{gathered}\left|e_{p, q}^1\right|=(1+2 r)\left|e_{p, q}^1\right|-\frac{r}{2}\left(4\left|e_{p, q}^1\,\,\right|\right) \\ \leq(1+2 r)\left|e_{p, q}^1\,\,\right|-\frac{r}{2}\left(\left|e_{p+1, q}^1\,\,\right|+\left|e_{p-1, q}^1\,\,\right|+\left|e_{p, q+1}^1\,\,\right|+\left|e_{p, q-1}^1 \,\,\right|\right) \\ \leq\left|(1+2 r) e_{p, q}^1-\frac{r}{2}\left(e_{p+1, q}^1+e_{p-1, q}^1+e_{p, q+1}^1+e_{p, q-1}^1\right)\right|\end{gathered}$                         (13)

From (11) and (13), we obtain

$\begin{aligned}\left|e_{p, q}^1\right| \leq & \left|(1-2 r) e_{p, q}^0+\frac{r}{2}\left(e_{p+1, q}^0+e_{p-1, q}^0+e_{p, q+1}^0+e_{p, q-1}^0\right)\right| \\ & +\delta\left|f\left(x_p, y_q, t_{\frac{1}{2}}, \frac{3}{2} \widetilde{u}_{p, q}^0\right)-f\left(x_p, y_q, t_{\frac{1}{2}}, \frac{3}{2} u_{p, q}^0\right)\right|\end{aligned}$

Since satisfies Lipchitz condition (3), we obtain

$\begin{gathered}e p,q 1 \leq 1-2 r E 0 \infty+2 r E 0 \infty\left|e_{p, q}^1\right| \leq(1-2 r)\left\|E^0\right\|_{\infty}+2 r\left\|E^0\right\|_{\infty}+\delta L\left|\frac{3}{2}\left(\tilde{u}_{p, q}^0-u_{p, q}^0\right)\right|. \\ \leq\left\|E^0\right\|_{\infty}+\frac{3}{2} \delta L\left\|E^0\right\|_{\infty}=\left(1+\frac{3}{2} \delta L\right)\left\|E^0\right\|_{\infty} \leq(1+2 \delta L)\left\|E^0\right\|_{\infty} .\end{gathered}$

Thus $\left\|E^1\right\|_{\infty} \leq(1+2 \delta L)\left\|E^0\right\|_{\infty}$

Now, supposes that $\left\|E^s\right\|_{\infty} \leq C\left\|E^0\right\|_{\infty}$, where $C=(1+2 \delta L), s=0,1,2, \ldots, n$, for $n+$, let $\left|e_{p, q}^{n+1}\right|=\operatorname{Max}_{\substack{1 \leq i \leq I-1 \\ 1 \leq j \leq I-1}}\left|e_{i, j}^{n+1}\,\,\right|$ we have

$\begin{gathered}\left|e_{p, q}^{n+1}\right|=(1+2 r)\left|e_{p, q}^{n+1}\right|-\frac{r}{2}\left[4\left|e_{p, q}^{n+1}\right|\right] . \\ \leq(1+2 r)\left|e_{p, q}^{n+1}\right|-\frac{r}{2}\left(\left|e_{p+1, q}^{n+1}\,\,\right|+\left|e_{p-1, q}^{n+1}\,\,\right|+\left|e_{p, q+1}^{n+1}\,\,\right|+\left|e_{p, q-1}^{n+1}\,\,\right|\right) . \\ \leq\left|(1+2 r) e_{p, q}^{n+1}-\frac{r}{2}\left(e_{p+1, q}^{n+1}+e_{p-1, q}^{n+1}+e_{p, q+1}^{n+1}+e_{p, q-1}^{n+1}\right)\right| .\end{gathered}$                          (14)

From (12) and (14), it follows that

$\begin{gathered}\left|e_{p, q}^{n+1}\right| \leq\left|\left(1-2^{1-\alpha} w_1-2 r\right) e_{p, q}^n+\frac{r}{2}\left(e_{p+1, q}^n+e_{p-1, q}^n+e_{p, q+1}^n+e_{p, q-1}^n\right)\right|+\left|2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right) e_{p, q}^s\right|+ \\ \delta\left|f\left(x_p, y_q, t_{n+\frac{1}{2}}, \frac{3}{2} \tilde{u}_{p, q}^n-\frac{1}{2} \tilde{u}_{p, q}^{n-1}\right)-f\left(x_p, y_q, t_{n+\frac{1}{2}}, \frac{3}{2} u_{p, q}^n-\frac{1}{2} u_{p, q}^{n-1}\right)\right|\end{gathered}$ .

Since satisfies Lipchitz condition (3), we obtain

$\begin{gathered}\left|e_{p, q}^{n+1}\right| \leq\left(1-2^{1-\alpha} w_1-2 r\right)\left|e_{p, q}^n\right|+\frac{r}{2}\left(\left|e_{p+1, q}^n\right|+\left|e_{p-1, q}^n\right|+\left|e_{p, q+1}^n\right|+\left|e_{p, q-1}^n\right|\right)+2^{1-\alpha} w_n\left|e_{p, q}^0\right|+ \\ 2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right)\left|e_{i, j}^s\right|+\delta L\left|\frac{3}{2}\left(\tilde{u}_{p, q}^n-u_{p, q}^n\right)+\frac{1}{2}\left(u_{p, q}^{n-1}-\tilde{u}_{p, q}^{n-1}\right)\right| .\end{gathered}$

$\begin{gathered}\leq\left(1-2^{1-\alpha} w_1-2 r\right)\left\|E^n\right\|_{\infty}+2 r\left\|E^n\right\|_{\infty}+2^{1-\alpha} w_n\left\|E^0\right\|_{\infty}+2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right)\left\|E^s\right\|_{\infty}+ \\ \delta L\left[\frac{3}{2}\left|e_{p, q}^n\right|+\frac{1}{2}\left|e_{p, q}^{n-1}\right|\right] .\end{gathered}$

$\begin{gathered}\leq\left(1-2^{1-\alpha} w_1\right) C^n\left\|E^n\right\|_{\infty}+2^{1-\alpha} w_n\left\|E^0\right\|+2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right) C^s\left\|E^0\right\|+\delta L\left[\frac{3}{2} C^n\left\|E^0\right\|+\frac{1}{2} C^{n-1}\left\|E^0\right\|\right] .\\ \quad \leq\left(1-2^{1-\alpha} w_1\right) C^n\left\|E^0\right\|+2^{1-\alpha} w_n C^n\left\|E^0\right\|+2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right) C^n\left\|E^0\right\|+2 \delta L \,\,C^n\left\|E^0\right\|\end{gathered}$

$\begin{gathered}=\left[1-2^{1-\alpha} w_1+2^{1-\alpha} w_n+2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right)+2 \delta L\right] C^n\left\|E^0\right\| . \\ =(1+2 \delta L) C^n\left\|E^0\right\|=(1+2 \delta L)^{n+1}\left\|E^0\right\| .\end{gathered}$

Thus

$\begin{aligned} & \left\|E^{n+1}\right\|_{\infty} \leq(1+2 \delta L)^{n+1}\left\|E^0\right\|. \\ \text { So, }\left\|E^{n+1}\right\|& \leq\left(1+2^{2-\alpha} k^\alpha \Gamma(2-\alpha) L\right)^{n+1}\left\|E^0\right\|_{\infty} \\ &\leq\left(1+2^{2-\alpha} T^\alpha \Gamma(2-\alpha) L\right)^{n+1}\left\|E^0\right\|_{\infty}\end{aligned}$

Hence, there exists C>0, such that $\left\|E^{n+1}\right\| \leq C\left\|E^0\right\|_{\infty}$ . 

3.3 Convergence analysis of Crank – Nicolson method

Let $u_{i, j}^n=u\left(x_i, y_j, t_n\right)$ be the exact solution of problem (2) at mesh point $\left(x_i, y_j, t_n\right), i, j=0,1,2, \ldots, I-1, n=0,1,2, \ldots, N$.

Definition 3.2

$e_{i, j}^n=u\left(x_i, y_j, t_n\right)-u_{i, j}^n, n=0,1,2, \ldots, N, i, j=0,1,2, \ldots ., I-1$

$\begin{aligned} E^n= & {\left[\begin{array}{c}E_1^n \\ E_2^n \\ \vdots \\ E_{I-1}^n\end{array}\right], E^0=\left[\begin{array}{c}0 \\ 0 \\ \vdots \\ 0\end{array}\right], E_i^n=\left[\begin{array}{llll}e_{i, 1}^n & e_{i, 2}^n & \cdots & e_{i, I-1}^n\end{array}\right]^T } \\ & i=0,1,2, \ldots, I-1, n=0,1,2, \ldots, N\end{aligned}$

By substitution $e_{i, j}^n$  in the C.N. Eq. (9), it follows that

For n=0,

$\begin{gathered}(1+2 r) e_{i, j}^1-\frac{r}{2}\left[e_{i+1, j}^1+e_{i-1, j}^1+e_{i, j+1}^1+e_{i, j-1}^1\right] \\ =\delta\left[f\left(x_i, y_j, t_{\frac{1}{2}}, \frac{3}{2} u\left(x_i, y_j, \frac{3}{2} u\left(x_i, y_j, t_0\right)\right)-f\left(x_i, y_j, t_{\frac{1}{2}}, \frac{3}{2} u_{i, j}^0\right)\right]+T_{i, j}^{\frac{1}{2}}\right.\end{gathered}$                    (15)

For n>0,

$\begin{gathered}(1+2 r) e_{i, j}^{n+1}-\frac{r}{2}\left[e_{i+1, j}^{n+1}+e_{i-1, j}^{n+1}+e_{i, j+1}^{n+1}+e_{i, j-1}^{n+1}\right]=\left(1-2^{1-\alpha} w_1-2 r\right) e_{i, j}^n+ \\ \frac{r}{2}\left[e_{i+1, j}^n+e_{i-1, j}^n+e_{i, j+1}^n+e_{i, j-1}^n\right]+2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right) e_{i, j}^s+ \\ \delta\left[f\left(x_i, y_j, t_{n+\frac{1}{2}}, \frac{3}{2} u\left(x_i, y_j, t_n\right)-\frac{1}{2} u\left(x_i, y_j, t_{n-1}\right)\right)-f\left(x_i, y_j, t_{n+\frac{1}{2}}, \frac{3}{2} u_{i, j}^n-\frac{1}{2} u_{i, j}^{n-1}\right)\right]+T_{i, j}^{n+1}\end{gathered}$                (16)

Theorem 3.3 There exists C>0 such that:

$\left|T_{i, j}\right| \leq C\left(k^{2-\alpha}+h^2\right)$

Proof: Let $u_{i, j}^n=u\left(x_i, y_j, t_n\right)$ be the exact solution of problem (2).

$\begin{aligned} T_{i, j}^{n+\frac{1}{2}}=\left(u_{i, j}^{n+1}-u_{i, j}^n\right)+ & 2^{1-\alpha} w_1 u_{i, j}^n-2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right) u_{i, j}^s-2^{1-\alpha} w_n \\ -\frac{r}{2}\left[u_{i+1, j}^{n+1}+u_{i-1, j}^{n+1}+u_{i, j+1}^{n+1}+\right. & \left.u_{i, j-1}^{n+1}-4 u_{i, j}^{n+1}\right]-\frac{r}{2}\left[u_{i+1, j}^n+u_{i-1, j}^n+u_{i, j+1}^n+u_{i, j-1}^n-4 u_{i, j}^n\right] \\ - & \delta f\left(x_i, y_j, t_{n+\frac{1}{2}}, \frac{3}{2} u_{i, j}^n-\frac{1}{2} u_{i, j}^{n-1}\right) .\end{aligned}$

Thus

$\begin{gathered}k^\alpha \Gamma(2-\alpha) 2^{1-\alpha}\left[u_t^\alpha\left(x_i, y_j, t_{n+\frac{1}{2}}\right)+O\left(k^{2-\alpha}\right)\right]-k^\alpha \Gamma(2-\alpha) 2^{1-\alpha}\left[u_{x x}\left(x_i, y_j, t_{n+1}\right)+u_{x x}\left(x_i, y_j, t_n\right)+O\left(h^2\right)\right] \\ -k^\alpha \Gamma\left[(2-\alpha) 2^{1-\alpha} f\left(x_i, y_j, t_{n+\frac{1}{2}}, u_{i, j}^{n+\frac{1}{2}}\right)+O\left(k^2\right)\right] \\ =k^\alpha \Gamma(2-\alpha) 2^{1-\alpha}\left[u_t-u_{x x}-u_{y y}-f\right]_{i, j}^{n+\frac{1}{2}}+k^\alpha \Gamma(2-\alpha) 2^{1-\alpha}\left[O\left(k^{2-\alpha}\right)+O\left(h^2\right)+O\left(k^2\right)\right]\end{gathered}$

Thus $\left|T_{i, j}^{n+\frac{1}{2}}\right| \leq C\left(k^{2-\alpha}+h^2\right), C>0$.

Remark 3.1 Based on Theorem 3.7, the C.N. formula (9) is consistent.

$\left|T_{i, j}^{n+\frac{1}{2}}\right| \rightarrow 0, h, k \rightarrow 0$.

Theorem 3.4 If $\left(1-2 r-2^{1-\alpha} w_1\right) \geq 0$, there exists $C>0$ such that

$|e_{i, j}^{n+1}| \leq C\left(k^{2-\alpha}+h^2\right), i, j=1,2, \ldots \ldots, I-1, \quad n=0,1, \ldots \ldots \ldots, N$

Proof: Define $\left\|E^n\right\|_{\infty}=\operatorname{Max}_{\substack{1 \leq i \leq I-1 \\ 1 \leq j \leq I-1}}\,\,\left|e_{i, j}^n\right|$.

By using mathematical induction, we can prove this theorem as follows:

For n=0, let $\left\|E^1\right\|_{\infty}=\left|e_{p, q}^1\right|=\operatorname{Max}_{\substack{1 \leq i \leq I-1 \\ 1 \leq j \leq I-1}}\,\,\left|e_{i, j}^1\right|$.

$\begin{aligned}\left|e_{i, j}^1\right| \leq\left|e_{p, q}^1\right| & =(1+2 r)\left|e_{p, q}^1\right|-\frac{r}{2}\left(4\left|e_{p, q}^1\right|\right). \\ \leq & (1+2 r)\,\,\left|e_{p, q}^1\,\,\right|-\frac{r}{2}\left(\left|e_{p+1, q}^1\,\,\right|+\left|e_{p-1, q}^1\,\,\right|+\left|e_{p, q+1}^1\,\,\right|+\left|e_{p, q-1}^1\,\,\right|\right). \\ & \leq\left|(1+2 r) e_{p, q}^1-\frac{r}{2}\left(e_{p+1, q}^1+e_{p-1, q}^1+e_{p, q+1}^1+e_{p, q-1}^1\right)\right|\end{aligned}$                       (17)

From (15) and (17), it follows that

$\left|e_{i, j}^1\right| \leq\left|\delta\left[f\left(x_p, y_q, t_{\frac{1}{2}}, \frac{3}{2} u\left(x_p, y_q, t_0\right)\right)-f\left(x_p, y_q, t_{\frac{1}{2}}, \frac{3}{2} u_{p, q}^0\right)\right]\right|+\left|T_{p, q}^{\frac{1}{2}}\right|$

Since satisfies Lipchitz condition (3), we obtain

$\left|e_{i, j}^1\right| \leq \delta \,\, L\left|u\left(x_p, y_q, t_n\right)-u_{p, q}^0\right|+\left|T_{p, q}^{\frac{1}{2}}\right|=\delta \,\, L\left|e_{p, q}^0\right|+\left|T_{p, q}^{\frac{1}{2}}\right|=\left|T_{p, q}^{\frac{1}{2}}\right| \leq C\left(k^{2-\alpha}+h^2\right)$

Thus $\left|e_{i, j}^1\right| \leq \frac{2 C}{\left(2^{1-\alpha} \,\,\,\, w_0\right)}\left(k^{2-\alpha}+h^2\right)$.

$\left\|E^1\right\| \leq\left(\frac{C_1}{2^{1-\alpha} \,\,\,\, w_0}\right)\left(k^{2-\alpha}+h^2\right), \quad C_1=2 C$

Now suppose that

$\left\|E^s\right\| \leq\left(\frac{C_1}{2^{1-\alpha} \,\,\,\, w_{s-1}}\right)\left(k^{2-\alpha}+h^2\right), C_s>0$, for $s=0,1,2, \ldots \ldots, n$

Set $M_s=\frac{C_s}{2^{1-\alpha} \,\,\,\, w_{s-1}}$, let $M=\operatorname{Max}\left\{C, C_1, C_2, \ldots \ldots \ldots, C_n\right\}$.

For $n+1$, Let $\left\|E^{n+1}\right\|_{\infty}=\operatorname{Max}_{\substack{1 \leq i \leq I-1 \\ 1 \leq j \leq I-1}}\left|e_{i, j}^{n+1}\right|=\left|e_{p, q}^{n+1}\right|$

$\begin{gathered}\left|e_{i, j}^{n+1}\right| \leq\left|e_{p, q}^{n+1}\right|=(1+2 r)\left|e_{p, q}^{n+1}\right|-\frac{r}{2}\left[4\left|e_{p, q}^{n+1}\right|\right] . \\ \leq(1+2 r)\left|e_{p, q}^{n+1}\right|-\frac{r}{2}\left[\left|e_{p+1, q}^{n+1}\right|+\left|e_{p-1, q}^{n+1}\right|+\left|e_{p, q+1}^{n+1}\right|+\left|e_{p, q-1}^{n+1}\right|\right] \\ \leq\left|(1+2 r) e_{p, q}^{n+1}-\frac{r}{2}\left(e_{p+1, q}^{n+1}+e_{p-1, q}^{n+1}+e_{p, q+1}^{n+1}+e_{p, q-1}^{n+1}\right)\right|\end{gathered}$                       (18)

From (16) and (18), it follows that

$\begin{gathered}\left|e_{i, j}^{n+1}\right| \leq\left|\left(1-2^{1-\alpha} w_1-2 r\right) e_{p, q}^n+\frac{r}{2}\left(e_{p+1, q}^n+e_{p-1, q}^n+e_{p, q+1}^n+e_{p, q-1}^n\right)\right|+2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right)\left|e_{i, j}^s\right|+ \\ \quad \delta\left|f\left(x_p, y_q, t_{n+\frac{1}{2}}, \frac{3}{2} u\left(x_p, y_q, t_n\right)-\frac{1}{2} u\left(x_p, y_q, t_{n-1}\right)\right)-f\left(x_p, y_q, t_{n+\frac{1}{2}}, \frac{3}{2} u_{p, q}^n-\frac{1}{2} u_{p, q}^{n-1}\right)\right|+\left|T_{p, q}^{n+\frac{1}{2}}\right|\end{gathered}$.

Since satisfies Lipchitz condition (3), we obtain

$\begin{gathered}\left|e_{i, j}^{n+1}\right| \leq\left(1-2^{1-\alpha} w_1-2 r\right)\left\|E^n\right\|_{\infty}+2 r\left\|E^n\right\|_{\infty}+2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right)\left\|E^s\right\|_{\infty}+\delta L\left[\frac{3}{2}\left|e_{p, q}^n\right|+\frac{1}{2}\left|e_{p, q}^{n-1}\right|\right]+\left|T_{p, q}^{n+\frac{1}{2}}\right|. \\ \leq\left(1-2^{1-\alpha} w_1\right) M_n\left(k^{2-\alpha}+h^2\right)+2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right) M_s\left(k^{2-\alpha}+h^2\right)+\delta L\left[\frac{3}{2} M_n\left(k^{2-\alpha}+h^2\right)_{\infty}+\right. \\ \left.\frac{1}{2} M_{n-1}\left(k^{2-\alpha}+h^2\right)\right]+C\left(k^{2-\alpha}+h^2\right) .\end{gathered}$

$\begin{gathered}\leq\left[1-2^{1-\alpha} w_1+2^{1-\alpha} \sum_{s=1}^{n-1}\left(w_{n-s}-w_{n-s+1}\right)+2 \,\,\delta L+2^{1-\alpha} w_n\right]\left(\frac{M}{2^{1-\alpha} w_n}\right)\left(k^{2-\alpha}+h^2\right). \\ =[1+2 \delta\,\, L]\left(\frac{1}{2^{1-\alpha} w_n}\right) M\left(k^{2-\alpha}+h^2\right) .\end{gathered}$

Thus $\left\|E^{n+1}\right\| \leq C_{n+1}\left(\frac{1}{2^{1-\alpha} \,\,\,\, w_n}\right)\left(k^{2-\alpha}+h^2\right)$, where $C_{n+1}=[1+2 \delta \,\,L] M$.