Numerical treatment of the fractional Rayleigh-Stokes problem using some orthogonal combinations of Chebyshev polynomials

1.   Introduction

A wide range of fields rely on Chebyshev polynomials (CPs). Some CPs are famously special polynomials of Jacobi polynomials (JPs). We can extract four kinds of CPs from JPs. They were employed in many applications; see ^[1,2,3,4]. However, others can be considered special types of generalized ultraspherical polynomials; see ^[5,6]. Some contributions introduced and utilized other specific kinds of generalized ultraspherical polynomials. In the sequence of papers ^[7,8,9,10], the authors utilized CPs of the fifth- and sixth-kinds to treat different types of differential equations (DEs). Furthermore, the eighth-kind CPs were utilized in ^[11,12] to solve other types of DEs.

Several phenomena that arise in different applied sciences can be better understood by delving into fractional calculus, which studies the integration and derivatives for non-integer orders. When describing important phenomena, fractional differential equations (FDEs) are vital. There are many examples of FDEs applications; see, for instance, ^[13,14,15]. Because it is usually not feasible to find analytical solutions for these equations, numerical methods are relied upon. Several methods were utilized to tackle various types of FDEs. Here are some techniques used to treat several FDEs: the Adomian decomposition method ^[16], a finite difference scheme ^[17], generalized finite difference method ^[18], Gauss collocation method ^[19], the inverse Laplace transform ^[20], the residual power series method ^[21], multi-step methods ^[22], Haar wavelet in ^[23], matrix methods in ^[24,25,26], collocation methods in ^{[27,28,29,30]}, Galerkin methods in ^[31,32,33], and neural networks method in ^[34].

Among the essential FDEs are the Rayleigh-Stokes equations. The fractional Rayleigh-Stokes equation is a mathematical model for the motion of fluids with fractional derivatives. This equation is used in many areas of study, such as non-Newtonian fluids, viscoelastic fluids, and fluid dynamics. Many contributions were devoted to investigating the types of Rayleigh-Stokes from a theoretical and numerical perspective. Theoretically, one can consult ^[35,36,37]. Several numerical approaches were followed to solve these equations. In ^[38], the authors used a finite difference method for the fractional Rayleigh-Stokes equation (FRSE). In ^[39], a computational method for two-dimensional FRSE is followed. The authors of ^[40] used a finite volume element algorithm to treat a nonlinear FRSE. In ^[41], a numerical method is applied to handle a type of Rayleigh-Stokes problem. Discrete Hahn polynomials treated variable-order two-dimensional FRSE in ^[42]. The authors of ^[43] numerically solved the FRSE.

The significance of spectral approaches in engineering and fluid dynamics has been better understood in recent years, and this trend is being further explored in the applied sciences ^[44,45,46]. In these techniques, approximations to integral and differential equations are assumed by expanding a variety of polynomials, which are frequently orthogonal. The three spectral techniques used most often are the collocation, tau, and Galerkin methods. The optimal spectral approach to solving the provided equation depends on the nature of the DE and the governing conditions that regulate it. The three spectral methods use distinct trial and test functions. In the Galerkin method, the test and trial functions are chosen so that each basis function member meets the given DE's underlying constraints; see ^[47,48]. The tau method is not limited to a specific set of basis functions like the Galerkin approach. This is why it solves many types of DEs; see ^[49]. Among the spectral methods, the collocation method is the most suitable; see, for example, ^[50,51].

In his seminal papers ^[52,53], Shen explored a new idea to apply the Galerkin method. He selected orthogonal combinations of Legendre and first-kind CPs to solve the second- and fourth-order DEs. The Galerkin approach was used to discretize the problems with their governing conditions. To address the even-order DEs, the authors of ^[54] employed a generalizing combination to solve even-order DEs.

This paper's main contribution and significance is the development of a new Galerkin approach for treating the FRSE. The suggested technique has the advantage that it yields accurate approximations by picking a small number of the retained modes of the selected Galerkin basis functions.

The current paper has the following structure. Section 2 presents some preliminaries and essential relations. Section 3 describes a Galerkin approach for treating FRSE. A comprehensive study on the convergence analysis is studied in Section 4. Section 5 is devoted to presenting some illustrative examples to show the efficiency and applicability of our proposed method. Section 6 reports some conclusions.

2.   Some fundamentals and formulas

This section defines the fractional Caputo derivative and reviews some of its essential properties. Next, we gather significant characteristics of the second-kind CPs. This paper will use some orthogonal combinations of the second-kind CPs to solve the FRSE.

2.1.   Caputo's fractional derivative

Definition 2.1. In Caputo's sense, the fractional-order derivative of the function $\xi(s)$ is defined as ^[55]

For $D^{\alpha}$ with $p-1 < \alpha < p, \; p\in \mathbb{N}$, the following identities are valid:

where ${\mathbb{N}} = \{1, 2, ...\}$ and ${\mathbb{N}}_0 = \{0, 1, 2, \ldots\}$, and $\lceil \alpha\rceil$ is the ceiling function.

2.2.   Shifted second-kind CPs

The shifted second-kind CPs $\mathbf{U}^{*}_{j}(t)$ are orthogonal regarding the weight function $\omega(t) = \sqrt{t\, (\tau-t)}$ in the interval $[0, \tau]$ and defined as ^[56,57]

where

with the following orthogonality relation ^[56]:

where

$\{\mathbf{U}^{*}_{m}(t)\}_{m\ge 0}$ can be generated by the recursive formula:

The following theorem that presents the derivatives of $\mathbf{U}^{*}_{m}(t)$ is helpful in what follows.

Theorem 2.1. ^[56] For all $j\ge n$, the following formula is valid:

The following particular formulas of (2.9) give expressions for the first- and second-order derivatives.

Corollary 2.1. The following derivative formulas are valid:

Proof. Special cases of Theorem 2.1. □

3.   Treatment for the FRSE

This section is devoted to analyzing a Galerkin approach to solve the following FRSE ^[38,58]:

governed by the following constraints:

where $a$ and $b$ are two positive constants and $\mathcal{S}(x, t)$ is a known smooth function.

Remark 3.1. The well-posedness and regularity of the fractional Rayleigh-Stokes problem are discussed in detail in ^[36].

3.1.   Selection of trial functions

We choose the trial functions to be

Due to (2.6), it can be seen that $\left\{\varphi_i(x)\right\}_{i\ge 0}$ satisfies the following orthogonality relation:

where

and $\hat{\omega}(x) = \frac{1}{x^{\frac{3}{2}}\, (\ell -x)^{\frac{3}{2}}}.$

Theorem 3.1. The second-derivative of $\varphi_i(x)$ can be expressed explicitly in terms of $\mathbf{U}^{*}_j(x)$ as

where

Proof. Based on the basis functions in (3.4), we can write

Using Corollary 2.1, Eq (3.9) may be rewritten as

With the aid of the recurrence relation (2.8), the following recurrence relation for $\mathbf{U}^{*}_{i}(x)$ holds:

Moreover, the last relation enables us to write the following relation:

If we insert relations (3.11) and (3.12) into relation (3.10), and perform some computations, then we get

where

This completes the proof.□

3.2.   Galerkin solution for FRSE with homogeneous boundary conditions

Consider the FRSE (3.1), governed by the conditions: $v(0, t) = v(\ell, t) = 0$.

Now, consider the following spaces:

where $\Omega = (0, \ell)\times(0, \tau].$

The approximate solution $\hat{v}(x, t)\in \mathbb{X}_{\mathcal{M}}(\Omega)$ may be expressed as

where

and $\boldsymbol{{C}} = (c_{ij})_{0\leq i, j\leq \mathcal{M}}$ is the unknown matrix to be determined whose order is $(\mathcal{M}+1)\times (\mathcal{M}+1)$.

The residual $\mathbf{\mathcal{R}}(x, t)$ of Eq (3.1) may be calculated to give

The philosophy in applying the Galerkin method is to find $\hat{v}(x, t)\in\mathbb{X}_{\mathcal{M}}(\Omega)$, such that

where $\bar{\omega}(x, t) = \hat{\omega}(x)\, \omega(t).$ The last equation may be rewritten as

In matrix form, Eq (3.19) can be written as

where

Moreover, (3.2) implies that

Now, Eq (3.20) along with (3.27) constitutes a system of algebraic equations of order $(\mathcal{M}+1)^{2}$, that may be solved using a suitable numerical procedure.

Now, the derivation of the formulas of the entries of the matrices $\mathit{\boldsymbol{\mathcal{A}}}$, $\mathit{\boldsymbol{\mathcal{B}}}$, $\mathit{\boldsymbol{\mathcal{H}}}$, $\mathit{\boldsymbol{\mathcal{K}}}$ and $\mathit{\boldsymbol{\mathcal{Q}}}$ are given in the following theorem.

Theorem 3.2. The following definite integral formulas are valid:

where $q_{j, s}$ and $a_{i, r}$ are given respectively in Eqs (2.6) and (3.6). Also, we have

and

where $\ _3\tilde{F}_2$ is the regularized hypergeometric function ^[59].

Proof. To find the elements $h_{i, r}$: Using Theorem 3.1, one has

Now, $\int_{0}^{\ell}\hat{\omega}(x)\, \mathbf{U}^{*}_j(x)\, \varphi_r(x)\, dx$ can be calculated to give the following result:

and therefore, we get the following desired result:

To find the elements $b_{j, s}$: Formula (2.10) along with the orthogonality relation (2.6) helps us to write

To find $k_{j, s}$: Using property (2.3) together with (2.4), one can write

If we note the following identity:

then, we get

Theorem 3.2 is now proved. □

Remark 3.2. The following algorithm shows our proposed Galerkin technique, which outlines the necessary steps to get the approximate solutions.

Remark 3.3. Based on the following substitution:

the FRSE (3.1) with non-homogeneous boundary conditions will convert to homogeneous ones $y(0, t) = y(\ell, t) = 0$.

4.   Error bound

In this section, we study the error bound for the two cases corresponding to the 1-D and 2-D Chebyshev-weighted Sobolev spaces.

Assume the following Chebyshev-weighted Sobolev spaces:

where $I_{1} = (0, \tau)$ and $I_{2} = (0, \ell)$ are quipped with the inner product, norm, and semi-norm

where $0 < \alpha < 1$ and $m\in \mathbb{N}.$

Also, assume the following two-dimensional Chebyshev-weighted Sobolev space:

equipped with the norm and semi-norm

where $0 < \alpha < 1$ and $r, s\in \mathbb{N}.$

Lemma 4.1. ^[60] For $n\in \mathbb{N},$ $n+r > 1,$ and $n+s > 1,$ where $r, s \in \mathbb{R}$ are any constants, we have

where

Theorem 4.1. Suppose $0 < \alpha < 1,$ and $\bar{\eta}(t) = \sum\limits_{j = 0}^{\mathcal{M}}\hat{\eta}_{j}\, \mathbf{U}^{*}_{j}(t)$ is the approximate solution of $\eta(t)\in\mathbf{H}^{\alpha, m}_{\omega(t)}(I_{1}).$ Then, for $0\leq k\leq m\leq \mathcal{M}+1,$ we get

where $\mathcal{A}\lesssim \mathcal{B}$ indicates the existence of a constant $\nu$ such that $\mathcal{A}\leq \nu\, \mathcal{B}.$

Proof. The definitions of $\eta(t)$ and $\hat{\eta}(t)$ allow us to have

To estimate the factor $\frac{||D_{t}^{\alpha+k}\, \mathbf{U}^{*}_{\mathcal{M}+1}(t)||^{2}_{L^{2}_{\omega(t)}}}{||D_{t}^{\alpha+m}\, \mathbf{U}^{*}_{\mathcal{M}+1}(t)||^{2}_{L^{2}_{\omega(t)}}},$ we first find $||D_{t}^{\alpha+k}\, \mathbf{U}^{*}_{\mathcal{M}+1}(t)||^{2}_{L^{2}_{\omega(t)}}.$

Equation (2.3) along with (2.4) allows us to write

and accordingly, we have

The following inequality can be obtained after applying the Stirling formula ^[44]:

By virtue of the Stirling formula ^[44] and Lemma 4.1, $||D_{t}^{\alpha+k}\, \mathbf{U}^{*}_{\mathcal{M}+1}(t)||^{2}_{L^{2}_{\omega(t)}}$ can be written as

where $\lambda^* = \max\limits_{0\leq r\leq \mathcal{M}+1}\left\{\frac{\lambda^{2}_{r, \mathcal{M}+1}\, \sqrt{\pi}}{2}\right\}$.

Similarly, we have

and accordingly, we have

Inserting Eq (4.16) into Eq (4.9), one gets

Therefore, we get the desired result. □

Theorem 4.2. Suppose $\bar{\zeta}(x) = \sum_{i = 0}^{\mathcal{M}} \hat{\zeta}_{i}\, \varphi_i(x),$ is the approximate solution of $\zeta_{\mathcal{M}}(x)\in\mathbf{Y}^{m}_{\hat{\omega}(x)}(I_{2}).$ Then, for $0\leq k\leq m\leq \mathcal{M}+1,$ we get

Proof. At first, based on the definitions of $\zeta(x)$ and $\bar{\zeta}(x),$ one has

Now, we have

and therefore, $||D_{x}^{k}\, \varphi_{\mathcal{M}+1}(x)||^{2}_{L^{2}_{\hat{\omega}(x)}}$ can be written as

The application of the Stirling formula ^[44] leads to

and hence, we get

Finally, we get the following estimation:

At the end, we get

□

Theorem 4.3. Given the following assumptions: $\alpha = 0,$ $0\leq p\leq r\leq \mathcal{M}+1,$ and the approximation to $v(x, t)\in{\mathbf{H}^{r, s}_{\ddot{\omega}}}(\Omega)$ is $\hat{v}(x, t)$. As a result, the estimation that follows is applicable:

Proof. According to the definitions of $v(x, t)$ and $\hat{v}(x, t),$ one has

Now, applying the same procedures as in Theorem 4.2, we obtain

□

Theorem 4.4. Given the following assumptions: $\alpha = 0,$ $0\leq q\leq s\leq \mathcal{M}+1,$ and the approximation to $v(x, t)\in{\mathbf{H}^{r, s}_{\ddot{\omega}}}(\Omega)$ is $\hat{v}(x, t)$. As a result, the estimation that follows is applicable:

Theorem 4.5. Let $\hat{v}(x, t)$ be the approximate solution of $v(x, t)\in{\mathbf{H}^{r, s}_{\bar{\omega}(x, t)}}(\Omega)$, and assume that $0 < \alpha < 1$. Consequently, for $0\leq p\leq r\leq \mathcal{M}+1,$ and $0\leq q\leq s\leq \mathcal{M}+1,$ we obtain

Proof. The proofs of Theorems 4.4 and 4.5 are similar to the proof of Theorem 4.3. □

Theorem 4.6. Let $\mathbf{\mathcal{R}}(x, t)$ be the residual of Eq (3.1), then $\left|\left|\mathbf{\mathcal{R}}(x, t)\right|\right|_{L^{2}_{\bar{\omega}(x, t)}}\rightarrow 0$ as $\mathcal{M}\rightarrow \infty$.

Proof. $\left|\left|\mathbf{\mathcal{R}}(x, t)\right|\right|_{L^{2}_{\bar{\omega}(x, t)}}$ of Eq (3.28) can be written as

Now, the application of Theorems 4.3–4.5 leads to

Therefore, it is clear that $\left|\left|\mathbf{\mathcal{R}}(x, t)\right|\right|_{L^{2}_{\bar{\omega}(x, t)}}\rightarrow 0$ as $\mathcal{M}\rightarrow \infty.$ □

5.   Illustrative examples

This section will compare our shifted second-kind Galerkin method (SSKGM) with other methods. Three test problems will be presented in this regard.

Example 5.1. ^[38] Consider the following equation:

where

governed by (3.2) and (3.3). Problem (5.1) has the exact solution: $u(x, t) = x^3\, (\ell -x)^3\, t^2.$

In Table 1, we compare the $L_{2}$ errors of the SSKGM with that obtained in ^[38] at $\ell = \tau = 1$. Table 2 reports the amount of time for which a central processing unit (CPU) was used for obtaining results in Table 1. These tables show the high accuracy of our method. Figure 1 illustrates the absolute errors (AEs) at different values of $\alpha$ at $\mathcal{M} = 4$ when $\ell = \tau = 1$. Figure 2 illustrates the AEs at different values of $\alpha$ at $\mathcal{M} = 4$ when $\ell = 3,$ and $\tau = 2$. Figure 3 shows the AEs at different $\alpha$ at $\mathcal{M} = 4$ when $\ell = 10,$ and $\tau = 5$.

Example 5.2. ^[38] Consider the following equation:

where

governed by (3.2) and (3.3). Problem (5.3) has the exact solution: $u(x, t) = t^2\, \sin (\pi \, x).$

Table 3 compares the $L_{2}$ errors of the SSKGM with those obtained by the method in ^[38] at $\ell = \tau = 1$. This table shows that our results are more accurate. Table 4 reports the CPU time used for obtaining results in Table 3. Moreover, Figure 4 sketches the AEs at different values $\mathcal{M}$ when $\alpha = 0.7$, and $\ell = \tau = 1$. Table 5 presents the maximum AEs at $\alpha = 0.8$ and $\mathcal{M} = 8$ when $\ell = \tau = 1$. Figure 5 sketches the AEs at different $\alpha$ for $\mathcal{M} = 10,$ $\ell = 3$ and $\tau = 1$.

Example 5.3. Consider the following equation:

where

governed by (3.2) and (3.3). The exact solution of this problem is: $u(x, t) = t^4\, \sin (2\, \pi\, x).$

Table 6 presents the maximum AEs at $\alpha = 0.5$ and $\mathcal{M} = 9$ when $\ell = \tau = 1$. Figure 6 sketches the AEs at different $\mathcal{M}$ and $\alpha = 0.9$ when $\ell = \tau = 1$.

Example 5.4. ^[61] Consider the following equation:

governed by the following constraints:

where

and the exact solution of this problem is: $u(x, t) = e^x\, t^{\gamma +2}.$ This problem is solved for the case $\gamma = 1$. In Table 7, we compare the $L_{2}$ errors of the SSKGM with that obtained in ^[61] at different values of $\alpha$. This table shows the high accuracy of our method. Figure 7 illustrates the AE (left) and the approximate solution (right) at $\mathcal{M} = 7$ when $\alpha = 0.6$.

6.   Conclusions

This study presented a Galerkin algorithm technique for solving the FRSE using orthogonal combinations of the second-kind CPs. The Galerkin method converts the FRSE with its underlying conditions into a matrix system whose entries are given explicitly. A suitable algebraic algorithm may be utilized to solve such a system, and by chance, the approximate solution can be obtained. We showcased the effectiveness and precision of the algorithm through a comprehensive study of the error analysis and by presenting multiple numerical examples. We think the proposed method can be applied to other types of FDEs. As an expected future work, we aim to employ this paper's developed theoretical results and suitable spectral methods to treat some other problems.

Authors contributions

W. M. Abd-Elhameed: Conceptualization, Methodology, Validation, Formal analysis, Funding acquisition, Investigation, Project administration, Supervision, Writing–Original draft, Writing–review & editing. A. M. Al-Sady: Methodology, Validation, Writing–Original draft; O. M. Alqubori: Methodology, Validation, Investigation; A. G. Atta: Conceptualization, Methodology, Validation, Formal analysis, Visualization, Software, Writing–Original draft, Writing–review & editing. All authors have read and agreed to the published version of the manuscript.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgement

This work was funded by the University of Jeddah, Jeddah, Saudi Arabia, under grant No. (UJ-23-FR-70). Therefore, the authors thank the University of Jeddah for its technical and financial support.

Conflict of interest

The authors declare that they have no competing interests.