# Construction of homoclinic and heteroclinic trajectories in mechanical systems with several equilibrium positions

Yuri Mikhlin, Gayana Manucharyan

► To cite this version:

Yuri Mikhlin, Gayana Manucharyan. Construction of homoclinic and heteroclinic trajectories in mechanical systems with several equilibrium positions. Chaos, Solitons & Fractals, 2003, 16 (2), pp. 299-309. 10.1016/S0960-0779(02)00404-6. hal-01348409

HAL Id: hal-01348409

https://hal.science/hal-01348409

Submitted on 23 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Construction of homoclinic and heteroclinic trajectories in mechanical systems with several equilibrium positions

Yu.V. Mikhlin *, G.V. Manucharyan

Department of Applied Mathematics, Kharkov Polytechnical University, 21 Frunze str., Kharkov 61002, Ukraine

## Abstract

A new approach for a construction of homo and heteroclinic trajectories of some principal non-linear dynamical systems is utilized here, namely the non-linear Schrodinger equation, non-autonomous Duffing equation and the equation of a parametrically excited damped pendulum are considered. Pade' and quasi-Pade' approximants and a convergence condition used earlier in the theory of non-linear normal vibration modes made possible to solve a boundary-value problems formulated for the orbits and to determine initial amplitude values of the trajectories with admissible precision. The approach proposed here is more exact than the generally accepted one because it is not necessary to use here separatrix trajectories of the corresponding autonomous equations.

## 1. Introduction

Homo and heteroclinic trajectories (HT) have been extensively studied in the literature [1,2]. A formation of HT is a criterion of a chaotic behavior beginning in dynamical systems.

Note that it is necessary to determine some important parameters for a construction of HT. For example, we have to know four parameters including two initial values to construct HT of the non-autonomous Duffing equation (we discuss the problem later, in Section 4). In most cases authors of last and recent publications on the HT construction make use of the well-known Melnikov condition of the trajectory formation [3–6] which gives us a single equation for a determination of all unknown parameters. As a result in the Melnikov condition, separatrix trajectories of the corresponding autonomous equations that is HT of zero approximation are utilized. A problem of effective analytic approximation of HT of non-autonomous system is difficult and it is not solved up to now.

Here a new approach for the HT construction in the non-linear dynamical systems with phase space of dimensions equal to two is proposed. Pade' approximants (PA) and quasi-Pade' approximants (QPA) are used for a representation both the HT in the dynamical system phase space and the corresponding time solution. (Note that QPA which contain powers of some parameter and exponential functions were considered in Ref. [7].) Convergence condition used earlier in the theory of non-linear normal vibration modes [8–10] as well conditions at infinity made possible to solve the boundary-value problem formulated for the HT and evaluate initial amplitude values with admissible precision. It is important that the HT criterion of the chaos beginning proposed here is more exact than the generally accepted Melnikov criterion, because it is not necessary to use here separatrix trajectories of the corresponding autonomous equations. Note that PA were successfully used earlier for the HT analytic approximation in non-linear Schrodinger equation and Lorenz system [11–14].
This work is structured as follow. First, the convergence condition is discussed in Section 2. Realization of the approach proposed here was made for the localized solution of the non-linear Schrodinger equation. Corresponding results are presented in Section 3. A boundary-value problem for homoclinic orbits and a construction of the trajectories for the non-autonomous Duffing equation is discussed in Section 4. A boundary-value problem and a construction of heteroclinic orbits was presented for the equation of a parametrically exited damped pendulum in Section 5. Comparisons between analytical estimates and numerical calculations are also presented with a discussion of the main features of the proposed approach.

## 2. Convergence condition

Let us assume that there are local expansions of solution obtained at small and large values of a parameter $c$ (for example, the parameter is an amplitude value or initial energy of the system). For small values of $c$ the local expansion can be determined as a power series in $c$, while for large values of $c$ it can be determined as a power series in $c^{-1}$:

$$y^{(0)} = \sum_{j=0}^{\infty} \alpha_j c^j, \quad y^{(\infty)} = \sum_{j=0}^{\infty} \beta_j c^{-j}. \qquad (1)$$

In order to join local expansions (1), fractional rational diagonal two-point PA [15] can be used. Let us consider

$$PA_s = \frac{\sum_{j=0}^{s} a_j c^j}{\sum_{j=0}^{s} b_j c^j} = \frac{\sum_{j=0}^{s} a_j c^{j-s}}{\sum_{j=0}^{s} b_j c^{j-s}} \quad (s = 1, 2, 3, \dots). \qquad (2)$$

By comparison expressions (1) and (2) and retaining only the terms with an order of $c^r(-s \le r \le s)$ we will obtain a system of $2(s+1)$ linear algebraic equations for the determination of coefficients $a_j, b_j$. Since generally the determinant of the system $\Delta_s$, is not equal to zero, the system has a single trivial exact solution. But we need in PA corresponding to the retaining terms in Eq. (1) having non-zero coefficients $a_j, b_j$. Without loss of generality it can be assumed that $b_0 = 1$. Now, the system of algebraic equations for determination $a_j, b_j$ becomes overdetermined. All of the unknown coefficients can be determined from $(2s+1)$ equations while the “residual” of this approximate solution can be obtained by substitution of all the coefficients into the remaining equation. Obviously, the residual (or “error”) is determined by the value of $\Delta_s$ since non-zero solutions and consequently exact PA will be obtained in the given approximation by $c$ only in the case when $\Delta_s = 0$. Hence the following is a necessary condition for convergence of the succession of $PA_s$ (2) at $s \to \infty$ to fractional rational function $P_\infty$ [8–10]. Namely,

$$\lim_{s \to \infty} \Delta_s = 0. \qquad (3)$$

It is possible to generalize the necessary condition for convergence (3) to QPA which contain powers of some unknown parameter and exponential functions. Besides, it is possible to utilize the condition (3) for obtaining some unknown parameters which is contained in local expansions.

**Example.** Consider a function $f(x) = \sqrt{(1+x)/(1+4x)}$. Local expansion of this function at small and large values of $x$ are the following:

$$f(x) \approx 1 - \frac{3}{2}x + \frac{39}{8}x^2 - \frac{267}{16}x^3 + \frac{7563}{128}x^4 + \dots \quad (x \to 0), \qquad (4)$$

$$f(x) \approx \frac{1}{2} + \frac{3}{16} \frac{1}{x} - \frac{21}{256} \frac{1}{x^2} + \dots \quad (x \to \infty). \qquad (5)$$

Let us introduce some unknown parameter $B$ in the second expansion (5):

$$f(x) \approx B + \frac{3}{16} \frac{1}{x} - \frac{21}{256} \frac{1}{x^2} + \dots \quad (x \to \infty). \qquad (6)$$

It is possible to obtain parameter $B$ constructing PA and using convergence condition (3).

So, let

$$f(x) \approx PA_1 = \frac{\alpha_0 + \alpha_1 x}{1 + \beta_1 x}.$$
Comparing it with expansions (4) and (6) one obtain four equations:

$$\alpha_0 - 1 = 0, \quad \beta_j - \frac{3}{2} - \alpha_1 = 0,$$

$$B\beta_1 - \alpha_1 = 0, \quad -\alpha_0 + B + \frac{3}{16}\beta_1 = 0.$$

Simple calculations show that an ambiguity of this approximation will be reduced if $B = 0.46967$.

If we will consider $PA_2 = (\alpha_0 + \alpha_1 x + \alpha_2 x^2)/(1 + \beta_1 x + \beta_2 x^2)$ then we obtain six equations:

$$\alpha_0 - 1 = 0, \quad \beta_1 - \frac{3}{2} - \alpha_1 = 0, \quad -\alpha_2 + \beta_2 - \frac{3}{2}\beta_1 + \frac{39}{8} = 0,$$

$$-\frac{3}{2}\beta_2 + \frac{39}{8}\beta_1 - \frac{267}{16} = 0, \quad B\beta_2 - \alpha_2 = 0, \quad -\alpha_1 + \frac{3}{16}\beta_2 + B\beta_1 = 0.$$

The value $B = 0.498595$ can be obtained from a condition $\Delta_2 = 0$ corresponding to Eq. (3) in this approximation. It is clear that sequence $B^{(i)}$ converges to the correct value $B = 0.5$.

### 3. Non-linear Schrodinger equation

In this section a realization of the approach proposed in this paper will be made for a localized solution of non-linear Schrodinger equation. The boundary-value problem considered here can be obtained in a problem of localized axially symmetric solutions of two-dimensional non-linear Schrodinger equation:

$$y''(x) + \frac{1}{x}y'(x) - y(x) + y^3(x) = 0, \quad y'(0) = y'(\infty) = 0, \quad y(\infty) = 0. \qquad (7)$$

A spectrum of the problem is discrete (a corresponding theorem of existence is contained in [16]). But all the solutions excepting one of them are unstable. Since the sought solutions are expected to be analytic functions of $x$, they can be expressed in Taylor series in vicinity of zero:

$$y(x) = a_0 + a_2x^2 + a_4x^4 + a_6x^6 + a_8x^8 + \dots,$$

where $a_0$ is an arbitrary constant, $a_2 = -(1/4)a_0^3 + (1/4)a_0$, $a_4 = -(1/16)a_0^3 + (1/64)a_0 + (3/64)a_0^5$, $a_6 = (13/768)a_0^5 - (19/2304)a_0^3 - (7/768)a_0^7 + (1/2304)a_0$, ...

To determine an initial value of the solution under consideration $a_0$, we will use the following procedure. We multiply the Eq. (7) by $y'(x)$ and integrate along the trajectory of a solution under consideration, substitute the expansion of $y(x)$ and integrate within limits from $x = 0$ to $x = \infty$: $0 = \int_0^\infty [y'(x) - y(x) + y^3(x)]y' dx$. It means that an energy along the integral trajectory saves on the average.

Substituting the expansion of $y(x)$ in the integral, one can obtain after integration the following:

$$0 = \int_{0}^{\infty} \left[ \frac{1}{x'} y' - y + y^3 \right] y' dx = \frac{1}{2} a_0^2 - \frac{1}{4} a_0^4 + (Ax^2 + Bx^4 + Cx^6 + \dots)|_{0}^{\infty}, \qquad (8)$$

$$\text{where } A = 2a_2^2, B = 4a_2a_4, C = (8a_4^2 + 12a_6a_2)/3, D = 4a_8a_2 + 6a_4a_6, \dots$$

Then, assuming that the obtained local expansion $Ax^2 + Bx^4 + Cx^6 + \dots$ represents an analytical function, we rebuild the expansion to the diagonal $PA_m$. This is an analytical continuation ad infinitum of the corresponding local expansions. A succession of the $PA_m$ is the following:

$$PA_4 = \frac{\alpha_2 x^2 + \alpha_4 x^4}{1 + \beta_2 x^2 + \beta_4 x^4}, \quad PA_6 = \frac{\alpha_2 x^2 + \alpha_4 x^4 + \alpha_6 x^6}{1 + \beta_2 x^2 + \beta_4 x^4 + \beta_6 x^6}, \quad PA_8 = \frac{\alpha_2 x^2 + \alpha_4 x^4 + \alpha_6 x^6 + \alpha_8 x^8}{1 + \beta_2 x^2 + \beta_4 x^4 + \beta_6 x^6 + \beta_8 x^8}, \dots$$

where unknown coefficients $\alpha_i, \beta_i$ are computed in terms of $a_0$ by comparing $PA_m$ and the local expansion from Eq. (8). Using conditions at infinity we obtain from (8) a succession of the following algebraic equations to determine $y(0) = a_0$:

$$\frac{\alpha_4}{\beta_4} + \frac{1}{2} a_0^2 - \frac{1}{4} a_0^4 = 0, \quad \frac{\alpha_6}{\beta_6} + \frac{1}{2} a_0^2 - \frac{1}{4} a_0^4 = 0, \quad \frac{\alpha_8}{\beta_8} + \frac{1}{2} a_0^2 - \frac{1}{4} a_0^4 = 0, \quad \text{etc.}$$

By solving these equations next values of $a_0$ were obtained:

$$PA_4 : a_0 = \pm 1.6563; \quad PA_6 : a_0 = \pm 2.1465; \quad PA_8 : a_0 = \pm 2.2380, \quad \text{etc.} \qquad (9)$$
The initial value corresponding to the decaying solution was numerically estimated by means of selection as $a_0^{\text{num}} \cong \pm 2.206$ [11]. It is clear that the succession (9) converges to the exact value.

### 4. Non-autonomous Duffing equation

Let us consider the well-known non-autonomous Duffing equation which can be obtained, for example, by discretization of dynamical models of non-linear rods, plates or shells:

$$y'' + \delta y' - \alpha y + \beta y^3 = f \cos \omega t, \quad (10)$$

where $y(t)$ is an unknown function; $\delta$ is a coefficient determining a friction; all coefficient are positive. Eq. (10) makes sense for a description of post-buckling dynamics of non-linear elastic systems.

Let $\alpha = 1, \beta = 1$. For a construction of HT we need information about initial point $(a_0, a_1)$, corresponding to time moment $t = 0$, and about some special relation of the system parameters, namely $\omega, f$ and $\delta$. Also it is necessary to know a shift of a saddle point of the autonomous Duffing equation $b_0$. Thus we should obtain a system of equations for determination of $a_0, a_1, b_0$, and the pointed out relation of parameters $\omega, f$ and $\delta$.

Let us do some assumptions for solving of this problem. We will use PA and QPA representations for the HT and for the corresponding time-solution. The convergence conditions (3) as well conditions at infinity give us a system of algebraic equations to obtain unknown parameters.

The first one is the condition at infinity. We assume that HT intersects an unstable equilibrium position, that is $(y, y') \underset{t \to \pm\infty}{\to} (b_0, 0)$.

Besides let us assume that $y(t)$ is an analytical function on the HT, so we can consider local expansion of $y(t)$ in Taylor series in vicinity of zero.

Using these assumptions, PA and QPA and convergence condition we will obtain the unknown $a_0, a_1, b_0$ and the relation between parameters $f, \delta$ corresponding to HT while a frequency $\omega$ fixed.

Let us consider a local expansion of $y(t)$ in Taylor series in vicinity of zero:

$$y = a_0 + a_1 t + a_2 t^2 + a_3 t^3 + a_4 t^4 + a_5 t^5 + a_6 t^6 + \dots, \quad (11)$$

where $a_0, a_1$ are the arbitrary constants, $a_2 = -(a_0^3 - a_0 - f + \delta a_1)/2$, $a_3 = (a_1 - 3a_0^2 a_1 - 2\delta a_2)/6$, $a_4 = (2a_2 - fw^2 - 6\delta a_3 - 6a_1^2 a_0 - 6a_2 a_0^2)/24, \dots$

Let us use now the procedure which was considered for the non-linear Schrodinger equation. If we multiply Eq. (10) by $y'(t)$ and integrate obtained equation within limits from $t = 0$ to arbitrary $t$ along any trajectory, we can obtain the following relation:

$$\int_0^t (y'' - y + y^3 + \delta y' - f \cos \omega t) y' dt = -\frac{y^2(t)}{2} + \frac{a_0^2}{2} + \frac{y^4(t)}{4} - \frac{a_0^4}{4} + \frac{\dot{y}^2(t)}{2} - \frac{a_1^2}{2} + \int_0^t (\delta y' - f \cos \omega t) y' dt = 0.$$

The last relation represents phase trajectories of the non-autonomous Duffing equation. Consider now HT which reaches an unstable equilibrium position $(b_0, 0)$ at infinity. If we integrate within limits from $t = 0$ to $t = \pm\infty$ we can obtain the following:

$$\int_0^{\pm\infty} (y'' - y + y^3 + \delta y' - f \cos \omega t) y' dt = -\frac{b_0^2}{2} + \frac{a_0^2}{2} + \frac{b_0^4}{4} - \frac{a_0^4}{4} - \frac{a_1^2}{2} + \int_0^{\pm\infty} (\delta y' - f \cos \omega t) y' dt = 0. \quad (12)$$

It is interesting that an integration along the closed HT within limits from $t = -\infty$ to $t = +\infty$ gives us the well-known Melnikov condition of the form:

$$\int_{-\infty}^{+\infty} (\delta y' - f \cos \omega t) y' dt = 0.$$

The equation can be considered as a condition of conservation of the non-autonomous dynamical system energy along the closed HT on the average. Note that here the Melnikov condition is not directly used.

Let us consider $\int (\delta y' - f \cos \omega t) y' dt$. After substitution of the local expansion of $y(t)$ and integration one obtains:

$$\int (\delta y' - f \cos \omega t) y' dt = At + Bt^2 + Ct^3 + Dt^4 + Et^5 + \dots, \quad (13)$$
where

$$A = (\delta a_1 - f)a_1, \quad B = (2(\delta a_1 - f)a_2 + 2\delta a_2 a_1)/2,$$

$$C = (3(\delta a_1 - f)a_3 + 4\delta a_2^2 + (fw^2/2 + 3\delta a_3)a_1)/3,$$

$$D = (4(\delta a_1 - f)a_4 + 4\delta a_4 a_1 + 6\delta a_2 a_3 + 2(fw^2/2 + 3\delta a_3)a_2)/4.$$

Now let us reorganize the expansion (13) to the PA:

$$PA_3^p = \frac{\alpha_1 t + \alpha_2 t^2 + \alpha_3 t^3}{1 + \beta_1 t + \beta_2 t^2 + \beta_3 t^3}. \qquad (14)$$

Comparing $PA_3^p$ in the form (14) and the local expansion (13) one has: $\alpha_1 = A$,

$$\alpha_2 = \frac{-ADC^2 - DA^2E + 2AD^2B + FA^2C - AFB^2 + B^3E - 2B^2DC + BC^3}{AEC - AD^2 - B^2E + 2BDC - C^3},$$

$$\alpha_3 = \frac{-A^2E^2 + 2AEC^2 + 2AEBD - 2ACD^2 - 2ACBF + DFA^2}{AEC - AD^2 - B^2E + 2BDC - C^3} + \frac{-2CB^2E + 3BDC^2 - D^2B^2 + FB^3 - C^4}{AEC - AD^2 - B^2E + 2BDC - C^3},$$

$$\beta_2 = \frac{AE^2 - EC^2 - EBD + CD^2 + CBF - DFA}{AEC - AD^2 - B^2E + 2BDC - C^3},$$

$$\beta_3 = \frac{-BDF + BE^2 - 2DCE + C^2F + D^3}{AEC - AD^2 - B^2E + 2BDC - C^3}.$$

Passing to the limit at infinity in the relation (13) and taking into account (14) we can obtain the following:

$$-\frac{b_0^2}{2} + \frac{a_0^2}{2} + \frac{b_0^4}{4} - \frac{a_0^4}{4} - \frac{a_1^2}{2} + \frac{\alpha_3}{\beta_3} = 0. \qquad (15)$$

The additional equation is obtained using a convergence condition (3) for $PA_3^p$ (14):

$$-2D^2BF + 2DBE^2 - 3D^2CE + 2DC^2F + D^4 - AE^3 + E^2C^2 - 2ECBF + 2EDFA - F^2AC + F^2B^2 + GAEC - GAD^2 - GB^2E + 2GBDC - GC^3 = 0. \qquad (16)$$

Now let us reorganize the local expansion (11) to the QPA which is similar to the separatrix solution representation in autonomous Duffing equation to obtain an analytical continuation of lie local expansion (16) ad infinitum. It is important that a correctness of the representation which contains exponential functions, is verified by computer calculations using the Runge-Kutta procedure (see, for example, Figs. 3 and 4 which are presented later).

$$y = a_0 + a_1 t + a_2 t^2 + a_3 t^3 + a_4 t^4 + \dots \approx e^{-t} \frac{\alpha_0 + \alpha_1 e^t + \alpha_2 e^{2t} + \alpha_3 e^{3t}}{1 + \beta_2 e^{2t}}. \qquad (17)$$

It follows from (17) that

$$b_0 = \frac{\alpha_3}{\beta_2}, \qquad (18)$$

where coefficients $\alpha_i, \beta_j$ are determined from the Eq. (17) by means of expansion of exponential functions in powers of $t$ and comparison of power coefficients.

Finally we can obtain a convergence condition in the form (3) for QPA in the form (17):

$$24a_5a_3 + 2a_5a_1 + 12a_5a_2 - 24a_4^2 - 4a_1a_4 - \frac{4}{15}a_1a_2 - \frac{7}{10}a_1a_3 - \frac{1}{10}a_1^2 - 8a_2a_4 + \frac{5}{6}a_2^2 + 4a_2a_3 + 6a_3^2 - 12a_4a_3 = 0. \qquad (19)$$

Eqs. (15), (16), (18) and (19) form a system of four non-linear algebraic equations for calculation $a_0, a_1, b_0$ and $f = f(\delta)$ while a frequency $\omega$ is fixed.

The algebraic system is solved by using Newton method while $\omega = 1$.

Table 1 presents values of $a_0, a_1, b_0$ and external amplitude $f$ as functions of the friction coefficient $\delta$ obtained from the algebraic system.
<table><caption>Table 1
Values of $a_0$, $a_1$, $b_0$ and external amplitude $f$ for the Duffing equation ($\omega = 1$)</caption><thead><tr><th>δ</th><th>a<sub>0</sub></th><th>a<sub>1</sub></th><th>f</th><th>b<sub>0</sub></th></tr></thead><tbody><tr><td>0.0001</td><td>1.21498</td><td>0.621807</td><td>8.76005e-005</td><td>5.81771e-005</td></tr><tr><td>0.0002</td><td>1.21499</td><td>0.621808</td><td>0.000175206</td><td>0.000116352</td></tr><tr><td>0.0005</td><td>1.21503</td><td>0.621803</td><td>0.000438083</td><td>0.00029088</td></tr><tr><td>0.0007</td><td>1.21505</td><td>0.621813</td><td>0.000613318</td><td>0.000407203</td></tr><tr><td>0.001</td><td>1.21508</td><td>0.621819</td><td>0.000876227</td><td>0.000581683</td></tr><tr><td>0.002</td><td>1.21519</td><td>0.621834</td><td>0.00175293</td><td>0.00116316</td></tr><tr><td>0.003</td><td>1.2153</td><td>0.621848</td><td>0.00263011</td><td>0.00174445</td></tr><tr><td>0.004</td><td>1.21541</td><td>0.621862</td><td>0.00350778</td><td>0.00232553</td></tr><tr><td>0.005</td><td>1.21553</td><td>0.621876</td><td>0.00438592</td><td>0.00290641</td></tr><tr><td>0.006</td><td>1.21564</td><td>0.621889</td><td>0.00526454</td><td>0.0034871</td></tr><tr><td>0.007</td><td>1.21575</td><td>0.621903</td><td>0.00614365</td><td>0.00406758</td></tr><tr><td>0.008</td><td>1.21586</td><td>0.621916</td><td>0.00702323</td><td>0.00464787</td></tr><tr><td>0.009</td><td>1.21597</td><td>0.62193</td><td>0.00790329</td><td>0.00522795</td></tr><tr><td>0.01</td><td>1.21609</td><td>0.621943</td><td>0.00878383</td><td>0.00580784</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>Fig. 1. Example of HT in phase space of the non-autonomous Duffing equation calculated by computer with initial amplitude values obtained from the algebraic system (12), (13), (15) and (16). Here $\delta = 0.001$, $a_0 = 1.21508$, $a_1 = 0.621819$, $b_0 = 0.00058$, $f = 0.00087$.</figcaption></figure>

If we use an initial point ($a_0, a_1$) obtained from the algebraic system, a numerical calculation of the non-autonomous Duffing equation HT by using the Runge-Kutta procedure shows a good correspondence of the analytical and numerical results.

Some examples of HT in phase space calculated by computer with initial amplitude values obtained from the algebraic system (15), (16), (18), (19) are shown in Fig. 1 (here $\delta = 0.001$, $a_0 = 1.21508$, $a_1 = 0.621819$, $b_0 = 0.00058$, $f = 0.00087$) and in Fig. 2 (here $\delta = 0.01$, $a_0 = 1.21609$, $a_1 = 0.621943$, $b_0 = 0.0058$, $f = 0.00878$). Fig. 3 (here $\delta = 0.001$, $a_0 = 1.21508$, $a_1 = 0.621819$, $b_0 = 0.00058$, $f = 0.00087$) and Fig. 4 (here $\delta = 0.01$, $a_0 = 1.21609$, $a_1 = 0.621943$, $b_0 = 0.0058$, $f = 0.00878$) represent a comparison of the trajectories constructed using the Runge-Kutta procedure and using the QPA (17). Fig. 5 shows boundaries of the chaotic behavior region in parameter place ($\delta, f$) obtained by the proposed here method (curve a) and by the Melnikov one (curve b).

It is essentially that the boundary obtained by procedure proposed in the present paper is closer to a real situation that the corresponding boundary obtained by the Melnikov method where a separatrix trajectory of the autonomous Duffing equation is used.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Example of HT in phase space of the non-autonomous Duffing equation calculated by computer with initial amplitude values obtained from the algebraic system (12), (13), (15) and (16). Here $\delta = 0.01$, $a_0 = 1.21609$, $a_1 = 0.621943$, $b_0 = 0.0058$, $f = 0.00878$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Comparison of the HT of the non-autonomous Duffing equation constructed using the Runge-Kutta procedure (curve a) and using the QPA (14) (curve b). Here $\delta = 0.001$, $a_0 = 1.21508$, $a_1 = 0.621819$, $b_0 = 0.00058$, $f = 0.00087$.</figcaption></figure>

## 5. Equation of a parametrically excited damped pendulum

A damped parametrically driven mathematical pendulum can be described by the following equation [17]:

$$x'' + \delta x' + (1 - f \cos \omega t) \sin x = 0. \quad (20)$$

Such an external perturbation which acts upon a parameter—the amplitude of the back-driving force—is easy to realize by a periodic vertical displacement of the suspension point of the pendulum.

Let us consider the following equation

$$y'' + \delta y' - (1 - f \cos \omega t) \sin y = 0 \quad (21)$$

obtained from Eq. (20) after a shift $y = x + \pi$.

In this case we apply proposed approach for a construction of heteroclinic trajectory corresponding to the saddle points $(2\pi n, 0)$ $(n \in Z)$ in the phase space.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Comparison of the HT of the non-autonomous Duffing equation constructed using the Runge-Kutta procedure (curve a) and using the QPA (14) (curve b). Here $\delta = 0.01$, $a_0 = 1.21609$, $a_1 = 0.621943$, $b_0 = 0.0058$, $f = 0.00878$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Boundaries of the chaotic behavior region in parameter space $(\delta, f)$ of the non-autonomous Duffing equation obtained by the proposed here method (curve a) and by the Melnikov one (curve b).</figcaption></figure>

Like for the non-autonomous Duffing equation we do the same assumptions and we obtain the first algebraic equation by multiplying by $y'$ and integrating within limits from $t = 0$ to $t = \infty$:

$$ 
\begin{aligned}
\int_{0}^{\infty} (y'' + \delta y' - (1 - f \cos wt) \sin y)y' dt &= -\frac{a_1^2}{2} + \cos(b_0) - \cos(a_0) + \int_{0}^{\infty} (\delta y' + f \cos wt \sin y)y' dt \\
&= -\frac{a_1^2}{2} + \cos(b_0) - \cos(a_0) + (At + Bt^2 + Ct^3 + Dt^4 + Et^5 + \dots)|_{0}^{\infty} \\
&= -\frac{a_1^2}{2} + \cos(b_0) - \cos(a_0) + \left( \frac{\alpha_1 t + \alpha_2 t^2 + \alpha_3 t^3}{1 + \beta_1 t + \beta_2 t^2 + \beta_3 t^3} \right)|_{0}^{\infty} \\
&= -\frac{a_1^2}{2} + \cos(b_0) - \cos(a_0) + \frac{\alpha_3}{\beta_3} = 0. \qquad (22)
\end{aligned}
$$
Thus the first equation is

$$ -\frac{a_1^2}{2} + \cos(b_0) - \cos(a_0) + \frac{\alpha_3}{\beta_3} = 0. \qquad (23) $$

Reorganization of the local expansion of the sought solution near zero $y = a_0 + a_1 t + a_2 t^2 + a_3 t^3 + a_4 t^4 + a_5 t^5 + a_6 t^6 + \dots$ to the QPA takes an analytical continuation of the expansion ad infinitum

$$ y = a_0 + a_1 t + a_2 t^2 + a_3 t^3 + a_4 t^4 + \dots \approx e^{-t} \frac{\alpha_0 + \alpha_1 e^t + \alpha_2 e^{2t} + \alpha_3 e^{3t}}{1 + \beta_1 e^t + \beta_2 e^{2t}} \qquad (24) $$

and the second equation can be obtained from here as a condition at infinity:

$$ b_0 = \frac{\alpha_3}{\beta_2}. \qquad (25) $$

Then one can obtain equations from the convergence condition (3) realized for PA from (22) and for QPA from (24):

$$ -2D^2BF + 2DBE^2 - 3D^2CE + 2DC^2F + D^4 - AE^3 + E^2C^2 - 2ECBF + 2EDFA - F^2AC + F^2B^2 + GAEC - GAD^2 - GB^2E + 2GBDC - GC^3 = 0, \qquad (26) $$

<table><caption>Table 2 Values of $a_0$, $a_1$, $b_0$ and external amplitude $f$ for the equation of a parametrically excited damped pendulum ($\omega = 1$)</caption><thead><tr><th>δ</th><th>a<sub>0</sub></th><th>a<sub>1</sub></th><th>f</th><th>b<sub>0</sub></th></tr></thead><tbody><tr><td>0.0001</td><td>1.26091725</td><td>-1.17885213</td><td>-0.00018061</td><td>0.00105889</td></tr><tr><td>0.0002</td><td>1.35528225</td><td>-1.25372099</td><td>-0.00040509</td><td>0.00171579</td></tr><tr><td>0.0005</td><td>1.51238141</td><td>-1.37218232</td><td>-0.00112148</td><td>0.00357714</td></tr><tr><td>0.0007</td><td>1.57849101</td><td>-1.41959572</td><td>-0.00159115</td><td>0.00476448</td></tr><tr><td>0.0009</td><td>1.63043627</td><td>-1.45579456</td><td>-0.00204328</td><td>0.00591774</td></tr><tr><td>0.001</td><td>1.65280848</td><td>-1.47109306</td><td>-0.00226211</td><td>0.00648265</td></tr><tr><td>0.0015</td><td>1.74179994</td><td>-1.53017077</td><td>-0.00328527</td><td>0.00920249</td></tr><tr><td>0.002</td><td>1.80733760</td><td>-1.57182762</td><td>-0.00420246</td><td>0.01176999</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Comparison of HT (curve a) obtained from the system (18) with initial amplitude values obtained from the algebraic system (20), (22), (23) and (24) with two other trajectories closed to it in initial time moment $t = 0$ (curves b and c). Here $\delta = 0.001$, $a_0 = 1.75019$, $a_1 = -1.53548$, $b_0 = 0.0098$, $f = -0.004073$.</figcaption></figure>
$$ 
\begin{aligned}
& 60a_2a_4^2 + 360a_4a_5a_2 + 60a_3a_5a_2 - 360a_3a_4^2 + \frac{4}{3}a_4a_1^2 - \frac{1}{4}a_1a_2^2 - 3a_1a_3^2 + 10a_6a_1^2 + a_3a_2^2 - 60a_4a_3^2 + 360a_5a_3^2 \\
& + 7a_2a_3^2 + 360a_6a_1a_4 - 30a_5a_2^2 - 60a_6a_2^2 + \frac{1}{60}a_1^3 - \frac{1}{4}a_2^3 + a_5a_1^2 - 360a_1a_5^2 + \frac{11}{180}a_2a_1^2 + \frac{1}{4}a_3a_1^2 - 720a_6a_3^2 \\
& - 720a_2a_5^2 + 3a_4a_2^2 + 30a_3a_4a_2 - 360a_3a_6a_2 - 720a_4^3 + 720a_6a_4a_2 + 1440a_3a_4a_5 - 13a_1a_5a_2 - 60a_5a_4a_1 \\
& + a_1a_4a_2 + 60a_3a_6a_1 - 7a_3a_4a_1 + 30a_6a_1a_2 - 30a_3a_5a_1 - \frac{13}{12}a_3a_1a_2 = 0.
\end{aligned} 
\quad (27) $$

Eqs. (23), and (25)–(27) form a system of four non-linear algebraic equations for calculation $a_0, a_1, b_0$ and $f = f(\delta)$ while $\omega$ is fixed.

Table 2 presents values of $a_0, a_1, b_0$ and external amplitude $f$ as functions of the friction coefficient $\delta$ obtained from the system using Newton method while $\omega = 1$.

Fig. 6 presents a comparison of HT (curve a) obtained from the system (21) with initial amplitude values obtained from the algebraic system (23), (25)–(27) with two other trajectories closed to it in initial time moment $t = 0$ (curves b and c). Here $\delta = 0.001, a_0 = 1.75019, a_1 = -1.53548, b_0 = 0.0098, f = -0.004073$.

## 6. Concluding remarks

The methodology presented in this work is sufficiently general to be applicable to other types of non-linear systems of dimensions equal to two or three or more.

To apply the approach one needs to continue local homoclinic or heteroclinic approximations at $t = 0$ to infinity by using quasi-Pade' approximations. Simultaneously we consider a phase trajectory of the solution which can be obtained if we multiply the original equation by $y'(t)$ and then integrate along the trajectory of the solution within limits from zero to infinity. In the relation obtained after integration, an analytical continuation of local expansion to infinity is made by using the Pade' approximation. Besides, the necessary conditions of convergence of the PA or QPA are utilized, and as a result it is obtained some system of algebraic equations to calculate relevant parameters of HT, including initial values. It is important that it is not necessary to use here a separatrix trajectory of the autonomous Duffing equation.

Checking numerical calculations of the HT wife initial amplitudes values obtained as solutions of the algebraic equations show an acceptable precision of the proposed analytical procedure. Moreover, boundaries of the chaotic behavior region in parameter place $(\delta, f)$ obtained by the proposed approach is closer to the known real situation than boundaries obtained by using the Melnikov method.

## Acknowledgements

The authors would like to thank Professor Leonid I. Manevitch of Chemical Physics Institute, Russian Academy of Science, Moscow, Russia and Dr. Komstantin V. Avramov of Kharkov Polytechnical University, Ukraine for valuable discussions and suggestions on the subject.

## References

[1] Guckenheimer J, Holmes P. Nonlinear oscillations, dynamical systems and bifurcations of vector fields. New York: Springer-Verlag; 1993.

[2] Wiggins S. Introduction to applied nonlinear dynamical systems and chaos. New York: Springer-Verlag; 1990.

[3] Melnikov VK. On the stability of the centre for time periodic perturbations. Trans Moscow Math Soc 1963;12:1–57.

[4] Holmes PJ. A nonlinear oscillator with a strange attractor. Philos Trans Roy Soc A 1978;292:419–48.

[5] Moon FC. Chaotic vibrations. New York: Wiley; 1987.

[6] Sanders JA. Melnikov's method and averaging. Celestial Mech 1982;28:171–81.

[7] Martin P, Baker GA. Two-point quasifractional approximant in physics, truncation error. J Math Phys 1991;32:1470.

[8] Manevich LI, Mikhlin YuV, Pilipchuk VN. The method of normal oscillation for essentially nonlinear systems. Moscow: Nauka; 1989.

[9] Mikhlin YuV. Matching of local expansions in the theory of nonlinear vibrations. J Sound Vibration 1995;182:577–88.
[10] Vakakis AF, Manevick LI, Mikhlin YuV, Pilipchuk VN, Zevin AA. Normal modes and localization in nonlinear systems. New York: Wiley; 1996.

[11] Emaci E, Vakakis AF, Andrianov IV, Mikhlin YuV. Study of two-dimensional axisymmetric breathers using Pade' approximants. Nonlinear Dyn 1997;13:327-38.

[12] Vakakis AF, Azeez FA. Analytic approximation of the homoclinic orbits of the Lorenz system. Nonlinear Dyn 1998;15:245-57.

[13] Mikhlin YuV. Analytical construction of homoclinic orbits of two- and three-dimensional dynamical systems. J Sound Vibration 2000;230(5):971-83.

[14] Mikhlin YuV, Manucharyan GV, Savenkova SN. Analysis of homoclinic trajectories in mechanical systems with several equilibrium positions, Proceedings of the Sixth Conference on Dynamical Systems. Theory and Applications, Lodz, 2001, p. 309-16.

[15] Baker GA, Graves-Morris P. Pade' approximants. London: Addison-Wesley; 1981.

[16] Molotkcv IA, Vakulenko SA. Localized nonlinear waves. Leningrad: Leningrad University Press; 1988.

[17] Leven RW, Koch BP. Chaotic behavior of a parametrically excited damped pendulum. Phys Lett 1981;86A(2):71-4.