An Implicit Numerical Approach for 2D Rayleigh Stokes Problem for a Heated Generalized Second Grade Fluid with Fractional Derivative

Abstract

In this research work, our aim is to use the fast algorithm to solve the Rayleigh–Stokes problem for heated generalized second-grade fluid (RSP-HGSGF) involving Riemann–Liouville time fractional derivative. We suggest the modified implicit scheme formulated in the Riemann–Liouville integral sense and the scheme can be applied to the fractional RSP-HGSGF. Numerical experiments will be conducted, to show that the scheme is stress-free to implement, and the outcomes reveal the ideal execution of the suggested technique. The Fourier series will be used to examine the proposed scheme stability and convergence. The technique is stable, and the approximation solution converges to the exact result. To demonstrate the applicability and viability of the suggested strategy, a numerical demonstration will be provided.

1. Introduction

Fractional calculus is associated with the study and application of arbitrary order derivatives and integrals. The presentations of fractional-order calculus in many areas of science and engineering, including geometric phenomena, have aroused great interest in this field. The discussion between Leibniz and L’ Hospital at the end of the seventeenth century had the first discussion about fractional calculus. The Erdelyi, Abel, Riemann, Laplace, Heaviside, Levy, Liouville, Riesz, Grünwald, Letnikov, and Fourier are the of the great mathematician who worked on it and contributed. The use of fractional-order integrals and derivatives plays a significant part in the solution of some chemical problems, and this field has received more attention since 1968. Initially, many books on fractional calculus and its applications have been authored by scholars such as Ross and Miller [1], Spanier and Oldham [2], Podlubny [3] is the most well-known book in the field of fractional calculus, and Samko et al. [4] explains the underlying theory of fractional calculus, as well as its applications and solutions.

Many researchers have solved such types of fractional-order problems with different methods, such as Shivanian and Jafferabadi [5] solved the RSP-HGSGF with fractional derivatives to find the numerical solution through spectral meshless radial point interpolation. The Riemann–Liouville sense has been used to define the time-fractional derivative. The Shape functions are constructed using the point interpolation approach and radial basis functions, which act as basic functions. Zuhang and Liu [6] describe an efficient numerical approach for approximating RSP-HGSGF in a bounded domain. They carried out the investigation of stability and convergence for the suggested scheme. Wu [7] used the numerical approach to solve SFP-HGSGF. The INAS for the SFP-HGSGF has been explored in terms of its stability, convergence, and consistency. Shen et al. [8] examined RSP-HGSGF in a flow on a heated flat plate and within a heated edge. The fractional calculus technique in the constitutive relationship model was utilized to describe such a viscoelastic fluid. The Fourier transform on fractional-order Laplace operator is utilized for the exact solution for the velocity and temperature fields. Yu et al. [9] solved the RSP-HGSGF by the Adomian decomposition approach. In general, such types of series solutions converge quickly without discretizing the problem, the Adomian decomposition approach yields very precise numerical solutions. Chen et al. [10] presented two numerical approaches in this study for solving a two-dimensional variable-order sub-diffusion anomalous problem. Fourier analysis was used to examine their stability, convergence, and solvability. Yu et al. [11] considered the numerical approximation for Riemann–Liouville fractional-order derivative for the fractional SFP-HGSGF. They solved the direct and inverse problem by the implicit scheme with Riemann–Liouville fractional derivative. Lin and Jiang [12] devised a simple approach for a RSP-HGSG with fractional derivative. They used kernel theory to originate the series of the exact solution to the problem and established the approximate solution of its fractional derivative, which are uniformly convergent using truncating series. In the meantime, their technique comes with an error estimate and stability analysis. Chen et al. [13] founded the solution of RSP-HGSGF of fractional-order. They suggested implicit and explicit techniques. The convergence, stability, and solvability were established. Chen et al. [14] discussed Stokes’ initial challenge attention in recent years. The variable order nonlinear RSP-HGSGF is researched and further discussed the fourth-order numerical technique. The Fourier approach is used to examine the theoretical analysis of the numerical schemes as stability, solvability, and convergence. Dehghan and Abbaszadeh [15] worked on 2D fractional-order RSP-HGSGF on rectangular domains including circular, L-shaped, and a unit square with circular holes and created a numerical solution. The fractional derivatives are calculated using the RL principle. They integrated the equation for the time variable and obtained a full discrete scheme based on the Galerkin FEM for the space direction. Finally, we compare Galerkin’s FEM findings to those of alternative numerical techniques. Nikan et al. [16] used the radial basis function and fractional derivatives, and the Rayleigh–Stokes’s problem for an edge in a generalized Old Royd–B fluid is solved. The finite difference technique is used to discretize the temporal derivative terms. However, the local RBF-FD is used to discretize the spatial derivative terms. They evaluate the distribution of data nodes inside the local support area to maintain a constant number of nodes. The suggested method’s stability and convergence analysis are also explored. On irregular domains, the RBF-FD results are compared to those of previous approaches, demonstrating the viability and efficiency of the novel methodology. Shen et al. [17] studied the flow of RSP-HGSGF on a heated flat plate and within a heated edge. A fractional calculus methodology in the constitutive relationship model was employed to describe such a viscoelastic fluid. The Fourier transform and the fractional Laplace operator were used to obtain a closed-form solution for the velocity and temperature fields. Another survey [18,19] studied the same model to describe a viscoelastic fluid. Guan et al. [20] provided an enhanced version of a nonlinear source term with a fractional RSP for the finite difference/finite element technique. To discretize the first-order time derivative by the backward difference formula and second-order Grünwald–Letnikov, respectively. In the space direction, they employ the Galerkin finite element approach to define a fully discrete strategy for the fractional RSP-HGSGF with a nonlinear source term. The level of accuracy in the $L_2$ norm are calculated in detail using a novel analytical technique. Ali et al. [21] used the modified implicit difference approximation for 2D modified anomalous fractional sub-diffusion equation. The Fourier series approach is used to examine the convergence and stability of the proposed scheme. The scheme is shown to be unconditionally stable, and the approximate solution converges to the exact solution. Bazhlekova et al. [22] investigated the RSP-HGSGF in time with RL fractional derivative, and the analysis of the problem in space semi-discrete, continuous, and completely discrete formulations has been presented. Mohebbi et al. [23] presented a comparison of a fourth-order approximation, and the meshless approach is used for 2D fractional RSP and generated a completely discrete implicit scheme. Sun et al. [24] contributed a review article related to important information on fractional calculus. They discussed the main real-world application and powerful mathematical tools. Nikan et al. [25] studied the numerical solution of the nonlinear fractional-order reaction-subdiffusion model. They used the radial base function-finite-difference for spatial discretization and time discretization process are carried based on the weighted discrete scheme. They discussed the theoretical analysis and, for the computational efficiency of the suggested scheme, tested two numerical examples that reported accurate solutions. In another survey [26], the author formulated an efficient meshless scheme for the fractional-order diffusion model. They eliminated the time derivative by integrating both sides of the proposed model, and for space derivatives, used the local hybridization of cubic and radial basis functions. Nikan and Avazzadeh [27] studied the local hybrid kernel meshless approach for the approximation of the fractional-order model. They utilized the central difference approximation and Gaussian kernels to approximate the time and space direction, respectively. They checked the validity of the suggested method through numerical examples, which is accurate and efficient. In the research [28,29,30,31,32,33,34], different researchers discussed various numerical approaches for the time and space fractional-order models.

The aim of this study is to propose a modified implicit scheme for fractional-order RSP-HGSGF, which is new for this model. It reduces the computational cost and can easily investigate theoretical analysis by any method for the resultant scheme. In the procedure, substituting the Riemann–Liouville derivative by the discretized form of the Riemann–Liouville integral operator with the first-order time derivative. Then eliminate the partial derivative with respect to time by backward difference approximation. Additionally, the stability and convergence criterion of the established method was investigated by the Fourier series method. In the end, numerical examples are presented and solved by the proposed method to confirm the accuracy and feasibility of the method. The numerical examples are coded in Maple 15.

The rest of the paper is arranged as follows: in Section 2, discussed the methodology of the proposed scheme followed by stability and convergence analysis in Section 2.1 and Section 2.2. Section 3 presents the numerical experiments and results, and discussed are explained in Section 4. In the last Section 5 the conclusion is discussed.

Here, we consider the following two-dimensional RSP-HGSGF with fractional derivative [22].

$$\frac{\partial w\left(x,y,t\right)}{\partial t}=D_t^{1-\gamma }\left(\frac{{\partial }^{2}w\left(x,y,t\right)}{\partial x^2}+\frac{{\partial }^{2}w\left(x,y,t\right)}{\partial y^2}\right)+\frac{{\partial }^{2}w\left(x,y,t\right)}{\partial x^2}+\frac{{\partial }^{2}w\left(x,y,t\right)}{\partial y^2}+F\left(x,y,t\right),$$

(1)

Initial and boundary conditions are:

$$w\left(x,y,t\right)=\phi \left(x,y\right),$$

(2)

$$\begin{gathered} w\left(0,y,t\right)={\phi }_1\left(y,t\right), w(L,y,t)={\phi }_2\left(y,t\right), \\ w\left(x,0,t\right)={\phi }_3(y,t), w\left(x,L,t\right)={\phi }_4\left(x,t\right), \end{gathered}$$

$$0\le x,y\le L,0\le t\le T.$$

(3)

where ${ }_0{D}_t^{1-\gamma }w\left(x,y,t\right)$ represent the frcational-order Riemann–Liouville derivative of order $1-\gamma$ defined by

$${ }_0{D}_t^{1-\gamma }w\left(x,y,t\right)=\frac{1}{\Gamma (\gamma )} \frac{\partial }{\partial t} {{\displaystyle \int }}_0^t\frac{w\left(x,n\right)}{{\left(t-ƞ\right)}^{1-\gamma }}dƞ=\frac{\partial }{\partial t}{I}_0^{\gamma }w\left(x,y,t\right),$$

(4)

$$\mathrm{and} I_0^{\gamma }w\left(x,y,t\right)=\frac{1}{\Gamma (\gamma )}{{\displaystyle \int }}_0^t\frac{w\left(x,n\right)}{{\left(t-ƞ\right)}^{1-\gamma }} dƞ,$$

(5)

Represent the fractional-order $\left(\gamma \ge 0\right)$ Riemann–Liouville integral.

Lemma 1.

The γ $(0<\gamma <1)$order Riemann–Liouville fractional integral of the function$w\left(x, y,t\right)$on [0,T] can be defined in discretized form as,

$$I_0^{\gamma }u\left(x,y,t_k\right)=\frac{{\tau }^{\gamma }}{\Gamma \left(\gamma +1\right)} {{\displaystyle \Sigma }}_{j=0}^{k-1} b_j^{(\gamma )}w\left(x,y,t_{k-j}\right).$$

(6)

Lemma 2.

The coefficients constant$b_k^{(\gamma )}\left(k=0, 1, 2,\dots .\right)$fulfil the following properties [29]:

(i) 

$b_0^{\gamma }=1, b_k^{\gamma }>0 , k=0,1,2\dots ,$

(ii) 

$b_{k-1}^{\gamma }>b_k^{\gamma }, k=1,2\dots .,$

(iii) 

There exist a positive constant$C>0$. Such that$\tau \le C{b}_k^{\gamma }{\tau }^{\gamma }, k=1,2,\dots ,$

(iv) 

${{\displaystyle \Sigma }}_{j=0}^k{b}_j^{(\gamma )}{\tau }^{\gamma }={\left(k+1\right)}^{\gamma }\le T^{\gamma }.$

2. Methodology of the Proposed Scheme

The 2D RSP-HGSGF in Equations (1)–(3) solve by the modified implicit scheme. Utilized the Riemann–Liouville approximation for time-fractional and central difference for space derivative and partition the bounded domain into sub-intervals of length $\Delta x$ and $\Delta y$. The space steps as $x_i=i\Delta x$, in the x-direction with $i=1,\dots ,M_x-1, \Delta x=\frac{L}{M_x} and y_j=j\Delta y$, in the y-direction with $j=1,\dots , M_y-1, \Delta y=\frac{L}{M_y}$. The time step is $t_k=k\tau , k=1,\dots ,N$ where $\tau =T/N$. Let $u_{i,j}^k$ be the numerical approximation to $w\left(x_i,y_j,t_k\right),$ by applying (4) to (1), we obtain

$$\frac{\partial w\left(x,y,t\right)}{\partial t}=\frac{\partial }{\partial t}{I}_0^{\gamma }\left(\frac{{\partial }^{2}w\left(x,y,t\right)}{\partial x^2}+\frac{{\partial }^{2}w\left(x,y,t\right)}{\partial y^2}\right)+\frac{{\partial }^{2}w\left(x,y,t\right)}{\partial x^2}+\frac{{\partial }^{2}w\left(x,y,t\right)}{\partial y^2}+F\left(x,y,t\right).$$

(7)

Applying Lemma 1 and backward difference approximation w.r.t time, we obtain

$$w_{i,j}^k-w_{i,j}^{k-1}=S_1{{\displaystyle \Sigma }}_{j=o}^{k-1} b_j^{(\gamma )}\left(\delta x^2{w}_{I,j}^{k-j}-\delta x^2{w}_{i,j}^{k-j-i}\right)+S_{2 }{{\displaystyle \Sigma }}_{j=o}^{k-1} b_j^{(\gamma )}\left(\delta x^2{w}_{I,j}^{k-j}-\delta x^2{w}_{i,j}^{k-j-i}\right)+S_3\left(\frac{w_{i,j}^k}{\Delta x^2}\right)+S_4\left(\frac{w_{i,j}^k}{\Delta y^2}\right)+F\left(x, y,t\right).$$

(8)

where

$$S_1=\frac{{\tau }^{\gamma }}{\Gamma \left(\gamma +1\right)\Delta x^2}, S_2=\frac{{\tau }^{\gamma }}{\Gamma \left(\gamma +1\right)\Delta y^2}, S_3=\frac{\tau }{\Delta x^2}, S_4=\frac{\tau }{\Delta y^2}.$$

And

$$\delta x^2{w}_{i,j}^k=w_{i+1,j}^k-2w_{i,j}^k+w_{i-1,j}^k.$$

(9)

Simplified form of the proposed scheme for 2D RSP-HGSGF (1)–(3) and the conditions are as follows.

$$w_{i,j}^k-w_{i,j}^{k-1}=S_{1 }\delta x^2{w}_{I,j}^k-S_{1 }{b}_{k-1}^{(\gamma )}\delta x^2{w}_{I,j}^0-S_{1 }{{\displaystyle \Sigma }}_{s=1}^{k-1}(b_{s-\gamma }^{(\gamma )}-b_s^{(\gamma )})\delta x^2{w}_{I,j}^{k-s}+S_{2 }\delta y^2{w}_{I,j}^k+S_2{{\displaystyle \Sigma }}_{s=1}^{k-1}(b_{s-\gamma }^{(\gamma )}-b_s^{(\gamma )}) \delta y^2{w}_{I,j}^{k-s}+S_{3 }\delta x^2{w}_{I,j}^k+S_4\delta y^2{w}_{I,j}^k+\tau F_{i,j}^k$$

(10)

where $i=1, 2,\dots ,M_x-1,j=1, 2,\dots , M_y-1, and k=1, 2,\dots ,N-1$

$$w_{I,j}^0=\phi \left(x_i, y_j\right),$$

(11)

$$\begin{gathered} w_{0,j}^k={\phi }_1\left(y_j,t_k\right), w_{i,0}^k={\phi }_{2 }\left(x_i,t_k\right), \\ {w}_{M_x,j}^k={\phi }_3\left(y_j,t_k\right), w_{i,M_y}^k={\phi }_4\left(x_i,t_k\right), \end{gathered}$$

$$0\le x,y\le L, 0\le t\le T.$$

(12)

2.1. Stability

To find the stability of the proposed scheme by Fourier technique. Let the approximate solution is $W_{i,j}^{k }$ for (10), we have

$$W_{i,j}^{k }-W_{i,j}^{k-1}=S_1\left(W_{i+1,j}^{k }-2W_{i,j}^{k }+W_{i-1,j}^{k }\right)-S_1{b}_{k-1}^{(\gamma )}\left(W_{i+1,j}^0-2W_{i,j}^{0 }+W_{i-1,j}^{0 }\right)-S_{1 }{{\displaystyle \Sigma }}_{s=1}^{k-1}\left(b_{s-1}^{(\gamma )}-b_s^{(\gamma )}\right)(W_{i+1,j}^{k-s }-2W_{i,j}^{k-s}+W_{i,j+1}^{k-s })+S_2\left(W_{i,j}^{k-s}-2W_{i,j}^{k-s}+W_{i,j-1}^{k-s }\right)+S_2{b}_{k-1}^{(\gamma )}\left(W_{i,j+1}^0-2W_{i,j}^{0 }+W_{i,j-1}^{0 }\right)-S_2{{\displaystyle \Sigma }}_{s=1}^{k-1}(b_{s-1}^{(\gamma )}-b_s^{(\gamma )})\left(W_{i,j+1}^{k-s }-2W_{i,j}^{k-s}+W_{i,j-1}^{k-s }\right)+S_3\left(W_{i+1,j}^k-2W_{i,j}^{k }+W_{i-1,j}^{k }\right)+S_4\left(W_{i,j+1}^k-2W_{i,j}^k+W_{i,j-1}^{k }\right).$$

(13)

Next, the error is defined as

$${\Psi }_{i,j}^k=w_{i,j}^k-W_{i,j}^k.$$

(14)

where ${\Psi }_{i,j}^k$ satisfies (13) and

$${\Psi }_{i,j}^{k }-{\Psi }_{i,j}^{k-1}=S_1\left({\Psi }_{i+1,j}^k-2{\Psi }_{i,j}^{k }+{\Psi }_{i-1,j}^{k }\right)-S_{1 }{b}_{k-1}^{(\gamma )}\left({\Psi }_{i+1,j}^0-2{\Psi }_{i,j}^{0 }+{\Psi }_{i-1,j}^{0 }\right)-S_{1 }{{\displaystyle \Sigma }}_{s=1}^{k-1}\left(b_{s-1}^{(\gamma )}-b_s^{(\gamma )}\right)({\Psi }_{i+1,j}^{k-s }-2{\Psi }_{i,j}^{k-s}+{\Psi }_{i,j+1}^{k-s })+S_2\left({\Psi }_{i,j+1}^{k-s}-2{\Psi }_{i,j}^{k-s}+{\Psi }_{i,j-1}^{k-s }\right)+S_2{b}_{k-1}^{(\gamma )}({\Psi }_{i,j+1}^0-2{\Psi }_{i,j}^{0 }+{\Psi }_{i,j-1}^{0 })-S_2{{\displaystyle \Sigma }}_{s=1}^{k-1}(b_{s-1}^{(\gamma )}-b_s^{(\gamma )})\left({\Psi }_{i,j+1}^{k-s }-2{\Psi }_{i,j}^{k-s}+{\Psi }_{i,j-1}^{k-s }\right)+S_3\left({\Psi }_{i+1,j}^k-2{\Psi }_{i,j}^{k }+{\Psi }_{i-1,j}^{k }\right)+S_4({\Psi }_{i,j+1}^{k }-2{\Psi }_{i,j}^k+{\Psi }_{i,j-1}^{k }).$$

(15)

The error initial and boundary conditions are given as

$${\Psi }_{0,j}^k={\Psi }_{M_x,j}^k={\Psi }_{i,0}^k={\Psi }_{i,M_y}^k={\Psi }_{i,j}^0=0.$$

(16)

By defining the following grid functions for $k=1, 2\dots , N$

$${\Psi }^k\left(x,y\right)=\left\{\begin{array}{c}{\Psi }_{i,j}^k, when x_{i-\frac{\Delta x}{2}}<x\le x_{i+\frac{\Delta x}{2}}, y_{j-\frac{\Delta y}{2}}<y\le y_{i+\frac{\Delta y}{2},}\\ 0, when 0\le x\le \frac{\Delta x}{2} or L-\frac{\Delta x}{2}\le x\le L,\\ 0, when 0\le y\le \frac{\Delta y}{2} or L-\frac{\Delta y}{2}\le y\le L.\end{array}\right\}$$

(17)

Then ${\Psi }^k{}^{ }\left(x,y\right)$ can be expanded in Fourier series such as

$${\Psi }^k\left(x,y\right)={{\displaystyle \Sigma }}_{l_1,l_2=-\infty }^{\infty }{\lambda }^k\left(l_1,l_2\right)e^{\sqrt[2]{-1}\pi (l_{1}x/L+l_{2}y/L)}.$$

(18)

where,

$${\lambda }^k\left(l_1, l_2\right)=\frac{1}{L} {{\displaystyle \int }}_0^L{{\displaystyle \int }}_0^L{\Psi }^k\left(x,y\right)e^{\sqrt[-2]{-1}\pi (l_{1}x/L+l_{2}y/L)}dxdy.$$

(19)

From the definition of $l^2$ norm and Parseval equality, we have

$$\Vert {\Psi }^k\Vert {}_{\infty }^2={{\displaystyle \Sigma }}_{i=1}^{M_x-1}{{\displaystyle \Sigma }}_{j=1}^{M_y-1}\Delta x\Delta y{\left|{\Psi }_{i,j}^k\right|}^2={{\displaystyle \Sigma }}_{l_1,l_2=-\infty }^{\infty }{\left|{\lambda }^k\left(l_1,l_2\right)\right|}^2.$$

(20)

Supposing that

$${\Psi }_{i,j}^k={\lambda }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}.$$

(21)

where ${\sigma }_1=2\pi l_1/L, {\sigma }_2=2\pi l_2/L$ and substituting (21) in (15)

$${\lambda }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}-{\lambda }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}=S_1\left({\lambda }^k{e}^{\sqrt{-1}\left({\sigma }_1\left(i+1\right)\Delta x+{\sigma }_{2}j\Delta y\right)}-2{\lambda }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}+{\lambda }^k{e}^{\sqrt{-1}\left({\sigma }_1\left(i-1\right)\Delta x+{\sigma }_{2}j\Delta y\right)}\right)-S_{1 }{b}_{k-1}^{(\gamma )}\left({\lambda }^0{e}^{\sqrt{-1}\left({\sigma }_1\left(i+1\right)\Delta x+{\sigma }_{2}j\Delta y\right)}-2 {\lambda }^0{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}+ {\lambda }^0{e}^{\sqrt{-1}\left({\sigma }_1\left(i-1\right)\Delta x+{\sigma }_{2}j\Delta y\right)}\right)-S_{1 }{{\displaystyle \Sigma }}_{s=1}^{k-1}\left(b_{s-1}^{(\gamma )}-b_s^{(\gamma )}\right)({\lambda }^{k-s}{e}^{\sqrt{-1}\left({\sigma }_1\left(i+1\right)\Delta x+{\sigma }_{2}j\Delta y\right)}-2 {\lambda }^{k-s}{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}+{\lambda }^{k-s}{e}^{\sqrt{-1}\left({\sigma }_1\left(i-1\right)\Delta x+{\sigma }_{2}j\Delta y\right)})+S_2({\lambda }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_2\left(j+1\right)\Delta y\right)}-2{\lambda }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}+{\lambda }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_2\left(j-1\right)\Delta y\right)}+S_2{b}_{k-1}^{(\gamma )}\left({\lambda }^0{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}-2{\lambda }^0{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}+{\lambda }^0{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}\right)-S_2{{\displaystyle \Sigma }}_{s=1}^{k-1}(b_{s-1}^{(\gamma )}-b_s^{(\gamma )})\left({\lambda }^{k-s}{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_2\left(j+1\right)\Delta y\right)}-2{\lambda }^{k-s}{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}+{\lambda }^{k-s}{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_2\left(j-1\right)\Delta y\right)}\right)+S_3\left({\lambda }^k{e}^{\sqrt{-1}\left({\sigma }_1\left(i+1\right)\Delta x+{\sigma }_{2}j\Delta y\right)}-2{\lambda }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}+{\lambda }^k{e}^{\sqrt{-1}\left({\sigma }_1\left(i-1\right)\Delta x+{\sigma }_{2}j\Delta y\right)}\right)+S_4\left({\lambda }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_2\left(j+1\right)\Delta y\right)}-2{\lambda }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}+{\lambda }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_2\left(j-1\right)\Delta y\right)}\right),$$

After simplifying we obtain

$${\lambda }^k\left[1+{\mu }_1+{\mu }_2\right]={\lambda }^{k-1}+{\lambda }^0{b}_{k-1}^{(\gamma )}{\mu }_1+{\mu }_1{{\displaystyle \Sigma }}_{s=1}^{k-1}(b_{s-1}^{(\gamma )}-b_s^{(\gamma )}){\lambda }^{k-s},$$

and

$${\lambda }^k=\frac{ {\lambda }^{k-1}+{\lambda }^0{b}_{k-1}^{(\gamma )}{\mu }_1+{\mu }_1{{\displaystyle \Sigma }}_{s=1}^{k-1}(b_{s-1}^{(\gamma )}-b_s^{(\gamma )}){\lambda }^{k-s}}{\left[1+{\mu }_1+{\mu }_2\right]}.$$

(22)

where

$${\mu }_1=[4S_1{\sin }^2\frac{{\sigma }_1\Delta x}{2}+4S_2{\sin }^2\frac{{\sigma }_2\Delta y}{2}],$$

$${\mu }_2=[4S_3{\sin }^2\frac{{\sigma }_1\Delta x}{2}+4S_4{\sin }^2\frac{{\sigma }_2\Delta y}{2}].$$

Proposition 1.

If${\lambda }^k \left(k=1, 2,\dots ,N\right)$satisfy (22), then$\left|{\lambda }^{k+1}\right|\le \left|{\lambda }^0\right|$.

Proof. 

By using mathematical induction, we take $k=1$ in (22)

$${\lambda }^1=\frac{(1+b_0^{(\gamma )}{\mu }_1){\lambda }^0}{\left(1+{\mu }_1+{\mu }_2\right)},$$

and as ${\mu }_1,{\mu }_2\ge 0, b_0^{(\gamma )}=1, \mathrm{then}$

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

Now, assuming that

$$\left|{\lambda }^m\right|\le \left|{\lambda }^0\right|; m=1, 2,\dots ,k-1,$$

and as 0 $<\gamma <1$, from (22) and Lemma 2, we obtain

$$\left|{\lambda }^k\right|\le \frac{\left|{\lambda }^{k-1}\right|+b_{k-1}^{(\gamma )}{\mu }_1\left|{\lambda }^0\right|+{\mu }_1{{\displaystyle \Sigma }}_{s=1}^{k-1}\left(b_{s-1}^{(\gamma )}-b_s^{(\gamma )}\right){\lambda }^{k-s}}{1+{\mu }_1+{\mu }_2},$$

$$\le \frac{1+b_{k-1}^{(\gamma )}{\mu }_1+{\mu }_1{{\displaystyle \Sigma }}_{s=1}^{k-1}\left(b_{s-1}^{(\gamma )}-b_s^{(\gamma )}\right)}{\left(1+{\mu }_1+{\mu }_2\right)}\left|{\lambda }^0\right|,$$

$$=\frac{1+b_{k-1}^{(\gamma )}{\mu }_1+{\mu }_1\left(1-b_{k-1}^{(\gamma )}\right)}{\left(1+{\mu }_1+{\mu }_2\right)}\left|{\lambda }^0\right|,$$

$$=\frac{1+{\mu }_1}{1+{\mu }_1+{\mu }_2}{\lambda }^0,$$

$$\left|{\lambda }^k\right|\le \left|{\lambda }^0\right|.$$

(23)

This completes the proof. □

Based on the above proof. It can summarize that the solution of Equation (10) satisfies the following inequality

$${\Vert {\lambda }^k\Vert }_2\le {\Vert {\lambda }^0\Vert }_2.$$

And demonstrated that the proposed scheme is unconditionally stable.

2.2. Convergence

Here, we use a similar method to examine the convergence of the scheme. Let $u\left(x_i,y_j,t_k\right)$ represent the exact solution, then the truncation error of the scheme is obtained as: from (8)

$$T_{i,j}^k=u\left(x_i,y_j,t_k\right)-u\left(x_i,y_j,t_{k-1}\right)-S_1{{\displaystyle \Sigma }}_{j=o}^{k-1}{b}_s^{(\gamma )}\delta x^2\left(u\left(x_i,y_j,t_{k-s}\right)-u\left(x_i,y_j,t_{k-s-1}\right)\right)+S_{2 }{{\displaystyle \Sigma }}_{j=o}^{k-1}{b}_s^{(\gamma )}\delta y^2\left(u\left(x_i,y_j,t_{k-s}\right)-u\left(x_i,y_j,t_{k-s-1}\right)\right)+S_3\delta x^{2}u\left(x_i,y_j,t_k\right)+S_4\delta y^{2}u\left(x_i,y_j,t_k\right)-\tau f\left(x_i,y_j,t_k\right).$$

(24)

From (1) we have,

$$T_{i,j}^k=\frac{u_{i,j}^k-u_{i,j}^{k-1}}{\tau }-\frac{\partial u\left(x_i,y_j,t_k\right)}{\partial t}+\left(\frac{{\partial }^{2}u\left(x_i,y_j,t_k\right)}{\partial x^2}\right)-S_1{{\displaystyle \Sigma }}_{s=0}^{k-1}{b}_s^{(\gamma )}\delta x^2 \left(u_{i,j}^{k-s}-u_{i,j}^{k-s-1}\right)+ \left(\frac{{\partial }^{2}u\left(x_i,y_j,t_k\right)}{\partial y^2}\right)-S_2{{\displaystyle \Sigma }}_{s=0}^{k-1}{b}_s^{(\gamma )}\delta y^2\left(u_{i,j}^{k-s}-u_{i,j}^{k-s-1}\right)+\left(\frac{{\partial }^{2}u\left(x_i,y_j,t_k\right)}{\partial x^2}\right)-S_3\delta x^2\left(u_{i,j}^k\right)+\left(\frac{{\partial }^{2}u\left(x_i,y_j,t_k\right)}{\partial x^2}\right)-S_4\delta y^2\left(u_{i,j}^k\right)=O\left(\tau +\tau {\left(\Delta x\right)}^2\right)+\tau {\left(\Delta y\right)}^2).$$

(25)

Since $i,j and k$ are finite, thus there is a positive constant $C_1$, for all $i,j and k$, which then have,

$$\left|T_{i,j}^k\right|\le C_1\left(\tau +\tau {\left(\Delta x\right)}^2\right)+\tau {\left(\Delta y\right)}^2).$$

(26)

The error is defined as

$${\phi }_{i,j}^k=u\left(x_{i },y_j,t_k\right)-u_{i,j}^k.$$

(27)

From (24), we have

$$u\left(x_i, y_j,t_k\right)=u\left(x_i, y_j,t_{k-1}\right)+S_1\left(u\left(x_i, y_j,t_k\right)-2u\left(x_i, y_j,t_k\right)+u\left(x_i, y_j,t_k\right)\right)-S_1{b}_{k-1}^{(\gamma )}((u\left(x_{i+1}, y_j,t_0\right)-2u\left(x_i, y_j,t_0\right)+u\left(x_{i-1}, y_j,t_0\right)-S_1{{\displaystyle \Sigma }}_{s=1}^{k-1} \left(b_{s-1}^{(\gamma )}-b_s^{(\gamma )}\right)\left(u\left(x_{i+1}, y_j,t_{k-s}\right)-2u\left(x_i, y_j,t_{k-s}\right)-u\left(x_{i-1}, y_j,t_{k-s}\right)\right)+S_2(u\left(x_i, y_{j+1},t_k\right)-2u\left(x_i, y_j,t_k\right)+ u\left(x_i, y_{j-1},t_k\right)-S_2{b}_{k-1}^{(\gamma )} ((u\left(x_i, y_{j+1},t_k\right)-2u\left(x_i, y_j,t_k\right)+ u\left(x_i, y_{j-1},t_k\right)-S_2{{\displaystyle \Sigma }}_{s=1}^{k-1} \left(b_{s-1}^{(\gamma )}-b_s^{(\gamma )}\right)\left(u\left(x_i, y_{j+1},t_{k-s}\right)-2u\left(x_i, y_j,t_{k-s}\right)+u\left(x_i, y_{j-1},t_{k-s}\right)\right)+S_3 \left(u\left(x_{i+1}, y_{j+1},t_k\right)-2u\left(x_i, y_{j+1},t_k\right)+u\left(x_{i-1}, y_{j+1},t_k\right)\right)+S_4 \left(u\left(x_i, y_{j+1},t_k\right)-2u\left(x_i, y_{j+1},t_k\right)+u\left(x_i, y_{j-1},t_k\right)\right)+\tau f\left(x_i,y_j,t_k\right).$$

(28)

To obtain the error equation, subtract (28) from (10) to obtain.

$${\phi }_{i,j}^k-{\phi }_{i,j}^{k-1}=S_1\left({\phi }_{i+1,j}^k-2{\phi }_{i,j}^{k-1}+{\phi }_{i-1,j}^k\right)-S_1{b}_{k-1}^{(\gamma )} ({\phi }_{i+1,j}^k-2{\phi }_{i,j}^{k-1}+{\phi }_{i-1,j}^k-S_1{{\displaystyle \Sigma }}_{s=1}^{k-1}\left(b_{s-1}^{(\gamma )}-b_s^{(\gamma )}\right)\left({\phi }_{i+1,j}^k-2{\phi }_{i,j}^{k-1}+{\phi }_{i-1,j}^k\right)+S_2\left({\phi }_{i,j}^k-2{\phi }_{i,j}^{k-1}+{\phi }_{i,j}^k\right)-S_2{b}_{k-1}^{(\gamma )}\left({\phi }_{i,j+1}^0-2{\phi }_{i,j}^0+{\phi }_{i,j-1}^0\right)-S_2{{\displaystyle \Sigma }}_{s=1}^{k-1} \left(b_{s-1}^{(\gamma )}-b_s^{(\gamma )}\right)\left({\phi }_{i,j+1}^{k-s}2-{\phi }_{i,j}^{k-s}+{\phi }_{i,j-1}^{k-s}\right)+S_3\left({\phi }_{i+1,j}^k-2{\phi }_{i,j}^{k-1}+{\phi }_{i-1,j}^k\right)+S_4\left({\phi }_{i,j+1}^{k-s}-2{\phi }_{i,j}^{k-s}+{\phi }_{i,j-1}^{k-s}\right)+ \tau T_{i,j}^k.$$

(29)

With error boundary conditions

$${\phi }_{0,j}^k={\phi }_{M_x,j}^k={\phi }_{0,j}^k={\phi }_{i,M_y}^k=0, k=1,2,\dots ,N.$$

And the initial condition

$${\phi }_{i,j}^0=0, i=1,2,\dots ,M_x-1, j=1,2,\dots ,M_y-1.$$

Next, we define the following grid functions for $k=1,2,\dots ,N$

$${\phi }^k\left(x,y\right)=\left\{\begin{array}{c}{\phi }_{i,j}^k , when x_{i-\frac{\Delta x}{2}}<x\le x_{i+\frac{\Delta x}{2}}, y_{j-\frac{\Delta y}{2}}<y\le y_{j+\frac{\Delta y}{2}},\\ 0, when 0\le x\le \frac{\Delta x}{2} or L-\frac{\Delta x}{2}\le x\le L,\\ 0, when 0\le y\le \frac{\Delta y}{2} or L-\frac{\Delta y}{2}\le y\le L.\end{array}\right\}$$

And

$$T^k\left(x,y\right)=\left\{\begin{array}{c}{T}_{i,j}^k, when x_{i-\frac{\Delta x}{2}}<x\le x_{i+\frac{\Delta x}{2}}, y_{j-\frac{\Delta y}{2}}<y\le y_{j+\frac{\Delta y}{2}},\\ 0, when 0\le x\le \frac{\Delta x}{2} or L-\frac{\Delta x}{2}\le x\le L,\\ 0, when 0\le y\le \frac{\Delta y}{2} or L-\frac{\Delta y}{2}\le y\le L.\end{array}\right\}$$

Here, the ${\phi }^k\left(x,y\right) and T^k\left(x,y\right)$ can be expanded in fourier series such as

$${\phi }^k\left(x,y\right)=\underset{l_1,l_2=-\infty }{\overset{\infty }{{\displaystyle \Sigma }}}{\xi }^k\left(l_1,l_2\right)e^{2\sqrt{-1}\pi \left(l_{1}x/L+l_{2}y/L\right)}, k=1,2,\dots ,N.$$

$$T^k\left(x,y\right)=\underset{l_1,l_2=-\infty }{\overset{\infty }{{\displaystyle \Sigma }}}{\Psi }^k\left(l_1,l_2\right)e^{2\sqrt{-1}\pi \left(l_{1}x/L+l_{2}y/L\right)}, k=1,2,\dots ,N.$$

where

$${\xi }^k\left(l_1,l_2\right)=\frac{1}{L} {{\displaystyle \int }}_0^L{{\displaystyle \int }}_0^L{\phi }^k\left(x,y\right)e^{2\sqrt{-1}\pi \left(l_{1}x/L+l_{2}y/L\right)}dxdy,$$

(30)

$${\Psi }^k\left(l_1,l_2\right)=\frac{1}{L} {{\displaystyle \int }}_0^L{{\displaystyle \int }}_0^L{\phi }^k\left(x,y\right)e^{2\sqrt{-1}\pi \left(l_{1}x/L+l_{2}y/L\right)}dxdy.$$

(31)

From the definition of $l^{2 }$ norm and the Parseval equality, we have

$$\Vert {\phi }^k\Vert {}_{l^2}{}^2={{\displaystyle \Sigma }}_{i=1}^{M_x-1}{{\displaystyle \Sigma }}_{j=1}^{M_y-1}\Delta x\Delta y {\left|e_{i,j}^k\right|}^2={{\displaystyle \Sigma }}_{l_1,l_2=-\infty }^{\infty }{\left|{\rho }^k\left(l_1,l_2\right)\right|}^2,$$

(32)

$$\Vert T^k\Vert {}_{l^2}{}^2={{\displaystyle \Sigma }}_{i=1}^{M_x-1}{{\displaystyle \Sigma }}_{j=1}^{M_y-1}\Delta x\Delta y {\left|e_{i,j}^k\right|}^2={{\displaystyle \Sigma }}_{l_1,l_2=-\infty }^{\infty }{\left|{\Psi }^k\left(l_1,l_2\right)\right|}^2.$$

(33)

Based on the above supposing that

$${\phi }_i^k={\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}i\Delta y\right)},$$

(34)

$$T_i^k={\Psi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}i\Delta y\right)}.$$

(35)

Respectively, where ${\sigma }_1=\frac{2\pi l_1}{L}, {\sigma }_2=\frac{2\pi l_2}{L}$, substituting (34) and (35) into (29)

$$\begin{array}{cc}\hfill {\xi }^k{e}^{\sqrt{-1}({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta }& {}^{y)}-{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}\hfill \\ & = S_{1 }\left({\xi }^k{e}^{\sqrt{-1}\left({\sigma }_1\left(i+1\right)\Delta x+{\sigma }_{2}j\Delta y\right)}-2{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}+{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_1\left(i-1\right)\Delta x+{\sigma }_{2}j\Delta y\right)}\right)\hfill \\ & -S_{1 }{b}_{k-1}^{(\gamma )}\left({\xi }^k{e}^{\sqrt{-1}\left({\sigma }_1\left(i+1\right)\Delta x+{\sigma }_{2}j\Delta y\right)}-2{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}+{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_1\left(i-1\right)\Delta x+{\sigma }_{2}j\Delta y\right)}\right)\hfill \\ & -S_1\underset{s=1}{\overset{k-1}{{\displaystyle \Sigma }}}\left(b_{s-1}^{(\gamma )}-b_s^{(\gamma )}\right)({\xi }^k{e}^{\sqrt{-1}\left({\sigma }_1\left(i+1\right)\Delta x+{\sigma }_{2}j\Delta y\right)}-2{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}+{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_1\left(i-1\right)\Delta x+{\sigma }_{2}j\Delta y\right)})\hfill \\ & +S_2({\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_2\left(j+1\right)\Delta y\right)}-2{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}+{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_2\left(j-1\right)\Delta y\right)}\hfill \\ & +S_2{b}_{k-1}^{(\gamma )} \left({\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}-2{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}+{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}\right)-S_2\underset{s=1}{\overset{k-1}{{\displaystyle \Sigma }}}(b_{s-1}^{(\gamma )}\hfill \\ & -b_s^{(\gamma )}) ({\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_2\left(j+1\right)\Delta y\right)}-2{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}+{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_2\left(j-1\right)\Delta y\right)})\hfill \\ & +S_3({\xi }^k{e}^{\sqrt{-1}\left({\sigma }_1\left(i+1\right)\Delta x+{\sigma }_{2}j\Delta y\right)}-2{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}+{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_1\left(i-1\right)\Delta x+{\sigma }_{2}j\Delta y\right)})\hfill \\ & +S_4({\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_2\left(j+1\right)\Delta y\right)}-2{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}+{\xi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_2\left(j-1\right)\Delta y\right)})\hfill \\ & +\tau \left({\Psi }^k{e}^{\sqrt{-1}\left({\sigma }_{1}i\Delta x+{\sigma }_{2}j\Delta y\right)}\right), \hfill \end{array}$$

Simplifying the above equation, we obtained

$${\xi }^k=\frac{ {\xi }^{k-1}+{\xi }^0{b}_{k-1}^{(\gamma )}{\mu }_1+{\mu }_1{{\displaystyle \Sigma }}_{s=1}^{k-1}(b_{s-1}^{(\gamma )}-b_s^{(\gamma )}){\xi }^{k-s}}{\left[1+{\mu }_1+{\mu }_2\right]},$$

(36)

where

$${\mu }_1=4S_1{\sin }^2\frac{{\sigma }_1\Delta x}{2}+4S_2{\sin }^2\frac{{\sigma }_2\Delta y}{2}],$$

$${\mu }_2=4S_3{\sin }^2\frac{{\sigma }_1\Delta x}{2}+4S_4{\sin }^2\frac{{\sigma }_2\Delta y}{2}].$$

Proposition 2.

Let ${\xi }^k\left(k=1,2,\dots ,N\right)$be the solution of (36), then there is a positive constant$C_2$thus that,

$$\left|{\xi }^k\right|\le C_{2}k\tau \left|{\Psi }^1\right|.$$

Proof. 

From ${\phi }^0=0$ and (30), we have

$${\xi }^o={\xi }^o\left(l_1,l_2\right)=0.$$

(37)

From (31) and (33), then there is a positive constant $C_2$, such that

$$\left|{\Psi }^k\right|\le C_2\left|{\Psi }^1\left(l_1,l_2\right)\right|.$$

(38)

Using mathematical induction, for $k=1$, then from (36) and (37), we obtain

$${\xi }^1=\frac{1}{1+{\mu }_1+{\mu }_2}\left(\tau {\Psi }^1\right).$$

Since ${\mu }_1,{\mu }_2\ge 0,$ from (3.38), we obtain

$$\left|{\xi }^1\right|\le \tau \left|{\Psi }^1\right|\le C_2\tau \left|{\Psi }^1\right|.$$

Now suppose that

$$\left|{\xi }^m\right|\le C_{2}m\tau \left|{\Psi }^1\right|, m=1,2,\dots ,k-1.$$

As 0 $<\gamma <1, {\mu }_1,{\mu }_2\ge 0$.

From (35), (37) and Lemma 2 we have

$$\left|{\xi }^k\right|=\frac{{|\xi |}^{k-1}+{\mu }_1{{\displaystyle \Sigma }}_{s=1}^{k-1}\left(b_{s-1}^{\gamma }-b_s^{\gamma }\right){|\xi |}^{k-s}+\tau {|\Psi |}^k}{\left(1+{\mu }_1+{\mu }_2\right)},$$

$$\left|{\xi }^k\right|=\frac{C_2\left(k-1\right)\tau {|\Psi |}^1+{\mu }_1{{\displaystyle \Sigma }}_{s=1}^{k-1}\left(b_{s-1}^{\gamma }-b_s^{\gamma }\right)C_2\left(k-s\right)\tau {|\Psi |}^1+C_2\tau {|\Psi |}^1}{\left(1+{\mu }_1+{\mu }_2\right)},$$

$$\le \left[\frac{\left(k-1\right)+{\mu }_1\left(k-1\right){{\displaystyle \Sigma }}_{s=1}^{k-1}\left(b_{s-1}^{\gamma }-b_s^{\gamma }\right)+1}{\left(1+{\mu }_1+{\mu }_2\right)}\right]C_2\tau {|\Psi |}^1,$$

$$=\left[\frac{k+{\mu }_1\left(k-1\right){{\displaystyle \Sigma }}_{s=1}^{k-1}\left(b_{s-1}^{\gamma }-b_s^{\gamma }\right)+1}{\left(1+{\mu }_1+{\mu }_2\right)}\right]C_2\tau {|\Psi |}^1,$$

$$=\left[\frac{k+{\mu }_1\left(k-1\right)+\left(1-b_{k-1}^{(\gamma )}\right)}{\left(1+{\mu }_1+{\mu }_2\right)}\right]C_2\tau {|\Psi |}^1,$$

$$\le k{C}_2\tau \left|{\Psi }^1\right|.$$

(39)

Proof completed via induction method. □

Theorem 1.

The modified implicit difference scheme$l^2$is convergent and the order of convergence is $O\left(\tau +\tau {\left(\Delta x\right)}^2+\tau {\left(\Delta y\right)}^2\right).$

Proof. 

From (25) and (33), we obtain

$$\Vert T^k\Vert \le \sqrt{M_x\Delta x}\sqrt{M_y\Delta y} C_1\left(\tau +\tau {\left(\Delta x\right)}^2+\tau {\left(\Delta y\right)}^2\right)=L{C}_1\left(\tau +\tau {\left(\Delta x\right)}^2+\tau {\left(\Delta y\right)}^2\right).$$

$$\Vert {\phi }^k\Vert {}_{l^2}\le k{C}_2\tau \Vert T^1\Vert \le C_1{C}_{2}k\tau L\left(\tau +\tau \left(\Delta x^2\right)+\tau \left(\Delta y^2\right)\right),$$

(40)

As $k\tau \le T,$ thus

$$\Vert {\phi }^k\Vert {}_{l^2}\le C_1{C}_{2}k\tau L\left(\tau +\tau \left(\Delta x^2\right)+\tau \left(\Delta y^2\right)\right),$$

where C =$C_1{C}_{2}TL.$ □

3. Numerical Experiment

Example 1.

Consider the following two-dimensional Rayleigh stokes problem for heated generalized second grade fluid with fractional derivative [22],

$$\frac{\partial u\left(x,y,t\right)}{\partial t}={ }_0{D}_t^{1-\gamma }\left(\frac{{\partial }^{2}u\left(x,y,t\right)}{\partial x^2}+\frac{{\partial }^{2}u(x,y,t}{\partial y^2}\right)+\frac{{\partial }^{2}u\left(x,y,t\right)}{\partial x^2}+\frac{{\partial }^{2}u\left(x,y,t\right)}{\partial y^2}+f\left(x,y,t\right),$$

$$0\le \gamma \le 1, 0\le t\le T,$$

(41)

with initial and boundary conditions,

$$u\left(x,y,0\right)=0, 0\le x, y\le 1,$$

(42)

$$u\left(0,y,t\right)=e^y{t}^{1+\gamma }, u\left(1,y,t\right)=e^{1+y}{t}^{1+\gamma },$$

$$u\left(x,0,t\right))=e^x{t}^{1+r}, u\left(x,1,t\right)=e^{1+x}{t}^{1+\gamma }, 0\le t\le T.$$

(43)

Here $f\left(x,y,t\right)=\left(\left(1+\gamma \right)t^{\gamma }-2\frac{\Gamma \left(2+\gamma \right)}{\Gamma \left(1+2\gamma \right)}{t}^{2\gamma }-2t^{1+\gamma }\right)e^{x+y}$ and the exact solution of (41) is given by

$$u\left(x,y,t\right)=e^{x+y}{t}^{1+\gamma }.$$

(44)

The error between the numerical solution and exact solution is defined as follows.

$$E_{\infty }=\underset{0\le i,j\le M,0\le k\le N}{\max }\left|u\left(x_i,y_j,t_k\right)-u_{i,j}^k\right|.$$

(45)

And rate of convergence for space variable can be define as

$$C_2-\mathrm{order}={\log }_2\left(\frac{\parallel E_{\infty }\left(16\tau , 2\Delta x\right)\parallel }{\parallel E_{\infty }\left(\tau , \Delta x\right)\parallel }\right).$$

(46)

The developed modified implicit scheme is applied to problems (41) to (44).

All following

Table 1. The error table for different values at $\tau ,\Delta x,\Delta y, \mathrm{and} \gamma$.

$\mathit{\tau }$	$\mathbf{\Delta }\mathit{x}=\mathbf{\Delta }\mathit{y}$	$\mathit{\gamma }=0.5$	$\mathit{\gamma }=0.6$	$\mathit{\gamma }=0.7$	$\mathit{\gamma }=0.8$	$\mathit{\gamma }=0.9$
$1⁄4$	1/2	0.013169	0.016017	0.018890	0.021807	0.024819
$1/8$	$1/8$	$0.001404$	$0.001638$	$0.008559$	$0.010150$	$0.001857$
$1/16$	$1/4$	$0.004573$	$0.005267$	$0.005934$	$0.006606$	$0.007312$
$1/64$	$1/8$	$0.001330$	$0.001484$	$0.001634$	$0.001789$	$0.001958$
$1/128$	$1/10$	$0.000770$	$0.000842$	$0.000914$	$0.000989$	$0.001073$
$1/144$	$1/12$	$0.000624$	$0.000688$	$0.000752$	$0.000819$	$0.000894$

,

Table 2. The error table for different values at $\tau ,\Delta x,\Delta y, \mathrm{and} \gamma$.

$\mathit{\tau }=\mathbf{\Delta }\mathit{x}=\mathbf{\Delta }\mathit{y}$	$\mathit{\gamma }=0.35$	$\mathit{\gamma }=0.65$	$\mathit{\gamma }=0.85$
1/2	$0.007989$	$0.024671$	$0.038181$
1/4	$0.004538$	$0.014258$	$0.021050$
$1/6$	$0.003387$	$0.009935$	$0.014272$
$1/8$	$0.002814$	$0.007770$	$0.010965$
$1/10$	$0.002428$	$0.006397$	$0.008915$

,

Table 3. The error table for different values at $\tau ,\Delta x,\Delta y$ and at fixed value of $\gamma =0.25$.

$\mathit{N}$	$\mathbf{\Delta }\mathit{x}=\mathbf{\Delta }\mathit{y}=1/5$	$\mathbf{\Delta }\mathit{x}=\mathbf{\Delta }\mathit{y}=1/10$	$\mathbf{\Delta }\mathit{x}=\mathbf{\Delta }\mathit{y}=1/15$	$\mathbf{\Delta }\mathit{x}=\mathbf{\Delta }\mathit{y}=1/20$
20	$0.0019274$	$0.0009539$	$0.0007670$	$0.0007014$
40	$0.0017467$	$0.0016060$	$0.0005798$	$0.0004302$
60	$0.0016504$	$0.0006685$	$0.0005453$	$0.0003295$

and

Table 4. The error $E_{\infty }$ and convergence order for space variable for Equation (41) for different values at $\tau ,\Delta x,\Delta y$ and at fixed value of $\gamma =0.5$.

	${\mathit{E}}_{\mathit{\infty }}$	${\mathit{C}}_2-\mathit{O}\mathit{r}\mathit{d}\mathit{e}\mathit{r}$	CPU-Time (s)
$\Delta x=\Delta y=\tau =\frac{1}{4}$	$0.0092878$	$-$	$0.04$
$\Delta x=\Delta y=\frac{1}{8}, \tau =\frac{1}{64}$	$0.0017302$	$2.425$	$499.68$
$\Delta x=\Delta y=\tau =\frac{1}{8}$	$0.0053546$	$-$	$2.03$
$\Delta x=\Delta y=\frac{1}{16},\tau =\frac{1}{128}$	$0.0009824$	$2.446$	$1095.67$

shows the errors $E_{\infty }$ for values of space step size ($\Delta x,\Delta y$) and $\tau$. where time step, $\tau$ is defined by $\tau =\frac{T}{N}$.

The above tables indicate that, as we reduce the time and space step size $\tau$ and ($\Delta x,\Delta y$), the error decrease for a fixed value of $\gamma$. This shows that the method converges to the exact solution.

4. Results and Discussion

A modified implicit scheme was developed and applied on RSP-HGSGF of fractional-order. A numerical example was given to support the theoretical study of stability and convergence analysis. The error between the exact and numerical solution was calculated using different values of $N$ and $M$. In addition, for different values of $\gamma$, Table 1, Table 2 and Table 3 shows the comparison of the numerical scheme with the exact solution in terms of maximum error. The error was calculated using Maple 15 software reduced as increase the number of space and time steps for different values of fractional-order $\gamma$. Table 4 represents the order of convergence for space variable performed in excellent agreement with theoretical analysis. Figure 1, Figure 2 and Figure 3 are plotted for different values of $M, N$ and fractional-order $\gamma$, such that Figure 1 is at $M=4, N=2, \gamma =0.25$, Figure 2 is at $M=N=6, \gamma =0.35$, and Figure 3 at $M=N=8, \gamma =0.85$ shows good agreement with the exact solution.

5. Conclusions

A modified implicit difference scheme is formulated for 2D RSP-HGSGF with derivative of fractional-order has been described in this research. The modified scheme based on Riemann–Liouville fractional integral operator, has been used because of the improvement of low computational cost and can easily be applied. The Fourier technique has been used for the theoretical analysis of stability and convergence with order $\left(\tau +{\left(\Delta x\right)}^2+{\left(\Delta y\right)}^2\right)$ is unconditionally stable and convergent. The numerical experiment is conducted, which shows that the modified implicit scheme is easy to implement, and the results show good performance of the proposed schemes. The results in Table 4 also support the effectiveness of the examined theoretical analysis. In the future we will extend this technique to the space fractional-order differential equations.