Application of Multistage Homotopy Perturbation Method to the Chaotic Genesio System

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.

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:d𝑥d𝑡=𝑦,d𝑦d𝑡=𝑧,d𝑧d𝑡=−𝑐𝑥−𝑏𝑦−𝑎𝑧+𝑥², (1.1)

where 𝑥, 𝑦, 𝑧 are state variables, and 𝑎, 𝑏, 𝑐 are positive real constants satisfying 𝑎𝑏<𝑐. For instance, the system is chaotic for the parameter 𝑎=1.2, 𝑏=2.92, 𝑐=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:d𝑢₁d𝑡+𝑔₁𝑡,𝑢₁,𝑢₂,…,𝑢_𝑚=𝑓₁(𝑡),d𝑢₂d𝑡+𝑔₂𝑡,𝑢₁,𝑢₂,…,𝑢_𝑚=𝑓₂(𝑡),⋮d𝑢_𝑚d𝑡+𝑔_𝑚𝑡,𝑢₁,𝑢₂,…,𝑢_𝑚=𝑓_𝑚(𝑡), (2.1) subject to the initial conditions:𝑢₁𝑡₀=𝑐₁,𝑢₂𝑡₀=𝑐₂,…,𝑢_𝑚𝑡₀=𝑐_𝑚. (2.2)

First we write system (2.1) in the operator form:𝐿𝑢₁+𝑁₁𝑢₁,𝑢₂,…,𝑢_𝑚−𝑓₁𝐿𝑢=0,₂+𝑁₂𝑢₁,𝑢₂,…,𝑢_𝑚−𝑓₂⋮𝐿𝑢=0,_𝑚+𝑁_𝑚𝑢₁,𝑢₂,…,𝑢_𝑚−𝑓_𝑚=0, (2.3)subject to the initial conditions (2.2), where 𝐿=d/d𝑡 is a linear operator and 𝑁₁,𝑁₂,…,𝑁_𝑚 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:𝐿𝑢₁−𝐿𝑣₁+𝑝𝐿𝑣₁+𝑝𝑁₁𝑢₁,𝑢₂,…,𝑢_𝑚−𝑓₁𝐿𝑢=0,₂−𝐿𝑣₂+𝑝𝐿𝑣₂+𝑝𝑁₂𝑢₁,𝑢₂,…,𝑢_𝑚−𝑓₂⋮𝐿𝑢=0,_𝑚−𝐿𝑣_𝑚+𝑝𝐿𝑣_𝑚+𝑝𝑁_𝑚𝑢₁,𝑢₂,…,𝑢_𝑚−𝑓_𝑚=0, (2.4)

where 𝑝∈[0,1] is an embedding parameter and 𝑣₁,𝑣₂,…,𝑣_𝑚 are initial approximations satisfying the given conditions. It is obvious that when the perturbation parameter 𝑝=0, (2.4) become a linear system and when 𝑝=1 we get the original nonlinear system.

Let us take the initial approximations as follows:𝑢₁(𝑡)=𝑢_1,0(𝑡)+𝑝𝑢_1,1(𝑡)+𝑝²𝑢_1,2(𝑡)+𝑝³𝑢_1,3(𝑡)𝑢+⋯,₂(𝑡)=𝑢_2,0(𝑡)+𝑝𝑢_2,1(𝑡)+𝑝²𝑢_2,2(𝑡)+𝑝³𝑢_2,3(𝑡)⋮𝑢+⋯,_𝑚(𝑡)=𝑢_𝑚,0(𝑡)+𝑝𝑢_𝑚,1(𝑡)+𝑝²𝑢_𝑚,2(𝑡)+𝑝³𝑢_𝑚,3(𝑡)𝑢+⋯,_1,0(𝑡)=𝑣₁(𝑡)=𝑢₁𝑡₀=𝑐₁,𝑢_2,0(𝑡)=𝑣₂(𝑡)=𝑢₂𝑡₀=𝑐₂,⋮𝑢_𝑚,0(𝑡)=𝑣_𝑚(𝑡)=𝑢_𝑚𝑡₀=𝑐_𝑚, (2.5)

where 𝑢_𝑖,𝑗, (𝑖=1,2,…,𝑚;𝑗=1,2,…) are functions yet to be determined. Substituting (2.5) into (2.4) and arranging the coefficients of the same powers of 𝑝, we get𝐿𝑢_1,1+𝐿𝑣₁+𝑁₁𝑢_1,0,𝑢_2,0,…,𝑢_𝑚,0−𝑓₁=0,𝑢_1,1𝑡₀𝐿𝑢=0,_2,1+𝐿𝑣₂+𝑁₂𝑢_1,0,𝑢_2,0,…,𝑢_𝑚,0−𝑓₂=0,𝑢_2,1𝑡₀⋮𝐿𝑢=0,_𝑚,1+𝐿𝑣_𝑚+𝑁_𝑚𝑢_1,0,𝑢_2,0,…,𝑢_𝑚,0−𝑓_𝑚=0,𝑢_𝑚,1𝑡₀𝐿𝑢=0,_1,2+𝑁₁𝑢_1,1,𝑢_2,1,…,𝑢_𝑚,1=0,𝑢_1,2𝑡₀𝐿𝑢=0,_2,2+𝑁₂𝑢_1,1,𝑢_2,1,…,𝑢_𝑚,1=0,𝑢_2,2𝑡₀⋮𝐿𝑢=0_𝑚,2+𝑁_𝑚𝑢_1,1,𝑢_2,1,…,𝑢_𝑚,1=0,𝑢_𝑚,2𝑡₀=0, (2.6)

and so forth. We solve the above systems of equations for the unknowns 𝑢_𝑖,𝑗(𝑖=1,2,…,𝑚;𝑗=1,2,…) by applying the inverse operator𝐿⁻¹(⋅)=^𝑡_𝑡0(⋅)d𝑡. (2.7)

Therefore, according to HPM the 𝑛-term approximations for the solutions of (2.3) can be expressed as𝜙_1,𝑛(𝑡)=𝑢₁(𝑡)=lim_𝑝→1𝑢₁(𝑡)=^𝑛−1_𝑘=0𝑢_1,𝑘(𝑡),𝜙_2,𝑛(𝑡)=𝑢₂(𝑡)=lim_𝑝→1𝑢₂(𝑡)=^𝑛−1_𝑘=0𝑢_2,𝑘(𝑡),⋮𝜙_𝑚,𝑛(𝑡)=𝑢_𝑚(𝑡)=lim_𝑝→1𝑢_𝑚(𝑡)=^𝑛−1_𝑘=0𝑢_𝑚,𝑘(𝑡). (2.8)

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 𝑡. A simple way of ensuring validity of the approximations for large 𝑡 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𝑢_1,0(𝑡)=𝑣₁(𝑡)=𝑢₁𝑡^∗=𝑐^∗₁,𝑢_2,0(𝑡)=𝑣₂(𝑡)=𝑢₂𝑡^∗=𝑐^∗₂,⋮𝑢_𝑚,0(𝑡)=𝑣_𝑚(𝑡)=𝑢_𝑚𝑡^∗=𝑐^∗_𝑚, (2.9)

where 𝑡^∗ is the left-end point of each subinterval.

Now we solve (2.6) for the unknowns 𝑢_𝑖,𝑗(𝑖=1,2,…,𝑚;𝑗=1,2,…) by applying the inverse linear operator𝐿⁻¹(⋅)=^𝑡_𝑡∗(⋅)d𝑡. (2.10)

In order to carry out the iterations in every subinterval of equal length Δ𝑡, [0,𝑡₁), [𝑡₁,𝑡₂), [𝑡₂,𝑡₃)⋯[𝑡_𝑗−1,𝑡), we need to know the values of the following:𝑢^∗_1,0(𝑡)=𝑢₁𝑡^∗,𝑢^∗_2,0(𝑡)=𝑢₂𝑡^∗,…,𝑢^∗_𝑚,0(𝑡)=𝑢_𝑚𝑡^∗. (2.11)

But, in general, we do not have these information at our clearance except at the initial point 𝑡^∗=𝑡₀. A simple way for obtaining the necessary values could be by means of the previous 𝑛-term approximations 𝜙_1,𝑛,𝜙_2,𝑛,…,𝜙_𝑚,𝑛 of the preceding subinterval given by (2.8), that is,𝑢^∗_1,0≃𝜙_1,𝑛𝑡^∗,𝑢^∗_2,0≃𝜙_2,𝑛𝑡^∗,…,𝑢^∗_𝑚,0≃𝜙_𝑚,𝑛𝑡^∗. (2.12)

3. Application

In this section,we will study the Genesio system (1.1) subject to the initial conditions:𝑥(0)=𝑐₁,𝑦(0)=𝑐₂,𝑧(0)=𝑐₃. (3.1)

According to the HPM, we can construct a homotopy which satisfies the following relation:𝑣^′₁−𝑥^′₀+𝑝𝑥^′₀−𝑣₂𝑣=0,_′2−𝑦^′₀+𝑝𝑦^′₀−𝑣₃𝑣=0,_′3−𝑧^′₀+𝑝𝑧^′₀+𝑐𝑣₁+𝑏𝑣₂+𝑎𝑣₃−𝑣²₁=0. (3.2)

We take the initial approximations as:𝑣₁(𝑡)=𝑣_1,0(𝑡)+𝑝𝑣_1,1(𝑡)+𝑝²𝑣_1,2(𝑡)+𝑝³𝑣_1,3(𝑡)𝑣+⋯,₂(𝑡)=𝑣_2,0(𝑡)+𝑝𝑣_2,1(𝑡)+𝑝²𝑣_2,2(𝑡)+𝑝³𝑣_2,3(𝑡)𝑣+⋯,₃(𝑡)=𝑣_3,0(𝑡)+𝑝𝑣_3,1(𝑡)+𝑝²𝑣_3,2(𝑡)+𝑝³𝑣_3,3(𝑡)𝑣+⋯,_1,0(𝑡)=𝑥₀(𝑡)=𝑥𝑡^∗=𝑐₁,𝑣_2,0(𝑡)=𝑦₀(𝑡)=𝑦𝑡^∗=𝑐₂,𝑣_3,0(𝑡)=𝑧₀(𝑡)=𝑧𝑡^∗=𝑐₃, (3.3)

where 𝑣_𝑖,𝑗,??𝑖,𝑗=1,2,3,… are functions yet to be determined. Substituting (3.3) into (3.2) and collecting terms the same powers of 𝑝, we have𝑣^′_1,1−𝑣_2,0𝑣=0,_′2,1−𝑣_3,0𝑣=0,_′3,1+𝑐𝑣_1,0+𝑏𝑣_2,0+𝑎𝑣_3,0−𝑣²_1,0𝑣=0,_′1,2−𝑣_2,1𝑣=0,_′2,2−𝑣_3,1𝑣=0,_′3,2+𝑐𝑣_1,1+𝑏𝑣_2,1+𝑎𝑣_3,1−2𝑣_1,0𝑣_1,1𝑣=0,_′1,3−𝑣_2,2𝑣=0,_′2,3−𝑣_3,2𝑣=0,_′3,3+𝑐𝑣_1,2+𝑏𝑣_2,2+𝑎𝑣_3,2−2𝑣_1,0𝑣_1,2−𝑣²_1,1=0. (3.4)

In order to obtain the unknowns 𝑣_𝑖,𝑗(𝑡),??𝑖,𝑗=1,2,3, we solve the above system taking the initial conditions 𝑣_𝑖,𝑗(0)=0,??𝑖,𝑗=1,2,3, we obtain,𝑣_1,1(𝑡)=𝑐₂𝑡−𝑡^∗,𝑣_2,1(𝑡)=𝑐₃𝑡−𝑡^∗,𝑣_3,1(𝑡)=𝑐²₁−𝑐𝑐₁−𝑏𝑐₂−𝑎𝑐₃𝑡−𝑡^∗,𝑣_1,2(𝑡)=12𝑐₃𝑡−𝑡^∗2,𝑣_2,2(𝑡)=12𝑐²₁−𝑐𝑐₁−𝑏𝑐₂−𝑎𝑐₃𝑡−𝑡^∗2,𝑣_3,2(𝑡)=122𝑐₁𝑐₂−𝑐𝑐₂−𝑏𝑐₃−𝑎𝑐²₁+𝑎𝑐𝑐₁+𝑎𝑏𝑐₂+𝑎²𝑐₃𝑡−𝑡^∗2,𝑣_1,3(𝑡)=16𝑐²₁−𝑐𝑐₁−𝑏𝑐₂−𝑎𝑐₃𝑡−𝑡^∗3,𝑣_2,3(𝑡)=162𝑐₁𝑐₂−𝑐𝑐₂−𝑏𝑐₃−𝑎𝑐²₁+𝑎𝑐𝑐₁+𝑎𝑏𝑐₂+𝑎²𝑐₃𝑡−𝑡^∗3,𝑣_3,3(𝑡)=162𝑐²₂+𝑏𝑐𝑐₁+2𝑎𝑏𝑐₃−𝑎²𝑐𝑐₁−𝑎²𝑏𝑐₂−2𝑎𝑐₁𝑐₂+𝑎𝑐𝑐₂+2𝑐₁𝑐₃−𝑏𝑐²₁+𝑏²𝑐₂+𝑎²𝑐²₁−𝑎³𝑐₃−𝑐𝑐₃𝑡−𝑡^∗3. (3.5)

Hence, the solution to the Genesio system (1.1) is:𝑥=^∞_𝑚=0𝑣_1,𝑚(𝑡),𝑦=∑^∞_𝑚=0𝑣_2,𝑚(𝑡),𝑧=∑^∞_𝑚=0𝑣_3,𝑚(𝑡). (3.6) To carry out the iterations on every subinterval of equal length Δ𝑡, we need to know the values of the following initial conditions:𝑐₁=𝑥𝑡^∗,𝑐₂=𝑦𝑡^∗,𝑐₃=𝑧𝑡^∗. (3.7)

In general, we do not have these information at our clearance except at the initial point 𝑡^∗=𝑡₀=0, but we can obtain these values following the MHPM as given in Section 2.2. We note that the 10-term approximations of 𝑥, 𝑦, and 𝑧 are denoted as 𝑥(𝑡)⋍𝜙₁₀(𝑡)=∑⁹_𝑖=0𝑣_1,𝑖, 𝑦(𝑡)⋍𝜓₁₀(𝑡)=∑⁹_𝑖=0𝑣_2,𝑖, and 𝑧(𝑡)⋍𝜉₁₀(𝑡)=∑⁹_𝑖=0𝑣_3,𝑖.

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 parameter 𝑎=1.2, 𝑏=2.92, and 𝑐=6 where the system exhibits chaotic behavior, alongside with its initial conditions 𝑥(0)=0.2, 𝑦(0)=−0.3, and 𝑧(0)=0.1. The time range studied in this work is 𝑡=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 Δ𝑡=0.001 and Δ𝑡=0.0001 is of the order of magnitude of 10⁻¹¹. So based on these observations we choose the RK4 solutions on the time step Δ𝑡=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 Δ𝑡=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 (Δ𝑡=0.001) and the RK4 solutions on Δ𝑡=0.001 are given in Table 2. In Table 2, we see that the 10-term MHPM solutions on the time step Δ𝑡=0.001 agree with the RK4 solutions at least 10 decimal places, while the 10-term classical HPM solutions are only valid for 𝑡≈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).

Table 2 

Differences between 10-term HPM and 10-term MHPM with RK4 solutions on Δ𝑡=0.001.

In Figure 1 we reproduce the well-known 𝑥-𝑦, 𝑥-𝑧, 𝑦-𝑧, and 𝑥-𝑦-𝑧 phase portraits of the chaotic Genesio system using the 10-term MHPM solutions on Δ𝑡=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 Δ𝑡=0.001) against the RK4 solutions on Δ𝑡=0.001. Both the 10-term MHPM solutions on Δ𝑡=0.001 and RK4 solutions on Δ𝑡=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 𝑡≈2 for 𝑥(𝑡), 𝑦(𝑡), and 𝑧(𝑡).

Figure 1 

Phase portraits using 10-term MHPM.

Figure 2 

Graphical comparisons between HPM, MHPM (Δ𝑡=0.001), and RK4 (Δ𝑡=0.001) for 𝑥(𝑡), 𝑦(𝑡), and 𝑧(𝑡).

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.

Acknowledgments

The authors would like to acknowledge the financial supports received from the Ministry of Higher Education via the FRGS Grant no. FRGS0409-106 and Universiti Putra Malaysia.