The combined Laplace transform and new homotopy perturbation methods for stiff systems of ODEs

Applied Mathematical Modelling 36 (2012) 3638–3644

Contents lists available at SciVerse ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

The combined Laplace transform and new homotopy perturbation methods for stiff systems of ODEs Hossein Aminikhah ⇑ Department of Applied Mathematics, School of Mathematical Sciences, University of Guilan, P.O. Box 1914, Rasht, Iran

a r t i c l e

i n f o

Article history: Received 16 July 2011 Received in revised form 30 September 2011 Accepted 11 October 2011 Available online 6 November 2011

a b s t r a c t In this paper, we propose new technique for solving stiff system of ordinary differential equations. This algorithm is based on Laplace transform and homotopy perturbation methods. The new technique is applied to solving two mathematical models of stiff problem. We show that the present approach is relatively easy, efficient and highly accurate. Ó 2011 Elsevier Inc. All rights reserved.

Keywords: Laplace transform method New homotopy perturbation method Stiff system of ordinary differential equations

1. Introduction Stiff differential equations arise in fluid mechanics, elasticity, electrical networks, chemical reactions, and many other areas of physical importance [1–5]. Computer simulation of dynamic systems very often leads to the solution of a set of stiff ordinary differential equations. The solution of this set of equations involves the eigenvalues of its Jacobian matrix. The greater the spread in eigenvalues, the more time consuming the solutions become when existing numerical methods are employed. Stiff differential equations can become a very serious problem for some systems, rendering accurate numerical solutions completely uneconomic. Several schemes have been developed for the numerical solution of differential equations. The homotopy perturbation method was proposed by He [6] in 1999. This method has been used by many mathematicians and engineers to solve various functional equations. Homotopy method was further developed and improved by He and applied to nonlinear oscillators with discontinuities [7], nonlinear wave equations [8], and boundary value problem [9]. It can be said that He’s HPM is a universal one and is able to solve various kinds of nonlinear functional equations. For example, it was applied to nonlinear Schrödinger equations [10], to nonlinear equations arising in heat transfer [11], to the quadratic Riccati differential equation [12], and to other equations [13–17]. In this method, the solution is considered as the summation of an infinite series which usually converges rapidly to exact solutions. In this paper we introduce a new form of homotopy perturbation and Laplace transform methods that have good stability properties. The new method is an efficient computational algorithm for linear and nonlinear stiff system of ordinary differential equations. 2. Analysis of the method To illustrate the basic ideas of this method, let us consider the following system of nonlinear differential equation

⇑ Tel./fax: +98 1313233509. E-mail addresses: [email protected], [email protected] 0307-904X/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2011.10.014

H. Aminikhah / Applied Mathematical Modelling 36 (2012) 3638–3644

AðUÞ  f ðrÞ ¼ 0;

r 2 X; U 2 Rn ;

3639

ð1Þ

with the following initial conditions

Uð0Þ ¼ a0 ;

U 0 ð0Þ ¼ a1 ; . . . ; U ðn1Þ ð0Þ ¼ an1 ;

ð2Þ

where A is a general differential operator and f(r) is a known analytical function. The operator A can be divided into two parts, L and N, where L is a linear and N is a nonlinear operator. Therefore Eq. (1) can be rewritten as

LðUÞ þ NðUÞ  f ðrÞ ¼ 0:

ð3Þ n

By the NHPM [18], we construct a homotopy Vðr; pÞ : X  ½0; 1 ! R , which satisfies

HðV; pÞ ¼ ð1  pÞ½LðVÞ  v 0  þ p½AðVÞ  f ðrÞ ¼ 0;

p 2 ½0; 1; r 2 X;

ð4Þ

or equivalently,

HðV; pÞ ¼ LðVÞ  v 0 þ pv 0 þ p½NðVÞ  f ðrÞ ¼ 0;

ð5Þ

where p 2 [0, 1] is an embedding parameter, v0 is an initial approximation of solution of Eq. (1). Clearly, we have from Eqs. (4) and (5)

HðV; 0Þ ¼ LðVÞ  v 0 ¼ 0;

ð6Þ

HðV; 1Þ ¼ AðVÞ  f ðrÞ ¼ 0:

ð7Þ

By applying Laplace transform on both sides of (5), we have

LfLðVÞ  v 0 þ pv 0 þ p½NðVÞ  f ðrÞg ¼ 0:

ð8Þ

Using the differential property of Laplace transform we have

sn LfVg  sn1 Vð0Þ  sn2 V 0 ð0Þ      V ðn1Þ ð0Þ ¼ Lfv 0  pv 0 þ p½NðVÞ  f ðrÞg;

ð9Þ

or

LfVg ¼

1 n1 fs Vð0Þ þ sn2 V 0 ð0Þ þ    þ V ðn1Þ ð0Þ þ Lfv 0  pv 0 þ p½NðVÞ  f ðrÞgg: sn

ð10Þ

By applying inverse Laplace transform on both sides of (10), we have

V ¼ L1

 n o 1 n1 ðn1Þ n2 0 s Vð0Þ þ s V ð0Þ þ    þ V ð0Þ þ Lf v  p v þ p½NðVÞ  f ðrÞg : 0 0 sn

ð11Þ

According to the HPM, we can first use the embedding parameter pas a small parameter, and assume that the solutions of Eq. (11) can be represented as a power series in p as

VðxÞ ¼

1 X

pn V n :

ð12Þ

n¼0

Now let us write Eq. (12) in the following form

( ( ( " ! #))) 1 X 1 n1 ðn1Þ n2 0 n : p Vn ¼ L s Vð0Þ þ s V ð0Þ þ    þ V ð0Þ þ L v 0  pv 0 þ p N p V n  f ðrÞ sn n¼0 n¼0

1 X

n

1

ð13Þ

Comparing coefficients of terms with identical powers of p, leads to

 1 n1 ðs Vð0Þ þ sn2 V 0 ð0Þ þ    þ V ðn1Þ ð0Þ þ Lfv 0 gÞ ; n s   1 p1 : V 1 ¼ L1 n ðLfNðV 0 Þ  v 0  f ðrÞgÞ ; s   1 1 2 ðLfNðV ; V ÞgÞ ; p : V2 ¼ L 0 1 sn   1 p3 : V 3 ¼ L1 n ðLfNðV 0 ; V 1 ; V 2 ÞgÞ ; s .. .   1 pj : V j ¼ L1 n ðLfNðV 0 ; V 1 ; V 2 ; . . . ; V j1 ÞgÞ ; s .. . p0 : V 0 ¼ L1



ð14Þ

Suppose that the initial approximation has the form V(0) = v0 = a0, V(0) = a1, . . . , V(n1)(0) = an1, therefore the exact solution may be obtained as following

3640

H. Aminikhah / Applied Mathematical Modelling 36 (2012) 3638–3644

U ¼ lim V ¼ V 0 þ V 1 þ V 2 þ   

ð15Þ

p!1

3. LTNHPM applied to stiff system of ODEs Consider an initial value problem for a large nonlinear stiff system of ordinary differential equations,

x0 ðtÞ ¼ FðxðtÞÞ þ ef ðxðtÞ; tÞ;

xð0Þ ¼ x0 ; T

ð16Þ T

where F(x) = (F1(x), F2(x), . . . , FN(x)) , x(t) = (x1(t), x2(t), . . . , xN(t)) , e is a small perturbation parameter and

f ðxðtÞ; tÞ ¼ ðf1 ðxðtÞ; tÞ; f2 ðxðtÞ; tÞ; . . . ; fN ðxðtÞ; tÞÞT : For solving system (16) by NHPM we construct the following homotopy

X 0 ðtÞ ¼ x0 ðtÞ  pðx0 ðtÞ  FðXðtÞÞ  ef ðXðtÞ; tÞÞ:

ð17Þ n

By the new homotopy technique, we construct a homotopy X : X  ½0; 1 ! R , which satisfies

HðXðtÞ; pÞ ¼ X 0 ðtÞ  x0 ðtÞ þ pðx0 ðtÞ  FðXðtÞÞ  ef ðXðtÞ; tÞÞ ¼ 0;

ð18Þ

where p 2 [0, 1] is an embedding parameter, x0(t) is an initial approximation of solution of Eq. (16). Clearly, we have from Eq. (18)

HðUðtÞ; 0Þ ¼ X 0 ðtÞ  x0 ðtÞ ¼ 0; HðUðtÞ; 1Þ ¼ X 0 ðtÞ  FðXðtÞÞ  ef ðXðtÞ; tÞ ¼ 0:

ð19Þ ð20Þ

By applying Laplace transform on both sides of (18), we have

LfX 0 ðtÞ  x0 ðtÞ þ pðx0 ðtÞ  FðXðtÞÞ  ef ðXðtÞ; tÞÞg ¼ 0:

ð21Þ

Using the differential property of Laplace transform we have

sLfXðtÞg  Xð0Þ ¼ Lfx0 ðtÞ  pðx0 ðtÞ  FðXðtÞÞ  ef ðXðtÞ; tÞÞg;

ð22Þ

or

LfXðtÞg ¼

1 fXð0Þ þ Lfx0 ðtÞ  pðx0 ðtÞ  FðXðtÞÞ  ef ðXðtÞ; tÞÞgg: s

ð23Þ

By applying inverse  Laplace transform on both sides of (23), we have 

XðtÞ ¼ L1

1 fXð0Þ þ Lfx0 ðtÞ  pðx0 ðtÞ  FðXðtÞÞ  ef ðXðtÞ; tÞÞgg : s

ð24Þ

According to the HPM, we use the embedding parameter p as a small parameter, and assume that the solutions of Eq. (24) can be represented as a power series in p as

XðtÞ ¼

1 X

pn X n ðtÞ:

ð25Þ

n¼0

Now let us write Eq. (24) in the following form

( ( ( ! !!))) 1 1 X X 1 n n Xð0Þ þ L x0 ðtÞ  p x0 ðtÞ  F p X n ðtÞ ¼ L p X n ðtÞ  ef p X n ðtÞ; t : s n¼0 n¼0 n¼0

1 X

n

1

ð26Þ

Comparing coefficients of terms with identical powers of p, leads to

( ( ( ! !!))) 1 1 X X 1 n n Xð0Þ þ L x0 ðtÞ  p x0 ðtÞ  F p X n ðtÞ ¼ L p X n ðtÞ  ef p X n ðtÞ; t ; s n¼0 n¼0 n¼0   1 ðXð0Þ þ Lfx0 ðtÞgÞ ; p0 : X 0 ðtÞ ¼ L1 s  1 p1 : X 1 ðtÞ ¼ L1  ðLfx0 ðtÞ  FðX 0 ðtÞÞ  ef ðX 0 ðtÞ; tÞgÞ ;  s  1 1 2 ðLfFðX 0 ðtÞ; X 1 ðtÞÞ þ ef ðX 0 ðtÞ; X 1 ðtÞ; tÞgÞ ; p : X 2 ðtÞ ¼ L s .. .   1 pj : X j ðtÞ ¼ L1 ðLfFðX 0 ðtÞ; X 1 ðtÞ; . . . ; X j1 ðtÞÞ þ ef ðX 0 ðtÞ; X 1 ðtÞ; . . . ; X j1 ðtÞ; tÞgÞ ; s .. . 1 X

n

1

ð27Þ

H. Aminikhah / Applied Mathematical Modelling 36 (2012) 3638–3644

3641

Suppose that the initial approximation has the form x0(t) = X(0) = x(0), therefore the exact solution may be obtained as following

xðtÞ ¼ lim XðtÞ ¼ X 0 ðtÞ þ X 1 ðtÞ þ X 2 ðtÞ þ   

ð28Þ

p!1

4. Examples Example 1. For first numerical experiment, we have chosen a particular case of the singularly perturbed test problem suggested by Kaps [5], namely, the initial value problem

x01 ðtÞ ¼ ðe1 þ 2Þx1 ðtÞ þ e1 ðx2 ðtÞÞ2 ; x02 ðtÞ ¼ x1 ðtÞ  x2 ðtÞð1 þ x2 ðtÞÞ; x1 ð0Þ ¼ 1;

ð29Þ

x2 ð0Þ ¼ 1:

The exact solution of this problem x1(t) = e2t, x2(t) = et does not depend on ~e. However, the problem becomes very stiff, as ~ The eigenvalues of the system are k1 = 1 and k2 = (e1 + 2) which enables its degree of stiffness to be regulated. e ! 0. To solve Eq. (29) by the LTNHPM, we construct the following homotopy

X 01 ðtÞ  x1;0 ðtÞ þ pðx1;0 ðtÞ þ ðe1 þ 2ÞX 1 ðtÞ  e1 ðX 2 ðtÞÞ2 Þ ¼ 0;

ð30Þ

X 02 ðtÞ  x2;0 ðtÞ þ pðx2;0 ðtÞ  X 1 ðtÞ þ X 2 ðtÞð1 þ X 2 ðtÞÞ ¼ 0: Applying Laplace transform on both sides of (30), we have

n o L X 01 ðtÞ  x1;0 ðtÞ þ pðx1;0 ðtÞ þ ðe1 þ 2ÞX 1 ðtÞ  e1 ðX 2 ðtÞÞ2 Þ ¼ 0;   L X 02 ðtÞ  x2;0 ðtÞ þ pðx2;0 ðtÞ  X 1 ðtÞ þ X 2 ðtÞð1 þ X 2 ðtÞÞ ¼ 0:

ð31Þ

Using the differential property of Laplace transform we have

sLfX 1 ðtÞg  X 1 ð0Þ ¼ Lfx1;0 ðtÞ  pðx1;0 ðtÞ þ ðe1 þ 2ÞX 1 ðtÞ  e1 ðX 2 ðtÞÞ2 Þg;

ð32Þ

sLfX 2 ðtÞg  X 2 ð0Þ ¼ Lfx2;0 ðtÞ  pðx2;0 ðtÞ  X 1 ðtÞ þ X 2 ðtÞð1 þ X 2 ðtÞÞg; or

1 fX 1 ð0Þ þ Lfx1;0 ðtÞ  pðx1;0 ðtÞ þ ðe1 þ 2ÞX 1 ðtÞ  e1 ðX 2 ðtÞÞ2 Þgg; s 1 LfX 2 ðtÞg ¼ fX 2 ð0Þ þ Lfx2;0 ðtÞ  pðx2;0 ðtÞ  X 1 ðtÞ þ X 2 ðtÞð1 þ X 2 ðtÞÞgg: s

LfX 1 ðtÞg ¼

ð33Þ

By applying inverse Laplace transform on both sides of (33), we have

 1 fX 1 ð0Þ þ Lfx1;0 ðtÞ  pðx1;0 ðtÞ þ ðe1 þ 2ÞX 1 ðtÞ  e1 ðX 2 ðtÞÞ2 Þgg ; s   1 1 fX 2 ð0Þ þ Lfx2;0 ðtÞ  pðx2;0 ðtÞ  X 1 ðtÞ þ X 2 ðtÞð1 þ X 2 ðtÞÞgg : X 2 ðtÞ ¼ L s

X 1 ðtÞ ¼ L1



ð34Þ

Suppose the solution of Eq. (34) to have the following form

X 1 ðtÞ ¼ X 1;0 ðtÞ þ pX 1;1 ðtÞ þ p2 X 1;2 ðtÞ þ    ;

ð35Þ

X 2 ðtÞ ¼ X 2;0 ðtÞ þ pX 2;1 ðtÞ þ p2 X 2;2 ðtÞ þ    ;

where Xi,j(t), i = 1, 2, j = 1, 2, . . . are unknown functions which should be determined. Substituting Eq. (35) into Eq. (19), collecting the same powers of p and equating each coefficient of p to zero, results in

( p0 :

X 1;0 ðtÞ ¼ L1

1

ðX 1 ð0Þ s  L1 1s ðX 2 ð0Þ

 þ Lfx1;0 ðtÞgÞ ;  þ Lfx2;0 ðtÞgÞ ;

X 2;0 ðtÞ ¼ 8 n o < X 1;1 ðtÞ ¼ L1  1 ðx1;0 ðtÞ þ ðe1 þ 2ÞX 1;0 ðtÞ  e1 ðX 2;0 ðtÞÞ2 Þ ; s 1 p : : X ðtÞ ¼ L1  1 ðx ðtÞ  X ðtÞ þ X ðtÞð1 þ X ðtÞÞ; 2;1 2;0 1;0 2;0 2;0 1 s 8    j1 P > > X 2;k ðtÞX 2;jk1 ðtÞ ; > X 1;j ðtÞ ¼ L1  1s ðe1 þ 2ÞX 1;0 ðtÞ  e1 < k¼0 pj :    j1 > P > 1 1 > : X 2;j ðtÞ ¼ L  X 1;j1 ðtÞ þ X 2;j1 ðtÞ þ X 2;k ðtÞX 2;jk1 ðtÞ ; s

k¼0

ð36Þ

j ¼ 2; 3; . . .

3642

H. Aminikhah / Applied Mathematical Modelling 36 (2012) 3638–3644

Assuming x1,0(t) = X1(0) = x1(0) = 1 and x2,0(t) = X2(0) = x2(0) = 1. Solving the above equation for X1,j(t), X2,j(t), j = 0, 1, . . . leads to the result

 (

X 1;0 ðtÞ ¼ 1 þ t; X 2;0 ðtÞ ¼ 1 þ t; X 1;1 ðtÞ ¼ 3t þ

1

2e 2

3  1 t2 þ 3t e ; 3

X 2;1 ðtÞ ¼ 2t  t  t3 ;    X 1;2 ðtÞ ¼ 3  21e t 2 þ 23  2e  61e2 t3  65e þ 121e2 t4  152 e t 5 ;   2 5 t ; X 2;2 ðtÞ ¼ 32 t2 þ 61e þ 2 t 3 þ 121 e þ 34 t4 þ 15 ( 1  3 5 1 2 43 17 7 X 1;3 ðtÞ ¼ 6e2 þ 3e  2 t þ 24e3 þ 3e2 þ 12e  13 t 4 þ    þ 315 et ;  1 11 3  1 5 37 4 17 7 X 2;3 ðtÞ ¼  6e þ 6 t  24e2 þ 8e þ 12 t þ     315 t ; .. . (

Therefore we gain the solution of Eq. (29) as 1 X 4 2 4 4 ð1Þn ð2tÞn t þ  ¼ x1 ðtÞ ¼ X 1;0 ðtÞ þ X 1;1 ðtÞ þ X 1;3 ðtÞ þ    ¼ 1  2t þ 2t 2  t 3 þ t 4  ¼ e2t ; 3 3 15 n! n¼0 1 X 1 1 1 4 1 5 ð1Þn t n t  t þ  ¼ x2 ðtÞ ¼ X 2;0 ðtÞ þ X 2;1 ðtÞ þ X 2;3 ðtÞ þ    ¼ 1  t þ t2  t 3 þ ¼ et ; 2 6 24 120 n! n¼0

which is exact solution. Example 2. Consider the following stiff problem taken from [4]

x01 ðtÞ ¼ 2x1 ðtÞ þ x2 ðtÞ þ 2 sin t; x02 ðtÞ ¼ ðe1 þ 2Þx1 ðtÞ þ ðe1 þ 1Þðx2 ðtÞ  cos t þ sin tÞ; x1 ð0Þ ¼ 2;

ð37Þ

x2 ð0Þ ¼ 3:

The exact solution, independent of e, is x1(t) = 2et + sin t, x2(t) = 2et + cos t. The eigenvalues of the system are k1 = 1 and k2 = e1 which enables its degree of stiffness to be regulated. For solving system (37), by LTNHPM, using the Taylor series of sin t and cos t, we construct the following homotopy

X 01 ðtÞ ¼ x1;0 ðtÞ  p x1;0 ðtÞ þ 2X 1 ðtÞ  X 2 ðtÞ  2

1 X n¼0

X 02 ðtÞ

1

¼ x2;0 ðtÞ  p x2;0 ðtÞ þ ðe

1

þ 2ÞX 1 ðtÞ  ðe

pn

! ð1Þn t 2nþ1 ; ð2n þ 1Þ!

 2nþ1 !! 1 X t t2n n n þ 1Þ X 2 ðtÞ þ ð1Þ p  : ð2n þ 1Þ! ð2nÞ! n¼0

Using the differential property of Laplace transform we have

(

1 X sLfX 1 ðtÞg  X 1 ð0Þ ¼ L x1;0 ðtÞ  p x1;0 ðtÞ þ 2X 1 ðtÞ  X 2 ðtÞ  2 ð1Þn pn n¼0

( 1

sLfX 2 ðtÞg  X 2 ð0Þ ¼ L x2;0 ðtÞ  p x2;0 ðtÞ þ ðe

1

þ 2ÞX 1 ðtÞ  ðe

t 2nþ1 ð2n þ 1Þ!

ð38Þ

!) ;

 2nþ1 !!) 1 X t t 2n n n þ 1Þ X 2 ðtÞ þ ð1Þ p  ; ð2n þ 1Þ! ð2nÞ! n¼0 ð39Þ

or

(

(

!))

1 X 1 t 2nþ1 ; ð1Þn pn X 1 ð0Þ þ L x1;0 ðtÞ  p x1;0 ðtÞ þ 2X 1 ðtÞ  X 2 ðtÞ  2 s ð2n þ 1Þ! n¼0 ( (  2nþ1 !!)) 1 X 1 t t 2n n n 1 1 LfX 2 ðtÞg ¼ X 2 ð0Þ þ L x2;0 ðtÞ  p x2;0 ðtÞ þ ðe þ 2ÞX 1 ðtÞ  ðe þ 1Þ X 2 ðtÞ þ ð1Þ p  : s ð2n þ 1Þ! ð2nÞ! n¼0

LfX 1 ðtÞg ¼

ð40Þ By applying inverse Laplace transform on both sides of (40), we have

( ( ( !))) 1 2nþ1 X 1 n n t ; X 1 ð0ÞþL x1;0 ðtÞp x1;0 ðtÞþ2X 1 ðtÞX 2 ðtÞ2 ð1Þ p X 1 ðtÞ¼L s ð2nþ1Þ! n¼0 ( ( (  2nþ1 !!))) 1 X t t 2n n n 1 1 1 1 X 2 ðtÞ¼L X 2 ð0ÞþL x2;0 ðtÞp x2;0 ðtÞþðe þ2ÞX 1 ðtÞðe þ1Þ X 2 ðtÞþ ð1Þ p  : s ð2nþ1Þ! ð2nÞ! n¼0 1

ð41Þ

3643

H. Aminikhah / Applied Mathematical Modelling 36 (2012) 3638–3644

Suppose the solutions of system (41) to be in the following form

X i ðtÞ ¼ X i;0 ðtÞ þ pX i;1 ðtÞ þ p2 X i;2 ðtÞ þ    ;

i ¼ 1; 2;

ð42Þ

where Xi,j(t), i = 1, 2 and j = 0, 1, 2, . . . are functions which should be determined. Substituting Eq. (42) into Eq. (41), collecting the same powers of p and equating each coefficient of p to zero, results in

( p0 : ( p1 :

X 1;0 ðtÞ ¼ L1 X 2;0 ðtÞ ¼ X 1;1 ðtÞ ¼

1



ðX 1 ð0Þ þ Lfx1;0 ðtÞgÞ ; s   L1 1s ðX 2 ð0Þ þ Lfx2;0 ðtÞgÞ ;   L1  1s ðx1;0 ðtÞ þ 2X 1;0 ðtÞ  X 2;0 ðtÞ  2tÞ ;  L1  1s ðx2;0 ðtÞ þ ð 1 þ 2ÞX 1;0 ðtÞ  ð 1 þ

 e e 1ÞðX 2;0 ðtÞ  1 þ tÞÞ ; X 2;1 ðtÞ ¼ 8 n

o t2j1 > < X 1;j ðtÞ ¼ L1  1s 2X 1;j1 ðtÞ  X 2;j1 ðtÞ þ 2ð1Þj ð2j1Þ! ; j n



o p : 1 > 1 1 1 t 2j2 t 2j1 : X 2;j ðtÞ ¼ L  s ðe þ 2ÞX 1;0 ðtÞ  ðe þ 1Þ X 2;0 ðtÞ þ ð1Þj ð2j2Þ!  ð2j1Þ! ;

ð43Þ j ¼ 2; 3; . . .

Assuming x1,0(t) = X1(0) = x1(0) = 2 and x2,0(t) = X2(0) = x2(0) = 3. Solving the above equation for X1,j(t), X2,j(t), j = 0, 1, . . . leads to the result

 (

X 1;0 ðtÞ ¼ 2ð1 þ tÞ; X 2;0 ðtÞ ¼ 3ð1 þ tÞ; X 1;1 ðtÞ ¼ 3t þ 12 t 2 ;

X 2;1 ðtÞ ¼ 5t þ 1e t 2 ;  1 4 X 1;2 ðtÞ ¼ 12 t 2 þ 31e  13 t 3  12 t ; 1 1 2  1 1 1 þ 241 e t 4 ; X 2;2 ðtÞ ¼ 2  e t þ 3e þ 3e2  16 t 3  24 (   1 5 1 1 6  120 X 1;3 ðtÞ ¼  16 þ 31e t3 þ 18 þ 121e2  121 e t4 þ 40 e t þ 360 t ; 1   5 1 6 1 1 1 1  120 X 2;3 ðtÞ ¼  3e þ 31e2 þ 16 t3 þ 18 þ 121e2 þ 241 e þ 121e3 t 4 þ 60 e  120e2 t þ 720 þ 720e t ;

(

.. . Therefore we gain the solution of Eq. (37) as

1 1 4 1 5 t  t þ  x1 ðtÞ ¼ X 1;0 ðtÞ þ X 1;1 ðtÞ þ X 1;3 ðtÞ þ    ¼ 2  t þ t2  t3 þ 2 12 120     n 1 1 X ð1Þ tn X 1 1 1 5 ð1Þn t2nþ1 ¼ 2 1  t þ t2     þ t  t3 þ t   ¼ 2 þ ¼ 2et þ sin t; 2 6 120 n! ð2n þ 1Þ! n¼0 n¼0 1 1 1 1 5 t þ  x2 ðtÞ ¼ X 2;0 ðtÞ þ X 2;1 ðtÞ þ X 2;3 ðtÞ þ    ¼ 3  2t þ t2  t3 þ t 4  2 3 8 60     n 1 1 X ð1Þ tn X 1 1 1 4 ð1Þn t2n ¼ 2 1  t þ t2     þ 1  t2 þ t   ¼ 2 þ ¼ 2et þ cos t; 2 2 24 n! ð2nÞ! n¼0 n¼0 which is exact solution.

5. Conclusion In this paper, we have introduced a combination of Laplace transform and new homotopy perturbation methods for stiff problem which we called LTNHPM. We have discussed the methodology for the construction of these schemes and studied their performance on several test problems. The solution is very rapidly convergent by utilizing the new homotopy perturbation method by modification of Laplace operator. References [1] [2] [3] [4]

J.E. Flaherty, R.E. OMalley, The numerical solution of boundary value problems for stiff differential equations, Math. Comput. 31 (1997) 66–93. T.D. Bui, T.R. Bui, Numerical methods for extremely stiff systems of ordinary differential equations, Appl. Math. Modell. 3 (1979) 355–358. J.C. Butcher, Numerical Methods for Ordinary Differential Equations, Wiley, New York, 2003. L.Gr. Ixaru, G. Vanden Berghe, H. De Meyer, Frequency evaluation in exponential fitting multistep algorithms for ODEs, J. Comput. Appl. Math. 140 (2000) 423–434. [5] P. Kaps, Rosenbrock-type methods, in: G. Dahlquist, R. Jeltsch, (Eds.), Numerical methods for stiff initial value problems, Berich Nr. 9, Inst. Fur Geometric und Practische Mathematik der RWTH Aachen, 1981. [6] J.H. He, Homotopy perturbation technique, Comput. Methods Appl. Mech. Eng. 178 (1999) 257–262. [7] J.H. He, The homotopy perturbation method for nonlinear oscillators with discontinuities, Appl. Math. Comput. 151 (2004) 287–292.

3644 [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18]

H. Aminikhah / Applied Mathematical Modelling 36 (2012) 3638–3644

J.H. He, Application of homotopy perturbation method to nonlinear wave equations, Chaos Soliton. Fract. 26 (2005) 695–700. J.H. He, Homotopy perturbation method for solving boundary value problems, Phys. Lett. A 350 (2006) 87–88. J. Biazar, H. Ghazvini, Exact solutions for nonlinear Schrödinger equations by He’s homotopy perturbation method, Phys. Lett. A 366 (2007) 79–84. D.D. Ganji, The application of He’s homotopy perturbation method to nonlinear equations arising in heat transfer, Phys. Lett. A 355 (2006) 337–341. Z. Odibat, S. Momani, Modified homotopy perturbation method: application to quadratic Riccati differential equation of fractional order, Chaos. Soliton. Fract. 36 (2008) 167–174. L. Cveticanin, Homotopy perturbation method for pure nonlinear differential equation, Chaos. Soliton. Fract. 30 (2006) 1221–1230. H. Aminikhah, M. Hemmatnezhad, An efficient method for quadratic Riccati differential equation, Commun. Nonlinear Sci. Numer. Simul. 15 (2010) 835–839. H. Aminikhah, J. Biazar, A new HPM for ordinary differential equations, Numer. Methods Partial Differ. Equat. 26 (2009) 480–489. A. Yıldırım, H. Koçak, Homotopy perturbation method for solving the space-time fractional advection–dispersion equation, Adv. Water Res. 32 (2009) 1711–1716. M.E. Berberler, A. Yıldırım, He’s homotopy perturbation method for solving the shock wave equation, Appl. Anal. 88 (2009) 997–1004. H. Aminikhah, An analytical approximation for solving nonlinear Blasius equation by NHPM, Numer. Methods Partial Differ. Equat. 26 (2010) 1291– 1299.