An effective method for solving nonlinear integral equations involving the Riemann-Liouville fractional operator

1.   Introduction

Many physical issues in science and technology that were previously solved by using the concepts of partial and ordinary differential equations can now be solved using more effective methods of integral equation theory. The generalization of differential and integral to an arbitrary order that isn't an integer is represented by fractional derivatives and fractional integrals respectively. In theoretical mechanics, applied mathematics, and mathematical physics, the theory of integral equations is a very useful tool ^[1,2,3,4]. In the recent decade, several numerical techniques have been established by various researchers to obtain approximate solutions of integral equations with integer or fractional order. For instance, an iterative numerical method based on Picard iteration has been presented by Micula ^[5] for linear fractional integral equations of second kind. Shiralashetti et al. ^[6] proposed a technique based on Fibonacci wavelets to solve nonlinear Volterra integral equations. Ali et al. ^[7] introduced a method that utilizes the hybrid orthonormal Bernstein and block-pulse functions wavelet method for the nonlinear Volterra integral equations with a weakly singular kernel. Bairwa et al. ^[8] presented a method named as q-homotopy analysis transform merthod for Abel's integral equations of the second kind. Bhat et al. ^[9] presented a method for the numerical solution of the third kind of Volterra integral equations based on Lagrange polynomial, modified Lagrange polynomial, and barycentric Lagrange polynomial. Hamdan et al. introduced the product integration and Haar Wavelet approaches for approximate solutions of the fractional Volterra integral equations of the second type ^[10]. Akg$\ddot{\mbox{u}}$l et al. ^[11] proposed a novel methodology based on reproducing kernels for solving the fractional order integro-differential transport model for a nuclear reactor. Khan et al. ^[12] investigated a competition model in the fractional derivative of the Caputo-Fabrizio type via Newton polynomial. A novel multi-parametric homotopy method for systems of linear and nonlinear Fredholm integral equations has been presented by Khan et al. ^[13]. Also, nice results were obtained for two-dimensional Fredholm integral equations via a multi-parametric homotopy approach ^[14]. More recent works on the solvability of integral equations have been introduced in references ^{[15,16,17,18,19]}.

Regarding the other perspective, stability is a significant factor for numerical applications and may be required to compare the outcome and performance of numerical methods. Both differential and integral equations have undergone several types of stability analysis. Recently, H-U-R and H-U stability have increasingly attained the interest of researchers. For instance, in the vast of practical problem scenarios, Lyapunov stability has been explored. Additionally, for a lot of topics, Mittag-Leffler and exponential type stabilities have been developed. Practical exponential stability of an impulsive stochastic food chain system with time-varying delays has been studied ^[20]. Ali et al. ^[21] studied on the Ulam-Hyers stability and the Ulam-Hyers-Rassias stability for nonlinear implicit fractional order differential equations. The Ulam-Hyers-Rassias stability of a nonlinear stochastic integral equation has been studied by Ngoc et al. ^[22]. Kumam et al. ^[23] developed some conditions on existence results and Hyers-Ulam stability for a class of nonlinear fractional order differential equations. Reinfelds et al. ^[24] presented the Hyers-Ulam stability of a nonlinear Volterra integral equation on time scales. Morales et al. ^[25] studied Hyers-Ulam and Hyers-Ulam-Rassias stability for nonlinear integral equations with delay. Subramanian et al. ^[26] studied the existence and U-H stability results for nonlinear coupled fractional differential equations with boundary conditions involving Riemann-Liouville and Erdélyi-Kober integrals. Li et al. ^[27] studied the Ulam stability of a fractional differential equation with multi-point boundary conditions and non-instantaneous integral impulse. Some examples and H-U-R stability for Volterra integral equations within weighted spaces have been presented by Castro et al. ^[28].

The existence requirements of the solutions to various types of integral equations are crucial elements to examine. We can use these criteria to determine under what conditions the solution to the problem exists. The fixed point approach is significant in this sense.

In this paper, our discussion will be on the nonlinear integral equation involving the Riemann-Liouville fractional operator, i.e.,

where $\lambda > 0$, $0 < \beta < \infty$, the functions $\mathcal{H}:C([0, \beta], \mathbb{R})\rightarrow \mathbb{R}$ and $\mathcal{N}, \psi, \xi :[0, \beta]\rightarrow \mathbb{R}$ are all continuous functions.

Sometimes, it is challenging to determine the exact solution to this equation explicitly. Thus, the current article's goal is to use the Laguerre polynomials to obtain approximate numerical solutions for nonlinear integral equations involving the RLFO. Remarkably, nonlinearity and singularity make the numerical procedure more challenging, so this is the strong motivation to study the numerical solution of Eq (1.1). In order to achieve this goal, we provide some conditions for the existence and uniqueness of solutions as well as a numerical scheme to solve the Eq (1.1). Moreover, we establish the conditions of Hyers-Ulam and Hyers-Ulam-Rassias stability for the considered Eq (1.1).

This paper is structured as follows: Section 2 contains some notations, and background information. In Section 3, we provides existence, uniqueness and stability criteria to the solutions of the Eq (1.1). The required method for solving the Eq (1.1) to obtain the approximate numerical solutions is presented in Section 4. In Section 5, accuracy of solutions of the proposed method is discussed. Some numerical experiments and graphs are shown in Section 6 to illustrate the effectiveness of the proposed technique. Finally, conclusions and future work are given in Section 7.

2.   Notations, definitions and auxiliary facts

In this section, we provide a few notations and background information which are useful for the presentation of our main results. Let $\mathbb{R}$ be the set of all real numbers and let $W = C([0, \beta], \mathbb{R})$ be the space of all continuous functions $\zeta: [0, \beta]\rightarrow \mathbb{R}$, where $0 < \beta < \infty$. Then $(W, \|.\|)$ is a Banach space with the norm $\|\zeta\| = \sup\{\lvert \zeta(y)\rvert : y\in[0, \beta]\}$.

Definition 2.1. ^[2] The Riemann-Liouville fractional integral of order $\lambda > 0$ of a function $\zeta(y)$ is described as

where $\Gamma(\lambda) = \int_{0}^{\infty}e^{-x}x^{\lambda-1}dx$, provided that the RHS is point-wise defined on $[0, \infty)$.

Definition 2.2. If for each function $\zeta(y)\in W$ satisfying

where $v(y)$ is a non-negative function, there exists a solution $\hat{\zeta}\in W$ of Eq (1.1) and a constant $\Lambda > 0$ independent of $\zeta$ and $\hat{\zeta}$ with the property $\lvert \zeta(y)-\hat{\zeta}(y)\rvert\leq\Lambda v(y)$, for all $y$, then we say that the Eq (1.1) has the H-U-R stability with respect to $v(y)$.

Definition 2.3. In particular, when $v(y)$ is a constant function in the above inequalities of Definition 2.2., we say that the Eq (1.1) has the H-U stability.

Theorem 2.1. (Banach's fixed point theorem ^[29]). Let $W$ be a Banach space. If a mapping $\mathcal{M} : W\rightarrow W$ is a contraction, then $\mathcal{M}$ has a unique fixed point in $W$.

Theorem 2.2. (Leray-Schauder alternative ^[26]). Let $\mathcal{M}:\mathcal{E}\rightarrow \mathcal{E}$ be a completely continuous operator. Let $\Omega = \{y\in\mathcal{E}:y = \eta\mathcal{M}(y),$ for some $0 < \eta < 1\}$. Then, either the set $\Omega$ is unbounded or $\mathcal{M}$ has at least one fixed point.

2.1.   Laguerre polynomials

The Laguerre polynomials $L_n(y) \; (n\ge 0)$, (^{[30,31,32,33]}) are defined by the generating function

Indeed, $u = L_n(y)$ is a solution of the following differential equation

From Eq (2.2),

Thus by Eq (2.3), we get

Moreover, one can also define the Laguerre polynomials by the recurrence relation as follows:

The first six Laguerre polynomials are given below:

$L_0(y) = 1$,

$L_1(y) = 1-y$,

$L_2(y) = \frac{1}{2}(2-4y+y^2)$,

$L_3(y) = \frac{1}{6}(6-18y+9y^2-y^3)$,

$L_4(y) = \frac{1}{24}(24-96y+72y^2-16y^3+y^4)$,

$L_5(y) = \frac{1}{120}(120-600y+600y^2-200y^3+25y^4-y^5)$.

2.2.   Approximation function and matrix representation for Laguerre polynomials

To determine the approximate numerical solutions of Eq (1.1), assume that the function $J_n(y)$ is approximated by the Laguerre polynomials as follows:

where the function $\left\{L_m(y)\right\}^{n}_{m = 0}$ denotes the basis of the Laguerre polynomials, defined in Eq (2.4). We have to determine the unknown Laguerre coefficients $s_m$, $(m = 0, 1, ..., n)$.

Now rewriting Eq (2.5) as,

Again, Eq (2.6) can be converted as,

where $\left\{\sigma_{mm}\right\}^{n}_{m = 0}$ are the parameters of the power basis, used to obtain the Laguerre polynomials, this matrix is upper triangular and is certainly invertible.

Now for $n = 1$, $2$ and $3$, the operational matrices are shown in the Eqs (2.8), (2.9) and (2.10) respectively as follows:

3.   Qualitative analysis

This section covers the existence and uniqueness of solutions of the Eq (1.1) via fixed point approach, as well as the H-U-R and H-U stability theorems are also given for the Eq (1.1).

Theorem 3.1. We consider the assumptions for Eq (1.1) as follows:

$(A_1)$ There exist constants $B_1 > 0$, $B_2 > 0$ and $B_3 > 0$ such that

$(A_2)$ There exists a constant $B > 0$ such that $\lvert\mathcal{H}(\zeta_1)-\mathcal{H}(\zeta_2)\rvert\leq B\lvert\zeta_1- \zeta_2 \rvert$, for all $\zeta_1, \zeta_2\in W$.

$(A_3)$ There exists a constant $H > 0$ such that $\lvert\mathcal{H}(\zeta(\mu))\rvert\leq H$, for all $\zeta(\mu)\in W$.

Then under assumptions $(A_1)-(A_3)$, Eq (1.1) has at least one solution defined on $[0, \beta]$.

Proof. Let us define the operator $\mathcal{M}:W\rightarrow W$ as

We will prove this result in four steps as follows:

Step 1. We need to prove that $\mathcal{M}$ is continuous.

Let $\{\zeta_n\}$ be a sequence in $W$ and $\zeta\in W$ such that $\lim_{n\rightarrow \infty}\|\zeta_n-\zeta\| = 0$.

Then for each $y\in[0, \beta]$,

which implies that $\|\mathcal{M}(\zeta_n)-\mathcal{M}(\zeta)\|\rightarrow 0$ as $n\rightarrow \infty$. Hence, $\mathcal{M}$ is continuous.

Step 2. Under the mapping $\mathcal{M}$, bounded sets of $W$ are mapped into bounded sets of $W$.

Let $B_\epsilon\subset W$ be bounded. Now, for $\zeta\in B_\epsilon$ and $y\in [0, \beta]$,

Step 3. $\mathcal{M}(B_\epsilon)$ is equicontinuous.

Let $\zeta\in B_\epsilon$ and $y_1, y_2\in[0, \beta]$ with $y_1 < y_2$, then

i.e.,

Now, it can be observed that

Then,

Thus, it follows from (3.2) and Eq (3.4) that, $\lvert\mathcal{M}(\zeta(y_2))-\mathcal{M}(\zeta(y_1))\rvert\rightarrow 0$ as $y_2\rightarrow y_1$. So, $\mathcal{M}(B_\epsilon)$ is equicontinuous.

Hence, combining Step 1, Step 2 and Step 3, the operator $\mathcal{M}$ is completely continuous by the consequence of Arzel$\grave{\mbox{a}}$-Ascoli theorem.

Step 4. Let $\Omega = \{\zeta\in W:\zeta = \eta\mathcal{M}(\zeta),$ for some $0 < \eta < 1\}$. We have to prove that the set $\Omega$ is bounded. Let $\zeta\in\Omega$, which implies that $\zeta = \eta\mathcal{M}(\zeta)$, for some $0 < \eta < 1$.

Now for $y\in [0, \beta]$, we have

This implies that the set $\Omega$ is bounded.

Thus, $\mathcal{M}$ has at least one fixed point according to the Leray-Schauder alternative, which is a solution of Eq (1.1) defined on $[0, \beta]$. □

Theorem 3.2. We consider the assumptions for Eq (1.1) as follows:

$(A_1)$ There exist constants $B_2 > 0$ and $B_3 > 0$ such that

$(A_2)$ There exists a constant $B$ with $0 < B < \frac{\Gamma(\lambda+1)}{B_{2}B_{3}\beta^{\lambda}}$ such that

$\lvert\mathcal{H}(\zeta_1)-\mathcal{H}(\zeta_2)\rvert\leq B\lvert\zeta_1- \zeta_2 \rvert$, for all $\zeta_1, \zeta_2\in W$.

Then under assumptions $(A_1)$ and $(A_2)$, Eq (1.1) has a unique solution defined on $[0, \beta]$.

Proof. Let us define the operator $\mathcal{M}:W\rightarrow W$ as

Now we have to prove that $\mathcal{M}$ is a contraction.

Let $\zeta_1, \zeta_2\in W$, then for all $y\in[0, \beta]$,

Thus, $\|\mathcal{M}(\zeta_1(y))-\mathcal{M}(\zeta_2(y))\|\leq\kappa\|\zeta_1-\zeta_2\|$, where $\kappa = \frac{B_{2}B_{3}B\beta^{\lambda}}{\Gamma(\lambda+1)}$.

By the assumption of $B$, $\kappa = \frac{B_{2}B_{3}B\beta^{\lambda}}{\Gamma(\lambda+1)} < 1$. So $\mathcal{M}$ is a contraction mapping. As a result, $\mathcal{M}$ has a unique fixed point, according to the Banach's fixed point theorem and hence, Eq (1.1) has a unique solution. □

Theorem 3.3. Assume that, the Eq (1.1) satisfying all assumptions of Theorem 3.2. Let $\zeta(y)\in W$ is such that satisfies the inequality

where $v(y)$ is a non-negative function. Then the Eq (1.1) has the H-U-R stability with respect to $v(y)$.

Proof. By the Theorem 3.2, $\exists$ a unique solution $\hat{\zeta}\in W$ of Eq (1.1).

Let $\zeta(y)\in W$ such that

According to Definition 2.2., we have to show that $\exists$ a constant $\Lambda > 0$ independent of $\zeta$ and $\hat{\zeta}$ such that, $\vert \zeta(y)-\hat{\zeta}(y)\rvert\leq\Lambda v(y)\rvert$. Now, for all $y\in[0, \beta]$,

Thus, $\lvert \zeta(y)-\hat{\zeta}(y)\rvert\leq\|\zeta-\hat{\zeta}\|\leq v(y)+\frac{B_2B_3B\beta^{\lambda}}{\Gamma(\lambda+1)} \|\zeta-\hat{\zeta}\|$, which implies that, $\lvert \zeta(y)-\hat{\zeta}(y)\rvert\leq\|\zeta-\hat{\zeta}\|\leq \Lambda v(y)$, where $\Lambda = \frac{1}{1-\kappa} > 0$, as $\kappa = \frac{B_{2}B_{3}B\beta^{\lambda}}{\Gamma(\lambda+1)} < 1$, by the assumption. Hence the Eq (1.1) has the H-U-R stability with respect to $v(y)$. □

Theorem 3.4. Assume that, the Eq (1.1) satisfying all assumptions of Theorem 3.2. Let $\zeta(y)\in W$ and $\epsilon > 0$ such that satisfies

then the Eq (1.1) has the H-U stability with respect to $\epsilon$.

Proof. We can prove this theorem in similar way by taking $v(y) = \epsilon$, in the Theorem 3.3, where $\epsilon > 0$. □

4.   Computational method

In this section, we will use the Laguerre polynomials to determine the approximate numerical solutions for the Eq (1.1).

Rewriting the Eq (1.1) as follows:

Now to approximate the unknown function in Eq (4.1), by using Eq (2.5), let

Putting the Eq (4.2) into the Eq (4.1), we get

Therefore, the Eq (4.3) becomes by using Eq (2.6) as follows:

Again, by applying Eq (2.7), the Eq (4.4) becomes as follows:

Thus, after finding the integration of the Eq (4.5) we have to calculate the values of unknown constants $s_m \; (m = 0, 1, 2, \dots, n)$, for this purpose we need $(n+1)$ equations. Now by choosing $y_i \; (i = 0, 1, 2, \dots, n)$ such that $0\leq y_i\leq\beta$, a system of $(n+1)$ equations can be obtained. After solving these equations the unknown coefficients $(s_0, s_1, \dots, s_n)$ are uniquely determined. Therefore, the approximate numerical solutions can be obtained by substituting the values of the coefficients into Eq (2.5).

4.1.   Algorithm for solutions

In this part, all steps of the proposed algorithm to find the approximate numerical solutions of Eq (1.1) are summarized as follows:

Step 1. Choose a value $n$ for $J_n(y)$ and then obtain the corresponding Laguerre polynomials

Step 2. Now we have to use Eqs (1.1), (2.7) and (4.2).

Step 3. Substitute Eq (4.2) into Eq (1.1) and use Eq (2.7).

Step 4. Compute all of the integration obtained in Step 3.

Step 5. Compute $(s_0, s_1, s_2, \dots, s_n)$, by choosing $y_i\in[0, \beta]$, where $i = 0, 1, 2, \dots, n$.

Step 6. Substitute the values of $(s_0, s_1, s_2, \dots, s_n)$ into the Eq (2.5) to get the approximate solution.

5.   Accuracy of solutions

Since Eq (4.3) and Eq (2.5) have the following forms given by Eq (5.1) and Eq (5.2) respectively:

and the unknown coefficients $(s_0, s_1, \dots, s_n)$ were determined by using Eq (4.5). Also, from Eq (4.2), we have

thus, Eq (5.3) is the unique approximate solution of Eq (5.1), and is substituted into Eq (5.1).

Now, assume that $y = y_k\in [0, \beta], k = 0, 1, 2, \dots$, then

and $\Omega(y_k)\leq 10^{-N_{p}}$, where $N_p$ is a positive integer. If $\max 10^{-N_{p}} = 10^{-N}$ ($N$ is a positive integer) prescribed, then $n$ is increased until the difference $\Omega(y_{k})$ at each of the points becomes smaller than the prescribed $10^{-N}$. For $\max 10^{-N_{p}}\neq 10^{-N}$, the error can be estimated by the following function:

If $\Omega_n(y)\rightarrow 0$, for sufficiently large $n$, then the error decreases.

6.   Numerical applications and discussions

In this section, five examples have been provided to demonstrate the efficiency and accuracy of the presented technique. The absolute errors are the values of $\lvert \zeta(y)-J_n(y)\rvert$ at selected points, where $\zeta(y)$ is the exact solution and the approximate solution is $J_n(y)$. Also, the convergence of solutions have been shown in the graphs.

Example 6.1. ^[34,35]. Consider the example as follows:

where the exact solution is $\zeta(y) = y$.

We can write this Eq (6.1) in the form of Eq (1.1) as,

Then by comparing with Eq (1.1), we get $\mathcal{N}(y) = y+\frac{4}{3}y^{\frac{3}{2}}$, $\psi(y) = -\Gamma(\frac{1}{2})$, $\lambda = \frac{1}{2}$, $\beta = 1$, $\xi(\mu) = 1$ and $\mathcal{H}(\zeta(\mu)) = \zeta(\mu)$. It's very easy to show that Eq (6.1) satisfies all assumption of Theorem 3.1. Hence, by Theorem 3.1, Eq (6.1) has at least one solution defined on $[0, 1]$.

Now applying the algorithm as, if $n = 1$, then by using Eq (2.8) and Eq (4.2) we have

this can be written as,

Therefore from Eq (6.3), after computing all integrations, we choose $y_{0} = 0.1$, $y_{1} = 0.2$ in the interval $[0, 1]$, and solving the algebraic system in MATLAB program, we get $s_0 = 1$ and $s_1 = -1$. Thus, the approximate solution for the Eq (6.1) is as follows:

This shows that the approximate solution of Eq (6.1) is same as the exact solution for $n = 1$, by this proposed method.

In ^[34], the maximum absolute error obtained for $n = 18$, nearly $5.0e-03$. In ^[35], the maximum absolute error obtained for $n = 14$, nearly $1.0e-10$. Therefore, our proposed method provides better result than the method given in ^[34,35]. The comparison of the approximate solution with the exact solution for the presented method is shown in Figure 1.

Also, for Eq (6.1), we see that $B_2 = \Gamma(\frac{1}{2})$, $B_3 = 1$, $\beta = 1$, $\lambda = \frac{1}{2}$, and for all $\zeta_1, \zeta_2\in W$, $\lvert\mathcal{H}(\zeta_1)-\mathcal{H}(\zeta_2)\rvert = \lvert\zeta_1-\zeta_2\rvert$, which implies $B = 1$. Now, $\frac{\Gamma(\lambda+1)}{B_{2}B_{3}\beta^{\lambda}} = 0.5$. So, the relation $0 < B < \frac{\Gamma(\lambda+1)}{B_{2}B_{3}\beta^{\lambda}}$ does not hold. Therefore, we can't say the H-U-R stability in this instance on the whole interval $[0, 1]$.

Now we will check the H-U-R stability on the subinterval $[0, 0.2]$, which contains the collocation points $y_0 = 0.1$ and $y_1 = 0.2$. Clearly, in this case, $\beta = 0.2$. Then, $\frac{\Gamma(\lambda+1)}{B_{2}B_{3}\beta^{\lambda}}\approx 1.1180$. So, the relation $0 < B < \frac{\Gamma(\lambda+1)}{B_{2}B_{3}\beta^{\lambda}}$ holds for $\beta = 0.2$. Thus, all conditions of Theorem 3.2 are satisfied on the interval [0, 0.2].

Now, if we simply choose $\zeta(y) = 0$, then

Therefore, by Theorem 3.3, we obviously have

Hence, we conclude that the Eq (6.1) has the H-U-R stability on $[0, 0.2]$ with respect to $v(y) = \frac{7y}{3}$.

Example 6.2. ^[5] Consider the example as follows:

where the exact solution is $\zeta(y) = y\Gamma(2/3)$.

Comparing with Eq (1.1), we get $\mathcal{N}(y) = y\Gamma(2/3)-\frac{1}{40}y^{8/3}$, $\psi(y) = \frac{1}{27}$, $\xi(\mu) = \mu$, $\mathcal{H}(\zeta(\mu)) = \zeta(\mu)$, $\lambda = \frac{2}{3}$, and $\beta = 1$. It's very easy to show that Eq (6.4) satisfies all assumption of Theorem 3.1. Hence, by Theorem 3.1, Eq (6.4) has at least one solution defined on $[0, 1]$. Now to find the approximate numerical solution we applying the proposed algorithm for $n = 1$, and by choosing $y_0 = 0.1$, $y_1 = 0.2$ in the given interval $[0, 1]$, then computing all steps, we get the approximate solution of Eq (6.4) as follows:

Table 1 shows the comparison of approximate, exact solutions and absolute errors obtained for $n = 1$. By using the proposed method, the maximum absolute error obtained for $n = 1$ is $1.7939e-05$. In ^[5], the maximum absolute error obtained for $n = 5$ ($5$ iterations) is $2.871932e-05$.

Thus, the proposed method described in Section 4, provides best approximation for less number of iterations compared to the method given in ^[5]. Moreover, Figure 2 shows the graph of approximate and exact solutions for $n = 1$.

Also, for Eq (6.4), we see that $B_2 = \frac{1}{27}$, $B_3 = 1$, $\beta = 1$, $\lambda = \frac{2}{3}$, and for all $\zeta_1, \zeta_2\in W$, $\lvert\mathcal{H}(\zeta_1)-\mathcal{H}(\zeta_2)\rvert = \lvert\zeta_1-\zeta_2\rvert$, which implies $B = 1$. Now, $\frac{\Gamma(\lambda+1)}{B_{2}B_{3}\beta^{\lambda}}\approx 24.3741$. So, the relation $0 < B < \frac{\Gamma(\lambda+1)}{B_{2}B_{3}\beta^{\lambda}}$ holds. Thus, all conditions of Theorem 3.2 are satisfied on the interval $[0, 1]$. Now, if we simply choose $\zeta(y) = 0$, then

Then, by Theorem 3.4, we obviously have

Hence, we conclude that the Eq (6.4) has the H-U stability on $[0, 1]$ with respect to $\epsilon = 1.3291$.

Example 6.3. ^[7,35] Consider the example as follows:

where the exact solution is $\zeta(y) = y^2$.

We can write this Eq (6.5) in the form of Eq (1.1) as,

Then by comparing with Eq (1.1), we get $\mathcal{N}(y) = y^2+\frac{16}{15}y^{5/2}$, $\psi(y) = -\Gamma(\frac{1}{2})$, $\xi(\mu) = 1$, $\mathcal{H}(\zeta(\mu)) = \zeta(\mu)$, $\lambda = \frac{1}{2}$ and $\beta = 1$. It's very easy to show that Eq (6.5) satisfies all assumption of Theorem 3.1. Hence, by Theorem 3.1, Eq (6.5) has at least one solution defined on $[0, 1]$.

Now to find the approximate numerical solution we applying the proposed algorithm for $n = 2$, and by choosing $y_0 = 0.1$, $y_1 = 0.15$, $y_2 = 0.2$ in the given interval $[0, 1]$, then computing all steps, we get the approximate solution of Eq (6.5) as follows:

This shows that the approximate solution of Eq (6.5) is same as the exact solution for $n = 2$, by this proposed method. In ^[7], the maximum absolute error found for $k = 6, M = 8$, i.e., for $256$ collocation points is $4.82e-08$. In ^[35], for $n = 14$, the maximum absolute error of order $1.0e-11$ was obtained, while we got the exact solution for $n = 2$, i.e., for $3$ collocation points only. Therefore, our proposed method described in Section 4 is better than ^[7,35]. Moreover, Figure 3 shows the graph of approximate and exact solutions for $n = 2$.

Also, for Eq (6.5), we see that $B_2 = \Gamma(\frac{1}{2})$, $B_3 = 1$, $\beta = 1$, $\lambda = \frac{1}{2}$, and for all $\zeta_1, \zeta_2\in W$, $\lvert\mathcal{H}(\zeta_1)-\mathcal{H}(\zeta_2)\rvert = \lvert\zeta_1-\zeta_2\rvert$, which implies $B = 1$. Now, $\frac{\Gamma(\lambda+1)}{B_{2}B_{3}\beta^{\lambda}} = 0.5$. So, the relation $0 < B < \frac{\Gamma(\lambda+1)}{B_{2}B_{3}\beta^{\lambda}}$ does not hold. Therefore, we can't say the H-U-R stability in this instance on the whole interval $[0, 1]$.

Now we will check the H-U-R stability on the subinterval $[0, 0.2]$, which contains the collocation points $y_0 = 0.1$, $y_1 = 0.15$, and $y_2 = 0.2$. Clearly, in this case $\beta = 0.2$. Then,

So, the relation $0 < B < \frac{\Gamma(\lambda+1)}{B_{2}B_{3}\beta^{\lambda}}$ holds for $\beta = 0.2$. Thus, all conditions of Theorem 3.2 are satisfied on the interval $[0, 0.2]$. Now, if we simply choose $\zeta(y) = 0$, then

Hence, by Theorem 3.3, we conclude that the Eq (6.5) has the H-U-R stability on $[0, 0.2]$ with respect to $v(y) = \frac{31y^2}{15}$.

Example 6.4. Consider the example as follows:

where $\zeta(y) = y^2+y$ is the exact solution.

Comparing with Eq (1.1), we get $\mathcal{N}(y) = y^2+y-\frac{\sin y}{3465\cdot\Gamma(\frac{5}{2})}\left(96 y^{\frac{11}{2}}+176 y^{\frac{9}{2}}\right)$, $\psi(y) = \sin y$, $\xi(\mu) = \mu$, $\mathcal{H}(\zeta(\mu)) = \zeta(\mu)$, $\lambda = \frac{5}{2}$, and $\beta = 1$. It's very easy to show that Eq (6.7) satisfies all assumption of Theorem 3.1. Hence, by Theorem 3.1, Eq (6.7) has at least one solution defined on $[0, 1]$.

Now to find the approximate numerical solution we applying the proposed algorithm for $n = 2$, and by choosing $y_0 = 0.1$, $y_1 = 0.2$, $y_2 = 0.3$ in the given interval $[0, 1]$, then computing all steps, we get the approximate solution of Eq (6.7) as follows:

This shows that the approximate solution of Eq (6.7) is same as the exact solution for $n = 2$, by this proposed method. Figure 4 shows the graph of approximate and exact solutions for $n = 2$.

Also, for Eq (6.7), we see that $B_2 = 1$, $B_3 = 1$, $\beta = 1$, $\lambda = \frac{5}{2}$, and for all $\zeta_1, \zeta_2\in W$, $\lvert\mathcal{H}(\zeta_1)-\mathcal{H}(\zeta_2)\rvert = \lvert\zeta_1-\zeta_2\rvert$, which implies $B = 1$. Now, $\frac{\Gamma(\lambda+1)}{B_{2}B_{3}\beta^{\lambda}}\approx3.3234$. So, the relation $0 < B < \frac{\Gamma(\lambda+1)}{B_{2}B_{3}\beta^{\lambda}}$ holds. Thus, all conditions of Theorem 3.2 are satisfied on the interval $[0, 1]$. Now, if we simply choose $\zeta(y) = 0$, then

Hence, by Theorem 3.3, we conclude that the Eq (6.7) has the H-U-R stability on $[0, 1]$ with respect to $v(y) = y^2+y$.

Example 6.5. ^[6] Consider the example as follows:

where the exact solution is $\zeta(y) = 1+3y$.

Comparing with Eq (1.1), we get $\mathcal{N}(y) = 1+3y-\frac{1}{2}y^2-y^3-\frac{3}{4}y^4$, $\psi(y) = \Gamma(2)$, $\xi(\mu) = 1$, $\mathcal{H}(\zeta(\mu)) = \zeta^2(\mu)$, $\lambda = 2$, and $0 < \beta < 1$. It's very easy to show that Eq (6.8) satisfies all assumption of Theorem 3.1. Thus, by Theorem 3.1, Eq (6.8) has at least one solution defined on $[0, \beta]$, where $0 < \beta < 1$.

Now to find the approximate numerical solution we applying the proposed algorithm for $n = 1$, and by choosing $y_0 = 0.1$, $y_1 = 0.2$ in the given interval $[0, \; \beta]$, then computing all steps, we get the approximate solution of Eq (6.8) as follows:

This shows that the approximate solution of Eq (6.8) is same as the exact solution for $n = 1$, by this proposed method.

In ^[6], the maximum absolute error found for $\hat{m} = 6$, i.e., for $6$ collocation points is $1.1424e-02$, while we got exact solution for $n = 1$, i.e., for $2$ collocation points only. Therefore, our proposed method described in Section 4 is better than the method given in ^[6]. Figure 5 shows the graph of exact and approximate solutions for $n = 1$.

Also, for Eq (6.8), we see that $B_2 = 1$, $B_3 = 1$, $\lambda = 2$, and for all $\zeta_1, \zeta_2\in W$, $\lvert\mathcal{H}(\zeta_1)-\mathcal{H}(\zeta_2)\rvert\leq B\lvert\zeta_1-\zeta_2\rvert$. As $0 < \beta < 1$, we choose the interval $[0, 0.9]$, which contains the collocation points $y_0 = 0.1$ and $y_1 = 0.2$. Now, $\frac{\Gamma(\lambda+1)}{B_{2}B_{3}\beta^{\lambda}}\approx 2.4691$. Thus, for those values of $B$ with $0 < B < 2.4691$ such that $\lvert\mathcal{H}(\zeta_1)-\mathcal{H}(\zeta_2)\rvert\leq B\lvert\zeta_1-\zeta_2\rvert$ holds, for all $\zeta_1, \zeta_2\in W$, then all conditions of Theorem 3.2 are satisfied on the interval $[0, 0.9]$.

Therefore, if we choose $\zeta(y) = 2$, then

Hence, by Theorem 3.3, we conclude that the Eq (6.8) has the H-U-R stability on $[0, 0.9]$ with respect to $v(y) = 1+y$.

7.   Conclusions and future work

In this study, an effective method based on Laguerre polynomials has been proposed to get approximate numerical solutions of Eq (1.1). We established the criteria for the existence and uniqueness of the solutions of Eq (1.1) by using two important theorems, namely, Banach's fixed point theorem and the Leray-Schauder alternative. Also, we established the conditions of H-U-R and H-U stability for the solutions of the considered Eq (1.1). Furthermore, we have provided five examples and compared our computational findings with the exact solutions and with the solutions obtained by some other numerical methods are given in Refs. ^{[5,6,7,34,35]}. For the comparison of numerical and exact solutions, we have given all relevant figures for all examples. According to these study, we see that our proposed method described in Section 4 has best accuracy and effectiveness in comparison to some referenced methods.

More venues can be studied in the future, including singular (Fredholm or Volterra) integral equations with kernels that are more complicated, such as kernels with modified (or delayed) inputs, or kernels satisfy requirements different from those utilized in this work.

Use of AI tools declaration

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

Conflict of interest

The authors declare no conflict of interest.