Hybrid Functions Approach via Nonlinear Integral Equations with Symmetric and Nonsymmetrical Kernel in Two Dimensions

Abstract

The second kind of two-dimensional nonlinear integral equation (NIE) with symmetric and nonsymmetrical kernel is solved in the Banach space $L_2[0,1]\times L_2[0,1]$. Here, the NIE’s existence and singular solution are described in this passage. Additionally, we use a numerical strategy that uses hybrid and block-pulse functions to obtain the approximate solution of the NIE in a two-dimensional problem. For this aim, the two-dimensional NIE will be reduced to a system of nonlinear algebraic equations (SNAEs). Then, the SNAEs can be solved numerically. This study focuses on showing the convergence analysis for the numerical approach and generating an estimate of the error. Examples are presented to prove the efficiency of the approach.

1. Introduction

Integral equations are used in many disciplines of applied mathematics to explore and solve problems. See [1,2,3,4,5,6] for more information on the topic of two-dimensional nonlinear integral equations, which have long been of growing interest in many fields, including medicine [7], biology [8], physics [9], geography and fuzzy control [10]. According to the references [11,12,13,14,15,16], many problems in engineering [17], applied mathematics and mathematical physics [18] can be reduced to two-dimensional nonlinear integral equations with a symmetric and nonsymmetrical kernel. The analytical solutions for these equations are typically difficult. Therefore, it is necessary to use numerical methods or semianalytical methods to obtain the solutions numerically. For example, Bernstein polynomial hybrids with functions of block-pulse form [19,20] and Legendre hybrids with functions of block-pulse form [21,22] have both recently been examined as computational approaches for solving two-dimensional nonlinear integral equations. Alhazmi et al., in [23], used the Lerch polynomial method for solving mixed integral equations in position and time with a strongly symmetric singular kernel. In [24], Azeem et al. conducted research by using the fractional derivative as a treatment for a partial cancer stem cell model. In [25], Attia et al. presented numerical solutions for the fractional differential equations, while in [26], Liaqat et al. presented Shehu transform and the Adomian decomposition technique in a novel algorithm form to establish approximate and exact solutions to quantum mechanics models. Electrical engineering was originally introduced to block-pulse functions by Harmuth; also, other academics have discussed the topic [27].

Recently, hybrid functions have been considered for solving numerous mathematical models, including [28,29,30]. Using block-pulse functions and Legendre polynomials, Maleknejad and Hashemizadeh described a method for solving mixed-type Hammerstein integral equations [31]. Also, hybrid functions have been applied for solving nonlinear Fredholm–Hammerstein systems. Hesameddini et al., in [32,33], obtained a numerical solution to partial differential equations with nonlocal integration terms. Rafiei et al., in [34] obtained the optimal solution of linear time-delay systems. In [35], Ray and Singh discussed the solution of stochastic Volterra–Fredholm integral equations, numerically using hybrid functions. Hashemzadeh and Maleknejad included the necessary definitions as well as some properties of Legendre polynomials and hybrid block-pulse functions [21].

The focus of this research is how to obtain the numerical solution of nonlinear integral equations with a symmetric and nonsymmetrical kernel in two dimensions. By applying special conditions, we utilize the Banach fixed point theorem to establish the existence and uniqueness of a solution for these equations [36,37,38,39,40]. Our investigation delves into the characteristics of hybrid functions, which are a combination of block-pulse functions and Legendre polynomials. We employ these hybrid functions to solve the integral equations, taking advantage of their useful properties. The approach involves transforming the integral equation into a system of algebraic equations, simplifying and improving the process for finding solutions. Overall, our study provides a novel technique for solving two-dimensional nonlinear integral equations with a symmetric and nonsymmetrical kernel using hybrid functions and demonstrates its effectiveness in obtaining accurate solutions.

The present study is structured into six sections, each focusing on a specific aspect of the numerical approach for solving a two-dimensional nonlinear integral equation with a symmetric and nonsymmetrical kernel. Section 2 delves into the existence and uniqueness of the solution for Equation (1), while Section 3 outlines a method for estimating the solution to this equation. In Section 4, we provide a detailed analysis of the convergence properties of our proposed method. The numerical results obtained from applying this approach are presented in Section 5, followed by concluding remarks in Section 6.

This study aims to present a numerical approach for solving the following two-dimensional nonlinear integral equation with a symmetric and nonsymmetrical kernel approximatively:

$$\begin{array}{ll}\gamma \psi (x,y)& =f(x,y)+{\lambda }_1{\int }_0^1{\int }_0^1\Phi (x,\tau ;y,v)\mu (\tau ,v,\psi (\tau ,v))\mathrm{d}v\mathrm{d}\tau \\ & +{\lambda }_2{\int }_0^x{\int }_0^{1}G(x,\tau ;y,v)\nu (\tau ,v,\psi (\tau ,v))\mathrm{d}v\mathrm{d}\tau ,\end{array}$$

(1)

where ${\lambda }_1$ and ${\lambda }_2$ are constant scalers having several physical meanings; the function $\psi (x,y)$ is unknown in the Banach spaces $L_2[0,1]\times L_2[0,1].$ The kernels $\Phi (x,\tau ;y,v)$ and $G(x,\tau ;y,v)$ are continuous in the same space, and the known function $f(x,y)$ is continuous in the space $L_2[0,1]\times L_2[0,1]$. In addition, the constant $\gamma$ defines the kind of nonlinear integral equations.

2. Existence of a Unique Solution for the Integral Equation with Symmetric and Nonsymmetrical Kernel

The existence of a unique solution of problem (1) is discussed and proved in this section using the Banach fixed point theorem. For this, we write Equation (1) in the form of an integral operator:

$$\overline{V}\psi (x,y)=\frac{1}{\gamma }f(x,y)+V\psi (x,y); \gamma \ne 0,$$

(2)

where

$$\begin{array}{cc}\hfill V\psi (x,y)& \hfill =\frac{{\lambda }_1}{\gamma }{\int }_0^1{\int }_0^1\Phi (x,\tau ;y,v)\mu (\tau ,v,\psi (\tau ,v))\mathrm{d}v\mathrm{d}\tau \hfill \\ \hfill & \hfill +\frac{{\lambda }_2}{\gamma }{\int }_0^x{\int }_0^{1}G(x,\tau ;y,v)\nu (\tau ,v,\psi (\tau ,v))\mathrm{d}v\mathrm{d}\tau .\hfill \end{array}$$

(3)

We assume the following conditions:

(i)

The kernels $\Phi (x,\tau ;y,v)$ and $G(x,\tau ;y,v)$ satisfy the following conditions: $\parallel \Phi (x,\tau ;y,v)\parallel \le A_1, \parallel G(x,\tau ;y,v))\parallel \le A_2,$ where $A_1$ and $A_2$ are two constants, assume $A=max\{A_1,A_2\}$.

(ii)

$\parallel f(x,y)\parallel ={\left[{\int }_0^1{\int }_0^1{\left|f(x,y)\right|}^2\mathrm{d}x\mathrm{d}y\right]}^{\frac{1}{2}}=D, D$ is a constant.

(iii)

The function $\mu (x,y,\psi (x,y\left)\right)$ satisfies the following conditions:

$$\begin{array}{cc}\hfill & \hfill \parallel \mu (x,y,\psi (x,y))\parallel ={\left[{\int }_0^1{\int }_0^1{\left|\mu (x,y,\psi (x,y))\right|}^2\mathrm{d}x\mathrm{d}y\right]}^{\frac{1}{2}}\le M_1\parallel \psi (x,y)\parallel \hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill & \hfill \parallel \mu (x,y,{\psi }_1(x,y))-\mu (x,y,{\psi }_2(x,y))\parallel \le M_2\parallel {\psi }_1(x,y)-{\psi }_2(x,y)\parallel ,\hfill \end{array}$$

(5)

where $M_1$ and $M_2$ are constants.

(iv)

The function $\nu (x,y,\psi (x,y\left)\right)$ is bounded and satisfies the following:

$$\begin{array}{cc}\hfill & \hfill \parallel \nu (x,y,\psi (x,y))\parallel ={\left[{\int }_0^1{\int }_0^1{\left|\nu (x,y,\psi (x,y))\right|}^2\mathrm{d}x\mathrm{d}y\right]}^{\frac{1}{2}}\le N_1\parallel \psi (x,y)\parallel \hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & \hfill \parallel \nu (x,y,{\psi }_1(x,y))-\nu (x,y,{\psi }_2(x,y))\parallel \le N_2\parallel {\psi }_1(x,y)-{\psi }_2(x,y)\parallel ,\hfill \end{array}$$

(7)

where, $N_1$ and $N_2$ are constants.

Theorem 1.

Assume that the conditions (i)–(iv) are satisfied. Then, Equation (1) has a unique solution $\psi (x,y)$ in the space $L_2[0,1]\times L_2[0,1]$, if the condition

$$\eta =A\left|\frac{\lambda }{\gamma }\right|[M+N]<1; \left(\lambda =max\{{\lambda }_1,{\lambda }_2\}, M=max\{M_1,M_2\}, N=max\{N_1,N_2\}\right)$$

(8)

is true.

To prove the theorem, the following two lemmas must be proven:

Lemma 1.

Under the conditions (i), (ii), (iii-4) and (iv-6), the operator $\overline{V}\psi (x,y)$, defined by Equation (2), maps the space $L_2[0,1]\times L_2[0,1]$ into itself.

Proof. 

In light of Formulas (2) and (3), we obtain

$$\begin{array}{cc}\hfill \parallel \overline{V}\psi (x,y)\parallel & \hfill \le \frac{1}{\left|\gamma \right|}\parallel f(x,y)\parallel +\left|\frac{{\lambda }_1}{\gamma }\right|∥{\int }_0^1{\int }_0^1|\Phi (x,\tau ;y,v)|\mu (\tau ,v,\psi (\tau ,v))\mathrm{d}v\mathrm{d}\tau ∥\hfill \\ \hfill & \hfill +\left|\frac{{\lambda }_2}{\gamma }\right|∥{\int }_0^x{\int }_0^1\left|G(x,\tau ;y,v)\right|\nu (\tau ,v,\psi (\tau ,v))\mathrm{d}v\mathrm{d}\tau ∥.\hfill \end{array}$$

From conditions (i) and (ii), we obtain

$$\begin{array}{cc}\hfill \parallel \overline{V}\psi (x,y)\parallel & \hfill \le \frac{D}{\left|\gamma \right|}+A\left|\frac{{\lambda }_1}{\gamma }\right|∥{\int }_0^1{\int }_0^1|\mu (\tau ,v,\psi (\tau ,v))\mathrm{d}v\mathrm{d}\tau ∥\hfill \\ \hfill & \hfill +A\left|\frac{{\lambda }_2}{\gamma }\right|∥{\int }_0^x{\int }_0^1\nu (\tau ,v,\psi (\tau ,v))\mathrm{d}v\mathrm{d}\tau ∥.\hfill \end{array}$$

Given conditions $(iii$-$4)$ and $(iv$-$6)$, the above inequality takes on the following form:

$$\begin{array}{cc}\hfill \parallel \overline{V}\psi (x,y)\parallel & \hfill \le \frac{D}{\left|\gamma \right|}+AM\parallel \psi (x,y)\parallel \left|\frac{{\lambda }_1}{\gamma }\right|∥{\int }_0^1{\int }_0^1\mathrm{d}v\mathrm{d}\tau ∥\hfill \\ \hfill & \hfill +AN\parallel \psi (x,y)\parallel \left|\frac{{\lambda }_2}{\gamma }\right|∥{\int }_0^x{\int }_0^1\mathrm{d}v\mathrm{d}\tau ∥,\hfill \end{array}$$

where ${max}_{0\le x\le 1}\left|x\right|=1$, so that the last inequality becomes

$$\begin{array}{c}\hfill \parallel \overline{V}\psi (x,y)\parallel \le \frac{D}{\left|\gamma \right|}+A\left|\frac{\lambda }{\gamma }\right|[M+N]\parallel \psi (x,y)\parallel ,\end{array}$$

since

$$\parallel \overline{V}\psi (x,y)\parallel \le \frac{D}{\left|\gamma \right|}+\eta \parallel \psi (x,y)\parallel ; \eta =A\left|\frac{\lambda }{\gamma }\right|[M+N]<1.$$

(9)

According to this inequality, the operator $\overline{V}$ maps the ball $B_r\subset L_2[0,1]\times L_2[0,1]$ into itself, where

$$r=\frac{D}{\left|\gamma \right|(1-\eta )}.$$

□

Lemma 2.

If the conditions (i), (iii-5) and (iv-7) are verified, then the operator $\overline{V}\psi (x,y)$ defined by Equation (2) is continuous in the space $L_2[0,1]\times L_2[0,1]$.

Proof. 

Suppose two functions, ${\Psi }_1(x,y)$ and ${\Psi }_2(x,y),$ satisfy Equation (2). Then,

$$\begin{array}{cc}\hfill \overline{V}{\psi }_1(x,y)-\overline{V}{\psi }_2(x,y)& \hfill =\frac{{\lambda }_1}{\gamma }{\int }_0^1{\int }_0^1\Phi (x,\tau ;y,v)[\mu (\tau ,v,{\psi }_1(\tau ,v))-\mu (\tau ,v,{\psi }_2(\tau ,v))]\mathrm{d}v\mathrm{d}\tau \hfill \\ \hfill & \hfill +\frac{{\lambda }_2}{\gamma }{\int }_0^x{\int }_0^{1}G(x,\tau ;y,v)[\nu (\tau ,v,{\psi }_1(\tau ,v))-\nu (\tau ,v,{\psi }_2(\tau ,v))]\mathrm{d}v\mathrm{d}\tau ,\hfill \end{array}$$

applying the properties of the norm, we obtain

$$\begin{array}{cc}\hfill \parallel \overline{V}{\psi }_1(x,y)-\overline{V}{\psi }_2(x,y)\parallel & \hfill \le \left|\frac{{\lambda }_1}{\gamma }\right|∥{\int }_0^1{\int }_0^1|\Phi (x,\tau ;y,v)||\mu (\tau ,v,{\psi }_1(\tau ,v))-\mu (\tau ,v,{\psi }_2(\tau ,v))|\mathrm{d}v\mathrm{d}\tau ∥\hfill \\ \hfill & \hfill +\left|\frac{{\lambda }_2}{\gamma }\right|∥{\int }_0^x{\int }_0^1\left|G(x,\tau ;y,v)\right||\nu (\tau ,v,{\psi }_1(\tau ,v))-\nu (\tau ,v,{\psi }_2(\tau ,v))|\mathrm{d}v\mathrm{d}\tau ∥.\hfill \end{array}$$

In view of the conditions $\left(i\right)$, $(iii$-$5)$ and $(iv$-$7)$, the above inequality becomes

$$\begin{array}{c}\hfill \parallel \overline{V}{\psi }_1(x,y)-\overline{V}{\psi }_2(x,y)\parallel \le A\left|\frac{\lambda }{\gamma }\right|[M+N]\parallel {\psi }_1(x,y)-{\psi }_2(x,y)\parallel ,\end{array}$$

since

$$\parallel \overline{V}{\psi }_1(x,y)-\overline{V}{\psi }_2(x,y)\parallel \le \eta \parallel {\psi }_1(x,y)-{\psi }_2(x,y)\parallel .$$

This inequality shows that $\overline{V}$ is a continuous operator in $L_2[0,1]\times L_2[0,1]$. Moreover, $\overline{V}$ is a contraction operator under the condition $\eta <1$.

The previous two Lemmas 1 and 2 show that the operator $\overline{V}$ defined by (2) is a contraction operator in the space $L_2[0,1]\times L_2[0,1]$. Hence, from the Banach fixed point theorem, $\overline{V}$ has a unique fixed point which is, of course, the unique solution of Equation (1). □

3. Method of Solution for the Main Problem

This section applies the collocation method, two-dimensional hybrid functions and the Gauss quadrature formula to transform the integral Equation (1) into nonlinear systems of equations. The following results are obtained by expanding the function $\Psi (x,y)$ in Equation (1) with respect to two-dimensional hybrid functions:

$$\Psi (x,y)=\sum _{m_1=1}^{\infty }\sum _{n_1=0}^{\infty }\sum _{m_2=1}^{\infty }\sum _{n_2=0}^{\infty }{c}_{m_1{n}_1{m}_2{n}_2}{h}_{m_1{n}_1{m}_2{n}_2}(x,y),$$

(10)

where the finite series in Equation (10) can be written as

$${\Psi }_{S,K}(x,y)=\sum _{m_1=1}^S\sum _{n_1=0}^{K-1}\sum _{m_2=1}^S\sum _{n_2=0}^{K-1}{c}_{m_1{n}_1{m}_2{n}_2}{h}_{m_1{n}_1{m}_2{n}_2}(x,y).$$

(11)

where $c_{m_1{n}_1{m}_2{n}_2}, m_1,m_2=1,2,\dots ,S, n_1,n_2=0,1,2,\dots ,K-1,$ and $S,K$ are the unknown hybrid coefficients.

Substituting Equation (11) into Equation (1) yields

$$\begin{array}{ll}\gamma {\psi }_{S,K}(x,y)& =f(x,y)+{\lambda }_1{\int }_0^1{\int }_0^1\Phi (x,\tau ;y,v)\mu (\tau ,v,{\psi }_{S,K}(\tau ,v))\mathrm{d}v\mathrm{d}\tau \\ & +{\lambda }_2{\int }_0^x{\int }_0^{1}G(x,\tau ;y,v)\nu (\tau ,v,{\psi }_{S,K}(\tau ,v))\mathrm{d}v\mathrm{d}\tau .\end{array}$$

(12)

Now, we discretize Equation (12) at the set of collocation nodes $(x_m,y_n)$ for $m,n=1,2,\dots ,SK$ as follows:

$$\begin{array}{ll}\gamma {\psi }_{S,K}(x_m,y_n)& =f(x_m,y_n)+{\lambda }_1{\int }_0^1{\int }_0^1\Phi (x_m,\tau ;y_n,v)\mu (\tau ,v,{\psi }_{S,K}(\tau ,v))\mathrm{d}v\mathrm{d}\tau \\ & +{\lambda }_2{\int }_0^{x_m}{\int }_0^{1}G(x_m,\tau ;y_n,v)\nu (\tau ,v,{\psi }_{S,K}(\tau ,v))\mathrm{d}v\mathrm{d}\tau ,\end{array}$$

(13)

where

$$x_m=\frac{1}{2}\left(cos\left(\frac{(2m-1)\pi }{2SK}\right)+1\right), m=1,2,\dots ,SK,$$

and

$$y_n=\frac{1}{2}\left(cos\left(\frac{(2n-1)\pi }{2SK}\right)+1\right), n=1,2,\dots ,SK,$$

the integral operators in Equation (13) are approximated by using the Gauss–Legendre quadrature formula. For this, we use the following transformations to convert the integrals over $[0,1]$ into the integral over $[-1,1]$,

$$\begin{array}{cc}\hfill \xi & \hfill =2\tau -1; \tau \in [0,1],\hfill \\ \hfill \varrho & \hfill =2v-1; v\in [0,1].\hfill \end{array}$$

The integral over $[0,x_m]$ must also be changed into the integral over $[-1,1]$, having the following form

$$\begin{array}{cc}\hfill \overline{\xi }& \hfill =\frac{2}{x_m}\tau -1; \tau \in [0,x_m].\hfill \end{array}$$

Then, Equation (13) is converted to

$$\begin{array}{ll} & \gamma {\psi }_{S,K}(x_m,y_n)=f(x_m,y_n)\\ +& \frac{{\lambda }_1}{4}{\int }_{-1}^1{\int }_{-1}^1\Phi (x_m,\frac{1}{2}(\xi +1);y_n,\frac{1}{2}(\varrho +1))\mu (\frac{1}{2}(\xi +1),\frac{1}{2}(\varrho +1),{\psi }_{S,K}(\frac{1}{2}(\xi +1),\frac{1}{2}(\varrho +1)))\mathrm{d}\varrho \mathrm{d}\xi \\ +& \frac{{\lambda }_2{x}_m}{4}{\int }_{-1}^1{\int }_{-1}^{1}G(x_m,\frac{x_m}{2}(\overline{\xi }+1);y_n,\frac{1}{2}(\varrho +1))\\ \times & \nu (\frac{x_m}{2}(\overline{\xi }+1),\frac{1}{2}(\varrho +1),{\psi }_{S,K}(\frac{x_m}{2}(\overline{\xi }+1),\frac{1}{2}(\varrho +1)))\mathrm{d}\varrho \mathrm{d}\overline{\xi }.\end{array}$$

The above equation can be expressed as follows using Gauss–Legendre quadrature:

$$\begin{array}{ll}\gamma & {\psi }_{S,K}(x_m,y_n)=f(x_m,y_n)\\ & +\frac{{\lambda }_1}{4}\sum _{j=1}^{{\ell }_1}\sum _{i=1}^{{\ell }_2}{w}_j{w}_i\Phi (x_m,\frac{1}{2}({\xi }_i+1);y_n,\frac{1}{2}({\varrho }_j+1))\\ & \times \mu (\frac{1}{2}({\xi }_i+1),\frac{1}{2}({\varrho }_j+1),{\psi }_{S,K}(\frac{1}{2}({\xi }_i+1),\frac{1}{2}({\varrho }_j+1)))\\ & +\frac{{\lambda }_2{x}_m}{4}\sum _{j=1}^{{\ell }_1}\sum _{i=1}^{{\overline{\ell }}_2}{w}_j{\overline{w}}_{i}G(x_m,\frac{x_m}{2}({\overline{\xi }}_i+1);y_n,\frac{1}{2}({\varrho }_j+1))\\ & \times \nu (\frac{x_m}{2}({\overline{\xi }}_i+1),\frac{1}{2}({\varrho }_j+1),{\psi }_{S,K}(\frac{x_m}{2}({\overline{\xi }}_i+1),\frac{1}{2}({\varrho }_j+1))),\\ & m=1,2,\dots ,SK, n=1,2,\dots ,SK,\end{array}$$

(14)

and $w_j,w_i$ and ${\overline{w}}_i$ are the corresponding weights.

This technique can be used to transform the two-dimensional nonlinear integral problem (1) into a solvable nonlinear system of algebraic equations.

4. Convergence Analysis

The aim of this section is to describe the uniform convergence of the hybrid functions expansion and to determine the maximum absolute truncation error of the function $\Psi$ based on hybrid functions.

Theorem 2.

If $\Psi \in C^4[0,1]$, then the function $\Psi (x,y)$ converges uniformly to the infinite sum of the hybrid functions of $\Psi (x,y)$ described by (10).

Proof. 

The hybrid coefficients are defined as

$$\begin{array}{ll}{c}_{m_1{n}_1{m}_2{n}_2}& =\frac{{\int }_0^1{\int }_0^1\Psi (x,y)h_{m_1{n}_1{m}_2{n}_2}(x,y)\mathrm{d}x\mathrm{d}y}{{\int }_0^1{\int }_0^1{h}_{m_1{n}_1{m}_2{n}_2}^2(x,y)\mathrm{d}x\mathrm{d}y}\\ & =\frac{{\int }_{\frac{m_2-1}{S}}^{\frac{m_2}{S}}{\int }_{\frac{m_1-1}{S}}^{\frac{m_1}{S}}\Psi (x,y)L_{n_1}(2Sx-2m_1+1)L_{n_2}(2Sy-2m_2+1)\mathrm{d}x\mathrm{d}y}{{\int }_{\frac{m_1-1}{S}}^{\frac{m_1}{S}}{L}_{n_1}^2(2Sx-2m_1+1)\mathrm{d}x{\int }_{\frac{m_2-1}{S}}^{\frac{m_2}{S}}{L}_{n_2}^2(2Sy-2m_2+1)\mathrm{d}y}.\end{array}$$

Suppose that $2m_1-1={\hat{m}}_1$ and $2Sx-{\hat{m}}_1=\Im$. Then,

$$\begin{array}{ll}{c}_{m_1{n}_1{m}_2{n}_2}& =\frac{{\int }_{\frac{m_2-1}{S}}^{\frac{m_2}{S}}\left({\int }_{-1}^1\Psi (\frac{{\hat{m}}_1+\Im }{2S},y)L_{n_1}(\Im )\mathrm{d}\Im \right)L_{n_2}(2Sy-2m_2+1)\mathrm{d}y}{{\int }_{-1}^1{L}_{n_1}^2(\Im )\mathrm{d}\Im {\int }_{\frac{m_2-1}{S}}^{\frac{m_2}{S}}{L}_{n_2}^2(2Sy-2m_2+1)\mathrm{d}y}\\ & =\frac{(2n_1+1)}{2}\frac{{\int }_{\frac{m_2-1}{S}}^{\frac{m_2}{S}}\left({\int }_{-1}^1\Psi (\frac{{\hat{m}}_1+\Im }{2S},y)L_{n_1}(\Im )\mathrm{d}\Im \right)L_{n_2}(2Sy-2m_2+1)\mathrm{d}y}{{\int }_{\frac{m_2-1}{S}}^{\frac{m_2}{S}}{L}_{n_2}^2(2Sy-2m_2+1)\mathrm{d}y}.\end{array}$$

From the technique of integration by parts with regard to ℑ and $(2n+1)L_n(\Im )=L_{n+1}^{\prime }(\Im )-L_{n-1}^{\prime }(\Im )$, we obtain

$$c_{m_1{n}_1{m}_2{n}_2}=-\frac{1}{2}\frac{{\int }_{\frac{m_2-1}{S}}^{\frac{m_2}{S}}\left({\int }_{-1}^1\frac{\partial }{\partial \Im }\Psi (\frac{{\hat{m}}_1+\Im }{2S},y)(L_{n_1+1}(\Im )-L_{n_1-1}(\Im ))\mathrm{d}\Im \right)L_{n_2}(2Sy-2m_2+1)\mathrm{d}y}{{\int }_{\frac{m_2-1}{S}}^{\frac{m_2}{S}}{L}_{n_2}^2(2Sy-2m_2+1)\mathrm{d}y}.$$

Once again, integration by parts of the above relation results

$$\begin{array}{l}{c}_{m_1{n}_1{m}_2{n}_2}=\\ \frac{1}{2}\frac{{\int }_{\frac{m_2-1}{S}}^{\frac{m_2}{S}}\left({\int }_{-1}^1\frac{{\partial }^2}{\partial {\Im }^2}\Psi (\frac{{\hat{m}}_1+\Im }{2S},y)[\frac{-L_{n_1}(\Im )}{2n_1+3}-\frac{L_{n_1}(\Im )-L_{n_1-2}(\Im )}{2n_1-1}]\mathrm{d}\Im \right)L_{n_2}(2Sy-2m_2+1)\mathrm{d}y}{{\int }_{\frac{m_2-1}{S}}^{\frac{m_2}{S}}{L}_{n_2}^2(2Sy-2m_2+1)\mathrm{d}y}.\end{array}$$

Now, we have

$$c_{m_1{n}_1{m}_2{n}_2}=\frac{1}{2(2n_1+3)(2n_1-1)}\frac{{\int }_{\frac{m_2-1}{S}}^{\frac{m_2}{S}}\left({\int }_{-1}^1\frac{{\partial }^2}{\partial {\Im }^2}\Psi (\frac{{\hat{m}}_1+\Im }{2S},y){\aleph }_{n_1}(\Im )\mathrm{d}\Im \right)L_{n_2}(2Sy-2m_2+1)\mathrm{d}y}{{\int }_{\frac{m_2-1}{S}}^{\frac{m_2}{S}}{L}_{n_2}^2(2Sy-2m_2+1)\mathrm{d}y},$$

where

$${\aleph }_{n_1}(\Im )=(2n_1-1)L_{n_1+2}(\Im )-(4n_1+2)L_{n_1}(\Im )+(2n_1+3)L_{n_1-2}(\Im ).$$

Similarly, changing the variable for y as $2m_2-1={\hat{m}}_2$, where $2Sy-{\hat{m}}_2=\wp$, and integrating by parts with respect to ℘, we obtain

$$\begin{array}{ll}{c}_{m_1{n}_1{m}_2{n}_2}& \hfill =\frac{1}{4(2n_1+3)(2n_1-1)(2n_2+3)(2n_2-1)}\hfill \\ & {\int }_{-1}^1{\int }_{-1}^1\frac{{\partial }^4}{\partial {\wp }^2\partial {\Im }^2}\Psi (\frac{{\hat{m}}_1+\Im }{2S},\frac{{\hat{m}}_2+\wp }{2S}){\aleph }_{n_1}(\Im ){\aleph }_{n_2}(\wp )\mathrm{d}\Im \mathrm{d}\wp ,\end{array}$$

where

$${\aleph }_{n_2}(\wp )=(2n_2-1)L_{n_2+2}(\wp )-(4n_2+2)L_{n_2}(\wp )+(2n_2+3)L_{n_2-2}(\wp ).$$

Using the chain derivatives and $\sigma ={max}_{(x,y)\in [0,1]}\left|\frac{{\partial }^4\Psi (x,y)}{\partial x^2\partial y^2}\right|$, it follows that

$$\begin{array}{ll}{c}_{m_1{n}_1{m}_2{n}_2}& \le \frac{\sigma }{64S^4(2n_1+3)(2n_1-1)(2n_2+3)(2n_2-1)}{\int }_{-1}^1{\int }_{-1}^1|{\aleph }_{n_1}(\Im )\left|\right|{\aleph }_{n_2}(\wp )|\mathrm{d}\Im \mathrm{d}\wp ,\\ & \le \frac{\sigma }{64m_1^2{m}_2^2(2n_1+3)(2n_1-1)(2n_2+3)(2n_2-1)}{\int }_{-1}^1{\int }_{-1}^1|{\aleph }_{n_1}(\Im )\left|\right|{\aleph }_{n_2}(\wp )|\mathrm{d}\Im \mathrm{d}\wp .\end{array}$$

(15)

However,

$$\begin{array}{cc}\hfill {\left({\int }_{-1}^1\left|{\aleph }_{n_1}(\Im )\right|\mathrm{d}\Im \right)}^2=& \hfill \left({\int }_{-1}^1|(2n_1-1)L_{n_1+2}(\Im )-(4n_1+2)L_{n_1}(\Im )+(2n_1+3)L_{n_1-2}(\Im )|\mathrm{d}\Im \right),\hfill \\ \hfill \le & \hfill 2{\int }_{-1}^1|{(2n_1-1)}^2{L}_{n_1+2}^2(\Im )+{(4n_1+2)}^2{L}_{n_1}^2(\Im )+{(2n_1+3)}^2{L}_{n_1-2}^2(\Im )|\mathrm{d}\Im ,\hfill \end{array}$$

from the Legendre orthogonal polynomial property, we obtain

$${\left({\int }_{-1}^1\left|{\aleph }_{n_1}(\Im )\right|\mathrm{d}\Im \right)}^2\le \frac{24{(2n_1+3)}^2}{2n_1-3},$$

(16)

thus

$${\int }_{-1}^1\left|{\aleph }_{n_1}(\Im )\right|\mathrm{d}\Im \le \frac{2\sqrt{6}(2n_1+3)}{\sqrt{2n_1-3}}$$

(17)

and

$${\int }_{-1}^1\left|{\aleph }_{n_2}(\wp )\right|\mathrm{d}\wp \le \frac{2\sqrt{6}(2n_2+3)}{\sqrt{2n_2-3}}.$$

(18)

By substituting (17) and (18) into (15), we obtain

$$\begin{array}{ll}{c}_{m_1{n}_1{m}_2{n}_2}& \le \frac{24\sigma }{64m_1^2{m}_2^2(2n_1-1)(2n_2-1)\sqrt{(2n_1-3)}\sqrt{(2n_2-3)}},\\ & \le \frac{3\sigma }{8m_1^2{m}_2^2{(2n_1-3)}^{\frac{3}{2}}{(2n_2-3)}^{\frac{3}{2}}}.\end{array}$$

(19)

Therefore, the following series is absolutely convergent:

$$\begin{array}{ll}|\Psi (x,y\left)\right|& =\left|\sum _{m_1=1}^{\infty }\sum _{n_1=0}^{\infty }\sum _{m_2=1}^{\infty }\sum _{n_2=0}^{\infty }{c}_{m_1{n}_1{m}_2{n}_2}{h}_{m_1{n}_1{m}_2{n}_2}(x,y)\right|\\ & \le \sum _{m_1=1}^{\infty }\sum _{n_1=0}^{\infty }\sum _{m_2=1}^{\infty }\sum _{n_2=0}^{\infty }\left|c_{m_1{n}_1{m}_2{n}_2}\right|\\ & <\infty ,\end{array}$$

and the series (10) converges to the function $\Psi (x,y)$ uniformly. □

Theorem 3.

The maximum absolute truncation error of the series solution (10) to nonlinear integral Equation (1) is

$$\parallel \Psi (x,y)-{\Psi }_{S,K}(x,y)\parallel \le \frac{3\sigma }{8S}{\left(\sum _{m_1=S+1}^{\infty }\frac{1}{m_1^4}\sum _{n_1=K}^{\infty }\frac{1}{{(2n_1-3)}^4}\sum _{m_2=S+1}^{\infty }\frac{1}{m_2^4}\sum _{n_2=K}^{\infty }\frac{1}{{(2n_2-3)}^4}\right)}^{\frac{1}{2}}.$$

Proof. 

$$\begin{array}{l}\parallel \Psi (x,y)-{\Psi }_{S,K}(x,y)\parallel \\ \hfill \le {\left(\sum _{m_1=S+1}^{\infty }\sum _{n_1=K}^{\infty }\sum _{m_2=S+1}^{\infty }\sum _{n_2=K}^{\infty }{c}_{m_1{n}_1{m}_2{n}_2}^2{\int }_0^1{\int }_0^1{h}_{m_1{n}_1{m}_2{n}_2}^2(x,y)\mathrm{d}x\mathrm{d}y\right)}^{\frac{1}{2}}.\hfill \end{array}$$

Using the orthogonality property of hybrid functions and taking relation (19) into consideration, we obtain

$$\begin{array}{ll}\parallel \Psi (x,y)-{\Psi }_{S,K}(x,y)\parallel & \le {\left(\sum _{m_1=S+1}^{\infty }\sum _{n_1=K}^{\infty }\sum _{m_2=S+1}^{\infty }\sum _{n_2=K}^{\infty }{c}_{m_1{n}_1{m}_2{n}_2}^2\frac{1}{S^2(2n_1+1)(2n_2+1)}\right)}^{\frac{1}{2}}\\ & \le \frac{3\sigma }{8S}{\left(\sum _{m_1=S+1}^{\infty }\frac{1}{m_1^4}\sum _{n_1=K}^{\infty }\frac{1}{{(2n_1-3)}^4}\sum _{m_2=S+1}^{\infty }\frac{1}{m_2^4}\sum _{n_2=K}^{\infty }\frac{1}{{(2n_2-3)}^4}\right)}^{\frac{1}{2}}.\end{array}$$

□

5. Application and Numerical Results

In order to show the accuracy and efficiency of the proposed method, some numerical examples are given in this section. We introduce the following notation to study the absolute values of this method’s errors:

$$R_{S,K}=|\Psi (x,y)-{\Psi }_{S,K}(x,y)|,$$

where $\Psi (x,y)$ and ${\Psi }_{S,K}(x,y)$ are the exact solution and the approximate solution of the integral equations, respectively.

Example 1.

Consider the following two-dimensional nonlinear integral equation with a symmetric and nonsymmetrical kernel:

$$16\psi (x,y)=f(x,y)+{\int }_0^1{\int }_0^1(x\tau +yv){\psi }^2(\tau ,v)\mathrm{d}v\mathrm{d}\tau +{\int }_0^x{\int }_0^1(x^2{\tau }^2+yv){\psi }^2(\tau ,v)\mathrm{d}v\mathrm{d}\tau ,$$

(20)

where

$$f(x,y)=\frac{-7}{24}-\frac{28yx}{45}+16(x^2+y^2)-\frac{x^2{y}^2}{360}(30y^4+72y^{5}x+45y^2{x}^2+80y^3{x}^3+30x^4+72y{x}^5).$$

Exact solution is $\psi (x,y)=x^2+y^2,$ using the proposed numerical technique, where $S=2$ and $K=2,4,6,8$ in the interval $[0,1)$.

In

Table 1. Absolute error of solution of Equation (20) by using the present method with $S=2$ and $K=2,4,6,8$.

$({\mathit{x}}_{\mathit{i}},{\mathit{y}}_{\mathit{i}})$	$\mathit{S}=2, \mathit{K}=2$	$\mathit{S}=2, \mathit{K}=4$	$\mathit{S}=2, \mathit{K}=6$	$\mathit{S}=2, \mathit{K}=8$
(0, 0)	5.62845$\times {10}^{-9}$	3.25447$\times {10}^{-10}$	2.36512$\times {10}^{-13}$	1.32654$\times {10}^{-16}$
(0.1, 0.1)	2.51405$\times {10}^{-7}$	2.36524$\times {10}^{-8}$	1.36524$\times {10}^{-10}$	6.32514$\times {10}^{-13}$
(0.2, 0.2)	5.62103$\times {10}^{-6}$	2.36985$\times {10}^{-7}$	5.36214$\times {10}^{-9}$	8.22551$\times {10}^{-12}$
(0.3, 0.3)	2.02154$\times {10}^{-4}$	3.58412$\times {10}^{-5}$	8.32541$\times {10}^{-8}$	6.32165$\times {10}^{-10}$
(0.4, 0.4)	4.58721$\times {10}^{-4}$	3.65413$\times {10}^{-4}$	2.21345$\times {10}^{-7}$	1.32114$\times {10}^{-9}$
(0.5, 0.5)	7.36212$\times {10}^{-4}$	2.23651$\times {10}^{-4}$	3.65221$\times {10}^{-7}$	2.36985$\times {10}^{-8}$
(0.6, 0.6)	1.36521$\times {10}^{-3}$	1.65214$\times {10}^{-4}$	7.32651$\times {10}^{-7}$	2.92541$\times {10}^{-8}$
(0.7, 0.7)	5.26512$\times {10}^{-3}$	1.36524$\times {10}^{-3}$	6.32541$\times {10}^{-6}$	6.32548$\times {10}^{-8}$
(0.8, 0.8)	5.62514$\times {10}^{-2}$	4.36210$\times {10}^{-3}$	8.36251$\times {10}^{-6}$	7.32614$\times {10}^{-8}$
(0.9, 0.9)	5.65214$\times {10}^{-2}$	6.25489$\times {10}^{-3}$	5.32658$\times {10}^{-5}$	1.36524$\times {10}^{-6}$

, we present the absolute error $|\Psi (x,y)-{\Psi }_{S,K}(x,y)|$, using the introduced numerical method with $S=2$ and $K=2,4,6,8$ in the interval $[0,1)$.

Table 2. The maximum error $R_{max}(x,y)$ for different values of $K=2,4,6,8$ and $S=2$ for Equation (20).

	$\mathit{S}=2, \mathit{K}=2$	$\mathit{S}=2, \mathit{K}=4$	$\mathit{S}=2, \mathit{K}=6$	$\mathit{S}=2, \mathit{K}=8$
$R_{max}$	6.2103$\times {10}^{-2}$	6.53210$\times {10}^{-3}$	5.32658$\times {10}^{-5}$	1.36524$\times {10}^{-6}$

shows the maximum absolute errors of the given method.

Moreover, in Figure 1, Figure 2, Figure 3 and Figure 4, we show a comparison between the exact solution and the approximate solution using the presented numerical technique with different values of $K=2,4,6,8$ with $S=2$ in the interval $[0,1)$.

Example 2.

Consider the nonlinear integral equation with a symmetric and nonsymmetrical kernel:

$$\begin{array}{ll}\psi (x,y)=f(x,y)& +0.003{\int }_0^1{\int }_0^1(x{\tau }^2+vcosy){\psi }^3(\tau ,v)\mathrm{d}v\mathrm{d}\tau \\ & +0.003{\int }_0^x{\int }_0^1(x^2\tau +yv)\psi (\tau ,v)\mathrm{d}v\mathrm{d}\tau ,\end{array}$$

(21)

where

$$\begin{array}{l}f(x,y)=\frac{1}{12}\left(27+16cos1-16cos2-18cos3+7cos4-12t^{2}x(2+cos1)sin{\left(\frac{1}{4}\right)}^4\right.\\ \left.-12sin1-36sin2+6sin3+6sin4\right)+xsiny+\frac{1}{2}{y}^3\left(-3x^2+x(2+3x)cosx-2sinx\right).\end{array}$$

The exact solution is $\psi (x,y)=xsiny,$ using the presented numerical technique with $S=2$ and $K=3,5,7,9$ in the interval $[0,1)$.

In

Table 3. Absolute error of solution of Equation (21) by using the present method with $S=2$ and $K=3,5,7,9$.

$({\mathit{x}}_{\mathit{i}},{\mathit{y}}_{\mathit{i}})$	$\mathit{S}=2, \mathit{K}=3$	$\mathit{S}=2, \mathit{K}=5$	$\mathit{S}=2, \mathit{K}=7$	$\mathit{S}=2, \mathit{K}=9$
(0, 0)	3.20514$\times {10}^{-5}$	5.32641$\times {10}^{-6}$	6.32141$\times {10}^{-9}$	2.36541$\times {10}^{-11}$
(0.1, 0.1)	3.25481$\times {10}^{-4}$	9.32541$\times {10}^{-5}$	5.32187$\times {10}^{-7}$	3.65874$\times {10}^{-8}$
(0.2, 0.2)	3.32541$\times {10}^{-3}$	3.21554$\times {10}^{-4}$	2.36414$\times {10}^{-6}$	7.36584$\times {10}^{-8}$
(0.3, 0.3)	4.32641$\times {10}^{-3}$	5.32654$\times {10}^{-4}$	5.32684$\times {10}^{-6}$	3.36241$\times {10}^{-7}$
(0.4, 0.4)	5.36854$\times {10}^{-3}$	6.36524$\times {10}^{-4}$	8.32546$\times {10}^{-6}$	6.32584$\times {10}^{-7}$
(0.5, 0.5)	6.93154$\times {10}^{-3}$	7.1.365$\times {10}^{-4}$	6.32541$\times {10}^{-5}$	8.65241$\times {10}^{-7}$
(0.6, 0.6)	1.32511$\times {10}^{-2}$	3.21547$\times {10}^{-3}$	9.99215$\times {10}^{-5}$	4.32516$\times {10}^{-6}$
(0.7, 0.7)	4.32658$\times {10}^{-2}$	4.36561$\times {10}^{-3}$	1.32154$\times {10}^{-4}$	8.69854$\times {10}^{-6}$
(0.8, 0.8)	5.32666$\times {10}^{-2}$	5.76524$\times {10}^{-3}$	2.34541$\times {10}^{-4}$	4.36215$\times {10}^{-5}$
(0.9, 0.9)	6.32541$\times {10}^{-2}$	7.96525$\times {10}^{-3}$	3.25456$\times {10}^{-4}$	1.05214$\times {10}^{-4}$

, we show the absolute error $|\Psi (x,y)-{\Psi }_{S,K}(x,y)|$, using the introduced numerical method with $S=2$ and $K=3,5,7,9$ in the interval $[0,1)$.

Table 4. The maximum error $R_{max}(x,y)$ for different values of $K=3,5,7,9$ and $S=2$ for Equation (21).

	$\mathit{S}=2, \mathit{K}=3$	$\mathit{S}=2, \mathit{K}=5$	$\mathit{S}=2, \mathit{K}=7$	$\mathit{S}=2, \mathit{K}=9$
$R_{max}$	6.32541$\times {10}^{-2}$	7.96525$\times {10}^{-3}$	3.25456$\times {10}^{-4}$	1.05214$\times {10}^{-4}$

shows the maximum absolute errors of the given method.

Furthermore, in Figure 5, Figure 6, Figure 7 and Figure 8, we present a comparison between the exact solution and the approximate solution using the introduced numerical method with different values of $K=3,5,7,9$ with $S=2$ in the interval $[0,1)$.

Example 3.

Consider the following two-dimensional nonlinear integral equation with a symmetric and nonsymmetrical kernel:

$$7\psi (x,y)=f(x,y)+0.01{\int }_0^1{\int }_0^1{(x-y)}^2\psi (\tau ,v)\mathrm{d}v\mathrm{d}\tau +0.01{\int }_0^x{\int }_0^1(\tau +v){\psi }^2(\tau ,v)\mathrm{d}v\mathrm{d}\tau ,$$

(22)

where

$$f(x,y)=7e^{t}x-0.00859141{(t-x)}^2-0.000416667x^3(2-3x+e^{2t}(-2+4t+3x)).$$

The exact solution is $\psi (x,y)=x{e}^y,$ using the proposed numerical technique, where $S=3$ and $K=2,3,6,7$ in the interval $[0,1)$.

In

Table 5. Absolute error of solution of Equation (22) by using the present method with $S=3$ and $K=2,3,6,7$.

$({\mathit{x}}_{\mathit{i}},{\mathit{y}}_{\mathit{i}})$	$\mathit{S}=3, \mathit{K}=2$	$\mathit{S}=3, \mathit{K}=3$	$\mathit{S}=3, \mathit{K}=6$	$\mathit{S}=3, \mathit{K}=7$
(0, 0)	9.32484$\times {10}^{-10}$	8.32544$\times {10}^{-11}$	5.21432$\times {10}^{-14}$	1.23698$\times {10}^{-16}$
(0.1, 0.1)	7.62514$\times {10}^{-9}$	3.52147$\times {10}^{-9}$	2.81647$\times {10}^{-11}$	5.21478$\times {10}^{-15}$
(0.2, 0.2)	5.03698$\times {10}^{-7}$	3.62149$\times {10}^{-8}$	3.89759$\times {10}^{-10}$	2.98547$\times {10}^{-13}$
(0.3, 0.3)	1.20584$\times {10}^{-5}$	5.20147$\times {10}^{-6}$	5.07896$\times {10}^{-10}$	4.69857$\times {10}^{-11}$
(0.4, 0.4)	2.96521$\times {10}^{-5}$	5.36987$\times {10}^{-5}$	4.36985$\times {10}^{-8}$	3.21458$\times {10}^{-10}$
(0.5, 0.5)	4.36514$\times {10}^{-5}$	5.58741$\times {10}^{-5}$	5.20142$\times {10}^{-8}$	3.69521$\times {10}^{-9}$
(0.6, 0.6)	4.96587$\times {10}^{-5}$	6.36925$\times {10}^{-5}$	1.25847$\times {10}^{-7}$	4.25871$\times {10}^{-9}$
(0.7, 0.7)	1.32548$\times {10}^{-3}$	2.01321$\times {10}^{-4}$	6.87452$\times {10}^{-7}$	3.02587$\times {10}^{-8}$
(0.8, 0.8)	3.02584$\times {10}^{-3}$	3.47585$\times {10}^{-4}$	1.01024$\times {10}^{-6}$	4.08754$\times {10}^{-8}$
(0.9, 0.9)	5.01478$\times {10}^{-3}$	1.75214$\times {10}^{-3}$	4.36954$\times {10}^{-6}$	5.31231$\times {10}^{-8}$

, we present the absolute error $|\Psi (x,y)-{\Psi }_{S,K}(x,y)|$, using the introduced numerical method with $S=3$ and $K=2,3,6,7$ in the interval $[0,1)$.

Table 6. The maximum error $R_{max}(x,y)$ for different values of $K=2,3,6,7$ and $S=3$ for Equation (22).

	$\mathit{S}=3, \mathit{K}=2$	$\mathit{S}=3, \mathit{K}=3$	$\mathit{S}=3, \mathit{K}=6$	$\mathit{S}=3, \mathit{K}=7$
$R_{max}$	5.01478$\times {10}^{-3}$	1.75214$\times {10}^{-3}$	4.36954$\times {10}^{-6}$	5.31231$\times {10}^{-8}$

shows the maximum absolute errors of the given method.

Moreover, in Figure 9, Figure 10, Figure 11 and Figure 12, we show a comparison between the exact solution and the approximate solution using the presented numerical technique with different values of $K=2,3,6,7$ with $S=3$ in the interval $[0,1)$.

6. Conclusions and Remarks

The following can be deduced from the above analysis and discussion:

• Under some conditions, Equation (1) has a unique solution $\Psi (x,y)$ in the space $L_2[0,1]\times L_2[0,1]$.
• After applying the proposed method, a two-dimensional integral equation with a symmetric and nonsymmetrical kernel of the second kind tends to result in an algebraic system of nonlinear equations.
• A nonlinear system of algebraic equations has a solution.
• Three illustrative examples are provided to evaluate and validate the effectiveness and dependability of the proposed method. Tables and Figures are used to show the numerical results. For example, Figure 1, Figure 5 and Figure 9 contained the numerical solution of Examples 1, 2 and 3, respectively, for different values of $x,y$ and K. Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12 formed the absolute errors of each example with different values of x and y.
• In Example 1, absolute errors in four cases of K are presented in Table 1 and Figure 1, Figure 2, Figure 3 and Figure 4. The error increases through x and y. When we take the Max. value error in Figure 1, it is $(6.2103\times {10}^{-2})$ at $S=2, K=2$. Also, the Min. error value in Figure 4 is $(1.36524\times {10}^{-6})$ at $S=2, K=8$ (see Table 2).
• In Example 2, from Table 4 at $S=2, K=3$, the error is as high as possible at point $x=y=0.9$, and its value is $(6.32541\times {10}^{-2})$. Likewise, the error begins to decrease, and when the value of $x=y=0$, its value is $(1.05214\times {10}^{-4})$.
• In Example 3, the error decreases as S and K increase, where the maximum value of the error at $x=y=0.9$ for $S=3, K=2$ is $(5.01478\times {10}^{-3})$, while for $S=3, K=7$, the minimum value of the error is $(5.31231\times {10}^{-8})$.
• In general, the error obtained by the proposed method decreases when the number of $\left(K\right)$ increases.

7. Future Work

We will consider the following two-dimensional nonlinear integral equation with phase-lag term:

$$\begin{array}{cc}\hfill \gamma \psi (x+\delta x,y+\delta y)& =f(x+\delta x,y+\delta y)+{\lambda }_1{\int }_0^1{\int }_0^1\Phi (x+\delta x,\tau ;y+\delta y,v)\mu (\tau ,v,\psi (\tau ,v))\mathrm{d}v\mathrm{d}\tau \hfill \\ \hfill & +{\lambda }_2{\int }_0^{x+\delta x}{\int }_0^{1}G(x+\delta x,\tau ;y+\delta y,v)\nu (\tau ,v,\psi (\tau ,v))\mathrm{d}v\mathrm{d}\tau .\hfill \end{array}$$