Numerical Methods for Solving Fredholm Integral Equations of Second Kind

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.

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. The B-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) where 𝑎 and 𝑏 are both constants. 𝑓(𝑥), 𝑔(𝑥), and 𝐾(𝑥,𝑡) are known functions while 𝑦(𝑥) is unknown function. 𝜆 (nonzero parameter) is called eigenvalue of the integral equation. The function 𝐾(𝑥,𝑡) is known as kernel of the integral equation.

2.1. Fredholm Integral Equation of First Kind

The linear integral equation is of form (by setting 𝑔(𝑥)=0 in (1)) 𝑓(𝑥)+𝜆^𝑏_𝑎𝐾(𝑥,𝑡)𝑦(𝑡)𝑑𝑡=0. (2) 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 setting 𝑔(𝑥)=1 in (1)) 𝑦(𝑥)=𝑓(𝑥)+^𝑏_𝑎𝐾(𝑥,𝑡)𝑦(𝑡)𝑑𝑡. (3) 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: ^𝑛_𝑗=1𝑔_𝑖,𝑗𝑦_𝑗(𝑥)=𝑓_𝑖(𝑥)+^𝑛_𝑗=1^𝑏_𝑎𝐾_𝑖,𝑗(𝑥,𝑡)𝑦_𝑗(𝑡)𝑑𝑡,𝑖=1,2,…,𝑛, (4) where 𝑓_𝑖(𝑥) and 𝐾_𝑖,𝑗(𝑥,𝑡) are known functions and 𝑦_𝑗(𝑥) are the unknown functions for 𝑖,𝑗=1,2,…,𝑛.

2.4. Nonlinear Fredholm-Hammerstein Integral Equation of Second Kind

Nonlinear Fredholm-Hammerstein integral equation of second kind is defined as follows: 𝑦(𝑥)=𝑓(𝑥)+^𝑏_𝑎𝐾(𝑥,𝑡)𝐹𝑦(𝑡)𝑑𝑡, (5) where 𝐾(𝑥,𝑡) is the kernel of the integral equation, 𝑓(𝑥) and 𝐾(𝑥,𝑡) are known functions, and 𝑦(𝑥) 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: ^𝑛_𝑗=1𝑔_𝑖,𝑗𝑦_𝑗(𝑥)=𝑓_𝑖(𝑥)+^𝑛_𝑗=1^𝑏_𝑎𝐾_𝑖,𝑗(𝑥,𝑡)𝐹_𝑖,𝑗𝑡,𝑦_𝑗(𝑡)𝑑𝑡,𝑖=1,2,…,𝑛, (6) where 𝑓_𝑖(𝑥) and 𝐾_𝑖,𝑗(𝑥,𝑡) are known functions and 𝑦_𝑗(𝑥) are the unknown functions for 𝑖,𝑗=1,2,…,𝑛.

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

Consider the following Fredholm integral equation of second kind defined in (3) 𝑦(𝑥)=𝑓(𝑥)+^𝑏_𝑎𝐾(𝑥,𝑡)𝑦(𝑡)𝑑𝑡,𝑎≤𝑥≤𝑏, (7) where 𝐾(𝑥,𝑡) and 𝑔(𝑥) are known functions and 𝑦(𝑥) 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 using B-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 𝐿²(𝑅).

Using this basis, an arbitrary function in 𝐿²(𝑅) 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.

Let 𝑚 and 𝑛 be two positive integers and let 𝑎=𝑥_−𝑚+1=⋯=𝑥₀<𝑥₁<⋯<𝑥_𝑛=𝑥_𝑛+1=⋯=𝑥_{𝑛+𝑚−1}=𝑏 (8) be an equally spaced knots sequence. The functions 𝐵_{𝑚,𝑗,𝑋}(𝑥)=𝑥−𝑥_𝑗𝑥_{𝑗+𝑚−1}−𝑥_𝑗𝐵_{𝑚−1,𝑗,𝑋}(𝑥)+𝑥_𝑗+𝑚−𝑥𝑥_𝑗+𝑚−𝑥_𝑗+1𝐵_{𝑚−1,𝑗+1,𝑋}(𝑥),𝐵𝑗=−𝑚+1,…,𝑛−1,_1,𝑗,𝑋(𝑥)=1,𝑥∈𝑥_𝑗,𝑥_𝑗+1,0,otherwise, (9) are called cardinal B-spline functions of order 𝑚≥2 for the knot sequence 𝑋={𝑥_𝑖}^{𝑛+𝑚−1}_{𝑖=−𝑚+1} and Supp𝐵_{𝑚,𝑗,𝑋}(𝑥)=[𝑥_𝑗,𝑥_𝑗+𝑚]∩[𝑎,𝑏].

By considering the interval [𝑎,𝑏]=[0,1], at any level 𝑗∈Ζ⁺, the discretization step is 2^−𝑗, and this generates 𝑛=2^𝑗 number of segments in [0,1] with knot sequence 𝑋^(𝑗)=𝑥^(𝑗)_−𝑚+1=⋯=𝑥₀^(𝑗)𝑥=0,_𝑘(𝑗)=𝑘2^𝑗𝑥,𝑘=1,…,𝑛−1,_𝑛(𝑗)=⋯=𝑥^(𝑗)_{𝑛+𝑚−1}=1. (10) Let 𝑗₀ be the level for which 2^𝑗0≥2𝑚−1; for each level, 𝑗≥𝑗₀, the scaling functions of order 𝑚 can be defined as follows in [2]: 𝜑_{𝑚,𝑗,𝑖}(𝑥)=𝐵_{𝑚,𝑗0,𝑖}2^𝑗−𝑗0𝑥𝐵𝑖=−𝑚+1,…,−1,_{𝑚,𝑗0,2𝑗−𝑚−𝑖}1−2^𝑗−𝑗0𝑥𝑖=2^𝑗−𝑚+1,…,2^𝑗𝐵−1,_𝑚,𝑗0,02^𝑗−𝑗0𝑥−2^−𝑗0𝑖𝑖=0,…,2^𝑗−𝑚. (11) And the two scale relations for the 𝑚-order semiorthogonal compactly supported B-wavelet functions are defined as follows: 𝜓_{𝑚,𝑗,𝑖−𝑚}=^{2𝑖+2𝑚−2}_𝑘=𝑖𝑞_𝑖,𝑘𝐵_{𝑚,𝑗,𝑘−𝑚}𝜓,𝑖=1,…,𝑚−1,_{𝑚,𝑗,𝑖−𝑚}=^{2𝑖+2𝑚−2}_{𝑘=2𝑖−𝑚}𝑞_𝑖,𝑘𝐵_{𝑚,𝑗,𝑘−𝑚}𝜓,𝑖=𝑚,…,𝑛−𝑚+1,_{𝑚,𝑗,𝑖−𝑚}=^{𝑛+𝑖+𝑚−1}_{𝑘=2𝑖−𝑚}𝑞_𝑖,𝑘𝐵_{𝑚,𝑗,𝑘−𝑚},𝑖=𝑛−𝑚+2,…,𝑛, (12) where 𝑞_𝑖,𝑘=𝑞_𝑘−2𝑖.

Hence, there are 2(𝑚−1) boundary wavelets and (𝑛−2𝑚+2) inner wavelets in the bounded interval [𝑎,𝑏]. Finally, by considering the level 𝑗 with 𝑗≥𝑗₀, the B-wavelet functions in [0,1] can be expressed as follows: 𝜓_{𝑚,𝑗,𝑖}(𝑥)=𝜓_{𝑚,𝑗0,𝑖}2^𝑗−𝑗0𝑥𝜓𝑖=−𝑚+1,…,−1,_{𝑚,2𝑗−2𝑚+1−𝑖,𝑖}1−2^𝑗−𝑗0𝑥𝑖=2^𝑗−2𝑚+2,…,2^𝑗𝜓−𝑚,_𝑚,𝑗0,02^𝑗−𝑗0𝑥−2^−𝑗0𝑖𝑖=0,…,2^𝑗−2𝑚+1. (13) The scaling functions 𝜑_{𝑚,𝑗,𝑖}(𝑥) occupy 𝑚 segments and the wavelet functions 𝜓_{𝑚,𝑗,𝑖}(𝑥) occupy 2𝑚−1 segments.

When the semiorthogonal wavelets are constructed from B-spline of order 𝑚, the lowest octave level 𝑗=𝑗₀ is determined in [19, 20] by 2^𝑗0≥2𝑚−1, (14) so as to have a minimum of one complete wavelet on the interval [0,1].

3.1.2. Function Approximation

A function 𝑓(𝑥) defined over [0,1] may be approximated by B-spline wavelets as [21, 22] 𝑓(𝑥)=^2𝑗0−1_{𝑘=1−𝑚}𝑐_𝑗0,𝑘𝜑_𝑗0,𝑘(𝑥)+^∞_𝑗=𝑗0^2𝑗−𝑚_{𝑘=1−𝑚}𝑑_𝑗,𝑘𝜓_𝑗,𝑘(𝑥). (15) If the infinite series in (15) is truncated at 𝑀, then (15) can be written as [2] 𝑓(𝑥)≅^2𝑗0−1_{𝑘=1−𝑚}𝑐_𝑗0,𝑘𝜑_𝑗0,𝑘(𝑥)+^𝑀_𝑗=𝑗0^2𝑗−𝑚_{𝑘=1−𝑚}𝑑_𝑗,𝑘𝜓_𝑗,𝑘(𝑥), (16) where 𝜑_2,𝑘 and 𝜓_𝑗,𝑘 are scaling and wavelets functions, respectively, and 𝐶 and Ψ are (2^𝑀+1+𝑚−1)×1 vectors given by 𝐶=𝑐_{𝑗0,1−𝑚},…,𝑐_𝑗0,2^𝑗₀₋₁,𝑑_{𝑗0,1−𝑚}𝑑,…,_{𝑗0,2𝑗0−𝑚},…,𝑑_{𝑀,1−𝑚},…,𝑑_𝑀,2^𝑀_−𝑚^𝑇,Ψ=𝜑_{𝑗0,1−𝑚},…,𝜑_{𝑗0,2𝑗0−1},𝜓_{𝑗0,1−𝑚}𝜓,…,_{𝑗0,2𝑗0−𝑚},…,𝜓_{𝑀,1−𝑚},…,𝜓_{𝑀,2𝑀−𝑚}^𝑇, (17)(18) with 𝑐_𝑗0,𝑘=¹₀𝑓(𝑥)𝜑_𝑗0,𝑘(𝑥)𝑑𝑥,𝑘=1−𝑚,…,2^𝑗0𝑑−1,_𝑗,𝑘=¹₀𝑓(𝑥)𝜓_𝑗,𝑘(𝑥)𝑑𝑥,𝑗=𝑗₀,…,𝑀,𝑘=1−𝑚,…,2^𝑀−𝑚, (19) where 𝜑_𝑗0,𝑘(𝑥) and 𝜓_𝑗,𝑘(𝑥) are dual functions of 𝜑_𝑗0,𝑘 and 𝜓_𝑗,𝑘, respectively. These can be obtained by linear combinations of 𝜑_𝑗0,𝑘, 𝑘=1−𝑚,…,2^𝑗0−1, and 𝜓_𝑗,𝑘, 𝑗=𝑗₀,…,𝑀, 𝑘=1−𝑚,…,2^𝑀−𝑚, as follows. Let Φ=𝜑_{𝑗0,1−𝑚},…,𝜑_{𝑗0,2𝑗0−1}^𝑇,Ψ=𝜓_{𝑗0,1−𝑚},…,𝜓_{𝑗0,2𝑗0−𝑚},…,𝜓_{𝑀,1−𝑚},…,𝜓_{𝑀,2𝑀−𝑚}^𝑇. (20)(21) Using (11), (20), (12)-(13), and (21), we get ¹₀ΦΦ^𝑇𝑑𝑥=𝑃₁,₁₀ΨΨ^𝑇𝑑𝑥=𝑃₂. (22) Suppose that Φ and Ψ are the dual functions of Φ and Ψ, respectively; then ¹₀ΦΦ^𝑇𝑑𝑥=𝐼₁,₁₀ΨΨ^𝑇𝑑𝑥=𝐼₂,Φ=𝑃₁⁻¹Φ,Ψ=𝑃₂⁻¹Ψ. (23)(24)

3.1.3. Application of B-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 𝑦(𝑥) as 𝑦(𝑥)=𝐶^𝑇Ψ(𝑥), (25) where Ψ(𝑥) is defined in (18) and 𝐶 is (2^𝑀+1+𝑚−1)×1 unknown vector defined similarly as in (17). We also expand 𝑓(𝑥) and 𝐾(𝑥,𝑡) by B-spline dual wavelets Ψ defined in (24) as 𝑓(𝑥)=𝐶₁^𝑇Ψ(𝑥),𝐾(𝑥,𝑡)=Ψ^𝑇(𝑡)ΘΨ(𝑥), (26) where Θ_𝑖,𝑗=¹₀¹₀𝐾(𝑥,𝑡)Ψ_𝑖(𝑡)𝑑𝑡Ψ_𝑗(𝑥)𝑑𝑥. (27) From (26) and (25), we get ¹₀𝐾(𝑥,𝑡)𝑦(𝑡)𝑑𝑡=¹₀𝐶^𝑇Ψ(𝑡)Ψ^𝑇(𝑡)ΘΨ(𝑥)𝑑𝑡=𝐶^𝑇ΘΨ(𝑥) (28) since ¹₀Ψ(𝑡)Ψ^𝑇(𝑡)𝑑𝑡=𝐼. (29) By applying (25)–(28) in (7) we have 𝐶^𝑇Ψ(𝑥)−𝐶^𝑇ΘΨ(𝑥)=𝐶^𝑇₁Ψ(𝑥). (30) By multiplying both sides of (30) with Ψ^𝑇(𝑥) from the right and integrating both sides with respect to 𝑥 from 0 to 1, we get 𝐶^𝑇𝑃−𝐶^𝑇Θ=𝐶^𝑇₁, (31) since ¹₀Ψ(𝑥)Ψ^𝑇(𝑥)𝑑𝑥=𝐼, (32) and 𝑃 is a (2^𝑀+1+𝑚−1)×(2^𝑀+1+𝑚−1) square matrix given by 𝑃=¹₀Ψ(𝑥)Ψ^𝑇(𝑥)𝑃𝑑𝑥=¹00𝑃₂. (33) Consequently, from (31), we get 𝐶^𝑇=𝐶^𝑇₁(𝑃−Θ)⁻¹. Hence, we can calculate the solution for 𝑦(𝑥)=𝐶^𝑇Ψ(𝑥).

3.2. Method of Moments

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

A set of subspaces {𝑉_𝑗}_𝑗∈𝑍 is said to be MRA of 𝐿²(𝑅) if it possesses the following properties: 𝑉_𝑗⊂𝑉_𝑗+1,∀𝑗∈Ζ,_𝑗∈Ζ𝑉_𝑗isdensein𝐿²(𝑅),_𝑗∈Ζ𝑉_𝑗𝑓(𝑥)=𝜙,∈𝑉_𝑗⟺𝑓(2𝑥)∈𝑉_𝑗+1,∀𝑗∈Ζ, (34)(35)(36)(37) where 𝑍 denotes the set of integers. Properties (34)–(36) state that {𝑉_𝑗}_𝑗∈𝑍 is a nested sequence of subspaces that effectively covers 𝐿²(𝑅). 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 𝑉_𝑗. A function 𝜑∈𝐿²(𝑅) is called a scaling function if it generates the nested sequence of subspaces 𝑉_𝑗 and satisfies the dilation equation; namely, 𝜑(𝑥)=_𝑘𝑝_𝑘𝜑(𝑎𝑥−𝑘), (38) with 𝑝_𝑘∈𝑙² and 𝑎 being any rational number.

For each scale 𝑗, since 𝑉_𝑗⊂𝑉_𝑗+1, there exists a unique orthogonal complementary subspace 𝑊_𝑗 of 𝑉_𝑗 in 𝑉_𝑗+1. This subspace 𝑊_𝑗 is called wavelet subspace and is generated by 𝜓_𝑗,𝑘=𝜓(2^𝑗𝑥−𝑘), where 𝜓∈𝐿² is called the wavelet. From the above discussion, these results follow easily: 𝑉_𝑗1∩𝑉_𝑗2=𝑉_𝑗2,𝑗₁>𝑗₂,𝑊_𝑗1∩𝑊_𝑗2=0,𝑗₁≠𝑗₂,𝑉_𝑗1∩𝑊_𝑗2=0,𝑗₁≤𝑗₂. (39) 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 𝑚 if ^∞_−∞𝑥^𝑝𝜓(𝑥)𝑑𝑥=0;𝑝=0,…,𝑚−1. (40) All wavelets must satisfy the previously mentioned condition for 𝑝=0.

(2) Semiorthogonality. The wavelets 𝜓_𝑗,𝑘 form a semiorthogonal basis if 𝜓_𝑗,𝑘,𝜓_𝑠,𝑖=0;𝑗≠𝑠;∀𝑗,𝑘,𝑠,𝑖∈Ζ. (41)

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: 𝑦(𝑥)≈^2𝑗0−1_𝑘=−1𝑐_𝑘𝜑_𝑗0,𝑘(𝑥)+^𝑀_𝑗=𝑗0^2𝑗−2_𝑘=−1𝑑_𝑗,𝑘𝜓_𝑗,𝑘(𝑥)=𝐶^𝑇Ψ(𝑥). (42) By substituting this expression into (7) and employing the Galerkin method, the following set of linear system of order (2^𝑀+1) is generated. The scaling and wavelet functions are used as testing and weighting functions: −−−−𝐶𝐷=𝜑,𝜑𝐾𝜑,𝜑𝜓,𝜑𝐾𝜓,𝜑𝐹𝜑,𝜓𝐾𝜑,𝜓𝜓,𝜓𝐾𝜓,𝜓𝐹¹₂, (43) where 𝐶=𝑐₋₁,𝑐₀,…,𝑐₃^𝑇,𝐷=𝑑_2,−1,…,𝑑_2,2,𝑑_3,−1,…,𝑑_3,6𝑑,…,_𝑀,−1,…,𝑑_𝑀,2^𝑀₋₂^𝑇,−=𝜑,𝜑𝐾𝜑,𝜑₁₀𝜑_𝑗0,𝑟(𝑥)𝜑_𝑗0,𝑖(𝑥)−𝑑𝑥₁₀𝜑_𝑗0,𝑟(𝑥)¹₀𝐾(𝑥,𝑡)𝜑_𝑗0,𝑖(𝑡)𝑑𝑡𝑑𝑥_𝑖,𝑟,−=𝜓,𝜑𝐾𝜓,𝜑₁₀𝜑_𝑗0,𝑟(𝑥)𝜓_𝑘,𝑗(𝑥)−𝑑𝑥₁₀𝜑_𝑗0,𝑟(𝑥)¹₀𝐾(𝑥,𝑡)𝜓_𝑘,𝑗(𝑡)𝑑𝑡𝑑𝑥_{𝑟,𝑘,𝑗},−=𝜑,𝜓𝐾𝜑,𝜓₁₀𝜓_𝑠,𝑙(𝑥)𝜑_𝑗0,𝑖(𝑥)−𝑑𝑥₁₀𝜓_𝑠,𝑙(𝑥)¹₀𝐾(𝑥,𝑡)𝜑_𝑗0,𝑖(𝑡)𝑑𝑡𝑑𝑥_{𝑖,𝑙,𝑠},−=𝜓,𝜓𝐾𝜓,𝜓₁₀𝜓_𝑠,𝑙(𝑥)𝜓_𝑘,𝑗(𝑥)−𝑑𝑥₁₀𝜓_𝑠,𝑙(𝑥)¹₀𝐾(𝑥,𝑡)𝜓_𝑘,𝑗(𝑡)𝑑𝑡𝑑𝑥_{𝑙,𝑠,𝑘,𝑗},𝐹₁=¹₀𝑓(𝑥)𝜑_𝑗0,𝑟(𝑥)𝐹𝑑𝑥,₂=¹₀𝑓(𝑥)𝜓_𝑠,𝑙(𝑥)𝑑𝑥, (44) and the subscripts 𝑖,𝑟,𝑘,𝑗,𝑙, and 𝑠 assume values as given below: 𝑖,𝑟=−1,…,2^𝑗0−1,𝑙,𝑘=𝑗₀,…,𝑀,𝑠,𝑗=−1,…,2^𝑀−2. (45) In fact, the entries with significant magnitude are in the ⟨𝐾𝜑,𝜑⟩−⟨𝜑,𝜑⟩ and ⟨𝐾𝜓,𝜓⟩−⟨𝜓,𝜓⟩ submatrices which are of order (2^𝑗0+1) and (2^𝑀+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 𝑥 yielding 𝑌^′(𝑥)=𝑓^′(𝑥)+¹₀𝐾^′(𝑥,𝑡)𝑦(𝑡)𝑑𝑡. (46) We apply variation iteration method for (46). According to this method, correction functional can be defined as 𝑦_𝑛+1(𝑥)=𝑦_𝑛(𝑥)+^𝑥₀𝜆(𝜉)𝑦^′_𝑛(𝜉)−𝑓^′(𝜉)−^𝑏_𝑎𝐾^′(𝜉,𝑡)𝑦_𝑛(𝑡)𝑑𝑡𝑑𝜉, (47) where 𝜆(𝜉) is a general Lagrange multiplier which can be identified optimally by the variational theory, the subscript 𝑛 denotes the 𝑛th order approximation, and 𝑦_𝑛 is considered as a restricted variation; that is, 𝛿𝑦_𝑛=0. The successive approximations 𝑦_𝑛(𝑥), 𝑛≥1 for the solution 𝑦(𝑥) can be readily obtained after determining the Lagrange multiplier and selecting an appropriate initial function 𝑦₀(𝑥). Consequently the approximate solution may be obtained by using 𝑦(𝑥)=lim_𝑛→∞𝑦_𝑛(𝑥). (48) To make the above correction functional stationary, we have 𝛿𝑦_𝑛+1(𝑥)=𝛿𝑦_𝑛(𝑥)+𝛿^𝑥₀𝜆(𝜉)𝑦^′_𝑛(𝜉)−𝑓^′(𝜉)−^𝑏_𝑎𝐾^′(𝜉,𝑡)𝑦_𝑛(𝑡)𝑑𝑡𝑑𝜉=𝛿𝑦_𝑛(𝑥)+^𝑥₀𝜆(𝜉)𝛿𝑦^′_𝑛(𝜉)𝑑𝜉=𝛿𝑦_𝑛(𝑥)+𝜆𝛿𝑦_{𝑛𝜉=𝑥}−^𝑥₀𝜆^′(𝜉)𝛿𝑦_𝑛(𝜉)𝑑𝜉. (49) Under stationary condition, 𝛿𝑦_𝑛+1=0 (50) implies the following Euler Lagrange equation: 𝜆^′(𝜉)=0, (51) with the following natural boundary condition: 1+𝜆(𝜉)_𝜉=𝑥=0. (52) Solving (51), along with boundary condition (52), we get the general Lagrange multiplier 𝜆=−1. (53) Substituting the identified Lagrange multiplier into (47) results in the following iterative scheme: 𝑦_𝑛+1(𝑥)=𝑦_𝑛(𝑥)−^𝑥₀𝑦^′_𝑛(𝜉)−𝑓^′(𝜉)−^𝑏_𝑎𝐾^′(𝜉,𝑡)𝑦_𝑛(𝑡)𝑑𝑡𝑑𝜉,𝑛≥0. (54) By starting with initial approximate function 𝑦₀(𝑥)=𝑓(𝑥) (say), we can determine the approximate solution 𝑦(𝑥) 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: ^𝑛_𝑗=1𝑦_𝑗(𝑥)=𝑓_𝑖(𝑥)+^𝑛_𝑗=1¹₀𝐾_𝑖,𝑗(𝑥,𝑡)𝑦_𝑗(𝑡)𝑑𝑡,𝑖=1,2,…,𝑛, (55) where 𝑓_𝑖(𝑥) and 𝐾_𝑖,𝑗(𝑥,𝑡) are known functions and 𝑦_𝑗(𝑥) are the unknown functions for 𝑖,𝑗=1,2,…,𝑛.

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 ℎ₀=1,0≤𝑥<1,0,others. (56) The mother wavelet function is defined as ℎ₁(𝑥)=11,0≤𝑥<2,1−1,2≤𝑥<1,0,others. (57) The family of wavelet functions generated by translation and dilation of ℎ₁(𝑥) are given by ℎ_𝑛(𝑥)=ℎ₁2^𝑗𝑥−𝑘, (58) where 𝑛=2^𝑗+𝑘, 𝑗≥0, 0≤𝑘<2^𝑗.

Mutual orthogonalities of all Haar wavelets can be expressed as ¹₀ℎ_𝑚(𝑥)ℎ_𝑛(𝑥)𝑑𝑥=2^−𝑗𝛿_𝑚𝑛=2^−𝑗,𝑚=𝑛=2^𝑗+𝑘,0,𝑚≠𝑛. (59)

4.1.2. Function Approximation

An arbitrary function 𝑦(𝑥)∈𝐿²[0,1) can be expanded into the following Haar series: 𝑦(𝑥)=^+∞_𝑛=0𝑐_𝑛ℎ_𝑛(𝑥), (60) where the coefficients 𝑐_𝑛 are given by 𝑐_𝑛=2^𝑗1₀𝑦(𝑥)ℎ_𝑛(𝑥)𝑑𝑥,𝑛=2^𝑗+𝑘,𝑗≥0,0≤𝑘<2^𝑗. (61) In particular, 𝑐₀=∫¹₀𝑦(𝑥)𝑑𝑥.

The previously mentioned expression in (60) can be approximately represented with finite terms as follows: 𝑦(𝑥)≈^𝑚−1_𝑛=0𝑐_𝑛ℎ_𝑛(𝑥)=𝐶^𝑇_(𝑚)ℎ_(𝑚)(𝑥), (62) where the coefficient vector 𝐶^𝑇_(𝑚) and the Haar function vector ℎ_(𝑚)(𝑥) are, respectively, defined as 𝐶^𝑇_(𝑚)=𝑐₀,𝑐₁,…,𝑐_𝑚−1,𝑚=2^𝑗,ℎ_(𝑚)(𝑥)=ℎ₀(𝑥),ℎ₁(𝑥),…,ℎ_𝑚−1(𝑥)^𝑇,𝑚=2^𝑗. (63) The Haar expansion for function 𝐾(𝑥,𝑡) of order 𝑚 is defined as follows: 𝐾(𝑥,𝑡)≈^𝑚−1_𝑢=0^𝑚−1_𝑣=0𝑎_𝑢𝑣ℎ_𝑣(𝑥)ℎ_𝑢(𝑡), (64) where 𝑎_𝑢𝑣=2^𝑖+𝑞∬¹₀𝐾(𝑥,𝑡)ℎ_𝑣(𝑥)ℎ_𝑢(𝑡)𝑑𝑥𝑑𝑡, 𝑢=2^𝑖+𝑗, 𝑣=2^𝑞+𝑟, 𝑖,𝑞≥0.

From (62) and (64), we obtain 𝐾(𝑥,𝑡)≈ℎ^𝑇_(𝑚)(𝑡)𝐾ℎ_(𝑚)(𝑥), (65) where 𝐾=𝑎_𝑢𝑣^𝑇_𝑚×𝑚. (66)

4.1.3. Operational Matrices of Integration

We define 𝐻_(𝑚)=ℎ_(𝑚)12𝑚,ℎ_(𝑚)32𝑚,…,ℎ_(𝑚)2𝑚−1,2𝑚 (67) where 𝐻₍₁₎=[1], 𝐻₍₂₎=¹¹₁₋₁.

Then, for 𝑚=4, the corresponding matrix can be represented as 𝐻₍₄₎=ℎ₍₄₎18,ℎ₍₄₎38,…,ℎ₍₄₎78=.111111−1−11−100001−1 (68) The integration of the Haar function vector ℎ_(𝑚)(𝑡) is ^𝑥₀ℎ_(𝑚)(𝑡)𝑑𝑡=𝑃_(𝑚)ℎ_(𝑚)(𝑥),𝑥∈[0,1), (69) where 𝑃_(𝑚) is the operational matrix of order 𝑚, and 𝑃₍₁₎=12,𝑃_(𝑚)=12𝑚2𝑚𝑃_(𝑚/2)−𝐻_(𝑚/2)𝐻⁻¹_(𝑚/2)0. (70) By recursion of the above formula, we obtain 𝑃₍₂₎=14,𝑃2−110₍₄₎=1,𝑃168−4−2−240−2211001−100₍₈₎=1.6432−16−8−8−4−4−4−4160−88−4−4444400−4400440000−441120000011−2000001−10200001−10−20000 (71) Therefore, we get 𝐻_(𝑚)⁻¹=1𝑚𝐻^𝑇_(𝑚)×diag1,1,2,2,2²,…,2²₂₂,…,2^𝛼−1,…,2^𝛼−1_2𝛼−1, (72) where 𝑚=2^𝛼 and 𝛼 is a positive integer.

The inner product of two Haar functions can be represented as ¹₀ℎ_(𝑚)(𝑡)ℎ^𝑇_(𝑚)(𝑡)𝑑𝑡=𝐷, (73) where 𝐷=diag1,1,1/2,1/2,1/2²,…,1/2²₂₂,…,1/2^𝛼−1,…,1/2^𝛼−1_2𝛼−1. (74)

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

Consider the following Fredholm integral equations system defined in (55): ^𝑚_𝑗=1𝑦_𝑗(𝑥)=𝑓_𝑖(𝑥)+^𝑚_𝑗=1¹₀𝐾_𝑖,𝑗(𝑥,𝑡)𝑦_𝑗(𝑡)𝑑𝑡,𝑖=1,2,…,𝑚. (75)

The Haar series of 𝑦_𝑗(𝑥) and 𝐾_𝑖,𝑗(𝑥,𝑡), 𝑖=1,2,…,𝑚; 𝑗=1,2,…,𝑚 are, respectively, expanded as 𝑦_𝑗(𝑥)≈𝑌^𝑇_𝑗ℎ_(𝑚)(𝑥)𝐾,𝑗=1,2,…,𝑚,_𝑖,𝑗(𝑥,𝑡)≈ℎ^𝑇_(𝑚)(𝑡)𝐾_𝑖,𝑗ℎ_(𝑚)(𝑥),𝑖,𝑗=1,2,…,𝑚. (76) Substituting (76) into (75), we get ^𝑚_𝑗=1𝑌^𝑇_𝑗ℎ_(𝑚)(𝑥)=𝑓_𝑖(𝑥)+^𝑚_𝑗=1¹₀𝑌^𝑇_𝑗ℎ_(𝑚)(𝑡)ℎ^𝑇_(𝑚)(𝑡)𝐾_𝑖,𝑗ℎ_(𝑚)(𝑥)𝑑𝑡,𝑖=1,2,…,𝑚. (77) From (77) and (73), we get ^𝑚_𝑗=1𝑌^𝑇_𝑗ℎ_(𝑚)(𝑥)=𝑓_𝑖(𝑥)+^𝑚_𝑗=1𝑌^𝑇_𝑗𝐷𝐾_𝑖,𝑗ℎ_(𝑚)(𝑥),𝑖=1,2,…,𝑚. (78) Interpolating 𝑚 collocation points, that is, {𝑥_𝑖}^𝑚_𝑖=1, in the interval [0,1] leads to the following algebraic system of equations: ^𝑚_𝑗=1𝑌^𝑇_𝑗ℎ_(𝑚)𝑥_𝑖=𝑓_𝑖𝑥_𝑖+^𝑚_𝑗=1𝑌^𝑇_𝑗𝐷𝐾_𝑖,𝑗ℎ_(𝑚)𝑥_𝑖,𝑖=1,2,…,𝑚. (79) Hence, 𝑌_𝑗, 𝑗=1,2,…,𝑚 can be computed by solving the above algebraic system of equations and consequently the solutions 𝑦_𝑗(𝑥)≈𝑌^𝑇_𝑗ℎ_(𝑚)(𝑥), 𝑗=1,2,…,𝑚.

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: 𝑦_𝑖(𝑥)=𝑓_𝑖(𝑥)+^𝑛_𝑗=1¹₀𝐾_𝑖,𝑗(𝑥,𝑡)𝑦_𝑗(𝑡)𝑑𝑡,𝑖=1,2,…,𝑛,0≤𝑥≤1. (80) A Taylor series expansion can be made for the solution of 𝑦_𝑗(𝑡) in the integral equation (80): 𝑦_𝑗(𝑡)=𝑦_𝑗(𝑥)+𝑦^′_𝑗(𝑥)+1)(𝑡−𝑥+⋯𝑦𝑚!^𝑗(𝑚)(𝑥)(𝑡−𝑥)^𝑚+𝐸(𝑡), (81) where 𝐸(𝑡) denotes the error between 𝑦_𝑗(𝑡) and its Taylor series expansion in (81).

If we use the first 𝑚 term of Taylor series expansion and neglect the term containing 𝐸(𝑡), that is, ∫¹₀∑^𝑛_𝑗=1𝐾_𝑖,𝑗(𝑥,𝑡)𝐸(𝑡)𝑑𝑡, then, substituting (81) for 𝑦_𝑗(𝑡) into the integral in (80), we have 𝑦_𝑖(𝑥)≈𝑓_𝑖(𝑥)+^𝑛_𝑗=1¹₀𝐾_𝑖,𝑗(𝑥,𝑡)^𝑚_𝑟=01()𝑟!𝑡−𝑥^𝑟𝑦_𝑗^(𝑟)(𝑥)𝑦𝑑𝑡,𝑖=1,2,…,𝑛,_𝑖(𝑥)≈𝑓_𝑖(𝑥)+^𝑛𝑚_𝑗=1𝑟=01𝑦𝑟!^𝑗(𝑟)(𝑥)¹₀𝐾_𝑖,𝑗(𝑥,𝑡)(𝑡−𝑥)^𝑟𝑦𝑑𝑡,𝑖=1,2,…,𝑛,_𝑖(𝑥)−^𝑛𝑚_𝑗=1𝑟=01𝑦𝑟!^𝑗(𝑟)(𝑥)¹₀𝐾_𝑖,𝑗(𝑥,𝑡)(𝑡−𝑥)^𝑟𝑑𝑡≈𝑓_𝑖(𝑥),𝑖=1,2,…,𝑛. (82)(83) 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 differentiate 𝑠 times both sides of (80) with respect to 𝑥; that is, 𝑦_𝑖^(𝑠)(𝑥)=𝑓_𝑖^(𝑠)(𝑥)+^𝑛_𝑗=1¹₀𝐾^(𝑠)_𝑖,𝑗(𝑥,𝑡)𝑦_𝑗(𝑡)𝑑𝑡,𝑖=1,2,…,𝑛,𝑠=1,2,…,𝑚, (84) where 𝐾^(𝑠)_𝑖,𝑗(𝑥,𝑡)=𝜕^(𝑠)𝐾_𝑖,𝑗(𝑥,𝑡)/𝜕𝑥^(𝑠), 𝑠=1,2,…,𝑚.

Applying the mean value theorem for integral in (84), we have 𝑦_𝑖^(𝑠)(𝑥)−^𝑛_𝑗=1¹₀𝐾^(𝑠)_𝑖,𝑗(𝑥,𝑡)𝑑𝑡𝑦_𝑗(𝑥)≈𝑓_𝑖^(𝑠)(𝑥),𝑖=1,2,…,𝑛,𝑠=1,2,…,𝑚. (85) 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].

An 𝑚-set of BPF is defined as follows: Φ_𝑖(𝑡)=()𝑇1,𝑖−1𝑚𝑇≤𝑡<𝑖𝑚,0,otherwise (86) with 𝑡∈[0,𝑇), 𝑇/𝑚=ℎ and 𝑖=1,2,…,𝑚.

4.3.1. Properties of BPF

(1) Disjointness. One hasΦ_𝑖(𝑡)Φ_𝑗(𝑡)=Φ_𝑖(𝑡),𝑖=𝑗;0,𝑖≠𝑗, (87), 𝑖,𝑗=1,2,…,𝑚. This property is obtained from definition of BPF.

(2) Orthogonality. One has^𝑇₀Φ_𝑖(𝑡)Φ_𝑗(𝑡)𝑑𝑡=ℎ,𝑖=𝑗;0,𝑖≠𝑗, (88)𝑡∈[0,𝑇), 𝑖,𝑗=1,2,…,𝑚. This property is obtained from the disjointness property.

(3) Completeness. For every 𝑓∈𝐿²,{Φ} is complete; if ∫Φ𝑓=0 then 𝑓=0 almost everywhere. Because of completeness of {Φ}, we have ^𝑇₀𝑓²(𝑡)𝑑𝑡=^∞_𝑖=1𝑓²_𝑖Φ_𝑖(𝑡)² (89) for every real bounded function 𝑓(𝑡) which is square integrable in the interval 𝑡∈[0,𝑇) and 𝑓_𝑖=(1/ℎ)𝑓(𝑡)Φ_𝑖(𝑡)𝑑𝑡.

4.3.2. Function Approximation

The orthogonality property of BPF is the basis of expanding functions into their Block-Pulse series. For every 𝑓(𝑡)∈𝐿²(𝑅), 𝑓(𝑡)=^𝑚_𝑖=1𝑓_𝑖Φ_𝑖(𝑡), (90) where 𝑓_𝑖 is the coefficient of Block-Pulse function, with respect to 𝑖th Block-Pulse function Φ_𝑖(𝑡).

The criterion of this approximation is that mean square error between 𝑓(𝑡) and its expansion is minimum 1𝜀=𝑇^𝑇₀𝑓(𝑡)−^𝑚_𝑗=1𝑓_𝑗Φ_𝑗(𝑡)²𝑑𝑡 (91) so that we can evaluate Block-Pulse coefficients. Now𝜕𝜀𝜕𝑓_𝑖2=−𝑇^𝑇₀𝑓(𝑡)−^𝑚_𝑗=1𝑓_𝑗Φ_𝑗(𝑡)Φ_𝑖(𝑡)𝑑𝑡=0,⟹𝑓_𝑖=1ℎ^𝑇₀𝑓(𝑡)Φ_𝑖(𝑡)𝑑𝑡usingorthogonalproperty. (92) In the matrix form, we obtain the following from (90) as follow: 𝑓(𝑡)=^𝑚_𝑖=1𝑓_𝑖Φ_𝑖(𝑡)=𝐹^𝑇Φ(𝑡)=Φ^𝑇𝐹where𝐹=𝑓₁,𝑓₂,…,𝑓_𝑚^𝑇,Φ(𝑡)=Φ₁(𝑡),Φ₂(𝑡),…,Φ_𝑚(𝑡)^𝑇. (93) Now let 𝐾(𝑡,𝑠) be two-variable function defined on 𝑡∈[0,𝑇) and 𝑠∈[0,1); then 𝐾(𝑡,𝑠) can be expanded to BPF as 𝐾(𝑡,𝑠)=Φ^𝑇(𝑡)𝐾Ψ(𝑠), (94) where Φ(𝑡) and Ψ(𝑠) are 𝑚₁ and 𝑚₂ dimensional Block-Pulse function vectors and 𝑘 is a 𝑚₁×𝑚₂ Block-Pulse coefficient matrix.

There are two different cases of multiplication of two BPF. The first case is Φ(𝑡)Φ^𝑇(𝑡)=Φ₁(𝑡)0⋯00Φ₂(𝑡)⋯0⋮⋮⋱⋮00⋯Φ_𝑚(𝑡). (95) It is obtained from disjointness property of BPF. It is a diagonal matrix with 𝑚 Block-Pulse functions.

The second case is Φ^𝑇(𝑡)Φ(𝑡)=1 (96) because ∑^𝑚_𝑖=1(Φ_𝑖(𝑡))²=∑^𝑚_𝑖=1Φ_𝑖(𝑡)=1.

Operational Matrix of Integration. BPF integration property can be expressed by an operational equation as ^𝑇₀Φ(𝑡)𝑑𝑡=𝑃Φ(𝑡), (97) where Φ(𝑡)=Φ₁(𝑡),Φ₂(𝑡),…,Φ_𝑚(𝑡)^𝑇. (98) A general formula for 𝑃_𝑚×𝑚 can be written as 1𝑃=2122⋯2012⋯2001⋯2⋮⋮⋮⋱⋮000⋯1. (99) By using this matrix, we can express the integral of a function 𝑓(𝑡) into its Block-Pulse series ^𝑡₀𝑓(𝑡)𝑑𝑡=^𝑡₀𝐹^𝑇Φ(𝑡)𝑑𝑡=𝐹^𝑇𝑃Φ(𝑡). (100)

4.3.3. Solution for Linear Integral Equations System

Consider the integral equations system from (55) as follows: ^𝑛_𝑗=1𝑦_𝑗(𝑥)=𝑓_𝑖(𝑥)+^𝑛_𝑗=1^𝛽_𝛼𝐾_𝑖,𝑗(𝑥,𝑡)𝑦_𝑗(𝑡)𝑑𝑡,𝑖=1,2,…,𝑛. (101) Block-Pulse coefficients of 𝑦_𝑗(𝑥), 𝑗=1,2,…,𝑛 in the interval 𝑥∈[𝛼,𝛽) can be determined from the known functions 𝑓_𝑖(𝑥), 𝑖=1,2,…,𝑛 and the kernels 𝐾_𝑖,𝑗(𝑥,𝑡), 𝑖,𝑗=1,2,…𝑛. Usually, we consider 𝛼=0 to facilitie the use of Block-Pulse functions. In case 𝛼≠0, we set 𝑋=((𝑥−𝛼)/(𝛽−𝛼))𝑇, where 𝑇=𝑚ℎ.

We approximate 𝑓_𝑖(𝑥), 𝑦_𝑗(𝑥), 𝐾_𝑖,𝑗(𝑥,𝑡) by its BPF as follows: 𝑓_𝑖(𝑥)≈𝐹^𝑇_𝑖Φ(𝑥),𝑦_𝑗(𝑥)≈𝑌^𝑇_𝑗Φ(𝑥),𝐾_𝑖,𝑗(𝑥,𝑡)≈Φ^𝑇(𝑡)𝐾_𝑖,𝑗Φ(𝑥), (102) where 𝐹_𝑖, 𝑌_𝑗, and 𝐾_𝑖,𝑗 are defined in Section 4.3.2, and substituting (102) into (101), we have ^𝑛_𝑗=1𝑌^𝑇_𝑗Φ(𝑥)=𝐹^𝑇_𝑖Φ(𝑥)+^𝑛_𝑗=10^𝑚ℎ𝑌^𝑇_𝑗Φ(𝑡)Φ^𝑇(𝑡)𝐾_𝑖,𝑗Φ(𝑥)𝑑𝑡,𝑖=1,2,…,𝑛,_𝑛𝑗=1𝑌^𝑇_𝑗Φ(𝑥)=𝐹^𝑇_𝑖Φ(𝑥)+^𝑛_𝑗=1𝑌^𝑇_𝑗ℎ𝐼𝐾_𝑖,𝑗Φ(𝑥),𝑖=1,2,…,𝑛, (103)(104) since ₀^𝑚ℎΦ(𝑡)Φ^𝑇(𝑡)𝑑𝑡=ℎ𝐼. (105) From (104), we get ^𝑛_𝑗=1𝐼−ℎ𝐾^𝑇_𝑖,𝑗𝑌_𝑗=𝐹_𝑖,𝑖=1,2,…,𝑛. (106) Set 𝐴_𝑖,𝑗=𝐼−ℎ𝐾^𝑇_𝑖,𝑗; then we have from (106) ^𝑛_𝑗=1𝐴_𝑖,𝑗𝑌_𝑗=𝐹_𝑖,𝑖=1,2,…,𝑛 (107) which is a linear system 𝐴₁₁𝐴₁₂...𝐴_1𝑛𝐴₂₁𝐴₂₂...𝐴_2𝑛𝐴............_𝑛1𝐴_𝑛2...𝐴_𝑛𝑛𝑌₁𝑌₂...𝑌_𝑛=𝐹₁𝐹₂...𝐹_𝑛. (108) After solving the above system we can find 𝑌_𝑗, 𝑗=1,2,…,𝑛 and hence obtain the solutions 𝑦_𝑗=Φ^𝑇𝑌_𝑗, 𝑗=1,2,…,𝑛.

5. Numerical Methods for Nonlinear Fredholm-Hammerstein Integral Equation

We consider the second kind nonlinear Fredholm integral equation of the following form: 𝑢(𝑥)=𝑓(𝑥)+¹₀𝐾(𝑥,𝑡)𝐹(𝑡,𝑢(𝑡))𝑑𝑡,0≤𝑥≤1, (109) where 𝐾(𝑥,𝑡) is the kernel of the integral equation, 𝑓(𝑥) and 𝐾(𝑥,𝑡) are known functions, and 𝑢(𝑥) 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 𝑦(𝑥)=𝐹(𝑥,𝑢(𝑥)),0≤𝑥≤1. (110) Now, from (16), we can approximate the functions 𝑢(𝑥) and 𝑦(𝑥) as 𝑢(𝑥)=𝐴^𝑇Ψ(𝑥),𝑦(𝑥)=𝐵^𝑇Ψ(𝑥), (111) where 𝐴 and 𝐵 are (2^𝑀+1+𝑚−1)×1 column vectors similar to 𝐶 defined in (17).

Again, by using dual of the wavelet functions, we can approximate the functions 𝑓(𝑥) and 𝐾(𝑥,𝑡) as follows: 𝐹(𝑥)=𝐷^𝑇Ψ(𝑥),𝐾(𝑥,𝑡)=Ψ^𝑇(𝑡)ΘΨ(𝑥), (112) where Θ_(𝑖,𝑗)=¹₀¹₀𝐾(𝑥,𝑡)Ψ_𝑖(𝑡)𝑑𝑡Ψ_𝑗(𝑥)𝑑𝑥. (113) From (110)–(112), we get ¹₀𝐾()𝐹((𝑡=𝑥,𝑡𝑡,𝑢))𝑑𝑡₁₀𝐵^𝑇Ψ(𝑡)Ψ^𝑇(𝑡)ΘΨ(𝑥)𝑑𝑡=𝐵^𝑇1₀Ψ(𝑡)Ψ^𝑇(𝑡)𝑑𝑡ΘΨ(𝑥)=𝐵^𝑇ΘΨ(𝑥),since¹₀Ψ(𝑡)Ψ^𝑇(𝑡)𝑑𝑡=𝐼. (114) Applying (110)–(114) in (109), we get 𝐴^𝑇Ψ(𝑥)=𝐷^𝑇Ψ(𝑥)+𝐵^𝑇ΘΨ(𝑥). (115) Multiplying (115) by Ψ^𝑇(𝑥) both sides from the right and integrating both sides with respect to 𝑥 from 0 to 1, we have 𝐴^𝑇𝑃=𝐷^𝑇+𝐵^𝑇𝐴Θ,_𝑇𝑃−𝐷^𝑇−𝐵^𝑇Θ=0, (116) where 𝑃 is a (2^𝑀+1+𝑚−1)×(2^𝑀+1+𝑚−1) square matrix given by 𝑃=¹₀Ψ(𝑥)Ψ^𝑇(𝑥)𝑃𝑑𝑥=¹𝑃₂₁₀,Ψ(𝑥)Ψ^𝑇(𝑥)𝑑𝑥=𝐼. (117) Equation (116) gives a system of (2^𝑀+1+𝑚−1) algebraic equations with 2(2^𝑀+1+𝑚−1) unknowns for 𝐴 and 𝐵 vectors given in (111).

To find the solution 𝑢(𝑥) in (111), we first utilize the following equation: 𝐹𝑥,𝐴^𝑇Ψ(𝑥)=𝐵^𝑇Ψ(𝑥), (118) with the collocation points 𝑥_𝑖=(𝑖−1)/(2^𝑀+1+𝑚−2), where 𝑖=1,2,…,2^𝑀+1+𝑚−1.

Equation (118) gives a system of (2^𝑀+1+𝑚−1) algebraic equations with 2(2^𝑀+1+𝑚−1) unknowns, for 𝐴 and 𝐵 vectors given in (111).

Combining (116) and (118), we have a total of 2(2^𝑀+1+𝑚−1) system of algebraic equations with 2(2^𝑀+1+𝑚−1) unknowns for 𝐴 and 𝐵. Solving those equations for the unknown coefficients in the vectors 𝐴 and 𝐵, we can obtain the solution 𝑢(𝑥)=𝐴^𝑇Ψ(𝑥).

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^𝑏_𝑎𝑓(𝑥)𝑑𝑥=^𝑛−1_𝑖=1^𝑥𝑖+1_𝑥𝑖−1𝑓(𝑥)=ℎ𝑑𝑥3𝑓(𝑎)+4ℎ3^𝑛/2_𝑖=1𝑓𝑥_2𝑖−1+2ℎ3^(𝑛−1)/2_𝑖=1𝑓𝑥_2𝑖+ℎ3𝑓(𝑏)−(𝑏−𝑎)ℎ180⁴𝑓⁽⁴⁾𝜂. (119)

5.2.2. Modified Trapezoid Rule

One has^𝑏_𝑎𝑓(𝑥)𝑑𝑥=^𝑛_𝑖=1^𝑥𝑖_𝑥𝑖−1𝑓(𝑥)=ℎ𝑑𝑥2𝑓(𝑎)+ℎ^𝑛−1_𝑖=1𝑓𝑥_𝑖+ℎ2𝑓(𝑏)+ℎ²𝑓12^′(𝑎)−𝑓^′(𝑏). (120) Consider the nonlinear Fredholm integral equation of second kind defined in (109) as follows: 𝑢(𝑥)=𝑓(𝑥)+^𝑏_𝑎𝐾(𝑥,𝑡)𝐹(𝑢(𝑡))𝑑𝑡,𝑎≤𝑥≤𝑏. (121) 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𝑢(𝑥)(𝑥)+ℎ=𝑓3𝐾𝑥,𝑡₀𝐹𝑢₀+4^𝑛/2_𝑗=1𝐾𝑥,𝑡_2𝑗−1𝐹𝑢_2𝑗−1+2^(𝑛/2)−1_𝑗=1𝐾𝑥,𝑡_2𝑗𝐹𝑢_2𝑗+𝐾𝑥,𝑡_𝑛𝐹𝑢_𝑛. (122) Hence, for 𝑥=𝑥₀,𝑥₁,…,𝑥_𝑛 and 𝑡=𝑡₀,𝑡₁,…,𝑡_𝑛 in (122), we have 𝑢𝑥_𝑖=𝑓𝑥_𝑖+ℎ3𝐾𝑥_𝑖,𝑡₀𝐹𝑢₀+4^𝑛/2_𝑗=1𝐾𝑥_𝑖,𝑡_2𝑗−1𝐹𝑢_2𝑗−1+2^(𝑛/2)−1_𝑗=1𝐾𝑥_𝑖,𝑡_2𝑗𝐹𝑢_2𝑗+𝐾𝑥_𝑖,𝑡_𝑛𝐹𝑢_𝑛. (123) Equation (123) is a nonlinear system of equations and, by solving (123), we obtain the unknowns 𝑢(𝑥_𝑖) for 𝑖=0,1,…,𝑛.

5.2.4. Modified Trapezoid Rule

One has𝑢(𝑥)(𝑥)+ℎ=𝑓2𝐾𝑥,𝑡₀𝐹𝑢₀+ℎ^𝑛−1_𝑗=1𝐾𝑥,𝑡_𝑗𝐹𝑢_𝑗+ℎ2𝐾𝑥,𝑡_𝑛𝐹𝑢_𝑛+ℎ²𝐽12𝑥,𝑡₀𝐹𝑢₀+𝐾𝑥,𝑡₀𝑢^′₀𝐹^′𝑢₀−𝐽𝑥,𝑡_𝑛𝐹𝑢_𝑛−𝐾𝑥,𝑡_𝑛𝑢^′_𝑛𝐹^′𝑢_𝑛, (124) where 𝐽(𝑥,𝑡)=𝜕𝐾(𝑥,𝑡)/𝜕𝑡.

For 𝑥=𝑥₀,𝑥₁,…,𝑥_𝑛 and 𝑡=𝑡₀,𝑡₁,…,𝑡_𝑛 in (124), we have 𝑢𝑥_𝑖=𝑓𝑥_𝑖+ℎ2𝐾𝑥_𝑖,𝑡₀𝐹𝑢₀+ℎ^𝑛−1_𝑗=1𝐾𝑥_𝑖,𝑡_𝑗𝐹𝑢_𝑗+ℎ2𝐾𝑥_𝑖,𝑡_𝑛𝐹𝑢_𝑛+ℎ²𝐽𝑥12^𝑖,𝑡₀𝐹𝑢₀+𝐾𝑥_𝑖,𝑡₀𝑢^′₀𝐹^′𝑢₀−𝐽𝑥_𝑖,𝑡_𝑛𝐹𝑢_𝑛−𝐾𝑥_𝑖,𝑡_𝑛𝑢^′_𝑛𝐹^′𝑢_𝑛, (125) for 𝑖=0,1,…,𝑛.

This is a system of (𝑛+1) equations and (𝑛+3) unknowns. By taking derivative from (121) and setting 𝐻(𝑥,𝑡)=𝜕𝐾(𝑥,𝑡)/𝜕𝑥, we obtain 𝑢^′(𝑥)=𝑓^′(𝑥)+^𝑏_𝑎𝐻(𝑥,𝑡)𝐹(𝑢(𝑡))𝑑𝑡,𝑎≤𝑥≤𝑏. (126) If 𝑢 is a solution of (121), then it is also solution of (126). By using trapezoid rule for (126) and replacing 𝑥=𝑥_𝑖, we get 𝑢^′𝑥_𝑖=𝑓^′𝑥_𝑖+ℎ2𝐻𝑥_𝑖,𝑡₀𝐹𝑢₀+ℎ^𝑛−1_𝑗=1𝐻𝑥_𝑖,𝑡_𝑗𝐹𝑢_𝑗+ℎ2𝐻𝑥_𝑖,𝑡_𝑛𝐹𝑢_𝑛, (127) for 𝑖=0,1,…,𝑛. In case of 𝑖=0,𝑛 from system (127), we obtain two equations.

Now (127) combined with (125) generates the nonlinear system of equations as follows: 𝑢𝑥_𝑖=ℎ2𝐾𝑥_𝑖,𝑡₀+ℎ²𝐽𝑥12^𝑖,𝑡₀𝐹𝑢₀+ℎ^𝑛−1_𝑗=1𝐾𝑥_𝑖,𝑡_𝑗𝐹𝑢_𝑗+ℎ2𝐾𝑥_𝑖,𝑡_𝑛−ℎ²𝐽𝑥12^𝑖,𝑡_𝑛𝐹𝑢_𝑛+ℎ²𝐾𝑥12^𝑖,𝑡₀𝑢^′₀𝐹^′𝑢₀−𝐾𝑥_𝑖,𝑡_𝑛𝑢^′_𝑛𝐹^′𝑢_𝑛,𝑢^′𝑥₀=𝑓^′𝑥₀+ℎ2𝐻𝑥₀,𝑡₀𝐹𝑢₀+ℎ^𝑛−1_𝑗=1𝐻𝑥₀,𝑡_𝑗𝐹𝑢_𝑗+ℎ2𝐻𝑥₀,𝑡_𝑛𝐹𝑢_𝑛,𝑢^′𝑥_𝑛=𝑓^′𝑥_𝑛+ℎ2𝐻𝑥_𝑛,𝑡₀𝐹𝑢₀+ℎ^𝑛−1_𝑗=1𝐻𝑥_𝑛,𝑡_𝑗𝐹𝑢_𝑗+ℎ2𝐻𝑥_𝑛,𝑡_𝑛𝐹𝑢_𝑛. (128) By solving this system with (𝑛+3) nonlinear equations and (𝑛+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 parameter 𝑎 and the translation parameter 𝑏𝜓_𝑎,𝑏(𝑡)=|𝑎|^−1/2𝜓𝑡−𝑏𝑎,𝑎,𝑏∈𝑅,𝑎≠0. (129) Legendre wavelets 𝜓_𝑚,𝑛(𝑡)=𝜓(𝑘,𝑛,𝑚,𝑡) have four arguments; 𝑘=2,3,…, 𝑛=2𝑛−1, 𝑛=1,2,…,2^𝑘−1, 𝑚 is the order for Legendre polynomials and 𝑡 is the normalized time.

Legendre wavelets are defined on [0, 1) by 𝜓_𝑚,𝑛(𝑡)=1𝑚+2^1/22^𝑘/2𝐿_𝑚2^𝑘𝑡−𝑛,𝑛−12^𝑘≤𝑡<𝑛+12^𝑘,0,otherwise, (130) where 𝐿_𝑚(𝑡) are the well-known Legendre polynomials of order m, which are orthogonal with respect to the weight function 𝑤(𝑡)=1 and satisfy the following recursive formula: 𝐿₀(𝑡)𝐿=1,₁(𝑡)𝐿=𝑡,_𝑚+1(𝑡)=2𝑚+1𝑚+1𝑡𝐿_𝑚(𝑡)−𝑚𝐿𝑚+1^𝑚−1(𝑡),𝑚=1,2,3,…. (131) The set of Legendre wavelets are an orthonormal set.

5.3.1. Function Approximation

A function 𝑓(𝑥)∈𝐿²[0,1] can be expanded as 𝑓(𝑥)=^∞∞_𝑛=1𝑚=0𝑐_𝑛,𝑚𝜓_𝑛,𝑚(𝑥), (132) where 𝑐_𝑛,𝑚=𝑓(𝑥),𝜓_𝑛,𝑚(𝑥). (133) If the infinite series in (132) is truncated, then (132) can be written as 𝑓(𝑥)≈^2𝑘−1_𝑛=1^𝑀−1_𝑚=0𝑐_𝑛,𝑚𝜓_𝑛,𝑚(𝑥)=𝐶^𝑇Ψ(𝑥), (134) where 𝐶 and Ψ(𝑥) are 2^𝑘−1𝑀×1 matrices given by 𝐶=𝑐_1,0,𝑐_1,1,…,𝑐_1,𝑀−1,𝑐_2,0𝑐,…,_2,𝑀−1,…,𝑐₂^𝑘−1_,0,…,𝑐₂^𝑘−1_,𝑀−1^𝑇,Ψ(𝑥)=𝜓_1,0(𝑥),…,𝜓_1,𝑀−1(𝑥),𝜓_2,0(𝑥),…,𝜓_2,𝑀−1(𝑥)𝜓,…,_2𝑘−1,0(𝑥),…,𝜓_{2𝑘−1,𝑀−1}(𝑥)^𝑇. (135)(136) Similarly, a function 𝑘(𝑥,𝑡)∈𝐿²([0,1]×[0,1]) can be approximated as 𝑘(𝑥,𝑡)≈Ψ^𝑇(𝑡)𝐾Ψ(𝑥), (137) where 𝐾 is (2^𝑘−1𝑀×2^𝑘−1𝑀) matrix, with 𝐾_𝑖,𝑗=𝜓_𝑖(𝑡),𝑘(𝑥,𝑡),𝜓_𝑗(𝑥). (138) Also, the integer power of a function can be approximated as 𝑦(𝑥)^𝑝=𝑌^𝑇Ψ(𝑥)^𝑝=𝑌^∗_𝑝^𝑇Ψ(𝑥), (139) where 𝑌^∗_𝑝 is a column vector, whose elements are nonlinear combinations of the elements of the vector 𝑌. 𝑌^∗_𝑝 is called the operational vector of the 𝑝th power of the function 𝑦(𝑥).

5.3.2. The Operational Matrices

The integration of the vector Ψ(𝑥) defined in (136) can be obtained as ^𝑡₀Ψ𝑡^′𝑑𝑡^′=𝑃Ψ(𝑡), (140) where 𝑃 is the (2^𝑘−1𝑀×2^𝑘−1𝑀) operational matrix for integration and is given in [23] as 𝑃=𝐿𝐻⋯𝐻𝐻0𝐿⋯𝐻𝐻⋮⋮⋱⋮⋮00⋯𝐿𝐻00⋯0𝐿. (141) In (141), 𝐻 and 𝐿 are (𝑀×𝑀) matrices given in [23] as1𝐻=2^𝑘,120⋯000⋯0⋮⋮⋱⋮00⋯0𝐿=2^𝑘11√3−√00⋯00330√33√5√0⋯000−55√30√55√7√⋯0000−77√5√0⋯00⋮⋮⋮⋮⋱⋮0000⋯02𝑀−3(2𝑀−3)√−√2𝑀−10000⋯2𝑀−1(2𝑀−1)√0.2𝑀−3 (142)The integration of the product of two Legendre wavelets vector functions is obtained as ¹₀Ψ(𝑡)Ψ^𝑇(𝑡)𝑑𝑡=𝐼, (143) where 𝐼 is an identity matrix.

The product of two Legendre wavelet vector functions is defined as Ψ(𝑡)Ψ^𝑇(𝑡)𝐶=𝐶^𝑇Ψ(𝑡), (144) where 𝐶 is a vector given in (135) and 𝐶 is (2^𝑘−1𝑀×2^𝑘−1𝑀) 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) where 𝑓∈𝐿²[0,1], 𝑘∈𝐿²([0,1]×[0,1]), 𝑦 is an unknown function, and 𝑝 is a positive integer.

We can approximate the following functions as 𝑓(𝑥)≈𝐹^𝑇Ψ(𝑥),𝑦(𝑥)≈𝑌^𝑇Ψ(𝑥),𝑘(𝑥,𝑡)≈Ψ^𝑇(𝑡)(𝑥),𝑦(𝑥𝐾Ψ)^𝑝≈𝑌^∗𝑇Ψ(𝑥). (146) Substituting (146) into (145), we have 𝑌^𝑇Ψ(𝑥)=𝐹^𝑇Ψ(𝑥)+¹₀𝑌^∗𝑇Ψ(𝑡)Ψ^𝑇(𝑡)𝐾Ψ(𝑥)𝑑𝑡=𝐹^𝑇Ψ(𝑥)+𝑌^∗𝑇1₀Ψ(𝑡)Ψ^𝑇(𝑡)𝑑𝑡𝐾Ψ(𝑥)=𝐹^𝑇Ψ(𝑥)+𝑌^∗𝑇(𝑥)=𝐹𝐾Ψ_𝑇+𝑌^∗𝑇𝐾Ψ(𝑥)⟹𝑌^𝑇−𝑌^∗𝑇𝐾−𝐹^𝑇=0. (147) Equation (147) is a system of algebraic equations. Solving (147), we can obtain the solution 𝑦(𝑥)≈𝑌^𝑇Ψ(𝑥).

5.4. Homotopy Perturbation Method

Consider the following nonlinear Fredholm integral equation of second kind of the form 𝑢(𝑥)=𝑓(𝑥)+¹₀𝐾(𝑥,𝑡)𝐹(𝑢(𝑡))𝑑𝑡,0≤𝑥≤1. (148) For solving (148) by Homotopy perturbation method (HPM) [14–16], we consider (148) as 𝐿(𝑢)=𝑢(𝑥)−𝑓(𝑥)−¹₀𝐾(𝑥,𝑡)𝐹(𝑢(𝑡))𝑑𝑡=0. (149) As a possible remedy, we can define 𝐻(𝑢,𝑝) by 𝐻()(𝑢),𝐻()(𝑢),𝑢,0=𝑁𝑢,1=𝐿 (150) where 𝑁(𝑢) is an integral operator with known solution 𝑢₀.

We may choose a convex homotopy by 𝐻𝑢,𝑝=1−𝑝𝑁(𝑢)+𝑝𝐿(𝑢)=0 (151) and continuously trace an implicitly defined curve from a starting point 𝐻(𝑢₀,0) to a solution function 𝐻(𝑈,1). The embedding parameter 𝑝 monotonically increases from zero to unit as the trivial problem 𝐿(𝑢)=0. The embedding parameter 𝑝∈(0,1] can be considered as an expanding parameter. The HPM uses the homotopy parameter 𝑝 as an expanding parameter; that is,𝑢=𝑢₀+𝑝𝑢₁+𝑝²𝑢₂+⋯. (152a) When 𝑝→1, (152a) corresponding to (151) become the approximate solution of (149) as follows: 𝑈=lim_𝑝→1𝑢=𝑢₀+𝑢₁+𝑢₂+⋯. (152b)The series in (152b) converges in most cases, and the rate of convergence depends on 𝐿(𝑢) [14].

Consider 𝑁(𝑢)=𝑢(𝑥)−𝑓(𝑥). (153) The nonlinear term 𝐹(𝑢) can be expressed in He polynomials [25] as 𝐹(𝑢)=^∞_𝑚=0𝑝^𝑚𝐻_𝑚𝑢₀,𝑢₁,…,𝑢_𝑚=𝐻₀𝑢₀+𝑝𝐻₁𝑢₀,𝑢₁+⋯+𝑝^𝑚𝐻_𝑚𝑢₀,𝑢₁,…𝑢_𝑚+⋯, (154) where 𝐻_𝑚𝑢₀,𝑢₁,…,𝑢_𝑚=1𝜕𝑚!𝜕𝑝^𝑚𝑚𝐹^𝑚_𝑘=0𝑝^𝑘𝑢_𝑘𝑝=0,𝑚≥0. (155) Substituting (152a), (153), and (154) into (151), we have 1−𝑝𝑢₀+𝑝𝑢₁(𝑥)𝑢+⋯−𝑓+𝑝₀+𝑝𝑢₁(𝑥)−+⋯−𝑓₁₀𝐾(𝑥,𝑡)^∞_𝑚=0𝑝^𝑚𝐻_𝑚𝑢₀,𝑢₁,…,𝑢_𝑚⟹𝑢𝑑𝑡=0₀+𝑝𝑢₁+⋯−𝑓(𝑥)−𝑝¹₀𝐾(𝑥,𝑡)^∞_𝑚=0𝑝^𝑚𝐻_𝑚𝑢₀,𝑢₁,…,𝑢_𝑚𝑑𝑡=0. (156) Equating the terms with identical power of 𝑝 in (156), we have 𝑝⁰∶𝑢₀(𝑥)−𝑓(𝑥)=0⟹𝑢₀(𝑥)(𝑥)𝑝=𝑓₁∶𝑢₁(𝑥)−¹₀𝐾(𝑥,𝑡)𝐻₀𝑑𝑡=0⟹𝑢₁(𝑥)=¹₀𝐾(𝑥,𝑡)𝐻₀𝑝𝑑𝑡₂∶𝑢₂(𝑥)−¹₀𝐾(𝑥,𝑡)𝐻₁𝑑𝑡=0⟹𝑢₂(𝑥)=¹₀𝐾(𝑥,𝑡)𝐻₁⋮𝑑𝑡 (157) and in general form we have 𝑢₀(𝑥)(𝑥),𝑢=𝑓_𝑛+1(𝑥)=¹₀𝐾(𝑥,𝑡)𝐻_𝑛𝑑𝑡,𝑛=0,1,2,…. (158) 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) where 𝐿(𝑢(𝑡)) and 𝑁(𝑢(𝑡)) are the linear and nonlinear terms, respectively.

The Adomian decomposition method (ADM) consists of representing 𝑢(𝑥) as a series 𝑢(𝑥)=^∞_𝑚=0𝑢_𝑚(𝑥). (160) In the view of ADM, the nonlinear term 𝑁𝑢 can be represented as 𝑁𝑢=^∞_𝑛=0𝐴_𝑛,where𝐴_𝑛=1𝜕𝑛!𝜕𝜆^𝑛𝑛𝑁^∞_𝑘=0𝜆^𝑘𝑢_𝑘𝜆=0. (161)(162) Now substituting (160) and (161) into (159), we have ^∞_𝑚=0𝑢_𝑚(𝑥)=𝑓(𝑥)+^𝑏_𝑎𝐾(𝑥,𝑡)𝐿^∞_𝑚=0𝑢_𝑚(𝑡)+^∞_𝑚=0𝐴_𝑚𝑑𝑡, (163) and, then, ADM uses the recursive relations 𝑢₀(𝑥)(𝑥),𝑢=𝑓_𝑚(𝑥)=^𝑏_𝑎𝐾(𝑥,𝑡)𝐿𝑢_𝑚−1(𝑡)+𝐴_𝑚−1(𝑡)𝑑𝑡,𝑚≥1, (164) where 𝐴_𝑚 is so-called Adomian polynomial.

Therefore, we obtain the 𝑛-terms approximate solution as 𝜑_𝑛=𝑢₀+𝑢₁+⋯+𝑢_𝑛 (165) with 𝑢(𝑥)=lim_𝑛→∞𝜑_𝑛. (166)

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 semiorthogonal B-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.