Utilizing Cubic B-Spline Collocation Technique for Solving Linear and Nonlinear Fractional Integro-Differential Equations of Volterra and Fredholm Types

Abstract

Fractional integro-differential equations (FIDEs) of both Volterra and Fredholm types present considerable challenges in numerical analysis and scientific computing due to their complex structures. This paper introduces a novel approach to address such equations by employing a Cubic B-spline collocation method. This method offers a robust and systematic framework for approximating solutions to the FIDEs, facilitating precise representations of complex phenomena. Within this research, we establish the mathematical foundations of the proposed scheme, elucidate its advantages over existing methods, and demonstrate its practical utility through numerical examples. We adopt the Caputo definition for fractional derivatives and conduct a stability analysis to validate the accuracy of the method. The findings showcase the precision and efficiency of the scheme in solving FIDEs, highlighting its potential as a valuable tool for addressing a wide array of practical problems.

1. Introduction

Efficient numerical methods play a vital role in approximating solutions to fractional integro-differential equations (FIDEs), which find applications in ac, physics, engineering, and biology. Among these methods, the Cubic B-spline collocation technique has emerged as a leading approach that is renowned for its accuracy and computational efficiency. Cubic B-splines, with their inherent flexibility and smoothness, excel at handling FIDEs by capturing the intricate behavior of functions involving fractional derivatives. Accurate numerical methods are indispensible for avoiding erroneous predictions in complex systems. The adaptability of the Cubic B-spline collocation technique has facilitated its successful application in various linear and nonlinear FIDEs, including those encountered in viscoelastic materials and population dynamics, underscoring its practical utility and versatility.

In recent years, significant efforts have been dedicated towards developing numerical methods capable of accurately approximating solutions to FIDEs. The Cubic B-spline collocation technique has emerged as a prominent approach among these methods. Cubic B-splines, known for their flexibility and smoothness properties, offer an excellent framework for approximating functions involving fractional derivatives. Benzahi et al. [1] delved into the numerical investigation of Fredholm FIDEs using the least squares method combined with a compact combination of shifted Chebyshev polynomials. This method, involving a series expansion in terms of Chebyshev polynomials, provides high accuracy, particularly in bounded domains. Yi and Huang [2] introduced a novel approach to tackle FIDEs with weakly singular kernels using the CAS (Chebyshev Adaptive Subinterval). This approach leverages the properties of wavelet functions to discretize the fractional derivative and integral operators efficiently while mantaining accuracy and stability. This method offers a promising alternative to existing numerical techniques, potentially enhancing the computational efficiency and accuracy for solving such equations. Mirzaee and Alipour [3] utilized Cubic B-spline approximation to address linear stochastic FIDEs examining the computational efficiency and accuracy of this method for stochastic modeling and analysis. Erfanian and Zeidabadi [4] proposed a Cubic B-spline finite element method in the complex plane to approximate solutions to linear Volterra IDEs, extending the applicability of the B-spline methods to problems involving complex-valued functions. Arshed [5] investigated the solution of FIDEs with weakly singular kernels using B-spline methods, emphasizing the effectiveness of B-spline approaches in addressing challenging problems. Qiao et al. [6] introduced an alternating direction implicit orthogonal spline collocation method for solving two-dimensional multi-term FIDEs, presenting a novel numerical technique for efficiently solving complex multi-dimensional fractional equations. Mohammed [7] explored a numerical solution approach for FIDEs using the least squares method and shifted Chebyshev polynomials, investigating the efficacy of least squares and Chebyshev polynomial methods in approximating solutions to FIDEs. Saeedi et al. [8] proposed a novel computational approach for solving nonlinear Fredholm FIDEs using a wavelet basis aiming to efficiently approximate the solutions of FIDEs. Awawdeh et al. [9] presented an analytic solution technique for FIDEs, focusing on equations with fractional derivatives and integrals employing Laplace transform and power series expansion methods to derive closed-form solutions for a class of FIDEs. Bhrawy and Alghamdi [10] introduced a shifted Jacobi–Gauss–Lobatto collocation method for solving nonlinear fractional Langevin equations involving two fractional orders in different intervals, presenting a novel collocation scheme based on shifted Jacobi–Gauss–Lobatto quadrature to approximate the solution of the fractional Langevin equation. Yang et al. [11] presented a spectral-collocation method for solving Fredholm FIDEs combining spectral methods with collocation techniques to approximate the solutions of Fredholm FIDEs. Bhrawy and Alofi [12] developed an operational matrix of fractional integration for shifted Chebyshev polynomials, offering a systematic approach to construct the operational matrix of fractional integration for shifted Chebyshev polynomials. Doha et al. [13] proposed efficient Chebyshev spectral methods for solving multi-term fractional order differential equations, introducing a spectral method based on Chebyshev polynomials. Irandoust et al. [14] applied Legendre wavelets to solving fractional differential equations, utilizing Legendre wavelets as basis functions to approximate the solution of fractional differential equations. Hashim et al. [15] introduced the homotopy analysis method (HAM) for solving fractional initial value problems (IVPs), which is a powerful analytical technique that constructs a homotopy between an IVP and a solvable auxiliary linear problem. Lio [16] presented the HAM as a general analytical technique for solving nonlinear problems, including FIDEs. Yang [17] proposed a numerical method for solving nonlinear Fredholm FIDEs, utilizing a hybrid approach combining block-pulse functions and Chebyshev polynomials to approximate the solution of these equations. Caputo and Fabrizio [18] explored the applications of new time and spatial fractional derivatives with exponential kernels, introducing fractional derivatives with exponential kernels that generalize the classical Caputo and Riemann–Liouville fractional derivatives. Yuan [19] investigated the use of Chebyshev wavelets to solving nonlinear fractional differential equations, employing Chebyshev wavelets as basis functions to approximate the solution of a specific nonlinear fractional differential equation. Li and Sun [20] introduced the Generalized Block-Pulse Operational Matrix (GBPOM) for numerically solving fractional differential equations, efficiently approximating the solutions of FDEs, with wide-ranging applications in numerical analysis. Araci [21] explored innovative identities concerning q-Genocchi numbers and polynomials with wide-ranging applications in numerical methodologies, particularly in addressing mathematical challenges like FIDEs. Gurbuz and Sezer [22] studied Laguerre polynomial solutions of initial and boundary value problems arising in science and engineering fields. Nazari and Shahmorad [23] applied the fractional differential transform method (FDTM) to FIDEs with non-local boundary conditions, obtaining analytical solutions for FIDEs. Ali et al. [24] presented a numerical solution for generalized nonlinear FIDEs with linear functional arguments using Chebyshev series, efficiently approximating generalized nonlinear FIDEs. Oa et al. [25] introduced numerical studies for solving FIDEs using the least squares method and Bernstein polynomials, providing accurate numerical solutions for FIDEs. Nigarchi and Nouri [26] proposed a numerical solution method for Volterra and Fredholm integral equations based on a special form of the Müntz Legendre polynomials. Beni [27] introduced the Legendre wavelet method combined with the Gauss quadrature rule for the numerical solution of FIDEs, obtaining reliable numerical solutions for FIDEs. Kumar et al. [28] conducted a comparative study of three numerical schemes for FIDEs, offering insights into the strengths and limitations of different numerical methods, including finite difference, finite element, and spectral methods. The book by Sun and Gao [29] delved into the theoretical foundations of fractional differential equations while emphasizing practical methodologies, particularly the finite difference method, for solving such equations.

Despite the significant strides made in obtaining numerical solutions for FIDEs, a literature gap persists in the development of versatile and precise numerical methods capable of handling both linear and nonlinear FIDEs featuring diverse integral operators. To bridge this gap, this paper suggests employing the Cubic B-spline collocation technique as a robust and adaptable approach to solve a broad spectrum of linear and nonlinear FIDEs incorporating Volterra and Fredholm integral operators. By harnessing the flexibility and computational efficiency of Cubic B-splines in conjunction with the collocation technique, this method aims to furnish accurate numerical solutions. This proposed approach holds promise in augmenting the efficiency and reliability of numerical investigations in analyzing intricate systems governed by FIDEs.

We examine the following linear and nonlinear FFIDEs:

$$D^{\alpha }\varphi \left(x\right)=f\left(x\right)+{\int }_a^{b}K(x,t)\varphi \left(t\right)dt,$$

(1)

$$D^{\alpha }\varphi \left(x\right)=f\left(x\right)+{\int }_a^{b}K(x,t){\left(\varphi \left(t\right)\right)}^{n}dt,\mathrm{where}n=2,3,4,\dots .$$

(2)

respectively. Additionally, we investigate the following VFIDE:

$$D^{\alpha }\varphi \left(x\right)=f\left(x\right)+{\int }_a^{x}K(x,t)\varphi \left(t\right)dt.$$

(3)

Equations (1)–(3) are subject to the following supplementary conditions:

• ${\varphi }^{\left(i\right)}={\delta }_i$.
• $n-1<\alpha \le n$, where n is a natural number ($n\in \mathbb{N}$).

In these equations, $f\left(x\right)$ and $K(x,t)$ represent given functions, and x and t are real variables ranging from 0 to 1. The term $D^{\alpha }\varphi \left(x\right)$ denotes the ${\alpha }^{th}$th Caputo fractional derivative of $\varphi \left(x\right)$, with $\varphi \left(x\right)$ being the function to be determined. If $\varphi \in C_{-1}^m$, where $1,m\in \mathbb{N}\cup 0$, then the Caputo fractional derivative of $\varphi \left(x\right)$ is given as follows:

$$\begin{array}{cc}\hfill D^{\alpha }\varphi \left(x\right)& =\left\{\begin{array}{cc}{J}^{\alpha }{\varphi }^m\left(x\right),\hfill & \mathrm{if}m-1<\alpha \le m,m\in \mathbb{N}\hfill \\ \frac{D^m\varphi \left(x\right)}{D{x}^m},\hfill & \mathrm{if}\alpha =m.\hfill \end{array}\right.\hfill \end{array}$$

(4)

This manuscript follows this structure: Section 2 provides a fundamental review of fractional derivatives for the readers’ reference. Section 3 details the derivation of numerical schemes for FFIDEs and VFIDEs. Moreover, the main steps of the scheme are described in algorithm in this section. Stability analysis is the focus of Section 4, while Section 5 delves into the convergence analysis of the proposed schemes. Section 6 presents the numerical findings, and concluding remarks are presented in Section 7.

2. Basic Definitions of Fractional Derivatives

In this section, we present the foundational definitions and properties of fractional derivatives.

Definition 1

([29]). A real function $\varphi \left(x\right)$, where $x>0$, belongs to the space $C_{\mu }$, where $\mu \in \mathbb{R}$, if there exists a real number $p>\mu$ such that $\varphi \left(x\right)=x^p{\varphi }_1\left(x\right)$, with ${\varphi }_1\left(x\right)\in C[0,1)$.

Definition 2

([29]). A function $\varphi \left(x\right)$, where $x>0$, is considered to be in the space $C_{\mu }^m$, where $m\in \mathbb{N}$, if ${\varphi }^{\left(m\right)}\in C_{\mu }$.

Definition 3

([29]). The left-sided Riemann–Liouville fractional integral operator of order $\alpha >0$ for a function $\varphi \in C_{\mu }$, where $\mu \ge -1$, is defined as

$$J^{\alpha }\varphi \left(x\right)=\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_0^x\frac{\varphi \left(t\right)}{{(x-t)}^{1-\alpha }}dt,\alpha >0,x>0.$$

For $\alpha =0$, we have $J^0\varphi \left(x\right)=\varphi \left(x\right).$ The following properties hold:

• $J^{\alpha }{J}^{\gamma }\varphi =J^{\alpha +\gamma }\varphi$.
• $J^{\alpha }{x}^{\gamma }=\frac{\mathsf{\Gamma }(\gamma +1)}{\mathsf{\Gamma }(\gamma -\alpha +1)}{x}^{\gamma -\alpha }$, for $\gamma >-1,\alpha >0,x>0$.
• $J^{\alpha }{D}^{\alpha }\varphi \left(x\right)=\varphi \left(x\right)-{\sum }_{k=0}^{m-1}{\varphi }^{\left(k\right)}\left(0\right)\frac{x^k}{k!}$, for $m-1<\alpha \le m,x>0$.
• $D^{\alpha }{J}^{\alpha }\varphi \left(x\right)=\varphi \left(x\right)$, for $m-1<\alpha \le m,x>0$.
• $D^{\alpha }C=0$.
•  

$$\begin{array}{cc}\hfill J^{\alpha }{x}^{\beta }& =\left\{\begin{array}{cc}0,\hfill & \mathrm{if}\beta \in {\mathbb{N}}_0,\beta <\lfloor \alpha \rfloor \hfill \\ \frac{\mathsf{\Gamma }(\beta +\alpha )}{\mathsf{\Gamma }(\beta -\alpha +1)}{x}^{\beta -\alpha },\hfill & \mathrm{if}\beta \ge \lfloor \alpha \rfloor .\hfill \end{array}\right.\hfill \end{array}$$

(5)

Here, ${\mathbb{N}}_0=0,1,2,\dots$, and $\lfloor \alpha \rfloor$ denotes the smallest integer greater than or equal to $\alpha$.

3. Numerical Scheme

In this section, we develop numerical techniques utilizing Cubic B-splines for solving FFIDEs and VFIDEs as expressed in Equations (1) and (3). The numerical method for the nonlinear FFIDE (2) is same as lthe inear FFIDE.

Let $h=\frac{b-a}{n}$ be the step size, where n is a positive integer. Setting $x_j=jh(0\le j\le n)$ partitions $a\le x\le b$ into n subintervals $[x_j,x_{j+1}]$ of equal length h, $j=0,1,2,\dots ,n-1$, where $a=x_0<x_1<\dots <x_{n-1}<x_n=b$. To obtain the approximate solution $\mathsf{\Phi }\left(x\right)$ to the exact solution $\varphi \left(x\right)$ of (1)–(3), we require

$$\mathsf{\Phi }\left(x\right)=\sum _{j=-1}^{n+1}{C}_j{B}_j\left(x\right),x\in [x_j,x_n].$$

(6)

Here, $C_j$s are the unknowns to be evaluated and $B_j\left(x\right)$ are Cubic B-spline basis functions.

$$\begin{array}{cc}\hfill B_j^3\left(x\right)& =\frac{1}{6h^3}\left\{\begin{array}{cc}{(x-x_j)}^3,\hfill & x\in [x_j,x_{j+1})\hfill \\ h^3+3h^2(s-s_{j+1})+3h{(x-x_{j+1})}^2-3{(x-x_{j+1})}^3,\hfill & x\in [x_{j+1},x_{j+2})\hfill \\ h^3+3h^2(x_{j+3}-x)+3h{(x_{j+3}-x)}^2-3{(x_{j+3}-x)}^3,\hfill & s\in [x_{j+2},x_{j+3})\hfill \\ {(x_{j+4}-x)}^3,\hfill & x\in [x_{j+3},x_{j+4})\hfill \\ 0,\hfill & otherwise.\hfill \end{array}\right.\hfill \end{array}$$

Only $B_{j-1}\left(x\right),B_j\left(x\right)$, and $B_{j+1}\left(x\right)$ survive at the grid point $x_j$ because of local support property of CuBS. Consequently, the approximation $\mathsf{\Phi }\left(x\right)$ at $x_j$ is given as

$$\mathsf{\Phi }\left(x_j\right)=\sum _{j=-1}^{n+1}{C}_j{B}_j\left(x_j\right),j=0,1,2,\dots ,n$$

(7)

The values of $\mathsf{\Phi }\left(x\right)$, ${\mathsf{\Phi }}_x\left(x\right)$, and ${\mathsf{\Phi }}_{xx}\left(x\right)$ at the grid point $x_j$ are given by

$$\left\{\begin{array}{c}{\displaystyle \mathsf{\Phi }\left(x\right)={\omega }_1{C}_{j-1}+{\omega }_2{C}_j+{\omega }_1{C}_{j+1},}\hfill \\ {\displaystyle {\mathsf{\Phi }}_x\left(x\right)=-{\nu }_1{C}_{j-1}+{\nu }_1{C}_{j+1},}\hfill \\ {\displaystyle {\mathsf{\Phi }}_{xx}\left(x\right)={\sigma }_1{C}_{j-1}+{\sigma }_2{C}_j+{\sigma }_1{C}_{j+1}.}\hfill \end{array}\right.$$

where ${\displaystyle {\omega }_1=\frac{1}{6}}$, ${\displaystyle {\omega }_2=\frac{4}{6},}$ ${\displaystyle {\nu }_1=\frac{1}{2h},}$${\displaystyle {\sigma }_1=\frac{1}{h^2},}$ and ${\displaystyle {\sigma }_2=-\frac{2}{h^2}}$.

Now, the Caputo fractional derivative of order $\alpha$, $0<\alpha <1$ is given by

$$D_x^{\alpha }\varphi \left(x\right)=\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^x\frac{{\varphi }^{\prime }\left(t\right)}{{(x-t)}^{\alpha }}dt.$$

For $n-1<\alpha <n$, we have

$$D_x^{\alpha }\varphi \left(x\right)=\frac{1}{\mathsf{\Gamma }(n-\alpha )}{\int }_0^x\frac{{\varphi }^{\left(n\right)}\left(t\right)}{{(x-t)}^{\alpha -(n-1)}}dt.$$

At $x=x_n$, we have

$$\begin{array}{cc}\hfill D_x^{\alpha }{\varphi \left(x\right)|}_{x=x_n}& =\frac{1}{\mathsf{\Gamma }(n-\alpha )}{\int }_0^{x_n}\frac{{\varphi }^{\left(n\right)}\left(t\right)}{{(x_n-t)}^{\alpha -(n-1)}}dt\hfill \\ \hfill & =\frac{1}{\mathsf{\Gamma }(n-\alpha )}\sum _{k=1}^n{\int }_{x_{k-1}}^{x_k}\frac{{\varphi }^{\left(n\right)}\left(t\right)}{{(x_n-t)}^{\alpha -(n-1)}}dt.\hfill \end{array}$$

(8)

Now, using Equation (8) in Equations (1) and (3), we obtain

$$\frac{1}{\mathsf{\Gamma }(n-\alpha )}\sum _{k=1}^n{\int }_{x_{k-1}}^{x_k}\frac{{\varphi }^{\left(n\right)}\left(t\right)}{{(x_n-t)}^{\alpha -(n-1)}}dt=f\left(x\right)+{\int }_a^{b}K(x,t)\varphi \left(t\right)dt$$

(9)

and

$$\frac{1}{\mathsf{\Gamma }(n-\alpha )}\sum _{k=1}^n{\int }_{x_{k-1}}^{x_k}\frac{{\varphi }^{\left(n\right)}\left(t\right)}{{(x_n-t)}^{\alpha -(n-1)}}dt=f\left(x\right)+{\int }_a^{x}K(x,t)\varphi \left(t\right)dt.$$

(10)

Since $\mathsf{\Phi }\left(x\right)$ is an approximation of $\varphi \left(x\right)$, substituting Equation (6) in Equations (9) and (10), we obtain

$$\frac{1}{\mathsf{\Gamma }(n-\alpha )}\sum _{k=1}^n{\int }_{x_{k-1}}^{x_k}\frac{{\displaystyle \sum _{j=-1}^{n+1}}{C}_j{B}_j^{\left(n\right)}\left(t\right)}{{(x_n-t)}^{\alpha -(n-1)}}dt=f\left(x\right)+{\int }_a^{b}K(x,t)\sum _{j=-1}^{n+1}{C}_j{B}_j\left(t\right)dt$$

(11)

and

$$\frac{1}{\mathsf{\Gamma }(n-\alpha )}\sum _{k=1}^n{\int }_{x_{k-1}}^{x_k}\frac{{\displaystyle \sum _{j=-1}^{n+1}}{C}_j{B}_j^{\left(n\right)}\left(t\right)}{{(x_n-t)}^{\alpha -(n-1)}}dt=f\left(x\right)+{\int }_a^{x}K(x,t)\sum _{j=-1}^{n+1}{C}_j{B}_j\left(t\right)dt.$$

(12)

Substituting $x=x_j$ in Equations (11) and (12), we obtain

$$\frac{1}{\mathsf{\Gamma }(n-\alpha )}\sum _{k=1}^n{\int }_{x_{k-1}}^{x_k}\frac{{\displaystyle \sum _{j=-1}^{n+1}}{C}_j{B}_j^{\left(n\right)}\left(t\right)}{{(x_n-t)}^{\alpha -(n-1)}}dt=f\left(x_j\right)+{\int }_a^{b}K(x_j,t)\sum _{j=-1}^{n+1}{C}_j{B}_j\left(t\right)dt$$

(13)

and

$$\frac{1}{\mathsf{\Gamma }(n-\alpha )}\sum _{k=1}^n{\int }_{x_{k-1}}^{x_k}\frac{{\displaystyle \sum _{j=-1}^{n+1}}{C}_j{B}_j^{\left(n\right)}\left(t\right)}{{(x_n-t)}^{\alpha -(n-1)}}dt=f\left(x_j\right)+{\int }_a^{x_j}K(x_j,t)\sum _{j=-1}^{n+1}{C}_j{B}_j\left(t\right)dt.$$

(14)

$j=0,1,2,3,\dots ,n$. Letting

$$G\left(x_j\right)=\frac{1}{\mathsf{\Gamma }(n-\alpha )}\sum _{k=1}^j{\int }_{x_{k-1}}^{x_k}\frac{{\displaystyle \sum _{j=-1}^{n+1}}{C}_j{B}_j^{\left(n\right)}\left(t\right)}{{(x_n-t)}^{\alpha -(n-1)}}dt,j=1,2,\dots ,n.$$

in Equations (13) and (14), we obtain

$$G\left(x_j\right)=f\left(x_j\right)+{\int }_a^{b}K(x_j,t)\sum _{j=-1}^{n+1}{C}_j{B}_j\left(t\right)dt,j=0,1,2,\dots ,n$$

(15)

and

$$G\left(x_j\right)=f\left(x_j\right)+{\int }_a^{x_j)}K(x_j,t)\sum _{j=-1}^{n+1}{C}_j{B}_j\left(t\right)dt,j=0,1,2,\dots ,n.$$

(16)

If we put $x=x_j$, we obtain $n+1$ equations in $n+3$ unknowns. To obtain two more equations, we introduce conditions given by

$${\mathsf{\Phi }}^{″}\left(x_0\right)=0,{\mathsf{\Phi }}^{″}\left(x_n\right)=0,$$

(17)

which are used to eliminate $C_{-1}$, $C_{n+1}$. Consequently, we obtain $n+1$ equations in $n+1$ unknowns, which can be written in matrix form as

$$AC=F,$$

(18)

where

$$C=\left[\begin{array}{c}{c}_0\\ c_1\\ ⋮\\ ⋮\\ c_{n-1}\\ c_n\end{array}\right],F=\left[\begin{array}{c}f\left(x_0\right)\\ f\left(x_1\right)\\ ⋮\\ ⋮\\ f\left(x_{n-1}\right)\\ f\left(x_n\right)\end{array}\right],$$

and

$$A=\left[\begin{array}{cccccc}\omega & 0& 0& 0& \dots & 0\\ {\lambda }_1& {\lambda }_2& {\lambda }_1& 0& \dots & 0\\ 0& {\lambda }_1& {\lambda }_2& {\lambda }_1& \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & 0\\ 0& \dots & 0& {\lambda }_1& {\lambda }_2& {\lambda }_1\\ 0& \dots & 0& 0& 0& \omega \end{array}\right],$$

(19)

where

$$\begin{array}{cc}\hfill {\lambda }_1& \hfill =a_{ij},fori\ne j,\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\lambda }_2& \hfill =a_{ij},fori=j,\hfill \end{array}$$

$$\begin{array}{cc}\hfill \omega & \hfill =-\frac{{\lambda }_1{\sigma }_2}{{\sigma }_1}+{\lambda }_2,\hfill \end{array}$$

and, for FFIDEs,

$$\begin{array}{cc}\hfill a_{ij}& \hfill =\frac{1}{\mathsf{\Gamma }(n-\alpha )}\sum _{k=1}^j{\int }_{x_{k-1}}^{x_k}\frac{{\displaystyle \sum _{j=-1}^{n+1}}{B}_j^{\left(n\right)}\left(t\right)}{{(x_n-t)}^{\alpha -(n-1)}}dt-{\int }_a^{b}K(x_j,t)\sum _{j=-1}^{n+1}{B}_j\left(t\right)dt\hfill \end{array}$$

for $j=1,...,n+1.$

For VFIDEs,

$$\begin{array}{cc}\hfill a_{ij}& \hfill =\frac{1}{\mathsf{\Gamma }(n-\alpha )}\sum _{k=1}^j{\int }_{x_{k-1}}^{x_k}\frac{{\displaystyle \sum _{j=-1}^{n+1}}{B}_j^{\left(n\right)}\left(t\right)}{{(x_n-t)}^{\alpha -(n-1)}}dt-{\int }_a^{x_j}K(x_j,t)\sum _{j=-1}^{n+1}{B}_j\left(t\right)dt\hfill \end{array}$$

for $j=1,\dots ,n+1.$ Once the matrices A and F are found, then $C_j$s can be evaluated easily. Algorithm 1 describing the presented scheme is given below:

Algorithm 1 Algorithm Describing the Proposed Scheme

Input $n;a;b;f\left(x\right),K(x,t),\alpha$: Output Solution vector C containing coefficients $C_j$: Procedure:    Initialize arrays C, F with size $n+1$.    Initialize matrix A with size $(n+1)\times (n+1)$.    Compute step size $h=\frac{b-a}{n}$.    Generate partition points $x_j=j\ast hforj=0,1,\dots ,n$.    Compute the elements of matrix A.    Compute the elements of vector F    Apply boundary conditions to modify A and F.    Solve the linear system $AC=F$ for the coefficients C using suitable numerical technique.    Output the solution vector C containing the coefficients $C_j$.

4. Stability Analysis

To investigate the accuracy of the proposed scheme, it is essential to examine whether it is stable or not. For this purpose, we prove the following:

Theorem 1.

Scheme (18) is stable.

Proof. 

System (18) is given by

$$AC=F.$$

Let us assume that ${\xi }_N$ and ${\varpi }_N$ are the modest alterations in A and F, respectively. Furthermore, let $\beta$ be the solution to system (18), so that

$$(A+{\xi }_N)\beta =(F+{\varpi }_N).$$

Since A is non-singular, $A+{\xi }_N$ is also non-singular. Now, let

$$\parallel {\xi }_N\parallel <\frac{1}{2\parallel A^{-1}\parallel }$$

so that

$$\parallel {(A+{\xi }_N)}^{-1}\parallel \le 2\parallel A^{-1}\parallel .$$

Equation (18) implies

$$\begin{array}{cc}\hfill C& =A^{-1}F,\hfill \\ \hfill \Rightarrow C-\beta & =A^{-1}F-\beta ,\hfill \\ \hfill \Rightarrow C-\beta & =A^{-1}F-\frac{F+{\varpi }_N}{A+{\xi }_N},\hfill \\ \hfill \Rightarrow (C-\beta )(A+{\xi }_N)& =A^{-1}F(A+{\xi }_N)-(F+{\varpi }_N),\hfill \\ \hfill \Rightarrow (C-\beta )(A+{\xi }_N)& =C(A+{\xi }_N)-(AC+{\varpi }_N),\hfill \\ \hfill \Rightarrow C-\beta & ={(A+{\xi }_N)}^{-1}(C{\xi }_N-{\varpi }_N),\hfill \end{array}$$

(20)

where we use the fact that $AC=CA$. As A is strictly diagonally dominant,

$$\parallel A^{-1}{\parallel }_{\infty }\le [\underset{0\le i\le N}{min}\left(\right|a_{ii}|-\sum _{i\ne j}|a_{ij}{\left|\right)]}^{-1}<\infty .$$

(21)

Letting

$$[\underset{0\le i\le N}{min}\left(\right|a_{ii}|-\sum _{i\ne j}|a_{ij}{\left|\right)]}^{-1}=q<\infty$$

(22)

so that, from (20), we have

$$\begin{array}{cc}\hfill {\parallel C-\beta \parallel }_{\infty }& =\parallel {(A+{\mathsf{\Omega }}_N)}^{-1}{\parallel }_{\infty }{\parallel (C{\xi }_N-{\varpi }_N)\parallel }_{\infty },\hfill \\ \hfill & \le 2\parallel A^{-1}{\parallel }_{\infty }{\parallel C\parallel }_{\infty }\parallel {\xi }_N{\parallel }_{\infty }+\parallel {\varpi }_N{)\parallel }_{\infty },\hfill \\ \hfill & \le {2q\parallel C\parallel }_{\infty }\parallel {\xi }_N{\parallel }_{\infty }+\parallel {\varpi }_N{)\parallel }_{\infty }.\hfill \end{array}$$

(23)

inequality (23) shows that the presented scheme is stable. □

5. Convergence Analysis

The suggested scheme’s convergence analysis is provided in this section. The detailed description is as follows:

Theorem 2.

If $\hat{\varphi }\left(x\right)$ is the exact solution of Equation (1) and $\hat{b}\left(x\right)$ is the B-spline collocation approximation to $\hat{\varphi }\left(x\right)$, the technique is then second-order convergent, and

$$\parallel \hat{\varphi }\left(x\right)-\hat{b}{\left(x\right)\parallel }_{\infty }\le \sigma h^2,$$

(24)

where $\sigma ={\kappa }_0Łh^2+R$ is a finite constant.

Proof. 

Assume that $\hat{\varphi }\left(x\right)$ is the exact solution of Equation (1), where $\hat{\varphi }\left(x\right)$ is approximated by $\hat{b}\left(x\right)$ so that

$$\hat{b}\left(x\right)=\sum _{j=-1}^{N+1}\widehat{C_j}{B}_j^3\left(x\right),$$

(25)

where $\hat{C}=(\widehat{C_{-}1},\widehat{C_0},\dots ,\widehat{C_{N+1}}).$ Moreover, let $\tilde{b}\left(x\right)$ be the evaluated Cubic B-spline-based collocation approximation of $\hat{b}\left(x\right)$, namely

$$\tilde{b}\left(x\right)=\sum _{j=-1}^{N+1}\widetilde{C_j}{B}_j^3\left(x\right),$$

(26)

where $\tilde{C}=(\widetilde{C_{-1}}$, $\widetilde{C_0},\dots ,\widetilde{C_{N+1}}).$ To approximate the errors, $\parallel \hat{\varphi }\left(x\right)-\hat{b}{\left(x\right)\parallel }_{\infty }$, we have to determine the errors $\parallel \hat{\varphi }\left(x\right)-\tilde{b}{\left(x\right)\parallel }_{\infty }$ and $\parallel \tilde{b}\left(x\right)-\hat{b}{\left(x\right)\parallel }_{\infty }$ separately. To compute $\tilde{b}\left(x\right)$ and $\hat{b}\left(x\right)$, the values of vectors $\hat{C}$ and $\tilde{C}$ must be computed from two linear equations,

$$A\hat{C}=\hat{F},$$

(27)

and

$$A\tilde{C}=\tilde{F}.$$

(28)

Now, by subtracting (28) from (27), we obtain

$$A(\tilde{C}-\hat{C})=\tilde{F}-\hat{F}.$$

(29)

The specification of matrix A in Equation (19) makes A strictly diagonally dominant, making it non-singular. Thus,

$$\begin{array}{c}\hfill (\tilde{C}-\hat{C})=A^{-1}(\tilde{F}-\hat{F}).\end{array}$$

Taking the infinity norm of the above equation, we obtain

$$\parallel (\tilde{C}-\hat{C}){\parallel }_{\infty }=\parallel A^{-1}{\parallel }_{\infty }{\parallel (\tilde{F}-\hat{F})\parallel }_{\infty }.$$

(30)

Let the sum of the $ith$ row of matrix $A={\left[a_{ij}\right]}_{(N+1)\times (N+1)}$ be ${\tau }_i(0\le i\le N)$. Then, we have

$$\begin{array}{c}\hfill {\tau }_0=\sum _{j=0}^N{a}_{0j}=\omega ,\end{array}$$

$$\begin{array}{c}\hfill {\tau }_i=\sum _{j=0}^N{a}_{ij}={\lambda }_2+2{\lambda }_{1}i=1,\dots ,N-1,\end{array}$$

$$\begin{array}{c}\hfill {\tau }_N=\sum _{j=0}^N{a}_{Nj}=\omega .\end{array}$$

It is well known in the theory of matrices that

$$\begin{array}{c}\hfill \sum _{j=0}^N{a}_{ij}^{-1}=1,i=0,1,\dots ,N.\end{array}$$

Here, $a_{ij}^{-1}$ represent the entries of $A^{-1}$. Now,

$$\parallel A^{-1}{\parallel }_{\infty }=\sum _{j=0}^N|a_{ij}^{-1}|\le \frac{1}{\tau },$$

(31)

where $\tau ={min}_{0\le i\le N}{\tau }_i=min(\omega ,{\lambda }_2+2{\lambda }_1).$ We substitute (31) into (30) to acquire

$$\parallel (\tilde{C}-\hat{C}){\parallel }_{\infty }\le \frac{1}{\tau }{\parallel (\tilde{F}-\hat{F})\parallel }_{\infty }.$$

(32)

To compute the upper bound of $\parallel (\tilde{F}-\hat{F}){\parallel }_{\infty }$, we have, from (16),

$$\begin{array}{cc}\hfill |\widetilde{F_i}-\widehat{F_i}|=G\left(x_j\right)-{\int }_a^{x_j)}K(x_j,t)\sum _{j=-1}^{n+1}{C}_j{B}_j\left(t\right)dt& \hfill \le |\widetilde{{\varphi }_i}-\widehat{{\varphi }_i}|.\hfill \end{array}$$

(33)

To simplify the RHS of (33), we use Theorem 1 from [30] to obtain

$$|\widetilde{{\varphi }_i}-\widehat{{\varphi }_i}|=|{\tilde{b}}_i\left(x\right)-{\hat{b}}_i\left(x\right)|\le {\kappa }_0Łh^4$$

(34)

where Ł is a constant. Thus, from (34), we can write (33) as

$$|\widetilde{F_i}-\widehat{F_i}|\le {\kappa }_0Łh^4.$$

(35)

Let ${\kappa }_0Łh^4=R_q$; then, (35) becomes

$$|\widetilde{F_i}-\widehat{F_i}|\le R_q.$$

(36)

Using (36) in (32), we obtain

$$\begin{array}{c}\hfill \parallel (\tilde{C}-\hat{C}){\parallel }_{\infty }\le \frac{1}{\tau }{R}_q=R{h}^2,\end{array}$$

(37)

where $R{h}^2=\frac{1}{\tau }{R}_q=max(\omega ,{\delta }_2+2{\delta }_1).$ To proceed further, we need the following from [31]: □

Theorem 3.

The B-splines $\{B_{-1},B_0,B_1,\cdots ,B_N,B_{N+1}\}$ satisfy $|{\sum }_{k=-1}^{N+1}{B}_k\left(x\right)|\le 1$.

Subtracting (26) from (25), we have

$$\begin{array}{c}\hfill \tilde{b}\left(x\right)-\hat{b}\left(x\right)=\sum _{j=-1}^{N+1}(\widetilde{C_j}-\widehat{C_j})B_j^3\left(x\right).\end{array}$$

Taking the infinity norm on both sides, we obtain

$$\begin{array}{cc}\hfill \parallel \tilde{b}\left(x\right)-\hat{b}{\left(x\right)\parallel }_{\infty }& =\parallel \sum _{j=-1}^{N+1}(\widetilde{C_j}-\widehat{C_j})B_j^3\left(x\right){\parallel }_{\infty },\hfill \\ & \le |\sum _{j=-1}^{N+1}{B}_j^3\left(x\right)| \parallel (\widetilde{C_j}-\widehat{C_j})B_j^3\left(x\right){\parallel }_{\infty },\hfill \\ & \le R{h}^2,\hfill \end{array}$$

that is,

$$\parallel \tilde{b}\left(x\right)-\hat{b}{\left(x\right)\parallel }_{\infty }\le R{h}^2.$$

(38)

Taking the norm of Equation (34), it is inferred that

$$\parallel \hat{\varphi }\left(x\right)-\tilde{b}{\left(x\right)\parallel }_{\infty }\le {\kappa }_0Łh^4,$$

(39)

so that, from (38) and (39), we have

$$\begin{array}{cc}\hfill \parallel \hat{\varphi }\left(x\right)-\hat{b}{\left(x\right)\parallel }_{\infty }& \le \parallel \hat{\varphi }\left(x\right)-\tilde{b}{\left(x\right)\parallel }_{\infty }+{\parallel \tilde{b}\left(x\right)-\hat{b}\left(x\right)\parallel }_{\infty },\hfill \\ \hfill & \hfill \le {\kappa }_0Łh^4+R{h}^2=\sigma h^2,\hfill \end{array}$$

(40)

where $\sigma ={\kappa }_0Łh^2+R$.

6. Numerical Findings and Discussion

The effectiveness and accuracy of the suggested schemes are validated in this section through the utilization of the $L_{\infty }$ error norm. The numerical results yielded by the proposed schemes are compared, with Mathematica 9 employed for obtaining both the numerical and graphical results.

Example 1

([25]). Consider the FFIDE

$$D^{\frac{1}{2}}\varphi \left(x\right)=\frac{\left(\frac{8}{3}\right)x^{\frac{3}{2}}-2x^{0.5}}{\sqrt{\pi }}+\frac{x}{12}+{\int }_0^{1}xt\varphi \left(t\right)dt,0\le x,t\le 1,$$

with

$$\varphi \left(0\right)=0$$

It has the exact solution $\varphi \left(x\right)=x^2-x$. The proposed method was utilized in the above example to obtain an approximate solution. Figure 1 illustrates a comparison between the approximate and exact solutions, while the error profile is depicted in Figure 2.

Table 1. Comparison of absolute errors for Example 1 with $n=3$.

x	Our Results	Results in [25]
0.1	$4.23273\times {10}^{-14}$	$1.6795\times {10}^{-5}$
0.2	$4.93217\times {10}^{-14}$	$2.7482\times {10}^{-6}$
0.3	$5.43732\times {10}^{-14}$	$7.2646\times {10}^{-6}$
0.4	$5.73153\times {10}^{-14}$	$1.3243\times {10}^{-5}$
0.5	$5.76206\times {10}^{-14}$	$1.5190\times {10}^{-5}$
0.6	$5.43454\times {10}^{-14}$	$1.3106\times {10}^{-5}$
0.7	$4.63518\times {10}^{-14}$	$6.9911\times {10}^{-6}$
0.8	$3.43614\times {10}^{-14}$	$3.1530\times {10}^{-6}$
0.9	$2.10942\times {10}^{-14}$	$1.7325\times {10}^{-5}$

presents a comparison of absolute errors and approximate solutions between the proposed method and those from [25].

The approximate solution for the case when $n=3$ is provided as follows:

$$\begin{array}{c}\hfill \mathsf{\Phi }\left(x\right)=\left\{\begin{array}{cc}5.55112\times {10}^{-17}-x+x^2-1.33227\times {10}^{-15}{x}^3,\hfill & x\in [0,\frac{1}{3})\hfill \\ -4.85723\times {10}^{-17}-x+x^2+8.32667\times {10}^{-16}{x}^3,\hfill & x\in [\frac{1}{3},\frac{2}{3})\hfill \\ 1.11022\times {10}^{-15}-x+x^2-2.22045\times {10}^{-15}{x}^3,\hfill & x\in [\frac{2}{3},1).\hfill \end{array}\right.\end{array}$$

The approximate solution when $n=10$ is given as

$$\begin{array}{c}\hfill \mathsf{\Phi }\left(x\right)=\left\{\begin{array}{cc}4.66988\times {10}^{-15}-x+x^2-1.53406\times {10}^{-11}{x}^3,\hfill & x\in [0,\frac{1}{10})\hfill \\ -2.71727\times {10}^{-14}-x+x^2+1.64988\times {10}^{-11}{x}^3,\hfill & x\in [\frac{1}{10},\frac{1}{5})\hfill \\ -2.95042\times {10}^{-13}-x+x^2+5.00009\times {10}^{-11}{x}^3,\hfill & x\in [\frac{1}{5},\frac{3}{10})\hfill \\ ⋮\hfill \\ ⋮\hfill \\ 2.87919\times {10}^{-10}-x+x^2-5.65507\times {10}^{-10}{x}^3,\hfill & x\in [\frac{7}{10},\frac{4}{5})\hfill \\ -7.76268\times {10}^{-12}-x+x^2+1.19798\times {10}^{-11}{x}^3,\hfill & x\in [\frac{4}{5},\frac{9}{10})\hfill \\ 3.25429\times {10}^{-12}-x+x^2-3.11928\times {10}^{-12}{x}^3,\hfill & x\in [\frac{9}{10},1).\hfill \end{array}\right.\end{array}$$

Example 2

([25]). Consider the FFIDE,

$$D^{\frac{5}{6}}\varphi \left(x\right)=-\frac{3}{91}\frac{x^{\frac{1}{6}}\mathsf{\Gamma }\left(\frac{5}{6}\right)(-91+216x^2)}{\pi }+(5-2e)x+{\int }_0^{1}x{e}^t\varphi \left(t\right)dt$$

with $\varphi \left(0\right)=0$.

The exact solution is given by $\varphi \left(x\right)=x-x^3$. We applied the proposed scheme to the aforementioned problem and obtained approximate results. Figure 3 plots the exact and approximate solutions, with tremendous agreement between the two solutions. The absolute error is depicted in Figure 4. Additionally,

Table 2. Comparison of absolute errors when $n=3$ for Example 2.

x	Present Method	Method in [25]
0.1	$1.0255\times {10}^{-14}$	$6.3036\times {10}^{-5}$
0.2	$9.5201\times {10}^{-15}$	$2.5659\times {10}^{-5}$
0.3	$8.7152\times {10}^{-15}$	$6.8668\times {10}^{-6}$
0.4	$7.7715\times {10}^{-15}$	$3.2130\times {10}^{-5}$
0.5	$4.6629\times {10}^{-15}$	$4.7716\times {10}^{-5}$
0.6	$3.6082\times {10}^{-15}$	$5.1213\times {10}^{-5}$
0.7	$1.9817\times {10}^{-14}$	$4.0208\times {10}^{-5}$
0.8	$4.1799\times {10}^{-14}$	$1.2286\times {10}^{-5}$
0.9	$5.8481\times {10}^{-14}$	$3.4964\times {10}^{-5}$

presents a comparison of the absolute errors and approximate solutions between the proposed scheme and those from [25].

The approximate solution for Example 2 when $n=3$ is given as

$$\begin{array}{c}\hfill \mathsf{\Phi }\left(x\right)=\left\{\begin{array}{cc}1.08663\times {10}^{-14}+x-1.77636\times {10}^{-15}{x}^2-x^3,\hfill & x\in [0,\frac{1}{3})\hfill \\ 1.17475\times {10}^{-14}+x+2.18159\times {10}^{-14}{x}^2-x^3,\hfill & x\in [\frac{1}{3},\frac{2}{3})\hfill \\ -1.33227\times {10}^{-15}+x-6.75016\times {10}^{-14}{x}^2-x^3,\hfill & x\in [\frac{2}{3},1).\hfill \end{array}\right.\end{array}$$

The approximate solution when $n=10$ is given as

$$\begin{array}{c}\hfill \mathsf{\Phi }\left(x\right)=\left\{\begin{array}{cc}2.90878\times {10}^{-14}+x-3.86358\times {10}^{-13}{x}^2-x^3,\hfill & x\in [0,\frac{1}{10})\hfill \\ 4.03289\times {10}^{-14}+x+2.98606\times {10}^{-12}{x}^2-x^3,\hfill & x\in [\frac{1}{10},\frac{1}{5})\hfill \\ -6.13398\times {10}^{-14}+x-4.61853\times {10}^{-12}{x}^2-x^3,\hfill & x\in [\frac{1}{5},\frac{3}{10})\hfill \\ ⋮\hfill \\ ⋮\hfill \\ -5.08038\times {10}^{-13}+x-3.21165\times {10}^{-12}{x}^2-x^3,\hfill & x\in [\frac{7}{10},\frac{4}{5})\hfill \\ 2.57216\times {10}^{-12}+x+1.13118\times {10}^{-11}{x}^2-x^3,\hfill & x\in [\frac{4}{5},\frac{9}{10})\hfill \\ -2.81375\times {10}^{-12}+x-8.66862\times {10}^{-12}{x}^2-x^3,\hfill & x\in [\frac{9}{10},1).\hfill \end{array}\right.\end{array}$$

Example 3

([27]). Consider the nonlinear FFIDE,

$$D^{\frac{3}{4}}\varphi \left(x\right)=\frac{4}{3}\mathsf{\Gamma }{\left(\frac{3}{4}\right)}^{-1}{x}^{\frac{3}{4}}-\frac{x}{4}+{\int }_0^1\left(xt\right){\left(\varphi \left(t\right)\right)}^{2}dt$$

with

$$\varphi \left(0\right)=0.$$

The analytical solution of this problem is $\varphi \left(x\right)=x$. We applied the proposed scheme to the above example and obtained numerical results. Figure 5 displays a comparison between the approximate and exact solutions for Example 3. The absolute error is also depicted in Figure 6. Additionally,

Table 3. Comparison of absolute errors when $n=13$ for Example 3.

x	Present Method	Method in [27]
0.1	$2.77556\times {10}^{-17}$	$6.53\times {10}^{-16}$
0.2	$4.16334\times {10}^{-16}$	$5.89\times {10}^{-15}$
0.3	$6.66134\times {10}^{-16}$	$5.44\times {10}^{-15}$
0.4	$2.22045\times {10}^{-16}$	$3.41\times {10}^{-15}$
0.5	$5.55112\times {10}^{-16}$	$4.07\times {10}^{-16}$
0.6	0	$2.61\times {10}^{-15}$
0.7	$1.11022\times {10}^{-16}$	$4.63\times {10}^{-15}$
0.8	$1.33227\times {10}^{-15}$	7
0.9	$3.33067\times {10}^{-16}$	$4.98\times {10}^{-15}$

presents a comparison of the absolute errors with those given in [27].

The approximate solution for Example 3 when $n=13$ is given as

$$\begin{array}{c}\hfill \mathsf{\Phi }\left(x\right)=\left\{\begin{array}{cc}-1.73472\times {10}^{-18}+x-8.88178\times {10}^{-16}{x}^2,\hfill & x\in [0,\frac{1}{13})\hfill \\ 1.249\times {10}^{-16}+x+6.39488\times {10}^{-14}{x}^2-2.55795\times {10}^{-13}{x}^3,\hfill & x\in [\frac{1}{13},\frac{2}{13})\hfill \\ -2.72005\times {10}^{-15}+x-2.98428\times {10}^{-13}{x}^2+5.25802\times {10}^{-13}{x}^3,\hfill & x\in [\frac{2}{13},\frac{3}{13})\hfill \\ ⋮\hfill \\ ⋮\hfill \\ 5.11591\times {10}^{-13}+x+2.50111\times {10}^{-12}{x}^2-1.08002\times {10}^{-12}{x}^3,\hfill & x\in [\frac{10}{13},\frac{11}{13})\hfill \\ -4.83169\times {10}^{-13}+x-1.81899\times {10}^{-12}{x}^2+6.82121\times {10}^{-13}{x}^3,\hfill & x\in [\frac{11}{13},\frac{12}{13})\hfill \\ 2.84217\times {10}^{-13}+x+7.95808\times {10}^{-13}{x}^2-3.41061\times {10}^{-13}{x}^3,\hfill & x\in [\frac{12}{13},1).\hfill \end{array}\right.\end{array}$$

The approximate solution when $n=16$ is given as

$$\begin{array}{c}\hfill \mathsf{\Phi }\left(x\right)=\left\{\begin{array}{cc}x,\hfill & \hfill x\in [0,\frac{1}{16})\hfill \\ x,\hfill & x\in [\frac{1}{16},\frac{1}{8})\hfill \\ 5.55112\times {10}^{-17}+x,\hfill & x\in [\frac{1}{8},\frac{3}{16})\hfill \\ ⋮\hfill \\ ⋮\hfill \\ 1-9.09495\times {10}^{-13}x2+2.27374\times {10}^{-13}{x}^3,\hfill & x\in [\frac{13}{16},\frac{7}{8})\hfill \\ 1.13687\times {10}^{-13}+x+9.09495\times {10}^{-13}{x}^2-2.27374\times {10}^{-13}{x}^3,\hfill & x\in [\frac{7}{8},\frac{15}{16})\hfill \\ -1.13687\times {10}^{-13}+x,\hfill & x\in [\frac{15}{16},1).\hfill \end{array}\right.\end{array}$$

Example 4

([27]). Consider the nonlinear FFIDE,

$$D^{\frac{3}{4}}\varphi \left(x\right)=\frac{64}{15}\mathsf{\Gamma }\left(\frac{3}{4}\right)\frac{\sqrt{2}{x}^{\frac{9}{4}}}{\pi }-\frac{1}{8}x+{\int }_0^1\left(xt\right){\left(\varphi \left(t\right)\right)}^{2}dt$$

with

$$\varphi \left(0\right)=0.$$

The analytical solution of this example is $\varphi \left(x\right)=x^3$. We applied the proposed method to the above problem and obtained numerical outcomes. A comparison between the approximate and exact solutions for Example 4 is exhibited in Figure 7. The absolute error is also shown in Figure 8.

Table 4. Comparison of absolute errors for $n=13$ for Example 4.

x	Present Method	Method in [27]
0.1	$1.86483\times {10}^{-17}$	$1.75\times {10}^{-13}$
0.2	$1.31839\times {10}^{-16}$	$4.05\times {10}^{-12}$
0.3	$4.85723\times {10}^{-17}$	$3.75\times {10}^{-12}$
0.4	$4.16334\times {10}^{-17}$	$2.21\times {10}^{-12}$
0.5	$4.85723\times {10}^{-16}$	$8.35\times {10}^{14}$
0.6	$1.38778\times {10}^{-16}$	$2.38\times {10}^{-12}$
0.7	$1.66533\times {10}^{-16}$	$3.92\times {10}^{-12}$
0.8	$3.33067\times {10}^{-16}$	$4.23\times {10}^{-12}$
0.9	$4.44089\times {10}^{-16}$	$3.11\times {10}^{-13}$

reports a comparison between the absolute errors and those mentioned in [27].

The approximate solution for Example 4 when $n=13$ is given as

$$\begin{array}{c}\hfill \mathsf{\Phi }\left(x\right)=\left\{\begin{array}{cc}{x}^3,\hfill & x\in [0,\frac{1}{13})\hfill \\ 1.11022\times {10}^{-16}{x}^2+x^3,\hfill & x\in [\frac{1}{13},\frac{2}{13})\hfill \\ -6.73073\times {10}^{-16}+1.29896\times {10}^{-14}x-8.52651\times {10}^{-14}{x}^2+x^3,\hfill & x\in [\frac{2}{13},\frac{3}{13})\hfill \\ ⋮\hfill \\ ⋮\hfill \\ -2.84217\times {10}^{-14}-1.7053\times {10}^{-13}x+2.84217\times {10}^{-13}{x}^2+x^3,\hfill & x\in [\frac{10}{13},\frac{11}{13})\hfill \\ 1.22213\times {10}^{-12}-4.26326\times {10}^{-12}x+4.77485\times {10}^{-12}{x}^2+x^3,\hfill & x\in [\frac{11}{13},\frac{12}{13})\hfill \\ 1.12834\times {10}^{-11}-3.62661\times {10}^{-11}x+4.01315\times {10}^{-11}{x}^2+x^3,\hfill & x\in [\frac{12}{13},1).\hfill \end{array}\right.\end{array}$$

The approximate solution when $n=8$ is given as

$$\begin{array}{c}\hfill \mathsf{\Phi }\left(x\right)=\left\{\begin{array}{cc}{x}^3,\hfill & x\in [0,\frac{1}{8})\hfill \\ -1.73472\times {10}^{-18}+5.55112\times {10}^{-17}x-4.44089\times {10}^{-16}{x}^2+x^3,\hfill & x\in [\frac{1}{8},\frac{1}{4})\hfill \\ 1.38778\times {10}^{-17}-2.22045\times {10}^{-16}x+8.88178\times {10}^{-16}{x}^2+x^3,\hfill & x\in [\frac{1}{4},\frac{3}{8})\hfill \\ ⋮\hfill \\ ⋮\hfill \\ x^3,\hfill & x\in [\frac{5}{8},\frac{3}{4})\hfill \\ x^3,\hfill & x\in [\frac{3}{4},\frac{7}{8})\hfill \\ x^3,\hfill & x\in [\frac{7}{8},1).\hfill \end{array}\right.\end{array}$$

Example 5

([23]). Consider the VFIDE,

$$D^{\frac{1}{3}}\varphi \left(x\right)=\frac{3x^{\frac{2}{3}}}{2\mathsf{\Gamma }\left(\frac{2}{3}\right)}-1+exp\left(x^2\right)-x^{2}exp\left(x^2\right)+{\int }_0^x{x}^{2}exp\left(xt\right)\varphi \left(t\right)dt$$

with

$$\varphi \left(0\right)=0,{\varphi }^{\prime }\left(0\right)=1.$$

The analytical solution to this problem is $\varphi \left(x\right)=x$. We employed the aforementioned scheme for the problem and obtained numerical outcomes. Figure 9 showcases a comparison between the approximate and exact solutions for Example 5, with the absolute error plotted in Figure 10. Furthermore,

Table 5. Absolute errors when $n=10$ for Example 5.

x	Present Method	Method in [21]
0	$1.24900\times {10}^{-16}$	$1.4247\times {10}^{-5}$
0.1	$2.60277\times {10}^{-12}$	$1.1250\times {10}^{-5}$
0.2	$1.80098\times {10}^{-12}$	$8.2823\times {10}^{-6}$
0.3	$3.91964\times {10}^{-13}$	$5.3691\times {10}^{-6}$
0.4	$1.74677\times {10}^{-12}$	$2.5271\times {10}^{-6}$
0.5	$4.84612\times {10}^{-14}$	$2.4902\times {10}^{-7}$
0.6	$1.04494\times {10}^{-12}$	$2.9990\times {10}^{-6}$
0.7	$1.95177\times {10}^{-13}$	$5.8130\times {10}^{-6}$
0.8	$5.46452\times {10}^{-13}$	$8.8544\times {10}^{-6}$
0.9	$4.91496\times {10}^{-13}$	$1.2387\times {10}^{-5}$
1	$7.51288\times {10}^{-13}$	$1.6811\times {10}^{-5}$

illustrates a comparison of the absolute errors and those reported in [23].

The approximate solution for Example 5 with $n=10$ is as follows:

$$\begin{array}{c}\hfill \mathsf{\Phi }\left(x\right)=\left\{\begin{array}{cc}-1.21431\times {10}^{-16}+x+1.77636\times {10}^{-14}{x}^2+7.09178\times {10}^{-10}{x}^3,\hfill & x\in [0,\frac{1}{10})\hfill \\ 8.50833\times {10}^{-13}+x+2.55303\times {10}^{-10}{x}^2-1.41789\times {10}^{-10}{x}^3,\hfill & x\in [\frac{1}{10},\frac{1}{5})\hfill \\ 1.69607\times {10}^{-11}+x+1.46355\times {10}^{-9}{x}^2-2.1555\times {10}^{-9}{x}^3,\hfill & x\in [\frac{1}{5},\frac{3}{10})\hfill \\ ⋮\hfill \\ ⋮\hfill \\ 2.83592\times {10}^{-10}+x+1.43797\times {10}^{-9}{x}^2-6.18314\times {10}^{-10}{x}^3,\hfill & x\in [\frac{7}{10},\frac{4}{5})\hfill \\ 2.08331\times {10}^{-11}+x+2.05944\times {10}^{-10}{x}^2-1.05032\times {10}^{-10}{x}^3,\hfill & x\in [\frac{4}{5},\frac{9}{10})\hfill \\ -2.42892\times {10}^{-10}+x-7.70342\times {10}^{-10}{x}^2+2.56705\times {10}^{-10}{x}^3,\hfill & x\in [\frac{9}{10},1).\hfill \end{array}\right.\end{array}$$

The approximate solution for $n=5$ is given as

$$\begin{array}{c}\hfill \mathsf{\Phi }\left(x\right)=\left\{\begin{array}{cc}-3.46945\times {10}^{-17}+x-2.66454\times {10}^{-15}{x}^2-1.91847\times {10}^{-13}{x}^3,\hfill & x\in [0,\frac{1}{5})\hfill \\ -1.18516\times {10}^{-14}+x-8.88178\times {10}^{-13}{x}^2+1.28431\times {10}^{-12}{x}^3,\hfill & x\in [\frac{1}{5},\frac{2}{5})\hfill \\ 2.07279\times {10}^{-13}+x+3.22586\times {10}^{-12}{x}^2-2.14229\times {10}^{-12}{x}^3,\hfill & x\in [\frac{2}{5},\frac{3}{5})\hfill \\ -4.3876\times {10}^{-13}+x-2.16716\times {10}^{-12}{x}^2+8.54428\times {10}^{-13}{x}^3,\hfill & x\in [\frac{3}{5},\frac{4}{5})\hfill \\ -1.04805\times {10}^{-13}+x-6.25278\times {10}^{-13}{x}^2+2.07834\times {10}^{-13}{x}^3,\hfill & x\in [\frac{4}{5},1).\hfill \end{array}\right.\end{array}$$

Example 6

([28]). Consider the VFIDE,

$$D^{\frac{5}{6}}\varphi \left(x\right)=\frac{-3\mathsf{\Gamma }\left(\frac{5}{6}\right)x^{\frac{1}{6}}(-91+216x^2)}{91\pi }+5x-xexp\left(x\right)(5-5x+3x^2-x^3)+{\int }_0^{x}xexp\left(t\right)\varphi \left(t\right)dt$$

with $\varphi \left(0\right)=0$.

The exact solution is given by $\varphi \left(x\right)=x-x^3$. We implemented the introduced scheme and obtained approximate solutions. Figure 11 illustrates the comparison between the approximate and exact solutions, with the absolute error depicted in Figure 12.

Table 6. Comparison of absolute errors for various values of h for Example 6.

h	[28] (S1)	[28] (S2)	[28] (S3)	Present Method
$\frac{1}{5}$	$2.73753\times {10}^{-1}$	$6.70569\times {10}^{-2}$	$7.26013\times {10}^{-2}$	$1.4439\times {10}^{-14}$
$\frac{1}{10}$	$1.15363\times {10}^{-1}$	$1.47491\times {10}^{-2}$	$1.61695\times {10}^{-2}$	$2.2417\times {10}^{-13}$
$\frac{1}{20}$	$5.05187\times {10}^{-2}$	$3.28518\times {10}^{-3}$	$3.65192\times {10}^{-3}$	$6.2883\times {10}^{-13}$
$\frac{1}{40}$	$2.24081\times {10}^{-2}$	$7.33189\times {10}^{-4}$	$8.26938\times {10}^{-4}$	$1.2594\times {10}^{-11}$

presents a comparison of the absolute errors and those in [28].

The approximate solution for Example 6 when $n=10$ is as follows:

$$\begin{array}{c}\hfill \mathsf{\Phi }\left(x\right)=\left\{\begin{array}{cc}2.24182\times {10}^{-13}+x+1.23457\times {10}^{-13}{x}^2-x^3,\hfill & x\in [0,\frac{1}{10})\hfill \\ 2.76688\times {10}^{-13}+x+1.58771\times {10}^{-11}{x}^2-x^3,\hfill & x\in [\frac{1}{10},\frac{1}{5})\hfill \\ -4.52582\times {10}^{-13}+x-3.88169\times {10}^{-11}{x}^2-x^3,\hfill & x\in [\frac{1}{5},\frac{3}{10})\hfill \\ ⋮\hfill \\ ⋮\hfill \\ \hfill 1.911\times {10}^{-11}+x+9.29532\times {10}^{-11}{x}^2-x^3,\hfill & x\in [\frac{7}{10},\frac{4}{5})\hfill \\ -2.35403\times {10}^{-11}+x-1.07008\times {10}^{-10}{x}^2-x^3,\hfill & x\in [\frac{4}{5},\frac{9}{10})\hfill \\ -7.70228\times {10}^{-11}+x-3.04937\times {10}^{-10}{x}^2-x^3,\hfill & x\in [\frac{9}{10},1).\hfill \end{array}\right.\end{array}$$

The approximate solution when $n=20$ is given by

$$\begin{array}{c}\hfill \mathsf{\Phi }\left(x\right)=\left\{\begin{array}{cc}3.13083\times {10}^{-14}+x+1.29674\times {10}^{-13}{x}^2-x^3,\hfill & x\in [0,\frac{1}{20})\hfill \\ 7.65464\times {10}^{-14}+x+5.44276\times {10}^{-11}{x}^2-x^3,\hfill & x\in [\frac{1}{20},\frac{1}{10})\hfill \\ -4.84529\times {10}^{-13}+x-1.13921\times {10}^{-10}{x}^2-x^3,\hfill & x\in [\frac{1}{10},\frac{3}{20})\hfill \\ ⋮\hfill \\ ⋮\hfill \\ -7.4283\times {10}^{-10}+x-2.88992\times {10}^{-9}{x}^2-x^3,\hfill & x\in [\frac{17}{20},\frac{9}{10})\hfill \\ 5.02723\times {10}^{-10}+x+1.72315\times {10}^{-9}{x}^2-x^3,\hfill & x\in [\frac{9}{10},\frac{19}{20})\hfill \\ -1.59787\times {10}^{-10}+x-4.79076\times {10}^{-10}{x}^2-x^3,\hfill & x\in [\frac{19}{20},1).\hfill \end{array}\right.\end{array}$$

7. Concluding Remarks

In conclusion, the utilization of the Cubic B-spline collocation method has evolved as a powerful tool for addressing the challenges posed by linear and nonlinear FFIDEs of Volterra and Fredholm types. By offering a durable and systematic technique for approximating solutions to these complex equations, our proposed method has exhibited exceptional accuracy and efficiency. Through systematic mathematical foundations and comparative analyses, we have highlighted the superiority of our approach over existing methods. Furthermore, the stability analysis conducted in this study provides additional confirmation of the precision and reliability of the Cubic B-spline collocation method. Overall, our research contributes to the advancement of computational techniques for solving FIDEs, offering promising solutions for addressing a wide range of practical problems across various fields of science and engineering. For future work, it would be valuable to explore the extension of this method to handle FIDEs with additional complexities, such as those involving time-varying parameters or non-local operators. Investigating the integration of adaptive strategies and refinement techniques could further enhance the method’s robustness and computational efficiency.