# Numerical Methods for Solving Fredholm Integral Equations of Second Kind

Authors: S. Saha Ray; P. K. Sahu
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

Integral equation has been one of the essential tools for various areas of applied mathematics. In this paper, we review different numerical methods for solving both linear and nonlinear Fredholm integral equations of second kind. The goal is to categorize the selected methods and assess their accuracy and efficiency. We discuss challenges faced by researchers in this field, and we emphasize the importance of interdisciplinary effort for advancing the study on numerical methods for solving integral equations.


---

## Body

## 1. Introduction

Integral equations occur naturally in many fields of science and engineering [1]. A computational approach to solve integral equation is an essential work in scientific research.Integral equation is encountered in a variety of applications in many fields including continuum mechanics, potential theory, geophysics, electricity and magnetism, kinetic theory of gases, hereditary phenomena in physics and biology, renewal theory, quantum mechanics, radiation, optimization, optimal control systems, communication theory, mathematical economics, population genetics, queuing theory, medicine, mathematical problems of radiative equilibrium, the particle transport problems of astrophysics and reactor theory, acoustics, fluid mechanics, steady state heat conduction, fracture mechanics, and radiative heat transfer problems. Fredholm integral equation is one of the most important integral equations.Integral equations can be viewed as equations which are results of transformation of points in a given vector spaces of integrable functions by the use of certain specific integral operators to points in the same space. If, in particular, one is concerned with function spaces spanned by polynomials for which the kernel of the corresponding transforming integral operator is separable being comprised of polynomial functions only, then several approximate methods of solution of integral equations can be developed.A computational approach to solving integral equation is an essential work in scientific research. Some methods for solving second kind Fredholm integral equation are available in the open literature. TheB-spline wavelet method, the method of moments based on B-spline wavelets by Maleknejad and Sahlan [2], and variational iteration method (VIM) by He [3–5] have been applied to solve second kind Fredholm linear integral equations. The learned researchers Maleknejad et al. proposed some numerical methods for solving linear Fredholm integral equations system of second kind using Rationalized Haar functions method, Block-Pulse functions, and Taylor series expansion method [6–8]. Haar wavelet method with operational matrices of integration [9] has been applied to solve system of linear Fredholm integral equations of second kind. Quadrature method [10], B-spline wavelet method [11], wavelet Galerkin method [12], and also VIM [13] can be applied to solve nonlinear Fredholm integral equation of second kind. Some iterative methods like Homotopy perturbation method (HPM) [14–16] and Adomian decomposition method (ADM) [16–18] have been applied to solve nonlinear Fredholm integral equation of second kind.


## 2. Fredholm Integral Equation

The general form of linear Fredholm integral equation is defined as follows:(1)g(x)y(x)=f(x)+λ∫abK(x,t)y(t)dt,
where a and b are both constants. f(x), g(x), and K(x,t) are known functions while y(x) is unknown function. λ (nonzero parameter) is called eigenvalue of the integral equation. The function K(x,t) is known as kernel of the integral equation.


### 2.1. Fredholm Integral Equation of First Kind

The linear integral equation is of form (by settingg(x)=0 in (1))
(2)f(x)+λ∫abK(x,t)y(t)dt=0.
Equation (2) is known as Fredholm integral equation of first kind.


### 2.2. Fredholm Integral Equation of Second Kind

The linear integral equation is of form (by settingg(x)=1 in (1))
(3)y(x)=f(x)+∫abK(x,t)y(t)dt.
Equation (3) is known as Fredholm integral equation of second kind.


### 2.3. System of Linear Fredholm Integral Equations

The general form of system of linear Fredholm integral equations of second kind is defined as follows:(4)∑j=1ngi,jyj(x)=fi(x)+∑j=1n∫abKi,j(x,t)yj(t)dt,kkkkkkkkkkkkkkkkkkkkkkkkkii=1,2,…,n,
where fi(x) and Ki,j(x,t) are known functions and yj(x) are the unknown functions for i,j=1,2,…,n.


### 2.4. Nonlinear Fredholm-Hammerstein Integral Equation of Second Kind

Nonlinear Fredholm-Hammerstein integral equation of second kind is defined as follows:(5)y(x)=f(x)+∫abK(x,t)F(y(t))dt,
where K(x,t) is the kernel of the integral equation, f(x) and K(x,t) are known functions, and y(x) is the unknown function that is to be determined.


### 2.5. System of Nonlinear Fredholm Integral Equations

System of nonlinear Fredholm integral equations of second kind is defined as follows:(6)∑j=1ngi,jyj(x)=fi(x)+∑j=1n∫abKi,j(x,t)Fi,j(t,yj(t))dt,kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkki=1,2,…,n,
where fi(x) and Ki,j(x,t) are known functions and yj(x) are the unknown functions for i,j=1,2,…,n.


## 2.1. Fredholm Integral Equation of First Kind

The linear integral equation is of form (by settingg(x)=0 in (1))
(2)f(x)+λ∫abK(x,t)y(t)dt=0.
Equation (2) is known as Fredholm integral equation of first kind.


## 2.2. Fredholm Integral Equation of Second Kind

The linear integral equation is of form (by settingg(x)=1 in (1))
(3)y(x)=f(x)+∫abK(x,t)y(t)dt.
Equation (3) is known as Fredholm integral equation of second kind.


## 2.3. System of Linear Fredholm Integral Equations

The general form of system of linear Fredholm integral equations of second kind is defined as follows:(4)∑j=1ngi,jyj(x)=fi(x)+∑j=1n∫abKi,j(x,t)yj(t)dt,kkkkkkkkkkkkkkkkkkkkkkkkkii=1,2,…,n,
where fi(x) and Ki,j(x,t) are known functions and yj(x) are the unknown functions for i,j=1,2,…,n.


## 2.4. Nonlinear Fredholm-Hammerstein Integral Equation of Second Kind

Nonlinear Fredholm-Hammerstein integral equation of second kind is defined as follows:(5)y(x)=f(x)+∫abK(x,t)F(y(t))dt,
where K(x,t) is the kernel of the integral equation, f(x) and K(x,t) are known functions, and y(x) is the unknown function that is to be determined.


## 2.5. System of Nonlinear Fredholm Integral Equations

System of nonlinear Fredholm integral equations of second kind is defined as follows:(6)∑j=1ngi,jyj(x)=fi(x)+∑j=1n∫abKi,j(x,t)Fi,j(t,yj(t))dt,kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkki=1,2,…,n,
where fi(x) and Ki,j(x,t) are known functions and yj(x) are the unknown functions for i,j=1,2,…,n.


## 3. Numerical Methods for Linear Fredholm Integral Equation of Second Kind

Consider the following Fredholm integral equation of second kind defined in (3)
(7)y(x)=f(x)+∫abK(x,t)y(t)dt,a≤x≤b,
where K(x,t) and g(x) are known functions and y(x) is unknown function to be determined.


### 3.1.B-Spline Wavelet Method

#### 3.1.1.B-Spline Scaling and Wavelet Functions on the Interval [0,1]

Semiorthogonal wavelets usingB-spline are specially constructed for the bounded interval and this wavelet can be represented in a closed form. This provides a compact support. Semiorthogonal wavelets form the basis in the space L2(R).Using this basis, an arbitrary function inL2(R) can be expressed as the wavelet series. For the finite interval [0,1], the wavelet series cannot be completely presented by using this basis. This is because supports of some basis are truncated at the left or right end points of the interval. Hence, a special basis has to be introduced into the wavelet expansion on the finite interval. These functions are referred to as the boundary scaling functions and boundary wavelet functions.Letm and n be two positive integers and let
(8)a=x-m+1=⋯=x0<x1<⋯<xn=xn+1=⋯=xn+m-1=b
be an equally spaced knots sequence. The functions
(9)Bm,j,X(x)=x-xjxj+m-1-xjBm-1,j,X(x)+xj+m-xxj+m-xj+1Bm-1,j+1,X(x),j=-m+1,…,n-1,B1,j,X(x)={1,x∈[xj,xj+1),0,otherwise,
are called cardinal B-spline functions of order m≥2 for the knot sequence X={xi}i=-m+1n+m-1 and SuppBm,j,X(x)=[xj,xj+m]∩[a,b].By considering the interval[a,b]=[0,1], at any level j∈Ζ+, the discretization step is 2-j, and this generates n=2j number of segments in [0,1] with knot sequence
(10)X(j)={x-m+1(j)=⋯=x0(j)=0,xk(j)=k2j,k=1,…,n-1,xn(j)=⋯=xn+m-1(j)=1.
Let j0 be the level for which 2j0≥2m-1; for each level, j≥j0, the scaling functions of order m can be defined as follows in [2]:
(11)φm,j,i(x)={Bm,j0,i(2j-j0x)i=-m+1,…,-1,Bm,j0,2j-m-i(1-2j-j0x)i=2j-m+1,…,2j-1,Bm,j0,0(2j-j0x-2-j0i)i=0,…,2j-m.
And the two scale relations for the m-order semiorthogonal compactly supported B-wavelet functions are defined as follows:
(12)ψm,j,i-m=∑k=i2i+2m-2qi,kBm,j,k-m,i=1,…,m-1,ψm,j,i-m=∑k=2i-m2i+2m-2qi,kBm,j,k-m,i=m,…,n-m+1,ψm,j,i-m=∑k=2i-mn+i+m-1qi,kBm,j,k-m,i=n-m+2,…,n,
where qi,k=qk-2i.Hence, there are2(m-1) boundary wavelets and (n-2m+2) inner wavelets in the bounded interval [a,b]. Finally, by considering the level j with j≥j0, the B-wavelet functions in [0,1] can be expressed as follows:
(13)ψm,j,i(x)={ψm,j0,i(2j-j0x)i=-m+1,…,-1,ψm,2j-2m+1-i,i(1-2j-j0x)i=2j-2m+2,…,2j-m,ψm,j0,0(2j-j0x-2-j0i)i=0,…,2j-2m+1.
The scaling functions φm,j,i(x) occupy m segments and the wavelet functions ψm,j,i(x) occupy 2m-1 segments.When the semiorthogonal wavelets are constructed fromB-spline of order m, the lowest octave level j=j0 is determined in [19, 20] by
(14)2j0≥2m-1,
so as to have a minimum of one complete wavelet on the interval [0,1].


#### 3.1.2. Function Approximation

A functionf(x) defined over [0,1] may be approximated by B-spline wavelets as [21, 22]
(15)f(x)=∑k=1-m2j0-1cj0,kφj0,k(x)+∑j=j0∞∑k=1-m2j-mdj,kψj,k(x).
If the infinite series in (15) is truncated at M, then (15) can be written as [2]
(16)f(x)≅∑k=1-m2j0-1cj0,kφj0,k(x)+∑j=j0M∑k=1-m2j-mdj,kψj,k(x),
where φ2,k and ψj,k are scaling and wavelets functions, respectively, and C and Ψ are (2M+1+m-1)×1 vectors given by
(17)C=[cj0,1-m,…,cj0,2j0-1,dj0,1-m,…,dj0,2j0-m,…,dM,1-m,…,dM,2M-m]T,(18)Ψ=[φj0,1-m,…,φj0,2j0-1,ψj0,1-m,…,ψj0,2j0-m,…,ψM,1-m,…,ψM,2M-m]T,
with
(19)cj0,k=∫01f(x)φ~j0,k(x)dx,k=1-m,…,2j0-1,kkkkkkkkkkdj,k=∫01f(x)ψ~j,k(x)dx,j=j0,…,M,k=1-m,…,2M-m,
where φ~j0,k(x) and ψ~j,k(x) are dual functions of φj0,k and ψj,k, respectively. These can be obtained by linear combinations of φj0,k, k=1-m,…,2j0-1, and ψj,k, j=j0,…,M, k=1-m,…,2M-m, as follows. Let
(20)Φ=[φj0,1-m,…,φj0,2j0-1]T,(21)Ψ-=[ψj0,1-m,…,ψj0,2j0-m,…,ψM,1-m,…,ψM,2M-m]T.
Using (11), (20), (12)-(13), and (21), we get
(22)∫01ΦΦTdx=P1,∫01Ψ-Ψ-Tdx=P2.
Suppose that Φ~ and Ψ-~ are the dual functions of Φ and Ψ-, respectively; then
(23)∫01Φ~ΦTdx=I1,∫01Ψ-~Ψ-Tdx=I2,(24)Φ~=P1-1Φ,Ψ-~=P2-1Ψ-.


#### 3.1.3. Application ofB-Spline Wavelet Method

In this section, linear Fredholm integral equation of the second kind of form (7) has been solved by using B-spline wavelets. For this, we use (16) to approximate y(x) as
(25)y(x)=CTΨ(x),
where Ψ(x) is defined in (18) and C is (2M+1+m-1)×1 unknown vector defined similarly as in (17). We also expand f(x) and K(x,t) by B-spline dual wavelets Ψ~ defined in (24) as
(26)f(x)=C1TΨ~(x),K(x,t)=Ψ~T(t)ΘΨ~(x),
where
(27)Θi,j=∫01[∫01K(x,t)Ψi(t)dt]Ψj(x)dx.
From (26) and (25), we get
(28)∫01K(x,t)y(t)dt=∫01CTΨ(t)Ψ~T(t)ΘΨ~(x)dt=CTΘΨ~(x)
since
(29)∫01Ψ(t)Ψ~T(t)dt=I.
By applying (25)–(28) in (7) we have
(30)CTΨ(x)-CTΘΨ~(x)=C1TΨ~(x).
By multiplying both sides of (30) with ΨT(x) from the right and integrating both sides with respect to x from 0 to 1, we get
(31)CTP-CTΘ=C1T,
since
(32)∫01Ψ~(x)ΨT(x)dx=I,
and P is a (2M+1+m-1)×(2M+1+m-1) square matrix given by
(33)P=∫01Ψ(x)ΨT(x)dx=(P100P2).
Consequently, from (31), we get CT=C1T(P-Θ)-1. Hence, we can calculate the solution for y(x)=CTΨ(x).


### 3.2. Method of Moments

#### 3.2.1. Multiresolution Analysis (MRA) and Wavelets [2]

A set of subspaces{Vj}j∈Z is said to be MRA of L2(R) if it possesses the following properties:
(34)Vj⊂Vj+1,∀j∈Ζ,(35)⋃j∈ΖVjisdenseinL2(R),(36)⋂j∈ΖVj=ϕ,(37)f(x)∈Vj⟺f(2x)∈Vj+1,∀j∈Ζ,
where Z denotes the set of integers. Properties (34)–(36) state that {Vj}j∈Z is a nested sequence of subspaces that effectively covers L2(R). That is, every square integrable function can be approximated as closely as desired by a function that belongs to at least one of the subspaces Vj. A function φ∈L2(R) is called a scaling function if it generates the nested sequence of subspaces Vj and satisfies the dilation equation; namely,
(38)φ(x)=∑kpkφ(ax-k),
with pk∈l2 and a being any rational number.For each scalej, since Vj⊂Vj+1, there exists a unique orthogonal complementary subspace Wj of Vj in Vj+1. This subspace Wj is called wavelet subspace and is generated by ψj,k=ψ(2jx-k), where ψ∈L2 is called the wavelet. From the above discussion, these results follow easily:
(39)Vj1∩Vj2=Vj2,j1>j2,Wj1∩Wj2=0,j1≠j2,Vj1∩Wj2=0,j1≤j2.
Some of the important properties relevant to the present analysis are given below [2, 19].(1) Vanishing Moment. A wavelet is said to have a vanishing moment of order m if
(40)∫-∞∞xpψ(x)dx=0;p=0,…,m-1.
All wavelets must satisfy the previously mentioned condition for p=0.(2) Semiorthogonality. The wavelets ψj,k form a semiorthogonal basis if
(41)〈ψj,k,ψs,i〉=0;j≠s;∀j,k,s,i∈Ζ.


#### 3.2.2. Method of Moments for the Solution of Fredholm Integral Equation

In this section, we solve the integral equation of form (7) in interval [0,1] by using linear B-spline wavelets [2]. The unknown function in (7) can be expanded in terms of the scaling and wavelet functions as follows:
(42)y(x)≈∑k=-12j0-1ckφj0,k(x)+∑j=j0M∑k=-12j-2dj,kψj,k(x)=CTΨ(x).
By substituting this expression into (7) and employing the Galerkin method, the following set of linear system of order (2M+1) is generated. The scaling and wavelet functions are used as testing and weighting functions:
(43)(〈φ,φ〉-〈Kφ,φ〉〈ψ,φ〉-〈Kψ,φ〉〈φ,ψ〉-〈Kφ,ψ〉〈ψ,ψ〉-〈Kψ,ψ〉)(CD)=(F1F2),
where
(44)C=[c-1,c0,…,c3]T,D=[d2,-1,…,d2,2,d3,-1,…,d3,6,…,dM,-1,…,dM,2M-2]T,〈φ,φ〉-〈Kφ,φ〉=(∫01φj0,r(x)φj0,i(x)dx-∫01φj0,r(x)∫01K(x,t)φj0,i(t)dtdx)i,r,〈ψ,φ〉-〈Kψ,φ〉=(∫01φj0,r(x)ψk,j(x)dx-∫01φj0,r(x)∫01K(x,t)ψk,j(t)dtdx)r,k,j,〈φ,ψ〉-〈Kφ,ψ〉=(∫01ψs,l(x)φj0,i(x)dx-∫01ψs,l(x)∫01K(x,t)φj0,i(t)dtdx)i,l,s,〈ψ,ψ〉-〈Kψ,ψ〉=(∫01ψs,l(x)ψk,j(x)dx-∫01ψs,l(x)∫01K(x,t)ψk,j(t)dtdx)l,s,k,j,F1=∫01f(x)φj0,r(x)dx,F2=∫01f(x)ψs,l(x)dx,
and the subscripts i,r,k,j,l, and s assume values as given below:
(45)i,r=-1,…,2j0-1,l,k=j0,…,M,s,j=-1,…,2M-2.
In fact, the entries with significant magnitude are in the 〈Kφ,φ〉-〈φ,φ〉 and 〈Kψ,ψ〉-〈ψ,ψ〉 submatrices which are of order (2j0+1) and (2M+1+1), respectively.


### 3.3. Variational Iteration Method [3–5]

In this section, Fredholm integral equation of second kind given in (7) has been considered for solving (7) by variational iteration method. First, we have to take the partial derivative of (7) with respect to x yielding
(46)Y′(x)=f′(x)+∫01K′(x,t)y(t)dt.
We apply variation iteration method for (46). According to this method, correction functional can be defined as
(47)yn+1(x)=yn(x)+∫0xλ(ξ)(yn′(ξ)-f′(ξ)-∫abK′(ξ,t)y~n(t)dt)dξ,
where λ(ξ) is a general Lagrange multiplier which can be identified optimally by the variational theory, the subscript n denotes the nth order approximation, and y~n is considered as a restricted variation; that is, δy~n=0. The successive approximations yn(x), n≥1 for the solution y(x) can be readily obtained after determining the Lagrange multiplier and selecting an appropriate initial function y0(x). Consequently the approximate solution may be obtained by using
(48)y(x)=limn→∞yn(x).
To make the above correction functional stationary, we have
(49)δyn+1(x)=δyn(x)+δ∫0xλ(ξ)(∫abK′(ξ,t)y~n(t)dtyn′(ξ)-f′(ξ)-∫abK′(ξ,t)y~n(t)dt)dξ=δyn(x)+∫0xλ(ξ)δ(yn′(ξ))dξ=δyn(x)+λδyn|ξ=x-∫0xλ′(ξ)δyn(ξ)dξ.
Under stationary condition,
(50)δyn+1=0
implies the following Euler Lagrange equation:
(51)λ′(ξ)=0,
with the following natural boundary condition:
(52)1+λ(ξ)|ξ=x=0.
Solving (51), along with boundary condition (52), we get the general Lagrange multiplier
(53)λ=-1.
Substituting the identified Lagrange multiplier into (47) results in the following iterative scheme:
(54)yn+1(x)=yn(x)-∫0x(yn′(ξ)-f′(ξ)-∫abK′(ξ,t)y~n(t)dt)dξ,n≥0.
By starting with initial approximate function y0(x)=f(x) (say), we can determine the approximate solution y(x) of (7).


## 3.1.B-Spline Wavelet Method

### 3.1.1.B-Spline Scaling and Wavelet Functions on the Interval [0,1]

Semiorthogonal wavelets usingB-spline are specially constructed for the bounded interval and this wavelet can be represented in a closed form. This provides a compact support. Semiorthogonal wavelets form the basis in the space L2(R).Using this basis, an arbitrary function inL2(R) can be expressed as the wavelet series. For the finite interval [0,1], the wavelet series cannot be completely presented by using this basis. This is because supports of some basis are truncated at the left or right end points of the interval. Hence, a special basis has to be introduced into the wavelet expansion on the finite interval. These functions are referred to as the boundary scaling functions and boundary wavelet functions.Letm and n be two positive integers and let
(8)a=x-m+1=⋯=x0<x1<⋯<xn=xn+1=⋯=xn+m-1=b
be an equally spaced knots sequence. The functions
(9)Bm,j,X(x)=x-xjxj+m-1-xjBm-1,j,X(x)+xj+m-xxj+m-xj+1Bm-1,j+1,X(x),j=-m+1,…,n-1,B1,j,X(x)={1,x∈[xj,xj+1),0,otherwise,
are called cardinal B-spline functions of order m≥2 for the knot sequence X={xi}i=-m+1n+m-1 and SuppBm,j,X(x)=[xj,xj+m]∩[a,b].By considering the interval[a,b]=[0,1], at any level j∈Ζ+, the discretization step is 2-j, and this generates n=2j number of segments in [0,1] with knot sequence
(10)X(j)={x-m+1(j)=⋯=x0(j)=0,xk(j)=k2j,k=1,…,n-1,xn(j)=⋯=xn+m-1(j)=1.
Let j0 be the level for which 2j0≥2m-1; for each level, j≥j0, the scaling functions of order m can be defined as follows in [2]:
(11)φm,j,i(x)={Bm,j0,i(2j-j0x)i=-m+1,…,-1,Bm,j0,2j-m-i(1-2j-j0x)i=2j-m+1,…,2j-1,Bm,j0,0(2j-j0x-2-j0i)i=0,…,2j-m.
And the two scale relations for the m-order semiorthogonal compactly supported B-wavelet functions are defined as follows:
(12)ψm,j,i-m=∑k=i2i+2m-2qi,kBm,j,k-m,i=1,…,m-1,ψm,j,i-m=∑k=2i-m2i+2m-2qi,kBm,j,k-m,i=m,…,n-m+1,ψm,j,i-m=∑k=2i-mn+i+m-1qi,kBm,j,k-m,i=n-m+2,…,n,
where qi,k=qk-2i.Hence, there are2(m-1) boundary wavelets and (n-2m+2) inner wavelets in the bounded interval [a,b]. Finally, by considering the level j with j≥j0, the B-wavelet functions in [0,1] can be expressed as follows:
(13)ψm,j,i(x)={ψm,j0,i(2j-j0x)i=-m+1,…,-1,ψm,2j-2m+1-i,i(1-2j-j0x)i=2j-2m+2,…,2j-m,ψm,j0,0(2j-j0x-2-j0i)i=0,…,2j-2m+1.
The scaling functions φm,j,i(x) occupy m segments and the wavelet functions ψm,j,i(x) occupy 2m-1 segments.When the semiorthogonal wavelets are constructed fromB-spline of order m, the lowest octave level j=j0 is determined in [19, 20] by
(14)2j0≥2m-1,
so as to have a minimum of one complete wavelet on the interval [0,1].


### 3.1.2. Function Approximation

A functionf(x) defined over [0,1] may be approximated by B-spline wavelets as [21, 22]
(15)f(x)=∑k=1-m2j0-1cj0,kφj0,k(x)+∑j=j0∞∑k=1-m2j-mdj,kψj,k(x).
If the infinite series in (15) is truncated at M, then (15) can be written as [2]
(16)f(x)≅∑k=1-m2j0-1cj0,kφj0,k(x)+∑j=j0M∑k=1-m2j-mdj,kψj,k(x),
where φ2,k and ψj,k are scaling and wavelets functions, respectively, and C and Ψ are (2M+1+m-1)×1 vectors given by
(17)C=[cj0,1-m,…,cj0,2j0-1,dj0,1-m,…,dj0,2j0-m,…,dM,1-m,…,dM,2M-m]T,(18)Ψ=[φj0,1-m,…,φj0,2j0-1,ψj0,1-m,…,ψj0,2j0-m,…,ψM,1-m,…,ψM,2M-m]T,
with
(19)cj0,k=∫01f(x)φ~j0,k(x)dx,k=1-m,…,2j0-1,kkkkkkkkkkdj,k=∫01f(x)ψ~j,k(x)dx,j=j0,…,M,k=1-m,…,2M-m,
where φ~j0,k(x) and ψ~j,k(x) are dual functions of φj0,k and ψj,k, respectively. These can be obtained by linear combinations of φj0,k, k=1-m,…,2j0-1, and ψj,k, j=j0,…,M, k=1-m,…,2M-m, as follows. Let
(20)Φ=[φj0,1-m,…,φj0,2j0-1]T,(21)Ψ-=[ψj0,1-m,…,ψj0,2j0-m,…,ψM,1-m,…,ψM,2M-m]T.
Using (11), (20), (12)-(13), and (21), we get
(22)∫01ΦΦTdx=P1,∫01Ψ-Ψ-Tdx=P2.
Suppose that Φ~ and Ψ-~ are the dual functions of Φ and Ψ-, respectively; then
(23)∫01Φ~ΦTdx=I1,∫01Ψ-~Ψ-Tdx=I2,(24)Φ~=P1-1Φ,Ψ-~=P2-1Ψ-.


### 3.1.3. Application ofB-Spline Wavelet Method

In this section, linear Fredholm integral equation of the second kind of form (7) has been solved by using B-spline wavelets. For this, we use (16) to approximate y(x) as
(25)y(x)=CTΨ(x),
where Ψ(x) is defined in (18) and C is (2M+1+m-1)×1 unknown vector defined similarly as in (17). We also expand f(x) and K(x,t) by B-spline dual wavelets Ψ~ defined in (24) as
(26)f(x)=C1TΨ~(x),K(x,t)=Ψ~T(t)ΘΨ~(x),
where
(27)Θi,j=∫01[∫01K(x,t)Ψi(t)dt]Ψj(x)dx.
From (26) and (25), we get
(28)∫01K(x,t)y(t)dt=∫01CTΨ(t)Ψ~T(t)ΘΨ~(x)dt=CTΘΨ~(x)
since
(29)∫01Ψ(t)Ψ~T(t)dt=I.
By applying (25)–(28) in (7) we have
(30)CTΨ(x)-CTΘΨ~(x)=C1TΨ~(x).
By multiplying both sides of (30) with ΨT(x) from the right and integrating both sides with respect to x from 0 to 1, we get
(31)CTP-CTΘ=C1T,
since
(32)∫01Ψ~(x)ΨT(x)dx=I,
and P is a (2M+1+m-1)×(2M+1+m-1) square matrix given by
(33)P=∫01Ψ(x)ΨT(x)dx=(P100P2).
Consequently, from (31), we get CT=C1T(P-Θ)-1. Hence, we can calculate the solution for y(x)=CTΨ(x).


## 3.1.1.B-Spline Scaling and Wavelet Functions on the Interval [0,1]

Semiorthogonal wavelets usingB-spline are specially constructed for the bounded interval and this wavelet can be represented in a closed form. This provides a compact support. Semiorthogonal wavelets form the basis in the space L2(R).Using this basis, an arbitrary function inL2(R) can be expressed as the wavelet series. For the finite interval [0,1], the wavelet series cannot be completely presented by using this basis. This is because supports of some basis are truncated at the left or right end points of the interval. Hence, a special basis has to be introduced into the wavelet expansion on the finite interval. These functions are referred to as the boundary scaling functions and boundary wavelet functions.Letm and n be two positive integers and let
(8)a=x-m+1=⋯=x0<x1<⋯<xn=xn+1=⋯=xn+m-1=b
be an equally spaced knots sequence. The functions
(9)Bm,j,X(x)=x-xjxj+m-1-xjBm-1,j,X(x)+xj+m-xxj+m-xj+1Bm-1,j+1,X(x),j=-m+1,…,n-1,B1,j,X(x)={1,x∈[xj,xj+1),0,otherwise,
are called cardinal B-spline functions of order m≥2 for the knot sequence X={xi}i=-m+1n+m-1 and SuppBm,j,X(x)=[xj,xj+m]∩[a,b].By considering the interval[a,b]=[0,1], at any level j∈Ζ+, the discretization step is 2-j, and this generates n=2j number of segments in [0,1] with knot sequence
(10)X(j)={x-m+1(j)=⋯=x0(j)=0,xk(j)=k2j,k=1,…,n-1,xn(j)=⋯=xn+m-1(j)=1.
Let j0 be the level for which 2j0≥2m-1; for each level, j≥j0, the scaling functions of order m can be defined as follows in [2]:
(11)φm,j,i(x)={Bm,j0,i(2j-j0x)i=-m+1,…,-1,Bm,j0,2j-m-i(1-2j-j0x)i=2j-m+1,…,2j-1,Bm,j0,0(2j-j0x-2-j0i)i=0,…,2j-m.
And the two scale relations for the m-order semiorthogonal compactly supported B-wavelet functions are defined as follows:
(12)ψm,j,i-m=∑k=i2i+2m-2qi,kBm,j,k-m,i=1,…,m-1,ψm,j,i-m=∑k=2i-m2i+2m-2qi,kBm,j,k-m,i=m,…,n-m+1,ψm,j,i-m=∑k=2i-mn+i+m-1qi,kBm,j,k-m,i=n-m+2,…,n,
where qi,k=qk-2i.Hence, there are2(m-1) boundary wavelets and (n-2m+2) inner wavelets in the bounded interval [a,b]. Finally, by considering the level j with j≥j0, the B-wavelet functions in [0,1] can be expressed as follows:
(13)ψm,j,i(x)={ψm,j0,i(2j-j0x)i=-m+1,…,-1,ψm,2j-2m+1-i,i(1-2j-j0x)i=2j-2m+2,…,2j-m,ψm,j0,0(2j-j0x-2-j0i)i=0,…,2j-2m+1.
The scaling functions φm,j,i(x) occupy m segments and the wavelet functions ψm,j,i(x) occupy 2m-1 segments.When the semiorthogonal wavelets are constructed fromB-spline of order m, the lowest octave level j=j0 is determined in [19, 20] by
(14)2j0≥2m-1,
so as to have a minimum of one complete wavelet on the interval [0,1].


## 3.1.2. Function Approximation

A functionf(x) defined over [0,1] may be approximated by B-spline wavelets as [21, 22]
(15)f(x)=∑k=1-m2j0-1cj0,kφj0,k(x)+∑j=j0∞∑k=1-m2j-mdj,kψj,k(x).
If the infinite series in (15) is truncated at M, then (15) can be written as [2]
(16)f(x)≅∑k=1-m2j0-1cj0,kφj0,k(x)+∑j=j0M∑k=1-m2j-mdj,kψj,k(x),
where φ2,k and ψj,k are scaling and wavelets functions, respectively, and C and Ψ are (2M+1+m-1)×1 vectors given by
(17)C=[cj0,1-m,…,cj0,2j0-1,dj0,1-m,…,dj0,2j0-m,…,dM,1-m,…,dM,2M-m]T,(18)Ψ=[φj0,1-m,…,φj0,2j0-1,ψj0,1-m,…,ψj0,2j0-m,…,ψM,1-m,…,ψM,2M-m]T,
with
(19)cj0,k=∫01f(x)φ~j0,k(x)dx,k=1-m,…,2j0-1,kkkkkkkkkkdj,k=∫01f(x)ψ~j,k(x)dx,j=j0,…,M,k=1-m,…,2M-m,
where φ~j0,k(x) and ψ~j,k(x) are dual functions of φj0,k and ψj,k, respectively. These can be obtained by linear combinations of φj0,k, k=1-m,…,2j0-1, and ψj,k, j=j0,…,M, k=1-m,…,2M-m, as follows. Let
(20)Φ=[φj0,1-m,…,φj0,2j0-1]T,(21)Ψ-=[ψj0,1-m,…,ψj0,2j0-m,…,ψM,1-m,…,ψM,2M-m]T.
Using (11), (20), (12)-(13), and (21), we get
(22)∫01ΦΦTdx=P1,∫01Ψ-Ψ-Tdx=P2.
Suppose that Φ~ and Ψ-~ are the dual functions of Φ and Ψ-, respectively; then
(23)∫01Φ~ΦTdx=I1,∫01Ψ-~Ψ-Tdx=I2,(24)Φ~=P1-1Φ,Ψ-~=P2-1Ψ-.


## 3.1.3. Application ofB-Spline Wavelet Method

In this section, linear Fredholm integral equation of the second kind of form (7) has been solved by using B-spline wavelets. For this, we use (16) to approximate y(x) as
(25)y(x)=CTΨ(x),
where Ψ(x) is defined in (18) and C is (2M+1+m-1)×1 unknown vector defined similarly as in (17). We also expand f(x) and K(x,t) by B-spline dual wavelets Ψ~ defined in (24) as
(26)f(x)=C1TΨ~(x),K(x,t)=Ψ~T(t)ΘΨ~(x),
where
(27)Θi,j=∫01[∫01K(x,t)Ψi(t)dt]Ψj(x)dx.
From (26) and (25), we get
(28)∫01K(x,t)y(t)dt=∫01CTΨ(t)Ψ~T(t)ΘΨ~(x)dt=CTΘΨ~(x)
since
(29)∫01Ψ(t)Ψ~T(t)dt=I.
By applying (25)–(28) in (7) we have
(30)CTΨ(x)-CTΘΨ~(x)=C1TΨ~(x).
By multiplying both sides of (30) with ΨT(x) from the right and integrating both sides with respect to x from 0 to 1, we get
(31)CTP-CTΘ=C1T,
since
(32)∫01Ψ~(x)ΨT(x)dx=I,
and P is a (2M+1+m-1)×(2M+1+m-1) square matrix given by
(33)P=∫01Ψ(x)ΨT(x)dx=(P100P2).
Consequently, from (31), we get CT=C1T(P-Θ)-1. Hence, we can calculate the solution for y(x)=CTΨ(x).


## 3.2. Method of Moments

### 3.2.1. Multiresolution Analysis (MRA) and Wavelets [2]

A set of subspaces{Vj}j∈Z is said to be MRA of L2(R) if it possesses the following properties:
(34)Vj⊂Vj+1,∀j∈Ζ,(35)⋃j∈ΖVjisdenseinL2(R),(36)⋂j∈ΖVj=ϕ,(37)f(x)∈Vj⟺f(2x)∈Vj+1,∀j∈Ζ,
where Z denotes the set of integers. Properties (34)–(36) state that {Vj}j∈Z is a nested sequence of subspaces that effectively covers L2(R). That is, every square integrable function can be approximated as closely as desired by a function that belongs to at least one of the subspaces Vj. A function φ∈L2(R) is called a scaling function if it generates the nested sequence of subspaces Vj and satisfies the dilation equation; namely,
(38)φ(x)=∑kpkφ(ax-k),
with pk∈l2 and a being any rational number.For each scalej, since Vj⊂Vj+1, there exists a unique orthogonal complementary subspace Wj of Vj in Vj+1. This subspace Wj is called wavelet subspace and is generated by ψj,k=ψ(2jx-k), where ψ∈L2 is called the wavelet. From the above discussion, these results follow easily:
(39)Vj1∩Vj2=Vj2,j1>j2,Wj1∩Wj2=0,j1≠j2,Vj1∩Wj2=0,j1≤j2.
Some of the important properties relevant to the present analysis are given below [2, 19].(1) Vanishing Moment. A wavelet is said to have a vanishing moment of order m if
(40)∫-∞∞xpψ(x)dx=0;p=0,…,m-1.
All wavelets must satisfy the previously mentioned condition for p=0.(2) Semiorthogonality. The wavelets ψj,k form a semiorthogonal basis if
(41)〈ψj,k,ψs,i〉=0;j≠s;∀j,k,s,i∈Ζ.


### 3.2.2. Method of Moments for the Solution of Fredholm Integral Equation

In this section, we solve the integral equation of form (7) in interval [0,1] by using linear B-spline wavelets [2]. The unknown function in (7) can be expanded in terms of the scaling and wavelet functions as follows:
(42)y(x)≈∑k=-12j0-1ckφj0,k(x)+∑j=j0M∑k=-12j-2dj,kψj,k(x)=CTΨ(x).
By substituting this expression into (7) and employing the Galerkin method, the following set of linear system of order (2M+1) is generated. The scaling and wavelet functions are used as testing and weighting functions:
(43)(〈φ,φ〉-〈Kφ,φ〉〈ψ,φ〉-〈Kψ,φ〉〈φ,ψ〉-〈Kφ,ψ〉〈ψ,ψ〉-〈Kψ,ψ〉)(CD)=(F1F2),
where
(44)C=[c-1,c0,…,c3]T,D=[d2,-1,…,d2,2,d3,-1,…,d3,6,…,dM,-1,…,dM,2M-2]T,〈φ,φ〉-〈Kφ,φ〉=(∫01φj0,r(x)φj0,i(x)dx-∫01φj0,r(x)∫01K(x,t)φj0,i(t)dtdx)i,r,〈ψ,φ〉-〈Kψ,φ〉=(∫01φj0,r(x)ψk,j(x)dx-∫01φj0,r(x)∫01K(x,t)ψk,j(t)dtdx)r,k,j,〈φ,ψ〉-〈Kφ,ψ〉=(∫01ψs,l(x)φj0,i(x)dx-∫01ψs,l(x)∫01K(x,t)φj0,i(t)dtdx)i,l,s,〈ψ,ψ〉-〈Kψ,ψ〉=(∫01ψs,l(x)ψk,j(x)dx-∫01ψs,l(x)∫01K(x,t)ψk,j(t)dtdx)l,s,k,j,F1=∫01f(x)φj0,r(x)dx,F2=∫01f(x)ψs,l(x)dx,
and the subscripts i,r,k,j,l, and s assume values as given below:
(45)i,r=-1,…,2j0-1,l,k=j0,…,M,s,j=-1,…,2M-2.
In fact, the entries with significant magnitude are in the 〈Kφ,φ〉-〈φ,φ〉 and 〈Kψ,ψ〉-〈ψ,ψ〉 submatrices which are of order (2j0+1) and (2M+1+1), respectively.


## 3.2.1. Multiresolution Analysis (MRA) and Wavelets [2]

A set of subspaces{Vj}j∈Z is said to be MRA of L2(R) if it possesses the following properties:
(34)Vj⊂Vj+1,∀j∈Ζ,(35)⋃j∈ΖVjisdenseinL2(R),(36)⋂j∈ΖVj=ϕ,(37)f(x)∈Vj⟺f(2x)∈Vj+1,∀j∈Ζ,
where Z denotes the set of integers. Properties (34)–(36) state that {Vj}j∈Z is a nested sequence of subspaces that effectively covers L2(R). That is, every square integrable function can be approximated as closely as desired by a function that belongs to at least one of the subspaces Vj. A function φ∈L2(R) is called a scaling function if it generates the nested sequence of subspaces Vj and satisfies the dilation equation; namely,
(38)φ(x)=∑kpkφ(ax-k),
with pk∈l2 and a being any rational number.For each scalej, since Vj⊂Vj+1, there exists a unique orthogonal complementary subspace Wj of Vj in Vj+1. This subspace Wj is called wavelet subspace and is generated by ψj,k=ψ(2jx-k), where ψ∈L2 is called the wavelet. From the above discussion, these results follow easily:
(39)Vj1∩Vj2=Vj2,j1>j2,Wj1∩Wj2=0,j1≠j2,Vj1∩Wj2=0,j1≤j2.
Some of the important properties relevant to the present analysis are given below [2, 19].(1) Vanishing Moment. A wavelet is said to have a vanishing moment of order m if
(40)∫-∞∞xpψ(x)dx=0;p=0,…,m-1.
All wavelets must satisfy the previously mentioned condition for p=0.(2) Semiorthogonality. The wavelets ψj,k form a semiorthogonal basis if
(41)〈ψj,k,ψs,i〉=0;j≠s;∀j,k,s,i∈Ζ.


## 3.2.2. Method of Moments for the Solution of Fredholm Integral Equation

In this section, we solve the integral equation of form (7) in interval [0,1] by using linear B-spline wavelets [2]. The unknown function in (7) can be expanded in terms of the scaling and wavelet functions as follows:
(42)y(x)≈∑k=-12j0-1ckφj0,k(x)+∑j=j0M∑k=-12j-2dj,kψj,k(x)=CTΨ(x).
By substituting this expression into (7) and employing the Galerkin method, the following set of linear system of order (2M+1) is generated. The scaling and wavelet functions are used as testing and weighting functions:
(43)(〈φ,φ〉-〈Kφ,φ〉〈ψ,φ〉-〈Kψ,φ〉〈φ,ψ〉-〈Kφ,ψ〉〈ψ,ψ〉-〈Kψ,ψ〉)(CD)=(F1F2),
where
(44)C=[c-1,c0,…,c3]T,D=[d2,-1,…,d2,2,d3,-1,…,d3,6,…,dM,-1,…,dM,2M-2]T,〈φ,φ〉-〈Kφ,φ〉=(∫01φj0,r(x)φj0,i(x)dx-∫01φj0,r(x)∫01K(x,t)φj0,i(t)dtdx)i,r,〈ψ,φ〉-〈Kψ,φ〉=(∫01φj0,r(x)ψk,j(x)dx-∫01φj0,r(x)∫01K(x,t)ψk,j(t)dtdx)r,k,j,〈φ,ψ〉-〈Kφ,ψ〉=(∫01ψs,l(x)φj0,i(x)dx-∫01ψs,l(x)∫01K(x,t)φj0,i(t)dtdx)i,l,s,〈ψ,ψ〉-〈Kψ,ψ〉=(∫01ψs,l(x)ψk,j(x)dx-∫01ψs,l(x)∫01K(x,t)ψk,j(t)dtdx)l,s,k,j,F1=∫01f(x)φj0,r(x)dx,F2=∫01f(x)ψs,l(x)dx,
and the subscripts i,r,k,j,l, and s assume values as given below:
(45)i,r=-1,…,2j0-1,l,k=j0,…,M,s,j=-1,…,2M-2.
In fact, the entries with significant magnitude are in the 〈Kφ,φ〉-〈φ,φ〉 and 〈Kψ,ψ〉-〈ψ,ψ〉 submatrices which are of order (2j0+1) and (2M+1+1), respectively.


## 3.3. Variational Iteration Method [3–5]

In this section, Fredholm integral equation of second kind given in (7) has been considered for solving (7) by variational iteration method. First, we have to take the partial derivative of (7) with respect to x yielding
(46)Y′(x)=f′(x)+∫01K′(x,t)y(t)dt.
We apply variation iteration method for (46). According to this method, correction functional can be defined as
(47)yn+1(x)=yn(x)+∫0xλ(ξ)(yn′(ξ)-f′(ξ)-∫abK′(ξ,t)y~n(t)dt)dξ,
where λ(ξ) is a general Lagrange multiplier which can be identified optimally by the variational theory, the subscript n denotes the nth order approximation, and y~n is considered as a restricted variation; that is, δy~n=0. The successive approximations yn(x), n≥1 for the solution y(x) can be readily obtained after determining the Lagrange multiplier and selecting an appropriate initial function y0(x). Consequently the approximate solution may be obtained by using
(48)y(x)=limn→∞yn(x).
To make the above correction functional stationary, we have
(49)δyn+1(x)=δyn(x)+δ∫0xλ(ξ)(∫abK′(ξ,t)y~n(t)dtyn′(ξ)-f′(ξ)-∫abK′(ξ,t)y~n(t)dt)dξ=δyn(x)+∫0xλ(ξ)δ(yn′(ξ))dξ=δyn(x)+λδyn|ξ=x-∫0xλ′(ξ)δyn(ξ)dξ.
Under stationary condition,
(50)δyn+1=0
implies the following Euler Lagrange equation:
(51)λ′(ξ)=0,
with the following natural boundary condition:
(52)1+λ(ξ)|ξ=x=0.
Solving (51), along with boundary condition (52), we get the general Lagrange multiplier
(53)λ=-1.
Substituting the identified Lagrange multiplier into (47) results in the following iterative scheme:
(54)yn+1(x)=yn(x)-∫0x(yn′(ξ)-f′(ξ)-∫abK′(ξ,t)y~n(t)dt)dξ,n≥0.
By starting with initial approximate function y0(x)=f(x) (say), we can determine the approximate solution y(x) of (7).


## 4. Numerical Methods for System of Linear Fredholm Integral Equations of Second Kind

Consider the system of linear Fredholm integral equations of second kind of the following form:(55)∑j=1nyj(x)=fi(x)+∑j=1n∫01Ki,j(x,t)yj(t)dt,kkkkkkkkkkkkkkkkkkkkkkki=1,2,…,n,
where fi(x) and Ki,j(x,t) are known functions and yj(x) are the unknown functions for i,j=1,2,…,n.


### 4.1. Application of Haar Wavelet Method [9]

In this section, an efficient algorithm for solving Fredholm integral equations with Haar wavelets is analyzed. The present algorithm takes the following essential strategy. The Haar wavelet is first used to decompose integral equations into algebraic systems of linear equations, which are then solved by collocation methods.


#### 4.1.1. Haar Wavelets

The compact set of scale functions is chosen as(56)h0={1,0≤x<1,0,others.
The mother wavelet function is defined as
(57)h1(x)={1,0≤x<12,-1,12≤x<1,0,others.
The family of wavelet functions generated by translation and dilation of h1(x) are given by
(58)hn(x)=h1(2jx-k),
where n=2j+k, j≥0, 0≤k<2j.Mutual orthogonalities of all Haar wavelets can be expressed as(59)∫01hm(x)hn(x)dx=2-jδmn={2-j,m=n=2j+k,0,m≠n.


#### 4.1.2. Function Approximation

An arbitrary functiony(x)∈L2[0,1) can be expanded into the following Haar series:
(60)y(x)=∑n=0+∞cnhn(x),
where the coefficients cn are given by
(61)kkkkkkicn=2j∫01y(x)hn(x)dx,n=2j+k,j≥0,0≤k<2j.
In particular, c0=∫01y(x)dx.The previously mentioned expression in (60) can be approximately represented with finite terms as follows:
(62)y(x)≈∑n=0m-1cnhn(x)=C(m)Th(m)(x),
where the coefficient vector C(m)T and the Haar function vector h(m)(x) are, respectively, defined as
(63)C(m)T=[c0,c1,…,cm-1],m=2j,h(m)(x)=[h0(x),h1(x),…,hm-1(x)]T,m=2j.
The Haar expansion for function K(x,t) of order m is defined as follows:
(64)K(x,t)≈∑u=0m-1∑v=0m-1auvhv(x)hu(t),
where auv=2i+q∬01K(x,t)hv(x)hu(t)dxdt, u=2i+j, v=2q+r, i,q≥0.From (62) and (64), we obtain
(65)K(x,t)≈h(m)T(t)Kh(m)(x),
where
(66)K=(auv)m×mT.


#### 4.1.3. Operational Matrices of Integration

We define(67)H(m)=[h(m)(12m),h(m)(32m),…,h(m)(2m-12m)],
where H(1)=[1], H(2)=[111-1].Then, form=4, the corresponding matrix can be represented as
(68)H(4)=[h(4)(18),h(4)(38),…,h(4)(78)]=[111111-1-11-100001-1].
The integration of the Haar function vector h(m)(t) is
(69)∫0xh(m)(t)dt=P(m)h(m)(x),kkkkkkkkkkkkkkix∈[0,1),
where P(m) is the operational matrix of order m, and
(70)P(1)=[12],P(m)=12m[2mP(m/2)-H(m/2)H(m/2)-10].
By recursion of the above formula, we obtain
(71)P(2)=14[2-110],P(4)=116[8-4-2-240-2211001-100],P(8)=164[32-16-8-8-4-4-4-4160-88-4-4444400-4400440000-441120000011-2000001-10200001-10-20000].
Therefore, we get
(72)H(m)-1=(1m)H(m)T×diag(1,1,2,2,22,…,22︸22,…,2α-1,…,2α-1︸2α-1),
where m=2α and α is a positive integer.The inner product of two Haar functions can be represented as(73)∫01h(m)(t)h(m)T(t)dt=D,
where
(74)D=diag(1,1,1/2,1/2,1/22,…,1/22︸22,…,1/2α-1,…,1/2α-1︸2α-1).


#### 4.1.4. Haar Wavelet Solution for Fredholm Integral Equations System [9]

Consider the following Fredholm integral equations system defined in (55):
(75)∑j=1myj(x)=fi(x)+∑j=1m∫01Ki,j(x,t)yj(t)dt,kkkkkkkkkkkkkkkkkkkkkkii=1,2,…,m.The Haar series ofyj(x) and Ki,j(x,t), i=1,2,…,m; j=1,2,…,m are, respectively, expanded as
(76)yj(x)≈YjTh(m)(x),j=1,2,…,m,Ki,j(x,t)≈h(m)T(t)Ki,jh(m)(x),kkkkkkkkkkkii,j=1,2,…,m.
Substituting (76) into (75), we get
(77)∑j=1mYjTh(m)(x)=fi(x)+∑j=1m∫01YjTh(m)(t)h(m)T(t)Ki,jh(m)(x)dt,kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkii=1,2,…,m.
From (77) and (73), we get
(78)∑j=1mYjTh(m)(x)=fi(x)+∑j=1mYjTDKi,jh(m)(x),kkkkkkkkkkkkkkkkkkkkkkkii=1,2,…,m.
Interpolating m collocation points, that is, {xi}i=1m, in the interval [0,1] leads to the following algebraic system of equations:
(79)∑j=1mYjTh(m)(xi)=fi(xi)+∑j=1mYjTDKi,jh(m)(xi),kkkkkkkkkkkkkkkkkkkkkkkkii=1,2,…,m.
Hence, Yj, j=1,2,…,m can be computed by solving the above algebraic system of equations and consequently the solutions yj(x)≈YjTh(m)(x), j=1,2,…,m.


### 4.2. Taylor Series Expansion Method

In this section, we present Taylor series expansion method for solving Fredholm integral equations system of second kind [7]. This method reduces the system of integral equations to a linear system of ordinary differential equation. After including boundary conditions, this system reduces to a system of equations that can be solved easily by any usual methods.Consider the second kind Fredholm integral equations system defined in (55) as follows:
(80)yi(x)=fi(x)+∑j=1n∫01Ki,j(x,t)yj(t)dt,kkkkkkkkkki=1,2,…,n,0≤x≤1.
A Taylor series expansion can be made for the solution of yj(t) in the integral equation (80):
(81)yj(t)=yj(x)+yj′(x)(t-x)+⋯+1m!yj(m)(x)(t-x)m+E(t),
where E(t) denotes the error between yj(t) and its Taylor series expansion in (81).If we use the firstm term of Taylor series expansion and neglect the term containing E(t), that is, ∫01∑j=1nKi,j(x,t)E(t)dt, then, substituting (81) for yj(t) into the integral in (80), we have
(82)yi(x)≈fi(x)+∑j=1n∫01Ki,j(x,t)∑r=0m1r!(t-x)ryj(r)(x)dt,kkkkkkkkkkkkkkkkkkkkkkii=1,2,…,n,yi(x)≈fi(x)+∑j=1n∑r=0m1r!yj(r)(x)∫01Ki,j(x,t)(t-x)rdt,kkkkkkkkkkkkkkkkkkkkkkii=1,2,…,n,(83)yi(x)-∑j=1n∑r=0m1r!yj(r)(x)[∫01Ki,j(x,t)(t-x)rdt]≈fi(x),i=1,2,…,n.
Equation (83) becomes a linear system of ordinary differential equations that we have to solve. For solving the linear system of ordinary differential equations (83), we require an appropriate number of boundary conditions.In order to construct boundary conditions, we first differentiates times both sides of (80) with respect to x; that is,
(84)yi(s)(x)=fi(s)(x)+∑j=1n∫01Ki,j(s)(x,t)yj(t)dt,kkkkkkkkkkki=1,2,…,n,s=1,2,…,m,
where Ki,j(s)(x,t)=∂(s)Ki,j(x,t)/∂x(s), s=1,2,…,m.Applying the mean value theorem for integral in (84), we have
(85)yi(s)(x)-[∑j=1n∫01Ki,j(s)(x,t)dt]yj(x)≈fi(s)(x),kkkkkkkkkkkkki=1,2,…,n,s=1,2,…,m.
Now (83) combined with (85) becomes a linear system of algebraic equations that can be solved analytically or numerically.


### 4.3. Block-Pulse Functions for the Solution of Fredholm Integral Equation

In this section, Block-Pulse functions (BPF) have been utilized for the solution of system of Fredholm integral equations [6].Anm-set of BPF is defined as follows:
(86)Φi(t)={1,(i-1)Tm≤t<iTm,0,otherwise
with t∈[0,T), T/m=h and i=1,2,…,m.


#### 4.3.1. Properties of BPF

(1) Disjointness. One has(87)Φi(t)Φj(t)={Φi(t),i=j;0,i≠j,, i,j=1,2,…,m. This property is obtained from definition of BPF.(2) Orthogonality.
One has(88)∫0TΦi(t)Φj(t)dt={h,i=j;0,i≠j,t∈[0,T), i,j=1,2,…,m. This property is obtained from the disjointness property.(3) Completeness.For every f∈L2,{Φ} is complete; if ∫Φf=0 then f=0 almost everywhere. Because of completeness of {Φ}, we have
(89)∫0Tf2(t)dt=∑i=1∞fi2∥Φi(t)∥2
for every real bounded function f(t) which is square integrable in the interval t∈[0,T) and fi=(1/h)f(t)Φi(t)dt.


#### 4.3.2. Function Approximation

The orthogonality property of BPF is the basis of expanding functions into their Block-Pulse series. For everyf(t)∈L2(R),
(90)f(t)=∑i=1mfiΦi(t),
where fi is the coefficient of Block-Pulse function, with respect to ith Block-Pulse function Φi(t).The criterion of this approximation is that mean square error betweenf(t) and its expansion is minimum
(91)ε=1T∫0T(f(t)-∑j=1mfjΦj(t))2dt
so that we can evaluate Block-Pulse coefficients. 										(92)Now∂ε∂fi=-2T∫0T(f(t)-∑j=1mfjΦj(t))Φi(t)dt=0,⟹fi=1h∫0Tf(t)Φi(t)dt(using  orthogonal  property).
In the matrix form, we obtain the following from (90) as follow:
(93)f(t)=∑i=1mfiΦi(t)=FTΦ(t)=ΦTFwhereF=[f1,f2,…,fm]T,Φ(t)=[Φ1(t),Φ2(t),…,Φm(t)]T.
Now let K(t,s) be two-variable function defined on t∈[0,T) and s∈[0,1); then K(t,s) can be expanded to BPF as
(94)K(t,s)=ΦT(t)KΨ(s),
where Φ(t) and Ψ(s) are m1 and m2 dimensional Block-Pulse function vectors and k is a m1×m2 Block-Pulse coefficient matrix.There are two different cases of multiplication of two BPF. The first case is(95)Φ(t)ΦT(t)=(Φ1(t)0⋯00Φ2(t)⋯0⋮⋮⋱⋮00⋯Φm(t)).
It is obtained from disjointness property of BPF. It is a diagonal matrix with m Block-Pulse functions.The second case is(96)ΦT(t)Φ(t)=1
because ∑i=1m(Φi(t))2=∑i=1mΦi(t)=1.Operational Matrix of Integration. BPF integration property can be expressed by an operational equation as
(97)∫0TΦ(t)dt=PΦ(t),
where
(98)Φ(t)=[Φ1(t),Φ2(t),…,Φm(t)]T.
A general formula for Pm×m can be written as
(99)P=12(122⋯2012⋯2001⋯2⋮⋮⋮⋱⋮000⋯1).
By using this matrix, we can express the integral of a function f(t) into its Block-Pulse series
(100)∫0tf(t)dt=∫0tFTΦ(t)dt=FTPΦ(t).


#### 4.3.3. Solution for Linear Integral Equations System

Consider the integral equations system from (55) as follows:
(101)∑j=1nyj(x)=fi(x)+∑j=1n∫αβKi,j(x,t)yj(t)dt,kkkkkkkkkkkikki=1,2,…,n.
Block-Pulse coefficients of yj(x), j=1,2,…,n in the interval x∈[α,β) can be determined from the known functions fi(x), i=1,2,…,n and the kernels Ki,j(x,t), i,j=1,2,…n. Usually, we consider α=0 to facilitie the use of Block-Pulse functions. In case α≠0, we set X=((x-α)/(β-α))T, where T=mh.We approximatefi(x), yj(x), Ki,j(x,t) by its BPF as follows:
(102)fi(x)≈FiTΦ(x),yj(x)≈YjTΦ(x),Ki,j(x,t)≈ΦT(t)Ki,jΦ(x),
where Fi, Yj, and Ki,j are defined in Section 4.3.2, and substituting (102) into (101), we have
(103)∑j=1nYjTΦ(x)=FiTΦ(x)+∑j=1n∫0mhYjTΦ(t)ΦT(t)Ki,jΦ(x)dt,kkkkkkkkkkkkkkkkkii=1,2,…,n,(104)∑j=1nYjTΦ(x)=FiTΦ(x)+∑j=1nYjThIKi,jΦ(x),kkkkkkkikkkki=1,2,…,n,
since
(105)∫0mhΦ(t)ΦT(t)dt=hI.
From (104), we get
(106)∑j=1n(I-hKi,jT)Yj=Fi,i=1,2,…,n.
Set Ai,j=I-hKi,jT; then we have from (106)
(107)∑j=1nAi,jYj=Fi,i=1,2,…,n
which is a linear system
(108)(A11A12...A1nA21A22...A2n............An1An2...Ann)(Y1Y2...Yn)=(F1F2...Fn).
After solving the above system we can find Yj, j=1,2,…,n and hence obtain the solutions yj=ΦTYj, j=1,2,…,n.


## 4.1. Application of Haar Wavelet Method [9]

In this section, an efficient algorithm for solving Fredholm integral equations with Haar wavelets is analyzed. The present algorithm takes the following essential strategy. The Haar wavelet is first used to decompose integral equations into algebraic systems of linear equations, which are then solved by collocation methods.


### 4.1.1. Haar Wavelets

The compact set of scale functions is chosen as(56)h0={1,0≤x<1,0,others.
The mother wavelet function is defined as
(57)h1(x)={1,0≤x<12,-1,12≤x<1,0,others.
The family of wavelet functions generated by translation and dilation of h1(x) are given by
(58)hn(x)=h1(2jx-k),
where n=2j+k, j≥0, 0≤k<2j.Mutual orthogonalities of all Haar wavelets can be expressed as(59)∫01hm(x)hn(x)dx=2-jδmn={2-j,m=n=2j+k,0,m≠n.


### 4.1.2. Function Approximation

An arbitrary functiony(x)∈L2[0,1) can be expanded into the following Haar series:
(60)y(x)=∑n=0+∞cnhn(x),
where the coefficients cn are given by
(61)kkkkkkicn=2j∫01y(x)hn(x)dx,n=2j+k,j≥0,0≤k<2j.
In particular, c0=∫01y(x)dx.The previously mentioned expression in (60) can be approximately represented with finite terms as follows:
(62)y(x)≈∑n=0m-1cnhn(x)=C(m)Th(m)(x),
where the coefficient vector C(m)T and the Haar function vector h(m)(x) are, respectively, defined as
(63)C(m)T=[c0,c1,…,cm-1],m=2j,h(m)(x)=[h0(x),h1(x),…,hm-1(x)]T,m=2j.
The Haar expansion for function K(x,t) of order m is defined as follows:
(64)K(x,t)≈∑u=0m-1∑v=0m-1auvhv(x)hu(t),
where auv=2i+q∬01K(x,t)hv(x)hu(t)dxdt, u=2i+j, v=2q+r, i,q≥0.From (62) and (64), we obtain
(65)K(x,t)≈h(m)T(t)Kh(m)(x),
where
(66)K=(auv)m×mT.


### 4.1.3. Operational Matrices of Integration

We define(67)H(m)=[h(m)(12m),h(m)(32m),…,h(m)(2m-12m)],
where H(1)=[1], H(2)=[111-1].Then, form=4, the corresponding matrix can be represented as
(68)H(4)=[h(4)(18),h(4)(38),…,h(4)(78)]=[111111-1-11-100001-1].
The integration of the Haar function vector h(m)(t) is
(69)∫0xh(m)(t)dt=P(m)h(m)(x),kkkkkkkkkkkkkkix∈[0,1),
where P(m) is the operational matrix of order m, and
(70)P(1)=[12],P(m)=12m[2mP(m/2)-H(m/2)H(m/2)-10].
By recursion of the above formula, we obtain
(71)P(2)=14[2-110],P(4)=116[8-4-2-240-2211001-100],P(8)=164[32-16-8-8-4-4-4-4160-88-4-4444400-4400440000-441120000011-2000001-10200001-10-20000].
Therefore, we get
(72)H(m)-1=(1m)H(m)T×diag(1,1,2,2,22,…,22︸22,…,2α-1,…,2α-1︸2α-1),
where m=2α and α is a positive integer.The inner product of two Haar functions can be represented as(73)∫01h(m)(t)h(m)T(t)dt=D,
where
(74)D=diag(1,1,1/2,1/2,1/22,…,1/22︸22,…,1/2α-1,…,1/2α-1︸2α-1).


### 4.1.4. Haar Wavelet Solution for Fredholm Integral Equations System [9]

Consider the following Fredholm integral equations system defined in (55):
(75)∑j=1myj(x)=fi(x)+∑j=1m∫01Ki,j(x,t)yj(t)dt,kkkkkkkkkkkkkkkkkkkkkkii=1,2,…,m.The Haar series ofyj(x) and Ki,j(x,t), i=1,2,…,m; j=1,2,…,m are, respectively, expanded as
(76)yj(x)≈YjTh(m)(x),j=1,2,…,m,Ki,j(x,t)≈h(m)T(t)Ki,jh(m)(x),kkkkkkkkkkkii,j=1,2,…,m.
Substituting (76) into (75), we get
(77)∑j=1mYjTh(m)(x)=fi(x)+∑j=1m∫01YjTh(m)(t)h(m)T(t)Ki,jh(m)(x)dt,kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkii=1,2,…,m.
From (77) and (73), we get
(78)∑j=1mYjTh(m)(x)=fi(x)+∑j=1mYjTDKi,jh(m)(x),kkkkkkkkkkkkkkkkkkkkkkkii=1,2,…,m.
Interpolating m collocation points, that is, {xi}i=1m, in the interval [0,1] leads to the following algebraic system of equations:
(79)∑j=1mYjTh(m)(xi)=fi(xi)+∑j=1mYjTDKi,jh(m)(xi),kkkkkkkkkkkkkkkkkkkkkkkkii=1,2,…,m.
Hence, Yj, j=1,2,…,m can be computed by solving the above algebraic system of equations and consequently the solutions yj(x)≈YjTh(m)(x), j=1,2,…,m.


## 4.1.1. Haar Wavelets

The compact set of scale functions is chosen as(56)h0={1,0≤x<1,0,others.
The mother wavelet function is defined as
(57)h1(x)={1,0≤x<12,-1,12≤x<1,0,others.
The family of wavelet functions generated by translation and dilation of h1(x) are given by
(58)hn(x)=h1(2jx-k),
where n=2j+k, j≥0, 0≤k<2j.Mutual orthogonalities of all Haar wavelets can be expressed as(59)∫01hm(x)hn(x)dx=2-jδmn={2-j,m=n=2j+k,0,m≠n.


## 4.1.2. Function Approximation

An arbitrary functiony(x)∈L2[0,1) can be expanded into the following Haar series:
(60)y(x)=∑n=0+∞cnhn(x),
where the coefficients cn are given by
(61)kkkkkkicn=2j∫01y(x)hn(x)dx,n=2j+k,j≥0,0≤k<2j.
In particular, c0=∫01y(x)dx.The previously mentioned expression in (60) can be approximately represented with finite terms as follows:
(62)y(x)≈∑n=0m-1cnhn(x)=C(m)Th(m)(x),
where the coefficient vector C(m)T and the Haar function vector h(m)(x) are, respectively, defined as
(63)C(m)T=[c0,c1,…,cm-1],m=2j,h(m)(x)=[h0(x),h1(x),…,hm-1(x)]T,m=2j.
The Haar expansion for function K(x,t) of order m is defined as follows:
(64)K(x,t)≈∑u=0m-1∑v=0m-1auvhv(x)hu(t),
where auv=2i+q∬01K(x,t)hv(x)hu(t)dxdt, u=2i+j, v=2q+r, i,q≥0.From (62) and (64), we obtain
(65)K(x,t)≈h(m)T(t)Kh(m)(x),
where
(66)K=(auv)m×mT.


## 4.1.3. Operational Matrices of Integration

We define(67)H(m)=[h(m)(12m),h(m)(32m),…,h(m)(2m-12m)],
where H(1)=[1], H(2)=[111-1].Then, form=4, the corresponding matrix can be represented as
(68)H(4)=[h(4)(18),h(4)(38),…,h(4)(78)]=[111111-1-11-100001-1].
The integration of the Haar function vector h(m)(t) is
(69)∫0xh(m)(t)dt=P(m)h(m)(x),kkkkkkkkkkkkkkix∈[0,1),
where P(m) is the operational matrix of order m, and
(70)P(1)=[12],P(m)=12m[2mP(m/2)-H(m/2)H(m/2)-10].
By recursion of the above formula, we obtain
(71)P(2)=14[2-110],P(4)=116[8-4-2-240-2211001-100],P(8)=164[32-16-8-8-4-4-4-4160-88-4-4444400-4400440000-441120000011-2000001-10200001-10-20000].
Therefore, we get
(72)H(m)-1=(1m)H(m)T×diag(1,1,2,2,22,…,22︸22,…,2α-1,…,2α-1︸2α-1),
where m=2α and α is a positive integer.The inner product of two Haar functions can be represented as(73)∫01h(m)(t)h(m)T(t)dt=D,
where
(74)D=diag(1,1,1/2,1/2,1/22,…,1/22︸22,…,1/2α-1,…,1/2α-1︸2α-1).


## 4.1.4. Haar Wavelet Solution for Fredholm Integral Equations System [9]

Consider the following Fredholm integral equations system defined in (55):
(75)∑j=1myj(x)=fi(x)+∑j=1m∫01Ki,j(x,t)yj(t)dt,kkkkkkkkkkkkkkkkkkkkkkii=1,2,…,m.The Haar series ofyj(x) and Ki,j(x,t), i=1,2,…,m; j=1,2,…,m are, respectively, expanded as
(76)yj(x)≈YjTh(m)(x),j=1,2,…,m,Ki,j(x,t)≈h(m)T(t)Ki,jh(m)(x),kkkkkkkkkkkii,j=1,2,…,m.
Substituting (76) into (75), we get
(77)∑j=1mYjTh(m)(x)=fi(x)+∑j=1m∫01YjTh(m)(t)h(m)T(t)Ki,jh(m)(x)dt,kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkii=1,2,…,m.
From (77) and (73), we get
(78)∑j=1mYjTh(m)(x)=fi(x)+∑j=1mYjTDKi,jh(m)(x),kkkkkkkkkkkkkkkkkkkkkkkii=1,2,…,m.
Interpolating m collocation points, that is, {xi}i=1m, in the interval [0,1] leads to the following algebraic system of equations:
(79)∑j=1mYjTh(m)(xi)=fi(xi)+∑j=1mYjTDKi,jh(m)(xi),kkkkkkkkkkkkkkkkkkkkkkkkii=1,2,…,m.
Hence, Yj, j=1,2,…,m can be computed by solving the above algebraic system of equations and consequently the solutions yj(x)≈YjTh(m)(x), j=1,2,…,m.


## 4.2. Taylor Series Expansion Method

In this section, we present Taylor series expansion method for solving Fredholm integral equations system of second kind [7]. This method reduces the system of integral equations to a linear system of ordinary differential equation. After including boundary conditions, this system reduces to a system of equations that can be solved easily by any usual methods.Consider the second kind Fredholm integral equations system defined in (55) as follows:
(80)yi(x)=fi(x)+∑j=1n∫01Ki,j(x,t)yj(t)dt,kkkkkkkkkki=1,2,…,n,0≤x≤1.
A Taylor series expansion can be made for the solution of yj(t) in the integral equation (80):
(81)yj(t)=yj(x)+yj′(x)(t-x)+⋯+1m!yj(m)(x)(t-x)m+E(t),
where E(t) denotes the error between yj(t) and its Taylor series expansion in (81).If we use the firstm term of Taylor series expansion and neglect the term containing E(t), that is, ∫01∑j=1nKi,j(x,t)E(t)dt, then, substituting (81) for yj(t) into the integral in (80), we have
(82)yi(x)≈fi(x)+∑j=1n∫01Ki,j(x,t)∑r=0m1r!(t-x)ryj(r)(x)dt,kkkkkkkkkkkkkkkkkkkkkkii=1,2,…,n,yi(x)≈fi(x)+∑j=1n∑r=0m1r!yj(r)(x)∫01Ki,j(x,t)(t-x)rdt,kkkkkkkkkkkkkkkkkkkkkkii=1,2,…,n,(83)yi(x)-∑j=1n∑r=0m1r!yj(r)(x)[∫01Ki,j(x,t)(t-x)rdt]≈fi(x),i=1,2,…,n.
Equation (83) becomes a linear system of ordinary differential equations that we have to solve. For solving the linear system of ordinary differential equations (83), we require an appropriate number of boundary conditions.In order to construct boundary conditions, we first differentiates times both sides of (80) with respect to x; that is,
(84)yi(s)(x)=fi(s)(x)+∑j=1n∫01Ki,j(s)(x,t)yj(t)dt,kkkkkkkkkkki=1,2,…,n,s=1,2,…,m,
where Ki,j(s)(x,t)=∂(s)Ki,j(x,t)/∂x(s), s=1,2,…,m.Applying the mean value theorem for integral in (84), we have
(85)yi(s)(x)-[∑j=1n∫01Ki,j(s)(x,t)dt]yj(x)≈fi(s)(x),kkkkkkkkkkkkki=1,2,…,n,s=1,2,…,m.
Now (83) combined with (85) becomes a linear system of algebraic equations that can be solved analytically or numerically.


## 4.3. Block-Pulse Functions for the Solution of Fredholm Integral Equation

In this section, Block-Pulse functions (BPF) have been utilized for the solution of system of Fredholm integral equations [6].Anm-set of BPF is defined as follows:
(86)Φi(t)={1,(i-1)Tm≤t<iTm,0,otherwise
with t∈[0,T), T/m=h and i=1,2,…,m.


### 4.3.1. Properties of BPF

(1) Disjointness. One has(87)Φi(t)Φj(t)={Φi(t),i=j;0,i≠j,, i,j=1,2,…,m. This property is obtained from definition of BPF.(2) Orthogonality.
One has(88)∫0TΦi(t)Φj(t)dt={h,i=j;0,i≠j,t∈[0,T), i,j=1,2,…,m. This property is obtained from the disjointness property.(3) Completeness.For every f∈L2,{Φ} is complete; if ∫Φf=0 then f=0 almost everywhere. Because of completeness of {Φ}, we have
(89)∫0Tf2(t)dt=∑i=1∞fi2∥Φi(t)∥2
for every real bounded function f(t) which is square integrable in the interval t∈[0,T) and fi=(1/h)f(t)Φi(t)dt.


### 4.3.2. Function Approximation

The orthogonality property of BPF is the basis of expanding functions into their Block-Pulse series. For everyf(t)∈L2(R),
(90)f(t)=∑i=1mfiΦi(t),
where fi is the coefficient of Block-Pulse function, with respect to ith Block-Pulse function Φi(t).The criterion of this approximation is that mean square error betweenf(t) and its expansion is minimum
(91)ε=1T∫0T(f(t)-∑j=1mfjΦj(t))2dt
so that we can evaluate Block-Pulse coefficients. 										(92)Now∂ε∂fi=-2T∫0T(f(t)-∑j=1mfjΦj(t))Φi(t)dt=0,⟹fi=1h∫0Tf(t)Φi(t)dt(using  orthogonal  property).
In the matrix form, we obtain the following from (90) as follow:
(93)f(t)=∑i=1mfiΦi(t)=FTΦ(t)=ΦTFwhereF=[f1,f2,…,fm]T,Φ(t)=[Φ1(t),Φ2(t),…,Φm(t)]T.
Now let K(t,s) be two-variable function defined on t∈[0,T) and s∈[0,1); then K(t,s) can be expanded to BPF as
(94)K(t,s)=ΦT(t)KΨ(s),
where Φ(t) and Ψ(s) are m1 and m2 dimensional Block-Pulse function vectors and k is a m1×m2 Block-Pulse coefficient matrix.There are two different cases of multiplication of two BPF. The first case is(95)Φ(t)ΦT(t)=(Φ1(t)0⋯00Φ2(t)⋯0⋮⋮⋱⋮00⋯Φm(t)).
It is obtained from disjointness property of BPF. It is a diagonal matrix with m Block-Pulse functions.The second case is(96)ΦT(t)Φ(t)=1
because ∑i=1m(Φi(t))2=∑i=1mΦi(t)=1.Operational Matrix of Integration. BPF integration property can be expressed by an operational equation as
(97)∫0TΦ(t)dt=PΦ(t),
where
(98)Φ(t)=[Φ1(t),Φ2(t),…,Φm(t)]T.
A general formula for Pm×m can be written as
(99)P=12(122⋯2012⋯2001⋯2⋮⋮⋮⋱⋮000⋯1).
By using this matrix, we can express the integral of a function f(t) into its Block-Pulse series
(100)∫0tf(t)dt=∫0tFTΦ(t)dt=FTPΦ(t).


### 4.3.3. Solution for Linear Integral Equations System

Consider the integral equations system from (55) as follows:
(101)∑j=1nyj(x)=fi(x)+∑j=1n∫αβKi,j(x,t)yj(t)dt,kkkkkkkkkkkikki=1,2,…,n.
Block-Pulse coefficients of yj(x), j=1,2,…,n in the interval x∈[α,β) can be determined from the known functions fi(x), i=1,2,…,n and the kernels Ki,j(x,t), i,j=1,2,…n. Usually, we consider α=0 to facilitie the use of Block-Pulse functions. In case α≠0, we set X=((x-α)/(β-α))T, where T=mh.We approximatefi(x), yj(x), Ki,j(x,t) by its BPF as follows:
(102)fi(x)≈FiTΦ(x),yj(x)≈YjTΦ(x),Ki,j(x,t)≈ΦT(t)Ki,jΦ(x),
where Fi, Yj, and Ki,j are defined in Section 4.3.2, and substituting (102) into (101), we have
(103)∑j=1nYjTΦ(x)=FiTΦ(x)+∑j=1n∫0mhYjTΦ(t)ΦT(t)Ki,jΦ(x)dt,kkkkkkkkkkkkkkkkkii=1,2,…,n,(104)∑j=1nYjTΦ(x)=FiTΦ(x)+∑j=1nYjThIKi,jΦ(x),kkkkkkkikkkki=1,2,…,n,
since
(105)∫0mhΦ(t)ΦT(t)dt=hI.
From (104), we get
(106)∑j=1n(I-hKi,jT)Yj=Fi,i=1,2,…,n.
Set Ai,j=I-hKi,jT; then we have from (106)
(107)∑j=1nAi,jYj=Fi,i=1,2,…,n
which is a linear system
(108)(A11A12...A1nA21A22...A2n............An1An2...Ann)(Y1Y2...Yn)=(F1F2...Fn).
After solving the above system we can find Yj, j=1,2,…,n and hence obtain the solutions yj=ΦTYj, j=1,2,…,n.


## 4.3.1. Properties of BPF

(1) Disjointness. One has(87)Φi(t)Φj(t)={Φi(t),i=j;0,i≠j,, i,j=1,2,…,m. This property is obtained from definition of BPF.(2) Orthogonality.
One has(88)∫0TΦi(t)Φj(t)dt={h,i=j;0,i≠j,t∈[0,T), i,j=1,2,…,m. This property is obtained from the disjointness property.(3) Completeness.For every f∈L2,{Φ} is complete; if ∫Φf=0 then f=0 almost everywhere. Because of completeness of {Φ}, we have
(89)∫0Tf2(t)dt=∑i=1∞fi2∥Φi(t)∥2
for every real bounded function f(t) which is square integrable in the interval t∈[0,T) and fi=(1/h)f(t)Φi(t)dt.


## 4.3.2. Function Approximation

The orthogonality property of BPF is the basis of expanding functions into their Block-Pulse series. For everyf(t)∈L2(R),
(90)f(t)=∑i=1mfiΦi(t),
where fi is the coefficient of Block-Pulse function, with respect to ith Block-Pulse function Φi(t).The criterion of this approximation is that mean square error betweenf(t) and its expansion is minimum
(91)ε=1T∫0T(f(t)-∑j=1mfjΦj(t))2dt
so that we can evaluate Block-Pulse coefficients. 										(92)Now∂ε∂fi=-2T∫0T(f(t)-∑j=1mfjΦj(t))Φi(t)dt=0,⟹fi=1h∫0Tf(t)Φi(t)dt(using  orthogonal  property).
In the matrix form, we obtain the following from (90) as follow:
(93)f(t)=∑i=1mfiΦi(t)=FTΦ(t)=ΦTFwhereF=[f1,f2,…,fm]T,Φ(t)=[Φ1(t),Φ2(t),…,Φm(t)]T.
Now let K(t,s) be two-variable function defined on t∈[0,T) and s∈[0,1); then K(t,s) can be expanded to BPF as
(94)K(t,s)=ΦT(t)KΨ(s),
where Φ(t) and Ψ(s) are m1 and m2 dimensional Block-Pulse function vectors and k is a m1×m2 Block-Pulse coefficient matrix.There are two different cases of multiplication of two BPF. The first case is(95)Φ(t)ΦT(t)=(Φ1(t)0⋯00Φ2(t)⋯0⋮⋮⋱⋮00⋯Φm(t)).
It is obtained from disjointness property of BPF. It is a diagonal matrix with m Block-Pulse functions.The second case is(96)ΦT(t)Φ(t)=1
because ∑i=1m(Φi(t))2=∑i=1mΦi(t)=1.Operational Matrix of Integration. BPF integration property can be expressed by an operational equation as
(97)∫0TΦ(t)dt=PΦ(t),
where
(98)Φ(t)=[Φ1(t),Φ2(t),…,Φm(t)]T.
A general formula for Pm×m can be written as
(99)P=12(122⋯2012⋯2001⋯2⋮⋮⋮⋱⋮000⋯1).
By using this matrix, we can express the integral of a function f(t) into its Block-Pulse series
(100)∫0tf(t)dt=∫0tFTΦ(t)dt=FTPΦ(t).


## 4.3.3. Solution for Linear Integral Equations System

Consider the integral equations system from (55) as follows:
(101)∑j=1nyj(x)=fi(x)+∑j=1n∫αβKi,j(x,t)yj(t)dt,kkkkkkkkkkkikki=1,2,…,n.
Block-Pulse coefficients of yj(x), j=1,2,…,n in the interval x∈[α,β) can be determined from the known functions fi(x), i=1,2,…,n and the kernels Ki,j(x,t), i,j=1,2,…n. Usually, we consider α=0 to facilitie the use of Block-Pulse functions. In case α≠0, we set X=((x-α)/(β-α))T, where T=mh.We approximatefi(x), yj(x), Ki,j(x,t) by its BPF as follows:
(102)fi(x)≈FiTΦ(x),yj(x)≈YjTΦ(x),Ki,j(x,t)≈ΦT(t)Ki,jΦ(x),
where Fi, Yj, and Ki,j are defined in Section 4.3.2, and substituting (102) into (101), we have
(103)∑j=1nYjTΦ(x)=FiTΦ(x)+∑j=1n∫0mhYjTΦ(t)ΦT(t)Ki,jΦ(x)dt,kkkkkkkkkkkkkkkkkii=1,2,…,n,(104)∑j=1nYjTΦ(x)=FiTΦ(x)+∑j=1nYjThIKi,jΦ(x),kkkkkkkikkkki=1,2,…,n,
since
(105)∫0mhΦ(t)ΦT(t)dt=hI.
From (104), we get
(106)∑j=1n(I-hKi,jT)Yj=Fi,i=1,2,…,n.
Set Ai,j=I-hKi,jT; then we have from (106)
(107)∑j=1nAi,jYj=Fi,i=1,2,…,n
which is a linear system
(108)(A11A12...A1nA21A22...A2n............An1An2...Ann)(Y1Y2...Yn)=(F1F2...Fn).
After solving the above system we can find Yj, j=1,2,…,n and hence obtain the solutions yj=ΦTYj, j=1,2,…,n.


## 5. Numerical Methods for Nonlinear Fredholm-Hammerstein Integral Equation

We consider the second kind nonlinear Fredholm integral equation of the following form:(109)u(x)=f(x)+∫01K(x,t)F(t,u(t))dt,kkkkkkkkkkkkkkkkkkkki0≤x≤1,
where K(x,t) is the kernel of the integral equation, f(x) and K(x,t) are known functions, and u(x) is the unknown function that is to be determined.


### 5.1.B-Spline Wavelet Method

In this section, nonlinear Fredholm integral equation of second kind of the form given in (109) has been solved by using B-spline wavelets [11].B-spline scaling and wavelet functions in the interval [0,1] and function approximation have been defined in Sections 3.1.1 and 3.1.2, respectively.First, we assume that(110)y(x)=F(x,u(x)),kkkkkk0≤x≤1.
Now, from (16), we can approximate the functions u(x) and y(x) as
(111)u(x)=ATΨ(x),y(x)=BTΨ(x),
where A and B are (2M+1+m-1)×1 column vectors similar to C defined in (17).Again, by using dual of the wavelet functions, we can approximate the functionsf(x) and K(x,t) as follows:
(112)F(x)=DTΨ~(x),K(x,t)=Ψ~T(t)ΘΨ~(x),
where
(113)Θ(i,j)=∫01[∫01K(x,t)Ψi(t)dt]Ψj(x)dx.
From (110)–(112), we get
(114)∫01K(x,t)F(t,u(t))dt=∫01BTΨ(t)Ψ~T(t)ΘΨ~(x)dt=BT[∫01Ψ(t)Ψ~T(t)dt]ΘΨ~(x)=BTΘΨ~(x),since∫01Ψ(t)Ψ~T(t)dt=I.
Applying (110)–(114) in (109), we get
(115)ATΨ(x)=DTΨ~(x)+BTΘΨ~(x).
Multiplying (115) by ΨT(x) both sides from the right and integrating both sides with respect to x from 0 to 1, we have
(116)ATP=DT+BTΘ,ATP-DT-BTΘ=0,
where P is a (2M+1+m-1)×(2M+1+m-1) square matrix given by
(117)P=∫01Ψ(x)ΨT(x)dx=[P1P2],∫01Ψ~(x)ΨT(x)dx=I.
Equation (116) gives a system of (2M+1+m-1) algebraic equations with 2(2M+1+m-1) unknowns for A and B vectors given in (111).To find the solutionu(x) in (111), we first utilize the following equation:
(118)F(x,ATΨ(x))=BTΨ(x),
with the collocation points xi=(i-1)/(2M+1+m-2), where i=1,2,…,2M+1+m-1.Equation (118) gives a system of (2M+1+m-1) algebraic equations with 2(2M+1+m-1) unknowns, for A and B vectors given in (111).Combining (116) and (118), we have a total of 2(2M+1+m-1) system of algebraic equations with 2(2M+1+m-1) unknowns for A and B. Solving those equations for the unknown coefficients in the vectors A and B, we can obtain the solution u(x)=ATΨ(x).


### 5.2. Quadrature Method Applied to Fredholm Integral Equation

In this section, Quadrature method has been applied to solve nonlinear Fredholm-Hammerstein integral equation [10].The quadrature methods like Simpson rule and modified trapezoid method are applied for solving a definite integral as follows.


#### 5.2.1. Simpson’s Rule

One has(119)∫abf(x)dx=∑i=1n-1∫xi-1xi+1f(x)dx=h3f(a)+4h3∑i=1n/2f(x2i-1)+2h3∑i=1(n-1)/2f(x2i)+h3f(b)-(b-a)180h4f(4)(η).


#### 5.2.2. Modified Trapezoid Rule

One has(120)∫abf(x)dx=∑i=1n∫xi-1xif(x)dx=h2f(a)+h∑i=1n-1f(xi)+h2f(b)+h212[f′(a)-f′(b)].
Consider the nonlinear Fredholm integral equation of second kind defined in (109) as follows:
(121)u(x)=f(x)+∫abK(x,t)F(u(t))dt,kkkkkkkkkkkkkkkkkkkia≤x≤b.
For solving (121), we approximate the right-hand integral of (121) with Simpson’s rule and modified trapezoid rule; then we get the following.


#### 5.2.3. Simpson’s Rule

One has(122)u(x)=f(x)+h3[∑j=1n/2K(x,t2j-1)F(u2j-1)K(x,t0)F(u0)kkk+4∑j=1n/2K(x,t2j-1)F(u2j-1)kkk+2∑j=1(n/2)-1K(x,t2j)F(u2j)kkk+K(x,tn)F(un)∑j=1n/2K(x,t2j-1)F(u2j-1)].
Hence, for x=x0,x1,…,xn and t=t0,t1,…,tn in (122), we have
(123)u(xi)=f(xi)+h3[∑j=1n/2K(x,t2j-1)F(u2j-1)K(xi,t0)F(u0)+4∑j=1n/2K(xi,t2j-1)F(u2j-1)+2∑j=1(n/2)-1K(xi,t2j)F(u2j)+K(xi,tn)F(un)∑j=1n/2K(x,t2j-1)F(u2j-1)].
Equation (123) is a nonlinear system of equations and, by solving (123), we obtain the unknowns u(xi) for i=0,1,…,n.


#### 5.2.4. Modified Trapezoid Rule

One has(124)u(x)=f(x)+h2K(x,t0)F(u0)+h∑j=1n-1K(x,tj)F(uj)+h2K(x,tn)F(un)+h212[(x,tn)un′F′(un)J(x,t0)F(u0)+K(x,t0)u0′F′(u0)-J(x,tn)F(un)-K(x,tn)un′F′(un)],
where J(x,t)=∂K(x,t)/∂t.Forx=x0,x1,…,xn and t=t0,t1,…,tn in (124), we have
(125)u(xi)=f(xi)+h2K(xi,t0)F(u0)+h∑j=1n-1K(xi,tj)F(uj)+h2K(xi,tn)F(un)+h212[K(xi,tn)un′F′(un)J(xi,t0)F(u0)+K(xi,t0)u0′F′(u0)-J(xi,tn)F(un)-K(xi,tn)un′F′(un)],
for i=0,1,…,n.This is a system of(n+1) equations and (n+3) unknowns. By taking derivative from (121) and setting H(x,t)=∂K(x,t)/∂x, we obtain
(126)u′(x)=f′(x)+∫abH(x,t)F(u(t))dt,kkkkkkkkkkkkkkkkkkkkka≤x≤b.
If u is a solution of (121), then it is also solution of (126). By using trapezoid rule for (126) and replacing x=xi, we get
(127)u′(xi)=f′(xi)+h2H(xi,t0)F(u0)+h∑j=1n-1H(xi,tj)F(uj)+h2H(xi,tn)F(un),
for i=0,1,…,n. In case of i=0,n from system (127), we obtain two equations.Now (127) combined with (125) generates the nonlinear system of equations as follows:
(128)u(xi)=(h2K(xi,t0)+h212J(xi,t0))F(u0)+h∑j=1n-1K(xi,tj)F(uj)+(h2K(xi,tn)-h212J(xi,tn))F(un)+h212(K(xi,t0)u0′F′(u0)-K(xi,tn)un′F′(un)),u′(x0)=f′(x0)+h2H(x0,t0)F(u0)+h∑j=1n-1H(x0,tj)F(uj)+h2H(x0,tn)F(un),u′(xn)=f′(xn)+h2H(xn,t0)F(u0)+h∑j=1n-1H(xn,tj)F(uj)+h2H(xn,tn)F(un).
By solving this system with (n+3) nonlinear equations and (n+3) unknowns, we can obtain the solution of (109).


### 5.3. Wavelet Galerkin Method

In this section, the continuous Legendre wavelets [12], constructed on the interval [0,1], are applied to solve the nonlinear Fredholm integral equation of the second kind. The nonlinear part of the integral equation is approximated by Legendre wavelets, and the nonlinear integral equation is reduced to a system of nonlinear equations.We have the following family of continuous wavelets with dilation parametera and the translation parameter b(129)ψa,b(t)=|a|-1/2ψ(t-ba),kkkkkkkka,b∈R,a≠0.
Legendre wavelets ψm,n(t)=ψ(k,n^,m,t) have four arguments; k=2,3,…, n^=2n-1, n=1,2,…,2k-1, m is the order for Legendre polynomials and t is the normalized time.Legendre wavelets are defined on [0, 1) by(130)ψm,n(t)={(m+12)1/22k/2Lm(2kt-n^),n^-12k≤t<n^+12k,0,otherwise,
where Lm(t) are the well-known Legendre polynomials of order m, which are orthogonal with respect to the weight function w(t)=1 and satisfy the following recursive formula:
(131)L0(t)=1,L1(t)=t,Lm+1(t)=2m+1m+1tLm(t)-mm+1Lm-1(t),m=1,2,3,….
The set of Legendre wavelets are an orthonormal set.


#### 5.3.1. Function Approximation

A functionf(x)∈L2[0,1] can be expanded as
(132)f(x)=∑n=1∞∑m=0∞cn,mψn,m(x),
where
(133)cn,m=〈f(x),ψn,m(x)〉.
If the infinite series in (132) is truncated, then (132) can be written as
(134)f(x)≈∑n=12k-1∑m=0M-1cn,mψn,m(x)=CTΨ(x),
where C and Ψ(x) are 2k-1M×1 matrices given by
(135)C=[c2k-1,0,…,c2k-1,M-1c1,0,c1,1,…,c1,M-1,c2,0,…,c2,M-1,…,c2k-1,0,…,c2k-1,M-1]T,(136)Ψ(x)=[ψ2k-1,0(x),…,ψ2k-1,M-1(x)ψ1,0(x),…,ψ1,M-1(x),ψ2,0(x),…,ψ2,M-1(x),…,ψ2k-1,0(x),…,ψ2k-1,M-1(x)]T.
Similarly, a function k(x,t)∈L2([0,1]×[0,1]) can be approximated as
(137)k(x,t)≈ΨT(t)KΨ(x),
where K is (2k-1M×2k-1M) matrix, with
(138)Ki,j=〈ψi(t),〈k(x,t),ψj(x)〉〉.
Also, the integer power of a function can be approximated as
(139)[y(x)]p=[YTΨ(x)]p=Yp*TΨ(x),
where Yp* is a column vector, whose elements are nonlinear combinations of the elements of the vector Y. Yp* is called the operational vector of the pth power of the function y(x).


#### 5.3.2. The Operational Matrices

The integration of the vectorΨ(x) defined in (136) can be obtained as
(140)∫0tΨ(t′)dt′=PΨ(t),
where P is the (2k-1M×2k-1M) operational matrix for integration and is given in [23] as
(141)P=[LH⋯HH0L⋯HH⋮⋮⋱⋮⋮00⋯LH00⋯0L].
In (141), H and L are (M×M) matrices given in [23] as(142)H=12k[20⋯000⋯0⋮⋮⋱⋮00⋯0],L=12k[11300⋯00-3303350⋯000-5530557⋯0000-7750⋯00⋮⋮⋮⋮⋱⋮0000⋯02M-3(2M-3)2M-10000⋯-2M-1(2M-1)2M-30].The integration of the product of two Legendre wavelets vector functions is obtained as
(143)∫01Ψ(t)ΨT(t)dt=I,
where I is an identity matrix.The product of two Legendre wavelet vector functions is defined as(144)Ψ(t)ΨT(t)C=C~TΨ(t),
where C is a vector given in (135) and C~ is (2k-1M×2k-1M) matrix, which is called the product operation of Legendre wavelet vector functions [23, 24].


#### 5.3.3. Solution of Fredholm Integral Equation of Second Kind

Consider the nonlinear Fredholm-Hammerstein integral equation of second kind of the form(145)y(x)=f(x)+∫01k(x,t)[y(t)]pdt,
where f∈L2[0,1], k∈L2([0,1]×[0,1]), y is an unknown function, and p is a positive integer.We can approximate the following functions as(146)f(x)≈FTΨ(x),y(x)≈YTΨ(x),k(x,t)≈ΨT(t)KΨ(x),[y(x)]p≈Y*TΨ(x).
Substituting (146) into (145), we have
(147)YTΨ(x)=FTΨ(x)+∫01Y*TΨ(t)ΨT(t)KΨ(x)dt=FTΨ(x)+Y*T(∫01Ψ(t)ΨT(t)dt)KΨ(x)=FTΨ(x)+Y*TKΨ(x)=(FT+Y*TK)Ψ(x)⟹YT-Y*TK-FT=0.
Equation (147) is a system of algebraic equations. Solving (147), we can obtain the solution y(x)≈YTΨ(x).


### 5.4. Homotopy Perturbation Method

Consider the following nonlinear Fredholm integral equation of second kind of the form(148)u(x)=f(x)+∫01K(x,t)F(u(t))dt,kkkkkkkkkkkkkkkkkkki0≤x≤1.
For solving (148) by Homotopy perturbation method (HPM) [14–16], we consider (148) as
(149)L(u)=u(x)-f(x)-∫01K(x,t)F(u(t))dt=0.
As a possible remedy, we can define H(u,p) by
(150)H(u,0)=N(u),H(u,1)=L(u),
where N(u) is an integral operator with known solution u0.We may choose a convex homotopy by(151)H(u,p)=(1-p)N(u)+pL(u)=0
and continuously trace an implicitly defined curve from a starting point H(u0,0) to a solution function H(U,1). The embedding parameter p monotonically increases from zero to unit as the trivial problem L(u)=0. The embedding parameter p∈(0,1] can be considered as an expanding parameter. The HPM uses the homotopy parameter p as an expanding parameter; that is,(152a)u=u0+pu1+p2u2+⋯.
When p→1, (152a) corresponding to (151) become the approximate solution of (149) as follows:
(152b)U=limp→1u=u0+u1+u2+⋯.The series in (152b) converges in most cases, and the rate of convergence depends on L(u) [14].Consider(153)N(u)=u(x)-f(x).
The nonlinear term F(u) can be expressed in He polynomials [25] as
(154)F(u)=∑m=0∞pmHm(u0,u1,…,um)=H0(u0)+pH1(u0,u1)+⋯+pmHm(u0,u1,…um)+⋯,
where
(155)Hm(u0,u1,…,um)=1m!∂m∂pm(F(∑k=0mpkuk))|p=0,m≥0.
Substituting (152a), (153), and (154) into (151), we have
(156)(1-p)((u0+pu1+⋯)-f(x))+p(∑m=0∞pmHm(u0,u1,⋯,um)(u0+pu1+⋯)-f(x)-∫01K(x,t)∑m=0∞pmHm(u0,u1,…,um)dt)=0⟹(u0+pu1+⋯)-f(x)-p∫01K(x,t)∑m=0∞pmHm(u0,u1,…,um)dt=0.
Equating the terms with identical power of p in (156), we have
(157)p0:u0(x)-f(x)=0⟹u0(x)=f(x)p1:u1(x)-∫01K(x,t)H0dt=0⟹u1(x)=∫01K(x,t)H0dtp2:u2(x)-∫01K(x,t)H1dt=0⟹u2(x)=∫01K(x,t)H1dt⋮
and in general form we have
(158)u0(x)=f(x),un+1(x)=∫01K(x,t)Hndt,n=0,1,2,….
Hence, we can obtain the approximate solution of aforesaid equation (148) from (152b).


### 5.5. Adomian Decomposition Method

Adomian decomposition method (ADM) [16–18] has been applied to a wide class of functional equations. This method gives the solution as an infinite series usually converging to an accurate solution. Let us consider the nonlinear Fredholm integral equation of second kind as follows:
(159)u(x)=f(x)+∫abK(x,t)(Lu(t)+Nu(t))dt,kkkkkkkkkkkkkkkkkkkkkkkkkka≤x≤b,
where L(u(t)) and N(u(t)) are the linear and nonlinear terms, respectively.The Adomian decomposition method (ADM) consists of representingu(x) as a series
(160)u(x)=∑m=0∞um(x).
In the view of ADM, the nonlinear term Nu can be represented as
(161)Nu=∑n=0∞An,(162)whereAn=1n!(∂n∂λnN(∑k=0∞λkuk))|λ=0.
Now substituting (160) and (161) into (159), we have
(163)∑m=0∞um(x)=f(x)+∫abK(x,t)(L(∑m=0∞um(t))+∑m=0∞Am)dt,
and, then, ADM uses the recursive relations
(164)u0(x)=f(x),um(x)=∫abK(x,t)(L(um-1(t))+Am-1(t))dt,kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkim≥1,
where Am is so-called Adomian polynomial.Therefore, we obtain then-terms approximate solution as
(165)φn=u0+u1+⋯+un
with
(166)u(x)=limn→∞φn.


## 5.1.B-Spline Wavelet Method

In this section, nonlinear Fredholm integral equation of second kind of the form given in (109) has been solved by using B-spline wavelets [11].B-spline scaling and wavelet functions in the interval [0,1] and function approximation have been defined in Sections 3.1.1 and 3.1.2, respectively.First, we assume that(110)y(x)=F(x,u(x)),kkkkkk0≤x≤1.
Now, from (16), we can approximate the functions u(x) and y(x) as
(111)u(x)=ATΨ(x),y(x)=BTΨ(x),
where A and B are (2M+1+m-1)×1 column vectors similar to C defined in (17).Again, by using dual of the wavelet functions, we can approximate the functionsf(x) and K(x,t) as follows:
(112)F(x)=DTΨ~(x),K(x,t)=Ψ~T(t)ΘΨ~(x),
where
(113)Θ(i,j)=∫01[∫01K(x,t)Ψi(t)dt]Ψj(x)dx.
From (110)–(112), we get
(114)∫01K(x,t)F(t,u(t))dt=∫01BTΨ(t)Ψ~T(t)ΘΨ~(x)dt=BT[∫01Ψ(t)Ψ~T(t)dt]ΘΨ~(x)=BTΘΨ~(x),since∫01Ψ(t)Ψ~T(t)dt=I.
Applying (110)–(114) in (109), we get
(115)ATΨ(x)=DTΨ~(x)+BTΘΨ~(x).
Multiplying (115) by ΨT(x) both sides from the right and integrating both sides with respect to x from 0 to 1, we have
(116)ATP=DT+BTΘ,ATP-DT-BTΘ=0,
where P is a (2M+1+m-1)×(2M+1+m-1) square matrix given by
(117)P=∫01Ψ(x)ΨT(x)dx=[P1P2],∫01Ψ~(x)ΨT(x)dx=I.
Equation (116) gives a system of (2M+1+m-1) algebraic equations with 2(2M+1+m-1) unknowns for A and B vectors given in (111).To find the solutionu(x) in (111), we first utilize the following equation:
(118)F(x,ATΨ(x))=BTΨ(x),
with the collocation points xi=(i-1)/(2M+1+m-2), where i=1,2,…,2M+1+m-1.Equation (118) gives a system of (2M+1+m-1) algebraic equations with 2(2M+1+m-1) unknowns, for A and B vectors given in (111).Combining (116) and (118), we have a total of 2(2M+1+m-1) system of algebraic equations with 2(2M+1+m-1) unknowns for A and B. Solving those equations for the unknown coefficients in the vectors A and B, we can obtain the solution u(x)=ATΨ(x).


## 5.2. Quadrature Method Applied to Fredholm Integral Equation

In this section, Quadrature method has been applied to solve nonlinear Fredholm-Hammerstein integral equation [10].The quadrature methods like Simpson rule and modified trapezoid method are applied for solving a definite integral as follows.


### 5.2.1. Simpson’s Rule

One has(119)∫abf(x)dx=∑i=1n-1∫xi-1xi+1f(x)dx=h3f(a)+4h3∑i=1n/2f(x2i-1)+2h3∑i=1(n-1)/2f(x2i)+h3f(b)-(b-a)180h4f(4)(η).


### 5.2.2. Modified Trapezoid Rule

One has(120)∫abf(x)dx=∑i=1n∫xi-1xif(x)dx=h2f(a)+h∑i=1n-1f(xi)+h2f(b)+h212[f′(a)-f′(b)].
Consider the nonlinear Fredholm integral equation of second kind defined in (109) as follows:
(121)u(x)=f(x)+∫abK(x,t)F(u(t))dt,kkkkkkkkkkkkkkkkkkkia≤x≤b.
For solving (121), we approximate the right-hand integral of (121) with Simpson’s rule and modified trapezoid rule; then we get the following.


### 5.2.3. Simpson’s Rule

One has(122)u(x)=f(x)+h3[∑j=1n/2K(x,t2j-1)F(u2j-1)K(x,t0)F(u0)kkk+4∑j=1n/2K(x,t2j-1)F(u2j-1)kkk+2∑j=1(n/2)-1K(x,t2j)F(u2j)kkk+K(x,tn)F(un)∑j=1n/2K(x,t2j-1)F(u2j-1)].
Hence, for x=x0,x1,…,xn and t=t0,t1,…,tn in (122), we have
(123)u(xi)=f(xi)+h3[∑j=1n/2K(x,t2j-1)F(u2j-1)K(xi,t0)F(u0)+4∑j=1n/2K(xi,t2j-1)F(u2j-1)+2∑j=1(n/2)-1K(xi,t2j)F(u2j)+K(xi,tn)F(un)∑j=1n/2K(x,t2j-1)F(u2j-1)].
Equation (123) is a nonlinear system of equations and, by solving (123), we obtain the unknowns u(xi) for i=0,1,…,n.


### 5.2.4. Modified Trapezoid Rule

One has(124)u(x)=f(x)+h2K(x,t0)F(u0)+h∑j=1n-1K(x,tj)F(uj)+h2K(x,tn)F(un)+h212[(x,tn)un′F′(un)J(x,t0)F(u0)+K(x,t0)u0′F′(u0)-J(x,tn)F(un)-K(x,tn)un′F′(un)],
where J(x,t)=∂K(x,t)/∂t.Forx=x0,x1,…,xn and t=t0,t1,…,tn in (124), we have
(125)u(xi)=f(xi)+h2K(xi,t0)F(u0)+h∑j=1n-1K(xi,tj)F(uj)+h2K(xi,tn)F(un)+h212[K(xi,tn)un′F′(un)J(xi,t0)F(u0)+K(xi,t0)u0′F′(u0)-J(xi,tn)F(un)-K(xi,tn)un′F′(un)],
for i=0,1,…,n.This is a system of(n+1) equations and (n+3) unknowns. By taking derivative from (121) and setting H(x,t)=∂K(x,t)/∂x, we obtain
(126)u′(x)=f′(x)+∫abH(x,t)F(u(t))dt,kkkkkkkkkkkkkkkkkkkkka≤x≤b.
If u is a solution of (121), then it is also solution of (126). By using trapezoid rule for (126) and replacing x=xi, we get
(127)u′(xi)=f′(xi)+h2H(xi,t0)F(u0)+h∑j=1n-1H(xi,tj)F(uj)+h2H(xi,tn)F(un),
for i=0,1,…,n. In case of i=0,n from system (127), we obtain two equations.Now (127) combined with (125) generates the nonlinear system of equations as follows:
(128)u(xi)=(h2K(xi,t0)+h212J(xi,t0))F(u0)+h∑j=1n-1K(xi,tj)F(uj)+(h2K(xi,tn)-h212J(xi,tn))F(un)+h212(K(xi,t0)u0′F′(u0)-K(xi,tn)un′F′(un)),u′(x0)=f′(x0)+h2H(x0,t0)F(u0)+h∑j=1n-1H(x0,tj)F(uj)+h2H(x0,tn)F(un),u′(xn)=f′(xn)+h2H(xn,t0)F(u0)+h∑j=1n-1H(xn,tj)F(uj)+h2H(xn,tn)F(un).
By solving this system with (n+3) nonlinear equations and (n+3) unknowns, we can obtain the solution of (109).


## 5.2.1. Simpson’s Rule

One has(119)∫abf(x)dx=∑i=1n-1∫xi-1xi+1f(x)dx=h3f(a)+4h3∑i=1n/2f(x2i-1)+2h3∑i=1(n-1)/2f(x2i)+h3f(b)-(b-a)180h4f(4)(η).


## 5.2.2. Modified Trapezoid Rule

One has(120)∫abf(x)dx=∑i=1n∫xi-1xif(x)dx=h2f(a)+h∑i=1n-1f(xi)+h2f(b)+h212[f′(a)-f′(b)].
Consider the nonlinear Fredholm integral equation of second kind defined in (109) as follows:
(121)u(x)=f(x)+∫abK(x,t)F(u(t))dt,kkkkkkkkkkkkkkkkkkkia≤x≤b.
For solving (121), we approximate the right-hand integral of (121) with Simpson’s rule and modified trapezoid rule; then we get the following.


## 5.2.3. Simpson’s Rule

One has(122)u(x)=f(x)+h3[∑j=1n/2K(x,t2j-1)F(u2j-1)K(x,t0)F(u0)kkk+4∑j=1n/2K(x,t2j-1)F(u2j-1)kkk+2∑j=1(n/2)-1K(x,t2j)F(u2j)kkk+K(x,tn)F(un)∑j=1n/2K(x,t2j-1)F(u2j-1)].
Hence, for x=x0,x1,…,xn and t=t0,t1,…,tn in (122), we have
(123)u(xi)=f(xi)+h3[∑j=1n/2K(x,t2j-1)F(u2j-1)K(xi,t0)F(u0)+4∑j=1n/2K(xi,t2j-1)F(u2j-1)+2∑j=1(n/2)-1K(xi,t2j)F(u2j)+K(xi,tn)F(un)∑j=1n/2K(x,t2j-1)F(u2j-1)].
Equation (123) is a nonlinear system of equations and, by solving (123), we obtain the unknowns u(xi) for i=0,1,…,n.


## 5.2.4. Modified Trapezoid Rule

One has(124)u(x)=f(x)+h2K(x,t0)F(u0)+h∑j=1n-1K(x,tj)F(uj)+h2K(x,tn)F(un)+h212[(x,tn)un′F′(un)J(x,t0)F(u0)+K(x,t0)u0′F′(u0)-J(x,tn)F(un)-K(x,tn)un′F′(un)],
where J(x,t)=∂K(x,t)/∂t.Forx=x0,x1,…,xn and t=t0,t1,…,tn in (124), we have
(125)u(xi)=f(xi)+h2K(xi,t0)F(u0)+h∑j=1n-1K(xi,tj)F(uj)+h2K(xi,tn)F(un)+h212[K(xi,tn)un′F′(un)J(xi,t0)F(u0)+K(xi,t0)u0′F′(u0)-J(xi,tn)F(un)-K(xi,tn)un′F′(un)],
for i=0,1,…,n.This is a system of(n+1) equations and (n+3) unknowns. By taking derivative from (121) and setting H(x,t)=∂K(x,t)/∂x, we obtain
(126)u′(x)=f′(x)+∫abH(x,t)F(u(t))dt,kkkkkkkkkkkkkkkkkkkkka≤x≤b.
If u is a solution of (121), then it is also solution of (126). By using trapezoid rule for (126) and replacing x=xi, we get
(127)u′(xi)=f′(xi)+h2H(xi,t0)F(u0)+h∑j=1n-1H(xi,tj)F(uj)+h2H(xi,tn)F(un),
for i=0,1,…,n. In case of i=0,n from system (127), we obtain two equations.Now (127) combined with (125) generates the nonlinear system of equations as follows:
(128)u(xi)=(h2K(xi,t0)+h212J(xi,t0))F(u0)+h∑j=1n-1K(xi,tj)F(uj)+(h2K(xi,tn)-h212J(xi,tn))F(un)+h212(K(xi,t0)u0′F′(u0)-K(xi,tn)un′F′(un)),u′(x0)=f′(x0)+h2H(x0,t0)F(u0)+h∑j=1n-1H(x0,tj)F(uj)+h2H(x0,tn)F(un),u′(xn)=f′(xn)+h2H(xn,t0)F(u0)+h∑j=1n-1H(xn,tj)F(uj)+h2H(xn,tn)F(un).
By solving this system with (n+3) nonlinear equations and (n+3) unknowns, we can obtain the solution of (109).


## 5.3. Wavelet Galerkin Method

In this section, the continuous Legendre wavelets [12], constructed on the interval [0,1], are applied to solve the nonlinear Fredholm integral equation of the second kind. The nonlinear part of the integral equation is approximated by Legendre wavelets, and the nonlinear integral equation is reduced to a system of nonlinear equations.We have the following family of continuous wavelets with dilation parametera and the translation parameter b(129)ψa,b(t)=|a|-1/2ψ(t-ba),kkkkkkkka,b∈R,a≠0.
Legendre wavelets ψm,n(t)=ψ(k,n^,m,t) have four arguments; k=2,3,…, n^=2n-1, n=1,2,…,2k-1, m is the order for Legendre polynomials and t is the normalized time.Legendre wavelets are defined on [0, 1) by(130)ψm,n(t)={(m+12)1/22k/2Lm(2kt-n^),n^-12k≤t<n^+12k,0,otherwise,
where Lm(t) are the well-known Legendre polynomials of order m, which are orthogonal with respect to the weight function w(t)=1 and satisfy the following recursive formula:
(131)L0(t)=1,L1(t)=t,Lm+1(t)=2m+1m+1tLm(t)-mm+1Lm-1(t),m=1,2,3,….
The set of Legendre wavelets are an orthonormal set.


### 5.3.1. Function Approximation

A functionf(x)∈L2[0,1] can be expanded as
(132)f(x)=∑n=1∞∑m=0∞cn,mψn,m(x),
where
(133)cn,m=〈f(x),ψn,m(x)〉.
If the infinite series in (132) is truncated, then (132) can be written as
(134)f(x)≈∑n=12k-1∑m=0M-1cn,mψn,m(x)=CTΨ(x),
where C and Ψ(x) are 2k-1M×1 matrices given by
(135)C=[c2k-1,0,…,c2k-1,M-1c1,0,c1,1,…,c1,M-1,c2,0,…,c2,M-1,…,c2k-1,0,…,c2k-1,M-1]T,(136)Ψ(x)=[ψ2k-1,0(x),…,ψ2k-1,M-1(x)ψ1,0(x),…,ψ1,M-1(x),ψ2,0(x),…,ψ2,M-1(x),…,ψ2k-1,0(x),…,ψ2k-1,M-1(x)]T.
Similarly, a function k(x,t)∈L2([0,1]×[0,1]) can be approximated as
(137)k(x,t)≈ΨT(t)KΨ(x),
where K is (2k-1M×2k-1M) matrix, with
(138)Ki,j=〈ψi(t),〈k(x,t),ψj(x)〉〉.
Also, the integer power of a function can be approximated as
(139)[y(x)]p=[YTΨ(x)]p=Yp*TΨ(x),
where Yp* is a column vector, whose elements are nonlinear combinations of the elements of the vector Y. Yp* is called the operational vector of the pth power of the function y(x).


### 5.3.2. The Operational Matrices

The integration of the vectorΨ(x) defined in (136) can be obtained as
(140)∫0tΨ(t′)dt′=PΨ(t),
where P is the (2k-1M×2k-1M) operational matrix for integration and is given in [23] as
(141)P=[LH⋯HH0L⋯HH⋮⋮⋱⋮⋮00⋯LH00⋯0L].
In (141), H and L are (M×M) matrices given in [23] as(142)H=12k[20⋯000⋯0⋮⋮⋱⋮00⋯0],L=12k[11300⋯00-3303350⋯000-5530557⋯0000-7750⋯00⋮⋮⋮⋮⋱⋮0000⋯02M-3(2M-3)2M-10000⋯-2M-1(2M-1)2M-30].The integration of the product of two Legendre wavelets vector functions is obtained as
(143)∫01Ψ(t)ΨT(t)dt=I,
where I is an identity matrix.The product of two Legendre wavelet vector functions is defined as(144)Ψ(t)ΨT(t)C=C~TΨ(t),
where C is a vector given in (135) and C~ is (2k-1M×2k-1M) matrix, which is called the product operation of Legendre wavelet vector functions [23, 24].


### 5.3.3. Solution of Fredholm Integral Equation of Second Kind

Consider the nonlinear Fredholm-Hammerstein integral equation of second kind of the form(145)y(x)=f(x)+∫01k(x,t)[y(t)]pdt,
where f∈L2[0,1], k∈L2([0,1]×[0,1]), y is an unknown function, and p is a positive integer.We can approximate the following functions as(146)f(x)≈FTΨ(x),y(x)≈YTΨ(x),k(x,t)≈ΨT(t)KΨ(x),[y(x)]p≈Y*TΨ(x).
Substituting (146) into (145), we have
(147)YTΨ(x)=FTΨ(x)+∫01Y*TΨ(t)ΨT(t)KΨ(x)dt=FTΨ(x)+Y*T(∫01Ψ(t)ΨT(t)dt)KΨ(x)=FTΨ(x)+Y*TKΨ(x)=(FT+Y*TK)Ψ(x)⟹YT-Y*TK-FT=0.
Equation (147) is a system of algebraic equations. Solving (147), we can obtain the solution y(x)≈YTΨ(x).


## 5.3.1. Function Approximation

A functionf(x)∈L2[0,1] can be expanded as
(132)f(x)=∑n=1∞∑m=0∞cn,mψn,m(x),
where
(133)cn,m=〈f(x),ψn,m(x)〉.
If the infinite series in (132) is truncated, then (132) can be written as
(134)f(x)≈∑n=12k-1∑m=0M-1cn,mψn,m(x)=CTΨ(x),
where C and Ψ(x) are 2k-1M×1 matrices given by
(135)C=[c2k-1,0,…,c2k-1,M-1c1,0,c1,1,…,c1,M-1,c2,0,…,c2,M-1,…,c2k-1,0,…,c2k-1,M-1]T,(136)Ψ(x)=[ψ2k-1,0(x),…,ψ2k-1,M-1(x)ψ1,0(x),…,ψ1,M-1(x),ψ2,0(x),…,ψ2,M-1(x),…,ψ2k-1,0(x),…,ψ2k-1,M-1(x)]T.
Similarly, a function k(x,t)∈L2([0,1]×[0,1]) can be approximated as
(137)k(x,t)≈ΨT(t)KΨ(x),
where K is (2k-1M×2k-1M) matrix, with
(138)Ki,j=〈ψi(t),〈k(x,t),ψj(x)〉〉.
Also, the integer power of a function can be approximated as
(139)[y(x)]p=[YTΨ(x)]p=Yp*TΨ(x),
where Yp* is a column vector, whose elements are nonlinear combinations of the elements of the vector Y. Yp* is called the operational vector of the pth power of the function y(x).


## 5.3.2. The Operational Matrices

The integration of the vectorΨ(x) defined in (136) can be obtained as
(140)∫0tΨ(t′)dt′=PΨ(t),
where P is the (2k-1M×2k-1M) operational matrix for integration and is given in [23] as
(141)P=[LH⋯HH0L⋯HH⋮⋮⋱⋮⋮00⋯LH00⋯0L].
In (141), H and L are (M×M) matrices given in [23] as(142)H=12k[20⋯000⋯0⋮⋮⋱⋮00⋯0],L=12k[11300⋯00-3303350⋯000-5530557⋯0000-7750⋯00⋮⋮⋮⋮⋱⋮0000⋯02M-3(2M-3)2M-10000⋯-2M-1(2M-1)2M-30].The integration of the product of two Legendre wavelets vector functions is obtained as
(143)∫01Ψ(t)ΨT(t)dt=I,
where I is an identity matrix.The product of two Legendre wavelet vector functions is defined as(144)Ψ(t)ΨT(t)C=C~TΨ(t),
where C is a vector given in (135) and C~ is (2k-1M×2k-1M) matrix, which is called the product operation of Legendre wavelet vector functions [23, 24].


## 5.3.3. Solution of Fredholm Integral Equation of Second Kind

Consider the nonlinear Fredholm-Hammerstein integral equation of second kind of the form(145)y(x)=f(x)+∫01k(x,t)[y(t)]pdt,
where f∈L2[0,1], k∈L2([0,1]×[0,1]), y is an unknown function, and p is a positive integer.We can approximate the following functions as(146)f(x)≈FTΨ(x),y(x)≈YTΨ(x),k(x,t)≈ΨT(t)KΨ(x),[y(x)]p≈Y*TΨ(x).
Substituting (146) into (145), we have
(147)YTΨ(x)=FTΨ(x)+∫01Y*TΨ(t)ΨT(t)KΨ(x)dt=FTΨ(x)+Y*T(∫01Ψ(t)ΨT(t)dt)KΨ(x)=FTΨ(x)+Y*TKΨ(x)=(FT+Y*TK)Ψ(x)⟹YT-Y*TK-FT=0.
Equation (147) is a system of algebraic equations. Solving (147), we can obtain the solution y(x)≈YTΨ(x).


## 5.4. Homotopy Perturbation Method

Consider the following nonlinear Fredholm integral equation of second kind of the form(148)u(x)=f(x)+∫01K(x,t)F(u(t))dt,kkkkkkkkkkkkkkkkkkki0≤x≤1.
For solving (148) by Homotopy perturbation method (HPM) [14–16], we consider (148) as
(149)L(u)=u(x)-f(x)-∫01K(x,t)F(u(t))dt=0.
As a possible remedy, we can define H(u,p) by
(150)H(u,0)=N(u),H(u,1)=L(u),
where N(u) is an integral operator with known solution u0.We may choose a convex homotopy by(151)H(u,p)=(1-p)N(u)+pL(u)=0
and continuously trace an implicitly defined curve from a starting point H(u0,0) to a solution function H(U,1). The embedding parameter p monotonically increases from zero to unit as the trivial problem L(u)=0. The embedding parameter p∈(0,1] can be considered as an expanding parameter. The HPM uses the homotopy parameter p as an expanding parameter; that is,(152a)u=u0+pu1+p2u2+⋯.
When p→1, (152a) corresponding to (151) become the approximate solution of (149) as follows:
(152b)U=limp→1u=u0+u1+u2+⋯.The series in (152b) converges in most cases, and the rate of convergence depends on L(u) [14].Consider(153)N(u)=u(x)-f(x).
The nonlinear term F(u) can be expressed in He polynomials [25] as
(154)F(u)=∑m=0∞pmHm(u0,u1,…,um)=H0(u0)+pH1(u0,u1)+⋯+pmHm(u0,u1,…um)+⋯,
where
(155)Hm(u0,u1,…,um)=1m!∂m∂pm(F(∑k=0mpkuk))|p=0,m≥0.
Substituting (152a), (153), and (154) into (151), we have
(156)(1-p)((u0+pu1+⋯)-f(x))+p(∑m=0∞pmHm(u0,u1,⋯,um)(u0+pu1+⋯)-f(x)-∫01K(x,t)∑m=0∞pmHm(u0,u1,…,um)dt)=0⟹(u0+pu1+⋯)-f(x)-p∫01K(x,t)∑m=0∞pmHm(u0,u1,…,um)dt=0.
Equating the terms with identical power of p in (156), we have
(157)p0:u0(x)-f(x)=0⟹u0(x)=f(x)p1:u1(x)-∫01K(x,t)H0dt=0⟹u1(x)=∫01K(x,t)H0dtp2:u2(x)-∫01K(x,t)H1dt=0⟹u2(x)=∫01K(x,t)H1dt⋮
and in general form we have
(158)u0(x)=f(x),un+1(x)=∫01K(x,t)Hndt,n=0,1,2,….
Hence, we can obtain the approximate solution of aforesaid equation (148) from (152b).


## 5.5. Adomian Decomposition Method

Adomian decomposition method (ADM) [16–18] has been applied to a wide class of functional equations. This method gives the solution as an infinite series usually converging to an accurate solution. Let us consider the nonlinear Fredholm integral equation of second kind as follows:
(159)u(x)=f(x)+∫abK(x,t)(Lu(t)+Nu(t))dt,kkkkkkkkkkkkkkkkkkkkkkkkkka≤x≤b,
where L(u(t)) and N(u(t)) are the linear and nonlinear terms, respectively.The Adomian decomposition method (ADM) consists of representingu(x) as a series
(160)u(x)=∑m=0∞um(x).
In the view of ADM, the nonlinear term Nu can be represented as
(161)Nu=∑n=0∞An,(162)whereAn=1n!(∂n∂λnN(∑k=0∞λkuk))|λ=0.
Now substituting (160) and (161) into (159), we have
(163)∑m=0∞um(x)=f(x)+∫abK(x,t)(L(∑m=0∞um(t))+∑m=0∞Am)dt,
and, then, ADM uses the recursive relations
(164)u0(x)=f(x),um(x)=∫abK(x,t)(L(um-1(t))+Am-1(t))dt,kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkim≥1,
where Am is so-called Adomian polynomial.Therefore, we obtain then-terms approximate solution as
(165)φn=u0+u1+⋯+un
with
(166)u(x)=limn→∞φn.


## 6. Conclusion and Discussion

In this work, we have examined many numerical methods to solve Fredholm integral equations. Using these methods except variational iteration method, the Fredholm integral equations have been reduced to a system of algebraic equations and this system can be easily solved by any usual methods. In this work, we have applied compactly supported semiorthogonalB-spline wavelets along with their dual wavelets for solving both linear and nonlinear Fredholm integral equations of second kind. The problem has been reduced to solve a system of algebraic equations. In order to increase the accuracy of the approximate solution, it is necessary to apply higher-order B-spline wavelet method. The method of moments based on compactly supported semiorthogonal B-spline wavelets via Galerkin method has been used to solve Fredholm integral equation of second kind. This method determines a strong reduction in the computation time and memory requirement in inverting the matrix. Variational iteration method has been successfully applied to find the approximate solution of Fredholm integral equation of both linear and nonlinear types. Taylor series expansion method reduces the system of integral equations to a linear system of ordinary differential equation. After including the required boundary conditions, this system reduces to a system of algebraic equations that can be solved easily. Block-Pulse functions and Haar wavelet method can be applied to the system of Fredholm integral equations by reducing into a system of algebraic equations. These methods give more accuracy if we increase their order. Quadrature method can be applied to solve the nonlinear Fredholm-Hammerstein integral equation of second kind by reducing it to a system of algebraic equations. Homotopy perturbation method (HPM) and Adomian decomposition method (ADM) can be also applied to approximate the solution of nonlinear Fredholm integral equation of second kind. The solutions obtained by HPM and ADM are applicable for not only weakly nonlinear equations, but also strong ones. The approximate solutions by these aforesaid methods highly agree with exact solutions.



---

*Source: 426916-2013-12-18.xml*