# Application of Multistage Homotopy Perturbation Method to the Chaotic Genesio System

Authors: M. S. H. Chowdhury; I. Hashim; S. Momani; M. M. Rahman
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

Finding accurate solution of chaotic system by using efficient existing numerical methods is very hard for its complex dynamical behaviors. In this paper, the multistage homotopy-perturbation method (MHPM) is applied to the Chaotic Genesio system. The MHPM is a simple reliable modification based on an adaptation of the standard
homotopy-perturbation method (HPM). The HPM is treated as an algorithm in a sequence of intervals for finding accurate approximate solutions to the Chaotic Genesio system. Numerical comparisons between the MHPM and the classical fourth-order Runge-Kutta (RK4) solutions are made. The results reveal that the new technique is a promising tool for the nonlinear chaotic systems of ordinary differential equations.


---

## Body

## 1. Introduction

Chaos is very interesting nonlinear phenomenon and has been intensively studied in the last three decades. The dynamical systems that exhibit chaotic behavior are sensitive to initial conditions. Chaotic behavior can be found in a variety of systems such as electrical circuits, lasers, fluid dynamics, mechanical devices, time evolution of the magnetic field of celestial bodies, population growth in ecology, the dynamics of molecular vibrations, and not forgetting the weather. The history of chaos theory has come a long way since Jacques Hadamard who in 1898 published a significant study of a free particle gliding frictionlessly on a surface of constant negative curvature which exhibits chaotic motion. The Genesio-Tesi system, proposed by Genesio and Tesi [1], is one of paradigms of chaos since it captures many features of chaotic systems. It includes a simple part and three simple ordinary differential equations (ODEs) that depend on three positive real parameters. The dynamic equation of the system is as follows:(1.1)dxdt=y,dydt=z,dzdt=-cx-by-az+x2,wherex, y, z are state variables, and a, b, c are positive real constants satisfying ab<c. For instance, the system is chaotic for the parameter a=1.2, b=2.92, c=6.In recent years, much attention has been devoted to the application of the HPM [2], to the solutions of various scientific models [3–6]. HPM yields rapidly convergent series solutions [7, 8]. Now, the application of HPM has been extended to the thermal problems [9]. The homotopy perturbation method admits some unknown parameters in the obtained series solutions. In [10], Ji-Huan suggested the least square method can be used to identify the unknown parameters involved in the series solutions.Recently, Chowdhury et al. [11] was the first to successfully apply the multistage homotopy-perturbation method (MHPM) to the chaotic and nonchaotic Lorenz system. The MHPM is also applied to solve the Chen system and a class of systems of ODEs [12–14].In this paper we are again interested in the accuracy of the MHPM for nonlinear systems of ODEs capable of exhibiting chaotic behavior. The system which is of interest to us is the chaotic Genesio system (1.1). Numerical comparisons with the HPM and fourth-order Runge-Kutta (RK4) solutions show that MHPM is accurate and efficient.


## 2. Solution Approaches

We consider a general system of first-order ODEs:(2.1)du1dt+g1(t,u1,u2,…,um)=f1(t),du2dt+g2(t,u1,u2,…,um)=f2(t),⋮dumdt+gm(t,u1,u2,…,um)=fm(t),
subject to the initial conditions:(2.2)u1(t0)=c1,u2(t0)=c2,…,um(t0)=cm.First we write system (2.1) in the operator form:(2.3)L(u1)+N1(u1,u2,…,um)-f1=0,L(u2)+N2(u1,u2,…,um)-f2=0,⋮L(um)+Nm(u1,u2,…,um)-fm=0,subject to the initial conditions (2.2), where L=d/dt is a linear operator and N1,N2,…,Nm are the nonlinear operators. We will next present the solution approaches for (2.3) based on the standard HPM and MHPM separately.


### 2.1.  Solution by HPM

According to HPM, we construct a homotopy for (2.3) which satisfies the following relations:(2.4)L(u1)-L(v1)+pL(v1)+p[N1(u1,u2,…,um)-f1]=0,L(u2)-L(v2)+pL(v2)+p[N2(u1,u2,…,um)-f2]=0,⋮L(um)-L(vm)+pL(vm)+p[Nm(u1,u2,…,um)-fm]=0,wherep∈[0,1] is an embedding parameter and v1,v2,…,vm are initial approximations satisfying the given conditions. It is obvious that when the perturbation parameter p=0, (2.4) become a linear system and when p=1 we get the original nonlinear system.Let us take the initial approximations as follows:(2.5)u1(t)=u1,0(t)+pu1,1(t)+p2u1,2(t)+p3u1,3(t)+⋯,u2(t)=u2,0(t)+pu2,1(t)+p2u2,2(t)+p3u2,3(t)+⋯,⋮um(t)=um,0(t)+pum,1(t)+p2um,2(t)+p3um,3(t)+⋯,u1,0(t)=v1(t)=u1(t0)=c1,u2,0(t)=v2(t)=u2(t0)=c2,⋮um,0(t)=vm(t)=um(t0)=cm,whereui,j, (i=1,2,…,m;j=1,2,…) are functions yet to be determined. Substituting (2.5) into (2.4) and arranging the coefficients of the same powers of p, we get(2.6)L(u1,1)+L(v1)+N1(u1,0,u2,0,…,um,0)-f1=0,u1,1(t0)=0,L(u2,1)+L(v2)+N2(u1,0,u2,0,…,um,0)-f2=0,u2,1(t0)=0,⋮L(um,1)+L(vm)+Nm(u1,0,u2,0,…,um,0)-fm=0,um,1(t0)=0,L(u1,2)+N1(u1,1,u2,1,…,um,1)=0,u1,2(t0)=0,L(u2,2)+N2(u1,1,u2,1,…,um,1)=0,u2,2(t0)=0⋮L(um,2)+Nm(u1,1,u2,1,…,um,1)=0,um,2(t0)=0,and so forth. We solve the above systems of equations for the unknownsui,j(i=1,2,…,m;j=1,2,…) by applying the inverse operator(2.7)L-1(⋅)=∫t0t(⋅)dt.Therefore, according to HPM then-term approximations for the solutions of (2.3) can be expressed as(2.8)ϕ1,n(t)=u1(t)=lim⁡p→1u1(t)=∑k=0n-1u1,k(t),ϕ2,n(t)=u2(t)=lim⁡p→1u2(t)=∑k=0n-1u2,k(t),⋮ϕm,n(t)=um(t)=lim⁡p→1um(t)=∑k=0n-1um,k(t).


### 2.2. Solution by MHPM

The approximate solutions (2.8) are generally, as will 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 (2.6) as an algorithm for approximating the solutions of (2.1) in a sequence of intervals choosing the initial approximations as(2.9)u1,0(t)=v1(t)=u1(t*)=c1*,u2,0(t)=v2(t)=u2(t*)=c2*,⋮um,0(t)=vm(t)=um(t*)=cm*,wheret* is the left-end point of each subinterval.Now we solve (2.6) for the unknowns ui,j(i=1,2,…,m;j=1,2,…) by applying the inverse linear operator(2.10)L-1(⋅)=∫t*t(⋅)dt.In order to carry out the iterations in every subinterval of equal lengthΔt, [0,t1), [t1,t2), [t2,t3)⋯[tj-1,t), we need to know the values of the following:(2.11)u1,0*(t)=u1(t*),u2,0*(t)=u2(t*),…,um,0*(t)=um(t*).But, in general, we do not have these information at our clearance except at the initial pointt*=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 given by (2.8), that is,(2.12)u1,0*≃ϕ1,n(t*),u2,0*≃ϕ2,n(t*),…,um,0*≃ϕm,n(t*).


## 2.1.  Solution by HPM

According to HPM, we construct a homotopy for (2.3) which satisfies the following relations:(2.4)L(u1)-L(v1)+pL(v1)+p[N1(u1,u2,…,um)-f1]=0,L(u2)-L(v2)+pL(v2)+p[N2(u1,u2,…,um)-f2]=0,⋮L(um)-L(vm)+pL(vm)+p[Nm(u1,u2,…,um)-fm]=0,wherep∈[0,1] is an embedding parameter and v1,v2,…,vm are initial approximations satisfying the given conditions. It is obvious that when the perturbation parameter p=0, (2.4) become a linear system and when p=1 we get the original nonlinear system.Let us take the initial approximations as follows:(2.5)u1(t)=u1,0(t)+pu1,1(t)+p2u1,2(t)+p3u1,3(t)+⋯,u2(t)=u2,0(t)+pu2,1(t)+p2u2,2(t)+p3u2,3(t)+⋯,⋮um(t)=um,0(t)+pum,1(t)+p2um,2(t)+p3um,3(t)+⋯,u1,0(t)=v1(t)=u1(t0)=c1,u2,0(t)=v2(t)=u2(t0)=c2,⋮um,0(t)=vm(t)=um(t0)=cm,whereui,j, (i=1,2,…,m;j=1,2,…) are functions yet to be determined. Substituting (2.5) into (2.4) and arranging the coefficients of the same powers of p, we get(2.6)L(u1,1)+L(v1)+N1(u1,0,u2,0,…,um,0)-f1=0,u1,1(t0)=0,L(u2,1)+L(v2)+N2(u1,0,u2,0,…,um,0)-f2=0,u2,1(t0)=0,⋮L(um,1)+L(vm)+Nm(u1,0,u2,0,…,um,0)-fm=0,um,1(t0)=0,L(u1,2)+N1(u1,1,u2,1,…,um,1)=0,u1,2(t0)=0,L(u2,2)+N2(u1,1,u2,1,…,um,1)=0,u2,2(t0)=0⋮L(um,2)+Nm(u1,1,u2,1,…,um,1)=0,um,2(t0)=0,and so forth. We solve the above systems of equations for the unknownsui,j(i=1,2,…,m;j=1,2,…) by applying the inverse operator(2.7)L-1(⋅)=∫t0t(⋅)dt.Therefore, according to HPM then-term approximations for the solutions of (2.3) can be expressed as(2.8)ϕ1,n(t)=u1(t)=lim⁡p→1u1(t)=∑k=0n-1u1,k(t),ϕ2,n(t)=u2(t)=lim⁡p→1u2(t)=∑k=0n-1u2,k(t),⋮ϕm,n(t)=um(t)=lim⁡p→1um(t)=∑k=0n-1um,k(t).


## 2.2. Solution by MHPM

The approximate solutions (2.8) are generally, as will 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 (2.6) as an algorithm for approximating the solutions of (2.1) in a sequence of intervals choosing the initial approximations as(2.9)u1,0(t)=v1(t)=u1(t*)=c1*,u2,0(t)=v2(t)=u2(t*)=c2*,⋮um,0(t)=vm(t)=um(t*)=cm*,wheret* is the left-end point of each subinterval.Now we solve (2.6) for the unknowns ui,j(i=1,2,…,m;j=1,2,…) by applying the inverse linear operator(2.10)L-1(⋅)=∫t*t(⋅)dt.In order to carry out the iterations in every subinterval of equal lengthΔt, [0,t1), [t1,t2), [t2,t3)⋯[tj-1,t), we need to know the values of the following:(2.11)u1,0*(t)=u1(t*),u2,0*(t)=u2(t*),…,um,0*(t)=um(t*).But, in general, we do not have these information at our clearance except at the initial pointt*=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 given by (2.8), that is,(2.12)u1,0*≃ϕ1,n(t*),u2,0*≃ϕ2,n(t*),…,um,0*≃ϕm,n(t*).


## 3. Application

In this section,we will study the Genesio system (1.1) subject to the initial conditions:(3.1)x(0)=c1,y(0)=c2,z(0)=c3.According to the HPM, we can construct a homotopy which satisfies the following relation:(3.2)v1′-x0′+p(x0′-v2)=0,v2′-y0′+p(y0′-v3)=0,v3′-z0′+p(z0′+cv1+bv2+av3-v12)=0.We take the initial approximations as:(3.3)v1(t)=v1,0(t)+pv1,1(t)+p2v1,2(t)+p3v1,3(t)+⋯,v2(t)=v2,0(t)+pv2,1(t)+p2v2,2(t)+p3v2,3(t)+⋯,v3(t)=v3,0(t)+pv3,1(t)+p2v3,2(t)+p3v3,3(t)+⋯,v1,0(t)=x0(t)=x(t*)=c1,v2,0(t)=y0(t)=y(t*)=c2,v3,0(t)=z0(t)=z(t*)=c3,wherevi,j,  i,j=1,2,3,… are functions yet to be determined. Substituting (3.3) into (3.2) and collecting terms the same powers of p, we have(3.4)v1,1′-v2,0=0,v2,1′-v3,0=0,v3,1′+cv1,0+bv2,0+av3,0-v1,02=0,v1,2′-v2,1=0,v2,2′-v3,1=0,v3,2′+cv1,1+bv2,1+av3,1-2v1,0v1,1=0,v1,3′-v2,2=0,v2,3′-v3,2=0,v3,3′+cv1,2+bv2,2+av3,2-2v1,0v1,2-v1,12=0.In order to obtain the unknownsvi,j(t),  i,j=1,2,3, we solve the above system taking the initial conditions vi,j(0)=0,  i,j=1,2,3, we obtain,(3.5)v1,1(t)=c2(t-t*),v2,1(t)=c3(t-t*),v3,1(t)=(c12-cc1-bc2-ac3)(t-t*),v1,2(t)=12c3(t-t*)2,v2,2(t)=12[c12-cc1-bc2-ac3](t-t*)2,v3,2(t)=12[2c1c2-cc2-bc3-ac12+acc1+abc2+a2c3](t-t*)2,v1,3(t)=16[c12-cc1-bc2-ac3](t-t*)3,v2,3(t)=16[2c1c2-cc2-bc3-ac12+acc1+abc2+a2c3](t-t*)3,v3,3(t)=16[2c22+bcc1+2abc3-a2cc1-a2bc2-2ac1c2+acc2+2c1c3-bc12+b2c2+a2c12-a3c3-cc3](t-t*)3.Hence, the solution to the Genesio system (1.1) is:(3.6)x=∑m=0∞v1,m(t),y=∑m=0∞v2,m(t),z=∑m=0∞v3,m(t).
To carry out the iterations on every subinterval of equal length Δt, we need to know the values of the following initial conditions:(3.7)c1=x(t*),c2=y(t*),c3=z(t*).In general, we do not have these information at our clearance except at the initial pointt*=t0=0, but we can obtain these values following the MHPM as given in Section 2.2. We note that the 10-term approximations of x, y, and z are denoted as x(t)⋍ϕ10(t)=∑i=09v1,i, y(t)⋍ψ10(t)=∑i=09v2,i, and z(t)⋍ξ10(t)=∑i=09v3,i.


## 4. Results and Discussions

The MHPM algorithm is coded in the computer algebra package Maple and we employ the Maple’s built-in fourth-order Runge-Kutta procedure (RK4). The Maple environment variable digits controlling the number of significant digits are set to 16 in all the calculations done in this paper. For the comparison, we set the parametera=1.2, b=2.92, and c=6 where the system exhibits chaotic behavior, alongside with its initial conditions x(0)=0.2, y(0)=-0.3, and z(0)=0.1. The time range studied in this work is t=0 to 20. In general, there is not a known exact solution for the Genesio system. Thus, the accuracy of the present method will be determined by its comparison to the numerical solution by Runge-kutta method. We determine the accuracy of RK4 for the solution of (1.1) for different time steps. From the results presented in Table 1 we see that the maximum difference between the RK4 solutions on time steps Δt=0.001 and Δt=0.0001 is of the order of magnitude of 10-11. So based on these observations we choose the RK4 solutions on the time step Δt=0.001 as the benchmark for our comparison purposes. So now we can compare the accuracy of the HPM and MHPM with the RK4 method on the chosen time step Δt=0.001. We choose this time step since a smaller one is computationally costly. We note that increasing the number of terms in the series solutions (1.1) improves the accuracy of the MHPM solutions but at the expense of increased computational efforts. In this work we fix the number of terms used to be ten. The details of the differences between the 10-term HPM solutions and 10-term MHPM solutions (Δt=0.001) and the RK4 solutions on Δt=0.001 are given in Table 2. In Table 2, we see that the 10-term MHPM solutions on the time step Δt=0.001 agree with the RK4 solutions at least 10 decimal places, while the 10-term classical HPM solutions are only valid for t≈1. For the chaotic Genesio system we observe that the MHPM has the advantage over the RK4 on achieving a good accuracy with a larger time step.Table 1
A determination of the accuracy of RK4 for the Genesio system (1.1).
tΔ=|RK40.01-RK40.001|Δ=|RK40.001-RK40.0001|ΔxΔyΔzΔxΔyΔz12.130E-101.525E-106.681E-102.12E-141.58E-146.68E-1423.299E-108.467E-109.038E-103.36E-148.40E-149.30E-1434.957E-101.443E-092.394E-094.85E-141.461E-132.362E-1341.533E-099.664E-105.188E-091.540E-139.27E-145.232E-1353.154E-104.876E-091.175E-093.51E-144.884E-131.044E-1363.474E-092.279E-091.324E-083.455E-132.380E-131.322E-1273.567E-098.028E-099.569E-093.631E-137.955E-139.820E-1383.430E-091.243E-081.969E-083.341E-131.257E-121.942E-1291.132E-085.660E-093.797E-081.135E-125.373E-133.824E-12101.655E-092.685E-082.350E-091.881E-132.688E-121.640E-13111.674E-081.348E-087.314E-081.667E-121.402E-127.308E-12122.241E-084.278E-086.208E-082.270E-124.234E-126.351E-12131.873E-085.370E-086.170E-081.822E-125.436E-126.050E-12143.733E-082.187E-081.333E-073.770E-122.083E-121.348E-11151.595E-081.228E-071.926E-081.679E-121.232E-111.565E-12167.856E-084.718E-082.002E-077.808E-124.950E-121.999E-11171.317E-087.821E-089.779E-081.510E-127.839E-121.008E-11182.585E-088.743E-081.959E-072.601E-129.132E-121.983E-11199.424E-086.053E-082.934E-079.655E-125.883E-123.032E-11205.455E-084.084E-081.401E-075.290E-124.470E-121.423E-11Table 2
Differences between 10-term HPM and 10-term MHPM with RK4 solutions onΔt=0.001.
tΔ=|HPM-RK40.001|Δ=|MHPM0.001-RK40.001|ΔxΔyΔzΔxΔyΔz11.856E-056.289E-050.00023552.12E-141.58E-146.68E-1420.016260.045390.22113.36E-148.40E-149.30E-1430.88151.91210.734.85E-141.461E-132.362E-13414.8328.01154.91.540E-139.27E-145.232E-135126.8236.311933.51E-144.885E-131.044E-136703.5137362933.456E-132.380E-131.322E-127291060632.573e+043.631E-137.956E-139.821E-13897902.174E+048.736E+043.341E-131.257E-121.942E-1292.827E+046.649E+042.570E+051.135E-125.373E-133.824E-12107.263E+041.795E+056.751E+051.881E-132.688E-121.640E-13111.700E+054.386E+051.617E+061.667E-121.403E-127.308E-12123.687E+059.877E+053.587E+062.270E-124.234E-126.352E-12137.510E+052.078E+067.464E+061.822E-125.436E-126.050E-12141.450E+064.130E+061.470E+073.771E-122.083E-121.348E-11152.675E+067.813E+062.763E+071.679E-121.232E-111.565E-12164.741E+061.417E+074.983E+077.809E-124.950E-121.999E-11178.117E+062.475E+078.667E+071.510E-127.840E-121.008E-11181.348E+074.184E+071.460E+082.601E-129.133E-121.983E-11192.177E+076.870E+072.391E+089.656E-125.884E-123.032E-11203.431E+071.099E+083.817E+085.291E-124.470E-121.424E-11In Figure1 we reproduce the well-known x-y, x-z, y-z, and x-y-z phase portraits of the chaotic Genesio system using the 10-term MHPM solutions on Δt=0.001. Having determined the benchmark timestep, we can now investigate the accuracy of the HPM and MHPM as compared to RK4. In Figure 2 we plot the 10-term HPM solutions and 10-term MHPM solutions (on Δt=0.001) against the RK4 solutions on Δt=0.001. Both the 10-term MHPM solutions on Δt=0.001 and RK4 solutions on Δt=0.001 seem to overlap on the scale used in figure, but the numerical results from the standard HPM start to stay away at about t≈2 for x(t), y(t), and z(t).Phase portraits using 10-term MHPM.Graphical comparisons between HPM, MHPM (Δt=0.001), and RK4 (Δt=0.001) for x(t), y(t), and z(t).


## 5. Conclusions

In this work, the MHPM was applied to the solutions of the well-known Genesio system. The MHPM is only a simple modification of the standard HPM. Comparisons between the HPM and MHPM solutions and the fourth-order Runge-Kutta (RK4) numerical solutions were made. For the chaotic Genesio system 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. We note that the MHPM solutions were computed via a simple algorithm with less amount of computations and without any need for perturbation techniques, special transformations, linearization, or discretization.



---

*Source: 974293-2012-05-13.xml*