Optimal Homotopy Asymptotic Method for an Anharmonic Oscillator: Application to the Chen System

Abstract

The aim of our work is to obtain the analytic solutions for a new nonlinear anharmonic oscillator by means of the Optimal Homotopy Asymptotic Method (OHAM), using only one iteration. The accuracy of the obtained results comes from the comparison with the corresponding numerical ones for specified physical parameters. Moreover, the OHAM method has a greater degree of flexibility than an iterative method as is presented in this paper. Based on these results, the analytically solutions of the Chen system were obtained for a special case (just one analytic first integral). The chaotic behaviors were excluded here. The provided solutions are usefully for many engineering applications.

1. Introduction

Many applications in electrical engineering, medicine, outdoor weather control application, secure communication techniques, and so on are based by the study of chaotic dynamic systems. Therefore, many techniques have been developed to study the perturbation of this simple periodic motion.

Even if the one dimensional harmonic oscillator potential is a suitable model for a series of physical problems, in many cases the concordance with the experimental values is not sufficient. A typical example is the diatomic molecule whose motions, especially vibrations, cannot be described sufficiently by using the one dimensional harmonic oscillator potential. Among many models of anharmonic potentials a privileged position for the potential energy of a diatomic molecule is occupied by the Morse oscillator potential due to its applications in quantum mechanics [1] to diatomic or polyatomic molecules, spectroscopy, and so on. This is a convenient model because it explicitly includes the effects of bond breaking, as well as the existence of unbound states, the anharmonicity of real bonds and bond dissociation. A modified version of anharmonic potential is used in different applications in nonlinear dynamical systems to solve the inverse scatting problem to derive so called soliton solutions [2].

Recently, He et al. [3] adopted a generalized Duffing oscillator using the homotopy perturbation method (HPM). Bel et al. [4] analyzed the synchronization and stability of coupled driven-damped Helmholtz-Duffing oscillators in bi-stability regimes. Synchronization of a nonlinear oscillator was explored by Vieira et al. [5]. Mariano et al. [6] considered a general class of nonlinear dynamical systems with memory. Furthermore, nonlinear oscillators equipped with fast varying periodic time delay feedback were developed by Liu et al. [7] using an analytical criterion. The mathematical techniques as periodic perturbations, stability, chaotic and asymptotic behaviors, and geometrical properties were developed in [8,9,10,11,12,13,14,15,16,17,18,19].

The Chen system [20] proposed in 1999, appears in a variety of studies, such as: Toda Lattice [21], Kowalevski top dynamics [22], Lotka–Volterra system [23], the battery model [24], Lagrange system [25], and others [26].

The relevance of the Chen chaotic system results from several approaches. For instance, an output feedback control algorithm for a single-input single-output variant was proposed in [27], the influence of the time delay was investigated in [28] showing that the single-scroll attractor is indeed chaotic, the global boundedness was explored in [29] using a suitable Lyapunov function, the equilibrium point stabilization problem by employing a simple linear feedback controller was discusses in [30], and others [31,32,33,34,35,36,37,38,39].

The aim of this work is to obtain the approximate closed-form solutions of the Chen system for a special case: just one analytic first integral known. For this case the Chen system could be reduced to an anharmonic oscillator which characterizes the behaviour of more physical systems of interest.

This paper is organized as follows. We focus on the analytical approaches of the approximate closed-form solutions for a special case of the Chen chaotic system in Section 2. In Section 3 we firstly proposed the differential equation which describes the nonlinear anharmonic oscillator. The OHAM method developed in [40,41,42] is applied for obtaining an approximate analytic solutions. Here, the study of the Chen chaotic system is reduced to a nonlinear anharmonic oscillator. Section 4 is devoted to the numerical simulations. There is an excellent agreement between the analytic approximate solution and the corresponding numerical solution which proved the validity of all obtained results. The last Section 5 is dedicated to some concluding remarks.

2. Approximate Closed-Form Solutions to the Chen’s System

The Chen system was analytically solved by mathematical techniques as the multistage homotopy–perturbation method in [33], the multistage homotopy analysis method [34], the differential transformation method (DTM) [36], the Adomian decomposition method (ADM), [37], the mechanical analysis [43], and the geometrical frame [44].

The Chen’s system has the form:

$$\left\{\begin{array}{c}\stackrel{.}{\mathrm{x}}=a\left(y-x\right)\hfill \\ \stackrel{.}{\mathrm{y}}=\left(c-a\right)x-xz+cy,\hfill \\ \stackrel{.}{\mathrm{z}}=-bz+xy\hfill \end{array}\right.$$

(1)

where $a,b,c$ are real parameters. The system is chaotic when $a=35,b=3,c=28$ ([45]).

Remark 1.

(a) 

If $a=0$ and $b=c$ the system (1) has the Hamilton–Poisson realization. The functional $H(x,y,z)=x$ is the Hamiltonian and the $C(x,y,z)=\frac{1}{2}(-2cxz+2yz+x{y}^2+x{z}^2)$ is functionally independent Casimir. Thus, the exact solution is written as an intersection between the surfaces $H=cst$ and $C=cst$.

(b) 

For $a\in R^{*},b=c=0$, the Hamiltonian function is $H(x,y,z)=y^2+z^2+2az$, but finding the Casimir functions remains an open problem. Therefore, it is impossible to write the exact solution as an intersection between the surfaces $H=cst$ and $C=cst$.

(c) 

Otherwise, the Chen’s system is chaotic.

In the following we focus just on the case $a\in R^{*},b=c=0$.

The system (1) becomes:

$$\left\{\begin{array}{c}\dot{x}=a(y-x)\hfill \\ \dot{y}=-ax-xz\hfill \\ \dot{z}=xy\hfill \end{array}\right..$$

(2)

Remark 2.

The considered system admits a symmetry with respect to $Oz$- axis, for $a\in R^{*},$$b=c=0$.

Considering the following representation:

$$\left\{\begin{array}{c}y=R·\frac{2·w}{1+w^2}\hfill \\ z=-a+R·\frac{1-w^2}{1+w^2}\hfill \end{array}\right.,$$

(3)

for the system (2) with $R^2={\left[y\left(0\right)\right]}^2+{(z\left(0\right)+a)}^2$, Equation (2${}_3$) gives

$$\begin{array}{c}x=-\frac{2·\dot{w}}{1+w^2}.\hfill \end{array}$$

(4)

Equation (2${}_2$) identically satisfies, and Equation (2${}_1$) yields:

$$\begin{array}{c}{\displaystyle \ddot{w}\left(t\right)+a·\dot{w}\left(t\right)+a·R·w\left(t\right)=2·\frac{w\left(t\right)}{1+w^2\left(t\right)}·{\left(\dot{w}\left(t\right)\right)}^2,t>0,}\hfill \end{array}$$

(5)

which describes the a nonlinear anharmonic oscillator presented in details in the following section.

For the unknown function w the initial conditions are:

$$\begin{array}{c}{\displaystyle w\left(0\right)=\frac{y\left(0\right)}{R+z\left(0\right)+a},\dot{w}\left(0\right)=-\frac{x\left(0\right)}{2}·\left(1+w{\left(0\right)}^2\right).}\hfill \end{array}$$

(6)

For the nonlinear differential problem (5) and (6), the first-order approximate solutions, taking account of the following section and Equations (3) and (4), become:

$$\left\{\begin{array}{c}\overline{x}=-\frac{2·\dot{\overline{w}}}{1+{\overline{w}}^2}\hfill \\ \overline{y}=R·\frac{2·\overline{w}}{1+{\overline{w}}^2}\hfill \\ \overline{z}=-a+R·\frac{1-{\overline{w}}^2}{1+{\overline{w}}^2}\hfill \end{array}\right..$$

(7)

3. Application to the Nonlinear Anharmonic Oscillator

The nonlinear anharmonic oscillator given by Equation (5) in our approach could be described by the following equation:

$$\begin{array}{c}{\displaystyle \ddot{w}\left(t\right)(1+w^2\left(t\right))+a·\dot{w}\left(t\right)(1+w^2\left(t\right))+a·R·w\left(t\right)(1+w^2\left(t\right))-2·w\left(t\right)·{\left(\dot{w}\left(t\right)\right)}^2=0,t>0,}\hfill \end{array}$$

(8)

subject to the initial conditions:

$$\begin{array}{c}{\displaystyle w\left(0\right)=A_1,\dot{w}\left(0\right)=B_1,}\hfill \end{array}$$

(9)

with $A_1=\frac{y\left(0\right)}{R+z\left(0\right)+a}\ne 0$, $B_1=-\frac{x\left(0\right)}{2}·\left(1+w{\left(0\right)}^2\right)\ne 0$, $R=\sqrt{{\left[y\left(0\right)\right]}^2+{(z\left(0\right)+a)}^2}$, $a\in \mathbf{R}$ given real numbers.

Taking into account of the initial conditions given by Equation (9), the approximate solutions, denoted $\overline{w}\left(t\right)$, of Equation (8) are deducted for the unknown function $w\left(t\right)$.

By using OHAM method (Marinca and Herisanu [40,41,42], for more details) to Equations (8) and (9) an analytic approximate solution, denoted ${\overline{w}}_{OHAM}\left(t\right)$, is obtained using only one iteration.

For an embedding parameter $p\in [0,1]$, the first-order approximate solution $\overline{w}$ for nonlinear differential problem given by Equations (8) and (9) could be written as:

$$\begin{array}{c}{\displaystyle \overline{w}\left(t\right)=w_0\left(t\right)+p·w_1(t,C_j),}\hfill \end{array}$$

(10)

with $w_0\left(t\right)$ the initial approximation and $w_1(t,C_j)$ the first approximation depending on the variable t and the several unknown optimal parameters $C_1$, $C_2$, $C_3$, …, $C_{N_{max}}$, with $N_{max}$ an arbitrary integer number.

The homotopic relation is given by [41]:

$$\begin{array}{c}\mathcal{H}\left[\mathcal{L}\left(\overline{w}\left(t\right)\right),H(t,C_i,p),\mathcal{N}\left(\overline{w}\left(t\right)\right)\right]=\hfill \\ \\ =(1-p)·\mathcal{L}\left(\overline{w}\left(t\right)\right)-H(t,C_i,p)·\left[\mathcal{L}\left(\overline{w}(t,C_i)\right)+\mathcal{N}\left(\overline{w}\left(t\right)\right)\right]=0,\hfill \end{array}$$

(11)

where the nonlinear operator $\mathcal{N}\left(\overline{w}\left(t\right)\right)$ has the form

$$\begin{array}{c}{\displaystyle \mathcal{N}\left(\overline{w}\left(t\right)\right)=N_0\left(w_0\left(t\right)\right)+\sum _{m\ge 1}{N}_m(w_0,w_1,\dots ,w_m)·p^m}\hfill \end{array}$$

(12)

and $H(t,C_i,p)=p{h}_1(t,C_i)+p^2{h}_2(t,C_i)+p^3{h}_3(t,C_i)+\dots$ is a known auxiliary function. Taking into account Equations (10) and (12) the homotopic relation from Equation (11) becomes:

$$\begin{array}{c}\mathcal{H}\left[\mathcal{L}\left(\overline{w}\left(t\right)\right),H(t,C_i),\mathcal{N}\left(\overline{w}\left(t\right)\right)\right]=\mathcal{L}\left(w_0\left(t\right)\right)+p\left[\mathcal{L}\left(w_1(t,C_i)\right)-H(t,C_i)\mathcal{N}\left(w_0\left(t\right)\right)\right]=0,\hfill \end{array}$$

(13)

where $H(t,C_i)\ne 0$ is an auxiliary convergence-control function depending of the variable t and of the parameters $C_1$, $C_2$, …, $C_{N_{max}}$. The linear operator $\mathcal{L}\left(w\right)$ has the form:

$$\mathcal{L}\left(w\right)\left(t\right)=\ddot{w}+2K·\dot{w}+(K^2-{\omega }_0^2)w,$$

(14)

where $K,{\omega }_0>0$ are unknown parameters at this moment.

Therefore, the form of the nonlinear operator $\mathcal{N}\left(w\right)$ corresponding to the unknown function w is obtained from Equation (8) by:

$$\begin{array}{c}\mathcal{N}\left(w\right)\left(t\right)=-2K·\dot{w}-(K^2-{\omega }_0^2)w+\ddot{w}\left(t\right)·w^2\left(t\right)-2·w\left(t\right)·{\left(\dot{w}\left(t\right)\right)}^2+\hfill \\ +a·\dot{w}\left(t\right)·\left(1+w^2\left(t\right)\right)+a·R·\left(1+w^2\left(t\right)\right)·w\left(t\right).\hfill \end{array}$$

(15)

The deformations problems are obtained by identifying the coefficients $p^0$ and $p^1$, respectively.

Some possibilities to choose the auxiliary function $H(t,C_i)$ could be:

$${\displaystyle H(t,C_i)=\sum _{k=1}^{N_{max}}{a}_k^{\left(2\right)}·cos(2k+1){\omega }_{0}t+b_k^{\left(2\right)}·sin(2k+1){\omega }_{0}t,}$$

(16)

where $C_i\in \left\{a_k^{\left(2\right)}\right|k=\overline{1,N_{max}}\}\cup \left\{b_k^{\left(2\right)}\right|k=\overline{1,N_{max}}\}$, or

$$H(t,C_i)=C_{1}cos\left({\omega }_{0}t\right)+B_{1}sin\left({\omega }_{0}t\right),$$

or

$$H(t,C_i)=C_{1}cos\left({\omega }_{0}t\right)+B_{1}sin\left({\omega }_{0}t\right)+C_{2}cos\left(3{\omega }_{0}t\right)+B_{2}sin\left(3{\omega }_{0}t\right),$$

or

$$\begin{array}{c}H(t,C_i)=C_{1}cos\left({\omega }_{0}t\right)+B_{1}sin\left({\omega }_{0}t\right)+C_{2}cos\left(3{\omega }_{0}t\right)+B_{2}sin\left(3{\omega }_{0}t\right)+\hfill \\ +C_{3}cos\left(5{\omega }_{0}t\right)+B_{3}sin\left(5{\omega }_{0}t\right),\hfill \end{array}$$

and so on.

3.1. The Zeroth-Order Deformation Problem

For the initial approximation $w_0$, Equation (14) becomes:

$$\ddot{w}+2K·\dot{w}+(K^2-{\omega }_0^2)w=0,w\left(0\right)=A_1,\dot{w}\left(0\right)=B_1$$

(17)

with the solution

$$w_0\left(t\right)=w\left(0\right)e^{-Kt}cos\left({\omega }_{0}t\right)+\frac{\dot{w}\left(0\right)}{{\omega }_0}{e}^{-Kt}sin\left({\omega }_{0}t\right).$$

(18)

3.2. The First-Order Deformation Problem

The nonlinear operator Equation (15) for the initial approximation $w_0\left(t\right)$ given by Equation (18), using Equation (15) becomes:

$$\begin{array}{c}\mathcal{N}\left(w_0\right)\left(t\right)=M_1·e^{-Kt}·cos\left({\omega }_{0}t\right)+M_2·e^{-3Kt}·cos\left({\omega }_{0}t\right)+M_3·e^{-3Kt}·cos\left(3{\omega }_{0}t\right)+\hfill \\ +P_1·e^{-Kt}·sin\left({\omega }_{0}t\right)+P_2·e^{-3Kt}·sin\left({\omega }_{0}t\right)+P_3·e^{-3Kt}·sin\left(3{\omega }_{0}t\right),\hfill \end{array}$$

(19)

where

$$\begin{array}{cc}{M}_1=a{B}_1-aK{A}_1-2K{B}_1+A_1{K}^2+aR{A}_1+A_1{\omega }_0^2,\hfill & \\ M_2=\frac{1}{4}a{A}_1^2{B}_1-\frac{5}{4}{A}_1{B}_1^2-\frac{3}{4}aK{A}_1^3+\frac{1}{2}K{A}_1^2{B}_1-\frac{3}{4}{K}^2{A}_1^3+\frac{3}{4}aR{A}_1^3+\hfill & \\ +\frac{a{B}_1^3}{4{\omega }_0^2}-\frac{3aK{A}_1{B}_1^2}{4{\omega }_0^2}+\frac{K{B}_1^3}{2{\omega }_0^2}-\frac{3K^2{A}_1{B}_1^2}{4{\omega }_0^2}+\frac{3aR{A}_1{B}_1^2}{4{\omega }_0^2}-\frac{5}{4}{A}_1^3{\omega }_0^2,\hfill & \\ M_3=\frac{3}{4}a{A}_1^2{B}_1-\frac{3}{4}{A}_1{B}_1^2-\frac{1}{4}aK{A}_1^3+\frac{3}{2}K{A}_1^2{B}_1-\frac{1}{4}{K}^2{A}_1^3+\frac{1}{4}aR{A}_1^3-\hfill & \\ -\frac{a{B}_1^3}{4{\omega }_0^2}+\frac{3aK{A}_1{B}_1^2}{4{\omega }_0^2}-\frac{K{B}_1^3}{2{\omega }_0^2}+\frac{3K^2{A}_1{B}_1^2}{4{\omega }_0^2}-\frac{3aR{A}_1{B}_1^2}{4{\omega }_0^2}+\frac{1}{4}{\omega }_0^2{A}_1^3\hfill & \\ P_1=-\frac{aK{B}_1}{{\omega }_0}+\frac{K^2{B}_1}{{\omega }_0}+\frac{aR{B}_1}{{\omega }_0}-a{\omega }_0{A}_1+B_1{\omega }_0+2K{\omega }_0{A}_1,\hfill & \\ P_2=-\frac{3aK{B}_1^3}{4{\omega }_0^3}-\frac{3K^2{B}_1^3}{4{\omega }_0^3}+\frac{3aR{B}_1^3}{4{\omega }_0^3}-\frac{a{A}_1{B}_1^2}{4{\omega }_0}-\frac{5B_1^3}{4{\omega }_0}-\frac{3aK{A}_1^2{B}_1}{4{\omega }_0}-\frac{K{A}_1{B}_1^2}{2{\omega }_0}-\hfill & \\ -\frac{3K^2{A}_1^2{B}_1}{4{\omega }_0}+\frac{3aR{A}_1^2{B}_1}{4{\omega }_0}-\frac{1}{4}a{\omega }_0{A}_1^3-\frac{5}{4}{\omega }_0{A}_1^2{B}_1-\frac{1}{2}K{\omega }_0{A}_1^3\hfill & \\ P_3=\frac{aK{B}_1^3}{4{\omega }_0^3}+\frac{K^2{B}_1^3}{4{\omega }_0^3}-\frac{aR{B}_1^3}{4{\omega }_0^3}+\frac{3a{A}_1{B}_1^2}{4{\omega }_0}-\frac{B_1^3}{4{\omega }_0}-\frac{3aK{A}_1^2{B}_1}{4{\omega }_0}+\frac{3K{A}_1{B}_1^2}{2{\omega }_0}-\hfill & \\ -\frac{3K^2{A}_1^2{B}_1}{4{\omega }_0}+\frac{3aR{A}_1^2{B}_1}{4{\omega }_0}-\frac{1}{4}a{\omega }_0{A}_1^3+\frac{3}{4}{\omega }_0{A}_1^2{B}_1-\frac{1}{2}K{A}_1^3\hfill & \end{array}$$

It depends on the elementary functions $e^{-Kt}·cos\left({\omega }_{0}t\right)$, $e^{-3Kt}·cos\left({\omega }_{0}t\right)$, $e^{-3Kt}·cos\left(3{\omega }_{0}t\right)$, $e^{-Kt}·sin\left({\omega }_{0}t\right)$, $e^{-3Kt}·sin\left({\omega }_{0}t\right)$, $e^{-3Kt}·sin\left(3{\omega }_{0}t\right)$. For $p=1$, from Equation (13) the first-order deformation problem becomes:

$$\begin{array}{c}\mathcal{L}\left(w_1(t,C_i)\right)=H(t,C_i)\mathcal{N}\left(w_0\left(t\right)\right)\hfill \end{array}$$

(20)

By integration of the first approximation $w_1(t,C_i)$, from Equation (20) and considering for $H(t,C_i)$ the expression given by Equation (16) yields:

$$\begin{array}{c}{\displaystyle w_1(t,C_i)=\sum _{k=1}^{N_{max}}{C}_k·e^{-3Kt}·cos\left((2k+1){\omega }_{0}t\right)+B_k·e^{-3Kt}·sin\left((2k+1){\omega }_{0}t\right),}\hfill \end{array}$$

(21)

where $C_i,B_i$ are unknown parameters, with ${\displaystyle \sum _{k=1}^{N_{max}}{C}_k=0}$ and ${\displaystyle \sum _{k=1}^{N_{max}}(2k+1)·B_k=0}$, respectively.

3.3. The First-Order Analytical Approximate Solution $\overline{w}$

Taking into account of Equations (18) and (21) was obtained the first-order approximate solution given by Equation (10), when $p=1$ in the form:

$$\begin{array}{c}{\displaystyle \overline{w}\left(t\right)=w_0\left(t\right)+w_1(t,C_i)=w\left(0\right)e^{-Kt}cos\left({\omega }_{0}t\right)+\frac{\dot{w}\left(0\right)}{{\omega }_0}{e}^{-Kt}sin\left({\omega }_{0}t\right)+}\hfill \\ {\displaystyle +\sum _{k=1}^{N_{max}}{C}_k·e^{-3Kt}·cos\left((2k+1){\omega }_{0}t\right)+B_k·e^{-3Kt}·sin\left((2k+1){\omega }_{0}t\right),}\hfill \end{array}$$

(22)

where the unknown parameters $C_i$, $B_i$, $i=1,2,3,\cdots$, are optimally identified.

4. Numerical Simulation

In order to demonstrate the accuracy and validity of the OHAM technique, a comparison between our approximate solutions with corresponding numerical results obtained via the fourth-order Runge–Kutta method is highlighted.

We consider the initial value problem given by (2) with initial conditions $x\left(0\right)=1$, $y\left(0\right)=1$, $z\left(0\right)=1$, $a=0.25$ and $b=c=0$. The convergence-control parameters K, ${\omega }_0$, $\left\{B_i\right\}{|}_{i=\overline{1,N_{max}}}$, $\left\{C_i\right\}{|}_{i=\overline{1,N_{max}}}$ are optimally determined by means of the least-square method.

Table 1. Comparison between the analytical approximate solution ${\overline{w}}_{OHAM}$ of Equation (22) and the corresponding numerical values for different values of the index number $N_{max}$.

t	${\mathit{N}}_{\mathit{max}}=10$	${\mathit{N}}_{\mathit{max}}=15$	${\mathit{N}}_{\mathit{max}}=20$	${\mathit{N}}_{\mathit{max}}=25$
	$w_{numerical}$
0	0.35078105
1.7	−0.40330489
3.4	−0.63904518
5.1	−0.29126219
6.8	0.12552955
8.5	0.29643558
10.2	0.16028326
11.9	−0.05820289
13.6	−0.15188552
15.3	−0.08407041
17	0.03107464
	${\overline{w}}_{OHAM}$ obtained from Equation (22) and
	Equation (A1)	Equation (A2)	Equation (A3)	Equation (A4)
0	0.35078105	0.35078105	0.35078105	0.35078105
1.7	−0.40841782	−0.40330994	−0.40330493	−0.40330493
3.4	−0.64462314	−0.63904979	−0.63904521	−0.63904521
5.1	−0.29646094	−0.29126751	−0.29126222	−0.29126223
6.8	0.12060783	0.12552535	0.12552950	0.12552951
8.5	0.28751514	0.29643084	0.29643554	0.29643554
10.2	0.14938850	0.16027866	0.16028322	0.16028322
11.9	−0.06828322	−0.05819707	−0.05820297	−0.05820293
13.6	−0.16709523	−0.15188883	−0.15188558	−0.15188556
15.3	−0.09821857	−0.08408416	−0.08407045	−0.08407044
17	0.01573297	0.03106425	0.03107455	0.03107461
	Relative errors $=|w_{numerical}-{\overline{w}}_{OHAM}|$
0	5.689 × ${10}^{-14}$	5.551 × ${10}^{-17}$	5.689 × ${10}^{-14}$	2.841 × ${10}^{-13}$
1.7	0.00511292	5.049 × ${10}^{-6}$	3.829 × ${10}^{-8}$	3.423 × ${10}^{-8}$
3.4	0.00557796	4.616 × ${10}^{-6}$	3.491 × ${10}^{-8}$	3.438 × ${10}^{-8}$
5.1	0.00519875	5.321 × ${10}^{-6}$	3.786 × ${10}^{-8}$	4.167 × ${10}^{-8}$
6.8	0.00492171	4.197 × ${10}^{-6}$	4.579 × ${10}^{-8}$	4.043 × ${10}^{-8}$
8.5	0.00892043	4.732 × ${10}^{-6}$	3.524 × ${10}^{-8}$	3.975 × ${10}^{-8}$
10.2	0.01089475	4.599 × ${10}^{-6}$	3.694 × ${10}^{-8}$	4.065 × ${10}^{-8}$
11.9	0.01008032	5.821 × ${10}^{-6}$	7.353 × ${10}^{-8}$	3.995 × ${10}^{-8}$
13.6	0.01520970	3.307 × ${10}^{-6}$	5.836 × ${10}^{-8}$	3.493 × ${10}^{-8}$
15.3	0.01414815	1.374 × ${10}^{-5}$	3.568 × ${10}^{-8}$	2.826 × ${10}^{-8}$
17	0.01534166	1.039 × ${10}^{-5}$	8.433 × ${10}^{-8}$	3.114 × ${10}^{-8}$

,

Table 2. Comparison between the analytic closed-form approximate solution $\overline{x}$ given by Equation (7)${}_1$ and corresponding numerical solution; the relative errors ${\epsilon }_x=|x_{numerical}-{\overline{x}}_{OHAM}|$.

t	${\mathit{x}}_{\mathit{numerical}}$	${\overline{\mathit{x}}}_{\mathit{OHAM}}$	Relative Errors ${\mathsf{\epsilon }}_{\mathit{x}}$
0	1	0.99999998	1.781 × ${10}^{-8}$
1.7	0.57062883	0.57062880	3.337 × ${10}^{-8}$
3.4	−0.10808510	−0.1080850	4.599 × ${10}^{-8}$
5.1	−0.49392956	−0.4939297	1.456 × ${10}^{-7}$
6.8	−0.38333158	−0.3833315	8.265 × ${10}^{-8}$
8.5	0.00446477	0.0044648	3.685 × ${10}^{-8}$
10.2	0.25485018	0.2548503	1.881 × ${10}^{-7}$
11.9	0.20964292	0.2096428	2.959 × ${10}^{-8}$
13.6	0.00265313	0.0026531	2.432 × ${10}^{-8}$
15.3	−0.13494391	−0.1349440	9.170 × ${10}^{-8}$
17	−0.11083968	−0.1108396	2.099 × ${10}^{-8}$

,

Table 3. Comparison between the analytic closed-form approximate solution $\overline{y}$ given by Equation (7)${}_2$ and corresponding numerical solution; the relative errors ${\epsilon }_y=|y_{\mathit{numerical}}-{\overline{y}}_{\mathit{OHAM}}|$.

t	${\mathit{y}}_{\mathit{numerical}}$	${\overline{\mathit{y}}}_{\mathit{OHAM}}$	Relative Errors ${\mathsf{\epsilon }}_{\mathit{y}}$
0	1	1.0000000000006326	6.326 × ${10}^{-13}$
1.7	−1.11056659	−1.11056670	1.180 × ${10}^{-7}$
3.4	−1.45269382	−1.45269367	1.479 × ${10}^{-7}$
5.1	−0.85957282	−0.85957335	5.275 × ${10}^{-7}$
6.8	0.39565568	0.39565592	2.411 × ${10}^{-7}$
8.5	0.87239595	0.87239585	9.511 × ${10}^{-8}$
10.2	0.50030350	0.50030353	3.028 × ${10}^{-8}$
11.9	−0.18571072	−0.18571121	4.917 × ${10}^{-7}$
13.6	−0.47530599	−0.47530611	1.210 × ${10}^{-7}$
15.3	−0.26726775	−0.26726775	8.200 × ${10}^{-9}$
17	0.09939116	0.09939132	1.591 × ${10}^{-7}$

and

Table 4. Comparison between the analytic closed-form approximate solution $\overline{z}$ given by Equation (7)${}_3$ and corresponding numerical solution; the relative errors ${\epsilon }_z=|z_{numerical}-{\overline{z}}_{OHAM}|$.

t	${\mathit{z}}_{\mathit{numerical}}$	${\overline{\mathit{z}}}_{\mathit{OHAM}}$	Relative Errors ${\mathsf{\epsilon }}_{\mathit{z}}$
0	1	0.999999999999494	5.060 × ${10}^{-13}$
1.7	0.90288415	0.90288402	1.293 × ${10}^{-7}$
3.4	0.42244420	0.42244411	8.863 × ${10}^{-8}$
5.1	1.10041934	1.10041980	4.580 × ${10}^{-7}$
6.8	1.30111457	1.30111456	1.648 × ${10}^{-8}$
8.5	1.09217208	1.09217192	1.641 × ${10}^{-7}$
10.2	1.27059023	1.27059079	5.585 × ${10}^{-7}$
11.9	1.33997205	1.33997212	6.118 × ${10}^{-8}$
13.6	1.27858895	1.27858892	3.430 × ${10}^{-8}$
15.3	1.32831159	1.32831174	1.498 × ${10}^{-7}$
17	1.34769247	1.34769251	3.643 × ${10}^{-8}$

emphasizes the accuracy of the OHAM technique by comparing the approximate analytic solutions $\overline{w}$, $\overline{x}$, $\overline{y}$, and $\overline{z}$, presented above with the corresponding numerical integration values (via the fourth-order Runge–Kutta method). These comparisons show the effectiveness, reliability, applicability, and efficiency of the OHAM.

All the convergence-control parameters, corresponding to the intervals $[0,20]$ are presented in Appendix A.

Figure 1, Figure 2 and Figure 3 present the comparisons between the analytical approximate solutions given by OHAM and numerical results provided by Runge–Kutta fourth steps integrator, for initial conditions $x\left(0\right)=1$, $y\left(0\right)=1$, $z\left(0\right)=1$, $a=0.25$ and for $N_{max}=25$. From these figures it follows that the behavior of the analytical approximate solutions and Runge–Kutta fourth steps integrator’s results are quite the same.

The corresponding relative errors are presented in detail in Appendix A.

Below we highlight the advantages of the OHAM method by comparison with an iterative method developed in [46]. By integration of the system (2) over the interval $[0,t]$, the following relations are obtained:

$$\begin{array}{c}{\displaystyle x\left(t\right)=x\left(0\right)+\underset{0}{\overset{t}{\int }}a(y\left(s\right)-x\left(s\right))ds}\hfill \\ {\displaystyle y\left(t\right)=y\left(0\right)+\underset{0}{\overset{t}{\int }}(-ax\left(s\right)-x\left(s\right)z\left(s\right))ds}\hfill \\ {\displaystyle z\left(t\right)=z\left(0\right)+\underset{0}{\overset{t}{\int }}x\left(s\right)y\left(s\right)ds}\hfill \end{array}.$$

(23)

The iterative algorithm is the following:

$$\begin{array}{c}{\displaystyle x_0\left(t\right)=x\left(0\right),x_1\left(t\right)=N_1(x_0,y_0,z_0)=\underset{0}{\overset{t}{\int }}a(y_0\left(s\right)-x_0\left(s\right))ds,}\hfill \\ {\displaystyle y_0\left(t\right)=y\left(0\right),y_1\left(t\right)=N_2(x_0,y_0,z_0)=\underset{0}{\overset{t}{\int }}(-a{x}_0\left(s\right)-x_0\left(s\right)z_0\left(s\right))ds,}\hfill \\ {\displaystyle z_0\left(t\right)=z\left(0\right),z_1\left(t\right)=N_3(x_0,y_0,z_0)=\underset{0}{\overset{t}{\int }}{x}_0\left(s\right)y_0\left(s\right)ds,}\hfill \\ \\ x_m\left(t\right)=N_1(x_0+x_1+\cdots +x_{m-1},y_0+y_1+\cdots +y_{m-1},z_0+z_1+\cdots +z_{m-1})-\hfill \\ -N_1(x_0+x_1+\cdots +x_{m-2},y_0+y_1+\cdots +y_{m-2},z_0+z_1+\cdots +z_{m-2}),\hfill \\ y_m\left(t\right)=N_2(x_0+x_1+\cdots +x_{m-1},y_0+y_1+\cdots +y_{m-1},z_0+z_1+\cdots +z_{m-1})-\hfill \\ -N_2(x_0+x_1+\cdots +x_{m-2},y_0+y_1+\cdots +y_{m-2},z_0+z_1+\cdots +z_{m-2}),\hfill \\ z_m\left(t\right)=N_3(x_0+x_1+\cdots +x_{m-1},y_0+y_1+\cdots +y_{m-1},z_0+z_1+\cdots +z_{m-1})-\hfill \\ -N_3(x_0+x_1+\cdots +x_{m-2},y_0+y_1+\cdots +y_{m-2},z_0+z_1+\cdots +z_{m-2}),m\ge 2.\hfill \end{array}$$

(24)

With the iterative method, the solutions of system (2) have the form:

$${\displaystyle x_{iter}\left(t\right)=\sum _{m=0}^{\infty }{x}_m\left(t\right),y_{iter}\left(t\right)=\sum _{m=0}^{\infty }{y}_m\left(t\right),z_{iter}\left(t\right)=\sum _{m=0}^{\infty }{z}_m\left(t\right).}$$

For the case using seven iterations, with the initial conditions $x\left(0\right)=1$, $y\left(0\right)=1$, $z\left(0\right)=1$ and the constants $a=0.25$, $b=c=0$, taking into account of the algorithm (24), the iterative solutions become:

$$\begin{array}{c}{\displaystyle x_{iter}\left(t\right)=\sum _{m=0}^7{x}_m\left(t\right)=1-0.15625t^2-0.02864583t^3+0.01888020t^4+0.00419108t^5-}\hfill \\ -0.00186326t^6-0.00066741t^7+0.00016160t^8+0.00008497t^9-9.153078·{10}^{-6}{t}^{10}-\hfill \\ -8.689161·{10}^{-6}{t}^{11}+6.523628·{10}^{-8}{t}^{12}+7.013146·{10}^{-7}{t}^{13}+5.223197·{10}^{-8}{t}^{14}-4.322356·{10}^{-8}{t}^{15},\hfill \\ \\ {\displaystyle y_{iter}\left(t\right)=\sum _{m=0}^7{y}_m\left(t\right)=1-1.25t-0.5t^2+0.2734375t^3+0.10270182t^4-0.04052734t^5-}\hfill \\ -0.0205508t^6+0.00508378t^7+0.0034984t^8-0.00048158t^9-0.0005136t^{10}+0.00002626t^{11}+\hfill \\ +0.0000650t^{12}+1.353350·{10}^{-6}{t}^{13}-7.134374·{10}^{-6}{t}^{14}-6.095126·{10}^{-7}{t}^{15}+6.773822·{10}^{-7}{t}^{16}+\hfill \\ +1.007179·{10}^{-7}{t}^{17}-5.542185·{10}^{-8}{t}^{18}-1.201667·{10}^{-8}{t}^{19},\hfill \\ \\ {\displaystyle z_{iter}\left(t\right)=\sum _{m=0}^7{z}_m\left(t\right)=1+t-0.625t^2-0.21874999t^3+0.11002604t^4+0.04710286t^5-}\hfill \\ -0.01472303t^6-0.00871044t^7+0.00159948t^8+0.00138481t^9-0.00011177t^{10}-0.00018913t^{11}-\hfill \\ -7.689201·{10}^{-7}{t}^{12}+0.00002208t^{13}+1.751434·{10}^{-6}{t}^{14}-2.202022·{10}^{-6}{t}^{15}-3.375254·{10}^{-7}{t}^{16}+\hfill \\ +1.858551·{10}^{-7}{t}^{17}+4.383605·{10}^{-8}{t}^{18}-1.298323·{10}^{-8}{t}^{19}-4.480401·{10}^{-9}{t}^{20}.\hfill \end{array}$$

(25)

A comparison between approximate closed-form solutions ${\overline{x}}_{OHAM}$, ${\overline{y}}_{OHAM}$, ${\overline{z}}_{OHAM}$ and the corresponding iterative solutions $x_{iter}$, $y_{iter}$, $z_{iter}$ given in (25) is shown in the Appendix A both graphically in Figure 4 and tabularly in

Table 5. Comparison between the approximate closed-form solution $\overline{x}$ given by Equation (7)${}_1$, corresponding numerical solution and the iterative solution $x_{iter}$ given by Equation (25).

t	${\mathit{x}}_{\mathit{numerical}}$	${\overline{\mathit{x}}}_{\mathit{OHAM}}$	${\mathit{x}}_{\mathit{iter}}$
0	1	0.99999998	1
1/2	0.95863428	0.95863428	0.95863421
1	0.83589390	0.83589396	0.83587428
3/2	0.65318706	0.65318710	0.65280614
2	0.44281996	0.44281983	0.44040998
5/2	0.23134451	0.23134437	0.22345920
3	0.03418218	0.03418216	0.01748463
7/2	−0.14118335	−0.14118327	−0.16734859
4	−0.28999110	−0.28999104	−0.32509324
9/2	−0.40661666	−0.40661673	−0.48802039
5	−0.48384683	−0.48384691	0.22295502

, respectively.

This comparison shows the precision and efficiency of the OHAM method (using just one iteration) against to the iterative method described in [46] (using seven iterations).

5. Conclusions

For a special case of the Chen system (just one analytic first integral) is shown that this system could be reduced to a nonlinear anharmonic oscillator. An analytic approximate solution for the anharmonic oscillator problem was been obtained by means of the Optimal Homotopy Asymptotic Method. The numerical outcomes contribute to a better knowledge of the accuracy and validity of the OHAM technique. The flexibility of this method results from the comparison with the corresponding iterative procedure. The results of our present study are useful in understanding of the behavior for the dynamical systems as complete synchronization or optimization of nonlinear system performance.