Modified homotopy perturbation method for solving Fredholm integral equations

Available online at www.sciencedirect.com

Chaos, Solitons and Fractals 37 (2008) 1528–1537 www.elsevier.com/locate/chaos

Modified homotopy perturbation method for solving Fredholm integral equations A. Golbabai *, B. Keramati Department of Applied Mathematics, Faculty of Science, Iran University of Science and Technology, Narmak, Tehran 16844, Iran Accepted 23 October 2006

Abstract In this paper we present a modification to homotopy perturbation method for solving linear Fredholm integral equations. Comparisons are made between the standard HPM and the modified one. The results reveal that the proposed method is very effective and simple and gives the exact solution. Ó 2006 Elsevier Ltd. All rights reserved.

1. Introduction Homotopy perturbation method has been recently intensively studied by scientists and engineers and used for solving nonlinear problems. This method [6,9] was proposed first by He which is, in fact, a coupling of the traditional perturbation method and homotopy in topology. HPM method yields a very rapid convergence of the solution series in most cases, usually only few iterations leading to very accurate solutions. This method has been applied to many problems [4,5,14,15], specially in integral equations [1–4,7,8]. Consider the Fredholm integral equation: Z b cðxÞ ¼ f ðxÞ þ kðx; tÞcðtÞdt; c 6 x 6 d; a

let LðuÞ ¼ uðxÞ  f ðxÞ 

Z

b

kðx; tÞuðtÞdt ¼ 0;

ð1Þ

a

with solution u(x) = c(x), we define the homotopy H(u, p) by H ðu; 0Þ ¼ F ðuÞ;

H ðu; 1Þ ¼ LðuÞ;

where F(u) is a functional operator with solution, say, u0, which can be obtained easily. We may choose a convex homotopy H ðu; pÞ ¼ ð1  pÞF ðuÞ þ pLðuÞ ¼ 0; *

Corresponding author. Tel.: +98 821 420 7982; fax: +98 218 791 094. E-mail address: [email protected] (A. Golbabai).

0960-0779/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2006.10.037

ð2Þ

A. Golbabai, B. Keramati / Chaos, Solitons and Fractals 37 (2008) 1528–1537

1529

and continuously trace an implicitly defined curve from a starting point H(u0, 0) to a solution H(c, 1). The embedding parameter p monotonically increases from 0 to 1 as the trivial problem F(u) = 0 continuously deformed to the original problem L(u) = 0. The parameter p can be considered as an expanding parameter [9,11,13]. In fact HPM uses the homotopy parameter p as an expanding parameter [10,13] to obtain u ¼ u0 þ pu1 þ p2 u2 þ    ;

ð3Þ

when p ! 1, (3) corresponds to (2) and gives an approximation to the solution of (1) as follows: c ¼ lim þu2 þ    :

ð4Þ

p!1

The series (4) converges in most cases, and the rate of convergence [12] depends on L(u). Taking F(u) = u(x)  f(x), and substituting (3) in (2) and equating the terms with identical power of p, we obtain p0 : u0  f ðxÞ ¼ 0 ) u0 ¼ f ðxÞ; Z b p1 : u1  kðx; tÞu0 ðtÞdt ¼ 0; a Z b u1 ¼ kðx; tÞu0 ðtÞdt; a

.. . and in general we have u0 ðxÞ ¼ f ðxÞ; Z b kðx; tÞun ðtÞdt; unþ1 ðxÞ ¼

n ¼ 1; 2; . . . ;

a

which is the standard Adomian’s decomposition method. Now consider Fredholm integral equation of the first kind Z b f ðxÞ ¼ kðx; tÞcðtÞdt; a

and the convex homotopy: H ðu; pÞ ¼ ð1  pÞuðxÞ þ p

Z

b

 kðx; tÞuðtÞdt  f ðxÞ ¼ 0;

ð5Þ

a

with the starting point H(0, 0) and the solution H(c, 1), by similar operations as above, we obtain u0 ¼ 0; unþ1

u1 ¼ f ðxÞ; Z b ¼ un  kðx; tÞun ðtÞdt; a

or unþ1 ¼ ðI  KÞun ðtÞ; ¼ ðI  KÞn u1 ðtÞ; where the operator K is defined as Z b kðx; tÞðÞdt; K¼ a

and for convergence we should have kI  Kk < 1: In this paper we introduce an accelerating parameter based on HPM which causes a rapid convergence and gives the exact solution.

2. The new modified HPM In this section we propose a scheme to Paccelerate the rate of convergence of HPM applied to linear Fredholm integral equations with kernels of the form Ni¼1 gi ðxÞhi ðtÞ. We define a new convex homotopy perturbation as follows:

1530

A. Golbabai, B. Keramati / Chaos, Solitons and Fractals 37 (2008) 1528–1537

H ðu; p; mÞ ¼ ð1  pÞF ðuÞ þ pLðuÞ þ pð1  pÞ

" N X

# mi gi ðxÞ ¼ 0;

ð6Þ

i¼1

where m = [mi] and mi, i = 1, 2, . . ., N are called the accelerating parameters, and for mi = 0, i = 1, 2, . . ., N, we define H ðu; p; 0Þ ¼ Hðu; pÞ; which is the standard HPM. It is easy to verify that the solutions H(u0, 0) and H(c, 1) also satisfy (6). 2.1. Application of H(u, p, m) to second kind FIE We first assume that k(x, t) = g(x)h(t), so for equation Z b kðx; tÞcðtÞdt; c 6 x 6 d; cðxÞ ¼ f ðxÞ þ a

we consider (6) as follows: H ðu; p; mÞ ¼ ð1  pÞF ðuÞ þ pLðuÞ þ pð1  pÞ½mgðxÞ ¼ 0; where F ðuÞ ¼ uðxÞ  f ðxÞ and

LðuÞ ¼ uðxÞ  f ðxÞ 

Z

b

gðxÞhðtÞuðtÞdt ¼ 0;

a

hence we can write

 Z ð1  pÞðu  f Þ þ p u  f 

b

 gðxÞhðtÞuðtÞdt þ mpð1  pÞgðxÞ ¼ 0;

a

or u  f  pgðxÞ

Z

b

hðtÞuðtÞdt þ mpgðxÞ  mp2 gðxÞ ¼ 0;

ð7Þ

a

now by using (3) in (7), and equating the terms with identical power of p, we obtain p0 :u0  f ðxÞ ¼ 0 ) u0 ¼ f ðxÞ; Z b hðtÞu0 ðtÞdt ¼ 0; p1 :u1 þ mgðxÞ  gðxÞ a Z b hðtÞf ðtÞdt; u1 ¼ ðc  mÞgðxÞ; where c ¼ a Z b p2 :u2  mgðxÞ  gðxÞ hðtÞu1 ðtÞdt ¼ 0; a Z b u2 ¼ ½m þ ðc  mÞagðxÞ; where a ¼ hðtÞgðtÞdt; a Z b hðtÞu2 ðtÞdt; p3 :u3 ¼ a

and in general Z unþ1 ¼

b

hðtÞun ðtÞdt;

n ¼ 2; 3; . . . ;

a

now we find m such that u2 = 0, since if u2 = 0 then u3 = u4 =    = 0, and the exact solution will be obtained as u(x) = u0(x) + u1(x), hence for all values of x we should have m þ ðc  mÞa ¼ 0; or m¼

ca ¼ a1

hR b a

ihR i b hðtÞf ðtÞdt a hðtÞgðtÞdt Rb hðtÞgðtÞdt  1 a

A. Golbabai, B. Keramati / Chaos, Solitons and Fractals 37 (2008) 1528–1537

1531

or Rb m ¼ Rb a

a

Z

kðt; tÞdt

kðt; tÞdt  1

b

hðtÞf ðtÞdt;

ð8Þ

a

provided that Z b kðt; tÞdt 6¼ 1: a

Now consider the general case kðx; tÞ ¼

N X

gi ðxÞhi ðtÞ;

i¼1

here we choose the convex homotopy as follows:

"

H ðu; p; mÞ ¼ ð1  pÞF ðuÞ þ pLðuÞ þ pð1  pÞ

N X

# mi gi ðxÞ ¼ 0;

i¼1

by doing similar manipulations, we obtain u0 ¼ f ; u1 ¼ 

N X

mi gi ðxÞ þ

i¼1

¼

i¼1

u2 ¼

b

mi gi ðxÞ þ

N X

Z N  X mi þ

gi ðxÞ

b

gi ðxÞ

Z

b

hi ðtÞu1 ðtÞdt;

a

 hi ðtÞu1 ðtÞdt gi ðxÞ

a

i¼1

unþ1 ¼

hi ðtÞu0 ðtÞdt;

a

a

i¼1

.. . N X

b

 hi ðtÞf ðtÞdt  mi gi ðxÞ

i¼1

¼

gi ðxÞ

Z

i¼1

N Z X

N X

N X

Z

b

hi ðtÞun ðtÞdt;

n ¼ 2; 3; . . . ;

a

i¼1

we try to find the parameters mi,i = 1, 2, . . ., N, such that u2 ¼ u3 ¼ u4 ¼    ¼ 0; hence from u2 we should have Z b mi þ hi ðtÞu1 ðtÞdt ¼ 0;

8 x 2 ½c; d;

ð9Þ

a

now by substituting u1 in (9), we obtain Z b  N Z b X hi ðtÞ hj ðtÞf ðtÞdt  mj gj ðtÞdt ¼ 0; mi þ j¼1

a

a

let cj ¼

Z

b

hj ðtÞf ðtÞdt;

bij ¼

a

Z

b

hi ðtÞgj ðtÞdt;

a

then mi þ

N X j¼1

bij ðcj  mj Þ ¼ 0;

i ¼ 1; 2; . . . N;

ð10Þ

1532

A. Golbabai, B. Keramati / Chaos, Solitons and Fractals 37 (2008) 1528–1537

under certain condition,the values of mi, i = 1, 2, . . ., N, can be obtained from the system of linear equations in (10). Let the matrix B and the vectors m and c be defined as follows: B ¼ ½bij ;

m ¼ ½mj ;

c ¼ ½cj ;

therefore from (10) we can write ðB  IÞm ¼ Bc; and if (B  I) is nonsingular then m ¼ ðB  IÞ1 Bc:

ð11Þ

In the case of non-degenerate kernels,by using Taylor expansion for functions of two variables, we can write k(x, t) (if possible) as follows: kðx; tÞ ¼

N X

gi ðxÞhi ðtÞ;

i¼1

and by applying H(u, p, m) method we can approximate the solution of the given integral equation. 2.2. Application of H(u,p,m) to first kind FIE Consider the integral equation: Z b kðx; tÞcðtÞdt; d 6 x 6 g; f ðxÞ ¼

ð12Þ

a

let k(x, t) = g(x)h(t), then from (12): f ðxÞ ¼ cgðxÞ; using (6) we can write ð1  pÞF ðuÞ þ pLðuÞ þ pð1  pÞ½mgðxÞ ¼ 0; where F ðuÞ ¼ u;

LðuÞ ¼

Z

b

gðxÞhðtÞdt  f ;

a

by similar operations in Section 2.1 we obtain p0 :u0 ¼ 0; p1 :u1 ¼ ðc  mÞgðxÞ; 2

p :u2 ¼ u1 þ mgðxÞ  gðxÞ

Z

b

hðtÞu1 ðtÞdt; Z where a ¼ a

¼ ½c  ðc  mÞagðxÞ; Z b hðtÞu2 ðtÞdt; u3 ¼ u2  gðxÞ

b

hðtÞgðtÞdt ¼ a

Z

b

kðt; tÞdt; a

a

and in general unþ1 ¼ un  gðxÞ

Z

b

hðtÞun ðtÞdt;

n ¼ 2; 3; . . . ;

a

for u2 = u3 = u4    = 0, from u2 we should have c  ðc  mÞa ¼ 0;

8x 2 ½c; d;

and consequently c m ¼ c  ; ða 6¼ 0Þ: a P Now assume kðx; tÞ ¼ Ni¼1 gi ðxÞhi ðtÞ, then it is easy to verify that

ð13Þ

A. Golbabai, B. Keramati / Chaos, Solitons and Fractals 37 (2008) 1528–1537

f ðxÞ ¼

N X

1533

ci gi ðxÞ;

i¼1

in this case as in Section 2.1, we choose the following convex homotopy "Z # " # N N b X X ð1  pÞu þ p gi ðxÞhi ðtÞdt  f þ pð1  pÞ mi gi ðxÞ ¼ 0; a

i¼1

ð14Þ

i¼1

by substituting (3) in (14) and equating the coefficients of the identical powers of p, we obtain u0 ¼ 0; N X ½ci  mi gi ðxÞ; u1 ¼ i¼1

u2 ¼

N X

" ci 

# N X ðcj  mj Þbij gi ðxÞ;

i¼1

where bij ¼

Z

b

hi ðtÞgj ðtÞdt a

j¼1

.. . unþ1 ¼ un 

N X

gi ðxÞ

Z

i¼1

b

hi ðtÞun ðtÞdt;

n ¼ 2; 3; . . . ;

a

as above for u2 = u3 = u4    = 0, parameters mi, i = 1, 2, . . ., N, should satisfy the following linear system of equations ci 

N X ðcj  mj Þbij ¼ 0;

i ¼ 1; 2; . . . ; N;

j¼1

or in matrix form we have Bm ¼ ðB  IÞc; where B ¼ ½bij ;

m ¼ ½mj ;

c ¼ ½cj ;

j ¼ 1; 2; . . . N ;

or m ¼ B1 ðB  IÞc;

ð15Þ

provided that (B) is not singular. For non-degenerate kernels k(x, t), using Taylor expansion, we can write kðx; tÞ ¼

N X

gi ðxÞhi ðtÞ;

f ðxÞ ¼

i¼1

N X

ci gi ðxÞ;

i¼1

and then by applying H(u, p, m) method, we can approximate the solution of the given problem.

3. Numerical comparisons

Example 1. Consider the equation Z 1  1 xtuðtÞdt; uðxÞ ¼ e3x  2e3 þ 1 x þ 9 0 with exact solution u(x) = e3x. We apply H(u, p) and H(u, p, m) methods to approximate the solution as follows H(u, p) method: We have u0 ðxÞ ¼ f ðxÞ; Z b kðx; tÞun ðtÞdt; unþ1 ðxÞ ¼ a

n ¼ 1; 2; . . . ;

1534

A. Golbabai, B. Keramati / Chaos, Solitons and Fractals 37 (2008) 1528–1537

so we obtain Z u1 ¼

1

xtu0 ðtÞdt ¼ x

0

Z

1 0

  1 t e3x  ð2e3 þ 1Þt dt; 9

  1 1 þ ð4e3  1Þ x; ¼ 9 27   1 1 1 þ ð4e3  1Þ x; u2 ¼ 3 9 27   1 1 1 u3 ¼ 2 þ ð4e3  1Þ x; 3 9 27   1 1 1 u4 ¼ 3 þ ð4e3  1Þ x; 3 9 27 .. .

Hence the solution will be as follows: uðxÞ ¼ u0 þ u1 þ u2 þ    ;

   1 1 1 1 1 1 ¼ e3x  ð2e3 þ 1Þx þ 1 þ þ 2 þ 3 þ    þ ð4e3  1Þ x ¼ e3x : 9 3 3 9 27 3

H(u, p, m) method: In this case we have: u0 ¼ f ðxÞ; u1 ¼ ðc  mÞgðxÞ;

where c ¼

Z

b

hðtÞf ðtÞdt;

a

and uðxÞ ¼ u0 þ u1 : For this problem g(x) = x and h(t) = t, hence 1 u0 ¼ e3x  ð2e3 þ 1Þx; 9 Z 1  1 u1 ¼ tðe3t  ð2e3 þ 1ÞtÞdt  m x; 9 0   4 3 2 e þ  m x; ¼ 27 27 and from (8): R1 m ¼ R1 0

0

t2 dt

t2 dt  1

Z

1 0

    1 1 4 3 2 2 3 1 t e3t  ð2e3 þ 1Þt dt ¼ 1 3 e þ ¼ e þ ; 9 27 27 27 27  1 3

hence we obtain    4 3 2 1 e þ þ x ¼ ð2e3 þ 1Þx; u1 ¼ 27 27 9 and the solution will be as follows: uðxÞ ¼ u0 þ u1 ¼ e3x : As we observe, after two terms the exact solution is obtained. Example 2. Approximate the solution of Z p ðcos xÞðcos tÞuðtÞdt: uðxÞ ¼ ð1  2pÞ cos x þ sin x þ 4 0

The exact solution is u(x) = cos x + sin x.

A. Golbabai, B. Keramati / Chaos, Solitons and Fractals 37 (2008) 1528–1537

1535

H(u, p) method: u0 ¼ ð1  2pÞ cos x þ sin x; u1 ¼ ð2pÞð1  2pÞ cos x; u2 ¼ ð2pÞ2 ð1  2pÞ cos x; .. . unþ1 ¼ ð2pÞn ð1  2pÞ cos x; hence uðxÞ ¼ ½1 þ ð2pÞ þ ð2pÞ2 þ ð2pÞ3 þ   ð1  2pÞ cos x þ sin x; as we observe the series is diverging and no convergent solution will be obtained by this method. H(u, p, m) method: ca gðxÞ ¼ cos x; hðtÞ ¼ 4 cos t; m ¼ ; a 1 Z p Z p c¼4 cos t½ð1  2pÞ cos tdt ¼ 2pð1  2pÞ; a ¼ 4 ðcos tÞ2 dt ¼ 2p; 0

0

hence m ¼ 4p2

) u1 ¼ ðc  mÞgðxÞ ¼ 2p cos x;

and the solution is as follows: uðxÞ ¼ u0 þ u1 ¼ cos x þ sin x: Example 3. Consider the equation Z 0

1 2

1 exþt uðtÞdt ¼ ðe  1Þex ; 2

1 06x6 ; 2

where the true solution is u(x) = ex. By using H(u, p) method, we obtain u0 ¼ 0; 1 u1 ¼ ðe  1Þex ; 2   1 e3 x ðe  1Þ e; u2 ¼ 2 2  2 1 e3 u3 ¼ ðe  1Þ ex ; 2 2 .. . hence uðxÞ ¼ u0 þ u1 þ u2 þ u4 þ    ; " #    2  3 1 e3 e3 e3 þ ¼ ðe  1Þ 1   þ    ex ; 2 2 2 2 1 1 ex ¼ ðe  1Þ 2 1 þ ðe3 Þ 2 ¼ ex ; in fact we could find the solution as a limit of (16). Now applying H(u, p, m), we can easily find the solution as following: gðxÞ ¼ ex ;

hðtÞ ¼ et ;

1 c ¼ ðe  1Þ; 2

ð16Þ

1536

A. Golbabai, B. Keramati / Chaos, Solitons and Fractals 37 (2008) 1528–1537

and from (13), we have

Z 12 c 1 e2t dt ¼ ðe  1Þ; m¼c ; a¼ a 2 0 so we obtain m ¼ 12 e  32, and consequently   1 1 3 u1 ¼ ðe  1Þ  e þ ex ¼ ex ; 2 2 2 and the solution is uðxÞ ¼ u0 þ u1 ¼ ex : Example 4. For integral equation Z

1 0

1 1 ðtex þ 1ÞuðtÞdt ¼ ex þ ; 3 2

ð17Þ

where the exact solution is u(x) = x, using H(u, p, m) method, we have 1 1 c1 ¼ ; c2 ¼ ; g1 ðxÞ ¼ ex ; g2 ðxÞ ¼ 1; h1 ðtÞ ¼ t; h2 ðtÞ ¼ 1; 3 2 by applying (15) and information of the problem, we obtain   1 1 2 B¼ e1 1 and 

1 1 2 e1 1



m1 m2





1 0 2 ¼ e1 0

"

" 1 # 3 1 2

¼

1 4 1 ðe 3

 1Þ

# ;

so the values of m1 and m2 are obtained as follows: m1 ¼

1 4

 16 ðe  1Þ ; 1  12 ðe  1Þ

m2 ¼

1

1 ðe  1Þ 6 ;  12 ðe  1Þ

consequently u1 ¼ ðc1  m1 Þg1 ðxÞ þ ðc2  m2 Þg2 ðxÞ; 1 6  4ðe  1Þ ¼ ex þ ; 12  6ðe  1Þ 12  6ðe  1Þ and finally uðxÞ ¼ u0 þ u1 ¼

1 6  4ðe  1Þ ex þ ; 12  6ðe  1Þ 12  6ðe  1Þ

which is in fact another solution of (17). It is easy to show that H(u, p) method can not give a convergent solution. 4. Conclusion In this paper we made some modifications to the homotopy perturbation method by introducing accelerating parameters for solving linear Fredholm integral equations. In contrast to the standard one, H(u, p, m) method is a simple and very effective tool for calculating the exact solutions and in most cases it can give solution to some problems which could not be solved by H(u, p) method.

References [1] Abbasbandy S. Application of He’s homotopy perturbation method to functional integral equations. Chaos, Solitons & Fractals 2007;31(5):1243–7.

A. Golbabai, B. Keramati / Chaos, Solitons and Fractals 37 (2008) 1528–1537

1537

[2] Abbasbandy S. Modified homotopy perturbation method for nonlinear equations and comparison with Adomian decomposition method. Appl Math Comput 2006;172:431–8. [3] Abbasbandy S. Numerical solution of integral equation: Homotopy perturbation method and Adomian’s decomposition method. Appl Math Comput 2006;173:493–500. [4] Chun C. Integration using He’s homotopy perturbation method. Chaos, Solitons & Fractals 2007;34(4):1130–4. [5] Cveticanin L. Homotopy perturbation method for pure nonlinear differential equation. Chaos, Solitons & Fractals 2006;30(5):1221–30. [6] He J-H. A coupling method of homotopy technique and perturbation technique for nonlinear problems. Int Nonlinear Mech 2000;35(1):37–43. [7] He J-H. The homotopy perturbation method for nonlinear oscillators with discontinuities. Appl Math Comput 2004;151:287–92. [8] He J-H. Comparison of homotopy perturbation method and homotopy analysis method. Appl Math Comput 2004;156:527–39. [9] He J-H. Homotopy perturbation technique. Comput Methods Appl Mech Eng 1999;173(3–4):257–62. [10] Hillermeier C. Generalized homotopy approach to multiobjective optimization. Int J Optim Theory Appl 2001;110(3):557–83. [11] Liao SJ. An approximate solution technique not depending on small parameters: a special example. Int J Non-Linear Mech 1995;30(3):371–80. [12] Liao SJ. Boundary element method for general nonlinear differential operators. Eng Anal Boundary Element 1997;20(2):91–9. [13] Nayef AH. Problems in perturbation. New york: John Wiley; 1985. [14] Siddiqui AM, Mahmood R, Ghori QK. Homotopy perturbation method for thin film flow of a third fluid down an inclined plane. Chaos, Solitons & Fractals 2008;35(1):140–7. [15] Ozis T, Yildirim A. A note on He’s homotopy perturbation method for van der Pol oscillator with very strong nonlinearity. Chaos, Solitons & Fractals 2007;34(3):989–91.