# Chaotic Behavior of the Biharmonic  Dynamics System

Authors: Vladimir S. Aslanov
Journal: International Journal of Mathematics and Mathematical Sciences
 (2009)

---

## Abstract

Motion of a biharmonic system under action of small periodic force and small damped force is studied. The biharmonic oscillator is a physical system acting under a biharmonic force likeasinθ+bsin2θ. The article contains biharmonic oscillator analysis, phase space research, and analytic solutions for separatrixes. The biharmonic oscillator performs chaotic motion near separatrixes under small perturbations. Melnikov method gives analytical criterion for heteroclinic chaos in terms of system parameters. A transition from chaotic to regular motion of the biharmonic oscillator was found as the heteroclinic chaos can be removed by increasing the coefficient of a damping force. The analytical results obtained using Melnikov method have been confirmed by a good match  with numeric research.


---

## Body

## 1. Introduction

### 1.1. Reduction of the Biharmonic Dynamics System to the Duffing Equation

As is well known,  nonlinear systems can perform a chaotic motion under the action of periodic forces [1–4]. Frequently the Duffing equation is used to illustrate chaos [1–4], and the chaotic behavior of various forms of the Duffing equation [5], some of which exhibit two-frequency excitation [6] as well as the chaotic motion of Duffing system with bounded noise [7] have been investigated. However the Duffing equation in the expanded (generalized) form has many mechanical applications and it can be interesting to researchers.In this paper, we will study periodically driven biharmonic dynamic system with a damping force:(1.1)θ¨=asinθ+bsin2θ+ε(a1sinθ+b1sin2θ+c)cosωt−δθ˙,  where ε and δ are assumed to be small positive parameters;  ω>0 is the frequency of the external force; and a, a1, b, b1, and c are coefficients. The terms in ε and in δ in (1.1) can be considered as small perturbations.Ifε=0 and δ=0, then the periodic and the damping force are absent, and we have the conservative system describing the motion of the undisturbed biharmonic oscillator as (1.2)θ¨=asinθ+bsin2θ. We will show that (1.1) is equivalent to the Duffing equation for small values of θ. It is known that (1.3)sinx=x−x33!+⋯(x=θ,2θ), and (1.1) can be written as (1.4)θ¨−a(θ−θ33!)−b[2θ−(2θ)33!]=ε{a(θ−θ33!)+b[2θ−(2θ)33!]+c}cosωt−δθ˙+⋯. If the variable θ is a small of order ε, then we have the Duffing equation (1.5)θ¨+λθ+μθ3=ε(c−λθ)cosωt−δθ˙+O(ε4), where λ=−(a+2b), μ=(a+8b)/3!.


### 1.2. Mechanical Applications

The disturbed system (1.1) has applications in space flight mechanics when studying a problem of a spacecraft motion about its center of mass in the atmosphere. Atmospheric re-entry is a critical phase for space vehicles. Dynamic stability issues play a crucial role for the success of their mission. For effective breaking while descending in the rarefied atmosphere of Mars or Titan we can compensate blunt-shaped spacecrafts by a small increase in length [8–11]. We will consider blunt-shaped spacecraft. An aerodynamic restoring moment strongly influences the motion of the spacecraft relative to the center of mass. The aerodynamic restoring moment coefficient for an axisymmetric rigid body can be written as [12, 13] (1.6)mθ(θ)=mO(θ)+cn(θ)x¯c, where θ is the spatial angle of attack  (θ is defined as the angle between symmetry axis and the velocity of the spacecraft V¯, Figure 1), mO is the pitching moment coefficient concerning the leading edge body, cn is the normal force coefficient, x¯c=xc/L is the nondimensional coordinate of the center of mass, and L is aerodynamic reference length (shield diameter). Such spacecrafts can have three positions of equilibrium according to the angle of attack θ: stable position at the points θ*=0 and θ*=π; and unstable in the third intermediate point θ*∈(0,π) [12, 13]. The restoring moment can be approximated as the biharmonic dependence (Figure 1) (1.7)mθ=asinθ+bsin2θ.
The presence of the second harmonic in (1.7) causes the possibility of appearance of an additional equilibrium position—saddle point on a phase portrait. For the considered spacecrafts position, θ=0 is stable; therefore, a derivative of the function mθ(θ) with respect to the angle of attack θ at this point is negative: (1.8)mθ′|θ=0=(acosθ+2bcos2θ)|θ=0<0 or (1.9)2b<−a, and if there exists an intermediate position of equilibrium inside the interval of (0,π), then (1.10)mθ(θ)=sinθ(a+2bcosθ)=0, which holds true, if  (1.11)|2b|>|a|. It is obvious that (1.9) and (1.10) are valid simultaneously when b<0. Note that the dependence of mθ(θ) given in Figure 1 satisfies conditions (1.9) and (1.10).Figure 1
A blunt-shaped spacecraft.The stable position occurs not only in the point ofθ=0 but also in the point of θ=π when (1.9) for the spacecraft is fulfilled. The motion of the spacecraft in a neighborhood of θ=π cannot be allowed, because in this case the back part of the spacecraft will move toward to an approach flow.The aerodynamic restoring moment can be written by(1.12)Mθ=mθ(θ)qSL, where S is reference surface and q is the dynamic pressure. Disturbance will simulate periodic change of position of the center of mass (1.13)x¯c=x¯c0+εcosωt, where x¯c0 is an initial position of the center of mass. Using expression (1.6), we represent the aerodynamic restoring moment as (1.14)Mθ=Asinθ+Bsin2θ+ε(Asinθ+Bsin2θ)cosωt. Thus, the planar motion of the spacecraft about the center of mass can be described by (1.15)Iθ¨=Asinθ+Bsin2θ+ε(Asinθ+Bsin2θ)cosωt+Md, where I is the transverse moment of inertia, Md=−δ1θ˙ is small damping moment. Obviously, the equation of spacecraft motion (1.15) corresponds to the biharmonic dynamics system (1.1) at c=0 up to constants.Let us observe, that the undisturbed equation (1.2) describes the motion of a known mechanical system—a heavy material point on a circle, rotating about a vertical axis [14] (1.16)θ¨=asinθ+bsin2θ(a=−gl,b=Ω2>0), where g is the gravitational acceleration, l is the radius of the circle and Ω is the angular velocity of the circle.


### 1.3. Aim and Structure of the Article

Let conditions (1.9) and (1.11) be satisfied, and let there be three positions of equilibrium. In this case, we find three regions in the phase portrait separated by separatrixes. Under the effect of disturbances the phase trajectory θ˙=θ˙(θ) can repeatedly intersect the separatrixes, thus moving from one an area to another, is accompanied by a jump change of the variable θ. We can observe chaos. It is an accepted fact in the theory of nonlinear dynamic systems that knowledge of the stable and unstable manifolds of hyperbolic equilibrium or hyperbolic periodic orbits may play a crucial role in understanding many issues of dynamics [1, 2]. For many dynamic systems, the only general way of studying such stable and unstable manifolds is computing them numerically. However, in some cases we can obtain analytic solutions.The aim of this paper is to analyze the motion of the disturbed system (1.1) near the undisturbed separatrixes and to define the boundaries of chaos. We will carry out the theoretical studies by means of the Melnikov method [15]. The Melnikov method is an analytical tool used to determine to first order, the existence of homo/heteroclinic intersections and so chaotic behavior. The Melnikov method allows us to obtain a necessary condition for the existence of chaos, therefore, numerical simulation is needed to confirm the predicted behavior and to give a deeper understanding of the global dynamics of the system. We will show that our theoretical results are in good agreement with the results of numerical calculations.The present paper is structured in the following way. Section2 gives the analysis of the unperturbed motion of the biharmonic dynamics system (1.1) and the phase portrait. On the phase portrait characteristic regions of the possible motions are found and for these areas the analytical solutions of the equation of the unperturbed motion are obtained. The main features of the phase space of the unperturbed system are defined. In Section 3 the Melnikov criterion for the perturbed system is analytically calculated for various areas of the phase portrait with help of the theory of residues. In Section 4 the Melnikov criterion is numerically calculated for various areas of a phase portrait for the case of disturbed motion of the spacecraft. By means of computer numerical simulations of the disturbed motion, we use several numerical techniques to check the validity of the analytical criterion for chaos obtained using Melnikov method.


## 1.1. Reduction of the Biharmonic Dynamics System to the Duffing Equation

As is well known,  nonlinear systems can perform a chaotic motion under the action of periodic forces [1–4]. Frequently the Duffing equation is used to illustrate chaos [1–4], and the chaotic behavior of various forms of the Duffing equation [5], some of which exhibit two-frequency excitation [6] as well as the chaotic motion of Duffing system with bounded noise [7] have been investigated. However the Duffing equation in the expanded (generalized) form has many mechanical applications and it can be interesting to researchers.In this paper, we will study periodically driven biharmonic dynamic system with a damping force:(1.1)θ¨=asinθ+bsin2θ+ε(a1sinθ+b1sin2θ+c)cosωt−δθ˙,  where ε and δ are assumed to be small positive parameters;  ω>0 is the frequency of the external force; and a, a1, b, b1, and c are coefficients. The terms in ε and in δ in (1.1) can be considered as small perturbations.Ifε=0 and δ=0, then the periodic and the damping force are absent, and we have the conservative system describing the motion of the undisturbed biharmonic oscillator as (1.2)θ¨=asinθ+bsin2θ. We will show that (1.1) is equivalent to the Duffing equation for small values of θ. It is known that (1.3)sinx=x−x33!+⋯(x=θ,2θ), and (1.1) can be written as (1.4)θ¨−a(θ−θ33!)−b[2θ−(2θ)33!]=ε{a(θ−θ33!)+b[2θ−(2θ)33!]+c}cosωt−δθ˙+⋯. If the variable θ is a small of order ε, then we have the Duffing equation (1.5)θ¨+λθ+μθ3=ε(c−λθ)cosωt−δθ˙+O(ε4), where λ=−(a+2b), μ=(a+8b)/3!.


## 1.2. Mechanical Applications

The disturbed system (1.1) has applications in space flight mechanics when studying a problem of a spacecraft motion about its center of mass in the atmosphere. Atmospheric re-entry is a critical phase for space vehicles. Dynamic stability issues play a crucial role for the success of their mission. For effective breaking while descending in the rarefied atmosphere of Mars or Titan we can compensate blunt-shaped spacecrafts by a small increase in length [8–11]. We will consider blunt-shaped spacecraft. An aerodynamic restoring moment strongly influences the motion of the spacecraft relative to the center of mass. The aerodynamic restoring moment coefficient for an axisymmetric rigid body can be written as [12, 13] (1.6)mθ(θ)=mO(θ)+cn(θ)x¯c, where θ is the spatial angle of attack  (θ is defined as the angle between symmetry axis and the velocity of the spacecraft V¯, Figure 1), mO is the pitching moment coefficient concerning the leading edge body, cn is the normal force coefficient, x¯c=xc/L is the nondimensional coordinate of the center of mass, and L is aerodynamic reference length (shield diameter). Such spacecrafts can have three positions of equilibrium according to the angle of attack θ: stable position at the points θ*=0 and θ*=π; and unstable in the third intermediate point θ*∈(0,π) [12, 13]. The restoring moment can be approximated as the biharmonic dependence (Figure 1) (1.7)mθ=asinθ+bsin2θ.
The presence of the second harmonic in (1.7) causes the possibility of appearance of an additional equilibrium position—saddle point on a phase portrait. For the considered spacecrafts position, θ=0 is stable; therefore, a derivative of the function mθ(θ) with respect to the angle of attack θ at this point is negative: (1.8)mθ′|θ=0=(acosθ+2bcos2θ)|θ=0<0 or (1.9)2b<−a, and if there exists an intermediate position of equilibrium inside the interval of (0,π), then (1.10)mθ(θ)=sinθ(a+2bcosθ)=0, which holds true, if  (1.11)|2b|>|a|. It is obvious that (1.9) and (1.10) are valid simultaneously when b<0. Note that the dependence of mθ(θ) given in Figure 1 satisfies conditions (1.9) and (1.10).Figure 1
A blunt-shaped spacecraft.The stable position occurs not only in the point ofθ=0 but also in the point of θ=π when (1.9) for the spacecraft is fulfilled. The motion of the spacecraft in a neighborhood of θ=π cannot be allowed, because in this case the back part of the spacecraft will move toward to an approach flow.The aerodynamic restoring moment can be written by(1.12)Mθ=mθ(θ)qSL, where S is reference surface and q is the dynamic pressure. Disturbance will simulate periodic change of position of the center of mass (1.13)x¯c=x¯c0+εcosωt, where x¯c0 is an initial position of the center of mass. Using expression (1.6), we represent the aerodynamic restoring moment as (1.14)Mθ=Asinθ+Bsin2θ+ε(Asinθ+Bsin2θ)cosωt. Thus, the planar motion of the spacecraft about the center of mass can be described by (1.15)Iθ¨=Asinθ+Bsin2θ+ε(Asinθ+Bsin2θ)cosωt+Md, where I is the transverse moment of inertia, Md=−δ1θ˙ is small damping moment. Obviously, the equation of spacecraft motion (1.15) corresponds to the biharmonic dynamics system (1.1) at c=0 up to constants.Let us observe, that the undisturbed equation (1.2) describes the motion of a known mechanical system—a heavy material point on a circle, rotating about a vertical axis [14] (1.16)θ¨=asinθ+bsin2θ(a=−gl,b=Ω2>0), where g is the gravitational acceleration, l is the radius of the circle and Ω is the angular velocity of the circle.


## 1.3. Aim and Structure of the Article

Let conditions (1.9) and (1.11) be satisfied, and let there be three positions of equilibrium. In this case, we find three regions in the phase portrait separated by separatrixes. Under the effect of disturbances the phase trajectory θ˙=θ˙(θ) can repeatedly intersect the separatrixes, thus moving from one an area to another, is accompanied by a jump change of the variable θ. We can observe chaos. It is an accepted fact in the theory of nonlinear dynamic systems that knowledge of the stable and unstable manifolds of hyperbolic equilibrium or hyperbolic periodic orbits may play a crucial role in understanding many issues of dynamics [1, 2]. For many dynamic systems, the only general way of studying such stable and unstable manifolds is computing them numerically. However, in some cases we can obtain analytic solutions.The aim of this paper is to analyze the motion of the disturbed system (1.1) near the undisturbed separatrixes and to define the boundaries of chaos. We will carry out the theoretical studies by means of the Melnikov method [15]. The Melnikov method is an analytical tool used to determine to first order, the existence of homo/heteroclinic intersections and so chaotic behavior. The Melnikov method allows us to obtain a necessary condition for the existence of chaos, therefore, numerical simulation is needed to confirm the predicted behavior and to give a deeper understanding of the global dynamics of the system. We will show that our theoretical results are in good agreement with the results of numerical calculations.The present paper is structured in the following way. Section2 gives the analysis of the unperturbed motion of the biharmonic dynamics system (1.1) and the phase portrait. On the phase portrait characteristic regions of the possible motions are found and for these areas the analytical solutions of the equation of the unperturbed motion are obtained. The main features of the phase space of the unperturbed system are defined. In Section 3 the Melnikov criterion for the perturbed system is analytically calculated for various areas of the phase portrait with help of the theory of residues. In Section 4 the Melnikov criterion is numerically calculated for various areas of a phase portrait for the case of disturbed motion of the spacecraft. By means of computer numerical simulations of the disturbed motion, we use several numerical techniques to check the validity of the analytical criterion for chaos obtained using Melnikov method.


## 2. The Unperturbed Solutions

Theε and δ terms in (1.1) are considered as small perturbations. If ε=0 and δ=0, then periodic and dissipation forces are absent, and we have the conservative system (1.2). It is obvious that if b=0 or a=0 (replacement of variables φ=2θ) we have the equation of a mathematical pendulum. However, if (1.9) and (1.11) are satisfied (b<0), then biharmonic dynamic system has a more complicated phase portrait in comparison with the mathematical pendulum or with the heavy material point on a circle (1.16). The equilibrium positions of the system (1.2) are defined from (1.10). If (1.9) and (1.11) are satisfied, then the undisturbed system (1.2) has four equilibrium positions at θ∈[−π,π]: two stable—center type  (2.1)θ=0,π and two unstable—saddle type (2.2)θ∗=±arccos(−a2b), where b<0. The center θ∗=−π coincides with the center θ∗=π. At θ∗→−π and at θ∗→π the speeds θ˙ coincide, therefore we can say, that phase trajectories are closed on a cylindrical phase space. We will consider the evolution of the cylindrical space in the range θ∈[−π,π]. We will separate two regions A0 and A1, divided by the two saddles s1 and s−1 (Figure 2). It is necessary to note, that the region A1 of the development of the cylinder undergoes a break at θ=π,−π. From (2.2) it follows, that if the coefficient a is equal 0, the saddle s1 is in the position: θ∗=π/2. At positive values of the coefficient a>0 the saddle s1 belongs to the interval: θ∗∈(0,π/2), and at negative values a<0 the saddle s1 belongs to the interval: θ∗∈(π/2,π) (Figure 2).The potential energyW(θ)=acosθ+bcos2θ and the phase space for a=1, b=−1.
(a)(b)The following energy integral corresponds to (1.2): (2.3)12θ˙2+acosθ+bcos2θ=E, where E is total energy. The biharmonic oscillator as well as the mathematical pendulum can perform oscillations and rotation.  The shape of the phase portrait depends on the potential energy: (2.4)W(θ)=acosθ+bcos2θ. The centers (2.1) correspond to the minimum of the potential energy (2.4), and the saddles (2.2)—to the maximum of the potential energy (2.4).  If E>W∗, where W∗=W(θ∗), then the motion is possible in the outer regions (Figure 2). In the opposite case (E<W∗) the motion can occur in any of the inner regions, depending on initial conditions. The equality E=W∗ corresponds to the motion along separatrixes. In this case, the two saddles s1 and s−1 are connected by four heteroclinic trajectories.First of all, we will consider the separatrixes, limiting the regionA0. Separating the variables in the energy integral (2.3), and taking into account (2.2) and (2.4), the equation of the motion on the separatrixes can be written as in the integrated form (2.5)t−t0=∫θ0θdθ2[W(θ∗)−acosθ−bcos2θ], where (2.6)W(θ∗)=acosθ∗+bcos2θ∗=−a24b. Substituting the variables: (2.7)x=tanθ2, we can rewrite (2.5) as (2.8)t−t0=2∫x0xdxψ(x),ψ(x)=Px4+Cx2+A, where P=2[W(θ∗)+a−b], C=4[W(θ∗)+b], A=2[W(θ∗)−a−b].The integral (2.8) can be simplified [16] (2.9)t−t0=2P∫x0xdxx∗2−x2=2Pln|x∗+xx∗−x||x0x, where x∗=tan(θ∗/2) and  P=2[W(θ∗)+a−b]=−(a−2b)2/2b>0 and conditions (1.9) and (1.11) are satisfied.Finally, using the change of variables (2.7) and (2.9), the solution of (1.2) for the heteroclinic orbits, for the region A0  (Figure 2), can be written as (2.10)θ+(t)=2arctan[tanθ∗2tanh(λt2)],σ+(t)=(θ˙)=λsinθ∗cosh(λt)+cosθ∗,[θ−(t),σ−(t)]=[−θ+(t),−σ+(t)], where (2.11)λ=tanθ∗2P=a2−4b22b. Now we will consider the region A1, including the center c1. Let us make a substitution to the new variable (2.12)β=π−θ in the equation of undisturbed motion (1.2) and obtain the following equation (2.13)β¨=−asinβ+bsin2β. Computing as in (2.3)–(2.9), we obtain the solution of this equation in the form (2.14)β=2arctan{tanβ∗2tanh[λ(t−t0)2]}, where β∗=π−θ∗. Then, coming back to the variable θ with help of (2.11), we will receive the equation of the heteroclinic orbits, bounding the region A1, including the center c1 (Figure 2): (2.15)θ+(t)=π−2arctan[cotθ∗2tanh(λt2)],σ+(t)=(θ˙)=λsinθ∗cosh(λt)−cosθ∗,[θ−(t),σ−(t)]=[2π−θ+(t),−σ+(t)].


## 3. Chaotic Motion. The Melnikov Criterion

### 3.1. General Positions

Now we set the stage for our study of the disturbed system (1.1). The stable and unstable manifolds do not necessarily coincide and it is possible that they can cross transversally leading to an infinite number of new heteroclinic points. Then, a heteroclinic tangle is generated. In this case, because of the perturbation, the motion of the system (1.1), near the unperturbed separatrices, becomes chaotic. Inside this chaotic layer small isolated regions of regular motion with periodic orbits can also appear. The existence of heteroclinic intersections may be proved by means of the Melnikov method [15]. We present a more convenient form the for application of Melnikov method to the disturbed nonautonomous equation of the second order (1.1) as three differential autonomous equations of the first order [2]: (3.1)θ˙=σ=f1+g1,σ˙=asinθ+bsin2θ+ε(a1sinθ+b1sin2θ+c)cosϕ−δσ=f2+g2,ϕ˙=ω, where (3.2)f1=σ,g1=0,f2=asinθ+bsin2θ,g2=ε(a1sinθ+b1sin2θ+c)cosϕ−δσ. The Melnikov function [2] for system (3.1) is given by (3.3)M±(t0,ϕ0)=∫−∞∞{f1[q±0(t)]g2[q±0(t),ωt+ωt0+ϕ0]−f2[q±0(t)]g1[q±0(t),ωt+ωt0+ϕ0]}dt=∫−∞∞{f1[q±0(t)]g2[q±0(t),ωt+ωt0+ϕ0]}dt, where q±0(t)=[θ±(t),σ±(t)] are the solutions of the undisturbed heteroclinic orbits (2.10) or (2.15)  for the areas  A0 or A1.


### 3.2. Case 1(a1=0,b1=0,c=1)

The disturbed system (3.1) in this case takes the form (3.4)θ˙=σ=f1+g1,σ˙=asinθ+bsin2θ+εcosϕ−δσ=f2+g2,ϕ˙=ω, where (3.5)f1=σ,g1=0,f2=asinθ+bsin2θ,g2=εcosϕ−δσ. Substituting (3.5) into (3.3) gives (3.6)M±(t0,ϕ0)=∫−∞∞σ±{εcos(ωt+ωt0+ϕ0)−δσ±}dt=ε∫−∞∞σ±cos(ωt+ωt0+ϕ0)dt−δ∫−∞∞(σ±)2dt=Mε+Mδ, where Mε and Mδ are the functions corresponding to both perturbations: the external periodic force (εcosωt) and the damping force (−δθ˙), respectively. The Melnikov function describes the splitting of the stable and unstable manifolds of the disturbed hyperbolic fixed points defined on the cross-section. Thus, there are transverse intersections between the stable and unstable trajectories, if M±(t0)=0.Firstly we consider the functionsMδ(0) and Mε(0) for the area A0, including the center c0 (Figure 2). Substituting (2.10) into (3.6) gives (3.7)Mδ(0)=−δ∫−∞∞(σ±)2dt=−δsin2θS∫−∞∞dt[cosh(λt)+cosθ∗]2,Mε(0)(t0,ϕ0)=ε∫−∞∞σ±cos(ωt+ωt0+ϕ0)dt=ελsinθS∫−∞∞cos(ωt+ωt0+ϕ0)cosh(λt)+cosθ∗dt. So using the tabulated integrals [16, 17], for the integrals (3.7) we obtain the following expressions (3.8)Mδ(0)=−2δλ(1−θ∗cotθ∗),Mε(0)(t0,ϕ0)=2επsinh(θ∗(ω/λ))λsin(θ∗)sinh(π(ω/λ))cos(ωt0+ϕ0)=Mεmax(0)cos(ωt0+ϕ0). Similar expressions can be obtained for Mδ(1) and Mε(1) for the region A1, including the center c1  (Figure 2), using the solutions (2.15) (3.9)Mδ(1)=−δ∫−∞∞(σ±)2dt=−δsin2θ∗∫−∞∞dt[cosh(λt)−cosθ∗]2,Mε(1)(t0,ϕ0)=ε∫−∞∞σ±cos(ωt+ωt0+ϕ0)dt=ελsinθS∫−∞∞cos(ωt+ωt0+ϕ0)cosh(λt)−cosθ∗dt, or (3.10)Mδ(1)=−2δλ[1+(π−θ∗)cotθ∗],Mε(1)(t0,ϕ0)=2επsinh[(π−θ∗)(ω/λ)]λsin(θ∗)sinh(π(ω/λ))cos(ωt0+ϕ0)=Mεmax(1)cos(ωt0+ϕ0), where Mεmax(1) and Mεmax(2) are measures of the maximum splitting of the stable and unstable manifolds, when the disturbed system (3.4) is only under the action of the one perturbation the external periodic force (εcosωt) for the regions A0 and A1, respectively.Obviously, ata=0 the undisturbed biharmonic oscillator (1.2) is transformed to the simpler system: θ¨=sin2θ. The regions A0 and A1 are equal. From (2.2), (2.10), and (2.15), we obtain (3.11)θ∗=π2,λ=−2b. Following the expressions (3.8)–(3.10), the Melnikov function becomes identical for the regions A0 and A1:
(3.12)M±(t0,ϕ0)=Mδ+Mε(t0,ϕ0)=−2δλ+επλsech(π2ωλ)cos(ωt0+ϕ0).From (3.8)–(3.10), it is easy to see that the conditions for the manifolds to intersect in terms of the parameters (δ,ε) is given by (3.13)δ<[πsinh(θ∗(ω/λ))λ2(1−θScotθ∗)sin(θ∗)sinh(π(ω/λ))]ε(fortheareaA0),δ<{πsinh[(π−θ∗)(ω/λ)]λ2[1+(π−θ∗)cotθ∗]sin(θ∗)sinh(π(ω/λ))}ε(fortheareaA1). Let us define a new parameter of the damping force, divided into amplitude of external force: (3.14)Δ=δε, then conditions (3.13) are given by (3.15)Δ<πsinh(θ∗(ω/λ))λ2(1−θ∗cotθ∗)sin(θ∗)sinh(π(ω/λ))=Δ0(fortheareaA0),Δ<πsinh[(π−θ∗)(ω/λ)]λ2[1+(π−θ∗)cotθ∗]sin(θ∗)sinh(π(ω/λ))=Δ1(fortheareaA1). Let us note, that θ∗ and λ, according to (2.2) and (2.11), depend on coefficients a and b, therefore criteria (3.15) are functions of the parameters a, b, and ω(3.16)Δj=Δj(a,b,ω),j=0,1. The criteria (3.16) define chaotic behavior of the perturbed system (3.4) in the regions A0 and A1. In Figure 3, we graph these criteria and the variable θ∗ as functions of parameter a for the fixed parameter b=−1 and ω=1. Figure 4 shows the criteria (3.16) as functions of the frequency ω.Figure 3
The criteriaΔj and the variable θ∗ of parameter a for b=−1 and ω=1.Figure 4
The criteriaΔj as functions of the frequency ω.


### 3.3. Case 2(a1≠0,b1≠0,c=0)

The disturbed system (3.1) in this case takes the form of (3.17)θ˙=σ=f1+g1,σ˙=asinθ+bsin2θ+ε(a1sinθ+b1sin2θ)cosϕ−δσ=f2+g2,ϕ˙=ω, where  f1=σ, g1=0, f2=asinθ+bsin2θ, g2=ε(a1sinθ+b1sin2θ)cosϕ−δσ, (3.18)M±(t0,ϕ0)=∫−∞∞σ±[ε(a1sinθ±+b1sin2θ±)cos(ωt+ωt0+ϕ0)−δσ±]dt or (3.19)M±(t0,ϕ0)=Mε+Mδ, where (3.20)Mε=ε∫−∞∞σ±(a1sinθ±+b1sin2θ±)cos(ωt+ωt0+ϕ0)dt,Mδ=−δ∫−∞∞(σ±)2dt. For the two regions A0 and A1 the functions (3.20) can be represented as (3.21)Mε(k)=−εI±(k)sin(ωt0+ϕ0),Mδ(k)=−δJ±(k),k=0,1, where (3.22)I±(k)=∫−∞∞σ±(k)(a1sinθ±(k)+b1sin2θ±(k))sin(ωt)dt,(3.23)J±(k)=∫−∞∞(σ±(k))2dt.  It is obvious that using (3.8) integrals (3.23) can be rewritten as
(3.24)J±(0)=2λ(1−θ∗cotθ∗),J±(1)=2λ[1+(π−θ∗)cotθ∗], The improper integral (3.22) in view of (2.10) and (2.15) is calculated numerically. For parameter (3.14) the conditions for the manifolds to intersect are given by (3.25)Δ<I±(0)J±(0)=Δ0(fortheareaA0),Δ<I±(1)J±(1)=Δ1(fortheareaA1). Criteria (3.25) define behavior of the perturbed system (3.17) in a vicinity of separatrixes.


## 3.1. General Positions

Now we set the stage for our study of the disturbed system (1.1). The stable and unstable manifolds do not necessarily coincide and it is possible that they can cross transversally leading to an infinite number of new heteroclinic points. Then, a heteroclinic tangle is generated. In this case, because of the perturbation, the motion of the system (1.1), near the unperturbed separatrices, becomes chaotic. Inside this chaotic layer small isolated regions of regular motion with periodic orbits can also appear. The existence of heteroclinic intersections may be proved by means of the Melnikov method [15]. We present a more convenient form the for application of Melnikov method to the disturbed nonautonomous equation of the second order (1.1) as three differential autonomous equations of the first order [2]: (3.1)θ˙=σ=f1+g1,σ˙=asinθ+bsin2θ+ε(a1sinθ+b1sin2θ+c)cosϕ−δσ=f2+g2,ϕ˙=ω, where (3.2)f1=σ,g1=0,f2=asinθ+bsin2θ,g2=ε(a1sinθ+b1sin2θ+c)cosϕ−δσ. The Melnikov function [2] for system (3.1) is given by (3.3)M±(t0,ϕ0)=∫−∞∞{f1[q±0(t)]g2[q±0(t),ωt+ωt0+ϕ0]−f2[q±0(t)]g1[q±0(t),ωt+ωt0+ϕ0]}dt=∫−∞∞{f1[q±0(t)]g2[q±0(t),ωt+ωt0+ϕ0]}dt, where q±0(t)=[θ±(t),σ±(t)] are the solutions of the undisturbed heteroclinic orbits (2.10) or (2.15)  for the areas  A0 or A1.


## 3.2. Case 1(a1=0,b1=0,c=1)

The disturbed system (3.1) in this case takes the form (3.4)θ˙=σ=f1+g1,σ˙=asinθ+bsin2θ+εcosϕ−δσ=f2+g2,ϕ˙=ω, where (3.5)f1=σ,g1=0,f2=asinθ+bsin2θ,g2=εcosϕ−δσ. Substituting (3.5) into (3.3) gives (3.6)M±(t0,ϕ0)=∫−∞∞σ±{εcos(ωt+ωt0+ϕ0)−δσ±}dt=ε∫−∞∞σ±cos(ωt+ωt0+ϕ0)dt−δ∫−∞∞(σ±)2dt=Mε+Mδ, where Mε and Mδ are the functions corresponding to both perturbations: the external periodic force (εcosωt) and the damping force (−δθ˙), respectively. The Melnikov function describes the splitting of the stable and unstable manifolds of the disturbed hyperbolic fixed points defined on the cross-section. Thus, there are transverse intersections between the stable and unstable trajectories, if M±(t0)=0.Firstly we consider the functionsMδ(0) and Mε(0) for the area A0, including the center c0 (Figure 2). Substituting (2.10) into (3.6) gives (3.7)Mδ(0)=−δ∫−∞∞(σ±)2dt=−δsin2θS∫−∞∞dt[cosh(λt)+cosθ∗]2,Mε(0)(t0,ϕ0)=ε∫−∞∞σ±cos(ωt+ωt0+ϕ0)dt=ελsinθS∫−∞∞cos(ωt+ωt0+ϕ0)cosh(λt)+cosθ∗dt. So using the tabulated integrals [16, 17], for the integrals (3.7) we obtain the following expressions (3.8)Mδ(0)=−2δλ(1−θ∗cotθ∗),Mε(0)(t0,ϕ0)=2επsinh(θ∗(ω/λ))λsin(θ∗)sinh(π(ω/λ))cos(ωt0+ϕ0)=Mεmax(0)cos(ωt0+ϕ0). Similar expressions can be obtained for Mδ(1) and Mε(1) for the region A1, including the center c1  (Figure 2), using the solutions (2.15) (3.9)Mδ(1)=−δ∫−∞∞(σ±)2dt=−δsin2θ∗∫−∞∞dt[cosh(λt)−cosθ∗]2,Mε(1)(t0,ϕ0)=ε∫−∞∞σ±cos(ωt+ωt0+ϕ0)dt=ελsinθS∫−∞∞cos(ωt+ωt0+ϕ0)cosh(λt)−cosθ∗dt, or (3.10)Mδ(1)=−2δλ[1+(π−θ∗)cotθ∗],Mε(1)(t0,ϕ0)=2επsinh[(π−θ∗)(ω/λ)]λsin(θ∗)sinh(π(ω/λ))cos(ωt0+ϕ0)=Mεmax(1)cos(ωt0+ϕ0), where Mεmax(1) and Mεmax(2) are measures of the maximum splitting of the stable and unstable manifolds, when the disturbed system (3.4) is only under the action of the one perturbation the external periodic force (εcosωt) for the regions A0 and A1, respectively.Obviously, ata=0 the undisturbed biharmonic oscillator (1.2) is transformed to the simpler system: θ¨=sin2θ. The regions A0 and A1 are equal. From (2.2), (2.10), and (2.15), we obtain (3.11)θ∗=π2,λ=−2b. Following the expressions (3.8)–(3.10), the Melnikov function becomes identical for the regions A0 and A1:
(3.12)M±(t0,ϕ0)=Mδ+Mε(t0,ϕ0)=−2δλ+επλsech(π2ωλ)cos(ωt0+ϕ0).From (3.8)–(3.10), it is easy to see that the conditions for the manifolds to intersect in terms of the parameters (δ,ε) is given by (3.13)δ<[πsinh(θ∗(ω/λ))λ2(1−θScotθ∗)sin(θ∗)sinh(π(ω/λ))]ε(fortheareaA0),δ<{πsinh[(π−θ∗)(ω/λ)]λ2[1+(π−θ∗)cotθ∗]sin(θ∗)sinh(π(ω/λ))}ε(fortheareaA1). Let us define a new parameter of the damping force, divided into amplitude of external force: (3.14)Δ=δε, then conditions (3.13) are given by (3.15)Δ<πsinh(θ∗(ω/λ))λ2(1−θ∗cotθ∗)sin(θ∗)sinh(π(ω/λ))=Δ0(fortheareaA0),Δ<πsinh[(π−θ∗)(ω/λ)]λ2[1+(π−θ∗)cotθ∗]sin(θ∗)sinh(π(ω/λ))=Δ1(fortheareaA1). Let us note, that θ∗ and λ, according to (2.2) and (2.11), depend on coefficients a and b, therefore criteria (3.15) are functions of the parameters a, b, and ω(3.16)Δj=Δj(a,b,ω),j=0,1. The criteria (3.16) define chaotic behavior of the perturbed system (3.4) in the regions A0 and A1. In Figure 3, we graph these criteria and the variable θ∗ as functions of parameter a for the fixed parameter b=−1 and ω=1. Figure 4 shows the criteria (3.16) as functions of the frequency ω.Figure 3
The criteriaΔj and the variable θ∗ of parameter a for b=−1 and ω=1.Figure 4
The criteriaΔj as functions of the frequency ω.


## 3.3. Case 2(a1≠0,b1≠0,c=0)

The disturbed system (3.1) in this case takes the form of (3.17)θ˙=σ=f1+g1,σ˙=asinθ+bsin2θ+ε(a1sinθ+b1sin2θ)cosϕ−δσ=f2+g2,ϕ˙=ω, where  f1=σ, g1=0, f2=asinθ+bsin2θ, g2=ε(a1sinθ+b1sin2θ)cosϕ−δσ, (3.18)M±(t0,ϕ0)=∫−∞∞σ±[ε(a1sinθ±+b1sin2θ±)cos(ωt+ωt0+ϕ0)−δσ±]dt or (3.19)M±(t0,ϕ0)=Mε+Mδ, where (3.20)Mε=ε∫−∞∞σ±(a1sinθ±+b1sin2θ±)cos(ωt+ωt0+ϕ0)dt,Mδ=−δ∫−∞∞(σ±)2dt. For the two regions A0 and A1 the functions (3.20) can be represented as (3.21)Mε(k)=−εI±(k)sin(ωt0+ϕ0),Mδ(k)=−δJ±(k),k=0,1, where (3.22)I±(k)=∫−∞∞σ±(k)(a1sinθ±(k)+b1sin2θ±(k))sin(ωt)dt,(3.23)J±(k)=∫−∞∞(σ±(k))2dt.  It is obvious that using (3.8) integrals (3.23) can be rewritten as
(3.24)J±(0)=2λ(1−θ∗cotθ∗),J±(1)=2λ[1+(π−θ∗)cotθ∗], The improper integral (3.22) in view of (2.10) and (2.15) is calculated numerically. For parameter (3.14) the conditions for the manifolds to intersect are given by (3.25)Δ<I±(0)J±(0)=Δ0(fortheareaA0),Δ<I±(1)J±(1)=Δ1(fortheareaA1). Criteria (3.25) define behavior of the perturbed system (3.17) in a vicinity of separatrixes.


## 4. Numerical Analysis

We have analyzed the evolution of the dynamical behavior of the disturbed system (1.1) as the parameters vary, studying the time histories of the variable θ and its derivative θ˙. Numerical techniques are based on the numerical integration of the equation of the disturbed motion (1.1) implementing a fixed-step fourth-order Runge-Kutta algorithm. For all numerical calculations the following biharmonic force parameters were used: a=1, b=−1 and the frequency of the perturbed force was ω=1. For the numerical analysis of the disturbed system (3.1) we use the P Poincaré cross-section method, examining manifolds with plane sections, perpendicular to the phase axis ϕ in the two-dimensional space (θ,θ˙), divided with an interval of 2π. It allows us to study the disturbed system (3.1) using a discrete phase instead of examining the continuous dynamics of the system. At ε=0, δ=0 the regular structure of phase space is observed, trajectories have no intersections, and Poincaré sections coincide with undisturbed phase portrait (Figure 5).Figure 5
Poincaré sections forε=0, δ=0.Disturbances (ε≠0) result in the complication of phase space and the occurrence of a chaotic layer near the undisturbed separatrixes (Figures 6–9). Figures 6 and 7 contain Poincaré sections for the case considered in Section 3.2, and in Figures 8 and 9—in Section 3.3. The growth of disturbances there leads to an increase in the width of the chaotic layer, and the new oscillatory modes determined by closed curves, uncharacteristic for the undisturbed case are observed.Figure 6
Poincaré sections in the case ofa1=0,b1=0,c=1 for ε=0.01, δ=0.Figure 7
Poincaré sections in the case ofa1=0,b1=0,c=1 for ε=0.02, δ=0.Figure 8
Poincaré sections in the case ofa1=1,b1=−1,c=0  for ε=0.01, δ=0.Figure 9
Poincaré sections in the case ofa1=1,b1=−1,c=0 for ε=0.02, δ=0.In the presence of damping phase trajectories eventually tended to reach steady positions of equilibrium of the undisturbed system (Figures10, 11).Figure 10
Poincaré sections in the case ofa1=0,b1=0,c=1 for ε=0.02, δ=0.0001.Figure 11
Poincaré sections in the case ofa1=1,b1=−1,c=0 for ε=0.02, δ=0.0001.In order to check in a quantitative way the validity of the analytic criteria (3.15) we focus on the evolution of the stable and unstable manifolds associated to the saddle fixed points. Parameters a=1, b=−1  and the frequency ω=1 in expressions (3.15) give a critical value for the regions A0 and A1: (4.1)Δ0=0.9115,Δ1=0.4530, or critical values of the coefficients of a damping force for ε=0.02: (4.2)δ0=εΔ0=0.01823,δ1=εΔ1=0.00906. Figure 12 shows numerical simulations of the phase space with initial conditions close to the undisturbed separatrix (θ0=−1.0572, θ˙0=0.01, ϕ0=π/10) for the region A0. Now, we reset the value of δ from δ0=0.01823 to greater ones (see Figure 12). It can be observed clearly that, for δ<δ0  (δ=0.018), the stable and unstable manifolds transversally intersect each other (Figure 12(a)). However, when δ>δ0  (δ=0.020), the invariant manifolds do not intersect (Figure 12(b)).Evolution of  stable and unstable manifolds as a function of the damping force coefficientsδ for a=1, b=−1, ε=0.02, ω=1 and two different values of δ close to  the critical value δ0=0.01823. (a) δ=0.018 (b) δ=0.020, with following initial conditions: θ0=−1.0572, θ˙0=0.01, ϕ0=π/10 for the area A0.
(a)(b)Similar results for the regionA1 are shown in Figure 13 (δ1=0.00906) for the following initial conditions: θ0=0.9472, θ˙0=0.2, ϕ0=π. Figure 13(a)  δ<δ1  (δ=0.0090).  Figure 13(b)  δ>δ1  (δ=0.0113). Thus the description, based on numerical simulations for some certain parameter values, makes a good match with the analytic criteria (3.15) provided by Melnikov method.Evolution of  stable and unstable manifolds as a function of the damping force coefficientsδ for a=1, b=−1, ε=0.02, ω=1 and two different values of δ close to the critical value δ1=0.00906. (a) δ=0.0090 (b) δ=0.0113, with following initial conditions: θ0=0.9472, θ˙0=0.2, ϕ0=π for the area A1.
(a)(b)


## 5. Conclusion

This work attempts to describe the transient cases occurring during a spacecraft descent in a planet atmosphere using methods of chaotic mechanics, in particular, the Melnikov method. We have suggested to introduce the concept of a biharmonic system (1.1) which reflects the behavior of the spacecrafts, and also, probably, more general mechanical systems. We have established the existence of transient heteroclinic chaos by means of the Melnikov method. Moreover, this method has provided an analytic criterion for the existence of chaotic behavior in terms of the system parameters. We have found a transition from chaotic to regular regime in the motion of the biharmonic oscillator, as the heteroclinic chaos can be removed by increasing the coefficient of a damping force. The analytic results given by the Melnikov method have been confirmed by a good match with the numeric research.



---

*Source: 319179-2009-03-26.xml*