# The Rayleigh–Stokes Problem for a Heated Generalized Second-Grade Fluid with Fractional Derivative: An Implicit Scheme via Riemann–Liouville Integral

Authors: Abdul Hamid Ganie; Abdulkafi Mohammed Saeed; Sadia Saeed; Umair Ali
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

The goal of this study is to use the fast algorithm to solve the Rayleigh–Stokes problem for heated generalized second-grade fluid (RSP-HGSGF) with Riemann–Liouville time fractional derivative using the fast algorithm. The modified implicit scheme, which is formulated by the Riemann–Liouville integral formula and applied to the fractional RSP-HGSGF, is proposed. Numerical experiments will be carried out to demonstrate that the scheme is simple to implement, and the results will reveal the best way to implement the suggested technique. The proposed scheme’s stability and convergence will be examined using the Fourier series. The method is stable, and the approximation solution approaches the exact solution. A numerical demonstration will be provided to demonstrate the applicability and viability of the suggested strategy.


---

## Body

## 1. Introduction

The study and application of arbitrary-order derivatives and integrals are associated with fractional calculus. The use of fractional-order calculus in a variety of fields of science and engineering, including geometric phenomena, has sparked a lot of interest in this area [1]. The first discussion of fractional calculus took place between Leibniz and L’Hospital at the end of the seventeenth century [2]. The great mathematicians Erdelyi, Abel, Riemann, Laplace, Heaviside, Levy, Liouville, Riesz, Gunwald, Letnikov, and Fourier worked on it and had contributed [3]. Fractional-order integrals and derivatives play an important role in solving some chemical problems, and this field has been paid much attention since 1968. The most well-known book in the field of fractional calculus, originally written by Ross and Miller and Ross [4], Spanier and Oldham [5], Podlubny [6], and Samko et al. [7], explains the underlying theory of fractional calculus as well as its applications and solutions.Many researchers have solved fractional-order problems using various methods. For example, Shivanian and Jafferabadi [8] used fractional derivatives to find the numerical solution for the RSP-HGSGF using spectral meshless radial point interpolation. The time-fractional derivative has been defined in the Riemann–Liouville sense. The Shape functions are created by using a point interpolation method and radial basis functions as basic functions. An efficient numerical approach for approximating RSP-HGSGF in a bounded domain is described by Liu et al. [9]. They investigated the proposed scheme’s stability and convergence. To solve SFP-HGSGF, Wu [10] used a numerical approach. The stability, convergence, and consistency of the INAS for the SFP HGSGF have been investigated. RSPHGSGF was studied in a flow on a heated flat plate and within a heated edge by Shen et al. [11]. A viscoelastic fluid was described using the fractional calculus technique in the constitutive relationship model. For the exact solution of the velocity and temperature fields, the Fourier transform on fractional-order Laplace operator is used. Yu et al. [12] used the Adomian decomposition method to solve the RSP-HGSGF. In general, without discretizing the problem, such series solutions converge quickly and the Adomian decomposition approach yields very precise numerical solutions. In this study, Chen et al. [13] presented two numerical methods for solving a two-dimensional variable-order subdiffusion anomalous problem. Their stability, convergence, and solvability were investigated using Fourier analysis. The numerical approximation for the Riemann–Liouville fractional-order derivative for the fractional SFP-HGSGF was studied by Yu et al. [14]. They used the implicit scheme with Riemann–Liouville fractional derivative to solve the direct and inverse problems. Lin and Jiang [15] devised a straightforward method for calculating the fractional derivative of an RSP-HGSG. They created the series of the exact solution to the problem using kernel theory and established the approximate solution of its fractional derivative using truncating series, which are uniformly convergent. Meanwhile, their method includes error estimation and stability analysis. Chen et al. [16] proposed the implicit and explicit techniques for solving the RSP-HGSGF of fractional order. The convergence, stability, and solvability of the problem have all been determined. In recent years, Chen et al. [17] discussed Stokes’ initial challenge attention. The variable-order nonlinear RSP-HGSGF is investigated, and the fourth-order numerical technique is discussed. The Fourier approach is used to investigate the numerical scheme's theoretical analysis. Dehghan and Abbas Zadeh [18] developed a numerical solution for 2D fractional-order RSP-HGSGF on rectangular domains such as circular, L-shaped, and a unit square with circular holes. The RL principle is used to calculate the fractional derivatives. They used the Galerkin FEM to obtain a fully discrete scheme for the space direction by integrating the equation for the time variable. Finally, we compare the results of Galerkin FEM to those of other numerical techniques. The Rayleigh–Stokes problem for an edge in a generalized Oldroyd-B fluid was solved by Nikan and Avazzadeh et al. [19] using the radial basis function and fractional derivatives. The temporal derivative terms are discretized using the finite difference technique, while the spatial derivative terms are discretized using the local RBF-FD.To maintain a constant number of nodes, they evaluate the distribution of data nodes within the local support area. The stability and convergence of the proposed method are also investigated. The RBF-FD results are compared to those of previous approaches on irregular domains, demonstrating the novel methodology’s viability and efficiency. RSP-HGSGF flow was investigated by Zhai et al. [20] on a heated flat plate and within a heated edge. To describe such a viscoelastic fluid, a fractional calculus methodology was used in the constitutive relationship model. The velocity and temperature fields were solved in closed form using the Fourier transform and the fractional Laplace operator. Another study looked at the same model to describe a viscoelastic fluid [21, 22]. For the finite difference/finite element technique, Guan et al. [23] provided an enhanced version of a nonlinear source term with a fractional RSP. The backward difference formula and second-order Grünwald–Letnikov derivative are used to discretize the first-order time derivative. They use the Galerkin finite element approach to define a fully discrete strategy for the fractional RSP-HGSGF with a nonlinear source term in the space direction. A novel analytical technique is used to calculate the level of accuracy in the L2 norm in great detail. For the 2D modified anomalous fractional subdiffusion equation, Ali et al. [24] used a modified implicit difference approximation. The proposed scheme’s convergence and stability are investigated using the Fourier series approach. It is shown that the scheme is unconditionally stable, and that the approximate solution converges to the exact solution. Bazhlekova et al. [25] investigated the RSP-HGSGF in time using the RL fractional derivative, and the problem was analysed in space using semidiscrete, continuous, and completely discrete formulations. Mohebbi et al. [26] compared the meshless approach to a fourth-order approximation for 2D fractional RSP and generated a completely discrete implicit scheme. Sun et al [27] contributed a review article on important fractional calculus information. They talked about the most important real-world applications as well as powerful mathematical tools. The numerical solution of a nonlinear fractional-order reaction-subdiffusion model was investigated by Nikan et al. [28]. For spatial discretization, they used the radial base function-finite difference method, and for time discretization, they used a weighted discrete scheme. They discussed theoretical analysis and tested two numerical examples for the computational efficiency of the proposed scheme, which yielded accurate results. In a separate study [29], the author proposed a meshless scheme for the fractional-order diffusion model. They eliminated the time derivative by integrating both sides of the proposed model and used local hybridization of cubic and radial basis functions for space derivatives. Nikan et al. [30] investigated the local hybrid kernel meshless approach for fractional-order model approximation. To approximate the time and space directions, they used the central difference approximation and Gaussian kernels, respectively. They verified the validity of the proposed method using numerical examples that are both accurate and efficient. Liu et al. [31] discussed the fractional dynamics modelled from the fractional-order PDEs. Fractional-order systems have importance in the field of electrochemistry, chaotic systems, biology etc. Ahmad et al. [32] formulated a new methodology named as variational iteration method I and successfully applied to a nonlinear model. They explained the compactness of the method and compared their results with the existed literature and found that the proposed method is more productive and reliable than others. Khan et al. [33] considered the numerical approach based on the collocation method for the inverse heat source problem and tasted the method both on regular and irregular domain. Different researchers discussed various numerical approaches for time and space fractional-order models in the research [27, 30, 34–38]. The goal of this research is to propose a new scheme for this model modified implicit scheme for fractional-order RSP-HGSGF. It lowers the computational cost and allows for easy theoretical analysis using any method for the final scheme. In the procedure, the discretized form of the Riemann–Liouville integral operator is used to replace the Riemann–Liouville derivative with the first-order time derivative. The partial derivative with respect to time is then eliminated using backward difference approximation. Additionally, we use the Fourier series method to investigate the established method's stability and convergence criterion. Finally, numerical examples are presented and solved using the proposed method to verify the method’s accuracy and feasibility. Maple 15 is used to code the numerical examples.The following is how the rest of the paper is organized: The methodology of the proposed scheme is discussed in Section2, followed by stability and convergence analysis in Sections 2.1 and 2.2. The numerical experiments and results are presented in Section 3 and discussed in Section 4. The conclusion is discussed in Section 5 of the report.The aim of this study is to propose a modified implicit scheme for fractional RSP-HGSGF based on the formulated Riemann–Liouville integral operator. The partial derivative w.r.t. time is eliminated by backward difference approximation. Additionally, we investigate the stability and convergence criterion of the established method by the Fourier series method.Here, we consider the following two-dimensional RSP-HGSGF with fractional derivative [22].(1)∂Yx,y,t∂t=Dt1−β∂2Yx,y,t∂x2+∂2Yx,y,t∂y2+∂2Yx,y,t∂x2+∂2Yx,y,t∂y2+Hx,y,t.Initial and boundary conditions are as follows:(2)Yx,y,t=φx,y,(3)Y0,y,t=Ω1y,t,YL,y,t=Ω2y,t,Yx,0,t=Ω3y,tYx,L,t=Ω4x,t,0≤x,y≤L,0≤t≤T,where 0Dt1−βYx,y,t represents the fractional-order Riemann–Liouville derivative of order 1−β.Lemma 1.
Theβ0〈β〈1-order Riemann–Liouville fractional integral of the function Yx,y,t on [0,T] can be defined in discretized form as(4)I0βYx,y,tm=τβГβ+1∑j=0m−1djβYx,y,tm−j.Lemma 2.
The coefficients constantdmβm=0,1,2,…. fulfils the following properties [29]:(i)
d0β=1,dmβ〉0,m=0,1,2…(ii)
dm−1β〉dmβ,m=1,2…(iii)
There exists a positive constantC〉0, such that τ≤Cdmβτβ,m=1,2,…(iv)
∑j=0mdjβτβ=m+1β≤Tβ


## 2. Methodology of the Proposed Scheme

The 2D RSP-HGSGF in equations (1)–(3) is solved by the modified implicit scheme. We utilized the Riemann–Liouville approximation for time-fractional and central difference for space derivative and partitioned the bounded domain into subintervals of lengths Δx and Δy. The space steps are xi=iΔx, in the x-direction with i=1,…,M1−1,Δx=L/M1,andyj=jΔy, in the y-direction with j=1,…,M2−1,Δy=L/M2. The time step is tm=mτ,m=1,…,N where τ=T/N. Let Yi,jm be the numerical approximation to Yxi,yj,tm; by applying (2) to (1), we obtain(5)∂Yx,y,t∂t=∂∂tI0β∂2Yx,y,t∂x2+∂2Yx,y,t∂y2+∂2Yx,y,t∂x2+∂2Yx,y,t∂y2+Hx,y,t.Applying Lemma1 and backward difference approximation w.r.t. time, we obtain(6)Yi,jm−Yi,jm−1=R1∑j=om−1djβδx2YI,jm−j−δx2Yi,jm−j−i+R2∑j=om−1djβδx2YI,jm−j−δx2Yi,jm−j−i+R3Yi,jmΔx2+R4Yi,jmΔy2+Hx,y,t,where(7)R1=τβГβ+1Δx2,R2=τβГβ+1Δy2,R3=τΔx2,R4=τΔy2.(8)δx2Yi,jm=Yi+1,jm−2Yi,jm+Yi−1,jm.The simplified form of the proposed scheme for 2D RSP-HGSGF (1)-(3) and the conditions are as follows:(9)Yi,jm−Yi,jm−1=R1δx2YI,jm−R1dm−1βδx2YI,j0−R1∑s=1m−1ds−ββ−dsβδx2YI,jm−s+R2δy2YI,jm+R2∑s=1m−1ds−ββ−dsβδy2YI,jm−s+R3δx2YI,jm+R4δy2YI,jm+τHi,jm,where i=1,2,…,M1−1,j=1,2,…,M2−1,andm=1,2,…,N−1.(10)YI,j0=Ωxi,yj,Y0,jm=Ω1yj,tm,Yi,0m=Ω2xi,tm,YM1,jm=Ω3yj,tm,Yi,M2m=Ω4xi,tm,0≤x,y≤L,0≤t≤T.


### 2.1. Stability

We find the stability of the proposed scheme by Fourier technique. Let the approximate solution beΨi,jm for (9); we have(11)Ψi,jm−Ψi,jm−1=R1Ψi+1,jm−2Ψi,jm+Ψi−1,jm−R1dm−1βΨi+1,j0−2Ψi,j0+Ψi−1,j0−R1∑s=1m−1ds−1β−dsβΨi+1,jm−s−2Ψi,jm−s+Ψi,j+1m−s+R2Ψi,jm−s−2Ψi,jm−s+Ψi,j−1m−s+R2dm−1βΨi,j+10−2Ψi,j0+Ψi,j−10−R2∑s=1m−1ds−1β−dsβΨi,j+1m−s−2Ψi,jm−s+Ψi,j−1m−s+R3Ψi+1,jm−2Ψi,jm+Ψi−1,jm+R4Ψi,j+1m−2Ψi,jm+Ψi,j−1m.Next, the error is defined as(12)φi,jm=Yi,jm−Ψi,jm.where φi,jm satisfies (11) and(13)φi,jm−φi,jm−1=R1φi+1,jm−2φi,jm+φi−1,jm−R1dm−1βφi+1,j0−2φi,j0+φi−1,j0−R1∑s=1m−1ds−1β−dsβφi+1,jm−s−2φi,jm−s+φi,j+1m−s+R2φi,j+1m−s−2φi,jm−s+φi,j−1m−s+R2dm−1βφi,j+10−2φi,j0+φi,j−10−R2∑s=1m−1ds−1β−dsβφi,j+1m−s−2φi,jm−s+φi,j−1m−s+R3φi+1,jm−2φi,jm+φi−1,jm+R4φi,j+1m−2φi,jm+φi,j−1m.The error initial and boundary conditions are given as(14)φ0,jm=φM1,jm=φi,0m=φi,M2m=φi,j0=0.Define the following grid functions form=1,2…,N:(15)φmx,y=φi,jm,whenxi−Δx2<x≤xi+Δx2,yj−Δy2<y≤yi+Δy2,0,when0≤x≤Δx2orL−Δx2≤x≤L,0,when0≤y≤Δy2orL−Δy2≤y≤L..Then,φmx,y can be expanded in Fourier series such as(16)φmx,y=∑l1,l2=−ααΧml1,l2e−12πl1xL+l2yL,where(17)Χml1,l2=1L∫0L∫0Lφmx,ye−1−2πl1xL+l2yLdxdy.From the definition ofl2 norm and Parseval equality, we have(18)φmα2=∑i=1M1−1∑j=1M2−1ΔxΔyφi,jm2=∑l1,l2=−ααΧml1,l22.Suppose that(19)φi,jm=Χme−1α1iΔx+α2jΔy,where α1=2πl1/Landα2=2πl2/L, and substituting (19) in (13), we getΧme−1α1iΔx+α2jΔy−Χme−1α1iΔx+α2jΔy=R1Χme−1α1i+1Δx+α2jΔy−2Χme−1α1iΔx+α2jΔy+Χme−1α1i−1Δx+α2jΔy−R1dm−1βΧ0e−1α1i+1Δx+α2jΔy−2Χ0e−1α1iΔx+α2jΔy+Χ0e−1α1i−1Δx+α2jΔy−R1∑s=1m−1ds−1β−dsβΧm−se−1α1i+1Δx+α2jΔy−2Χm−se−1α1iΔx+α2jΔy+Χm−se−1α1i−1Δx+α2jΔy+R2Χme−1α1iΔx+α2j+1Δy−2Χme−1α1iΔx+α2jΔy+Χme−1α1iΔx+α2j−1Δy+R2dm−1βΧ0e−1α1iΔx+α2jΔy−2Χ0e−1α1iΔx+α2jΔy+Χ0e−1α1iΔx+α2jΔy−R2∑s=1m−1ds−1β−dsβΧm−se−1α1iΔx+α2j+1Δy−2Χm−se−1α1iΔx+α2jΔy+Χm−se−1α1iΔx+α2j−1Δy+R3Χme−1α1i+1Δx+α2jΔy−2Χme−1α1iΔx+α2jΔy+Χme−1α1i−1Δx+α2jΔy+R4Χme−1α1iΔx+α2j+1Δy−2Χme−1α1iΔx+α2jΔy+Χme−1α1iΔx+α2j−1Δy.After simplifying, we get(20)Χm1+ν1+ν2=Χm−1+Χ0dm−1βν1+ν1∑s=1m−1ds−1β−dsβΧm−s,Χm=Χm−1+Χ0dm−1βν1+ν1∑s=1m−1ds−1β−dsβΧm−s1+ν1+ν2,where ν1=4R1sin²α1Δx/2+4R2sin²α2Δy/2 and ν2=4R3sin²α1Δx/2+4R4sin²α2Δy/2.Proposition 1.
IfΧmm=1,2,…,N satisfies (20), then Χm+1≤Χ0.Proof: .
By using mathematical induction, we takem=1 in (20).(21)Χ1=1+d0βν1Χ01+ν1+ν2,and as ν1,ν2≥0,b0β=1,then(22)Χ1≤Χ0.
Now, assume that(23)Χn≤Χ0;n=1,2,…,m−1,and as 0 <β<1, from (20) and Lemma 2, we obtain(24)Χm≤Χm−1+bm−1βν1Χ0+ν1∑s=1m−1ds−1β−dsβΧm−s1+ν1+ν2,≤1+dm−1βν1+ν1∑s=1m−1ds−1β−dsβ1+ν1+ν2Χ0,=1+dm−1βν1+ν21−dm−1β1+ν1+ν2Χ0,=1+ν11+ν1+ν2Χ0,Χm≤Χ0.
This completes the proof.
Based on the above proof, it can be summarized that the solution of (5) satisfies the following inequality:
Χm₂≤Χ0₂.
And, we demonstrated that the proposed scheme is unconditionally stable.


### 2.2. Convergence

Here, we use a similar method to examine the convergence of the scheme. LetYxi,yj,tm represent the exact solution; then, the truncation error of the scheme is obtained as follows: from (3),(25)Ti,jm=Yxi,yj,tm−Yxi,yj,tm−1−R1∑j=om−1dsβδx2Yxi,yj,tm−s−Yxi,yj,tm−s−1,+R2∑j=ok−1dsβδy2Yxi,yj,tm−s−Yxi,yj,tm−s−1+R3δx2Yxi,yj,tm+R4δy2Yxi,yj,tm−τhxi,yj,tm.From (1), we have(26)Ti,jm=Yi,jm−Yi,jm−1τ−∂Yxi,yj,tm∂t+∂2Yxi,yj,tm∂x2−R1∑s=0m−1dsβδx2Yi,jm−s−Yi,jm−s−1+∂2Yxi,yj,tm∂y2−R2∑s=0m−1bsβδy2Yi,jm−s−Yi,jm−s−1+∂2Yxi,yj,tm∂x2−R3δx2Yi,jm+∂2Yxi,yj,tm∂x2−R4δy2Yi,jm=Oτ+τΔx²+τΔy².Sincei,j,andm are finite, there is a positive constant C1, for all i,j,andm, which then have(27)Ti,jm≤C1τ+τΔx²+τΔy².The error is defined as(28)ϕi,jm=Yxi,yj,tm−Yi,jm.From (25), we have(29)Yxi,yj,tm=Yxi,yj,tm−1+R1Yxi,yj,tm−2Yxi,yj,tm+Yxi,yj,tm−R1dm−1βYxi+1,yj,t0−2Yxi,yj,t0+Yxi−1,yj,t0−R1∑s=1m−1ds−1β−dsβYxi+1,yj,tm−s−2Yxi,yj,tm−s−Yxi−1,yj,tm−s+R2Yxi,yj+1,tm−2Yxi,yj,tm+Yxi,yj−1,tm−R2dm−1βYxi,yj+1,tm−2Yxi,yj,tm+Yxi,yj−1,tm−R2∑s=1m−1ds−1β−dsβYxi,yj+1,tm−s−2Yxi,yj,tm−s+Yxi,yj−1,tm−s+R3Yxi+1,yj+1,tm−2Yxi,yj+1,tm+Yxi−1,yj+1,tm+R4Yxi,yj+1,tm−2Yxi,yj+1,tm+Yxi,yj−1,tm+τhxi,yj,tm.To obtain the error equation, subtract (29) from (5) to obtain(30)ϕi,jm−ϕi,jm−1=R1ϕi+1,jm−2ϕi,jm−1+ϕi−1,jm−R1dm−1βϕi+1,jm−2ϕi,jm−1+ϕi−1,jm−R1∑s=1m−1ds−1β−dsβϕi+1,jm−2ϕi,jm−1+ϕi−1,jm+R2ϕi,jm−2ϕi,jm−1+ϕi,jm−R2dm−1βϕi,j+10−2ϕi,j0+ϕi,j−10−R2∑s=1m−1ds−1β−dsβϕi,j+1m−s2−ϕi,jm−s+ϕi,j−1m−s+R3ϕi+1,jm−2ϕi,jm−1+ϕi−1,jm+R4ϕi,j+1m−s−2ϕi,jm−s+ϕi,j−1m−s+τTi,jm.With error boundary conditions,(31)ϕ0,jm=ϕM1,jm=ϕ0,jm=ϕi,M2m=0,m=1,2,…,N.And, the initial condition(32)ϕi,j0=0,i=1,2,…,M1−1,j=1,2,…,M2−1.Next, we define the following grid functions form=1,2,…,N:(33)ϕmx,y=ϕi,jm,whenxi−Δx2<x≤xi+Δx2,yj−Δy2y≤yj+Δy2,0,when0≤x≤Δx2orL−Δx2≤x≤L,0,when0≤y≤Δy2orL−Δy2≤y≤L..(34)Tmx,y=Ti,jm,whenxi−Δx2<x≤xi+Δx2,yj−Δy2<y≤yj+Δy2,0,when0≤x≤Δx2orL−Δx2≤x≤L,0,when0≤y≤Δy2orL−Δy2≤y≤L..Here, theϕmx,yandTmx,y can be expanded in Fourier series such as(35)ϕmx,y=∑l1,l2=−ααζml1,l2e2−1πl1x/L+l2y/L,m=1,2,…,N,Tmx,y=∑l1,l2=−ααφml1,l2e2−1πl1x/L+l2y/L,m=1,2,…,N,where(36)ζml1,l2=1L∫0L∫0Lϕmx,ye2−1πl1xL+l2yLdxdy,(37)φml1,l2=1L∫0L∫0Lϕmx,ye2−1πl1x/L+l2y/Ldxdy.From the definition ofl2 norm and the Parseval equality, we have(38)ϕm=l22∑i=1M1−1∑j=1M2−1ΔxΔyei,jm²=∑l1,l2=−ααρml1,l22,(39)Tm=l22∑i=1M1−1∑j=1M2−1ΔxΔyei,jm²=∑l1,l2=−ααφml1,l22.Based on the previous equations, suppose that(40)φim=ζme−1α1iΔx+α2iΔy,(41)Tim=φme−1α1iΔx+α2iΔy.Respectively, we haveα1=2πl1/Landα2=2πl2/L; substitute (40) and (41) into (30) and we getζme−1α1iΔx+α2jΔy−ζme−1α1iΔx+α2jΔy=R1ζme−1α1i+1Δx+α2jΔy−2ζme−1α1iΔx+α2jΔy+ζme−1α1i−1Δx+α2jΔy−R1dm−1βζme−1α1i+1Δx+α2jΔy−2ζme−1α1iΔx+α2jΔy+ζme−1α1i−1Δx+α2jΔy−R1∑s=1m−1ds−1β−dsβζme−1α1i+1Δx+α2jΔy−2ζme−1α1iΔx+α2jΔy+ζme−1α1i−1Δx+α2jΔy+R2ζme−1α1iΔx+α2j+1Δy−2ζme−1α1iΔx+α2jΔy+ζme−1α1iΔx+α2j−1Δy+R2dm−1βζme−1α1iΔx+α2jΔy−2ζme−1α1iΔx+α2jΔy+ζme−1α1iΔx+α2jΔy−R2∑s=1m−1ds−1β−dsβζme−1α1iΔx+α2j+1Δy−2ζme−1α1iΔx+α2jΔy+ζme−1α1iΔx+α2j−1Δy+R3ζme−1α1i+1Δx+α2jΔy−2ζme−1α1iΔx+α2jΔy+ζme−1α1i−1Δx+α2jΔy+R4ζme−1α1iΔx+α2j+1Δy−2ζme−1α1iΔx+α2jΔy+ζme−1α1iΔx+α2j−1Δy+τφme−1α1iΔx+α2jΔy,Simplifying the previous equation, we obtain(42)ζm=ζm−1+ζ0bm−1βν1+ν1∑s=1m−1ds−1β−dsβζm−s1+ν1+ν2,where(43)ν1=4R1sin²α1Δx2+4R2sin²α2Δy2,ν2=4R3sin²α1Δx2+4R4sin²α2Δy2.Proposition 2.
Letζmm=1,2,…,N be the solution of (42); then, there is a positive constant C2 so thatζm≤C2mτφ1.Proof: .
Fromϕ0=0 and (36), we have(44)ζo=ζol1,l2=0.
From (37) and (39), then there is a positive constant C2, such that(45)φm≤C2φ1l1,l2.
Using mathematical induction, form=1, then from (42) and (44), we obtain(46)ζ1=1/1+ν1+ν2τφ1.
Sinceν1,ν2≥0, from (45), we get(47)ζ1≤τφ1≤C2τφ1.
Now, suppose that(48)ζm≤C2mτφ1,n=1,2,…,m−1.
As 0<β<1,ν1,ν2≥0.
From (41) and (44) and Lemma 2, we have(49)ζm=ζm−1+ν1∑s=1m−1ds−1β−dsβζm−s+τφm1+ν1+ν2,ζm=C2m−1τφ1+ν1∑s=1m−1ds−1β−dsβC2m−sτφ1+C2τφ11+ν1+ν2,≤m−1+ν1m−1∑s=1m−1bs−1β−bsβ+11+ν1+ν2C2τφ1,=m+ν1m−1∑s=1m−1bs−1β−bsβ+11+ν1+ν2C2τφ1,=m+ν1m−1+1−bm−1β1+ν1+ν2C2τφ1,≤mC2τφ1.
The proof is completed via the induction method.Theorem 1.
The modified implicit difference schemel2 is convergent, and the order of convergence is Oτ+τΔx2+τΔy2.Proof: .
From (27) and (39), we obtain(50)Tk≤MxΔxMyΔyC1τ+τΔx2+τΔy2=LC1τ+τΔx2+τΔy2φkl2≤kC2τT1≤C1C2kτLτ+τΔx2+τΔy2.
Askτ≤T, thus(51)φkl2≤C1C2kτLτ+τΔx2+τΔy2,where C=C1C2TL.


## 2.1. Stability

We find the stability of the proposed scheme by Fourier technique. Let the approximate solution beΨi,jm for (9); we have(11)Ψi,jm−Ψi,jm−1=R1Ψi+1,jm−2Ψi,jm+Ψi−1,jm−R1dm−1βΨi+1,j0−2Ψi,j0+Ψi−1,j0−R1∑s=1m−1ds−1β−dsβΨi+1,jm−s−2Ψi,jm−s+Ψi,j+1m−s+R2Ψi,jm−s−2Ψi,jm−s+Ψi,j−1m−s+R2dm−1βΨi,j+10−2Ψi,j0+Ψi,j−10−R2∑s=1m−1ds−1β−dsβΨi,j+1m−s−2Ψi,jm−s+Ψi,j−1m−s+R3Ψi+1,jm−2Ψi,jm+Ψi−1,jm+R4Ψi,j+1m−2Ψi,jm+Ψi,j−1m.Next, the error is defined as(12)φi,jm=Yi,jm−Ψi,jm.where φi,jm satisfies (11) and(13)φi,jm−φi,jm−1=R1φi+1,jm−2φi,jm+φi−1,jm−R1dm−1βφi+1,j0−2φi,j0+φi−1,j0−R1∑s=1m−1ds−1β−dsβφi+1,jm−s−2φi,jm−s+φi,j+1m−s+R2φi,j+1m−s−2φi,jm−s+φi,j−1m−s+R2dm−1βφi,j+10−2φi,j0+φi,j−10−R2∑s=1m−1ds−1β−dsβφi,j+1m−s−2φi,jm−s+φi,j−1m−s+R3φi+1,jm−2φi,jm+φi−1,jm+R4φi,j+1m−2φi,jm+φi,j−1m.The error initial and boundary conditions are given as(14)φ0,jm=φM1,jm=φi,0m=φi,M2m=φi,j0=0.Define the following grid functions form=1,2…,N:(15)φmx,y=φi,jm,whenxi−Δx2<x≤xi+Δx2,yj−Δy2<y≤yi+Δy2,0,when0≤x≤Δx2orL−Δx2≤x≤L,0,when0≤y≤Δy2orL−Δy2≤y≤L..Then,φmx,y can be expanded in Fourier series such as(16)φmx,y=∑l1,l2=−ααΧml1,l2e−12πl1xL+l2yL,where(17)Χml1,l2=1L∫0L∫0Lφmx,ye−1−2πl1xL+l2yLdxdy.From the definition ofl2 norm and Parseval equality, we have(18)φmα2=∑i=1M1−1∑j=1M2−1ΔxΔyφi,jm2=∑l1,l2=−ααΧml1,l22.Suppose that(19)φi,jm=Χme−1α1iΔx+α2jΔy,where α1=2πl1/Landα2=2πl2/L, and substituting (19) in (13), we getΧme−1α1iΔx+α2jΔy−Χme−1α1iΔx+α2jΔy=R1Χme−1α1i+1Δx+α2jΔy−2Χme−1α1iΔx+α2jΔy+Χme−1α1i−1Δx+α2jΔy−R1dm−1βΧ0e−1α1i+1Δx+α2jΔy−2Χ0e−1α1iΔx+α2jΔy+Χ0e−1α1i−1Δx+α2jΔy−R1∑s=1m−1ds−1β−dsβΧm−se−1α1i+1Δx+α2jΔy−2Χm−se−1α1iΔx+α2jΔy+Χm−se−1α1i−1Δx+α2jΔy+R2Χme−1α1iΔx+α2j+1Δy−2Χme−1α1iΔx+α2jΔy+Χme−1α1iΔx+α2j−1Δy+R2dm−1βΧ0e−1α1iΔx+α2jΔy−2Χ0e−1α1iΔx+α2jΔy+Χ0e−1α1iΔx+α2jΔy−R2∑s=1m−1ds−1β−dsβΧm−se−1α1iΔx+α2j+1Δy−2Χm−se−1α1iΔx+α2jΔy+Χm−se−1α1iΔx+α2j−1Δy+R3Χme−1α1i+1Δx+α2jΔy−2Χme−1α1iΔx+α2jΔy+Χme−1α1i−1Δx+α2jΔy+R4Χme−1α1iΔx+α2j+1Δy−2Χme−1α1iΔx+α2jΔy+Χme−1α1iΔx+α2j−1Δy.After simplifying, we get(20)Χm1+ν1+ν2=Χm−1+Χ0dm−1βν1+ν1∑s=1m−1ds−1β−dsβΧm−s,Χm=Χm−1+Χ0dm−1βν1+ν1∑s=1m−1ds−1β−dsβΧm−s1+ν1+ν2,where ν1=4R1sin²α1Δx/2+4R2sin²α2Δy/2 and ν2=4R3sin²α1Δx/2+4R4sin²α2Δy/2.Proposition 1.
IfΧmm=1,2,…,N satisfies (20), then Χm+1≤Χ0.Proof: .
By using mathematical induction, we takem=1 in (20).(21)Χ1=1+d0βν1Χ01+ν1+ν2,and as ν1,ν2≥0,b0β=1,then(22)Χ1≤Χ0.
Now, assume that(23)Χn≤Χ0;n=1,2,…,m−1,and as 0 <β<1, from (20) and Lemma 2, we obtain(24)Χm≤Χm−1+bm−1βν1Χ0+ν1∑s=1m−1ds−1β−dsβΧm−s1+ν1+ν2,≤1+dm−1βν1+ν1∑s=1m−1ds−1β−dsβ1+ν1+ν2Χ0,=1+dm−1βν1+ν21−dm−1β1+ν1+ν2Χ0,=1+ν11+ν1+ν2Χ0,Χm≤Χ0.
This completes the proof.
Based on the above proof, it can be summarized that the solution of (5) satisfies the following inequality:
Χm₂≤Χ0₂.
And, we demonstrated that the proposed scheme is unconditionally stable.


## 2.2. Convergence

Here, we use a similar method to examine the convergence of the scheme. LetYxi,yj,tm represent the exact solution; then, the truncation error of the scheme is obtained as follows: from (3),(25)Ti,jm=Yxi,yj,tm−Yxi,yj,tm−1−R1∑j=om−1dsβδx2Yxi,yj,tm−s−Yxi,yj,tm−s−1,+R2∑j=ok−1dsβδy2Yxi,yj,tm−s−Yxi,yj,tm−s−1+R3δx2Yxi,yj,tm+R4δy2Yxi,yj,tm−τhxi,yj,tm.From (1), we have(26)Ti,jm=Yi,jm−Yi,jm−1τ−∂Yxi,yj,tm∂t+∂2Yxi,yj,tm∂x2−R1∑s=0m−1dsβδx2Yi,jm−s−Yi,jm−s−1+∂2Yxi,yj,tm∂y2−R2∑s=0m−1bsβδy2Yi,jm−s−Yi,jm−s−1+∂2Yxi,yj,tm∂x2−R3δx2Yi,jm+∂2Yxi,yj,tm∂x2−R4δy2Yi,jm=Oτ+τΔx²+τΔy².Sincei,j,andm are finite, there is a positive constant C1, for all i,j,andm, which then have(27)Ti,jm≤C1τ+τΔx²+τΔy².The error is defined as(28)ϕi,jm=Yxi,yj,tm−Yi,jm.From (25), we have(29)Yxi,yj,tm=Yxi,yj,tm−1+R1Yxi,yj,tm−2Yxi,yj,tm+Yxi,yj,tm−R1dm−1βYxi+1,yj,t0−2Yxi,yj,t0+Yxi−1,yj,t0−R1∑s=1m−1ds−1β−dsβYxi+1,yj,tm−s−2Yxi,yj,tm−s−Yxi−1,yj,tm−s+R2Yxi,yj+1,tm−2Yxi,yj,tm+Yxi,yj−1,tm−R2dm−1βYxi,yj+1,tm−2Yxi,yj,tm+Yxi,yj−1,tm−R2∑s=1m−1ds−1β−dsβYxi,yj+1,tm−s−2Yxi,yj,tm−s+Yxi,yj−1,tm−s+R3Yxi+1,yj+1,tm−2Yxi,yj+1,tm+Yxi−1,yj+1,tm+R4Yxi,yj+1,tm−2Yxi,yj+1,tm+Yxi,yj−1,tm+τhxi,yj,tm.To obtain the error equation, subtract (29) from (5) to obtain(30)ϕi,jm−ϕi,jm−1=R1ϕi+1,jm−2ϕi,jm−1+ϕi−1,jm−R1dm−1βϕi+1,jm−2ϕi,jm−1+ϕi−1,jm−R1∑s=1m−1ds−1β−dsβϕi+1,jm−2ϕi,jm−1+ϕi−1,jm+R2ϕi,jm−2ϕi,jm−1+ϕi,jm−R2dm−1βϕi,j+10−2ϕi,j0+ϕi,j−10−R2∑s=1m−1ds−1β−dsβϕi,j+1m−s2−ϕi,jm−s+ϕi,j−1m−s+R3ϕi+1,jm−2ϕi,jm−1+ϕi−1,jm+R4ϕi,j+1m−s−2ϕi,jm−s+ϕi,j−1m−s+τTi,jm.With error boundary conditions,(31)ϕ0,jm=ϕM1,jm=ϕ0,jm=ϕi,M2m=0,m=1,2,…,N.And, the initial condition(32)ϕi,j0=0,i=1,2,…,M1−1,j=1,2,…,M2−1.Next, we define the following grid functions form=1,2,…,N:(33)ϕmx,y=ϕi,jm,whenxi−Δx2<x≤xi+Δx2,yj−Δy2y≤yj+Δy2,0,when0≤x≤Δx2orL−Δx2≤x≤L,0,when0≤y≤Δy2orL−Δy2≤y≤L..(34)Tmx,y=Ti,jm,whenxi−Δx2<x≤xi+Δx2,yj−Δy2<y≤yj+Δy2,0,when0≤x≤Δx2orL−Δx2≤x≤L,0,when0≤y≤Δy2orL−Δy2≤y≤L..Here, theϕmx,yandTmx,y can be expanded in Fourier series such as(35)ϕmx,y=∑l1,l2=−ααζml1,l2e2−1πl1x/L+l2y/L,m=1,2,…,N,Tmx,y=∑l1,l2=−ααφml1,l2e2−1πl1x/L+l2y/L,m=1,2,…,N,where(36)ζml1,l2=1L∫0L∫0Lϕmx,ye2−1πl1xL+l2yLdxdy,(37)φml1,l2=1L∫0L∫0Lϕmx,ye2−1πl1x/L+l2y/Ldxdy.From the definition ofl2 norm and the Parseval equality, we have(38)ϕm=l22∑i=1M1−1∑j=1M2−1ΔxΔyei,jm²=∑l1,l2=−ααρml1,l22,(39)Tm=l22∑i=1M1−1∑j=1M2−1ΔxΔyei,jm²=∑l1,l2=−ααφml1,l22.Based on the previous equations, suppose that(40)φim=ζme−1α1iΔx+α2iΔy,(41)Tim=φme−1α1iΔx+α2iΔy.Respectively, we haveα1=2πl1/Landα2=2πl2/L; substitute (40) and (41) into (30) and we getζme−1α1iΔx+α2jΔy−ζme−1α1iΔx+α2jΔy=R1ζme−1α1i+1Δx+α2jΔy−2ζme−1α1iΔx+α2jΔy+ζme−1α1i−1Δx+α2jΔy−R1dm−1βζme−1α1i+1Δx+α2jΔy−2ζme−1α1iΔx+α2jΔy+ζme−1α1i−1Δx+α2jΔy−R1∑s=1m−1ds−1β−dsβζme−1α1i+1Δx+α2jΔy−2ζme−1α1iΔx+α2jΔy+ζme−1α1i−1Δx+α2jΔy+R2ζme−1α1iΔx+α2j+1Δy−2ζme−1α1iΔx+α2jΔy+ζme−1α1iΔx+α2j−1Δy+R2dm−1βζme−1α1iΔx+α2jΔy−2ζme−1α1iΔx+α2jΔy+ζme−1α1iΔx+α2jΔy−R2∑s=1m−1ds−1β−dsβζme−1α1iΔx+α2j+1Δy−2ζme−1α1iΔx+α2jΔy+ζme−1α1iΔx+α2j−1Δy+R3ζme−1α1i+1Δx+α2jΔy−2ζme−1α1iΔx+α2jΔy+ζme−1α1i−1Δx+α2jΔy+R4ζme−1α1iΔx+α2j+1Δy−2ζme−1α1iΔx+α2jΔy+ζme−1α1iΔx+α2j−1Δy+τφme−1α1iΔx+α2jΔy,Simplifying the previous equation, we obtain(42)ζm=ζm−1+ζ0bm−1βν1+ν1∑s=1m−1ds−1β−dsβζm−s1+ν1+ν2,where(43)ν1=4R1sin²α1Δx2+4R2sin²α2Δy2,ν2=4R3sin²α1Δx2+4R4sin²α2Δy2.Proposition 2.
Letζmm=1,2,…,N be the solution of (42); then, there is a positive constant C2 so thatζm≤C2mτφ1.Proof: .
Fromϕ0=0 and (36), we have(44)ζo=ζol1,l2=0.
From (37) and (39), then there is a positive constant C2, such that(45)φm≤C2φ1l1,l2.
Using mathematical induction, form=1, then from (42) and (44), we obtain(46)ζ1=1/1+ν1+ν2τφ1.
Sinceν1,ν2≥0, from (45), we get(47)ζ1≤τφ1≤C2τφ1.
Now, suppose that(48)ζm≤C2mτφ1,n=1,2,…,m−1.
As 0<β<1,ν1,ν2≥0.
From (41) and (44) and Lemma 2, we have(49)ζm=ζm−1+ν1∑s=1m−1ds−1β−dsβζm−s+τφm1+ν1+ν2,ζm=C2m−1τφ1+ν1∑s=1m−1ds−1β−dsβC2m−sτφ1+C2τφ11+ν1+ν2,≤m−1+ν1m−1∑s=1m−1bs−1β−bsβ+11+ν1+ν2C2τφ1,=m+ν1m−1∑s=1m−1bs−1β−bsβ+11+ν1+ν2C2τφ1,=m+ν1m−1+1−bm−1β1+ν1+ν2C2τφ1,≤mC2τφ1.
The proof is completed via the induction method.Theorem 1.
The modified implicit difference schemel2 is convergent, and the order of convergence is Oτ+τΔx2+τΔy2.Proof: .
From (27) and (39), we obtain(50)Tk≤MxΔxMyΔyC1τ+τΔx2+τΔy2=LC1τ+τΔx2+τΔy2φkl2≤kC2τT1≤C1C2kτLτ+τΔx2+τΔy2.
Askτ≤T, thus(51)φkl2≤C1C2kτLτ+τΔx2+τΔy2,where C=C1C2TL.


## 3. Numerical Experiment

Example 1.
Consider the following two-dimensional Rayleigh–Stokes problem for heated generalized second-grade fluid with the fractional derivative [22]:(52)∂Yx,y,t∂t=0Dt1−β∂2Yx,y,t∂x2+∂²Yx,y,t∂y2+∂2Yx,y,t∂x2+∂2Yx,y,t∂y2+hx,y,t,0≤β≤1,0≤t≤T,with initial and boundary conditions(53)Yx,y,0=0,0≤x,y≤1,Y0,y,t=eyt1+β,Y1,y,t=e1+yt1+β,Yx,0,t=ext1+β,Yx,1,t=e1+xt1+β,0≤t≤T.
Here,hx,y,t=1+βtβ−2Γ2+β/Γ1+2βt2β−2t1+βex+y and the exact solution of (52) is given by(54)Yx,y,t=ex+yt1+β.
The error between the numerical solution and exact solution is defined as follows:(55)E∞=max0≤i,j≤M,0≤m≤NYxi,yj,tm−Yi,jm.
And, the rate of convergence for space variable can be defined as(56)=−order=log2E∞16τ,2ΔxE∞τ,Δx.
The developed modified implicit scheme is applied to problems (52) to (54).
Tables1–4 show the errors E∞ for values of space step size Δx,Δy and τ. Here, time step τ is defined by τ=T/N.
Tables1–4 indicate that, as we reduce the time and space step size τ and Δx,Δy, the error decreases for a fixed value of γ.This shows that the method converges to the exact solution.Example 2.
Consider the following Cable equation:(57)∂Yx,t∂t=D0t1−ρ1K∂2Yx,t∂x2−μ2D0t1−ρ2Yx,t+2t+π2t1+ρ1Γ2+ρ1+t1+ρ2Γ2+ρ2sinπx,with initial and boundary conditions(58)Yx,0=0,0≤x≤L,Y0,t=t2sinπx,YL,t=β2t,0t≤T.The exact solution is Yx,t=t2sinπx.


## 4. Results and Discussion

A modified implicit scheme is developed and applied on RSP-HGSGF. Numerical example is given to support theoretical study. The error between the exact and numerical solution is calculated using different values ofN and M. Also, at different values of γ, Tables 1–4 are created to show the comparison of the numerical scheme with the exact solution in terms of maximum error. In example 2, we solved the fractional-order Cable equation, and the numerical results are shown in Table 5 for various values of space and time step size. The values of ρ1andρ2 are also changed, and the obtained results are converging with reduced step sizes. Here, the error is calculated using Maple15 software with the increase in the number of space and time steps. Figures 1–3 are plotted for different values of MandN, and fractional order γ shows good agreement with the exact solution.Table 1
The error table for different values atτ,Δx,Δy,andγ.
τΔx=Δyγ=0.5γ=0.6γ=0.7γ=0.8γ=0.91/41/22.180E − 21.316E − 21.601E − 21.889E − 22.481E − 21/161/44.573E − 35.267E − 35.934E − 36.606E − 27.312E − 21/641/81.330E − 31.484E − 31.634E − 31.789E − 31.958E − 31/1281/107.701E − 48.426E − 49.141E − 49.898E − 41.073E − 4Table 2
The error table for different values atτ,Δx,Δy,andγ.
τ=Δx=Δyγ=0.35γ=0.65γ=0.851/27.9899E − 32.4671E − 23.8181E − 2¼4.5383E − 31.4258E − 22.1050E − 21/63.3873E − 39.9350E − 31.4272E − 21/82.8140E − 37.7703E − 31.0965E − 21/102.4289E − 36.3977E − 38.9159E − 3Table 3
The error table for different values atτ,Δx,Δy,andγ.
τΔx=Δyγ=0.5γ=0.6γ=0.7γ=0.8γ=0.91/161/45.0794E − 35.9521E − 36.7049E − 30.007463528.2610E − 31/81/81.4043E − 31.6387E − 38.5596E − 31.0150E − 31.8575E − 31/1441/126.2458E − 46.8873E − 47.5239E − 48.1996E − 48.9478E − 4Table 4
The error table for different values atτ,Δx,andΔy and at a fixed value of γ=0.25.
NΔx=Δy=1/5Δx=Δy=1/10Δx=Δy=1/15Δx=Δy=1/20201.9274E − 39.5392E − 47.6704E − 47.0147E − 4401.7467E − 31.6060E − 35.7985E − 44.3020E − 4601.6504E − 36.6854E − 45.4537E − 43.2953E − 4Table 5
Numerical results of example 2 of the proposed scheme for various values ofρ1,ρ2, N, and Δx.
ΔxNρ1,ρ2=0.25ρ1,ρ2=0.5ρ1,ρ2=0.951/10401.173E − 28.131E − 26.552E − 3808.032E − 36.452E − 35.762E − 31105.912E − 35.528E − 35.171E − 31/20408.146E − 36.312E − 33.272E − 3804.537E − 32.989E − 32.409E − 31102.428E − 32.015E − 39.146E − 41/40407.294E − 33.781E − 32.491E − 3803.639E − 32.173E − 31.591E − 31101.578E − 31.257E − 39.015E − 4Figure 1
Comparing equations (52) and (54) at M=4,N=2, and γ=0.25.Figure 2
Comparing equations (52) and (54) at M=6,N=6, and γ=0.35.Figure 3
Comparing equations (52) and (54) at M=8,N=8, and γ=0.85.


## 5. Conclusion

A modiﬁed implicit diﬀerence scheme is formulated for 2D RSP-HGSGF, and a derivative of fractional order has been described in this paper. The modiﬁed scheme has the improvement of low computational cost and can be easily applied. The Fourier technique has been used for the theoretical analysis stability, and convergence with orderτ+Δx²+Δy² is unconditionally stable and convergent. The numerical experiment for the 2D RSP-HGSGF and 1D Cable equation is conducted, which shows that the modified implicit scheme is easy to implement, and the results show good performance of the proposed schemes [39].



---

*Source: 6948461-2022-05-06.xml*