An efficient two-level factored method for advection-dispersion problem with spatio-temporal coefficients and source terms

1.   Introduction

The study of various phenomena in science and engineering often leads to the formulation of systems of fractional integro-differential equations (FIDEs). Examples of such systems can be found in electric circuit analysis, the activity of interacting inhibitory and excitatory neurons ^[1], glass-forming processes, non-hydrodynamics, drop-wise condensation, and wind ripple formation in deserts ^[2,3]. Due to the presence of fractional derivative operators, obtaining analytical solutions for these functional equations is generally infeasible, like the fractional derivative model of viscoelasticity, the so-called fractional Zener model ^[4,5]. Consequently, researchers have refined existing methods to provide semi-analytical or numerical solutions for FIDEs. For instance, Wang et al. proposed the use of Bernoulli wavelets and operational matrices to solve coupled systems of nonlinear fractional-order integro-differential equations ^[6]. Dief and Grace developed a new technique based on iterative refinement to approximate the analytical solution of a linear system of FIDEs ^[7]. In ^[3], the single term Walsh series (STWS) method was employed to handle second-order Volterra integro-differential equations. The Muntz-Legendre wavelets were applied to find approximate solutions for systems of fractional integro-differential Volterra-Fredholm equations ^[8]. Heydari et al. presented a Chebyshev wavelet method for solving a class of nonlinear singular fractional Volterra integro-differential equations ^[9]. Mohammed and Malik applied a modified series algorithm to solve systems of linear FIDEs ^[10]. In ^[11], Genocchi polynomials combined with the collocation method were utilized to numerically solve a system of Volterra integro-differential equations. Youbi et al. introduced an iterative reproducing kernel algorithm to investigate approximate solutions for fractional systems of Volterra integro-differential equations in the Caputo-Fabrizio operator sense ^[12]. Akbar et al. extended the optimal homotopy asymptotic method to systems of fractional order integro-differential equations ^[13]. Wang et al. combined a mixed element method with the second-order backward difference scheme to numerically solve a class of two-dimensional nonlinear fourth-order partial differential equations ^[14].

Spectral methods offer semi-analytic approximate solutions to various functional equations by expressing them as linear combinations of basis functions. The three main spectral methods are collocation, Galerkin, and tau methods. Orthogonal polynomials serve as fundamental basis functions in spectral methods. Classic polynomials such as Jacobi, Hermite, and Laguerre polynomials are commonly employed for numerical solutions of diverse equations. For instance, Legendre polynomials and Chebyshev polynomials of the first to sixth kind, which are special cases of Jacobi polynomials, are applied in spectral methods for solving multidimensional partial Volterra integro-differential equations, nonlinear fractional delay systems, distributed-order fractional differential equations, fractional-order wave equations, population balance equations, fractional-order diffusion equations, and Volterra-Fredholm integral equations ^{[15,16,17,18,19,20,21]}.

While Gegenbauer polynomials, a specific form of Jacobi polynomials, are less commonly used compared to their orthogonal counterparts, they have found applications in various areas. Usman et al. utilized Gegenbauer polynomials to approximate solutions for multidimensional fractional-order delay problems arising in mathematical physics and engineering ^[22]. Shifted Gegenbauer polynomials have been employed in extracting features from color images ^[23]. Alkhalissi et al. introduced the generalized Gegenbauer-Humbert polynomials for solving fractional differential equations ^[24]. Faheem and Khan proposed a wavelet collocation method based on Gegenbauer polynomials for solving fourth-order time-fractional integro-differential equations with a weakly singular kernel ^[25]. In ^[26], a Gegenbauer wavelets method was utilized to solve FIDEs.

In this study, we focus on solving a system of Caputo fractional-order Volterra integro-differential equations with variable coefficients in the following form:

where $x \in \mathbf{I} = [0, 1]$ and $r = 1, 2$ and the initial conditions are

In (1.1), $\xi _{r, k} \in C(\mathbb{R}), \xi _{r, 0}(x) = 1, r = 1, 2, k = 0, 1, \ldots, n$, $\mu_{r, k}, \gamma _{r, j} \in (0, 1]$, such that $\mu_{r, n} > \mu_{t, n-1} > \ldots > \mu _{r, 1} > 0, \gamma_{r, 1} > \gamma_{r, 0} > 0, \; M_r = \max\{ p_{r, j}, q_{r, k}\}, q_{r, k}-1 < \mu _{r, k}\leq q_{r, k}$, $p_{r, j}-1 < \gamma _{r, j}\leq p_{r, j}, r = 1, 2, j = 1, 2, k = 1, 2, \ldots, n$, $y_r(x), r = 1, 2$ are real continuous functions, $\mathsf{f}_r \in C[0, 1]$, $\kappa _{r, j} \in C([0, 1] \times [0, 1])$ are known functions, $\nu _{r, j} \in \mathbb{R}, r, j = 1, 2$, and ${_0^C}\mathsf{D}_x^{\mu}$ is the Caputo fractional derivative operator. A rigid plate submerged in Newtonian fluid can be modeled using FIDEs like (1.1) ^[27]. Among the important special cases of (1.1), the Bagley-Torvik equation with fractional-order derivative describes the motion of physical systems in Newtonian fluids ^[28]. Structural dynamics most frequently use the Kelvin-Voigt model, which is another special case ^[29].

The use of Gegenbauer polynomials in addressing various types of fractional equations, particularly FIDEs, has been relatively limited (readers can refer to ^{[22,23,24,25,26]}). In this study, our objective is to develop a combined approach using the tau method and Gegenbauer polynomials. Initially, we investigate the existence and uniqueness of solutions to these equations by leveraging Krasnoselskii's fixed point theorem. Subsequently, we derive integral operational matrices of both integer and fractional orders, as well as the operational matrix of the product, specifically tailored to Gegenbauer polynomials. Through appropriate approximations utilizing these operational matrices, the original system is transformed into an algebraic system. By employing the tau spectral method and the inner product by basis functions, we eliminate the independent variable $x$ and obtain a system of $2(N+1)$ algebraic equations that determine the coefficients of the series solutions. Solving this system allows us to approximate the solutions of the main system. Additionally, we estimate an error bound for the residual function within a Gegenbauer-weighted Sobolev space, which reveals that increasing the number of basis functions in the series solutions leads to smaller errors.

The objectives of this paper can be summarized as follows:

● Investigating the existence and uniqueness of solutions for the system of FIDEs with variable coefficients.

● Expressing the approximate solutions of the model (1.1) as linear combinations of shifted Gegenbauer polynomials.

● Developing operational matrices for integration and product operations associated with Gegenbauer polynomials.

● Estimating error bounds for the approximate solutions and residual functions within a Gegenbauer-weighted Sobolev space.

The structure of this paper is as follows: Section 2 provides a brief review of relevant definitions in fractional calculus. The existence and uniqueness of solutions to the system (1.1) are investigated in Section 3. In Section 4, the shifted Gegenbauer polynomials are introduced, their connection to Jacobi polynomials is stated, and operational matrices for integration and product operations are derived. Section 5 outlines the necessary approximations for the functions involved in the system (1.1). Error bounds for the approximate solutions and residual functions are estimated in Section 6. The effectiveness of the proposed method is demonstrated through two numerical examples in Section 7. Finally, Section 8 provides concluding remarks.

2.   Preliminaries

In this section, some definitions and properties of the fractional derivative and integral operators are recalled.

Definition 2.1. The well-known non-integer derivative operator in the Caputo sense of the differentiable function $g$ of the order $\mu \in (0, 1)$ is defined below ^[30]:

The following properties are achieved directly from Definition 2.1:

1) ${_0^C}\mathsf{D}_x^{\mu} \varrho = 0, \quad \varrho \in \mathbb{R};$

2) ${_0^C}\mathsf{D}_x^{\mu} (\lambda _1 \, g_1(x)+\lambda _2 \, g_2(x)) = \lambda _1 \, {_0^C}\mathsf{D}_x^{\mu} g_1(x)+ \lambda _2 \, {_0^C}\mathsf{D}_x^{\mu} g_2(x), \quad \lambda _1, \lambda _2 \in \mathbb{R};$

3) ${_0^C}\mathsf{D}_x^{\mu} x^{\nu} = \left\{\begin{array}{ll} \frac{\Gamma(\nu +1)}{\Gamma(\nu -\mu +1)} x^{\nu -\mu}, & \lfloor \mu \rfloor > \nu, \\ 0, & \text{otherwise} \end{array} \right.$

Definition 2.2. The Riemann-Liouville integral operator of the function $g \in C[0, 1]$ of the order $\mu$ is defined below ^[30]:

The following properties follow from Definition 2.2:

1) ${_0^{\mathsf{RL}}}\mathsf{I}_x^{\mu _1} ({_0^{\mathsf{RL}}}\mathsf{I}_x^{\mu _2}g(x)) = {_0^{\mathsf{RL}}}\mathsf{I}_x^{\mu _1+\mu _2}g(x), \quad \mu _1, \mu _2 \in \mathbb{R};$

2) ${_0^{\mathsf{RL}}}\mathsf{I}_x^{\mu}g(x) = {_0^{\mathsf{RL}}}\mathsf{I}_x^{s-\mu}(D^sg(x)), \quad s = \lceil \mu \rceil, \, \, D^s = \frac{d^s}{dx^s};$

3) ${_0^{\mathsf{RL}}}\mathsf{I}_x^{\mu}x^{\nu} = \left\{\begin{array}{ll} \frac{\Gamma(\nu +1)}{\Gamma(\nu +\mu +1)} x^{\nu +\mu}, & \nu, \mu \in \mathbb{R}^+, \\ 0, & \text{otherwise}; \end{array} \right.$

4) ${_0^{\mathsf{RL}}}\mathsf{I}_x^{\mu}({_0^C}\mathsf{D}_x^{\mu}g(x)) = g(x)-g(0)$

3.   Existence and uniqueness

This section deals with the existence of unique solutions to system (1.1). By applying the Riemann-Liouville integral operator of order $\mu _{r, n} \in (0, 1)$ to Eq (1.1), the following equivalent equation is achieved:

Now, suppose that $C(\mathbf{I}, X)$ is a Banach space of real-valued continuous functions from $\mathbf{I} = [0, 1]$ into $X \subseteq \mathbb{R}$ equipped with the following norm

and the following assumptions hold for any $x \in \mathbf{I}$ and $(x, t) \in \mathbf{I} \times \mathbf{I}$:

Theorem 3.1. Suppose that the assumptions in (3.2) and the following inequality hold

then, problems (1.1)–(1.2) have a unique solution on $C(\mathbf{I}, X)$.

Proof. Let $D_q = \{ z \in C(\mathbf{I}, X)\, |\, \Vert z \Vert _C \leq q\}$ subject to

Then, $D_q$ is a closed, bounded, and convex subset of $C(\mathbf{I}, X)$. The operators $\mathcal{B}_1$ and $\mathcal{B}_2$ are defined as shown below:

It's necessary to be shown that $\mathcal{B}_1+\mathcal{B}_2$ has a fixed point in $D_q$. The process of the proof is divided into four stages:

Stage 1. It is shown that $\mathcal{B}_1 y_r(x)+\mathcal{B}_2 u_r(x) \in D_q$ for every $y_r, u_r \in D_q$. Using (3.2) and (3.4), one obtains the following:

Therefore, $\Vert \mathcal{B}_1 y_r+\mathcal{B}_2 u_r \Vert _C \leq q$, which implies that $\mathcal{B}_1 y_r+\mathcal{B}_2 u_r \in D_q$ for any $y_r, u_r \in D_q$.

Stage 2. It is shown that the operator $\mathcal{B}_1$ is a contraction mapping on $D_q$. For each $y_r, u_r \in D_q$ and each $x \in \mathbf{I}$, one gets the following

From (3.3), this shows that $\mathcal{B}_1$ is a contraction mapping on $D_q$.

Stage 3. It's shown that the operator $\mathcal{B}_2$ is compact and continuous. For $y_r \in D_q$, one has

This shows that $\mathcal{B}_2$ is uniformly bounded on $D_q$. It remains to prove the compactness of the operator $\mathcal{B}_2$. For $x_1, x_2 \in \mathbf{I}$ such that $x_1 < x_2$ and $y_r \in D_q$, one has

By taking norm, one gets

Using a change of variable, one gets:

As $x_2$ tends to $x_1$, the righthand side of the above inequality tends to zero and $\mathcal{B}_2$ is equicontinuous. From Stages 1–3 along with the Arzela-Ascoli theorem, one deduces that the operator $\mathcal{B}_2$ is compact and continuous ^[31]. Finally, by Krasnoselskii's fixed-point theorem ^[32], problem (1.1) has at least a set of solutions as $(y_1(x), y_2(x))$ on $\mathbf{I}$.

Stage 4. Define $\mathcal{E}y_r(x) = \mathcal{B}_1 y_r(x)+\mathcal{B}_2 y_r(x)$. For any $y_r, u_r \in C(\mathbf{I}, X)$ and $x \in \mathbf{I}$ one has:

Imposing the hypothesis in (3.3) implies that $\mathcal{E}$ is a contraction mapping. It follows that $\mathcal{E}$ has a unique fixed point, which is a solution of problems (1.1)–(1.2). □

4.   Shifted Gegenbauer polynomials and their operational matrices

4.1.   Shifted Gegenbauer polynomials

The shifted Gegenbauer polynomials (SGPs) $\mathcal{G}_i^{\sigma}(x), i = 0, 1, 2, \cdots$ are orthogonal polynomials on the interval $[0, 1]$

textcolorred, with respect to the weight function $\omega ^{\sigma}(x) = x^{\sigma -\frac{1}{2}} (1-x)^{\sigma -\frac{1}{2}}, \sigma > -\frac{1}{2}$, that can be defined by the following recurrence relation:

where $\mathbf{I} = [0, 1]$. If $\mathcal{G}_i^{\sigma}(x)$ and $\mathcal{G}_j^{\sigma}(x)$ are the Gegenbauer polynomials of degrees $i$ and $j$, respectively, then they satisfy the following relation:

where $\mathsf{h}_i^{\sigma}$ is the normalizing factor as follows:

These polynomials can be represented in series form as

where the coefficients $\rho _{k, i}^{\sigma}, k = 0, 1, \cdots, i$ are computed below

Remark 4.1. The shifted Gegenbauer polynomials are classified as an especial case of the classic Jacobi polynomials. The relation between these polynomials is as ^[34]:

where $\mathcal{J}_i^{\alpha, \beta}(x)$ is the shifted Jacobi polynomials of the degree $i$. The normalizing factor of the shifted Jacobi polynomials regarding the weight function $w^{\alpha, \beta}(x) = x^{\beta}(1-x)^{\alpha}$ is

and the $l$-th derivative of $\mathcal{J}_i^{\alpha, \beta}(x)$ w.r.t. $x$ is

The square-integrable function $y\in L^2(\mathbb{I})$ can be considered as a linear combination of the SGPs, that is

where the coefficients $\mathsf{y}_k, k = 0, 1, \cdots$ are calculated by the following relation:

To use the shifted Gegenbauer polynomials as basis functions and present approximations in matrix form, a finite form of the series in (4.8) is considered:

where $Y$ and $\Lambda (x)$ are $(N+1) \times 1$ vectors as follows

4.2.   Operational matrices

To deduce computational time and avoid multiplying or integrating the basis functions, introducing and using operational matrices are recommended.

4.2.1.   Operational matrix of the integration with the integer order

Consider the vector $\Lambda ^{\sigma} (x)$ in (4.11); the integral of this vector can be represented in a matrix form. To compute the integral operational matrix, the $i$-th component of $\Lambda ^{\sigma} (x)$ is first considered. Using the series given by (4.3), the integral of $\mathcal{G}_i^{\sigma}(x)$ can be computed as follows:

Now, $x^{k+1}$ is approximated in terms of the Gegenbauer polynomials:

where

Therefore, the integral in (4.12) will be as

Thus, (4.13) can be written in matrix form as

where $\mathbf{\Theta}^{(\sigma)}$ is called the operational matrix of the integration and its entries are calculated below:

4.2.2.   Operational matrix of the integration with the fractional order

Similar to what was stated in the previous subsection, the Riemann-Liouville integral of the $i$-th Gegenbauer polynomial, $\mathcal{G}_i^{\sigma}(x)$, can be calculated as follows:

The function $x^{k+\mu}$ can be approximated as follows

where

So, (4.15) is written as

Therefore, one gets

where $\mathbf{\Theta}^{(\mu, \sigma)}$ is called the operational matrix of the integration of the order $\sigma$ and its entries are calculated as

4.2.3.   Operational matrix of the product

After substituting appropriate approximations into Eq (1.1), terms like $\int _0^x \Lambda ^{\sigma}(t)\, {\Lambda ^{\sigma}}^T(t)\, V\, dt$ are appeared where $V$ is an $(N\times 1)$ vector. To reduce the computational time of the product of the vectors $\Lambda ^{\sigma}(x)$ and ${\Lambda ^{\sigma}}^T(x)$, an algorithm is suggested and this product is accomplished in a matrix form:

where $V^*$ is the $(N\times N)$ operational matrix of the product corresponding to the vector $V$. As an approximate way of calculating its entries, the following steps can be proceeded:

Step 1. If $\mathcal{G}_j^{\sigma}(x)$ and $\mathcal{G}_k^{\sigma}(x)$ are the Gegenbauer polynomials of degrees $j$ and $k$, respectively, compute their products as follows:

where coefficients $\pi _{r}^{(j, k)}, r = 0, 1, \cdots, j+k$ are calculated as:

  if $j\geq k$ then

    for $r = 0..k+j$ do

      if $r > j$ then

        $\pi^{(j, k)}_{r} = \sum\limits _{l = r-j}^{k}\rho_{r-l, j}^{\sigma}\, \rho_{l, k}^{\sigma}$

      else

        $r^* = \min\{r, k\}$

        $\pi^{(j, k)}_{r} = \sum\limits _{l = 0}^{r^*}\rho_{r-l, j}^{\sigma}\, \rho_{l, k}^{\sigma}$

      end if

    end for

  else

    for $r = 0..k+j$ do

      if $r\leq j$ then

        $r^* = \min\{r, j\}$

        $\pi^{(j, k)}_{r} = \sum\limits _{l = 0}^{r^*}\rho_{r-l, j}^{\sigma}\, \rho_{l, k}^{\sigma}$

      else

        $r^{**} = \min\{r, k\}$

        $\pi^{(j, k)}_{r} = \sum\limits _{l = r-j}^{r^{**}}\rho_{r-l, j}^{\sigma}\, \rho_{l, k}^{\sigma}$

      end if

    end for

  end if

Step 2. Using (4.18), compute the integral of the product of three Gegenbauer polynomials $\mathcal{G}_i^{\sigma}(x)\, \mathcal{G}_j^{\sigma}(x),$ and $\mathcal{G}_k^{\sigma}(x)$ as follows

Step 3. The entries of $V^*$ are calculated below:

For more details, please refer to ^[33].

Remark 4.2. In the process of approximating the integral parts in system (1.1), the following integral may be appeared:

where $K$ is an $(N \times N)$ known matrix, $W^*$ is the $(N \times N)$ product operational matrix corresponding to the given $W$, and $\mathbf{\Theta}^{(\sigma)}$ is the integral operational matrix introduced in (4.14). The above integral can be approximated as follows:

where $\Delta$ is an $(N \times 1)$ vector with the following components:

Remark 4.3. If $\Lambda(x)$ is the basis vector in (4.11), one has

where $\mathcal{Q}^{\sigma}$ is the $(N+1) \times (N+1)$ diagonal matrix such that $\mathcal{Q}^{\sigma} _{i, i} = \mathsf{h}_i^{\sigma}, \, \, i = 0, 1, \cdots, N$.

5.   Solution process

In this section, we proceed to approximate the functions and integral components of Eq (1.1) by utilizing the operational matrices derived in Section 4. To this end, we apply the proposed methodology to two distinct systems of fractional Volterra integro-differential equations featuring variable coefficients.

5.1.   System I

Consider the following system of fractional Volterra integro-differential equations with variable coefficients:

with the initial conditions $y_0(0) = 0, y_1(0) = 1$. According to the highest orders of derivatives in Eq (5.1), the following approximations are considered:

where $\mathcal{Y}_0$ and $\mathcal{Y}_1$ are the vectors of unknown coefficients and must be determined. Integrating Eq (5.2) along with the initial conditions leads to the following approximations:

For approximating other terms in (5.1), approximations in Eq (5.2) are written as follows:

so, one has

Using the representation (5.4), one can get the following expressions:

Similarly, one has

Now, using (4.12) and (5.3) yields the following approximations for System I coefficients:

where $\mathcal{W}_1^*$ and $\mathcal{Z}_2^*$ are $(N+1) \times (N+1)$ operational matrices of the product corresponding to the vectors $\mathcal{W}_1$ and $\mathcal{Z}_2$, respectively. The kernels and integral parts can be approximated as follows:

where $\mathbf{\Theta }^{(\sigma)}$ is the operational matrix of the integration in (4.14) and $\mathcal{W}_3^*$ and $\mathcal{W}_4^*$ are the operational matrices of the product. To use the tau method, the sources functions $\mathsf{f}_1$ and $\mathsf{f}_2$ should be approximated in terms of the Gegenbauer polynomials:

where vectors $\mathcal{F}_3$ and $\mathcal{F}_4$ are obtained using Eq (4.9). By substituting approximations (5.2)–(5.11) into system (5.1), the following algebraic system is achieved:

The inner product of algebraic system (5.12) by the basis vector $\Lambda ^{\sigma}(x)$ leads to the following system of algebraic equations:

where $\Delta _i, i = 1, 2$, and $\mathcal{Q}$ are $(N+1)\times (N+1)$ matrices in (4.21) and (4.23), respectively.

5.2.   System II

Now, consider the following system of fractional Volterra integro-differential equations with variable coefficients as the second illustrative example:

with the initial conditions $y_0(0) = y_1(0) = 0$. Based on the highest orders of derivatives in Eq (5.14), the following approximations are set:

Integrating the approximations in (5.15) along with the initial conditions leads to the following approximations:

For approximating other terms in (5.14), approximations in Eq (5.15) are written as follows:

Using (4.12) and (5.16), the coefficients of System II can be approximated as follows:

where $\mathcal{W}_1^*, \mathcal{Z}_1^*$, and $\mathcal{Z}_2^*$ are $(N+1) \times (N+1)$ operational matrices of the product corresponding to the vectors $\mathcal{W}_1\mathcal{Z}_1$ and $\mathcal{Z}_2$, respectively. Now, the kernels and integral parts can be approximated as follows:

where $\mathbf{\Theta }^{(\sigma)}$ is the operational matrix of the integration in (4.14) and $\mathcal{W}_3^*, \mathcal{Z}_1^*$, and $\mathcal{Y}_0^*$ are the operational matrices of the product corresponding to the vectors $\mathcal{W}_3, \mathcal{Z}_1, \mathcal{Y}_0$. To use the tau method, the sources functions $\mathsf{f}_1$ and $\mathsf{f}_2$ should be approximated in terms of the Gegenbauer polynomials:

where vectors $\mathcal{F}_4$ and $\mathcal{F}_5$ are obtained using Eq (4.9). By substituting approximations (5.15)–(5.20) into system (5.14), the following algebraic system is achieved:

The inner product of algebraic system (5.21) by the basis vector $\Lambda ^{\sigma}(x)$ leads to the following system of algebraic equations:

where $\Delta _i, i = 1, 2, 3$, and $\mathcal{Q}$ are $(N+1)\times (N+1)$ matrices in (4.21) and (4.23).

The systems in Eqs (5.13) and (5.22) involve $2(N+1)$ algebraic equations in $2(N+1)$ unknown variables $\mathsf{y}_{0, i}, \mathsf{y}_{1, i}, i = 0, 1, \cdots, N$. By solving the resulted systems, the values of these variables are estimated and approximate solutions derived from Eqs (5.3) and (5.16).

6.   Error bounds

In this section, we aim to compute error bounds for the obtained approximate solutions within a Gegenbauer-weighted Sobolev space. To begin, we present the definitions of this space.

Definition 6.1. Suppose that $m \in \mathbb{N}, \mathbf{I} = [0, 1]$, then the Gegenbauer-weighted Sobolev space $\mathcal{A}_{\omega ^{\sigma}}^m(\mathbf{I})$ is defined as below:

This space is equipped with the following norm and semi-norm:

where $\Vert \, .\, \Vert _{ \omega_k ^{\sigma}}$ denotes the $L^2_{ \omega_k ^{\sigma}}$-norm and $\omega_k ^{\sigma}(x) = x^{\sigma +k-\frac{1}{2}}(1-x)^{\sigma +k-\frac{1}{2}}$.

Definition 6.2. The following inequality holds for any $s \in \mathbb{R}$ and $z \in \mathcal{A}_{\omega ^{\sigma}}^m(\mathbf{I})$

where $s = [s]+\iota, 0 < \iota < 1$. Inequality (6.1) is known as the Gagliardo-Nirenberg inequality ^[35].

Definition 6.3. If $\langle., . \rangle _{m, \omega ^{\sigma}}$ and $\Vert. \Vert _{m, \omega ^{\sigma}}$ are the inner and the norm in the space $\mathcal{A}_{\omega ^{\sigma}}^m(\mathbf{I})$, respectively, then the following inequality holds for any two functions $\phi, \varphi \in \mathcal{A}_{\omega ^{\sigma}}^m(\mathbf{I})$ ^[36]:

Theorem 6.1. Assume that $y \in \mathcal{A}_{\omega ^{\sigma}}^m(\mathbf{I}), m \in \mathbb{N}, 0\leq \eta \leq m,$ and $y_N(x)$ is the Gegenbauer approximation to $y(x)$. An estimation of the error bound can be obtained as follows:

where $\mathcal{C}_0$ is a positive constant.

Proof. Since $y_N(x)$ is the Gegenbauer approximation to $y(x)$, subtracting the series in (4.13) from (4.11) and differentiating it leads to

According to the relation between the Gegenbauer and Jacobi polynomials, mentioned in Remark 4.1, one can get

So, from (6.4) and (4.6), one gets

Similarly, one can obtain

with the aid of the Stirling formula ^[35], the following inequality is obtained:

where $\mathcal{C}_1$ is a positive constant. Based on the definition of the semi-norm and using (6.5)–(6.7), one gets

then, the following inequality is acquired

Using the Gagliardo–Nirenberg inequality leads to the desired bound for $\eta = [\eta]+\eta _0, 0 < \eta _0 < 1$:

□

Corollary 6.2. Using Theorem 6.1, an error bound for $\frac{dy}{dx}-\frac{dy_N}{dx}$ can be estimated as follows:

Therefore, one has

where $\mathcal{C}_1$ is a positive constant.

Theorem 6.3. If $y\in \mathcal{A}_{\omega ^{\sigma}}^m(\mathbf{I}), m \in \mathbb{N}, 0\leq \eta \leq m, 0 < \mu < 1,$ and $y_N(x)$ is the Gegenbauer approximation to $y(x)$, then an error bound can be estimated for ${_0^C}\mathsf{D}_x^{\mu} y(x)-{_0^C}\mathsf{D}_x^{\mu} y_N(x)$ as follows:

Proof. By noting the relation between the classical derivative and Riemann-Liouville integral operators (the second property in Definition 2.2), one has

where the star sign denotes the convolution of $x^{-\mu}$ and $(y^{\prime}(x) -y^{\prime}_N(x))$. Applying Young's convolution inequality $\Vert u_1 *u_2 \Vert _q \leq \delta _0\, \Vert u_1 \Vert _1 \, \Vert u_2 \Vert _q$ and Corollary 6.2 to (6.11) leads to the desired result:

□

Now, suppose that $y_N(x)$ is the obtained solutions from the suggested method, so one has the following approximate system:

where $\mathcal{R}_r(x), r = 1, 2$ are the residual functions. By subtracting Eq (6.12) from Eq (1.1), the following error system is achieved:

where $e_{j, N} = y_{j}(x)-y_{j, N}(x)$. To obtain an error bound for the residual function in (6.13), this function is multiplied by $e_{r, N}$; therefore, using the inequality in (6.2) yields:

where $\omega_0^{\sigma}(x, t) = \omega ^{\sigma}(x)\, \omega ^{\sigma}(t)$ and $\mathbf{I}^2 = \mathbf{I} \times \mathbf{I}$. So, one has

From the righthand side of (6.14), increasing the value of $N$ leads to a smaller error bound for the residual function.

7.   Numerical examples

In Section 5, we implemented the proposed approach on two systems of Volterra-type FIDEs. Systems I and II were transformed into corresponding systems of algebraic equations, as shown in Eqs (5.13) and (5.22). We solved these systems for various values of $N$ and $\sigma$ and calculated the maximum absolute errors (MAEs) and least square errors (LSEs). The obtained results are presented in the forms of figures and tables. All computations were performed using Maple 16 software. Also, the convergence rate $C_R(y_j), j = 0, 1$ is computed by the following formula:

Example 1. Consider system (5.1) with the exact solutions $(y_0(x), y_1(x)) = (-x^2, 1-x^3)$, initial conditions $(y_0(0), y_1(0)) = (0, 1)$, and source functions

Table 1 displays the values of approximate solutions at equally spaced points $x_i = 0.2 i, i = 0, 1, \cdots, 5$ as well as the corresponding least square errors obtained from the proposed tau-Gegenbauer method (for $N = 7$, $\sigma = 1$) and the block-by-block approach combined with a finite difference approximation ^[37] (for $N = 10$, $h = 0.1$). It is evident that the tau-Gegenbauer method yields results with significantly smaller errors compared to those reported in ^[37]. MAEs for various values of $N$ (ranging from 2 to 7) and $\sigma = 1$ are presented in Table 2. As expected, increasing the number of Gegenbauer polynomials in the series solutions (represented by $N$) leads to a reduction in errors. Table 4 reports the MAEs of the approximate solutions for $N = 5$ and different values of $\sigma$.

In Figure 1, we present the exact and approximate solutions, as well as the corresponding absolute error functions, for $N = 7$ and $\sigma = 1$. The plots demonstrate a good agreement between the approximate solutions and the exact ones. Furthermore, Figure 2 illustrates the absolute error functions for various values of $\sigma$ while keeping $N = 5$ constant. These plots further highlight the accuracy of the proposed method. The results substantiate the effectiveness and reliability of the tau-Gegenbauer method in providing accurate approximate solutions. The convergence rate of approximate solutions have been computed and listed in Table 2 for different values of $N$. In adition, the convergence rates at $x = 0.5, 0.7$ are calculated for different values of $N$ and $\sigma = 1$, which can be seen in Table 3.

Example 2. Consider system (5.14) as the second example with the exact solutions $(y_0(x), y_1(x)) = (x^4, x(1-x))$, initial conditions $(y_0(0), y_1(0)) = (0, 0)$, and source functions

A comparison between exact and approximate solutions at equally spaced points $x_i = 0.2 i, i = 0, 1, \cdots, 5$ is seen in Table 5 for $N = 7$, $\sigma = 1$. The LSEs of the results obtained from the proposed method are smaller than those obtained from the block-by-block approach in ^[37] (for $N = 10$, $h = 0.1$). MAEs have been computed for $N = 2, 3, \cdots, 8$ and $\sigma = 1$ and listed in Table 6. In Table 7, MAEs of the approximate solutions are computed for $N = 5$ and diverse values of $\sigma$. Figures of exact and approximate solutions and absolute error functions are plotted in Figure 3 for $N = 7$ and $\sigma = 1$. Plots of absolute error functions are seen in Figure 4 for various values of $\sigma$ and $N = 5$. The convergence rate of approximate solutions have been computed and listed in Table 6 for different values of $N$. Error bounds of $y_{0, N}(x)$ and $y_{1, N}(x)$ in (6.3) are $8.2270\times 10^{-3}$ and $2.0428\times 10^{-3}$ for $N = 7, m = 2, \eta = 1, \sigma = 1$, respectively. Based on the convergenc rate in Table 6, the absolute error of $y_{0, 7}(x)$ is proportional to $\frac{1}{N^{2.1527}} = 1.5162\times 10^{-2}$ and the absolute error of $y_{1, 7}(x)$ is proportional to $\frac{1}{N^{1.7767}} = 3.1515\times 10^{-2}$.

8.   Conclusions

A matrix-based tau spectral method utilizing Gegenbauer polynomials was developed to numerically solve a specific class of FIDE systems. A fractional-order integral operational matrix was constructed for this purpose. In ^[22], differential equations with time-fractional delays were analyzed using spectral operational matrices. Operational matrices of the fractional-order derivative related to shifted Gegenbauer polynomials were derived. To solve fractional differential equations based on the Gegenbauer-Humbert polynomials, a collocation wavelet method was proposed in ^[24]. The operational matrices of the fractional derivatives were derived. Additional to the main equation, derivative operational matrices were also used to approximate initial and boundary conditions. In contrast, the authors of the current paper used integral operational matrices without requiring conditions to be approximated. The authors have observed in previous works that operational matrices reduce error. The derived error bound for the residual function, evaluated within a Gegenbauer-weighted Sobolev space, demonstrates that selecting a sufficiently large value for the parameter $N$ leads to a sufficiently small error. Notably, the Gegenbauer polynomials were dependent on the parameter $\sigma$, and altering its value yields different versions of these polynomials. The impact of varying $\sigma$ can be observed in Tables 4 and 7, as well as Figures 2 and 4. The most favorable outcomes were obtained for $\sigma = 0.5$ and $1$. A comparison between the results obtained from the proposed method and those reported in ^[37] (employing a block-by-block approach combined with a finite difference method) was presented in Tables 1 and 5. The numerical results obtained through the tau-Gegenbauer scheme exhibited better agreement with the exact solutions. The authors intend to apply the proposed approach to systems of fractional integral equations with variable orders and generalize it to solve two-dimensional integro-partial differential equations of the fractional order.

Use of AI tools declaration

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

Acknowledgments

The authors are thankful to dear referees for the valuable comments which have improved the final version of this work. The authors extend their appreciation to the Deputyship for Research and Innovation, Ministry of Education in Saudi Arabia for funding this research work through project number: 445-9-676. Furthermore, the authors would like to extend their appreciation to Taibah University for its supervision support.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.