Applied Mathematics Letters 24 (2011) 1634–1639
Contents lists available at ScienceDirect
Applied Mathematics Letters journal homepage: www.elsevier.com/locate/aml
Approximate analytical solutions of reaction–diffusion equations with exponential source term: Homotopy perturbation method (HPM) S.O. Ajadi ∗ , M. Zuilino Department of Mathematics, Ben-Gurion University of Negev, PO Box 653 Beer-Sheva, 84105, Israel
article
info
Article history: Received 14 September 2009 Received in revised form 1 April 2011 Accepted 5 April 2011 Keywords: HPM Exponential source term Reaction–diffusion Combustion
abstract In this letter, the solutions of some nonlinear differential equations have been obtained by means of the homotopy perturbation method (HPM). Applications of the homotopy method to some nonlinear reaction–diffusion equations with exponential source term show rapid convergence of the sequence constructed by this method to the exact solutions. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction In reality, closed-form solutions for nonlinear problems are difficult to come-by and this has resulted in the development of numerical and other approximate solutions [1]. Recently, the application of the homotopy perturbation method (HPM) to nonlinear initial boundary value problems shows that there is a rapid convergence of the sequence constructed by this method to the exact solutions [2,3]. In this paper, we consider the initial value problem (IVP) dθ dx
+ δ eθ = 0,
θ (0) = 0,
(1) (2)
which describes the temperature equation in a pressure driven porous media combustion with weak internal thermal diffusion [4]. In the absence of reactant consumption, the general heat balance equation for the one-step reaction system can be written as [1]
∂θ = 1θ + δ eθ = 0, in Ω ∂t ∂θ + Bi θ = 0, ∂ Ω , ∂n
(3) (4)
where 1 is an operator and δ is related to the characteristic chemical reaction, while ∂∂n is the outward normal derivative on the boundary and Bi is the Biot number, a parameter which determines whether or not the temperatures inside a body will vary significantly in space.
∗
Corresponding address: Department of Mathematics, Obafemi Awolowo University, Ile Ife, 220005, Nigeria. Tel.: +234 8030545990. E-mail addresses:
[email protected],
[email protected] (S.O. Ajadi).
0893-9659/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.aml.2011.04.003
S.O. Ajadi, M. Zuilino / Applied Mathematics Letters 24 (2011) 1634–1639
1635
2. Homotopy perturbation technique In line with [5,6,3,7–11], we illustrate the homotopy perturbation method, we consider the nonlinear equation A(u) − f (r ) = 0,
r ∈ Ω,
(5)
with the boundary conditions:
B u,
∂u ∂n
= 0,
(r ∈ ∂ Ω ),
(6)
where A is a general differential operator, B is a boundary operator, f (r ) is a known analytic function, and ∂ Ω is the boundary of domain Ω . The operator A is generally divided into two parts; L and N, where L and N are linear and nonlinear parts of A, respectively. Therefore, (3) may be written as L(u) + N (u) − f (r ) = 0.
(7)
We construct a homotopy v(r , p) : Ω × [0, 1] → ℜ which satisfies H (v, p) = [L(v) − L(u0 )] + p[L(u0 )] + p[N (v) − f (r )] = 0,
(8)
H (v, p) = [L(u) − L(u0 )] + p[A(u) − f (r )] = 0
(9)
or
where p ∈ [0, 1] is called the homotopy parameter and u0 is an initial approximation of (5) which satisfies the specified boundary conditions. When p = 0 or p = 1, we have H (v, 0) = L(u) − L(u0 ) = 0,
H (v, 1) = A(u) − f (r ) = 0.
(10)
On the other hand, if p ∈ (0, 1), then the homotopy H (v, p) deforms from L(u) − L(u0 ) to A(u) − f (r ). Thus, the solution of (5)–(6) may be expressed as
v = v0 + pv1 + p2 v2 + p3 v3 + · · · .
(11)
Eventually, at p = 1, the system takes the original form of the equation and the final stage of deformation gives the desired solution. Thus taking limit u = lim v = v0 + v1 + v2 + . p→1
(12)
3. Some particular examples 3.1. Example I We start by applying the homotopy technique to an initial value problem of the form dθ dx
+ δ eθ = 0.
(13)
θ (0) = 0.
(14)
In line with [1], we define homotopy as dθ dx
−
dy0 dx
+p
dy0 dx
+ δe
θ
= 0.
(15)
Suppose that the solution of (13)–(14) takes the form
θ = v0 + pv1 + p2 v2 + p3 v3 + p4 v4 + p5 v5 · · · ,
(16)
with an initial approximation
v0 (x) = y0 (x) = c ,
(17)
where c is to be determined. Eq. (15) may be expressed as [7] dθ dx
−
dy0 dx
[ +p
] θ2 θ3 + δ ev0 1 + θ1 + 1 + 1 · · · = 0, dx 2! 3!
dy0
(18)
1636
S.O. Ajadi, M. Zuilino / Applied Mathematics Letters 24 (2011) 1634–1639
where θ1 = pv1 + p2 v2 + p3 v3 + p4 v4 + p5 v5 · · ·. After substituting (16) into (18), and collecting terms in powers of p, we obtain dv 1 dx dv 2
+
dy0 dx
+ δ ev0 = 0,
(19)
+ δ ev0 v1 = 0, v2 dv 3 + δ ev0 v2 + 1 = 0,
(20)
dx
(21)
dx
2
dv 4
v3 + δ ev0 v3 + v1 v2 + 1 = 0,
dx
(22)
6
dv 5 dx
v 2 v2 v2 v4 + δ ev0 v4 + v1 v3 + 1 + 2 + 1 = 0, 2
2
(23)
24
v1 (0) = v2 (0) = v3 (0) = v4 (0) = v5 (0) · · · = 0.
(24)
The solutions of (19)–(23) are
v1 = −δ ec x,
v2 =
δ2 2
e2c x2 ,
v3 = −
δ3 3
e3c x3 ,
v4 =
δ4 4
e4c x4 ,
and
v5 = −
δ5 5
e5c x5 .
(25)
The fifth order approximation is given by
θ = v0 + v1 + v2 + v3 + v4 + v5 = c − δ ec x +
δ2 2
e2c x2 −
δ3 3
e3c x3 +
δ4 4
e4c x4 −
δ5 5
e5c x5 .
(26)
By using (14) on (26), we obtain c = 0, and the solution corresponds to the exact solution
θ=
∞ − (−1)n δ n xn = − ln(1 + δ x). n n! n =1
(27)
3.2. Example II The steady state formulations of (3)–(4) due to Frank-Kamenetskii have been developed and studied by many others [12,1,13,14] d2 θ dx2 dθ dx
+
j dθ x dx
+ δ x−β eθ = 0,
(0) = 0,
θ (1) = 0,
(28) (29)
where j is related to the geometry and β is a numerical exponent. Thus after substituting (16) into (28)–(29) and collecting terms in p, d2 v 1 dx2 d2 v 2 dx2 d2 v 3 dx2 d2 v 4 dx2 d2 v 5 dx2
+ + + + +
j dv1 x dx j dv2 x dx j dv3 x dx j dv4 x dx j dv5 x dx
+
d2 y0 dx2
+ δ x−β ev0 = 0,
+ δ x−β ev0 v1 = 0, v2 + δ x−β ev0 v2 + 1 = 0, 2! v13 −β v0 + δx e v3 + v1 v2 + = 0, 3! v2 v 2 v2 v4 + δ x−β ev0 v4 + v1 v3 + 2 + 1 + 1 = 0. 2! 2! 4!
(30) (31) (32)
(33)
(34)
S.O. Ajadi, M. Zuilino / Applied Mathematics Letters 24 (2011) 1634–1639
1637
The solutions of (30)–(34) are
δ 2 e2c x4−2β δ ec x2−β , v2 = , 2 (2 − β)(j + 1 − β) 2(2 − β) (j + 1 − β)(j + 3 − 2β) δ 3 e3c (2j + 4 − 3β) v3 = − x6−3β , 2 2(2 − β) (j + 1 − β)2 (j + 3 − 2β)(j + 5 − 3β) v0 = c ,
v1 = −
δ 4 e4c (3j2 + 16j − 11jβ − 25β + 17 + 9β 2 ) x8−4β , 12(2 − β)4 (j + 1 − β)3 (j + 3 − 2β)(j + 5 − 3β)(j + 7 − 4β) 4 24j + 316j3 − 206j2 β + 609j2 β 2 − 1890j2 β −δ 5 e5c +1416j2 − 745jβ 3 + 3468jβ 2 − 5190jβ + 2516j x10−5β +324β 4 − 1991β 3 + 4419β 2 − 4234β + 1488 , v5 = 5 120(2 − β) (j + 1 − β)4 (j + 3 − 2β)2 (j + 5 − 3β)(j + 7 − 4β)(j + 9 − 5β) v4 =
(35)
where β ̸= 2, which is also satisfied by the corresponding closed-form solution [14]. The value of c is determined by using the Dirichlet condition in the boundary conditions (29). For example, when j = β = 0 and δ = 0.86, c = 0.8363. 3.3. Example III An unsteady 1D form of (3)–(4) is
∂θ ∂ 2θ = 2 + δ eθ = 0, (x, t ) ∈ Ω ⊂ ℜ2 , ∂t ∂x θ (x, 0) = 1/2x(1 − x), dθ (0, t ) = 0, θ (1, t ) = 0.
(36) (37) (38)
dx
In line with [3], we construct homotopy in the form
2 ∂θ ∂v0 ∂ θ ∂v0 θ − =p + e − , ∂t ∂t ∂ x2 ∂t
(x, t ) ∈ Ω ⊂ ℜ2 .
(39)
After substituting Eq. (16) into Eq. (39) and collecting terms in p,
∂v1 ∂ 2 v0 dv0 = − + ev 0 , ∂t ∂ x2 dt ∂v2 ∂ 2 v1 = + ev0 v1 , ∂t ∂ x2 ∂v3 ∂ 2 v2 v12 v0 = + e v + , 2 ∂t ∂ x2 2! ∂v4 ∂ 2 v3 v0 +e = v3 + v1 v2 + ∂t ∂ x2 ∂v5 ∂ 2 v4 v0 + e = v4 + v1 v3 + ∂t ∂ x2
(40) (41) (42)
v13 3!
,
(43)
v22 v 2 v2 v4 + 1 + 1 2! 2! 4!
.
(44)
In line with [1], it is convenient to choose v0 = 1/2x(1 − x), then
1 v1 = e 2 x(1−x) − 1 t , v2 = v3 =
t2
e
2! t3 3!
1 x(1−x) 2
+
(45) 1 2
5 − 6e
1 x(1−x) 2
2 −x
+ 5e
1
− 2 e 2 x(1−x) ,
1 x(1−x) 2
1
−7
2
(46)
2 −x
x(1−x)
+ 2e
+
1 2
4 −x
1
e 2 x(1−x) ,
(47)
1638
S.O. Ajadi, M. Zuilino / Applied Mathematics Letters 24 (2011) 1634–1639
Fig. 1. Plot of θ(x) against x (HPM).
Fig. 2. Plot of θ(x) against x (Exact).
Fig. 3. 3D plot of θ(x, t ) against x (HPM).
v4 =
t4
4!
−20 + 6e
1
3 x(1−x) 2
2
x(1−x)
− 20e
1 x(1−x) 2
+ 31e
1 x(1−x) 2
4
1
+ 34
1 x(1−x) 2
1 2
2 −x
1
ex(1−x)
2
− 59 −x e + 21 −x e + 52 −x 2 2 2 4 6 1 1 1 1 1 −8 −x + − x − 20 − x e 2 x(1−x) e 2 x(1−x) . 2
2
(48)
2
] 30 ∗ 2 ∂ 2 v4∗ 15 ∗ 2 ∗ 1 ∗4 1 1 ∗ ∗ x(1−x) ∗ 2 v5 = +e v4 + 4v1 v3 + v2 + v1 v2 + v1 e 2 x(1−x) , 5! ∂ x2 2! 2! 4! t5
[
(49)
where
v1∗ = v1 /t ,
v2∗ =
2! t2
e
−1 x(1−x) 2
v2 ,
v3∗ =
3! t3
3!e
−1 x(1−x) 2
v3 ,
v4∗ =
4! t4
3!e
−1 x(1−x) 2
v4 .
(50)
4. Conclusion It is observed that Fig. 1. obtained by the fifth order approximation (HPM) converges to the profiles given by the exact solution (Fig. 2). Similarly, in the unsteady state, Figs. 3 and 4 show that the HPM solution is reliable as it confirms the numerical solutions obtained by NDSolve in the MATHEMATICA package. Although the approximate homotopy perturbation
S.O. Ajadi, M. Zuilino / Applied Mathematics Letters 24 (2011) 1634–1639
1639
Fig. 4. 3D plot of θ(x, t ) against x (NDSolve).
method is found to work extremely well in the examples considered, the approach may be less effective and accurate in the presence of more complicated nonlinear source terms. References [1] E. Balakrishnan, A. Swift, G.C. Wake, Critical values for some non-class a geometries in thermal ignition theory, Mathematical and Computer Modelling 24 (1996) 1–10. [2] D. Lesnic, The decomposition method for Cauchy reaction–diffusion problems, Applied Mathematics Letters 20 (2007) 412–418. [3] A. Yildirim, Application of He’s homotopy pertubation method for solving the Cauchy reaction–diffusion problem, Computers & Mathematics with Applications 57 (2009) 612–618. [4] V. Bykov, Mathematical problems of multiphase combustion, Doctoral Thesis, Ben-Gurion Univerisity of Negev, 2004, pp. 63–64. [5] J.H. He, Homotopy perturbation method: a new nonlinear analytical technique, Applied Mathematics and Computation 135 (2003) 73–79. [6] J.H. He, Homotopy perturbation method for solving boundary value problems, Physics Letters A 350 (2006) 87–88. [7] A. Yildirim, H. Kocak, Homotopy pertubation method for solving the space–time fractional advection dispersion equation, Advances in Water Resources 32 (2009) 1711–1716. [8] I. Ates, A. Yildirim, Application variational iteration method to fractional initial-value problems, International Journal of Nonlinear Sciences and Numerical Simulation 10 (2009) 877–883. [9] A. Konuralp, C. Konuralp, A. Yildirim, Numerical solution to Van Der Pol equation with fractional damping, Physica Scripta T 136 (2009) 014034. 5pp. [10] A. Yildirim, Y. Gulkanat, Analytical approach to fractional Zakharov–Kuznetsov equations by He’s homotopy pertubation method, Communications in Theoretical Physics 53 (2010) 1005–1010. [11] S. Momani, A. Yildirim, Analytical approximate solutions of the fractional convection–diffusion equation with nonlinear source term by He’s homotopy pertubation method, International Journal of Computer Mathematics 57 (2010) 1055–1065. [12] S.O. Ajadi, A note on the thermal stability of a reactive non-Newtonian flow in a cylindrical pipe, International Communications in Heat and Mass Transfer 36 (2009) 63–68. [13] S.S. Okoya, S.O. Ajadi, A.B. Kolawole, Thermal runaway for a reaction–diffusion equation, International Communications in Heat and Mass Transfer 30 (2003) 845–850. [14] S.S. Okoya, S.O. Ajadi, Critical parameters for thermal conduction equation, Journal of Mechanics Research Communications 26 (1999) 363–370.