High-order compact scheme for the two-dimensional fractional Rayleigh–Stokes problem for a heated generalized second-grade fluid

In this article, an unconditionally stable compact high-order iterative finite difference scheme is developed on solving the two-dimensional fractional Rayleigh–Stokes equation. A relationship between the Riemann–Liouville (R–L) and Grunwald–Letnikov (G–L) fractional derivatives is used for the time-fractional derivative, and a fourth-order compact Crank–Nicolson approximation is applied for the space derivative to produce a high-order compact scheme. The stability and convergence for the proposed method will be proven; the proposed method will be shown to have the order of convergence O(τ+h4)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$O(\tau + h^{4})$\end{document}. Finally, numerical examples are provided to show the high accuracy solutions of the proposed scheme.


Introduction
The interest in fractional calculus has increased in recent years because of its applicability in various fields of science and engineering [1][2][3][4][5][6][7][8][9][10]. Many physical phenomena from fluid mechanics, viscoelasticity, physics and other sciences can be modeled mathematically with the help of fractional differential equations (FDEs) which represents the physical phenomena more appropriately than ordinary differential equations. In most cases, it is difficult to find the analytical solution of the FDE due to its complexity. Therefore researchers resorted to different numerical methods with a variety of stability and convergence properties to solve various FDEs [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. The Rayleigh-Stokes problem (RSP) is one of the most extensively researched FDEs over the last few years. This equation plays an important role in representing the behavior of some non-Newtonian fluids and the fractional derivative is also used in this model problem to capture the visco-elastic behavior of the flow [26,27]. Also, fractional calculus has described the visco-elastic model more accurately than the non-fractional model, which is why we considered the Rayleigh-Stokes problem with fractional derivative.
In this paper, we considered the two-dimensional (2D) RSP with fractional derivative of the form: ∂w(x, y, t) ∂t = 0 D 1-γ t ∂ 2 w(x, y, t) ∂x 2 + ∂ 2 w(x, y, t) ∂y 2 + ∂ 2 w(x, y, t) ∂x 2 + ∂ 2 w(x, y, t) having initial and boundary conditions w(x, y, t) = g(x, y, t), (x, y) ∈ ∂Ω, where 0 < γ < 1, Ω = {(x, y) | 0 ≤ x ≤ L, 0 ≤ y ≤ L} and 0 D 1-γ t represents the Riemann-Liouville fractional derivative of order, which is define as The 2D RSP defined in (1) has been solved by several researchers such as Chen et al. [28] solved the RSP using implicit and explicit of finite differences methods; also Fourier analysis is used for the convergence and stability of the proposed schemes. The convergence of both schemes were found to be of order O(τ + x 2 + y 2 ). Zhuang and Liu [29] solved the same problem by an implicit numerical approximation scheme, and its stability and convergence were also established with the order of convergence O(τ + x 2 + y 2 ). Mohebbi et al. [30] used a higher-order implicit finite difference scheme for (1), and they discussed its convergence and stability by Fourier analysis. The convergence order of their scheme was shown to be O(τ + x 4 + y 4 ). Although available, higher-order schemes for solving the two-dimensional Rayleigh-Stokes with better accuracy and simplicity are still very limited. Motivated by this, we try to further formulate another high-order stable numerical method for the Rayleigh-Stokes problem for a heated generalized second-grade fluid with fractional derivative (1) with given boundary and initial conditions (2)-(3) but with better accuracy than the schemes in [30]. This paper aims to propose an unconditionally stable compact iterative scheme for solving (1) having order of convergence O(τ + x 4 + y 4 ). A fourth-order Crank-Nicolson difference scheme is applied for the discretization of the space derivative and a relationship between the Grunwald-Letnikov, and the Riemann-Liouville fractional derivative is used for the fractional time derivative. The stability and convergence analysis of this proposed method will be established using Fourier analysis. The paper is organized as follows: in Sect. 2, we discuss the formulation of the proposed scheme, followed by the stability of the proposed scheme in Sect. 3. The convergence analysis is discussed in Sect. 4. In Sect. 5, some numerical examples are presented with discussion, and finally, the conclusion is presented in Sect. 6.

Stability of the proposed scheme
Let w k i,j (i, j = 0, 1, 2, . . . , n, k = 0, 1, 2, . . . , l) be the approximate solution and W k i,j be the exact solution of (1), then the error ε k i,j = W k i,jw k i,j will also satisfy (1), so from (12) where . . , k + 1. Assume the error function at initial and boundary value are defined by and the error function for the grid function when k = 0, 1, 2, . . . , N -1 is defined as Then we can express ε k (x, y) by a Fourier series as where The relationship between the Parseval equality and the L 2 norm is Now suppose that (13), and after simplification we get where λ 1 = Cos(β x) + Cos(α y) and λ 2 = Cos(β x) Cos(α y).

Theorem 1 The C-N high-order compact finite difference scheme
Proof Using (17) and Proposition 1, we have This shows C-N high-order compact finite difference scheme (12) is unconditionally stable.
Hence, we proved the proposition.

Numerical results
In this section, we will show the effectiveness of our proposed method by performing three numerical experiments on the two-dimensional fractional Rayleigh-Stokes problem with the help of Core i7 Duo 3.40 GHz, Window 7 and 4 GB RAM using Mathematica software. We used Successive Over-Relaxation (SOR) as the acceleration technique, also for the convergence criterion a tolerance ζ = 10 -5 is used for the maximum error (L ∞ ).
The convergence order for the space variable and time variable is calculated with the help of the C 2 -order of convergence and C 1 -order of convergence, respectively [31], where h is the space step, τ is the time step and L ∞ is the maximum error = e l ∞ = max 1≤i,j≤N-1 |W k i,jw k i,j |. The following test problems were used for the numerical experiments.

Test Problem 3 ([28])
∂w(x, y, t) where 0 < x, y < 1, and having initial and boundary conditions w(x, y, 0) = 0, with exact solution We solved Test problem 1, Test problem 2 and Test problem 3 using our proposed method for the different values of γ , τ , h and L = T = 1. Table 1 and Table 2 shows the C 1 -order of convergence for Test problem 1 and Test problem 3, respectively. Similarly,     Table 3 and Table 4 show C 2 -order of convergence for the Test problem 1 and Test problem 3, respectively, where Max error is L ∞ error. From Table 1 to Table 4, it can be observed that the computational order of convergence is in agreement with the theoretical order of convergence, i.e. Order of convergence is O(τ + h 4 ). Table 5, Table 6 and Table 7 show the comparison between the Implicit compact method, the Radial basis function Meshless Method (RMM) and our proposed method for (τ = 1 50 , γ = 0.25, ν = 0.1), (h = 1 20 , τ = 1 50 , ν = 1 30 ) and (γ = 0.25, ν = 1 30 ), respectively, from which it can be seen that our proposed    Figure 1 The convergence of the proposed method compared with the implicit method and RMM for Test problem 2 method gives better results as compared to the Implicit compact method and the Radial basis function Meshless Method (RMM), also it can be seen in Fig. 1 and Fig. 2. Similarly, Tables 8-11 show the number of iterations, execuation time and errors (maximum and average error) for the Test problem 1 and Test problem 2. Figures 3-5 are the 3-D graphs of L ∞ error, the approximate solution and exact solution for the Test problem 1, respectively, when (h = τ = 1 30 , γ = 0.5), and similarly Figs. 6-8, represent the 3-D graphs of L ∞ error, the approximate solution and the exact solution for the Test problem 1, respectively, when (h = τ = 1 25 , γ = 0.5, ν = 0.1). From Fig. 4, Fig. 5, Fig. 7 and Fig. 8, it can be seen that Fig. 5 and Fig. 8 (the approximate solutions) are good approximations to the exact solutions ( Fig. 4 and Fig. 7).

Conclusion
In this article, we solved the two-dimensional fractional RSP using a compact high-order finite difference method. We proved that the Crank-Nicolson compact method is better in accuracy as compared to the implicit compact and radial basis function meshless method. The proposed method is unconditionally stable and convergent. The C 1 -order of convergence and C 2 -order of convergence show that the theoretical and experimental order of convergence agree for the time and space variables, respectively, i.e. O(τ + h 4 ).