Legendre spectral collocation method for solving nonlinear fractional Fredholm integro-differential equations with convergence analysis

1.   Introduction

For the last three decades, researchers have been drawn to fractional calculus, which has numerous applications in engineering and physics. Fractional calculus finds utility in various fields like image processing ^[1], dynamical system control theory, signal processing, electrical networks, optics, probability, statistics ^[2], and chemical physics. Fractional integro-differential equations, both linear and nonlinear, play a pivotal role in natural sciences and engineering, solving mathematical modeling challenges in spatiotemporal developments, biological and physical problems, epidemic modeling, and boundary value problems^[3,4,5]. Given that analytical solutions are often elusive, the development of numerical methods for approximating solutions has become imperative.

The author in ^[6] has provided a solution to the dynamic model of atmospheric $CO_{2}$ concentration via a fractional mathematical model of the non-linear nature, but ^[7] computed the solution of some Cauchy problems and diffusion equations modeled with the Hilfer-Prabhakar fractional derivative via the Kharrat-Toma transform. The authors in ^[8], investigated a fractional extension of Lienard's equation by using a fractional operator with an exponential kernel.

In this paper, we present the Legendre-Gauss-Lobatto collocation (L-G-LC) method to solve the non-linear fractional Fredholm integro-differential equations (non-FFIDEs), given by:

with the initial conditions;

where $D^{\alpha_{1}}$ denotes the fractional derivative of order $\alpha_{1}$, and $0<\alpha_{1}<2$.

Physical processes such as neutron transport ^[9], neural networks ^[10], population models ^[11], filtering and scattering ^[12], inverse problems ^[13], and disease spread ^[14] are all modeled using Fredholm integro-differential equations (FIDEs). In ^[15], the Legendre collocation method was used for fractional Volterra integro-differential equations. Whereas, in this work, we will solve fractional Fredholm integro-differential equations and provide the convergence analysis. Many papers published in recent years (e.g., ^[16,17]) have been devoted to discussing the FIDEs. A linear fractional Fredholm differential equation with variable coefficients has been approximated via the Taylor matrix method ^[18]. While ^[19] introduces a robust algorithm for mathematically resolving a family of two-dimensional FIDEs, the authors in ^[20] utilized the reproducing kernel method for approximating the solution of a nonlinear FIDE, and the fractional derivative is given in the Caputo sense in this paper. ^[21] used a computational method for solving a class of nonlinear Fredholm integro-differential equations of fractional order based on the second kind of Chebyshev wavelet. ^[22] constructed a method based on Haar wavelet approximation. The authors in ^[23] used an accurate numerical approach based on Wilson wavelets and the collocation method, as well as the Kumar and Sloan scheme, to numerically solve the nonlinear fractional integro-differential equations of Fredholm-Hammerstein. ^[24] used Taylor expansion to solve the linear fractional integro-differential equations of the Fredholm and Volterra types. It entails the $n^{th}$-order Taylor expansion of the unknown function at any point, yielding a system of equations for the unknown function and its m derivatives. ^[25] utilized a novel iterative algorithm for solving Volterra partial integro-differential problems with a weakly singular kernel. The author in ^[26] solves integro-differential equations (IDEs) by utilizing Jacobi-Gauss quadrature, while ^[27] solves fractional differential and integro-differential equations by utilizing an operational matrix method. The authors in ^[28] utilized the shifted Jacobi spectral collocation method to solve IDEs, while in ^[29] the author used the Taylor polynomial for solving non-linear IDEs. ^[30] solved IDEs with weakly singular kernels by utilizing spline collocation, while in ^[31,32] the authors find a solution to FIDEs by applying Legendre multiwavelets and a shifted Lucas polynomial, respectively. There are several other numerical methods to consider ^[33,34]. Differential collocation schemes can be obtained by either directly approximating fractional derivative operators or by recasting the governing differential equation into an equivalent integral equation ^[35]. Because of numerical differentiation, the integral collocation scheme has better stability properties than its differential counterpart, whereas numerical integration is inherently stable ^[36]. However, differential collocation schemes must use efficient integration preconditioners to overcome ill-conditioning issues when solving differential equations, which worsen as the number of collocation nodes increases.

In the last four decades, spectral methods ^{[37,38,39,40,41]} have been widely used in a variety of fields. Initially, Fourier-expanded spectral techniques were used in a few contexts, such as periodic boundary conditions and simple geometric areas. They have recently advanced theoretically and been used as powerful techniques to solve a variety of problems. When compared to other numerical techniques, spectral techniques have a superior character based on thoroughness and exponential averages of convergence. The fundamental step in all spectral techniques is to express the problem solution as a finite series of several functions. Spectral methods include many types, such as collocation ^{[42,43,44,45]}, tau ^[46], Galerkin ^[47], and Petrov-Galerkin ^[48]. The coefficients will then be chosen to minimize the absolute error. During this time, the numerical solution in the spectral collocation technique will be implemented to nearly satisfy IDEs (see ^[49,50,51]). On the other hand, at selection points, the residuals may be permitted to be zero. The collocation approach has been successfully applied in a wide range of scientific and engineering areas due to its obvious advantages. Because their global nature fits well with the nonlocal definition of fractional operators, spectral collocation methods are promising candidates for solving fractional differential equations.

The main goal of the paper is to use the L-G-LC method for approximating non-FFIDEs with L-GL interpolation nodes. We estimate the residuals of the aforementioned problem by using a finite expansion of a Legendre polynomial for independent variables and L-GL quadrature points to approximate the solution of an equation. When these equations are combined with the initial conditions, they yield an algebraic system of (N+1) equations that can be solved. In addition, we investigate the convergence of approximation solutions. To demonstrate the method's accuracy, numerical simulations of some non-FFIDEs are presented.

This paper is structured as follows: Section 2 contains some preliminary information about FDEs. In Section 3, we utilized the L-G-LC method for building a technique to solve nonlinear FFIDEs. Section 4 discusses the spectral collocation method's convergence analysis. Numerical simulations are presented in Section 5 to ensure the effectiveness of the proposed method. Some observations and conclusions are provided in Section 6.

2.   Fractional calculus

This section introduces the main definitions used in the following section as the left and right Caputo (RL-C) definitions. ^[43] gives the RL-C derivative $D^\alpha_{1}$ of order $\alpha_{1}$ as

The operator $D_{\pm}^{\alpha_{1}}$ satisfies the following properties:

where $D^{\alpha_{1}}{\pm}$ and $I{\pm}^{\alpha_{1}}$ are operators of left and right Caputo differential and integral, respectively.

where $\lfloor\alpha_{1} \rfloor$ and $\lceil\alpha_{1} \rceil$ denote the floor and ceiling functions, respectively, while $\mathbb{N} = \{1, 2, \ldots\}$ and $\mathbb{N}_0 = \{0, 1, 2, \ldots\}$.

For $\alpha_{1} > 0$, the fractional integrals of order $\alpha_{1}$ ^[43], both left-sided and right-sided, are specified as follows:

3.   The spectral collocation method

In the next section, we will propose a Legendre spectral collocation (LSC) method for solving Eq (1.1).

3.1.   L-GL interpolation

The Legendre polynomials $\mathcal{L}_{\kappa}(\varrho),$ $\kappa = 0, 1\ldots,$ follow the Rodrigues formula ^[52]:

Furthermore, $\mathcal{L}_{\kappa}(\varrho)$ aligns with a polynomial of degree $\kappa$, resulting in the $p$th derivative of $\mathcal{L}_{\kappa}(\varrho)$ as:

where

Orthogonality is attained through the execution of the subsequent procedures:

where $\chi(\varrho) = 1$, $h_{\kappa} = \frac{2}{2\kappa+1}.$

The L-GL quadrature method was employed to effectively compute the integrals presented earlier. For any $\mathcal{Z}\in S_{2\nu_{1}-1}$ on the interval $[-1, 1]$, we can express:

For instance, contemplate the discrete inner product:

In the case of L-GL, we ascertain that ^[53] $\varrho_{\nu_{1}, 0} = -1$, $\quad \varrho_{\nu_{1}, \nu_{1}} = 1$, $\quad \varrho_{\nu_{1}, \varepsilon}$ and $(\varepsilon = 1, \cdots, \nu_{1}-1)$ are the zeros of $(l_{\nu_{1}}(\varrho))^{'}$, and $\varpi_{\nu_{1}, \varepsilon} = \frac{2}{\nu_{1}(\nu_{1}+1)(\mathcal{L}_{\nu_{1}}(\varrho_{\nu_{1}, \varepsilon}))^{2}},$ where $\varpi_{\nu_{1}, \varepsilon}$ (with $0\leq \varepsilon\leq \nu_{1}$) and $\varrho_{\nu_{1}, \varepsilon}$ (with $0\leq \varepsilon\leq \nu_{1}$) serve as the Christoffel numbers and nodes within the interval $[-1, 1]$, respectively. To apply these polynomials in the range $\varrho\in (0, l_1)$, we introduce a shifted Legendre polynomials (SLPs) by utilizing $\varrho = \dfrac{2\varrho}{l_1}-1$.

If we denote by $\mathcal{L}_{l_1, \rho}{(\varrho)}$ the SLP $\mathcal{L}_\rho{\left(\dfrac{2\varrho}{l_1}-1\right)}$, then $\mathcal{L}_{l_1, \rho}{(\varrho)}$ can be obtained as ^[52]:

The analytical representation of the SLP denoted as $\mathcal{L}_{l_1, \rho}{(\varrho)}$ of degree $\rho$ is expressed as follows:

The condition of orthogonality is expressed as:

where $w_{l_1}(\varrho) = 1\ \text{and}\ \hbar_\kappa = \dfrac{l_1}{2\kappa+1}.$

In terms of SLPs, a square integrable function $\mathcal{Z}(\varrho)$ in the interval $(0, l)$ can be expressed as:

where the coefficients $e_\varepsilon$ are:

The initial $(\nu_{1}+1)$ terms of the SLPs find application in practical scenarios. Consequently, $\mathcal{Z}(\varrho)$ is formulated as follows:

3.2.   LSC scheme

In this section we solve Eq (1.1) by using the following transformations $t=2\lambda-1$, $s=2\varrho-1$, $\mathcal{Y}(2\varrho-1)=\mathcal{Z}(\varrho)$, $\phi(2\varrho-1)=\phi(\varrho)$, $2G(\mathcal{Y}(2\lambda-1)=F(\mathcal{Z}(\lambda))$, and $\eta(2\varrho-1,2\lambda-1)=\sigma(\varrho,\lambda)$, we obtain,

with the initial conditions

We used the L-G-LC method to solve non-FFIDEs with initial conditions. The LSC method for Eq (3.12) is to explore the approximate solution in the form,

As a result, inserting Eq (3.14) into Eq (3.12),

where $I_{\varrho, \nu_{1}}$, $I_{\lambda, \nu_{1}}$ are Legendre-Gauss-Lobatto interpolation operators.

Now, we describe how we implemented our form Eq (3.15) by utilizing the Legendre-Gauss interpolation, which serves as the foundation for our scheme, setting,

we can obtain

where

The fractional derivative of $\mathcal{Z}_{\nu_{1}}(\varrho)$ is then evaluated as

Utilizing Eqs (3.19) and (3.17), we can express Eq (3.12) as:

where $\nu_{1}^{\varepsilon} = \sum_{\varepsilon = 0}^{\nu_{1}}\phi(\varrho)\varpi_{\varepsilon} L_{\varepsilon}(\varrho)$.

Merging Eqs (3.14) and (3.13), we get

Conversely, we can formulate

Equations (3.21)–(3.23) are equivalent to a discretized system of $(\nu_{1}+1)$ algebraic equations with the unknowns $e_{\varepsilon}, \, \rho = 0, \cdots, \nu_{1}$,

At last, the system of $(\nu_{1}+1)$ algebraic equations generated are solved. In our implementation, this system has been solved using the Mathematica function. As a result, $\mathcal{Z}_{\nu{1}}(\varrho)$ can be calculated in closed form.

4.   Convergence analysis

A discussion of error analysis is included, as well as some useful lemmas.

Assume$P_{N} : L^{2}(I)\rightarrow \varrho_N$ to be the $L^2$ orthogonal projection, as defined by ^[53]

Here are some definitions of weighted Hilbert spaces. For a nonnegative integer $\eta$, define ^[52,53]

whereas $\partial^{\imath}_{\varrho}\mathcal{Z}(\varrho) = \frac{\partial^{\imath}\mathcal{Z}(\varrho)}{\partial{\varrho}^{\imath}}$, related to the norm and semi-norm as follows:

Lemma 1. Assume $I_{N}$ denotes the polynomial of degree $N$ interpolating $\mathcal{Z}$ at one of these point sets, and that the interpolation error is estimated as $\mathcal{Z}-I_{N}\mathcal{Z}$ in the norms of the Sobolev spaces $H^{1}(-1, 1)$. When $u\in H^{\eta}(-1, 1)$ and $\eta\geq 1,$ the relationship holds, as indicated by ^[53].

Lemma 2. Let $\mathcal{Z}\in H^{\eta}(I), \quad I\equiv(-1, 1)$. The following estimates are satisfied by the interpolation of $\mathcal{Z}$ ($I_{N}\xi$) computed at any point of Jacobi Gauss points (Gauss-Radau points, Gauss-lobatto, or Gauss) ^[53]:

Lemma 3. s Consider $e(x) = \mathcal{Z}(\varrho)-\mathcal{Z} N(\varrho)$ to represent the error function of the solution. The subsequent inequality is applicable in this context:

where

$B_{1}=I_{\varrho,N}D^{\alpha_{1}}\mathcal{Z}(\varrho)-D^{\alpha_{1}}\mathcal{Z}(\varrho)$

$B_{2} = I_{\varrho, N}\int_{-1}^{1}(I-I_{\lambda, N})\bigg[\sigma(\varrho, \lambda)F(\mathcal{Z}(\lambda))\bigg]d\lambda$

$B_{3} = I_{\varrho, N}\int\limits_{-1}^{1}I_{\lambda, N}\bigg[\sigma(\varrho, \lambda)F(\mathcal{Z}(\lambda))-\sigma(\varrho, \lambda)F(\mathcal{Z}_{N}(\lambda))\bigg]d\lambda$.

Proof. By using the Caputo definition, we write the equation of non-FFIDEs as follows:

and when utilizing the approximate solution we have,

Subtracting (4.5) from (4.4) yields

hence

The desired result can be obtained directly from the above.

4.1.   Error analysis

Theorem 1. Let $I_{N}\mathcal{Z}(\varpi)$ be the spectral approximate and let $\mathcal{Z}(\varpi)$ be the exact solution of the equation of non-FFIDEs and, consequently we obtain

Proof. In the $\varrho$-direction, the interpolation operator is formally defined as $I_{\varrho, N}:C(-1, 1)\rightarrow P_{N}$ for any $\mathcal{Z}(\varrho)\in C(-1, 1)$

The Jacobi polynomial is reduced to the Legender polynomial $L_{n}(\varrho)$ in the special case on condition that $\nu_{1} = \nu = 0$. We can write $\varrho_{\ell} = \varrho_{\ell}^{0, 0}$, $\phi_{\ell} = \phi_{\ell}^{0, 0},$ and $I_{\varrho, N} = I_{\varrho, N}^{0, 0}$. Equation (4.10) can be formulated utilizing Lemma 3 and the Gronwall inequality:

We compute $B_1$ by using Lemma 3, Lemma (3-3) in ^[54],

The quantity $| B_{2} |$ is subsequently approximated as follows:

By using the Cauchy inequality, we can get

Hence,

We have now determined an estimate for the term $|B_{3}|$. The Legendre-Gauss integration formula (3.3) has been applied to achieve this result.

We obtain it by using the Cauchy-Schwarz inequality

By using the Lipschitz condition, we can write

which $L$ is Lipschitz condition, and $Max|\sigma(\varrho,\lambda)|\leq M$ and $L<1/M$.

Furthermore, by utilizing the triangle inequality, we derive that

Moreover, we can infer from Lemma 4 that

Consequently, the combination of (4.11), (4.14), and (4.19) results in the desired conclusion of this theorem.

5.   Numerical results

We examine several examples to validate the proposed methodology's effectiveness and accuracy. The absolute error (AE) is defined as the difference between the exact and measured values of an approximate solution. The programs used in this work are run on a PC with an Intel(R) Core(TM) i7-10510U CPU running 1.80 GHz and 2.30 GHz, with 2.00 GB of RAM, and Mathematica version 12 running the code.

The definition of the absolute error (AEs) is as follows:

Where $\mathcal{Z}_{Approx}(\varrho)$ and $\mathcal{Z}(\varrho)$ are the approximate and exact solutions at $z$, respectively.

Example 1. First, we present non-FFIDEs.

While $F(\varrho)$ is derived from the exact solution $\mathcal{Z}(\varrho) = \varrho^{3}-\varrho^{5}$ with $\alpha_{1} = 0.9$.

The methodology outlined in Section 3.2 is applied to address the problem. The outcomes presented in Table 1 depict the AE corresponding to various choices of $\nu_{1}$. The AE for Example 1 is illustrated in Figure 1 for $\nu_{1} = 6$, while the AE is illustrated in Figure 2 for $\nu_{1} = 12$. Additionally, Figure 3 displays the numerical approximation of $\kappa_{\nu_{1}}(\varrho)$ for various fractional orders $(\alpha_{1})$. From Table 1 and Figures 1 and 2, it is evident that the proposed algorithm yields a highly accurate approximation of the exact solution, even with a minimal number of collocation points. The results highlight the superior accuracy of our method. Taking $\alpha_{1} = 0.5$ and $\nu_{1} = 12$, we obtain $\mathcal{Z}_{\nu_{1}}(\varrho)$ of Example 1 as

Taking $\alpha_{1} = 0.7$ and $\nu_{1} = 12$, we obtain $\mathcal{Z}_{\nu_{1}}(\varrho)$ of Example 1 as

Taking $\alpha_{1} = 0.9$ and $\nu_{1} = 12$, we obtain $\mathcal{Z}_{\nu_{1}}(\varrho)$ of Example 1 as

Example 2. The following non-FFIDEs are considered ^[21]:

where $\alpha_{1} = 0.5$ and the exact solution $\mathcal{Z}(\varrho) = \varrho^{2}-\varrho$. We examine the convergence and computational processing time of our approach. It is evident that our proposed method outperforms Chebyshev wavelet ^[21]. Table 2 presents the root mean square errors between Chebyshev wavelet ^[21] and the current method across different $\nu_{1}$ values, along with the corresponding CPU time in seconds. Furthermore, in Figure 4, we represent the logarithmic graphs of $M_E$ (i.e., $\log_{10} M_E$) obtained by the proposed method with different values of $\nu_{1}$. It is demonstrated from the results of this example that the present scheme provides very highly accurate approximation of the solution for the problems and yields accurate convergence rates.

Example 3. Consider the following non-FFIDEs:

where $F$ is from exact solution $\mathcal{Z}(\varrho) = \varrho^{3}-\varrho^{2}$.

The outcomes of the L-G-LC method for various $\nu_{1}$ values are presented in Table 3. In Figure 5, we graph the AE curve of Example 3 for $\nu_{1} = 5$ and $\nu_{1} = 8$. The concordance between the curves of $\mathcal{Z}(\varrho)$ and $\mathcal{Z}_{Approx}(\varrho)$ from Example 3 is illustrated in Figure 6. From the results, we verify that our scheme reveals superior accuracy, even for just a few points.

Example 4. Consider the following non-FFIDEs

where $F$ is from exact solution $\mathcal{Z}(\varrho) = e^{-0.5\varrho}sin(\pi \varrho)$.

The outcomes of the L-G-LC method for various $\nu_{1}$ values are displayed in Table 4. Furthermore, in Figure 7, we depict the logarithmic graphs of $M_E$ (i.e., $\log_{10} M_E$) obtained by the proposed method for various values of $\nu_{1}$ and $\alpha_{1} = 0.2, 0.5, 0.9$. Thus, we have illustrated that the present method offers a precise approximation for problems characterized by nonlinearity and lack of smoothness.

Example 5. Finally, we present non-FFIDEs.

where the exact solution $\mathcal{Z}(\varrho) = \varrho$ for $\alpha_{1} = 1$.

Table 5 exhibits the root mean square errors with $\alpha_{1} = 1$. The exact and approximate solutions are graphed in Figure 8 for $\alpha_{1} = 1$. Figure 9–12 show the AE curve versus $\varrho$ in Example 5 for $\nu_{1} = 12$ with various $\alpha_{1} = 0.25, 0.5, 0.75, 1$, respectively. Moreover, in the absence of exact solutions, we have plotted the approximate solutions for various values of $\alpha_{1}$. The results show that our technique achieves greater accuracy, even for a few points.

Example 6. The following non-FFIDE with non-smooth solution is considered ^[53]:

where the exact solution $\mathcal{Z}(\varrho) = \sqrt{\varrho}$.

When the solution to a problem is not sufficiently smooth, the performance of a numerical scheme, particularly its order of convergence, can degrade significantly. This means that the accuracy of the approximation does not improve as quickly as expected when the grid or step size is refined. However, this challenge can be addressed effectively by incorporating fractional-order Legendre functions $\mathcal{L}_{\varepsilon}(\lambda^{\gamma})$ into the numerical approach. These functions extend the traditional Legendre polynomials to fractional orders, providing greater flexibility and adaptability. By tailoring the basis functions to better match the problem's irregularities, fractional-order Legendre functions can enhance the approximation's accuracy, maintaining or even improving the order of convergence despite the lack of smoothness in the solution. Figure 13 presents the AE of Example 6 for $\gamma=\frac{1}{2}$ and $\nu_{1}$. Figure 14 compares the exact and approximate solutions. We applied our technique by using fractional-order in this example with a non-smooth solution and we can see from the outcomes that the suggested scheme yields superior accuracy. Additionally, it should be noted that good approximations can be made with only a few points.

Example 7. Consider the following non-FFIDE with non-smooth solution ^[53]:

where $F(\varrho)$ is obtained from the exact solution $\mathcal{Z}(\varrho) = \varrho^{\frac{1}{2}}-\frac{1}{3!}\varrho^{\frac{3}{2}}+\frac{1}{5!}\varrho^{\frac{5}{2}}$.

For solving the nonsmooth solution, we used fractional-order Legendre $\mathcal{L}_{\varepsilon}(\lambda^{\gamma})$. For our algorithm, we get the maximum errors shown in Table 6, for various values of $\nu_{1}$ and $\gamma$. In Table 6 we compare our method with the approach presented in ^[55], in terms of the MAE. From the results, we verify that our scheme reveals superior accuracy, even for just a few points.

6.   Conclusions

In this study, we introduce a precise and efficient numerical algorithm based on the L-G-LC method to solve non-FFIDEs with initial value problems. The resolution of the nonlinear algebraic equations system was employed to solve the problem. Utilizing the LGL points as collocation nodes in the approximate solution preserves spectral convergence for the spatial variable. To illustrate the efficacy of the derived numerical algorithm, numerical examples were presented. The algorithm demonstrates efficiency, applicability to various operators, and extensibility to multi-dimensional problems, laying the foundation for future research. In subsequent investigations, we aim to address fractional integro-differential equations involving generalized formulations.

In conclusion, based on the theoretical formulation and numerical illustrations, we confirm that:

(1) The presented method yields accurate and reliable solutions, when compared with other approaches.

(2) The error decays exponentially as $\nu_{1}\rightarrow \infty$ in the case of smooth solutions.

(3) An upper bound for the absolute error in the approximate solution can be determined.

(4) If the solution is not smooth, then the order of convergence of the numerical scheme may deteriorate. However, this can be prevented by using fractional order Legendre functions.

Furthermore, we will study in the future a variable-order and stochastic fractional integro-differential equation with a non-smooth solution by utilizing a combination of the finite difference and spectral methods.

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 extend their appreciation to the Deanship of Scientific Research at King Khalid University for funding this work through Large Research Group project under grant number (R.G.P. 2/81/44).

Conflict of interest

The authors declare no conflicts of interest.