A tau-Gegenbauer spectral approach for systems of fractional integro-differential equations with the error analysis

Khadijeh Sadri; David Amilo; Kamyar Hosseini; Evren Hinçal; Aly R. Seadawy; Khadijeh Sadri; David Amilo; Kamyar Hosseini; Evren Hinçal; Aly R. Seadawy

doi:10.3934/math.2024190

AIMS Mathematics

2024, Volume 9, Issue 2: 3850-3880. doi: 10.3934/math.2024190

Previous Article Next Article

Research article

A tau-Gegenbauer spectral approach for systems of fractional integro-differential equations with the error analysis

1.
Department of Mathematics, Near East University TRNC, Mersin 10, Nicosia 99138, Turkey
2.
Mathematics Research Center, Near East University TRNC, Mersin 10, Nicosia 99138, Turkey
3.
Department of Computer Science and Mathematics, Lebanese American University, Beirut, Lebanon
4.
Mathematics Department, Faculty of Science, Taibah University, Al-Madinah Al-Munawarah 41411, Kingdom of Saudi Arabia

Received: 21 November 2023 Revised: 20 December 2023 Accepted: 27 December 2023 Published: 11 January 2024
MSC : 34K05, 34K28, 65R20

This research paper focused on the solution of systems of fractional integro-differential equations (FIDEs) of the Volterra type with variable coefficients. The proposed approach combined the tau method and shifted Gegenbauer polynomials in a matrix form. The investigation of the existence and uniqueness of solutions for these systems was carried out using Krasnoselskii's fixed point theorem. The equations employed Caputo-style derivative operators, and to minimize computational operations involving derivatives and multiplications, integral and product operational matrices were derived. By introducing suitable polynomial approximations and employing the tau spectral method, the original system of FIDE was transformed into an algebraic system. Solving this algebraic system provided approximate solutions to the main system. Error bounds were computed in the Gegenbauer-weighted Sobolev space. The proposed algorithm was implemented and tested on two systems of integro-fractional differential equations to demonstrate its efficiency and simplicity. By varying the parameter $\sigma$ in the Gegenbauer polynomials, the impact of this variation on the approximate solutions can be observed. A comparison with another method utilizing the block-by-block approach was also presented.

Keywords:

shifted Gegenbauer polynomials,
system of fractional integro-differential equations (FIDEs),
existence and uniqueness,
error bound

Citation: Khadijeh Sadri, David Amilo, Kamyar Hosseini, Evren Hinçal, Aly R. Seadawy. A tau-Gegenbauer spectral approach for systems of fractional integro-differential equations with the error analysis[J]. AIMS Mathematics, 2024, 9(2): 3850-3880. doi: 10.3934/math.2024190

Related Papers:

[1]	Chuanhua Wu, Ziqiang Wang . The spectral collocation method for solving a fractional integro-differential equation. AIMS Mathematics, 2022, 7(6): 9577-9587. doi: 10.3934/math.2022532
[2]	Ali H. Tedjani, Sharifah E. Alhazmi, Samer S. Ezz-Eldien . An operational approach for one- and two-dimension high-order multi-pantograph Volterra integro-differential equation. AIMS Mathematics, 2025, 10(4): 9274-9294. doi: 10.3934/math.2025426
[3]	Obaid Algahtani, M. A. Abdelkawy, António M. Lopes . A pseudo-spectral scheme for variable order fractional stochastic Volterra integro-differential equations. AIMS Mathematics, 2022, 7(8): 15453-15470. doi: 10.3934/math.2022846
[4]	Xiaojun Zhou, Yue Dai . A spectral collocation method for the coupled system of nonlinear fractional differential equations. AIMS Mathematics, 2022, 7(4): 5670-5689. doi: 10.3934/math.2022314
[5]	Mahmoud A. Zaky, Weam G. Alharbi, Marwa M. Alzubaidi, R.T. Matoog . A Legendre tau approach for high-order pantograph Volterra-Fredholm integro-differential equations. AIMS Mathematics, 2025, 10(3): 7067-7085. doi: 10.3934/math.2025322
[6]	Chuanli Wang, Biyun Chen . An $ hp $-version spectral collocation method for fractional Volterra integro-differential equations with weakly singular kernels. AIMS Mathematics, 2023, 8(8): 19816-19841. doi: 10.3934/math.20231010
[7]	Yumei Chen, Jiajie Zhang, Chao Pan . Numerical approximation of a variable-order time fractional advection-reaction-diffusion model via shifted Gegenbauer polynomials. AIMS Mathematics, 2022, 7(8): 15612-15632. doi: 10.3934/math.2022855
[8]	Kanagaraj Muthuselvan, Baskar Sundaravadivoo, Kottakkaran Sooppy Nisar, Suliman Alsaeed . Discussion on iterative process of nonlocal controllability exploration for Hilfer neutral impulsive fractional integro-differential equation. AIMS Mathematics, 2023, 8(7): 16846-16863. doi: 10.3934/math.2023861
[9]	Saowaluck Chasreechai, Sadhasivam Poornima, Panjaiyan Karthikeyann, Kulandhaivel Karthikeyan, Anoop Kumar, Kirti Kaushik, Thanin Sitthiwirattham . A study on the existence results of boundary value problems of fractional relaxation integro-differential equations with impulsive and delay conditions in Banach spaces. AIMS Mathematics, 2024, 9(5): 11468-11485. doi: 10.3934/math.2024563
[10]	Chang Phang, Abdulnasir Isah, Yoke Teng Toh . Poly-Genocchi polynomials and its applications. AIMS Mathematics, 2021, 6(8): 8221-8238. doi: 10.3934/math.2021476

Abstract

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:

$\begin{equation} {_0^C}\mathsf{D}_x^{\mu _{r, n}}y_r(x)+\sum\limits _{k = 1}^{n-1} \xi _{r, k}(x)\, {_0^C}\mathsf{D}_x^{\mu _{r, n-k}}y_r(x)+\xi _{r, n}y_r(x)-\sum \limits _{j = 1}^{2}\nu _{r, j} \int _0^x \kappa _{r, j}(x, t)\, {_0^C}\mathsf{D}_t^{\gamma _{r, 2-j}}\, y_j(t)\, dt -\mathsf{f}_r(x) = 0 \end{equation}$

(1.1)

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

$\begin{equation} \frac{d^{m_r} }{dx}y_r(0) = y_{r, m_{r}}, \quad m_r = 0, 1, \ldots, M_{r-1}, \, r = 1, 2 \end{equation}$

(1.2)

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]:

$\begin{equation} {_0^C}\mathsf{D}_x^{\mu} g(x) = \frac{1}{\Gamma(1-\mu)} \int _0^x (x-s)^{-\mu} g^{\prime}(s) ds, \quad 0 < \mu < 1 \end{equation}$

(2.1)

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]:

$\begin{equation} {_0^{\mathsf{RL}}}\mathsf{I}_x^{\mu} g(x) = \frac{1}{\Gamma(\mu)} \int _0^x (x-s)^{\mu -1} g(s) ds, \quad \mu > 1 \end{equation}$

(2.2)

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:

$\begin{equation} \begin{aligned} y_r(x) = y_{r, 0} &-\sum\limits _{k = 1}^{n-1} \frac{1}{\Gamma(\mu _{r, n})} \int _0^x (x-t)^{\mu _{r, n}-1} \xi _{r, k}(t)\, {_0^C}\mathsf{D}_t^{\mu _{r, n-k}}y_r(t)\, dt \\ & - \frac{1}{\Gamma(\mu _{r, n})}\int _{0}^x (x-t)^{\mu _{r, n}-1}\xi _{r, n}(t)\, y_r(t)\, dt \\ & +\sum\limits _{j = 1}^{2} \frac{\nu _{r, j}}{\Gamma(\mu _{r, n})} \int _0^x (x-t)^{\mu _{r, n}-1}\int _0^t \kappa _{r, j}(t, s)\, {_0^C}\mathsf{D}_s^{\gamma _{r, 2-j}}\, y_j(s)\, ds\, dt \\ & +\frac{1}{\Gamma(\mu _{r, n})} \int _0^x (x-t)^{\mu _{r, n}-1}\, \mathsf{f}_r(t)\, dt \end{aligned} \end{equation}$

(3.1)

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

$\begin{equation*} \Vert z \Vert _C = \max \{\sup\limits_{x \in \mathbf{I}} \vert z(x) \vert, \sup\limits_{x \in \mathbf{I}} \vert {_0^C}\mathsf{D}_x^{\mu} z(x) \vert \}, \quad \forall z \in C(\mathbf{I}, X) \end{equation*}$

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

$\begin{equation} \begin{aligned} & \mathcal{M}_{r, j} = \sup\limits_{(x, t)\in \mathbf{I}} \vert \kappa _{r, j}(x, t) \vert, \, \, r, j = 1, 2, \\ & \mathcal{P}_{r, k} = \sup\limits_{(x)\in \mathbf{I}} \vert \xi_{r, k}(x) \vert, \, \, r = 1, 2, k = 1, \cdots, n, \\ & \mathcal{F}_r = \sup\limits_{(x)\in \mathbf{I}} \vert \mathsf{f}_{r}(x) \vert, \, \, r = 1, 2 \end{aligned} \end{equation}$

(3.2)

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

$\begin{equation} \frac{\mathcal{P}_{r, n}}{\Gamma(\mu _{r, n}+1)} < 1, \quad \sum\limits_{k = 1}^{n}\frac{\mathcal{P}_{r, k}}{\Gamma(\mu _{r, n}+1)}+\sum\limits_{j = 1}^{2}\frac{ \vert \nu _{r, j} \vert \, \mathcal{M}_{r, j}}{\Gamma(\mu _{r, n}+2)} < 1 \end{equation}$

(3.3)

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

$\begin{equation} q\geq \frac{ \vert y_{r, 0} \vert +\frac{\mathcal{F}_r}{\Gamma(\mu _{r, n}+1)}}{1-\bigg \{ \sum\limits_{k = 1}^n \frac{\mathcal{P}_{r, k}}{\Gamma(\mu _{r, n}+1)}+\sum\limits_{j = 1}^2 \frac{\vert \nu _{r, j}\vert \, \mathcal{M}_{r, j}}{\Gamma(\mu _{r, n}+2)}\bigg\}} \end{equation}$

(3.4)

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:

$\begin{equation*} \mathcal{B}_1 y_r(x) = y_{r, 0}-\frac{1}{\Gamma(\mu _{r, n})}\int _{0}^x (x-t)^{\mu _{r, n}-1}\xi _{r, n}(t)\, y_r(t)\, dt+ \frac{1}{\Gamma(\mu _{r, n})} \int _0^x (x-t)^{\mu _{r, n}-1}\, \mathsf{f}_r(t)\, dt \end{equation*}$

$\begin{align*} \mathcal{B}_2 y_r(x)& = \sum\limits _{j = 1}^{2} \frac{\nu _{r, j}}{\Gamma(\mu _{r, n})} \int _0^x (x-t)^{\mu _{r, n}-1}\int _0^t \kappa _{r, j}(t, s)\, {_0^C}\mathsf{D}_s^{\gamma _{r, 2-j}}\, y_j(s)\, ds\, dt\\ & \quad -\sum\limits _{k = 1}^{n-1} \frac{1}{\Gamma(\mu _{r, n})} \int _0^x (x-t)^{\mu _{r, n}-1} \xi _{r, k}(t)\, {_0^C}\mathsf{D}_t^{\mu _{r, n-k}}y_r(t)\, dt \end{align*}$

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:

$\begin{align*} \Vert \mathcal{B}_1 y_r+\mathcal{B}_2 u_r \Vert _C & \leq \vert y_{r, 0} \vert +\frac{\mathcal{P}_{r, n}\, \Vert y_r \Vert _C}{\Gamma(\mu _{r, n}+1)}+\frac{\mathcal{F}_{r}}{\Gamma(\mu _{r, n}+1)}+\sum\limits_{j = 1}^2\frac{\vert \nu _{r, j} \vert \, \mathcal{M}_{r, j}\, \Vert u_r \Vert _C}{\Gamma(\mu _{r, n}+2)}+\sum\limits_{k = 1}^{n-1}\frac{\mathcal{P}_{r, k}\, \Vert u_r \Vert _C}{\Gamma(\mu _{r, n}+1)}\\ & \leq \vert y_{r, 0} \vert +\frac{\mathcal{F}_{r}}{\Gamma(\mu _{r, n}+1)}+\bigg\{ \frac{\mathcal{P}_{r, n}\, \Vert y_r \Vert _C}{\Gamma(\mu _{r, n}+1)}+\sum\limits_{k = 1}^{n-1}\frac{\mathcal{P}_{r, k}\, \Vert u_r \Vert _C}{\Gamma(\mu _{r, n}+1)}+\sum\limits_{j = 1}^2\frac{\vert \nu _{r, j} \vert \, \mathcal{M}_{r, j}\, \Vert u_r \Vert _C}{\Gamma(\mu _{r, n}+2)}\bigg\}q\\ & \leq q \end{align*}$

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

$\begin{align*} \Vert \mathcal{B}_1 y_r- \mathcal{B}_1 u_r \Vert _C & = \Vert \frac{1}{\Gamma(\mu _{r, n})}\int _0^x (x-t)^{\mu _{r, n}-1}\xi _{r, n}(t)\, (u_r(t)-y_r(t))\, dt \Vert _C\\ &\leq \frac{\mathcal{P}_{r, n}}{\Gamma(\mu_{r, n}+1)} \Vert y_r-u_r \Vert _C \end{align*}$

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

$\begin{align*} \Vert \mathcal{B}_2 y_r \Vert _C& = \bigg \Vert \sum\limits _{j = 1}^{2} \frac{\nu _{r, j}}{\Gamma(\mu _{r, n})} \int _0^x (x-t)^{\mu _{r, n}-1}\int _0^t \kappa _{r, j}(t, s)\, {_0^C}\mathsf{D}_s^{\gamma _{r, 2-j}}\, y_j(s)\, ds\, dt\\ & \quad -\sum\limits _{k = 1}^{n-1} \frac{1}{\Gamma(\mu _{r, n})} \int _0^x (x-t)^{\mu _{r, n}-1} \xi _{r, k}(t)\, {_0^C}\mathsf{D}_t^{\mu _{r, n-k}}y_r(t)\, dt \bigg \Vert _C\\ & \leq \sum\limits_{j = 1}^2 \frac{\vert \nu _{r, j} \vert\, \mathcal{M}_{r, j}}{\Gamma(\mu_{r, n}+1)}\Vert y_j \Vert _C+\sum\limits_{k = 1}^{n-1} \frac{\mathcal{P}_{r, k}}{\Gamma(\mu _{r, n}+2)}\Vert y_r \Vert _C\\ & \leq \bigg\{ \sum\limits_{j = 1}^2 \frac{\vert \nu _{r, j} \vert \, \mathcal{M}_{r, j}}{\Gamma(\mu_{r, n}+1)}+\sum\limits_{k = 1}^{n-1} \frac{\mathcal{P}_{r, k}}{\Gamma(\mu _{r, n}+2)}\bigg \} q \end{align*}$

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

$\begin{align*} \mathcal{B}_2 y_r(x_2)-\mathcal{B}_2 y_r(x_1)& = \sum\limits_{j = 1}^2 \frac{\nu _{r, j}}{\Gamma(\mu _{r, n})}\int _0^{x_2}(x_2-t)^{\mu _{r, n}-1}\int _0^t \kappa _{r, j}(t, s)\, {_0^C}\mathsf{D}_s^{\gamma _{r, 2-j}}\, y_j(s)\, ds\, dt \\ & \quad -\sum\limits _{k = 1}^{n-1} \frac{1}{\Gamma(\mu _{r, n})} \int _0^{x_2} (x_2-t)^{\mu _{r, n}-1} \xi _{r, k}(t)\, {_0^C}\mathsf{D}_t^{\mu _{r, n-k}}y_r(t)\, dt\\ & \quad - \sum\limits_{j = 1}^2 \frac{\nu _{r, j}}{\Gamma(\mu _{r, n})}\int _0^{x_1}(x_1-t)^{\mu _{r, n}-1}\int _0^t \kappa _{r, j}(t, s)\, {_0^C}\mathsf{D}_s^{\gamma _{r, 2-j}}\, y_j(s)\, ds\, dt \\ & \quad +\sum\limits _{k = 1}^{n-1} \frac{1}{\Gamma(\mu _{r, n})} \int _0^{x_1} (x_1-t)^{\mu _{r, n}-1} \xi _{r, k}(t)\, {_0^C}\mathsf{D}_t^{\mu _{r, n-k}}y_r(t)\, dt\\ & = \sum\limits_{j = 1}^2 \frac{\nu _{r, j}}{\Gamma(\mu _{r, n})}\int _0^{x_2}\bigg [(x_2-t)^{\mu _{r, n}-1}-(x_1-t)^{\mu _{r, n}-1}\bigg] \int _0^t \kappa _{r, j}(t, s)\, {_0^C}\mathsf{D}_s^{\gamma _{r, 2-j}}\, y_j(s)\, ds\, dt \\ & \quad +\sum\limits_{j = 1}^2 \frac{\nu _{r, j}}{\Gamma(\mu _{r, n})}\int _{x_1}^{x_2}(x_1-t)^{\mu _{r, n}-1}\int _0^t \kappa _{r, j}(t, s)\, {_0^C}\mathsf{D}_s^{\gamma _{r, 2-j}}\, y_j(s)\, ds\, dt \\ & \quad -\sum\limits _{k = 1}^{n-1} \frac{1}{\Gamma(\mu _{r, n})} \int _0^{x_2} \bigg [(x_2-t)^{\mu _{r, n}-1}-(x_1-t)^{\mu _{r, n}-1} \bigg] \xi _{r, k}(t)\, {_0^C}\mathsf{D}_t^{\mu _{r, n-k}}y_r(t)\, dt\\ & \quad - \sum\limits _{k = 1}^{n-1} \frac{1}{\Gamma(\mu _{r, n})} \int _{x_1}^{x_2} (x_1-t)^{\mu _{r, n}-1} \xi _{r, k}(t)\, {_0^C}\mathsf{D}_t^{\mu _{r, n-k}}y_r(t)\, dt \end{align*}$

By taking norm, one gets

$\begin{align*} \Vert \mathcal{B}_2 y_r(x_2)-\mathcal{B}_2 y_r(x_1) \Vert _C &\leq \sum\limits_{j = 1}^2 \frac{\vert \nu _{r, j}\vert \, \mathcal{M}_{r, j}\, q}{\Gamma(\mu _{r, n})}\int _0^{x_2}\bigg [(x_2-t)^{\mu _{r, n}-1}-(x_1-t)^{\mu _{r, n}-1}\bigg] \, t \, dt\\ & \quad + \sum\limits_{j = 1}^2 \frac{\vert \nu _{r, j}\vert \, \mathcal{M}_{r, j}\, q}{\Gamma(\mu _{r, n})}\int _0^{x_2} (x_1-t)^{\mu _{r, n}-1}\, t \, dt\\ & \quad +\sum\limits _{k = 1}^{n-1} \frac{\mathcal{P}_{r, k}\, q}{\Gamma(\mu _{r, n})} \int _0^{x_2} \bigg [(x_2-t)^{\mu _{r, n}-1}-(x_1-t)^{\mu _{r, n}-1} \bigg] \, dt \\ & \quad +\sum\limits _{k = 1}^{n-1} \frac{\mathcal{P}_{r, k}\, q}{\Gamma(\mu _{r, n})} \int _{x_1}^{x_2} (x_1-t)^{\mu _{r, n}-1}\, dt \end{align*}$

Using a change of variable, one gets:

$\begin{align*} \Vert \mathcal{B}_2 y_r(x_2)-\mathcal{B}_2 y_r(x_1) \Vert _C &\leq \sum\limits_{j = 1}^2 \frac{\vert \nu _{r, j}\vert \, \mathcal{M}_{r, j}\, q}{\Gamma(\mu _{r, n})} \bigg \{-\int _{x_2}^0 (x_2-u)\, u^{\mu _{r, n}-1}\, du + \int _{x_1}^{x_1-x_2} (x_1-u)\, u^{\mu _{r, n}-1}\, du\\ & \quad -\int _{0}^{x_1-x_2} (x_1-u)\, u^{\mu _{r, n}-1}\, du \bigg\}+\sum\limits _{k = 1}^{n-1} \frac{\mathcal{P}_{r, k}\, q}{\Gamma(\mu _{r, n})} \bigg \{ \int _{x_2}^{0} u^{\mu _{r, n}-1}\, du\\ & \quad -\int _{x_1}^{x_1-x_2} u^{\mu _{r, n}-1}\, du +\int _{0}^{x_1-x_2} u^{\mu _{r, n}-1}\, du \bigg \} \\ & = \sum\limits_{j = 1}^2 \frac{\vert \nu _{r, j}\vert \, \mathcal{M}_{r, j}\, q}{\Gamma(\mu _{r, n})} \bigg ( \frac{x_2^{\mu _{r, n}+1}}{\mu _{r, n}\, (\mu _{r, n}+1)}+\frac{x_1\, (x_1-x_2)^{\mu _{r, n}}}{\mu _{r, n}}-\frac{(x_1-x_2)^{\mu _{r, n}+1}}{\mu _{r, n}+1}\\ & \quad -\frac{x_1^{\mu _{r, n}+1}}{\mu _{r, n}\, (\mu _{r, n}+1)}-\frac{x_1\, (x_1-x_2)^{\mu _{r, n}}}{\mu _{r, n}}+\frac{(x_1-x_2)^{\mu _{r, n}+1}}{\mu _{r, n}+1}\bigg )\\ & \quad +\sum\limits _{k = 1}^{n-1} \frac{\mathcal{P}_{r, k}\, q}{\Gamma(\mu _{r, n})} \bigg ( \frac{x_2^{\mu _{r, n}}}{\mu _{r, n}} +\frac{(x_1-x_2)^{\mu _{r, n}}}{\mu _{r, n}}-\frac{x_1^{\mu _{r, n}}}{\mu _{r, n}}-\frac{(x_1-x_2)^{\mu _{r, n}}}{\mu _{r, n}} \bigg )\\ & = \sum\limits_{j = 1}^2 \frac{\vert \nu _{r, j} \vert \, \mathcal{M}_{r, j}\, q}{\Gamma(\mu _{r, n}+2)}\bigg( x_2^{\mu _{r, n}+1}-x_1^{\mu _{r, n}+1}\bigg) +\sum\limits_{k = 1}^{n-1}\frac{\mathcal{P}_{r, k}\, q}{\Gamma(\mu _{r, n}+1)}\bigg(x_1^{\mu _{r, n}}-x_2^{\mu _{r, n}}\bigg) \end{align*}$

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:

$\begin{align*} \Vert \mathcal{E}y_r-\mathcal{E}u_r \Vert _C & = \bigg \Vert \sum\limits_{k = 1}^{n-1} \frac{1}{\Gamma(\mu _{r, n})}\int _0^x (x-t)^{\mu _{r, n-1}} \, \xi _{r, k}(t)\, {_0^C}\mathsf{D}_s^{\mu _{r, n-k}}\, (u_r(t)-y_r(t))\, dt \\ & \quad +\frac{1}{\Gamma(\mu _{r, n})}\int _0^x(x-t)^{\mu _{r, n}-1}\, \xi _{r, n}(t)\, (u_r(t)-y_r(t))\, dt \\ & \quad +\sum\limits _{j = 1}^2 \frac{\nu _{r, j}}{\Gamma(\mu _{r, n})}\int _0^x (x-t)^{\mu _{r, n}-1} \int _0^t \kappa _{r, j}(t, s) \, {_0^C}\mathsf{D}_s^{\gamma _{r, 2-j}}\, (u_j(s)-y_j(s))\, ds\, dt \bigg\Vert _C \\ & \leq \sum\limits_{k = 1}^{n-1} \frac{\mathcal{P}_{r, k}}{\Gamma(\mu _{r, n}+1)}\Vert y_r - u_r \Vert _C+ \frac{\mathcal{P}_{r, n}}{\Gamma(\mu _{r, n}+1)}\Vert y_r - u_r \Vert _C+\sum\limits_{j = 1}^{2} \frac{\vert \nu _{r, j} \vert \, \mathcal{M}_{r, j}}{\Gamma(\mu _{r, n}+2)} \Vert y_r - u_r \Vert _C \end{align*}$

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:

$\begin{align*} & \mathcal{G}_{i+1}^{\sigma}(x) = \frac{2(i+\sigma )}{i+1}(2x-1) \mathcal{G}_{i}^{\sigma}(x)-\frac{i+2\sigma -1}{i+1}\mathcal{G}_{i-1}^{\sigma}(x), \quad i = 1, 2, \cdots, \quad x \in\mathbf{I}, \\ & \mathcal{G}_{0}^{\sigma}(x) = 1, \quad \mathcal{G}_{1}^{\sigma}(x) = 2 \sigma (2x-1) \end{align*}$

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:

$\begin{equation} \int _0^1 \mathcal{G}_i^{\sigma}(x)\, \mathcal{G}_j^{\sigma}(x)\, \omega ^{\sigma}(x)\, dx = \left \{ \begin{array}{ll} 0, & i \neq j, \\ \mathsf{h}_i^{\sigma}, & i = j \end{array} \right. \end{equation}$

(4.1)

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

$\begin{equation} \mathsf{h}_i^{\sigma} = \frac{\Gamma(\sigma +\frac{1}{2})^2 \, \Gamma(i+2 \sigma)}{\Gamma(2\sigma)^2\, (2i+2\sigma)\, i!}, \quad i = 0, 1, 2, \cdots \end{equation}$

(4.2)

These polynomials can be represented in series form as

$\begin{equation} \mathcal{G}_i^{\sigma}(x) = \sum \limits _{k = 0}^i \rho _{k, i}^{\sigma} \, x^k, \quad x\in \mathbf{I} \end{equation}$

(4.3)

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

$\begin{equation} \rho _{k, i}^{\sigma} = \frac{(-1)^{i-k}\, \Gamma(\sigma +\frac{1}{2})\, \Gamma(i+k+2\sigma)}{\Gamma(2 \sigma)\, \Gamma(i+\sigma +\frac{1}{2})\, (i-k)!\, k!}, \quad i = 0, 1, \cdots, \, \, k = 0, 1, \cdots, i \end{equation}$

(4.4)

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]:

$\begin{equation} \mathcal{G}_i^{\sigma}(x) = \frac{i!\, \Gamma(\sigma +\frac{1}{2})}{\Gamma(i+\sigma +\frac{1}{2})} \mathcal{J}_i^{\sigma -\frac{1}{2}, \sigma -\frac{1}{2}}(x) \end{equation}$

(4.5)

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

$\begin{equation} \varsigma _{i}^{\alpha, \beta} = \frac{\Gamma(i+\alpha +1)\, \Gamma(i+\beta +1)}{(2i +\alpha +\beta +1)\, \Gamma(i+1)\, \Gamma(i+\alpha +\beta +1)}, \quad i = 0, 1, 2, \cdots \end{equation}$

(4.6)

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

$\begin{equation} \frac{d^l \, \mathcal{J}_i^{\alpha, \beta}(x)}{d\, x^l} = \frac{\Gamma(l+i+\alpha +\beta)}{\Gamma(i+\alpha +\beta +1)}\mathcal{J}_{i-l}^{\alpha +l, \beta +l}(x) \end{equation}$

(4.7)

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

$\begin{equation} y(x) = \sum\limits_{k = 0}^{\infty} \mathsf{y}_k \, \mathcal{G}_k^{\sigma}(x) \end{equation}$

(4.8)

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

$\begin{equation} \mathsf{y}_k = \frac{1}{\mathsf{h}_k^{\sigma}} \int _0 ^ 1 y(x)\, \mathcal{G}_k^{\sigma}(x)\, \omega ^{\sigma}(x)\, dx. \end{equation}$

(4.9)

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:

$\begin{equation} y(x)\approx y_N(x) = \sum\limits_{k = 0}^{N} \mathsf{y}_k \, \mathcal{G}_k^{\sigma}(x) = Y^T \, \Lambda ^{\sigma} (x) = {\Lambda ^{\sigma} }^T (x)\, Y, \end{equation}$

(4.10)

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

$\begin{equation} Y = [\mathsf{y}_0, \mathsf{y}_1, \mathsf{y}_2, \cdots, \mathsf{y}_N]^T, \quad \Lambda ^{\sigma} (x) = [\mathcal{G}_0^{\sigma}(x), \mathcal{G}_1^{\sigma}(x), \mathcal{G}_2^{\sigma}(x), \cdots, \mathcal{G}_N^{\sigma}(x)]^T \end{equation}$

(4.11)

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:

$\begin{equation} \int _0^x \mathcal{G}_i^{\sigma}(s)\, ds = \sum \limits_{k = 0}^i \rho _{k, i}^{\sigma} \, \int _0^x s^k\, ds = \sum \limits_{k = 0}^i \rho _{k, i}^{\sigma} \frac{x^{k+1}}{k+1} \end{equation}$

(4.12)

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

$\begin{equation*} x^{k+1}\approx \sum \limits_{j = 0}^N d_{j, k+1} ^{\sigma} \, \mathcal{G}_j^{\sigma}(x) \end{equation*}$

where

$\begin{align*} d_{l, k+1} ^{\sigma}& = \frac{1}{\mathsf{h}_j ^{\sigma}} \int _0^1 x^{k+1}\, \mathcal{G}_j^{\sigma}(x)\, \omega ^{\sigma}(x)\, dx\\ & = \frac{1}{\mathsf{h}_j ^{\sigma}}\sum\limits_{l = 0}^j \rho_{l, j}^{\sigma} \int _0^1 x^{k+1}\, x^l\, x^{\sigma -\frac{1}{2}}(1-x)^{\sigma -\frac{1}{2}}\, dx\\ & = \frac{1}{\mathsf{h}_j ^{\sigma}}\sum\limits_{l = 0}^j \rho_{l, j}^{\sigma} \frac{\Gamma(k+l+\sigma +\frac{3}{2})\, \Gamma(\sigma +\frac{1}{2})}{\Gamma(k+l+2\sigma +2)} \end{align*}$

Therefore, the integral in (4.12) will be as

$\begin{equation} \int _0^x \mathcal{G}_i^{\sigma}(s)\, ds \approx \sum \limits_{j = 0}^N \bigg\{ \sum \limits_{k = 0}^i \sum\limits_{l = 0}^j \frac{\rho _{k, i}^{\sigma}\, \rho _{l, j}^{\sigma}\, \Gamma(k+l+\sigma +\frac{3}{2})\Gamma(\sigma +\frac{1}{2})}{\mathsf{h}_j ^{\sigma} \, (k+1)\, \Gamma(k+l+2\sigma +2) } \bigg\} \mathcal{G}_j^{\sigma}(x), \quad i = 0, 1, \cdots, N \end{equation}$

(4.13)

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

$\begin{equation} \int _0^x \Lambda ^{\sigma}(s)\, ds \approx \mathbf{\Theta}^{(\sigma)}\, \Lambda ^{\sigma}(x) \end{equation}$

(4.14)

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

$\begin{equation*} \mathbf{\Theta}^{(\sigma)}_{i, j} = \sum \limits_{k = 0}^i \sum\limits_{l = 0}^j \frac{\rho _{k, i}^{\sigma}\, \rho _{l, j}^{\sigma}\, \Gamma(k+l+\sigma +\frac{3}{2})\Gamma(\sigma +\frac{1}{2})}{\mathsf{h}_j ^{\sigma} \, (k+1)\, \Gamma(k+l+2\sigma +2) }, \quad i, j = 0, 1, \cdots, N \end{equation*}$

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:

$\begin{equation} {_0^{\mathsf{RL}}}\mathsf{I}_x^{\mu} (\mathcal{G}_i^{\sigma}(x)) = \sum\limits_{k = 0}^i \rho _{k, i}^{\sigma} \, {_0^{\mathsf{RL}}}\mathsf{I}_x^{\mu} (x^k) = \sum\limits_{k = 0}^i \frac{\rho _{k, i}^{\sigma}\, \Gamma(k+1)}{\Gamma(k+\mu +1)} x^{k+\mu} \end{equation}$

(4.15)

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

$\begin{equation*} x^{k+\mu}\approx \sum\limits_{j = 0}^N b_{j, k+\mu}^{\sigma}\, \mathcal{G}_j^{\sigma}(x) \end{equation*}$

where

$\begin{align*} b_{j, k+\mu}^{\sigma}& = \frac{1}{\mathsf{h}_j ^{\sigma}}\int _0^1 x^{k+\mu}\, \mathcal{G}_j^{\sigma}(x)\, \omega ^{\sigma}(x)\, dx\\ & = \frac{1}{\mathsf{h}_j ^{\sigma}} \sum\limits _{l = 0}^j \rho _{l, j}^{\sigma} \int _0^1 x^{k+l+\mu +\sigma -\frac{1}{2}}\, (1-x)^{\sigma -\frac{1}{2}}\, dx\\ & = \frac{1}{\mathsf{h}_j ^{\sigma}} \sum\limits _{l = 0}^j \frac{\rho _{l, j}^{\sigma} \, \Gamma(k+l+\mu +\sigma +\frac{1}{2})\, \Gamma(\sigma +\frac{1}{2})}{\Gamma(k+l+\mu +2\sigma +1)} \end{align*}$

So, (4.15) is written as

$\begin{equation*} {_0^{\mathsf{RL}}}\mathsf{I}_x^{\mu} (\mathcal{G}_i^{\sigma}(x))\approx \sum\limits_{j = 0}^N \bigg \{ \sum\limits_{k = 0}^i \sum\limits_{l = 0}^j \frac{\rho _{k, i}^{\sigma}\, \rho _{l, j}^{\sigma} \, \Gamma(k+1)\, \Gamma(k+l+\mu +\sigma +\frac{1}{2})\, \Gamma(\sigma +\frac{1}{2})}{\mathsf{h}_j ^{\sigma}\, \Gamma(k+\mu +1)\, \Gamma(k+l+\mu +2\sigma +1)}\bigg\} \, \mathcal{G}_j^{\sigma}(x), \quad i = 0, 1, \cdots, N \end{equation*}$

Therefore, one gets

$\begin{equation} {_0^{\mathsf{RL}}}\mathsf{I}_x^{\mu} (\Lambda ^{\sigma}(x))\approx \mathbf{\Theta}^{(\mu, \sigma)}\, \Lambda ^{\sigma}(x) \end{equation}$

(4.16)

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

$\begin{equation*} \mathbf{\Theta}^{(\mu, \sigma)}_{i, j} = \sum\limits_{k = 0}^i \sum\limits_{l = 0}^j \frac{\rho _{k, i}^{\sigma}\, \rho _{l, j}^{\sigma} \, \Gamma(k+1)\, \Gamma(k+l+\mu +\sigma +\frac{1}{2})\, \Gamma(\sigma +\frac{1}{2})}{\mathsf{h}_j ^{\sigma}\, \Gamma(k+\mu +1)\, \Gamma(k+l+\mu +2\sigma +1)}, \quad i, j = 0, 1, \cdots, N \end{equation*}$

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:

$\begin{equation} \Lambda ^{\sigma}(x)\, {\Lambda ^{\sigma}}^T(x)\, V \approx V^*\, \Lambda ^{\sigma}(x) \end{equation}$

(4.17)

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:

$\begin{equation} \mathcal{G}_j^{\sigma}(x)\, \mathcal{G}_k^{\sigma}(x) = \sum\limits_{r = 0}^{j+k} \pi _{r}^{(j, k)}\, x^r \end{equation}$

(4.18)

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

$\begin{equation} \begin{aligned} \mathsf{q}_{i, j, k}& = \int _0^1 \mathcal{G}_i^{\sigma}(x)\, \mathcal{G}_j^{\sigma}(x)\, \mathcal{G}_k^{\sigma}(x)\, \omega ^{\sigma}(x)\, dx \\ & = \sum\limits _{r = 0}^{j+k} \pi^{(j, k)}_{r} \int _0^1 x^r\, \mathcal{G}_i^{\sigma}(x) \, \omega ^{\sigma}(x)\, dx\\ & = \sum\limits _{r = 0}^{j+k} \sum\limits _{l = 0}^{i} \frac{\pi^{(j, k)}_{r}\, \rho _{l, i}^{\sigma}\, \Gamma(r+l+\sigma +\frac{1}{2})\, \Gamma(\sigma +\frac{1}{2})}{\Gamma(r+l+2\sigma +1)} \end{aligned} \end{equation}$

(4.19)

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

$\begin{equation} V^*_{j, k} = \frac{1}{\mathsf{h}_k^{\sigma}} \sum\limits _{i = 0}^N V_i\, \mathsf{q}_{i, j, k}, \quad j, k = 0, 1, \cdots, N \end{equation}$

(4.20)

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:

$\begin{equation*} \int _0^1 \Lambda ^T(x)\, K\, W^*\, \mathbf{\Theta}^{(\sigma)}\, \Lambda (x)\, dx \end{equation*}$

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:

$\begin{equation} \int _0^1\Lambda ^T(x)\, K\, W^*\, \mathbf{\Theta}^{(\sigma)}\, \Lambda (x)\, dx \approx \Delta \end{equation}$

(4.21)

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

$\begin{equation} \begin{aligned} &\Delta _s = \sum\limits _{m = 1}^{N+1} \sum\limits _{n = 1}^{N+1} \mathsf{a}_{m, n}\, \mathsf{q}_{s-1, m-1, n-1}, \quad s = 1, 2, \cdots, N+1, \\ & \mathsf{a}_{m, n} = \sum\limits _{i = 1}^{N+1} \sum\limits _{j = 1}^{N+1} K_{m, j}\, W^*_{j, i}\, \mathbf{\Theta}^{(\sigma)}_{i, n}, \quad m, n = 1, 2, \cdots, N+1 \end{aligned} \end{equation}$

(4.22)

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

$\begin{equation} \int _0^1 \Lambda(x) \, \Lambda ^T(x)\, \omega^{\sigma}(x)\, dx = \mathcal{Q}^{\sigma} \end{equation}$

(4.23)

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:

$\begin{equation} \begin{aligned} & {_0^C}\mathsf{D}_x^{0.9} y_0(x)+{_0^C}\mathsf{D}_x^{0.6} y_0(x)+2x\, y_0(x)-\int _0^x 5\, \exp(x)\, {_0^C}\mathsf{D}_t^{0.3} y_0(t)\, dt -\mathsf{f}_1(x) = 0, \\ & {_0^C}\mathsf{D}_x^{0.8} y_1(x)+x\, {_0^C}\mathsf{D}_x^{0.5} y_1(x)-\int _0^x (\sin(x)-t)\, {_0^C}\mathsf{D}_t^{0.8} y_0(t)\, dt -\mathsf{f}_2(x) = 0 \end{aligned} \end{equation}$

(5.1)

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:

$\begin{equation} {_0^C}\mathsf{D}_x^{0.9} y_0(x) \approx \mathcal{Y}_0 ^T\, \Lambda ^{\sigma}(x), \quad {_0^C}\mathsf{D}_x^{0.8} y_1(x) \approx \mathcal{Y}_1 ^T\, \Lambda ^{\sigma}(x) \end{equation}$

(5.2)

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:

$\begin{equation} \begin{aligned} & y_0(x) \approx \mathcal{Y}_0 ^T\, \mathbf{\Theta}^{(0.9, \sigma)}\, \Lambda^{\sigma}(x)+y_0(0) = \mathcal{W}_1^T \, \Lambda ^{\sigma}(x), \, \, \mathcal{W}_1 = {\mathbf{\Theta}^{(0.9, \sigma)}}^T\, \mathcal{Y}_0, \\ & y_1(x) \approx \mathcal{Y}_1 ^T\, \mathbf{\Theta}^{(0.8, \sigma)}\, \Lambda^{\sigma}(x)+y_1(0) = \mathcal{Y}_1 ^T\, \mathbf{\Theta}^{(0.8, \sigma)}\, \Lambda^{\sigma}(x)+E_1^T\, \Lambda(x) = \mathcal{Z}_1^T \, \Lambda ^{\sigma}(x), \, \, \mathcal{Z}_1 = E_1+{\mathbf{\Theta}^{(0.8, \sigma)}}^T\, \mathcal{Y}_1 \end{aligned} \end{equation}$

(5.3)

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

$\begin{equation} {_0^C}\mathsf{D}_x^{0.9} y_0(x) = {_0^C}\mathsf{D}_x^{0.3}(\mathsf{D}_x^{0.6} y_0(x)) \approx \mathcal{Y}_0 ^T\, \Lambda ^{\sigma}(x) \end{equation}$

(5.4)

so, one has

$\begin{equation} \begin{aligned} {_0^C}\mathsf{D}_x^{0.6} y_0(x) \approx \mathcal{Y}_0 ^T\, \mathbf{\Theta}^{(0.3, \sigma)}\, \Lambda ^{\sigma}(x)+{_0^C}\mathsf{D}_x^{0.6} y_0(0) = \mathcal{W}_2^T\, \Lambda ^{\sigma}(x), \, \, \mathcal{W}_2 = {\mathbf{\Theta}^{(0.3, \sigma)}}^T\, \mathcal{Y}_0 \end{aligned} \end{equation}$

(5.5)

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

$\begin{equation} {_0^C}\mathsf{D}_x^{0.9} y_0(x) = {_0^C}\mathsf{D}_x^{0.6}(\mathsf{D}_x^{0.3} y_0(x)) \approx \mathcal{Y}_0 ^T\, \Lambda ^{\sigma}(x) \quad \Longrightarrow \quad \mathsf{D}_x^{0.3} y_0(x) \approx \mathcal{W}_3^T\, \Lambda ^{\sigma}(x), \, \, \mathcal{W}_3 = {\mathbf{\Theta}^{(0.6, \sigma)}}^T\, \mathcal{Y}_0, \end{equation}$

(5.6)

$\begin{equation} {_0^C}\mathsf{D}_x^{0.9} y_0(x) = {_0^C}\mathsf{D}_x^{0.1}(\mathsf{D}_x^{0.8} y_0(x)) \approx \mathcal{Y}_0 ^T\, \Lambda ^{\sigma}(x) \quad \Longrightarrow \quad \mathsf{D}_x^{0.8} y_0(x) \approx \mathcal{W}_4^T\, \Lambda ^{\sigma}(x), \, \, \mathcal{W}_4 = {\mathbf{\Theta}^{(0.1, \sigma)}}^T\, \mathcal{Y}_0 \end{equation}$

(5.7)

Similarly, one has

$\begin{equation} {_0^C}\mathsf{D}_x^{0.8} y_1(x) = {_0^C}\mathsf{D}_x^{0.3}(\mathsf{D}_x^{0.5} y_1(x)) \approx \mathcal{Y}_1 ^T\, \Lambda ^{\sigma}(x) \quad \Longrightarrow \quad \mathsf{D}_x^{0.5} y_1(x) \approx \mathcal{Z}_2^T\, \Lambda ^{\sigma}(x), \, \, \mathcal{Z}_2 = {\mathbf{\Theta}^{(0.3, \sigma)}}^T\, \mathcal{Y}_1 \end{equation}$

(5.8)

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

$\begin{equation} \begin{aligned} & 2\, x \approx \mathcal{F}_1\, \Lambda ^{\sigma}(x), \quad x \approx\mathcal{F}_2\, \Lambda ^{\sigma}(x), \\ & 2\, x \, y_0(x) \approx \mathcal{F}_1 \, \Lambda ^{\sigma}(x) \, {\Lambda ^{\sigma} }^T(x)\, \mathcal{W}_1 \approx \mathcal{F}_1 \, \mathcal{W}_1^*\, \Lambda ^{\sigma}(x), \\ & x\, \mathsf{D}_x^{0.5} y_1(x) \approx \mathcal{F}_2 \, \Lambda ^{\sigma}(x) \, {\Lambda ^{\sigma} }^T(x)\, \mathcal{Z}_2 \approx \mathcal{F}_2 \, \mathcal{Z}_2^*\, \Lambda ^{\sigma}(x) \end{aligned} \end{equation}$

(5.9)

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:

$\begin{equation} \begin{aligned} & 5\, \exp(x) \approx {\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_1\, \Lambda ^{\sigma}(t), \quad \sin(x)-t \approx {\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_2\, \Lambda ^{\sigma}(t), \\ & \int _0^x 5\, \exp(x)\, \mathsf{D}_t^{0.3} y_0(t)\, dt \approx {\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_1\, \mathcal{W}_3^*\, \int _0^x \Lambda ^{\sigma}(t) dt \approx {\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_1\, \mathcal{W}_3^*\, \mathbf{\Theta }^{(\sigma)}\, \Lambda ^{\sigma}(x), \\ & \int _0^x (\sin(x)-t)\, \mathsf{D}_t^{0.8} y_0(t)\, dt \approx {\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_2\, \mathcal{W}_4^*\, \int _0^x \Lambda ^{\sigma}(t) dt \approx {\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_2\, \mathcal{W}_4^*\, \mathbf{\Theta} ^{(\sigma)}\, \Lambda ^{\sigma}(x) \end{aligned} \end{equation}$

(5.10)

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:

$\begin{equation} \mathsf{f}_1(x) \approx \mathcal{F}_3^T\, \Lambda ^{\sigma}(x), \quad \mathsf{f}_2(x) \approx \mathcal{F}_4^T\, \Lambda ^{\sigma}(x) \end{equation}$

(5.11)

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:

$\begin{equation} \begin{aligned} &\mathcal{R}_1(x) = \mathcal{Y}_0 ^T\, \Lambda ^{\sigma}(x)+ \mathcal{W}_2 ^T\, \Lambda ^{\sigma}(x)+\mathcal{F}_1^T \, \mathcal{W}_1^* \, \Lambda ^{\sigma}(x)-{\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_1\, \mathcal{W}_3^*\, \mathbf{\Theta}^{(\sigma)} \, \Lambda ^{\sigma}(x)-\mathcal{F}_3^T\, \Lambda ^{\sigma}(x) \approx 0, \\ & \mathcal{R}_2(x) = \mathcal{Y}_1 ^T\, \Lambda ^{\sigma}(x)+ \mathcal{F}_2^T \, \mathcal{Z}_2^* \, \Lambda ^{\sigma}(x)-{\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_2\, \mathcal{W}_4^*\, \mathbf{\Theta}^{(\sigma)} \, \Lambda ^{\sigma}(x)-\mathcal{F}_4^T\, \Lambda ^{\sigma}(x) \approx 0 \end{aligned} \end{equation}$

(5.12)

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

$\begin{equation} \left\{ \begin{array}{l} \mathcal{Y}_0 ^T\, \mathcal{Q}+ \mathcal{W}_2 ^T\, \mathcal{Q}+\mathcal{F}_1^T \, \mathcal{W}_1^* \, \mathcal{Q}-\Delta _1-\mathcal{F}_3^T\, \mathcal{Q} \approx 0, \\ \mathcal{Y}_1 ^T\, \mathcal{Q}+ \mathcal{Z}_2 ^T\, \mathcal{Q}-\Delta _2-\mathcal{F}_4^T\, \mathcal{Q} \approx 0 \end{array} \right. \end{equation}$

(5.13)

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:

$\begin{equation} \begin{aligned} & {_0^C}\mathsf{D}_x^{0.7} y_0(x)+{_0^C}\mathsf{D}_x^{0.5} y_0(x)-2\, y_0(x)-\int _0^x \bigg(x\, {_0^C}\mathsf{D}_t^{0.3} y_0(t)+(x-t)\, y_1(t)\bigg)\, dt -\mathsf{f}_1(x) = 0, \\ & {_0^C}\mathsf{D}_x^{0.6} y_1(x)+\exp(x)\, {_0^C}\mathsf{D}_x^{0.4} y_1(x)+x\, y_1(x)-\int _0^x (x-2t)\, {_0^C}\mathsf{D}_t^{0.8} y_0(t)\, dt -\mathsf{f}_2(x) = 0 \end{aligned} \end{equation}$

(5.14)

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:

$\begin{equation} {_0^C}\mathsf{D}_x^{0.8} y_0(x) \approx \mathcal{Y}_0 ^T\, \Lambda ^{\sigma}(x), \quad {_0^C}\mathsf{D}_x^{0.6} y_1(x) \approx \mathcal{Y}_1 ^T\, \Lambda ^{\sigma}(x) \end{equation}$

(5.15)

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

$\begin{equation} \begin{aligned} & y_0(x) \approx \mathcal{Y}_0 ^T\, \mathbf{\Theta}^{(0.8, \sigma)}\, \Lambda^{\sigma}(x)+y_0(0) = \mathcal{W}_1^T \, \Lambda ^{\sigma}(x), \, \, \mathcal{W}_1 = {\mathbf{\Theta}^{(0.8, \sigma)}}^T\, \mathcal{Y}_0, \\ & y_1(x) \approx \mathcal{Y}_1 ^T\, \mathbf{\Theta}^{(0.6, \sigma)}\, \Lambda^{\sigma}(x)+y_1(0) = \mathcal{Z}_1^T \, \Lambda ^{\sigma}(x), \, \, \mathcal{Z}_1 = {\mathbf{\Theta}^{(0.6, \sigma)}}^T\, \mathcal{Y}_1 \end{aligned} \end{equation}$

(5.16)

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

$\begin{equation} \begin{aligned} &{_0^C}\mathsf{D}_x^{0.8} y_0(x) = {_0^C}\mathsf{D}_x^{0.1}(\mathsf{D}_x^{0.7} y_0(x)) \approx \mathcal{Y}_0 ^T\, \Lambda ^{\sigma}(x), \\ &\quad \Longrightarrow \, \, \mathsf{D}_x^{0.7} y_0(x) \approx \mathcal{Y}_0 ^T\, \mathbf{\Theta}^{(0.1, \sigma)}\, \Lambda(x) = \mathcal{W}_2^T\, \Lambda(x), \quad \mathcal{W}_2 = {\mathbf{\Theta}^{(0.1, \sigma)}}^T\, \mathcal{Y}_0, \\ &{_0^C}\mathsf{D}_x^{0.8} y_0(x) = {_0^C}\mathsf{D}_x^{0.5}(\mathsf{D}_x^{0.3} y_0(x)) \approx \mathcal{Y}_0 ^T\, \Lambda ^{\sigma}(x), \\ & \quad \Longrightarrow \, \, \mathsf{D}_x^{0.3} y_0(x) \approx \mathcal{Y}_0 ^T\, \mathbf{\Theta}^{(0.5, \sigma)}\, \Lambda(x) = \mathcal{W}_3^T\, \Lambda(x), \quad \mathcal{W}_3 = {\mathbf{\Theta}^{(0.5, \sigma)}}^T\, \mathcal{Y}_0, \\ & {_0^C}\mathsf{D}_x^{0.8} y_0(x) = {_0^C}\mathsf{D}_x^{0.3}(\mathsf{D}_x^{0.5} y_0(x)) \approx \mathcal{Y}_0 ^T\, \Lambda ^{\sigma}(x), \\ &\quad \Longrightarrow \, \, \mathsf{D}_x^{0.5} y_0(x) \approx \mathcal{Y}_0 ^T\, \mathbf{\Theta}^{(0.3, \sigma)}\, \Lambda(x) = \mathcal{W}_4^T\, \Lambda(x), \quad \mathcal{W}_4 = {\mathbf{\Theta}^{(0.3, \sigma)}}^T\, \mathcal{Y}_0, \\ & {_0^C}\mathsf{D}_x^{0.6} y_1(x) = {_0^C}\mathsf{D}_x^{0.2}(\mathsf{D}_x^{0.4} y_1(x)) \approx \mathcal{Y}_1 ^T\, \Lambda ^{\sigma}(x), \\ & \quad \Longrightarrow \, \, \mathsf{D}_x^{0.4} y_1(x) \approx \mathcal{Y}_1 ^T\, \mathbf{\Theta}^{(0.2, \sigma)}\, \Lambda(x) = \mathcal{Z}_2 ^T\, \Lambda(x), \quad \mathcal{Z}_2 = {\mathbf{\Theta}^{(0.2, \sigma)}}^T\, \mathcal{Y}_1 \end{aligned} \end{equation}$

(5.17)

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

$\begin{equation} \begin{aligned} & 2 \approx \mathcal{F}_1\, \Lambda ^{\sigma}(x), \quad \exp(x) \approx\mathcal{F}_2\, \Lambda ^{\sigma}(x), \quad x \approx\mathcal{F}_3\, \Lambda ^{\sigma}(x), \\ & 2 \, y_0(x) \approx \mathcal{F}_1 \, \Lambda ^{\sigma}(x) \, {\Lambda ^{\sigma} }^T(x)\, \mathcal{W}_1 \approx \mathcal{F}_1 \, \mathcal{W}_1^*\, \Lambda ^{\sigma}(x), \\ & \exp(x)\, \mathsf{D}_x^{0.4} y_1(x) \approx \mathcal{F}_2 \, \Lambda ^{\sigma}(x) \, {\Lambda ^{\sigma} }^T(x)\, \mathcal{Z}_2 \approx \mathcal{F}_2 \, \mathcal{Z}_2^*\, \Lambda ^{\sigma}(x), \\ & x \, y_1(x) \approx \mathcal{F}_3 \, \Lambda ^{\sigma}(x) \, {\Lambda ^{\sigma} }^T(x)\, \mathcal{Z}_1 \approx \mathcal{F}_3 \, \mathcal{Z}_1^*\, \Lambda ^{\sigma}(x) \end{aligned} \end{equation}$

(5.18)

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:

$\begin{equation} \begin{aligned} & x \approx {\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_1\, \Lambda ^{\sigma}(t), \quad x-t \approx {\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_2\, \Lambda ^{\sigma}(t), \quad x-2t \approx {\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_3\, \Lambda ^{\sigma}(t), \\ & \int _0^x x\, \mathsf{D}_t^{0.3} y_0(t)\, dt \approx {\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_1\, \mathcal{W}_3^*\, \int _0^x \Lambda ^{\sigma}(t) dt \approx {\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_1\, \mathcal{W}_3^*\, \mathbf{\Theta }^{(\sigma)}\, \Lambda ^{\sigma}(x), \\ & \int _0^x (x-t)\, y_1(t)\, dt \approx {\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_2\, \mathcal{Z}_1^*\, \int _0^x \Lambda ^{\sigma}(t) dt \approx {\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_2\, \mathcal{Z}_1^*\, \mathbf{\Theta} ^{(\sigma)}\, \Lambda ^{\sigma}(x), \\ & \int _0^x (x-2t)\, \mathsf{D}_t^{0.8} y_0(t)\, dt \approx {\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_3\, \mathcal{Y}_0^*\, \int _0^x \Lambda ^{\sigma}(t) dt \approx {\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_3\, \mathcal{Y}_0^*\, \mathbf{\Theta }^{(\sigma)}\, \Lambda ^{\sigma}(x) \end{aligned} \end{equation}$

(5.19)

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:

$\begin{equation} \mathsf{f}_1(x) \approx \mathcal{F}_4^T\, \Lambda ^{\sigma}(x), \quad \mathsf{f}_2(x) \approx \mathcal{F}_5^T\, \Lambda ^{\sigma}(x) \end{equation}$

(5.20)

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:

$\begin{equation} \begin{aligned} &\mathcal{R}_1(x) = \mathcal{W}_2 ^T\, \Lambda ^{\sigma}(x)+ \mathcal{W}_4 ^T\, \Lambda ^{\sigma}(x)-\mathcal{F}_1^T \, \mathcal{W}_1^* \, \Lambda ^{\sigma}(x)-{\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_1\, \mathcal{W}_3^*\, \mathbf{\Theta}^{(\sigma)} \, \Lambda ^{\sigma}(x)\\ & \quad -{\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_2\, \mathcal{Z}_1^*\, \mathbf{\Theta}^{(\sigma)} \, \Lambda ^{\sigma}(x)-\mathcal{F}_4^T\, \Lambda ^{\sigma}(x) \approx 0, \\ & \mathcal{R}_2(x) = \mathcal{Y}_1 ^T\, \Lambda ^{\sigma}(x)+ \mathcal{F}_2^T \, \mathcal{Z}_2^* \, \Lambda ^{\sigma}(x)+\mathcal{F}_3^T \, \mathcal{Z}_1^* \, \Lambda ^{\sigma}(x) -{\Lambda ^{\sigma}}^T(x)\, \mathcal{K}_3\, \mathcal{Y}_0^*\, \mathbf{\Theta}^{(\sigma)} \, \Lambda ^{\sigma}(x)-\mathcal{F}_5^T\, \Lambda ^{\sigma}(x) \approx 0 \end{aligned} \end{equation}$

(5.21)

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

$\begin{equation} \left\{ \begin{array}{l} \mathcal{W}_2 ^T\, \mathcal{Q}+ \mathcal{W}_4 ^T\, \mathcal{Q}-\mathcal{F}_1^T \, \mathcal{W}_1^* \, \mathcal{Q}-\Delta _1-\Delta _2-\mathcal{F}_4^T\, \mathcal{Q} \approx 0, \\ \mathcal{Y}_1 ^T\, \mathcal{Q}+\mathcal{F}_2^T\, \mathcal{Z}_2 ^*\, \mathcal{Q}+\mathcal{F}_3^T\, \mathcal{Z}_1 ^*\, \mathcal{Q}-\Delta _3-\mathcal{F}_5^T\, \mathcal{Q} \approx 0 \end{array} \right. \end{equation}$

(5.22)

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:

$\begin{equation*} \mathcal{A}_{\omega ^{\sigma}}^m(\mathbf{I}) = \{ z\, | \, z \text{ is measurable and } \Vert z \Vert _{k, \omega ^{\sigma}} < \infty, k = 0, 1, \cdots, m \} \end{equation*}$

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

$\begin{equation*} \Vert z \Vert _{m, \omega ^{\sigma}} = \bigg ( \sum\limits_{k = 0}^{m} \bigg \Vert \frac{d^k z(x)}{dx^k}\bigg \Vert _{\omega _k^{\sigma}}^2\bigg) ^{\frac{1}{2}}, \quad \vert z \vert _{m, \omega ^{\sigma}} = \bigg \Vert \frac{d^m z(x)}{dx^m}\bigg \Vert _{\omega _m^{\sigma}} \end{equation*}$

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})$

$\begin{equation} \Vert z \Vert _{s, \omega ^{\sigma}} \leq \Vert z \Vert _{[s]+1, \omega ^{\sigma}}^{\iota}\, \Vert z \Vert _{[s], \omega ^{\sigma}}^{1-\iota} \end{equation}$

(6.1)

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]:

$\begin{equation} \langle \phi, \varphi \rangle _{m, \omega ^{\sigma}}\leq \frac{1}{2} \bigg(\Vert \phi \Vert _{m, \omega ^{\sigma}}^2+\Vert \varphi \Vert _{m, \omega ^{\sigma}}^2\bigg) \end{equation}$

(6.2)

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:

$\begin{equation} \Vert y-y_N \Vert _{\eta, \omega ^{\sigma}} \leq \mathcal{C}_0\, (N+2\sigma +1)^{\eta -m}(N+2)^{\frac{\eta -m}{2}} \vert y \vert _{m, \omega ^{\sigma}} \end{equation}$

(6.3)

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

$\begin{equation*} \frac{d^l}{dx^l}(y(x)-y_N(x)) = \frac{d^l}{dx^l} \sum\limits_{k = N+1}^{\infty} g _k^{\sigma}\, \mathcal{G}_k^{\sigma}(x) = \sum\limits_{k = N+1}^{\infty} g _k^{\sigma} \frac{d^l \mathcal{G}_k^{\sigma}(x)}{dx^l}, \quad l\leq m \end{equation*}$

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

$\begin{equation} \begin{aligned} \frac{d^l \mathcal{G}_k^{\sigma}(x)}{dx^l}& = \frac{\Gamma(k+1)\, \Gamma(\sigma +\frac{1}{2})}{\Gamma(k+\sigma +\frac{1}{2})} \frac{d^l \mathcal{J}_k^{\sigma -\frac{1}{2}, \sigma -\frac{1}{2}}(x)}{dx^l}\\ & = \frac{\Gamma(k+1)\, \Gamma(\sigma +\frac{1}{2})\, \Gamma(l+k+2\sigma)}{\Gamma(k+\sigma +\frac{1}{2})\, \Gamma(k+2\sigma)}\mathcal{J}_{k-l}^{l+\sigma -\frac{1}{2}, l+\sigma -\frac{1}{2}}(x) \end{aligned} \end{equation}$

(6.4)

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

$\begin{equation} \bigg \Vert \frac{d^l}{dx^l}(y(x)-y_N(x)) \bigg \Vert _{\omega _l^{\sigma}}^2 = \sum\limits_{k = N+1}^{\infty} {g _k^{\sigma}}^2 \, \frac{\Gamma ^2(k+1)\, \Gamma ^2(\sigma +\frac{1}{2})\, \Gamma ^2 (k+l+2\sigma)}{2\, \Gamma ^2(k+2\sigma)\, \Gamma(l+2\sigma)\, \Gamma(k-l+1)\, (k+\sigma)} \end{equation}$

(6.5)

Similarly, one can obtain

$\begin{equation} \bigg \Vert \frac{d^m y}{dx^m} \bigg \Vert _{\omega _m^{\sigma}}^2 = \sum\limits_{k = m+1}^{\infty} {g _k^{\sigma}}^2 \, \frac{\Gamma ^2(k+1)\, \Gamma ^2(\sigma +\frac{1}{2})\, \Gamma ^2 (k+m+2\sigma)}{2\, \Gamma ^2(k+2\sigma)\, \Gamma(m+2\sigma)\, \Gamma(k-m+1)\, (k+\sigma)} \end{equation}$

(6.6)

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

$\begin{equation} \frac{\Gamma ^2(k+l+2\sigma)\, \Gamma(m+2\sigma)\, \Gamma(k-m+1)}{\Gamma ^2(k+m+2\sigma)\, \Gamma(l+2\sigma)\, \Gamma(k-l+1)}\leq \mathcal{C}_1\, (2\sigma)^{m-l}\, (k+2\sigma)^{2(l-m)}\, (k+1)^{l-m} \end{equation}$

(6.7)

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

$\begin{align*} \vert y-y_N \vert _{l, \omega ^{\sigma}}^2& = \bigg\Vert \frac{d^l}{dx^l}(y-y_N) \bigg \Vert _{\omega _l^{\sigma}}^2\\ & = \sum\limits _{k = N+1}^{\infty} \bigg\{\frac{2 {g_k^{\sigma}}^2\Gamma^2(k+1)\Gamma^2(\sigma +\frac{1}{2})\Gamma^2(k+l+2\sigma)\Gamma^2(k+2\sigma)\Gamma(m+2\sigma)\Gamma(k-m+1)(k+\sigma)}{2{g_k^{\sigma}}^2\Gamma(l+2\sigma)\Gamma(k-l+1)\Gamma^2(k+1)\Gamma^2(\sigma+\frac{1}{2})\Gamma^2(k+m+2\sigma)(k+\sigma)}\\ & \quad \quad \times {g_k^{\sigma}}^2\frac{\Gamma^2(k+1)\Gamma^2(\sigma+\frac{1}{2})\Gamma^2(k+m+2\sigma)}{2\Gamma^2(k+2\sigma)\Gamma(m+2\sigma)\Gamma(k-m+1)(k+\sigma)}\bigg\}\\ & \leq \sum\limits_{k = N+1}^{\infty} \mathcal{C}_1(2\sigma)^{m-l} {g_k^{\sigma}}^2 (k+2\sigma)^{2(l-m)} (k+1)^{l-m}\frac{\Gamma^2(k+1) \Gamma^2(\sigma +\frac{1}{2})\Gamma^2(k+m+2\sigma)}{2\Gamma^2(k+2\sigma)\Gamma(m+2\sigma)\Gamma(k-m+1)(k+\sigma)}\\ & \leq \sum\limits_{k = m+1}^{\infty} \mathcal{C}_1(2\sigma)^{m-l} {g_k^{\sigma}}^2 (N+2\sigma +1)^{2(l-m)} (N+2)^{l-m}\\ & \quad \quad \times \frac{\Gamma^2(k+1) \Gamma^2(\sigma +\frac{1}{2})\Gamma ^2(k+m+2\sigma)}{2\Gamma^2(k+2\sigma)\Gamma(m+2\sigma)\Gamma(k-m+1)(k+\sigma)}\\ & \leq \mathcal{C}_1(2\sigma)^{m-l} {g_k^{\sigma}}^2 (N+2\sigma +1)^{2(l-m)} (N+2)^{l-m} \vert y \vert _{m, \omega ^{\sigma}}^2 \end{align*}$

then, the following inequality is acquired

$\begin{equation} \Vert y-y_N \Vert _{m, \omega ^{\sigma}} \leq \mathcal{C}_2 (2\sigma)^{\frac{m-l}{2}}(N+2\sigma +1)^{l-m}(N+2)^{\frac{l-m}{2}} \vert y \vert _{m, \omega ^{\sigma}}, \quad l \leq m \end{equation}$

(6.8)

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

$\begin{align*} \Vert y-y_N \Vert _{\eta, \omega ^{\sigma}} &\leq \Vert y-y_N \Vert _{[\eta]+1, \omega ^{\sigma}}^{\eta _0} \Vert y-y_N \Vert _{[\eta], \omega ^{\sigma}}^{1-\eta _0} \\ & \leq \mathcal{C}_0 (N+2\sigma +1)^{\eta -m}(N+2)^{\frac{\eta -m}{2}} \vert y \vert _{m, \omega ^{\sigma}} \end{align*}$

□

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

$\begin{align*} \frac{d^l}{dx^l}\bigg(\frac{dy(x)}{dx}-\frac{dy_N(x)}{dx}\bigg)& = \frac{d^{l+1}}{dx^{l+1}}(y(x)-y_N(x)) = \sum\limits_{k = N+1}^{\infty} g_{k}^{\sigma}\, \frac{d^{l+1}\mathcal{G}_k^{\sigma}(x)}{dx^{l+1}}\\ & = \sum\limits_{k = N+1}^{\infty} g_{k}^{\sigma} \frac{\Gamma(k+1)\, \Gamma(\sigma +\frac{1}{2})}{\Gamma(k+\sigma +\frac{1}{2})} \frac{d^{l+1}}{dx^{l+1}}\mathcal{J}_{k}^{\sigma -\frac{1}{2}, \sigma -\frac{1}{2}}(x)\\ & = \sum\limits_{k = N+1}^{\infty} g_{k}^{\sigma} \frac{\Gamma(k+1)\, \Gamma(\sigma +\frac{1}{2})\, \Gamma(l+k+2\sigma +1)}{\Gamma(k+\sigma +\frac{1}{2})\, \Gamma(k+2\sigma)} \mathcal{J}_{k-l-1}^{\sigma +l+\frac{1}{2}, \sigma +l+\frac{1}{2}}(x) \end{align*}$

Therefore, one has

$\begin{equation} \bigg \Vert \frac{dy}{dx}- \frac{dy_N}{dx}\bigg\Vert _{\eta, \omega _{1}} \leq \mathcal{C}_1\, (N+2\sigma +2)^{\frac{\eta -m}{2}}\, (N+2\sigma +4)^{2(\eta -m)}\, (N+1)^{\frac{\eta -m}{2}}\, \vert y \vert _{m, \omega ^{\sigma}} \quad l\leq m, \end{equation}$

(6.9)

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:

$\begin{equation} \bigg \Vert {_0^C}\mathsf{D}_x^{\mu} y-{_0^C}\mathsf{D}_x^{\mu} y_N \bigg \Vert _{\eta, \omega _1^{\sigma}} \leq \delta _0\, \mathcal{C}_1 \frac{\Gamma(\sigma -\mu +\frac{1}{2})\, \Gamma(\sigma +\frac{1}{2})}{\Gamma(1-\mu)\, \Gamma(2\sigma -\mu +1)}(N+2\sigma +2)^{\frac{\eta -m}{2}} (N+2\sigma +4)^{2(\eta -m)} (N+1)^{\frac{\eta -m}{2}} \vert y \vert _{m, \omega ^{\sigma}} \end{equation}$

(6.10)

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

$\begin{equation} \begin{aligned} {_0^C}\mathsf{D}_x^{\mu} y(x)-{_0^C}\mathsf{D}_x^{\mu} y_N(x)& = {_0^{\mathsf{RL}}}\mathsf{I}_x^{1-\mu}(y^{\prime}(x)-y_N^{\prime}(x))\\ & = \frac{1}{\Gamma(1-\mu)}\int _0^x (x-\xi)^{-\mu} (y^{\prime}(\xi)-y_N^{\prime}(\xi))\, d\xi \\ & = \frac{1}{\Gamma(1-\mu)} \bigg\{ x^{-\mu} * (y^{\prime}(x) -y^{\prime}_N(x)) \bigg\} \end{aligned} \end{equation}$

(6.11)

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:

$\begin{align*} \Vert {_0^C}\mathsf{D}_x^{\mu} y-{_0^C}\mathsf{D}_x^{\mu} y_N \Vert _{\eta, \omega _1^{\sigma}}& \leq \frac{\delta _0}{\Gamma(1-\mu)} \Vert x^{-\mu} \Vert _{L^1_{\omega ^{\sigma}}(\mathbf{I})} \Vert y^{\prime}-y_N^{\prime} \Vert _{\eta, \omega _1^{\sigma}} \\ & = \frac{\delta _0}{\Gamma(1-\mu)}\frac{\Gamma(\sigma -\mu +\frac{1}{2})\, \Gamma(\sigma +\frac{1}{2})}{\Gamma(2\sigma -\mu +1)}\Vert y^{\prime} -y_N^{\prime} \Vert _{\eta, \omega _1^{\sigma}}\\ & \leq \delta _0\, \mathcal{C}_1 \frac{\Gamma(\sigma -\mu +\frac{1}{2})\, \Gamma(\sigma +\frac{1}{2})}{\Gamma(1-\mu)\, \Gamma(2\sigma -\mu +1)}(N+2\sigma +2)^{\frac{\eta -m}{2}} (N+2\sigma +4)^{2(\eta -m)}(N+1)^{\frac{\eta -m}{2}} \vert y \vert_{m, \omega ^{\sigma}} \end{align*}$

□

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

$\begin{equation} \begin{aligned} {_0^C}\mathsf{D}_x^{\mu _{r, n}}y_{r, N}(x)&+\sum\limits _{k = 1}^{n-1} \xi _{r, k}(x)\, {_0^C}\mathsf{D}_x^{\mu _{r, n-k}}y_{r, N}(x)+\xi _{r, n}y_{r, N}(x)\\ &-\sum \limits _{j = 1}^{2}\nu _{r, j} \int _0^x \kappa _{r, j}(x, t)\, {_0^C}\mathsf{D}_t^{\gamma _{r, 2-j}}\, y_{j, N}(t)\, dt -\mathsf{f}_r(x)+\mathcal{R}_r(x) = 0 \end{aligned} \end{equation}$

(6.12)

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:

$\begin{equation} \begin{aligned} \mathcal{R}_r(x) = {_0^C}\mathsf{D}_x^{\mu _{r, n}}e_{r, N}(x)&+\sum\limits _{k = 1}^{n-1} \xi _{r, k}(x)\, {_0^C}\mathsf{D}_x^{\mu _{r, n-k}}e_{r, N}(x)+\xi _{r, n}e_{r, N}(x)\\ &- \sum \limits _{j = 1}^{2}\nu _{r, j} \int _0^x \kappa _{r, j}(x, t)\, {_0^C}\mathsf{D}_t^{\gamma _{r, 2-j}}\, e_{j, N}(t)\, dt \end{aligned} \end{equation}$

(6.13)

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:

$\begin{align*} \langle \mathcal{R}_r(x), e_{r, N}(x) \rangle _{\eta, \omega ^{\sigma}} & = \langle {_0^C}\mathsf{D}_x^{\mu _{r, n}}e_{r, N}(x), e_{r, N}(x) \rangle _{\eta, \omega ^{\sigma}}+\langle\sum\limits _{k = 1}^{n-1} \xi _{r, k}(x)\, {_0^C}\mathsf{D}_x^{\mu _{r, n-k}}e_{r, N}(x), e_{r, N}(x)\rangle _{\eta, \omega ^{\sigma}}\\ &+\langle \xi _{r, n}e_{r, N}(x), e_{r, N}(x) \rangle _{\eta, \omega ^{\sigma}}-\langle \sum \limits _{j = 1}^{2}\nu _{r, j} \int _0^x \kappa _{r, j}(x, t)\, {_0^C}\mathsf{D}_t^{\gamma _{r, 2-j}}\, e_{j, N}(t)\, dt, e_{r, N}(x) \rangle _{\eta, \omega ^{\sigma}}\\ & \leq \frac{1}{2} \Vert {_0^C}\mathsf{D}_x^{\mu _{r, n}}e_{r, N} \Vert _{\eta, \omega _1^{\sigma}}^2 +\frac{1}{2} \Vert e_{r, N} \Vert _{\eta, \omega ^{\sigma}}^2+\frac{1}{2} \sum\limits_{k = 1}^{n-1}\Vert \xi _{r, k}\Vert _{L^2_{\omega^{\sigma}}(\mathbf{I})}^2 \, \Vert {_0^C}\mathsf{D}_x^{\mu _{r, n-k}}e_{r, N} \Vert _{\eta, \omega _1^{\sigma}}^2\\ &\, \, +\frac{1}{2} \Vert e_{r, N} \Vert _{\eta, \omega ^{\sigma}}^2+\frac{1}{2} \Vert \xi _{r, N}\Vert _{L^2_{\omega^{\sigma}}(\mathbf{I})}^2 \, \Vert e_{r, N} \Vert _{\eta, \omega ^{\sigma}}^2+\frac{1}{2} \Vert e_{r, N} \Vert _{\eta, \omega ^{\sigma}}^2 \\ & \, \, +\frac{1}{2} \sum\limits_{j = 1}^2 \vert \nu _{r, j} \vert ^2 \, \Vert \kappa_{r, j} \Vert _{L^2_{ \omega_0^{\sigma} }(\mathbf{I}^2)}^2\, \Vert {_0^C}\mathsf{D}_x^{\gamma _{r, 2-j}}e_{j, N} \Vert _{\eta, \omega _1^{\sigma}}^2+\frac{1}{2} \Vert e_{r, N} \Vert _{\eta, \omega ^{\sigma}}^2 \end{align*}$

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

$\begin{equation} \begin{aligned} \Vert \mathcal{R}_r \Vert _{\eta, \omega ^{\sigma}}^2 &\leq \Vert {_0^C}\mathsf{D}_x^{\mu _{r, n}}e_{r, N} \Vert _{\eta, \omega _1^{\sigma}}^2+\sum\limits_{k = 1}^{n-1}\Vert \xi _{r, k}\Vert _{L^2_{\omega^{\sigma}}(\mathbf{I})}^2 \, \Vert {_0^C}\mathsf{D}_x^{\mu _{r, n-k}}e_{r, N} \Vert _{\eta, \omega _1^{\sigma}}^2+3 \Vert e_{r, N} \Vert _{\eta, \omega ^{\sigma}}^2\\ & \, \, +\sum\limits_{j = 1}^2 \vert \nu _{r, j} \vert ^2 \, \Vert \kappa_{r, j} \Vert _{L^2_{ \omega_0^{\sigma} }(\mathbf{I}^2)}^2\, \Vert {_0^C}\mathsf{D}_x^{\gamma _{r, 2-j}}e_{j, N} \Vert _{\eta, \omega _1^{\sigma}}^2+\Vert \xi _{r, N}\Vert _{L^2_{\omega^{\sigma}}(\mathbf{I})}^2 \, \Vert e_{r, N} \Vert _{\eta, \omega ^{\sigma}}^2\\ & \leq \delta _0^2\, \mathcal{C}_1^2\, \frac{\Gamma ^2(\sigma -\mu _{r, n}+\frac{1}{2})\, \Gamma ^2(\sigma +\frac{1}{2})}{\Gamma ^2(1-\mu _{r, n})\, \Gamma ^2(2\sigma -\mu _{r, n}+1)}\\ & \quad \quad \times (N+2\sigma +2)^{\eta -m}(N+2\sigma +4)^{4(\eta -m)}(N+1)^{\eta -m}\vert y \vert _{m, \omega ^{\sigma}}^2\\ & + \delta _0^2\, \mathcal{C}_1^2\, \sum\limits _{k = 1}^{n-1}\Vert \xi _{r, k} \Vert _{L^2_{\omega ^{\sigma}}(\mathbf{I})} ^2\frac{\Gamma^2(\sigma -\mu _{r, n-k}+\frac{1}{2})\, \Gamma^2(\sigma +\frac{1}{2})}{\Gamma^2(1-\mu _{r, n-k})\, \Gamma^2(2\sigma -\mu _{r, n}+1)}\\ & \quad \quad \times (N+2\sigma +2)^{\eta -m}\, (N+2\sigma +4)^{4(\eta -m)} (N+1)^{\eta -m}\, \vert y \vert_{m, \omega ^{\sigma}}^2\\ & +3 \mathcal{C}_1^2\, (N+2\sigma +1)^{2(\eta-m)} (N+2)^{\eta -m}\vert y \vert _{m, \omega ^{\sigma}}^2\\ & +\delta _0^2\, \mathcal{C}_1^2\, \sum\limits _{j = 1}^{2} \vert \nu _{r, j} \vert ^2 \, \Vert \kappa_{r, j} \Vert _{L^2_{\omega _0^{\sigma}}(\mathbf{I}^2)}^2 \frac{\Gamma ^2(\sigma -\gamma_{r, 2-j}+\frac{1}{2})\, \Gamma^2(\sigma +\frac{1}{2})}{\Gamma^2(1-\gamma_{r, 2-j})\, \Gamma^2(2\sigma -\gamma _{r, 2-j})}\\ & \quad \quad \times (N+2\sigma +2)^{\eta -m}(N+2\sigma +4)^{4(\eta -m)} (N+1)^{\eta -m} \vert y \vert _{m, \omega^{\sigma}}\\ & +\mathcal{C}_1^2\, \Vert \mathcal{\xi}_{r, N} \Vert _{L^2_{\omega ^{\sigma}}(\mathbf{I})}^2(N+2\sigma +1)^{2(\eta -m)}(N+2)^{\eta -m} \vert y \vert _{m, \omega ^{\sigma}}^2, \quad r = 1, 2 \end{aligned} \end{equation}$

(6.14)

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:

$\begin{equation*} C_R(y_j) = \frac{|\ln(e_{i+1}^j/e_i^j)|}{\ln(i+1/i)}, \quad e_i^j = |y_j-y_{j, N}|, \quad j = 0, 1, \quad i = 1, 2, \ldots, 8 \end{equation*}$

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

$\begin{align*} & \mathsf{f}_1(x) = -\frac{2}{\Gamma(2.1)}x^{1.1}-\frac{2}{\Gamma(2.4)}x^{1.4}-2x^3+\frac{10}{\Gamma(3.7)} \exp(x)\, x^{2.7}, \\ & \mathsf{f}_2(x) = -\frac{6}{\Gamma(3.2)}x^{2.2}-\frac{6}{\Gamma(3.5)}x^{3.5}+\frac{2}{\Gamma(3.2)}\sin(x) \, x^{2.2}-\frac{2}{3.2 \Gamma(2.2)} x^{3.2} \end{align*}$

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 . As expected, increasing the number of Gegenbauer polynomials in the series solutions (represented by $N$ ) leads to a reduction in errors. reports the MAEs of the approximate solutions for $N = 5$ and different values of $\sigma$ .

Table 1. Values of exact and approximate solutions at selected points for

$N = 7$ and

$\sigma = 1$ for Example 1.

	Exact values		Proposed scheme		Method in ^[37]
$x_i$	${y_0}_{Exact}$	${y_1}_{Exact}$	$y_{0, N}$	$y_{1, N}$	$y_{0, N}$	$y_{1, N}$
$0.0$	$0.0000$	$1.0000$	$3.2144\times10^{-5}$	$0.9999946673$	$0.0000$	$1.0000$
$0.2$	$-0.0400$	$0.9920$	$-0.0399856254$	$0.9919993594$	$-0.04387683293$	$0.99071904024$
$0.4$	$-0.1600$	$0.9360$	$-0.1599851160$	$0.9360018667$	$-0.1667659284$	$0.93075651777$
$0.6$	$-0.3600$	$0.7840$	$-0.3599697321$	$0.7840028442$	$-0.36880182427$	$0.77152149574$
$0.8$	$-0.6400$	$0.4880$	$-0.6399544939$	$0.4880059282$	$-0.64932899451$	$0.46536622778$
$1.0$	$-1.0000$	$0.0000$	$-0.9999322792$	$1.2027 \times 10^{-5}$	$-1.00653583547$	$-0.0349831482$
LSE			$3.7817\times 10^{-5}$	$5.0342\times 10^{-6}$	$3.3483\times 10^{-4}$	$3.1115\times 10^{-3}$

| Show Table

DownLoad: CSV

Table 2. MAEs and convergence rate of approximate solutions for various values of

$N$ and

$\sigma = 1$ of Example 1.

$N$	MAE $(y_0)$	$C_R(y_0)$	MAE $(y_1)$	$C_R(y_1)$
$2$	$2.2605 \times 10^{-2}$	$--$	$6.6949\times 10^{-2}$	$--$
$3$	$2.3660 \times 10^{-3}$	$5.5664$	$4.3205\times 10^{-4}$	$12.4379$
$4$	$9.4636 \times 10^{-4}$	$3.1852$	$9.8083\times 10^{-5}$	$5.1540$
$5$	$2.8943 \times 10^{-4}$	$5.3092$	$4.9989\times 10^{-5}$	$3.0205$
$6$	$1.9745 \times 10^{-4}$	$2.0975$	$1.3829\times 10^{-5}$	$7.0482$
$7$	$7.7834 \times 10^{-5}$	$6.0389$	$1.2027\times 10^{-5}$	$0.9057$
$8$	$4.2115 \times 10^{-5}$	$4.5995$	$8.6863\times 10^{-6}$	$2.4369$

| Show Table

DownLoad: CSV

Table 3. Converegence rates at

$x = 0.5, 0.7$ for various values of

$N$ and

$\sigma = 1$ of Example 1.

	$x=0.5$				$x=0.7$
$N$	Error $(y_0)$	$C_R(y_0)$	Error $(y_1)$	$C_R(y_1)$	Error $(y_0)$	$C_R(y_0)$	Error $(y_1)$	$C_R(y_1)$
$2$	$4.9384 \times 10^{-3}$	$--$	$1.8672\times 10^{-4}$	$--$	$9.6234 \times 10^{-3}$	$--$	$1.6530 \times 10^{-2}$	$--$
$4$	$2.2761 \times 10^{-4}$	$4.4394$	$6.3837\times 10^{-7}$	$8.1923$	$3.3882 \times 10^{-4}$	$4.8280$	$3.8006 \times 10^{-5}$	$8.7646$
$6$	$3.5898 \times 10^{-5}$	$4.5551$	$5.0491\times 10^{-6}$	$5.1004$	$7.2318 \times 10^{-5}$	$3.8089$	$6.4580 \times 10^{-6}$	$4.3713$
$8$	$4.6590 \times 10^{-5}$	$0.9062$	$4.3682\times 10^{-8}$	$16.5114$	$5.4926 \times 10^{-5}$	$0.9562$	$1.4408 \times 10^{-6}$	$5.2145$
$10$	$3.7557 \times 10^{-4}$	$9.3530$	$1.9654\times 10^{-5}$	$27.3774$	$7.2617 \times 10^{-5}$	$1.2513$	$2.9547 \times 10^{-5}$	$13.5374$

| Show Table

DownLoad: CSV

Table 4. LSEs of approximate solutions for various values of

$\sigma$ and

$N = 5$ in Example 1.

$\sigma$	$0.5$	$1$	$1.25$	$1.5$	$2$	$2.5$
LSE $(y_0)$	$2.6457 \times 10^{-5}$	$1.4897 \times 10^{-4}$	$2.3482\times 10^{-4}$	$3.2302 \times 10^{-4}$	$4.9850\times 10^{-4}$	$6.6611\times 10^{-4}$
LSE $(y_1)$	$8.6608 \times 10^{-5}$	$2.0616 \times 10^{-5}$	$2.8292\times 10^{-5}$	$3.6091 \times 10^{-5}$	$5.1284\times 10^{-5}$	$6.5460\times 10^{-5}$

| Show Table

DownLoad: CSV

In , 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, 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 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.

Figure 1. Plots of (a) Exact solution

$y_0(x)$ , (b) Approximate solution

$y_{0, N}(x)$ , (c) Exact solution

$y_1(x)$ , (d) Approximate solution

$y_{1, N}(x)$ , (e) Absolute error function

$\vert y_0(x)-y_{0, N}(x) \vert$ , (f) Absolute error function

$\vert y_1(x)-y_{1, N}(x) \vert$ , for

$N = 7$ and

$\sigma = 1$ for Example 1.

DownLoad: Full-Size Img PowerPoint

Figure 2. Plots of absolute error functions for

$N = 5$ and

$\sigma = 0.5, 1.5, 2.5$ for Example 1.

DownLoad: Full-Size Img PowerPoint

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

$\begin{align*} & \mathsf{f}_1(x) = \frac{\Gamma(5)}{\Gamma(4.3)}x^{3.3}+\frac{\Gamma(5)}{\Gamma(3.5)}x^{4.5}-\frac{23}{12}x^4-\frac{1}{6} x^{3}-\frac{\Gamma(5)}{\Gamma(5.7)}x^{5.7}, \\ & \mathsf{f}_2(x) = \frac{1}{\Gamma(1.4)}x^{0.4}-\frac{2}{\Gamma(2.4)}x^{1.4}+\frac{1}{\Gamma(1.6)}\, x^{0.6}\exp(x)-\frac{2}{ \Gamma(2.6)} x^{1.6}\exp(x)\\ & \quad +x^2-x^3+\bigg ( \frac{2\Gamma(5)}{5.2 \Gamma(4.2)}- \frac{5}{\Gamma(5.2)}\bigg) x^{5.2} \end{align*}$

A comparison between exact and approximate solutions at equally spaced points $x_i = 0.2 i, i = 0, 1, \cdots, 5$ is seen in 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 . In , 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 for $N = 7$ and $\sigma = 1$ . Plots of absolute error functions are seen in for various values of $\sigma$ and $N = 5$ . The convergence rate of approximate solutions have been computed and listed in 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 , 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}$ .

Table 5. Values of exact and approximate solutions at selected points for

$N = 7$ and

$\sigma = 1$ for Example 2.

	Exact		Proposed scheme		Method in ^[37]
$x_i$	${y_0}_{Exact}$	${y_1}_{Exact}$	$y_{0, N}$	$y_{1, N}$	$y_{0, N}$	$y_{1, N}$
$0.0$	$0.0000$	$0.0000$	$-2.3182\times10^{-6}$	$0.0018380780$	$0.0000$	$0.0000$
$0.2$	$0.0016$	$0.1600$	$0.0016007312$	$0.1603567253$	$0.0017309978$	$0.1600362527$
$0.4$	$0.0256$	$0.2400$	$0.0256076988$	$0.2400674525$	$0.0271385452$	$0.2390378788$
$0.6$	$0.1296$	$0.2400$	$0.1296217631$	$0.2402228353$	$0.1359682357$	$0.2321624668$
$0.8$	$0.4096$	$0.1600$	$0.4096477585$	$0.1599868537$	$0.4262711145$	$0.1302339480$
$1.0$	$1.0000$	$0.0000$	$1.0000904470$	$-0.0009318649$	$1.0325094318$	$-0.0805029267$
LSE			$3.8913\times 10^{-5}$	$6.4800\times 10^{-4}$	$7.1677\times 10^{-3}$	$1.0402\times 10^{-2}$

| Show Table

DownLoad: CSV

Table 6. MAEs of approximate solutions for various values of

$N$ and

$\sigma = 1$ of Example 2.

$N$	MAE $(y_0)$	$C_R(y_0)$	MAE $(y_1)$	$C_R(y_1)$
$2$	$1.3064 \times 10^{-1}$	$--$	$4.6931\times 10^{-3}$	$--$
$3$	$1.9922 \times 10^{-2}$	$4.6382$	$9.3530\times 10^{-3}$	$1.7008$
$4$	$2.0981 \times 10^{-4}$	$15.8278$	$4.8217\times 10^{-3}$	$2.3031$
$5$	$1.9955 \times 10^{-4}$	$0.2247$	$3.3191\times 10^{-3}$	$1.6735$
$6$	$1.2604 \times 10^{-4}$	$2.5201$	$2.4172\times 10^{-3}$	$1.7391$
$7$	$9.0447 \times 10^{-5}$	$2.1527$	$1.8381\times 10^{-3}$	$1.7767$
$8$	$6.5017 \times 10^{-5}$	$2.4742$	$8.4361\times 10^{-4}$	$5.8323$

| Show Table

DownLoad: CSV

Table 7. LSEs of approximate solutions for various values of

$\sigma$ and

$N = 5$ in Example 2.

$\sigma$	$0.5$	$1$	$1.25$	$1.5$	$2$	$2.5$
LSE $(y_0)$	$4.5374 \times 10^{-6}$	$8.3138 \times 10^{-4}$	$1.2031\times 10^{-4}$	$1.5449 \times 10^{-4}$	$2.1416\times 10^{-4}$	$2.6366\times 10^{-4}$
LSE $(y_1)$	$8.1582 \times 10^{-4}$	$1.1886 \times 10^{-3}$	$1.4490\times 10^{-3}$	$1.7085 \times 10^{-3}$	$2.1904\times 10^{-3}$	$2.6115\times 10^{-3}$

| Show Table

DownLoad: CSV

Figure 3. Plots of (a) Exact solution

$y_0(x)$ , (b) Approximate solution

$y_{0, N}(x)$ , (c) Exact solution

$y_1(x)$ , (d) Approximate solution

$y_{1, N}(x)$ , (e) Absolute error function

$y_{0, N}(x)$ , (f) Absolute error function

$y_{1, N}(x)$ , for

$N = 7$ and

$\sigma = 1$ for Example 2.

DownLoad: Full-Size Img PowerPoint

Figure 4. Plots of absolute error functions for

$N = 5$ and

$\sigma = 0.5, 1.25, 1.5$ for Example 2.

DownLoad: Full-Size Img PowerPoint

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 , as well as and . 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.

References

[1]	M. M. Khader, N. H. Sweilam, On the approximate solutions for system of fractional integro-differential equations using Chebyshev pseudo-spectral method, Appl. Math. Model., 37 (2013), 9819–9828. https://doi.org/10.1016/j.apm.2013.06.010 doi: 10.1016/j.apm.2013.06.010
[2]	H. Aminikhah, A new analytical method for solving systems of linear integro-differential equations, J. King Saud Univ. Sci., 23 (2011), 349–353. https://doi.org/10.1016/j.jksus.2010.07.016 doi: 10.1016/j.jksus.2010.07.016
[3]	R. C. G. Sekar, K. Murugesan, A new analytical method for solving systems of linear integro-differential equations, Appl. Math. Comput., 273 (2016), 484–492. https://doi.org/10.1016/j.amc.2015.09.092 doi: 10.1016/j.amc.2015.09.092
[4]	F. Saedpanah, Well-posedness of an integro-differential equation with positive type kernels modeling fractional order viscoelasticity, Eur. J. Mech. A-Solid., 44 (2014), 201–211. https://doi.org/10.1016/j.euromechsol.2013.10.014 doi: 10.1016/j.euromechsol.2013.10.014
[5]	S. Larson, M. Racheva, F. Saedpanah, Discontinuous Galerkin method for an integro-differential equation modeling dynamic fractional order viscoelasticity, Comput. Method. Appl. M., 283 (2015), 196–209. https://doi.org/10.1016/j.cma.2014.09.018 doi: 10.1016/j.cma.2014.09.018
[6]	J. Wang, T. Z. Xu, Y. Q. Wei, J. Q. Xie, Numerical simulation for coupled systems of nonlinear fractional order integro-differential equations via wavelets method, Appl. Math. Comput., 324 (2018), 36–50. https://doi.org/10.1016/j.amc.2017.12.010 doi: 10.1016/j.amc.2017.12.010
[7]	S. A. Deif, S. R. Grace, Iterative refinement for a system of linear integro-differential equations of fractional type, J. Comput. Appl. Math., 294 (2016), 138–150. https://doi.org/10.1016/j.cam.2015.08.008 doi: 10.1016/j.cam.2015.08.008
[8]	F. Saemi, H. Ebrahimi, M. Shafiee, An effective scheme for solving system of fractional Volterra-Fredholm integro-differential equations based on the Muntz-Legendre wavelets, J. Comput. Appl. Math., 374 (2020), 112773. https://doi.org/10.1016/j.cam.2020.112773 doi: 10.1016/j.cam.2020.112773
[9]	M. H. Heydari, M. R. Hooshmandasl, F. Mohammadi, C. Cattani, Wavelets method for solving systems of nonlinear singular fractional Volterra integro-differential equations, Commun. Nonlinear Sci., 19 (2014), 37–48. https://doi.org/10.1016/j.cnsns.2013.04.026 doi: 10.1016/j.cnsns.2013.04.026
[10]	O. H. Mohammed, A. M. Malik, A modified computational algorithm for solving systems of linear integro-differential equations of fractional order, J. King Saud Univ. Sci., 31 (2019), 946–955. https://doi.org/10.1016/j.jksus.2018.09.005 doi: 10.1016/j.jksus.2018.09.005
[11]	J. R. Loh, C. Phang, A new numerical scheme for solving system of Volterra integro-differential equation, Alex. Eng. J., 57 (2018), 1117–1124. https://doi.org/10.1016/j.aej.2017.01.021 doi: 10.1016/j.aej.2017.01.021
[12]	F. Youbi, S. Momani, S. Hasan, M. Al-Smadi, Effective numerical technique for nonlinear Caputo-Fabrizio systems of fractional Volterra integro-differential equations in Hilbert space, Alex. Eng. J., 61 (2022), 1778–1786. https://doi.org/10.1016/j.aej.2021.06.086 doi: 10.1016/j.aej.2021.06.086
[13]	M. Akbar, R. Nawaz, S. Ahsan, K. S. Nisar, A. H. Abdel-Aty, H. Eleuch, New approach to approximate the solution for the system of fractional order Volterra integro-differential equations, Results Phys., 19 (2020), 103453. https://doi.org/10.1016/j.rinp.2020.103453 doi: 10.1016/j.rinp.2020.103453
[14]	D. Wang, Y. Liu, H. Li, Z. Fang, Second-order time stepping scheme combined with a mixed element method for a 2D nonlinear fourth-order fractional integro-differential equations, Fractal Fract., 6 (2022), 201. https://doi.org/10.3390/fractalfract6040201 doi: 10.3390/fractalfract6040201
[15]	H. R. Marzban, A generalization of Muntz-Legendre polynomials and its implementation in optimal control of nonlinear fractional delay systems, Chaos Soliton. Fract., 158 (2022), 112093. https://doi.org/10.1016/j.chaos.2022.112093 doi: 10.1016/j.chaos.2022.112093
[16]	M. Pourbabaee, A. Saadatmandi, A new operational matrix based on Müntz-Legendre polynomials for solving distributed order fractional differential equations, Math. Comput. Simulat., 194 (2022), 210–235. https://doi.org/10.1016/j.matcom.2021.11.023 doi: 10.1016/j.matcom.2021.11.023
[17]	A. A. El-Sayed, P. Agarwal, Spectral treatment for the fractional-order wave equation using shifted Chebyshev orthogonal polynomials, J. Comput. Appl. Math., 424 (2023), 114933. https://doi.org/10.1016/j.cam.2022.114933 doi: 10.1016/j.cam.2022.114933
[18]	I. Terghini, A. Hasseine, D. Caccavo, H. J. Bart, Solution of the population balance equation for wet granulation using second kind Chebyshev polynomials, Chem. Eng. Res. Des., 189 (2023), 262–271. https://doi.org/10.1016/j.cherd.2022.11.028 doi: 10.1016/j.cherd.2022.11.028
[19]	D. S. Mohamed, R. A. Taher, Comparison of Chebyshev and Legendre polynomials methods for solving two dimensional Volterra-Fredholm integral equations, Journal of the Egyptian Mathematical Society, 25 (2017), 302–307. https://doi.org/10.1016/j.joems.2017.03.002 doi: 10.1016/j.joems.2017.03.002
[20]	N. H. Sweilam, A. M. Nagy, A. A. El-Sayed, Second kind shifted Chebyshev polynomials for solving space fractional order diffusion equation, Chaos Soliton. Fract., 73 (2015), 141–147. https://doi.org/10.1016/j.chaos.2015.01.010 doi: 10.1016/j.chaos.2015.01.010
[21]	W. Zheng, Y. Chen, J. Zhou, A Legendre spectral method for multidimensional partial Volterra integro-differential equations, J. Comput. Appl. Math., 436 (2024), 115302. https://doi.org/10.1016/j.cam.2023.115302 doi: 10.1016/j.cam.2023.115302
[22]	M. Usman, M. Hamid, T. Zubair, R. U. Haq, W. Wang, M. B. Liu, Novel operational matrices-based method for solving fractional-order delay differential equations via shifted Gegenbauer polynomials, Appl. Math. Comput., 372 (2020), 124985. https://doi.org/10.1016/j.amc.2019.124985 doi: 10.1016/j.amc.2019.124985
[23]	K. M. Hosny, M. M. Darwish, Novel quaternion discrete shifted Gegenbauer moments of fractional-orders for color image analysis, Appl. Math. Comput., 421 (2022), 126926. https://doi.org/10.1016/j.amc.2022.126926 doi: 10.1016/j.amc.2022.126926
[24]	J. H. S. Alkhalissi, I. Emiroglu, M. Bayram, A. Secer, F. Tasci, A new operational matrix of fractional derivative based on the generalized Gegenbauer-Humbert polynomials to solve fractional differential equations, Alex. Eng. J., 60 (2021), 3509–3519. https://doi.org/10.1016/j.aej.2021.02.012 doi: 10.1016/j.aej.2021.02.012
[25]	M. Faheem, A. Khan, A wavelet collocation method based on Gegenbauer scaling function for solving fourth-order time-fractional integro-differential equations with a weakly singular kernel, Appl. Numer. Math., 184 (2023), 197–218. https://doi.org/10.1016/j.apnum.2022.10.003 doi: 10.1016/j.apnum.2022.10.003
[26]	G. Ozaltun, A. Konuralp, S. Gumgum, Gegenbauer wavelet solutions of fractional integro-differential equations, J. Comput. Appl. Math., 420 (2023), 114830. https://doi.org/10.1016/j.cam.2022.114830 doi: 10.1016/j.cam.2022.114830
[27]	T. D. Phung, D. T. Duc, V. K. Tuan, Multi-term fractional oscillation integro-differential equations, Fract. Calc. Appl. Anal., 25 (2022), 1713–1733. https://doi.org/10.1007/s13540-022-00074-8 doi: 10.1007/s13540-022-00074-8
[28]	P. J. Torvik, R. L. Bagley, On the appearance of the fractional derivative in the behavior of real materials, J. Appl. Mech., 51 (1984), 294–298. https://doi.org/10.1115/1.3167615 doi: 10.1115/1.3167615
[29]	Y. A. Rossikhin, M. V. Shitikova, Application of fractional calculus for dynamic problems of solid mechanics: novel trends and recent results, Appl. Mech. Rev., 63 (2010), 010801. https://doi.org/10.1115/1.4000563 doi: 10.1115/1.4000563
[30]	I. Podlubny, Fractional differential equations, San Diego: Academic Press, 1998.
[31]	Q. Lei, A. Chen, Arzela-Ascoli theorem for demi-linear mappings, J. Funct. Space., 2014 (2014), 679825. https://doi.org/10.1155/2014/679825 doi: 10.1155/2014/679825
[32]	Y. Liu, Z. Li, Krasnoselskii type fixed point theorem and applications, P. Am. Math. Soc., 136 (2008), 1213–1220. https://doi.org/10.1090/S0002-9939-07-09190-3 doi: 10.1090/S0002-9939-07-09190-3
[33]	K. Sadri, K. Hosseini, D. Baleanu, S. Salahshour, C. Park, Designing a matrix collocation method for fractional delay integro-differential equations with weakly singular kernels based on Vieta-Fibonacci polynomials, Fractal Fract., 6 (2022), 2. https://doi.org/10.3390/fractalfract6010002 doi: 10.3390/fractalfract6010002
[34]	E. H. Doha, A. H. Bhrawy, S. S. Ezz-Eldien, A new Jacobi operational matrix: An application for solving fractional differential equations, Appl. Math. Model., 36 (2012), 4931–4943. https://doi.org/10.1016/j.apm.2011.12.031 doi: 10.1016/j.apm.2011.12.031
[35]	B. Y. Guo, L. L. Wang, S. S. Ezz-Eldien, Jacobi approximations in non-uniformly Jacobi-weighted Sobolev spaces, J. Approx. Theory, 128 (2004), 1–41. https://doi.org/10.1016/j.jat.2004.03.008 doi: 10.1016/j.jat.2004.03.008
[36]	E. Kreyszig, Introductory functional analysis with applications, Canada: John Wiley & Sons. Inc., 1978.
[37]	S. S. Ahmed, S. A. Hamasalih, Solving a system of Caputo fractional-order Volterra integro-differential equations with variable coefficients based on the finite difference approximation via the block-by-block method, Symmetry, 15 (2023), 607. https://doi.org/10.3390/sym15030607 doi: 10.3390/sym15030607

This article has been cited by:

1.	Ravikiran A. Mundewadi, Raju B. Jummannaver, Hosoya polynomial method for the numerical solution of Volterra integral equations, 2024, 2731-6734, 10.1007/s43994-024-00191-5
2.	Khadijeh Sadri, David Amilo, Muhammad Farman, Evren Hinçal, Bivariate Jacobi polynomials depending on four parameters and their effect on solutions of time-fractional Burgers’ equations, 2024, 83, 18777503, 102450, 10.1016/j.jocs.2024.102450
3.	Aly R. Seadawy, Bayan A. Alsaedi, Soliton solutions of nonlinear Schrödinger dynamical equation with exotic law nonlinearity by variational principle method, 2024, 56, 1572-817X, 10.1007/s11082-024-06367-x
4.	Syed T. R. Rizvi, Aly R. Seadawy, Bazgha Mustafa, Discussion on chirped periodic and optical soliton solutions for chiral nonlinear Schrödinger dynamical equation with Bohm potential, 2024, 56, 1572-817X, 10.1007/s11082-024-06658-3
5.	Hoorieh Fakhari, Akbar Mohebbi, Galerkin spectral and finite difference methods for the solution of fourth-order time fractional partial integro-differential equation with a weakly singular kernel, 2024, 70, 1598-5865, 5063, 10.1007/s12190-024-02173-6
6.	Waleed Mohamed Abd-Elhameed, Omar Mazen Alqubori, Ahmed Gamal Atta, A collocation procedure for the numerical treatment of FitzHugh–Nagumo equation using a kind of Chebyshev polynomials, 2025, 10, 2473-6988, 1201, 10.3934/math.2025057
7.	Saman Bagherbana, Jafar Biazar, Hossein Aminikhah, Convergence analysis of a novel fractional product integration–spectral collocation scheme for nonlinear weakly singular Volterra integro–differential systems, 2025, 0020-7160, 1, 10.1080/00207160.2025.2465770
8.	Ramy Mahmoud Hafez, Hany Mostafa Ahmed, Omar Mazen Alqubori, Amr Kamel Amin, Waleed Mohamed Abd-Elhameed, Efficient Spectral Galerkin and Collocation Approaches Using Telephone Polynomials for Solving Some Models of Differential Equations with Convergence Analysis, 2025, 13, 2227-7390, 918, 10.3390/math13060918
9.	Waleed Mohamed Abd-Elhameed, Omar Mazen Alqubori, Naher Mohammed A. Alsafri, Amr Kamel Amin, Ahmed Gamal Atta, A Matrix Approach by Convolved Fermat Polynomials for Solving the Fractional Burgers’ Equation, 2025, 13, 2227-7390, 1135, 10.3390/math13071135
10.	David Amilo, Khadijeh Sadri, Evren Hincal, Comparative Analysis of Caputo and Variable-Order Fractional Derivative Algorithms Across Various Applications, 2025, 11, 2349-5103, 10.1007/s40819-025-01887-w

Reader Comments

Your name:*

Email:*
© 2024 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

AIMS Mathematics

1.8 3.4

Metrics

Article views(1417) PDF downloads(78) Cited by(10)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Figures and Tables

Figures(4) / Tables(7)

AIMS Mathematics

A tau-Gegenbauer spectral approach for systems of fractional integro-differential equations with the error analysis

Related Papers:

Abstract

1. Introduction

2. Preliminaries

3. Existence and uniqueness

4. Shifted Gegenbauer polynomials and their operational matrices

4.1. Shifted Gegenbauer polynomials

4.2. Operational matrices

4.2.1. Operational matrix of the integration with the integer order

4.2.2. Operational matrix of the integration with the fractional order

4.2.3. Operational matrix of the product

5. Solution process

5.1. System I

5.2. System II

6. Error bounds

7. Numerical examples

8. Conclusions

Use of AI tools declaration

Acknowledgments

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Catalog

AIMS Mathematics

A tau-Gegenbauer spectral approach for systems of fractional integro-differential equations with the error analysis

Related Papers:

Abstract

1. Introduction

2. Preliminaries

3. Existence and uniqueness

4. Shifted Gegenbauer polynomials and their operational matrices

4.1. Shifted Gegenbauer polynomials

4.2. Operational matrices

4.2.1. Operational matrix of the integration with the integer order

4.2.2. Operational matrix of the integration with the fractional order

4.2.3. Operational matrix of the product

5. Solution process

5.1. System I

5.2. System II

6. Error bounds

7. Numerical examples

8. Conclusions

Use of AI tools declaration

Acknowledgments

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog