The multistage homotopy-perturbation method: A powerful scheme for handling the Lorenz system

Available online at www.sciencedirect.com

Chaos, Solitons and Fractals 40 (2009) 1929–1937 www.elsevier.com/locate/chaos

The multistage homotopy-perturbation method: A powerful scheme for handling the Lorenz system M.S.H. Chowdhury a, I. Hashim a, S. Momani a

b,*

School of Mathematical Sciences, National University of Malaysia, 43600 Bangi Selangor, Malaysia b Department of Mathematics and Physics, College of Arts and Sciences, Qatar University, Qatar Accepted 17 September 2007

Communicated by J.-H. He

Abstract In this paper, a new reliable algorithm based on an adaptation of the standard homotopy-perturbation method (HPM) is presented. The HPM is treated as an algorithm in a sequence of intervals (i.e. time step) for finding accurate approximate solutions to the famous Lorenz system. Numerical comparisons between the multistage homotopy-perturbation method (MHPM) and the classical fourth-order Runge–Kutta (RK4) method reveal that the new technique is a promising tool for the nonlinear systems of ODEs. Ó 2007 Elsevier Ltd. All rights reserved.

1. Introduction In this paper, we consider the famous Lorenz system [16], dx ¼ rðy  xÞ; dt dy ¼ Rx  y  xz; dt dz ¼ xy þ bz; dt

ð1Þ ð2Þ ð3Þ

where x, y and z are respectively proportional to the convective velocity, the temperature difference between descending and ascending flows, and the mean convective heat flow, and r, b and the so-called bifurcation parameter R are real constants. The motivation of this paper is to extend the application of the analytic homotopy-perturbation method (HPM) [6,11–15,23] to solve the Lorenz system (1)–(3). The homotopy perturbation method (HPM) was first proposed by Chinese mathematician He [10–15]. The essential idea of this method is to introduce a homotopy parameter, say p,

*

Corresponding author. E-mail address: [email protected] (S. Momani).

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

1930

M.S.H. Chowdhury et al. / Chaos, Solitons and Fractals 40 (2009) 1929–1937

which takes the values from 0 to 1. When p ¼ 0, the system of equations usually reduces to a sufficiently simplified form, which normally admits a rather simple solution. As p gradually increases to 1, the system goes through a sequence of ‘deformation’, the solution of each of which is ‘close’ to that at the previous stage of ‘deformation’. Eventually at p ¼ 1, the system takes the original form of the equation and the final stage of ‘deformation’ gives the desired solution. The HPM has been employed to solve a large variety of linear and nonlinear problems [11,12,15,3–5]. In particular, Gill et al. [7] studied the ion-acoustic solitons in a weakly relativistic electron–positron–ion plasma using HPM. In [19], Noor and Mohyud-Din presented an algorithm based on HPM for solving a boundary-value problem. Shakeri and Dehghan [22] applied HPM to an inverse problem of diffusion equation. Recently, the application of the method has been extended to linear and nonlinear differential equations of fractional order [17,18,21]. The HPM yields a very rapid convergence of the solution series in most cases, usually only a few iterations leading to very accurate solutions. Thus He’s HPM is a universal one which can solve various kinds of nonlinear equations. To the best of our knowledge, this work is the first attempt to apply HPM to systems of nonlinear ODEs capable of showing chaotic behaviours such as the Lorenz system. It was shown that the hybrid numeric–analytic procedure of the Adomian decomposition method (ADM) [2] is a powerful scheme for the chaotic Lorenz and Chen systems [9,20,1]. The success the ADM has had thus far motivates us to develop a new hybrid numeric–analytic procedure based on the HPM. In this paper, the HPM is treated as an algorithm in a sequence of intervals (i.e. time step) for finding accurate approximate solutions to the Lorenz system. We shall call this technique as the multistage HPM (for short MHPM). Comparison with the classical fourth-order Runge-Kutta (RK4) shall be made.

2. Solution approaches Consider a general systems of first-order ODEs du1 þ g1 ðt; u1 ; u2 ; . . . ; um Þ ¼ f1 ðtÞ; dt du2 þ g2 ðt; u1 ; u2 ; . . . ; um Þ ¼ f2 ðtÞ; dt .. . dum þ gm ðt; u1 ; u2 ; . . . ; um Þ ¼ fm ðtÞ dt

ð4Þ ð5Þ

ð6Þ

subject to the initial conditions u1 ðt0 Þ ¼ c1 ; u2 ðt0 Þ ¼ c2 ; . . . ; um ðt0 Þ ¼ cm :

ð7Þ

First write system (4)–(6) in the operator form Lðu1 Þ þ N 1 ðu1 ; u2 ; . . . ; um Þ  f1 ¼ 0; Lðu2 Þ þ N 2 ðu1 ; u2 ; . . . ; um Þ  f2 ¼ 0; .. . Lðum Þ þ N m ðu1 ; u2 ; . . . ; um Þ  fm ¼ 0

ð8Þ ð9Þ

ð10Þ

subject to the initial conditions (7), where L ¼ d=dt is a linear operator and N 1 ; N 2 ; . . . ; N m are the nonlinear operators. We shall next present the solution approaches of (8)–(10) based on the standard HPM and MHPM separately. 2.1. Solution by HPM According to HPM, we construct a homotopy for (8)–(10) which satisfies the following relations: Lðu1 Þ  Lðv1 Þ þ pLðv1 Þ þ p½N 1 ðu1 ; u2 ; . . . ; um Þ  f1  ¼ 0; Lðu2 Þ  Lðv2 Þ þ pLðv2 Þ þ p½N 2 ðu1 ; u2 ; . . . ; um Þ  f2  ¼ 0; .. .

ð11Þ ð12Þ

Lðum Þ  Lðvm Þ þ pLðvm Þ þ p½N m ðu1 ; u2 ; . . . ; um Þ  fm  ¼ 0;

ð13Þ

M.S.H. Chowdhury et al. / Chaos, Solitons and Fractals 40 (2009) 1929–1937

1931

where p 2 ½0; 1 is an embedding parameter and v1 ; v2 ; . . . ; vm are initial approximations which satisfying the given conditions. It is obvious that when the perturbation parameter p ¼ 0, Eqs. (11)–(13) become a linear system of equations and when p ¼ 1 we get the original nonlinear system of equations. Let us take the initial approximations as follows: u1;0 ðtÞ ¼ v1 ðtÞ ¼ u1 ðt0 Þ ¼ c1 ; u2;0 ðtÞ ¼ v2 ðtÞ ¼ u2 ðt0 Þ ¼ c2 ; .. . um;0 ðtÞ ¼ vm ðtÞ ¼ um ðt0 Þ ¼ cm

ð14Þ ð15Þ

ð16Þ

and u1 ðtÞ ¼ u1;0 ðtÞ þ pu1;1 ðtÞ þ p2 u1;2 ðtÞ þ p3 u1;3 ðtÞ þ    ; 2

3

ð17Þ

u2 ðtÞ ¼ u2;0 ðtÞ þ pu2;1 ðtÞ þ p u2;2 ðtÞ þ p u2;3 ðtÞ þ    ; .. .

ð18Þ

um ðtÞ ¼ um;0 ðtÞ þ pum;1 ðtÞ þ p2 um;2 ðtÞ þ p3 um;3 ðtÞ þ    ;

ð19Þ

where ui;j , ði ¼ 1; 2; . . . ; m; j ¼ 1; 2; . . .Þ are functions yet to be determined. Substituting (14)–(19) into (11)–(13) and arranging the coefficients of the same powers of p, we get Lðu1;1 Þ þ Lðv1 Þ þ N 1 ðu1;0 ; u2;0 ; . . . ; um;0 Þ  f1 ¼ 0;

u1;1 ðt0 Þ ¼ 0;

ð20Þ

Lðu2;1 Þ þ Lðv2 Þ þ N 2 ðu1;0 ; u2;0 ; . . . ; um;0 Þ  f2 ¼ 0; u2;1 ðt0 Þ ¼ 0; .. . Lðum;1 Þ þ Lðvm Þ þ N m ðu1;0 ; u2;0 ; . . . ; um;0 Þ  fm ¼ 0; um;1 ðt0 Þ ¼ 0;

ð21Þ

Lðu1;2 Þ þ N 1 ðu1;1 ; u2;1 ; . . . ; um;1 Þ ¼ 0; Lðu2;2 Þ þ N 2 ðu1;1 ; u2;1 ; . . . ; um;1 Þ ¼ 0; .. .

u1;2 ðt0 Þ ¼ 0; u2;2 ðt0 Þ ¼ 0;

ð23Þ ð24Þ

Lðum;2 Þ þ N m ðu1;1 ; u2;1 ; . . . ; um;1 Þ ¼ 0;

um;2 ðt0 Þ ¼ 0;

ð25Þ

ð22Þ

etc. We solve the above systems of equations for the unknowns ui;j ði ¼ 1; 2; . . . ; m; j ¼ 1; 2; . . .Þ by applying the inverse operator Z t 1 ðÞdt: ð26Þ L ðÞ ¼ 0

Therefore, according to HPM the n-term approximations for the solutions of (8)–(10) can be expressed as /1;n ðtÞ ¼ u1 ðtÞ ¼ lim u1 ðtÞ ¼ p!1

/2;n ðtÞ ¼ u2 ðtÞ ¼ lim u2 ðtÞ ¼ p!1

n1 X

u1;k ðtÞ;

ð27Þ

u2;k ðtÞ;

ð28Þ

k¼0 n1 X k¼0

.. . /m;n ðtÞ ¼ um ðtÞ ¼ lim um ðtÞ ¼ p!1

n1 X

um;k ðtÞ:

ð29Þ

k¼0

2.2. Solution by MHPM The approximate solutions (27)–(29) are generally, as shall be shown in the numerical experiments of this paper, not valid for large t. A simple way of ensuring validity of the approximations for large t is to treat (20)–(25) as an algorithm for approximating the solutions of (4)–(7) in a sequence of intervals choosing the initial approximations as

1932

M.S.H. Chowdhury et al. / Chaos, Solitons and Fractals 40 (2009) 1929–1937

u1;0 ðtÞ ¼ v1 ðtÞ ¼ u1 ðt Þ ¼ c1 ; u2;0 ðtÞ ¼ v2 ðtÞ ¼ u2 ðt Þ ¼ c2 ; .. . um;0 ðtÞ ¼ vm ðtÞ ¼ um ðt Þ ¼ cm :

ð30Þ ð31Þ

ð32Þ

Now we solve (20)–(25) for the unknowns ui;j ði ¼ 1; 2; . . . ; m; j ¼ 1; 2; . . .Þ by applying the inverse linear operator Z t L1 ðÞ ¼ ðÞdt: ð33Þ t

In order to carry out the iterations in every subinterval of equal length Dt, ½0; t1 Þ, ½t1 ; t2 Þ, ½t2 ; t3 Þ . . . ½tj1 ; tÞ, we would need to know the values of the following: u1;0 ðtÞ ¼ u1 ðt Þ; u2;0 ðtÞ ¼ u2 ðt Þ; . . . ; um;0 ðtÞ ¼ um ðt Þ:

ð34Þ 

But, in general, we do not have these information at our clearance except at the initial point t ¼ t0 . A simple way for obtaining the necessary values could be by means of the previous n-term approximations /1;n ; /2;n ; . . . ; /m;n of the preceding subinterval, i.e. u1;0 ’ /1;n ðt Þ; u2;0 ’ /2;n ðt Þ; . . . ; um;0 ’ /m;n ðt Þ:

ð35Þ

3. Application In this section, we shall study the Lorenz system (1)–(3) subject to the initial conditions xðt Þ ¼ c1 ;

yðt Þ ¼ c2 ;

zðt 0Þ ¼ c3 :

ð36Þ

According to the HPM, we can construct a homotopy which satisfies the following relations: v01  x00 þ pðx00  rv2 þ rv1 Þ ¼ 0;

ð37Þ

v02  y 00 þ pðy 00  Rv2 þ v2 þ v1 v3 Þ ¼ 0; v03  z00 þ pðz00  bv3  v1 v2 Þ ¼ 0:

ð38Þ ð39Þ

The initial approximations we take are v1;0 ðtÞ ¼ x0 ðtÞ ¼ xðt Þ ¼ c1 ; v2;0 ðtÞ ¼ y 0 ðtÞ ¼ yðt Þ ¼ c2 ; v3;0 ðtÞ ¼ z0 ðtÞ ¼ zðt Þ ¼ c3 ;

ð40Þ ð41Þ ð42Þ

and v1 ðtÞ ¼ v1;0 ðtÞ þ pv1;1 ðtÞ þ p2 v1;2 ðtÞ þ p3 v1;3 ðtÞ þ    ; 2

3

ð43Þ

v2 ðtÞ ¼ v2;0 ðtÞ þ pv2;1 ðtÞ þ p v2;2 ðtÞ þ p v2;3 ðtÞ þ    ;

ð44Þ

v3 ðtÞ ¼ v3;0 ðtÞ þ pv3;1 ðtÞ þ p2 v3;2 ðtÞ þ p3 v3;3 ðtÞ þ    ;

ð45Þ

where vi;j , i; j ¼ 1; 2; 3; . . . are functions yet to be determined. Substituting (40) and (43) into (37) and collecting terms the same powers of p, we have v01;1  rv2;0 þ rv1;0 ¼ 0;

ð46Þ

v02;1 v03;1 v01;2 v02;2 v03;2 v01;3 v02;3 v03;3

ð47Þ

 Rv1;0 þ v2;0 þ v1;0 v3;0 ¼ 0;  bv3;0  v1;0 v2;0 ¼ 0;

ð48Þ

 rv2;1 þ rv1;1 ¼ 0;

ð49Þ

 Rv1;1 þ v2;1 þ v1;0 v3;1 þ v1;1 v3;0 ¼ 0;

ð50Þ

 bv3;1  v1;0 v2;1  v1;1 v2;0 ¼ 0;

ð51Þ

 rv2;2 þ rv1;2 ¼ 0;

ð52Þ

 Rv1;2 þ v2;2 þ v1;2 v3;0 þ v1;1 v3;1 þ v1;0 v3;2 ¼ 0;

ð53Þ

 bv3;2  v1;1 v2;1  v1;2 v2;0  v1;0 v2;2 ¼ 0:

ð54Þ

M.S.H. Chowdhury et al. / Chaos, Solitons and Fractals 40 (2009) 1929–1937

1933

In order to obtain the nine unknowns vi;j ðtÞ, i; j ¼ 1; 2; 3; we solve the above system which includes nine equations with nine unknowns, considering the initial conditions vi;j ð0Þ ¼ 0, i; j ¼ 1; 2; 3, we obtain v1;1 ðtÞ ¼ ðrc1 þ rc2 Þðt  t Þ; v2;1 ðtÞ ¼ ðRc1  c2  c1 c3 Þðt  t Þ; v3;1 ðtÞ ¼ ðbc3 þ c1 c2 Þðt  t Þ; 1 v1;2 ðtÞ ¼ ½rðrc2  rc1 Þ þ rðRc1  c2  c3 c1 Þðt  t Þ2 ; 2 1 v2;2 ðtÞ ¼ ½Rðrc2  rc1 Þ  ðRc1  c2  c1 c3 Þ  c1 ðbc3 þ c1 c2 Þ  ðrc2  rc1 Þc3 ðt  t Þ2 ; 2 1 v3;2 ðtÞ ¼ ½bðbc3 þ c1 c2 Þ þ c1 ðRc1  c2  c1 c3 Þ þ c2 ðrc2  rc1 Þðt  t Þ2 ; 2 1 v1;3 ðtÞ ¼ ½rðrRc1  rc2  rc1 c3  r2 c2 þ r2 c1 Þ þ rðrc3 c2 þ rc3 c1 þ rRc2  rRc1  Rc1 6 þ c2 þ c3 c1  bc3 c1  c21 c2 Þðt  t Þ3 ; 1 v2;3 ðtÞ ¼ ½RðRrc1  rc2  rc3 c1  r2 c2 þ r2 c1 Þ  ðrc3 c1  rc3 c2 þ rRc2  rRc1  Rc1 þ c2 þ c1 c3 Þ 6  bc1 c3  c21 c2  c1 ðb2 c3 þ bc1 c2 þ Rc21  c1 c2  c21 c3 þ rc22  rc1 c2 Þ  2ðbrc3 c2 þ rc1 c22  rbc3 c1

ð55Þ ð56Þ ð57Þ ð58Þ ð59Þ ð60Þ

ð61Þ

 rc21 c2 Þ  c3 ðrRc1  rc2  rc3 c1  r2 c2 þ r2 c1 Þðt  t Þ3 ; ð62Þ 1 v3;3 ðtÞ ¼ ½bðb2 c3 þ bc1 c2 þ Rc21  c1 c2  c21 c3 þ rc22  rc2 c1 Þ þ c1 ðrc1 c3  rc3 c2 þ rRc2  rRc1  Rc1 þ c2 þ c1 c3 6  bc1 c3  c21 c2 Þ þ c2 ðrRc1  rc2  rc1 c3  r2 c2 þ r2 c1 Þ þ 2ðRrc1 c2  rc22  rc1 c3 c2  Rrc21 þ rc1 c2 þ rc21 c3 Þðt  t Þ3 :

ð63Þ

Hence, the explicit solution to the Lorenz system (1)–(3) is 1 X v1;m ðtÞ; x¼ y¼ z¼

m¼0 1 X m¼0 1 X

ð64Þ

v2;m ðtÞ;

ð65Þ

v3;m ðtÞ:

ð66Þ

m¼0

To carry out the iterations in every subinterval of equal length Dt, ½0; t1 Þ, ½t1 ; t2 Þ, ½t2 ; t3 Þ . . . ½tn1 ; tÞ, we would need to know the values of the following initial conditions: c1 ¼ xðt Þ;

c2 ¼ yðt Þ;

c3 ¼ zðt Þ:

ð67Þ 

In general, we do not have these information at our clearance except at the initial point t ¼ t0 ¼ 0 but we can obtain these values following the MHPM as given in Section approximations of x, y and z are P P (2.2). We note that the 10-term P denoted as xðtÞ ’ /10 ðtÞ ¼ 9i¼0 v1;i , yðtÞ ’ w10 ðtÞ ¼ 9i¼0 v2;i and zðtÞ ’ n10 ðtÞ ¼ 9i¼0 v3;i . For practical computations, a finite number of terms in the series solution are used in a time step procedure just outlined. 4. Results and discussion For a direct comparison with Guellal et al. [8] and Hashim et al. [9], we set r ¼ 10, b ¼ 8=3 and take the initial conditions xð0Þ ¼ 15:8, yð0Þ ¼ 17:48 and zð0Þ ¼ 35:64. The time range studied in this work is ½0; 20 as considered in [8,9]. In this analysis, we attempt to demonstrate the accuracy of the MHPM for the solutions of both non-chaotic and chaotic systems. In this work we fix the number of terms used to be ten and we select the time step size Dt ¼ 0:01 for non-chaotic and time step size Dt ¼ 0:001 for chaotic in MHPM solutions and choose time step Dt ¼ 0:001 for the RK4 method [9]. 4.1. Non-chaotic solutions First we consider the non-chaotic case R ¼ 23:5 with the other parameters as given above. The 10-term HPM series solutions to the Lorenz system (1)–(3) for this case are given by

1934

M.S.H. Chowdhury et al. / Chaos, Solitons and Fractals 40 (2009) 1929–1937

x ¼ 15:8  16:80t þ 1130:460t2 þ 993:025333t3  35589:992422t4  153747:229879t5 þ 1540403:337229t6 þ 7267052:894670t7  51567839:729863t8  362718526:968872t9 ; 2

3

ð68Þ

4

y ¼ 17:48 þ 209:2920t þ 1428:36760t  13242:971635t  112463:607362t þ 770494:772458t

5

þ 6627340:363498t6  33987218:889220t7  378014514:001849t8 þ 1312791002:118400t9 ; 2

3

4

z ¼ 35:64 þ 181:1440t  1748:100133t  13727:718041t þ 110271:943941t þ 830001:233062t

ð69Þ 5

 5135201:506905t6  50598363:252419t7 þ 220874715:920757t8 þ 2763659177:836651t9 :

ð70Þ

In Table 1 we present the absolute errors between the 10-term HPM solutions and the 10-term MHPM solutions and the RK4 solutions on time step Dt ¼ 0:01 for R ¼ 23:5. The 10-term MHPM solutions on the time step Dt ¼ 0:01 agree with the RK4 solutions to at least seven decimal places, while the 10-term classical HPM solutions are only valid for t  1. For the non-chaotic case we see that the MHPM has the advantage over the RK4 on achieving a good accuracy with a larger time step.

Table 1 Differences between 10-term HPM and 10-term MHPM with RK4 solutions for R ¼ 23:5 D ¼ jHPM  RK40:001 j

t

2 4 6 8 10 12 14 16 18 20

D ¼ jMHPM0:01  RK40:001 j

Dx

Dy

Dz

Dx

Dy

Dz

1.979E+11 9.834E+13 3.740E+15 4.953E+16 3.678E+17 1.893E+18 7.569E+18 2.515E+19 7.251E+19 1.870E+20

5.715E+11 3.188E+14 1.259E+16 1.698E+17 1.275E+18 6.610E+18 2.656E+19 8.858E+19 2.562E+20 6.624E+20

1.465E+12 7.381E+14 2.821E+16 3.745E+17 2.785E+18 1.435E+19 5.742E+19 1.909E+20 5.506E+20 1.421E+21

1.808E09 3.190E09 2.676E09 1.916E10 5.986E09 1.391E08 1.595E08 2.814E09 1.469E08 2.171E08

1.854E09 4.960E09 6.124E09 5.001E09 1.478E09 8.112E09 2.482E08 2.692E08 8.032E10 2.308E08

2.481E09 1.430E09 3.585E09 1.003E08 1.646E08 1.929E08 7.172E09 2.453E08 4.246E08 2.961E08

40

15 10

35

5 30

y

z

0 -5

25

-10 20

-15 -20 -20

15 -15

-10

-5

0

5

10

-20

15

-15

-10

-5

0

5

10

15

x

x 40 35

z 40 35 30 25 20 15

z

30 25 20 15 -20

-15

-10

-5

0

y

5

10

15

-20 -15 -10 -5

x

0

5 10 15

Fig. 1. Phase portraits using 10-term MHPM on Dt ¼ 0:01 for R ¼ 23:5.

-20

-10-5 -15

0 5

y

15 10

M.S.H. Chowdhury et al. / Chaos, Solitons and Fractals 40 (2009) 1929–1937

1935

The x  y, x  z, y  z and x  y  z phase portraits obtained using the 10-term MHPM solutions on Dt ¼ 0:01 are shown in Fig. 1. 4.2. Chaotic solutions Finally, we consider the case R ¼ 28 with the other parameters as given above which exhibits chaotic solutions. The 10-term HPM series solution to the Lorenz system (1)–(3) for this case can be written as x ¼ 15:8  16:80t þ 774:960t2 þ 2170:525333t3  24633:707422t4  128128:338879t5 þ 676062:199037t6 þ 5675321:096379t7  14338709:761281t8  225821492:890099t9 ;

ð71Þ

Table 2 Differences between 10-term HPM and 10-term MHPM with RK4 solutions for R ¼ 28 D ¼ jHPM  RK40:001 j

y

2 4 6 8 10 12 14 16 18 20

D ¼ jMHPM0:001  RK40:001 j

Dx

Dy

Dz

Dx

Dy

Dz

1.185E+11 6.004E+13 2.298E+15 3.054E+16 2.272E+17 1.171E+18 4.686E+18 1.558E+19 4.495E+19 1.160E+20

9.993E+10 3.676E+13 1.229E+15 1.514E+16 1.073E+17 5.351E+17 2.089E+18 6.813E+18 1.937E+19 4.937E+19

7.318E+11 3.753E+14 1.443E+16 1.921E+17 1.431E+18 7.383E+18 2.956E+19 9.831E+19 2.838E+20 7.324E+20

4.413E09 2.041E08 3.098E08 1.454E07 1.173E06 4.896E07 3.992E05 6.830E05 5.619E04 7.503E03

5.700E09 1.926E09 5.767E08 2.386E07 9.988E06 5.279E06 1.576E05 7.706E05 8.503E05 1.243E02

8.405E09 2.860E08 7.697E09 5.281E08 7.549E06 3.714E06 1.190E04 9.285E05 9.133E04 6.807E03

25 20 15 10 5 0 -5 -10 -15 -20 -25 -20

45 40 35 30

z

t

25 20 15 10

-15

-1

-5

0

x

5

10

15

5 -20

20

-15

-10

-5

0

x

5

10

15

20

45 40

z

35

z

30

45 35 25 15 5

25 20 15 10 5 -25 -20 -15 -10

-5

0

y

5

10

15

20

25

-20 -15 -10 -5

0

x

5 10 15 20

Fig. 2. Phase portraits using 10-term MHPM on Dt ¼ 0:01 for R ¼ 28.

-25

-15

-5

5

y

15

25

1936

M.S.H. Chowdhury et al. / Chaos, Solitons and Fractals 40 (2009) 1929–1937

y ¼ 17:48 þ 138:1920t þ 1426:11760t2  7682:957635t3  88697:876862t4 þ 277508:980544t5 þ 4648786:966503t6  5795646:712647t7  217578053:362370t8  85529221:224559t9 ;

ð72Þ

z ¼ 35:64 þ 181:1440t  1186:410133t2  11745:608041t3 þ 49476:516475t4 þ 646858:660665t5  1440422:065397t6  32048523:232191t7 þ 15278390:313161t8 þ 1429750324:578864t9 :

ð73Þ

The differences between the 10-term HPM solutions and 10-term MHPM solutions Dt ¼ 0:001 and the RK4 solutions on Dt ¼ 0:001 are given in Table 2. We observe that the MHPM solutions are accurate enough as compared with the RK4 solutions, but the same conclusion as in the non-chaotic case applies to the standard HPM. In Fig. 2 we reproduce the well-known x  y, x  z, y  z and x  y  z phase portraits of the chaotic Lorenz system using the 10-term MHPM solutions on Dt ¼ 0:01.

5. Conclusions In this work, the MHPM was applied to the solutions of the well known Lorenz system. Comparisons between the decomposition solutions and the fourth-order Runge–Kutta (RK4) numerical solutions were made. For the non-chaotic case studied we found that the 10-term MHPM solutions on a larger time step achieved comparable accuracy compared with the RK4 solutions on a much smaller time step. Whereas for the chaotic case, both the 10-term MHPM and the RK4 solutions are of comparable accuracy on the same time step Dt ¼ 0:001. The MHPM was shown to be a simple, yet powerful analytic-numeric scheme for handling chaotic systems.

Acknowledgments The authors acknowledge the financial supports received from the Academy of Sciences Malaysia under the SAGA Grant No. P24c (STGL-011-2006), the Malaysian Technical Cooperation Program and the International Islamic University Chittagong, Bangladesh.

References [1] Abdulaziz O, Noor NFM, Hashim I, Noorani MSM. Further accuracy tests on Adomian decomposition method for chaotic systems. Chaos Solitons & Fractals 2008;36(5):1405–1411. [2] Adomian G. Nonlinear stochastics systems theory and application to physics. Dordrecht: Reidel; 1988. [3] Chowdhury MSH, Hashim I. Solutions of a class of singular second-order IVPs by homotopy-perturbation method. Phys Lett A 2007;365:439–47. [4] Chowdhury MSH, Hashim I, Abdulaziz O. Application of homotopy-perturbation method to nonlinear population dynamics models. Phys Lett A 2007;368:251–8. [5] Chowdhury MSH, Hashim I. Solutions of time-dependent Emden–Fowler type equations by homotopy-perturbation method. Phys Lett A, doi:10.1016/j.physleta.2007.04.020. [6] Chowdhury MSH, Hashim I. Application of homotopy-perturbation method to Klein–Gordon and sine-Gordon equations. Chaos Solitons & Fractals, in press. [7] Gill TS, Singh A, Kaur H, Saini NS, Bala P. Ion-acoustic solitons in weakly relativistic plasma containing electron–positron and ion. Phys Lett A 2007;361(4–5):364–7. [8] Guellal S, Grimalt P, Cherruault Y. Numerical study of Lorenz’s equation by the Adomian method. Comput Math Appl 1997;33:25–9. [9] Hashim I, Noorani MSM, Ahmad R, Bakar SA, Ismail ES, Zakaria AM. Accuracy of the Adomian decomposition method applied to the Lorenz system. Chaos Solitons & Fractals 2006;28:1149–58. [10] He JH. Approximate analytical solution of Blasius’ equation. Commun Nonlinear Sci Numer Simul 1998;3:260–3. [11] He JH. A coupling method of homotopy technique and perturbation technique for nonlinear problems. Int J Non-Linear Mech 2000;35:37–43. [12] He JH. Application of homotopy perturbation method to nonlinear wave equations. Chaos Solitons & Fractals 2005;26:695–700. [13] He JH. Non-perturbative methods for strongly nonlinear problems. Germany: Die Deutsche bibliothek; 2006.

M.S.H. Chowdhury et al. / Chaos, Solitons and Fractals 40 (2009) 1929–1937 [14] [15] [16] [17] [18] [19] [20] [21] [22] [23]

1937

He JH. Some asymptotic methods for strongly nonlinear equations. Int J Mod Phys B 2006;20(10):1141–99. He JH. Homotopy perturbation method for solving boundary value problems. Phys Lett A 2006;350:87–8. Lorenz EN. Deterministic nonperiodic flow. J Atmospheric Sci 1963;20(2):130–41. Momani S, Odibat Z. Comparison between homotopy perturbation method and the variational iteration method for linear fractional partial differential equations. Computers Math Appl, in press. Momani S, Odibat Z. Homotopy perturbation method for nonlinear partial differential equations of fractional order. Phys Lett A 2007;365:345–50. Noor MA, Mohyud-Din ST. An efficient algorithm for solving fifth-order boundary value problems. Math Comput Model 2007;45(7-8):954–64. Noorani MSM, Hashim I, Ahmad R, Bakar SA, Ismail ES, Zakaria AM. Comparing numerical methods for the solutions of the Chen system. Chaos Solitons & Fractals 2007;32:1296–304. Odibat Z, Momani S. Modified homotopy perturbation method: Application to quadratic Riccati differential equation of fractional order. Chaos Solitons & Fractals, doi:10.1016/j.chaos.2006.06.041. Shakeri F, Dehghan M. Inverse problem of diffusion equation by He’s homotopy perturbation method. Phys Scr 2007;75(4):551–6. Zhang LN, He JH. Homotopy perturbation method for the solution of the electrostatic potential differential equation. Math Prob Eng 2006:83878.