# A Novel and Efficient Numerical Algorithm for Solving 2D Fredholm Integral Equations

Authors: H. Bin Jebreen
Journal: Journal of Mathematics
 (2020)

---

## Abstract

A novel and efficient numerical method is developed based on interpolating scaling functions to solve 2D Fredholm integral equations (FIE). Using the operational matrix of integral for interpolating scaling functions, FIE reduces to a set of algebraic equations that one can obtain an approximate solution by solving this system. The convergence analysis is investigated, and some numerical experiments confirm the accuracy and validity of the method. To show the ability of the proposed method, we compare it with others.


---

## Body

## 1. Introduction

In this work, we focus on developing the wavelet Galerkin scheme for 2D Fredholm integral equations (2D-FIE), given by(1)wx,y−∫01∫01kx,y,s,tgws,tdsdt=fx,y,x,y∈D,where the functions kx,y,s,t and fx,y are continuous functions on D≔Ω×Ω and S with S≔x,y,s,t:x,s∈Ω,y,t∈Ω, respectively, as prescribed before. Function g is given and can be linear or nonlinear. The unknown wx,y is sought.Most of the mathematical models for physical phenomena use integral equations. Consequently, the application of integral equations is a momentous subject within applied mathematics. Also, this equation plays a significant role in reformulations of other mathematical problems. Principally, the Fredholm integral equation can be obtained from converting boundary value problems [1]. Equation (1) is an effective tool to model the problems that come from electromagnetic scattering, computer graphics, and aerodynamics [2, 3].Many papers studied the Fredholm integral equations numerically. Among them, we can mention an effective numerical method based on block-pulse functions [3]. Effati et al. [4] utilized a neural network scheme for solving this type of integral equation. A sparse representation based on multiwavelets has been developed and utilized an efficient algorithm for solving the system of FIE [5]. In [2], numerous numerical methods are presented and investigated for solving the numerical solution of FIE. For some related works, we refer the readers to [6–8].There are so many papers on developing and analyzing the numerical methods for solving two-dimensional FIE. In [9], a novel algorithm is developed to solve the problem using the integral mean value theorem. Kazemi et al. [10] proposed the iterative method based on quadrature formula to solve 2D-FEI. A numerical method based on developing the two-dimensional triangular orthogonal functions is applied in [11]. Aziz et al. [12] used Haar wavelet method to solve problem (1). The principal advantage of the method is that, unlike others, it does not apply numerical integration. In [13], an Euler-type method is proposed to solve considering problems. Singh et al. [14] proposed a numerical method based on Legendre scaling functions to solve the multidimensional Fredholm integral equations. In [15], an accurate scheme is developed for solving 2D Fredholm integral equation by using the cosine-trigonometric functions. In [16], an iterative algorithm based on sparse least-squares is used to solve the problem. For some related works, we refer the readers to [17–20].Multiwavelets are found an interesting basis for solving a variety of equations [21–23]. Multiwavelets have some properties of wavelets, such as orthogonality, vanishing moments, and compact support. In contrast to the wavelets and biorthogonal wavelets, they can have high smoothness and high approximate order coupled with short support [24]. Contrary to biorthogonal wavelets, multiwavelets can have the high vanishing moments without enlarging their support [25]. In general, the multiwavelets are a very powerful tool for expressing a variety of operators. At the present work, we apply the interpolating scaling functions constructed in [6, 26].An outline of the remaining part of the paper is as follows. In Section2, we devote the properties of interpolating scaling functions and related projections. In Section 3, we define the wavelet Galerkin method and investigate the convergence analysis. In Section 4, we solve some numerical examples to confirm the accuracy and efficiency of the scheme and show its ability comparing with other methods. Section 5 contains a few concluding remarks.


## 2. Interpolating Scaling Functions

LetJ∈ℤ+∪0. We consider the uniform finite discretizations Ω≔0,1=∪b∈ℬXJ,b, where the subintervals XJ,b≔xb,xb+1 are determined by the point xb≔b/2J, with ℬ≔0,…,2J−1. For k∈ℛ≔0,1,…,r−1, we introduce the subspace VJr as a space of piecewise polynomial bases of degree less than multiplicity parameter r that is spanned by(2)VJr≔Spanϕj,bk≔D2jTbϕk,b∈ℬj,k∈ℛ⊂L2Ω,r≥0,where T and D are the translation and dilation operators, respectively, and ϕkk∈ℛ are the primal interpolating scaling bases introduced by Alpert et al. [26]. Given nodes τkk∈ℛ, which are the roots of Legendre polynomial of degree r, the interpolating scaling bases are defined as(3)ϕkx=2ωkLk2x−1,x∈Ω,0,o.w,where Lktk∈ℛ are the Lagrange interpolating polynomials at the point τkk∈ℛ and ωkk∈ℛ are the Gauss-Legendre quadrature weights [22, 26]. These bases form orthonormal bases on Ω with respect to the L2-inner product. According to the definition of the space VJr, the spaces Vjrj∈ℤ+∪0 have dimension N≔2Jr and obviously are nested(4)V0r⊂V1r⊂⋯⊂VJr⊂⋯⊂L2Ω.Every functionp∈L2Ω can be represented in the form(5)p≈PJrp=∑b∈ℬJ∑k∈ℛpJ,bkϕJ,bk,where PJr is the orthogonal projection that maps L2Ω onto the subspace VJr. To find the coefficients pJ,bk that are determined by p,ϕJ,bk=∫XJ,bfxϕJ,bkxdx, we shall compute the integrals. We apply the r-point Gauss-Legendre quadrature by a suitable choice of the weights ωk and nodes τk for k∈ℛ to avoid these integrals [5, 26] via(6)pJ,bk≈2−J/2ωk2p2−Jτk+12+b,k∈ℛ,b∈ℬJ.Convergence analysis of the projectionPJrp is investigated for the r-times continuously differentiable function p∈ℂrΩ [26].(7)PJrp−p≤2−Jr24rr!supx∈0,1prx.For the full proof of this approximation and further details, we refer the readers to [6]. Thus, we can conclude that PJrp converges to p with rate of convergence O2−Jr.Assume that the vector functionΦJr≔Φr,J,0,…,Φr,J,2J−1T with Φr,J,b≔ϕJ,b0,…,ϕJ,br−1 includes the scaling functions and called multiscaling function. Approximation (5) may be rewritten using the vector P that includes entries Pbr+k+1≔pJ,bk as follows:(8)PJrp=PTΦJr,where P is an N-dimensional vector. The building blocks of these bases construction can be applied to approximate a higher-dimensional function. To this end, one can introduce the two-dimensional subspace VJr,2≔VJr×VJr⊂L2Ω×Ω that is spanned by(9)ϕJ,bkϕJ,b′k′:b,b′∈ℬJ,k,k′∈ℛ.Thus, by this assumption, to derive an approximation of the functionp∈L2Ω×Ω by the projection operator PJr, we get(10)ps,t≈PJrps,t=ΦJrTsPΦJrt,where components of the square matrix P of order N are obtained by(11)Prb+k+1,rb′+k′+1≈2−Jωk2ωk′2p2−Jτ^k+b,2−Jτ^k′+b′,where τ^k=τk+1/2. Consider the 2r-th partial derivatives of p:Ω×Ω⟶ℝ are continuous. Utilizing this assumption, the error of this approximation can be bounded as follows:(12)PJrp−p≤ℳmax21−rJ4rr!2+21−Jr4rr!,where ℳmax is a constant [27, 28].


## 3. Wavelet Galerkin Scheme (WGS)

For the completeness of this paper, we describe the wavelet Galerkin scheme for two-dimensional Fredholm integral equation (1). In the operator form, equation (1) can be written as(13)I−Kw=f,where I is the identity operator and the Fredholm operator(14)Kwx,y≔∫01∫01kx,y,s,tgws,tdsdt,is assumed to be compact on L2D.Suppose that the unique solution of (1) can be approximated as an expansion of the interpolating scaling functions via(15)wx,y≈PJrwx,y=ΦJrTxWΦJry,where W is an N×N matrix whose elements must be found and the projection operator PJr maps L2D on the subspace VJr. Also, we can write the same expansion for function f. These are substituted into (13), and the coefficients W are determined by using the Galerkin method. For later use, introduce the residual as(16)rJr≔I−λKJrwJr−fJr,where wJr≔PJrw, fJr≔PJrf, and KJr≔PJrK. Due to the Galerkin method, it is necessitous that rJr,ΦJriΦJrj=0, for i,j=1:N. To derive the approximate solution, we solve the system of algebraic equations(17)I−λKJrwJr,ΦJriΦJrj=0.This system can be linear or nonlinear and it depends on the functiong. The critical step of this method is how to approximate Kw. For this purpose, assume that(18)kx,y,s,tgws,t≔ΦJrsΓx,yΦJrt,x,y,s,t∈S,where Γx,y is an N×N matrix whose entries are obtained by (11). Substituting (18) into (14), one can write(19)Kwx,y≔∫01∫01kx,y,s,tgws,tdsdt=∫01∫01ΦJrsΓx,yΦJrtdsdt=IϕΦJr1TΓx,yIϕΦJr1=ΦJrxGΦJry,where G is an N×N matrix and Iϕ is the operational matrix of integral for interpolating scaling functions introduced in [22]. Due to this approximation, system (17) may be written as(20)W−G=F,where fx,y≔ΦJrTxFΦJr.


### 3.1. Convergence Analysis

Theorem 1.
Assume thatK is a compact operator, I−K is injective, and the sequence KJr:L2S⟶L2S is pointwise convergent to K and collectively compact. Then, I−KJr−1 exists and is uniformly bounded and the solution of (13) satisfies the error estimate(21)w−wJ≤I−KJr−11−εJ0I−KJr−1w−PJrw+PJrK−PJrKJrw.Proof.
Motivated by the hypothesis of the theorem, one can conclude that(22)ℱJr−ℱℱJr⟶0,asJ⟶∞.
Using (22), for all sufficiently large J, we can write(23)I−K−1KJr−KKJr<1,and as a consequence of this, I−KJr is reversible.
PickJ≥J0 such that(24)εJ0≡supJ≥J0KJr−PJrKJr<1I−KJr−1.
Due to geometric series theorem,I+I−KJr−1KJr−PJrKJr−1 exists and is uniformly bounded, i.e.,(25)I+I−KJr−1KJr−PJrKJr−1≤11−εJ0I−KJr−1.
The investigation of convergence analysis is founded by the approximation ofI−KJr by I−PJrKJr,(26)I−PJrKJr=I−KJr+KJr−PJrKJr=I−KJrI+I−KJr−1KJr−PJrKJr.
Taking norm of (26) and using (25), one can obtain(27)I−PJrKJr−1≤I−KJr−11−εJ0I−KJr−1.
Applying the operatorPJr to both sides of (13) and then rearrange to get(28)I−PJrKJrw=PJrf+w−PJrw+PJrK−PJrKJrw.
Subtracting (28) of the original equation (13), we have(29)I−PJrKJrw−wJ=w−PJrw+PJrK−PJrKJrw.
Taking norm and using (27),(30)w−wJ≤I−KJr−11−εJ0I−KJr−1w−PJrw+PJrK−PJrKJrw.
It is straightforward to show thatPJrK−PJr‖Jrw⟶0 as J⟶∞. Assume that wJ is a sequence of continuous functions so that wJ⟶w as J⟶∞. Since the orthonormal projection PJr satisfies PJr=1, we can obtain(31)w−PJrw≤w−wJ+wJ−PJrwJ+PJrw−wJ≤2w−wJ+wJ−PJrwJ.
Thus,∀ε>0, there exists a number J0 so that, for any J0≤J, one can write w−wJ≤ε/4. This then implies that(32)u−PJrw≤ε2+wJ−PJrwJ.
This implies thatw−PJrw≤ε, for sufficiently large value of J and because ϵ is arbitrary, consequently, PJrw⟶w as J⟶∞.


## 3.1. Convergence Analysis

Theorem 1.
Assume thatK is a compact operator, I−K is injective, and the sequence KJr:L2S⟶L2S is pointwise convergent to K and collectively compact. Then, I−KJr−1 exists and is uniformly bounded and the solution of (13) satisfies the error estimate(21)w−wJ≤I−KJr−11−εJ0I−KJr−1w−PJrw+PJrK−PJrKJrw.Proof.
Motivated by the hypothesis of the theorem, one can conclude that(22)ℱJr−ℱℱJr⟶0,asJ⟶∞.
Using (22), for all sufficiently large J, we can write(23)I−K−1KJr−KKJr<1,and as a consequence of this, I−KJr is reversible.
PickJ≥J0 such that(24)εJ0≡supJ≥J0KJr−PJrKJr<1I−KJr−1.
Due to geometric series theorem,I+I−KJr−1KJr−PJrKJr−1 exists and is uniformly bounded, i.e.,(25)I+I−KJr−1KJr−PJrKJr−1≤11−εJ0I−KJr−1.
The investigation of convergence analysis is founded by the approximation ofI−KJr by I−PJrKJr,(26)I−PJrKJr=I−KJr+KJr−PJrKJr=I−KJrI+I−KJr−1KJr−PJrKJr.
Taking norm of (26) and using (25), one can obtain(27)I−PJrKJr−1≤I−KJr−11−εJ0I−KJr−1.
Applying the operatorPJr to both sides of (13) and then rearrange to get(28)I−PJrKJrw=PJrf+w−PJrw+PJrK−PJrKJrw.
Subtracting (28) of the original equation (13), we have(29)I−PJrKJrw−wJ=w−PJrw+PJrK−PJrKJrw.
Taking norm and using (27),(30)w−wJ≤I−KJr−11−εJ0I−KJr−1w−PJrw+PJrK−PJrKJrw.
It is straightforward to show thatPJrK−PJr‖Jrw⟶0 as J⟶∞. Assume that wJ is a sequence of continuous functions so that wJ⟶w as J⟶∞. Since the orthonormal projection PJr satisfies PJr=1, we can obtain(31)w−PJrw≤w−wJ+wJ−PJrwJ+PJrw−wJ≤2w−wJ+wJ−PJrwJ.
Thus,∀ε>0, there exists a number J0 so that, for any J0≤J, one can write w−wJ≤ε/4. This then implies that(32)u−PJrw≤ε2+wJ−PJrwJ.
This implies thatw−PJrw≤ε, for sufficiently large value of J and because ϵ is arbitrary, consequently, PJrw⟶w as J⟶∞.


## 4. Numerical Examples

To illustrate the accuracy and efficiency of the scheme, some equations of form (1) solve by the proposed method. In this section, the maximum absolute value error,(33)eJ=w−wJ,and L2-error are used to show the accuracy of the proposed method. All numerical computations are carried out simultaneously using Maple and MATLAB software.Example 1.
Let us consider the following 2D linear Fredholm integral equation:(34)kx,y,s,t≔x8+y1+s+t,fx,y≔1+x+y−2−x48+6y.
The exact solution iswx,y=1/1+x+y2 [9, 11].
The proposed method is compared with integral mean value theorem [9] and 2D-TFs method [11] in Table 1. The results in Table 1 show that the proposed method has better accuracy than other methods. Figure 1 shows the surfaces of the approximate solution and the absolute error when r=4 and J=2. In Figure 2, we show the effect of the parameters r and J on L2-error. It is obvious that increasing the parameters r and J reduces the error. The order of convergence is presented numerically for Example 1 in Table 2. As we expected, the order of convergence must approach multiplicity r.Table 1
Comparison of maximum absolute value error with others for Example1.
rWGS[11][9]r=4,J=3r=3,J=4m=n=32m=n=320.2,0.21.12e−066.21e−069.54e−032.58e−060.4,0.43.85e−071.40e−069.03e−035.05e−060.6,0.61.44e−075.83e−071.55e−047.39e−060.8,0.81.51e−083.65e−072.62e−049.63e−06Figure 1
The surfaces of the approximate solution and the absolute error takingr=4 and J=2 for Example 1.
(a)(b)Figure 2
Effects of the multiplicity parameterr and refinement level J for Example 1.
(a)(b)Table 2
Order of convergence for Example1.
r=2r=3r=41.8699394602.7756460623.6716741341.9646187852.9385994553.8998856621.9903239412.9808911773.920293300↓↓↓P234Example 2.
Let us now consider the following 2D linear FIE:(35)kx,y,s,t≔xyes+t,fx,y≔e−x−y.
The exact solutionwx,y=e−x−y−1/2xy is given in [9].
Table3 illustrates a comparison among the numerical results of the proposed method, integral mean value theorem [9], and 2D-TFs method [11]. The results show that the proposed method is better than others. Figure 3 shows the surfaces of the approximate solution and the absolute error taking r=4 and J=2. Table 4 shows the effects of the multiplicity parameter r and refinement level J for Example 2. The order of convergence is presented numerically for Example 2 in Table 5.Table 3
Comparison of maximum absolute value errors with others for Example2.
rWGS[11][9]r=4,J=3m=n=32m=n=320.2,0.29.18e−079.36e−036.49e−060.4,0.42.23e−061.48e−022.60e−050.6,0.64.21e−066.68e−055.84e−050.8,0.86.56e−064.58e−031.04e−05Figure 3
The surfaces of the approximate solution and the absolute error takingr=4 and J=2 for Example 2.
(a)(b)Table 4
Effects of the multiplicity parameterr and refinement level J for Example 2.
rJ=1J=2J=3J=421.46e−024.20e−031.13e−032.95e−0431.16e−031.54e−042.07e−052.69e−0645.53e−053.64e−062.45e−071.59e−08Table 5
Order of convergence for Example2.
r=2r=3r=41.7975071362.9256100003.9256100001.8940665552.8858289803.8930847961.9375359132.9439526903.950893907↓↓↓P234Example 3.
Let us consider the nonlinear 2D-FIE given in [10]:(36)wx,y=cos13xy236−cos1xy212+xy218+siny−∫01∫01xy2s6w3s,tdsdt.
The exact solutionwx,y=siny is given in [10].
In Table6, we compare the proposed method and quadrature formula proposed in [10]. It shows that the proposed method is much more flexible than the method presented in [10] and offers better accuracy using fewer computations. Figure 4 shows the plots of the approximate solution and the absolute error when r=4 and J=2. Figure 5 shows the effects of the multiplicity parameter r and refinement level J for Example 3. The order of convergence is presented numerically for Example 3 in Table 7.Table 6
Comparison of maximum absolute value errors with the proposed method in [10] for Example 3.
0.1,0.10.3,0.30.5,0.50.7,0.70.9,0.9WGM (r=7,J=1)1.37e−109.60e−122.98e−101.85e−093.82e−09[10](n=100)8.18e−102.21e−081.02e−072.81e−075.97e−07Figure 4
The plots of the approximate solution and the absolute error takingr=7 and J=1 for Example 3.
(a)(b)Figure 5
Effects of the multiplicity parameterr and refinement level J for Example 3.
(a)(b)Table 7
Order of convergence for Example3.
r=2r=3r=41.9821519993.0085620143.8537555461.9952775902.9456100003.9456100001.9975359132.9630556103.972803421↓↓↓P234Example 4.
Consider the following nonlinear 2D-FIE given in [12, 18]:(37)wx,y−∫01∫01x1−s2y+1t2+11−e−ws,tdsdt=−log1+xyy2+1+x16y+1.
The exact solution is reported in [12] via(38)wx,y=−log1+xyy2+1.
Table8 illustrates a comparison among the numerical results of the proposed method, Nystrom method [18], and Haar wavelet method [12]. The approximate solution and the absolute error are plotted in Figure 6 when r=4 and J=2.Table 8
Comparison of maximum absolute value errors with others for Example4.
J12345WGS7.6e−054.7e−062.9e−071.8e−081.1e−09Haar wavelet method [12]2.6e−037.3e−041.9e−045.0e−051.3e−05Nystrom method [18]1.4e−023.4e−038.3e−042.1e−045.2e−05Figure 6
The plots of the approximate solution and the absolute error takingr=4 and J=2 for Example 4.
(a)(b)


## 5. Conclusion

Interpolating scaling functions are applied to obtain the approximate solution of 2D Fredholm integral equations. The convergence analysis is investigated and the numerical results are compared with existing methods. It shows that the proposed method is much more flexible than others and the results show that the proposed method has better accuracy than other methods. To increase the accuracy, it is not necessary to increase the degree of polynomials (unlike bases such as Chebyshev and Legendre). By increasing the level of refinementJ, we can increase the accuracy. Also, these bases have the interpolation property and this helps us to avoid the direct integral for finding the coefficients and this reduces the computational cost.



---

*Source: 6662721-2020-12-15.xml*