Jacobi Spectral Galerkin Method for Distributed-Order Fractional Rayleigh–Stokes Problem for a Generalized Second Grade Fluid

Distributed-order fractional differential operators provide a powerful tool for mathematical modeling of multiscale multiphysics processes, where the differential orders are distributed over a range of values rather than being just a fixed fraction. In this work, we consider the Rayleigh-Stokes problem for a generalized second-grade fluid which involves the distributed-order fractional derivative in time. We develop a spectral Galerkin method for this model by employing Jacobi polynomials as temporal and spatial basis/test functions. The suggested approach is based on a novel distributed order fractional differentiation matrix for Jacobi polynomials. Numerical results for one- and two-dimensional examples are presented illustrating the performance of the algorithm. The results show that our scheme can achieve the spectral accuracy for the problem under consideration with smooth solution and allows a great flexibility to deal with multi-dimensional temporally-distributed order fractional Rayleigh-Stokes problems as the global behavior of the solution is taken into account.


INTRODUCTION
Let ⊂ R d (d = 1, 2) be a bounded convex domain with a polygonal boundary ∂ , and T > 0 be a fixed time. The purpose of this paper is to extend the approach in Hafez and Zaky [1] to the distributed-order time-fractional Rayleigh-Stokes problem for a generalized second-grade fluid in one and two space dimensions. The mathematical model is given by with initial-boundary conditions: U(x, t) x∈∂ = g(x, t), t ∈ (0, T ], where is the Laplacian operator, ∂ is the boundary of , H(x, t), f (x), and g(x, t) are given smooth functions on . The distributed-order fractional derivative 0D w(ν) t is defined by 0 where w(ν) is a non-negative smooth weight function satisfying the conditions ∂ξ dξ , 0 < ν < 1.
(5) The fractional-order derivatives appearing in (4) is defined in the Caputo sense. Unlike Riemann-Liouville derivatives, Caputo derivatives are not singular on the domain boundaries. That feature makes them particularly appealing for non-local numerical methods, like the spectral methods. Zhang et al. [2] have also confirmed that Riemann-Liouville derivatives could cause mass-balance errors on bounded domains. The time derivative is thus defined in the Caputo sense.
Recently, the fractional-order Rayleigh-Stokes problem has received considerable attention in recent years. This model plays an important role in describing the behavior of some non-Newtonian fluids [3,4]. In order to gain insights into the behavior of the solution of this model, there has been substantial interest in deriving a closed form solution for the special case w(ν) = δ(ν − 1 + β), where δ(·) is the Dirac delta function and β ∈ (0, 1). For example, Shen et al. [4] derived exact solutions of this model using the Fourier transform and the fractional Laplace transform. Girault and Saadouni [5] analyzed the existence and uniqueness of a weak solution of a closely related time-dependent grade-two fluid model. Zhao and Yang [6] obtained exact solutions using the eigenfunction expansion. Xue and Nie [7] obtained also closed form solutions of this model in a porous half-space using both Fourier and fractional Laplace transforms. The exact solutions obtained in these studies involve infinite series and special functions, e.g., generalized Mittag-Leffler functions, and thus are inconvenient for numerical evaluation. Further, closed-form solutions are available only for a restricted class of problem settings. Hence, it is imperative to develop efficient and optimally accurate numerical algorithms for problem (1). Wu [8] developed an implicit numerical scheme by transforming the above mentioned problem into an integral equation. Lin and Jiang [9] introduced a numerical sachem based on the reproducing kernel space. Mohebbi et al. [10] compared a fourth-order compact scheme with radial basis functions meshless approach. Recently, Bazhlekova et al. [11] developed two fully discrete schemes based on the backward difference method and backward Euler method and a semidiscrete scheme based on the Galerkin finite element method. Abdelkawy and Alqahtani [12] proposed spectral collocation techniques for solving fractional Stokes' first problem for a heated generalized second grade fluid. Bhrawy et al. [13] developed two shifted Jacobi-Gauss collocation schemes. More recently, Dehghan and Abbaszadeh [14] developed a finite element method for two-dimensional fractional Rayleigh-Stokes model on complex geometries. Shivanian and Jafarabadi [15] developed a spectral meshless radial point nterpolation technique. Zaky [16] develop efficient algorithms based on the Legendre-tau approximation for one-and two-dimensional fractional Rayleigh-Stokes problems for a generalized second-grade fluid. Yang and Jiang [17] proposed a numerical algorithm based on the L1 finite difference scheme for the temporal direction while the Legendre spectral method for the spatial direction.
Theoretical studies on numerical methods for fractional differential equations involving distributed-order derivatives have received considerable attention in the last decade [18][19][20]. A general form of distributed-order fractional ordinary differential equation is solved in Katsikadelis [21], Mashayekhi and Razzaghi [22], and Zaky et al. [23]. Numerical methods for solving distributed-order time-fractional diffusion equations are presented in Morgado and Rebelo [24], Abdelkawy et al. [25], and Zaky and Machado [26]. Numerical methods for distributedorder space-fractional diffusion equations are provided in Abbaszadeh [27], Kazmi and Khaliq [28], and Fan and Liu [29]. Numerical methods for multi-dimensional distributedorder generalized Schrödinger equations are provided in Bhrawy and Zaky [30]. For numerical methods for solving distributedorder fractional optimal control problems, see [31,32]. In this paper, we use a non-local representation of the solution of the distributed-order time-fractional Rayleigh-Stokes problem to introduce spectral solutions. The spectral and pseudospectral methods are well-known for their high accuracy and have been used extensively in scientific computation, see [33][34][35][36][37][38][39][40][41] and the references therein. The main contribution of this paper is to develop Jacobi-Galerkin algorithms for solving the multidimensional distributed-order time-fractional Rayleigh-Stokes problem (1).
The outline of this work is as follows. In section 2, we introduce the distributed-order fractional differentiation matrix of the shifted Jacobi polynomials. In section 3, we derive a timespace discretization for the one-dimensional distributed-order fractional Rayleigh-Stokes problem. In section 4, we consider the numerical solution of the two-dimensional case. In section 5, we present various numerical results exhibiting the efficiency of our numerical schemes. We end this paper with a few concluding remarks in section 6. with J (γ ,η) 0 .
Using the linear map z = 2x L − 1, the interval [−1, 1] can be rescaled onto [0, L]. Hence, the set of shifted Jacobi polynomials J (γ ,η) can be generated by: The terminal values of the shifted Jacobi polynomials satisfy The J (γ ,η) L,j (x) satisfies the following orthogonality relation where w The following Jacobi-Gauss quadrature rule is commonly used to approximate the previous integrals For the shifted Jacobi-Gauss: The weights are given by where C (γ ,η) and the nodes x [42]) The q times differentiation of the Jacobi polynomials J (γ ,η) L,j (x) are given by where where (·) i is the Pochhammer symbol, λ = γ + η + 1, and 3 F 2 is the generalized hypergeometric function.
Lemma 2.2. (see [43]) The Caputo fractional derivative of order ν ∈ (0, 1) of the shifted Jacobi polynomials is given by . (14) Lemma 2.3. Let x (0,0) G,L,r and ̟ (0,0) G,L,r (0 ≤ r ≤ N) be the set of nodes and the weights of the Legendre-Gauss quadrature formula given in (10). Then, the distributed-order fractional derivative of the shifted Jacobi polynomials is given by

ONE-DIMENSIONAL CASE
We consider the following one-dimensional distributed-order fractional Rayleigh-Stokes problem: Let P N ( ) and P M (I) denote the set of polynomials of degree N in space and M in time, respectively. Since we have the homogeneous initial and boundary conditions, we choose appropriate basis for the space ansatz from as well as for time For the sake of clarity, we define the multiindex L = (N, M) and To simplify the notation, we introduce the following integral operators for the Jacobi-Galerkin spectral formulation: The spectral-Galerkin approximation of the solution U ∈ W L is given by where θ (γ ,η) T ,j (t) ∈ P t M (I), and U ij are the unknown coefficients. The key idea behind the Galerkin approximation is to finedÛ ∈ W L such that The actual linear system for (23) depends on the choice of basis functions of W L . We shall construct below suitable spectral basis functions θ (γ ,η) L,i (x) and ϑ (γ ,η) T ,j (t) of W L . Therefore, we construct basis functions using compact combinations of the Jacobi polynomials. In this case, we define θ (γ ,η) where the parameters µ j , κ i and λ i are chosen to satisfy the boundary conditions in (17). Such basis functions are referred to as modal basis functions. (11). Lemma 3.1. (see [1]) For all i, j ≥ 0, there exists a unique set of verify the boundary conditions in (17).
The set of basis functions θ (γ ,η) L,i (x) ∈ P s N+2 ( ) and ϑ (γ ,η) T ,j (t) ∈ P t M+1 (I) are linearly independent. Hence, by dimension argument, we obtain It is clear that the Galerkin formulation of (23) is equivalent to following discrete discretization θ (γ ,η) for 0 ≤ j ≤ M − 1 and 0 ≤ i ≤ N − 2. Throughout this paper, we assume that the indices j and s vary between 0 and M − 1 and that the indices i and r vary between 0 and N − 2. Moreover, we assume that repeated indices are summed over. Thus, the matrix form of (27) becomes Frontiers in Physics | www.frontiersin.org Let us denote where H is the matrix whose entries are H ij = θ (γ ,η) T ,j (t) , U is the matrix of unknown coefficients, and the non-zero elements of the matrices A, B, C, D ω and E are given explicitly by Theorem 3.2. The previous integrals can be computed using the Jacobi-Gauss quadrature rule (9). The discretization of the distributed-order fractional Rayleigh-Stokes problem (16) is equivalent to the following matrix equation In order to be able to solve (29), we shall recast it in a more convenient form. To do so, we make use of the Kronecker product (represented by ⊗). If we consider the matrices H ∈ R n,m and G ∈ R q,p , then the Kronecker product of H and G is defined as follows Let h i ∈ R n denote the columns of H ∈ R n,m so that H = The Kronecker product has the useful property that for any three matrices H, G and H for which the matrix product is defined, we have: where T denotes the transpose. Equation (29) can finally be expressed in matrix form as follows:

Theorem 3.2. Let
Then the non-zero elements A ir , B ir , C js , D ω js , and E js are given by where and h (γ ,η) L,i is given by (8).
Proof: We shall only determine the non-zero elements of A as the proof for the non-zero entries of the other matrices can easily be obtained similarly. From (25), we have Using the orthogonality relation (7), we obtain Frontiers in Physics | www.frontiersin.org In particular, the special cases for the shifted Legendre basis can be obtained directly by taking γ = η = 0, and for the shifted Chebyshev basis of the first and second kinds can be obtained directly by taking γ = η = − 1 2 and γ = η = 1 2 , respectively. Using a proper transformation the problems with nonhomogeneous initial-boundary conditions can be transformed into problems with homogeneous initial-boundary conditions. Let whereŨ is an unknown function satisfying the modified problem with the homogeneous initial and boundary conditions where while U e (x, t) is an arbitrary function satisfying the original non-homogeneous boundary conditions.

TWO-DIMENSIONAL CASE
In this section, we consider the following two-dimensional distributed-order fractional Rayleigh-Stokes problem: with homogeneous initial and boundary conditions, where 2 = × . The two-dimensional Galerkin approximation can be written as T ,j (t).
(46) Then the spectral Jacobi-Galerkin scheme (28) in the two-dimensional case can be expressed in the following form θ (γ ,η) Let us denote where H is the matrix whose entries are T ,j (t) and U is the matrix of unknown coefficients. The Jacobi-Galerkin discretization (47) is equivalent to the following matrix equation For computational convenience, we recast Equation (50) using the Kronecker product in the following matrix form Using a suitable iterative method the above linear system can be solved to obtain the numerical solution (46). In our implementation, the Mathematica function FindRoot with zero initial approximation has been used to solve this system.

NUMERICAL RESULTS AND COMPARISONS
In this section, we present numerical results to verify the efficiency of the spectral Galerkin algorithms. We consider The initial condition, the boundary conditions and the function g(x, t) are selected such as the continuous problem has an exact non-smooth solution in time direction U(x, t) = e x t κ+2 .
In Tables 1, 2, we compare the L ∞ -errors of the present method with the compact finite difference approximation (CFDA) [44], implicit numerical approximation scheme (INAS) [8] and reproducing kernel method (RKM) [9]. We see in these tables that the results are accurate for even small choices of N and M. These results are in perfect agreement with what was expected for a spectral method. Also, this result indicates that the Jacobi-Galerkin method can converge reasonably well for problem (52) with non-smooth data. In Table 3, we list 1the L ∞ -errors of the present method for case III.
Example 5.2. Consider the following two-dimensional problem: ∂U(x, y, t) ∂t = 1 + 0 D w(ν) t ∂ 2 U(x, y, t) ∂x 2 + ∂ 2 U(x, y, t) ∂y 2 + H(x, y, t), The initial condition, the boundary conditions and the right side function H are selected such as the continuous problem has an exact non-smooth solution in time direction U(x, y, t) = e x+y t α+1 .

CONCLUSION
We have presented a Galerkin technique for solving the distributed-order time-fractional Rayleigh-Stokes problem for a generalized second-grade fluid with Jacobi polynomials that is efficient, adaptable to different operators, and easily generalizes to multiple dimensions. By expanding the model solution in terms of Jacobi polynomials in both time and space, we were able to derive adaptable schemes those easily accommodate the distributed fractional-order differential operator. All the calculations can be performed numerically with reasonable accuracy and with relatively small number of degrees of freedom. It should be pointed out that the proposed method to discretize the model equation could also accommodate other numerical methods. For instance, if the model solution is not smooth in time, the order of convergence of the spectral Galerkin schemes may be deteriorated. That could be prevented by simply replacing the Jacobi basis functions in time by fractional order Jacobi functions or using smoothing transformations and deriving the corresponding mass and diffusion matrices by following the same procedure as we described in this study.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/supplementary material.