Dynamical Tangles in Third-Order Oscillator with Single Jump Function

Research Article Dynamical Tangles in Third-Order Oscillator with Single Jump Function JilíPetrDela,1Tomas Gotthans,1and Milan Guzan2 1DepartmentofRadioElectronics,BrnoUniversityofTechnology,Purkynova118,61200Brno,CzechRepublic 2Department of Theoretical Electrotechnics and Electrical Measurement, Technical University of Kosice, Letna 9, 042 00 Kosice, Slovakia Correspondence should be addressed to Tomas Gotthans; go[email protected]tbr.cz Received 16 July 2014; Revised 18 September 2014; Accepted 18 September 2014; Published 3 December Academic Editor: Esteban Tlelo-Cuautle Copyright © 2014 Jiˇ r´ ıPetr ˇ zela et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. This contribution brings a deep and detailed study of the dynamical behavior associated with nonlinear oscillator described by a single third-order differential equation with scalar jump nonlinearity. The relative primitive geometry of the vector field allows making an exhaustive numerical analysis of its possible solutions, visualizations of the invariant manifolds, and basins of attraction as well as proving the existence of chaotic motion by using the concept of both Shilnikov theorems. The aim of this paper is also to complete, carry out and link the previous works on simple Newtonian dynamics, and answer the question how individual types of the phenomenon evolve with time via understandable notes. 1. Introduction It is well known that a majority of the real physical systems can be modeled by the system of the first-order differential equations with some sort of nonlinearity. In the case of the systems with at least three degrees of freedom the solution is notrestrictedtostableequilibriumorlimitcyclesbutthereis a certain chance to observe a much more complicated motion like chaos or hyperchaos [1]. This is a long-term unpredictable behavior caused by the so-called folding and stretching mechanism; first is responsible for solution bounded in finite state space volume and second for extreme sensitivity to the tiny changes of the initial conditions. Looking at this signal in time and frequency domain it resembles noise in many aspects. In reality, the individual waveforms combined together give rise to the strange attractors with fractal dimension [2] characterized by density, ergodicity, and mixing property. For chaotic attractor produced by third-order dynamical system value of geometrical dimension belongs to the range between two and three. Since chaos is a robust steady state dynamical motion, it should be somehow distinguished from chaotic transients [3]. The rigorous mathematical tool proving its existence can be picked as one of the famous Shilnikov theorems (ST) [4]. Roughly speaking, if there hold certain conditions for the eigenvalues and the strategic orbits associated with thesameequilibriumisdiscovered,theso-calledShilnikov’s chaos can be observed. As will be clarified later, additional informationmustbeobtainedbeforethestartofthesearching procedure, such as location of the fixed points, eigenspaces, boundary planes, attraction sets, and corresponding basins. The description of procedure solving this problem for the famous Chuas equations [5]canbefoundinpublication[6]. Many associated problems like vector field geometry of the so-called double-hook or dual double-scroll attractors [7] aresolvedintheinterestingbook[8]. Also chaos evolution principlesforsimpledrivensystemscanbefoundhere. This paper is organized as follows. The second section introduces the mathematical model of the nonlinear oscillator and brings its brief linear analysis. The third section focuses on linear topological conjugacy LTC [9] and presents equivalent dynamical systems. In other words the question if the mathematical model under inspection forms an entire class of the dynamical systems will be answered using similar approach as demonstrated in [10]. The fourth section exhibits one possible approach to find two mirrored homoclinic Hindawi Publishing Corporation e Scientific World Journal Volume 2014, Article ID 239407, 12 pages http://dx.doi.org/10.1155/2014/239407 2The Scientific World Journal orbits or heteroclinic connection between two fixed points. These trajectories are confirmed numerically together with associated chaotic behavior. The visualization of the basins of attraction and different manifolds is the core of the next section. Illustration of the structural stability of the chaoticattractors[11] by calculation of the largest Lyapunov exponents (LE) in the neighborhood of the nominal system parameters is a content of a next section. Such form of stability is essential from the viewpoint of physical construction of chaotic oscillator, for example, as electronic circuit. Since values of the circuit elements are functions of mathematical model parameters the sensitivity with respect to chaos deformation or destruction will be calculated. Finally concluding remarks, further research suggestions, and future topics are provided. 2. Mathematical Model Assume the following dynamical system [12,13] described by a single third-order differential equation, where the individual state variables can be interpreted as position, velocity, and acceleration, which belongs to the task from classical Newtonian dynamics: ... 𝑥+𝜑1𝑥+𝜑2𝑥=𝜉1𝑥+𝜉2sign (𝑥), 𝜑𝐴 1=0.6, 𝜑𝐴 2=1, 𝜑𝐵 1=0.7, 𝜑𝐵 2=0.8, (1) where dots represent derivatives with respect to the independent variable, time. There are, in fact, two different dynamical systems, one for each sign combination inside right-hand side function. Let these variants denote by symbols 𝐴,𝐵and let the nominal set of the parameters lead to the strange attractors: 𝜉𝐴 1=−1.2, 𝜉𝐴 2=2, 𝜉𝐵 1=1.2, 𝜉𝐵 2=−0.6. (2) Note that there is a significant degree of vector field symmetry. Both systems have simply defined boundary plane BP :{𝑥=0,𝑦∈R1,𝑧∈R1}separating the vector field into two outer affine segments and single virtual inner region. In further text, these segments will be denoted as 𝐷+,𝐷0,and 𝐷−.Therearejustthreeequilibriums,onepereachvector field segment, located at x𝑒0 =(000)𝑇,x±=(±𝜉2 𝜉100)𝑇.(3) It is evident that the characteristic polynomial and set of theeigenvaluesisthesameforbothfixedpoints.Tobemore specific it is cubic polynomial: 𝜆3+𝜑1𝜆2+𝜑2𝜆∓𝜉1=3 ∏ 𝑖=1 (𝜆−𝜆𝑖)=0. (4) For the system case A,aftersubstituting(2) in (1) we get apairofthecomplexconjugatedandasinglerealeigenvalue, namely, configuration R3∈𝐸𝑠 1⊕𝐸𝑢 2: 𝜆1,2 =𝜆󸀠±𝜆󸀠󸀠𝑗=0.162±1.128𝑗, 𝜆3=−0.924. (5) Thus fixed point is saddle-focus with stability index one. The vector field geometry resembles double scroll attractor generatedbywell-knownChua’sequationwithsuppressed inner segment [5]. The system case 𝐵for the substitution (3) in (1) has the reverse stability index two, in detail R3∈𝐸𝑢 1⊕𝐸𝑠 2: 𝜆1,2 =𝜆󸀠±𝜆󸀠󸀠𝑗=−0.594±1.160𝑗, 𝜆3=0.588. (6) The fixed point located at the origin cannot be treated as a regular one; it acts more likely as a virtual equilibrium. The inner segment 𝐷0is very narrow and cannot be easily observed even in the case of conventional numerical analysis. Following the rules of linear algebra we cannot derive additional useful information about dynamical system global motion. That is why further analysis is restricted to the existing numerical methods, that is, approaches which utilize the numerical integration process. 3. Linear Transform of Coordinates Suppose original dynamical system (1) xandassociatedsystem after linear transformation of the coordinates xwritten in the compact matrix form [14] x=Ax +bsign (w𝑇x), x=T−1ATx+T−1bsign (w𝑇Tx), (7) where Tis regular square matrix (3×3)oftherealnumbers. Note that such transforms can shift, rotate, or linearly stretch and compress the state space volume in some direction while leaving eigenvalues unchanged. New system can be advantageous from the viewpoint of circuitry implementation, symbolic analysis, understanding underlying dynamics, attractor visualization, and so forth. Since (1) expressed in terms of (7) is already in normal form A𝑅=(010 001 𝜉1−𝜑2−𝜑1), b𝑅=(00𝜉2), w𝑅=(001), (8) first example of LTC is the generation of state matrix A𝐽 in Jordan form [15].Thisisfundamentalconversionwith transformation matrix T𝐽with columns composed of the real and imaginary part of the complex eigenvector 𝜆𝑐and real eigenvector 𝜆𝑟: T𝐽=(Re (𝜆1 𝑐)Im (𝜆1 𝑐)𝜆1 𝑟 Re (𝜆2 𝑐)Im (𝜆2 𝑐)𝜆2 𝑟 Re (𝜆3 𝑐)Im (𝜆3 𝑐)𝜆3 𝑟), The Scientific World Journal 3 20 −2 −4 2 0 −2 −4 2 0 −2 −4 4 y x z (a) 1 0 −1 −5 −5 0 0 5 5 yx z (b) 44 0 0 0 2 22 −2 −2 −2 −4 −4 yx z (c) 4 0 02 2 −2 −2 −4 −4 yx z0 5 −5 (d) Figure 1: System configuration 𝐴in normal form (a), Jordan form (b), and first and second equivalent. T−1 𝐽A𝑅T𝐽=A𝐽=(𝜆󸀠𝜆󸀠󸀠 0 −𝜆󸀠󸀠 𝜆󸀠0 00𝜆 3), (9) where state matrix A𝑅directly represents system (1) and upper index of 𝜆𝑐and 𝜆𝑟denotes the corresponding element oftheeigenvector.Theknowledgeofsuchsystemcanbe useful for motion analysis and return map diagnosis since eigenspaces in each region of the state space are orthogonal. Using the concept of LTC the concrete form of the state matrix Aand vector wcanbeprescribedandtransformation between desired system and equivalent system in the normal formcanbeestablishedaccordinglyto[16]. The only restriction is that the resulting transformation must be regular square matrix. For the first equivalent system in the sense [17] the transformation T𝐼is A𝐼=(−𝜑1−1 0 𝜑20−1 −𝜉100 ), w𝐼=(100), T𝐼=(w𝑇 𝐼 w𝑇 𝐼A𝐼 w𝑇 𝐼A2 𝐼)=( 100 −𝜑1−10 𝜑2 1−𝜑2𝜑11). (10) 4The Scientific World Journal y 0 0 −1 1 0.5 0.5 0 −0.5 −0.5 −1 x z (a) y 0 0 1 1 0 −0.5 −1 −1 x z (b) y 1 1 1 2 2 0 0 0 −1 −1 −1 −2 −2 x z (c) y1 0.5 0.5 0.5 0 00 −0.5 −0.5 −0.5 x z (d) Figure 2: System configuration 𝐵in normal form (a), Jordan form (b), and first and second equivalent. Analogically for the second equivalent system the matrix T𝐼𝐼 becomes A𝐼𝐼 =(−𝜑1𝜑2−𝜉1 −1 0 0 0−10), w𝐼𝐼 =(100), T𝐼𝐼 =(w𝑇 𝐼𝐼 w𝑇 𝐼𝐼A𝐼𝐼 w𝑇 𝐼𝐼A2 𝐼𝐼)=( 100 −𝜑1𝜑2−𝜉1 𝜑2 1−𝜑2−𝜑1𝜑2+𝜉1𝜑1𝜉1).(11) The three-dimensional perspective views on the chaotic state space attractors together with plane projections for each dynamical system mentioned above and obtained by using Mathcad with build-in fourth-order Runge-Kutta numerical integration method are shown in Figures 1and 2,respectively. For these simulations final time equals 𝑡max =1000with time step 𝑡step =0.01. The initial conditions were set: x0𝐴 = (0.1,0,0)𝑇and x0𝐵 =(0.1,0,0)𝑇. Note that LTC operation is demonstrated by means of Figure 3. 4. Strategic Trajectories Before starting with description of the procedure for finding some strategic orbit in the sense of Shilnikov, the results (5) The Scientific World Journal 5 Figure 3: Geometrical visualization of LTC between equivalent systems; see text. and (6) prompted that both essential conditions desired by ST are satisfied simultaneously for system of classes 𝐴and 𝐵: 𝜆󸀠𝜆3<0∧󵄨󵄨󵄨󵄨𝜆3󵄨󵄨󵄨󵄨>󵄨󵄨󵄨󵄨󵄨𝜆󸀠󵄨󵄨󵄨󵄨󵄨.(12) Remember that this condition itself does not guarantee the presence of chaotic behavior; ST requires also strategic orbit associated with some fixed point. The problem with searching for specific state trajectory can be effectively converted into optimization task [18]. The basic form of the vector field with de facto two linear segments significantly simplifies the algebraic set of the fitness functions for optimization. The process of derivation of such set is described in step-by-step manner in paper [19]. The optimization routine seeks through parameter space, alters the eigenvalues, and rotates the corresponding eigenspaces simultaneously. The goal function value is minimized and penalized if geometry of the vector field or desired property of the system changes. It is preserved by choosing the suitable guess values as well as by the restrictions on the parameter space under inspection. The geometric structures like points, lines, and planes important for optimization are defined in Figure 4. 4.1. Homoclinic Tangle. By definition, the homoclinic orbit is forwards and backwards asymptotic to the same equilibrium point. Due to the vector filed symmetry there is always a pair of the homoclinic tangles. Therefore we can focus on homoclinic orbit associated with fixed point x−located in segment 𝐷−. Following full integration method described in [20]tends to be very time consuming and, as a part of an optimization routine, it can diverge even in the case if homoclinic orbit exists. Our problem should be addressed separately for dynamical system cases 𝐴and 𝐵. In case 𝐵let further investigation be focused on homoclic orbit associated with fixed point x−.Sincewehaveunstable eigenvector we can start with numerical integration in the intersection of this eigenvector and boundary plane. Dynamical flow became a part of optimization procedure which will bestoppedassoonastrajectoryleavessegment𝐷+.This corresponds to the unique mapping which should be satisfied: (0,𝜆3𝜉2 𝜉1,𝜆2 3𝜉2 𝜉1)𝑇Φ𝑡>0 󳨀󳨀󳨀→(0,𝛼,(𝜉2 𝜉1−𝜔1𝛼−(𝜉2/𝜉1)(𝜎2/𝜎1) 𝜔2−𝜔1(𝜎2/𝜎1))𝜎3 𝜎1+𝜔3𝛼−(𝜉2/𝜉1)(𝜎2/𝜎1) 𝜔2−𝜔1(𝜎2/𝜎1))𝑇,(13) where 𝛼is arbitrary value. In the case 𝐴, the situation is analogical, except that backward integration instead of standard is performed. To get homoclinic connection there must exist a mapping (13) but for times 𝑡<0. Final trajectory will consist of three pieces; free-motion in one segment of the vector field and forced motions along unstable eigenspace (forward integration) and stable eigenspace (backward integration) in the other segment. Utilizing this concept the new set of the parameters for system case 𝐴has been found as 𝜑1=0.63339561,𝜑2=1.00676348, 𝜉1=−1.29333691,𝜉2=1.98006640,initialconditionsx0𝐴 = (−𝜉2/𝜉1,0,0)𝑇+(1⋅10−10,0,0)𝑇, fourth order Runge-Kutta build-in MATLAB function integration with variable step (initial step ℎ=1⋅10−6) and with maximal step ℎ=1⋅10−2, and final value of fitness function 1.3⋅10−3. Similarly for dynamical system case B we get 𝜑1= 0.715499177650,𝜑2= 0.577154986868,𝜉1= 1.284935534298,𝜉2=−0.631662627700,initialconditions x0𝐵 =(𝜉2/𝜉1,0,0)𝑇−(1⋅10−8,0,0)𝑇, fourth order RungeKutta build-in MATLAB function integration with variable step (initial step ℎ=1⋅10−6) and with maximal step ℎ=1⋅10−2,andfinalvalueoffitnessfunction7.62⋅10−4.The numerically integrated trajectories are provided in Figure 5. A careful reader can raise an objection that integration as a part of optimization can be removed by solving a system of the linear differential equations. This is possible although the resulting analytical formulas are quite complicated. By considering known initial conditions and substituting desired line-type intersection it is possible to separate internal system parameters 𝜑1,𝜑2,𝜉1,and𝜉2.Thusthisapproachalsocannot provide a time required for trajectory to leave affine segment under inspection. 4.2. Heteroclinic Tangle. The heteroclinic orbit can be considered as a generalization of a saddle loop. In fact, searching for the heteroclinic connection is a two-dimensional problem. 6The Scientific World Journal (a) (b) Figure 4: The geometric structures important for optimization, starting situation, and systems 𝐴(a) and 𝐵(b), where boundary planes are white, fixed points are black dots, eigenvectors are green, eigenplanes are gray, intersections of eigenvectors and boundary planes are red crosses, and intersections of eigenplanes and boundary plane are blue lines. (a) (b) Figure 5: Trajectory close to the homoclinic orbit for dynamical system case 𝐴(a) and case 𝐵(b). (a) (b) Figure 6: Trajectory close to the heteroclinic orbit for dynamical system cases 𝐴(a) and 𝐵(b). Since objective functions are known analytically MATLAB and build-in gradient-based technique has been utilized. The fitness function for this geometric structure is minimization of distance between point 𝑋and line 𝑌such that 𝑧−(𝜉2 𝜉1−𝜔1𝑦−(𝜉2/𝜉1)(𝜎2/𝜎1) 𝜔2−𝜔1(𝜎2/𝜎1))𝜎3 𝜎1 −𝜔3𝑦−(𝜉2/𝜉1)(𝜎2/𝜎1) 𝜔2−𝜔1(𝜎2/𝜎1)=0, (14) where 𝜔is an imaginary part of the complex eigenvector and 𝜎represents its real part and auxiliary constants. Although theanalyticsolutioncanbederivedthetimenecessaryforthis unique mapping is unknown. Thus a numerical integration process should be a part of the optimization with the initial conditions set in vicinity of some fixed point. For the first sign case of the differential equations (1) one can found the following set of the parameters: 𝜑1=−0.6261, 𝜑2=0.9536,𝜉1=−1.2655,and𝜉2=1.9937. Analogically for case 𝐵sign variant we get 𝜑1=−0.6378,𝜑2=0.6011, 𝜉1=1.3156,and𝜉2=−2.1697. The strategic orbits are unstable such that for numerical integration the initial conditions must be chosen carefully. For the strategic orbits, the process of numerical integration mustbegininthecloseneighborhoodofthefixedpoint. Using time-forward integration the state point goes into the The Scientific World Journal 7 10 10 0 0 −10 10 0 −10 −10 10 10 0 0 −10 10 0 −10 10 0 −10 10 0 −10 10 0 −10 10 0 −10 10 0 −10 −10 10 0−10 10 0−10 10 0−10 10 0−10 10 0−10 10 0−10 10 10 0 0 −10 −10 10 10 0 0 −10 −10 10 10 0 0 −10 10 0 −10 10 0 −10 −10 10 0−10 10 0−10 10 0−10 10 0 −10 10 0 −10 10 0−10 10 0 −10 y = 0.5 y = 5 y = 5.5 y = 6.5 y = 8 y=2 y=3 y=4 y=−4 y=−3 y=−2 y = −0.5 y = −8.0 y = −6.5 y = −5.5 y = −5 Figure 7: The graphical illustration of the basin of attraction for chaotic attractor, system case 𝐴. These slices are 𝑥-𝑧projections of crosssections with the 𝑦-axis. Black region represents initial condition inside basin of attraction. All other trajectories tend to infinity. directionofunstablemanifoldtowardstheoppositeequilibrium. The saddle loop is finished by using time-backward numerical integration. The resulting state trajectories for both sets of the parameters are given in Figure 6.Despitebeing structurally unstable, these trajectories can be constructed and destructed via a manipulation with vector field geometry, that is, by changing the internal system parameters. Recently it has been verified that inverse approach can be used; starting with the satisfaction of one ST the original mathematical model can be derived [21]. 5. Overall Numerical Analysis Although studied dynamical system is algebraically quite simple, it is nonlinear and chaotic. Therefore there is no closed-form analytic solution. Thus the analysis is restricted to the numerical procedures mostly based on integration of the state space trajectory. This process is also involved in the routine for calculation of the spectrum of the LE [22]. These real numbers measure the average ratio of exponential separation between trajectories in the state space. For chaotic 8The Scientific World Journal 0 0 7 7 −7 −7 0 0 7 7 −7 0 7 −7 0 7 −7 0 7 −7 0 7 −7 0 7 −7 0 7 −7 −7 07−7 07−7 07−7 07−7 0 0 7 7 −7 −7 07 −7 07−7 0 0 7 7 −7 0 7 −7 0 7 −7 −7 07−7 07−7 y = 0.7 y=1 y = 1.4 y=2 y = 0.1 y = 0.4 y=−2 y = −0.4 y = −0.1 y=−1.4 y=−1 y = −0.7 Figure 8: The graphical illustration of the basin of attraction for chaotic attractor, system case 𝐵. These slices are cross-sections of the 𝑦-axis where 𝑦={−2,−1.4,−1,−0.7,−0.4,−0.1,0.1,0.4,0.7,1,1.4,2}(left to right). Black and green regions mark basins of attraction for case 𝐴and case 𝐵, respectively. Blue and red regions denote unbounded solutions. motion it is necessary to have one positive LE; the sum of all LE must be negative since the dynamics is dissipative due to the parameter 𝜑1>0. This concept has been usedtoprovethattheregionforchaosiswideenough to preserve some structural stability [23] of the desired strange attractor; such form of stability allows developing the electronic circuits with the same dynamics [24]. The routine for LE spectrum computation can be effectively utilized also forthepurposeofdetectingchaosinthegeneralclassof arbitrary-order dynamical systems; for further details see [25]. Due to discontinuity in the system the derivations in Jacobi matrix have to be threated in order to determine Lyapunov exponents correctly. We obtained more precise results by using sigmoid function with extreme value of 𝛾 as defined in [26]. One other possible demonstration dealing withdiscontinuityintheforcedsystemcanbefoundin[27]. The remaining question which should be answered is the following: how many attractors are available for the nominal set of the parameters and show the associated basins of attractions. The most straightforward approach to visualize these subspaces is by means of repeated integrations. The grid for the initial condition was a cube with edge lengths 𝛿 ∈ (−2,2)with 400points. Due to the symmetry of thevectorfieldtwomirroredstrangeattractorsarehighly expected. It eventually turns out that the system case 𝐴has only unbounded solution or chaotic attractor; trivial fixed point solution is out of question. These results are visible in Figure 7.Forthesystemcase𝐵two chaotic attractors have been found; see Figure 8. The graphical illustration of the sensitivity to initial conditionscanbeseeninFigure 9.Thetotalnumberof random generations is 10⋅103with standard deviation 𝜎= 0.01around initial point 𝑥0𝐴𝐵 =(0,0,0)𝑇,totaltime𝑇end = 100, and integration step ℎ=0.01. 6. Circuitry Implementation In order to evaluate geometrical structural stability of equationsanewcircuit(notpresentedsofarfor(1)) was assembled and measured. The circuit synthesis methods dedicated to modeling the nonlinear dynamical systems are well-known The Scientific World Journal 9 (a) (b) Figure 9: The graphical illustration of the sensitivity to initial conditions. Where the total number of random generations is 10⋅103with standard deviation 𝜎=0.01around initial point 𝑥0𝐴𝐵 =(0,0,0)𝑇,totaltime𝑇end =100, and integration step ℎ=0.01. On the left there is caseAandontherightthereiscase𝐵. The green points represent initial conditions, red color is a final point, and gray is original attractor. C1 C2 C3+ + − − Ra Rc R2 R1 Figure 10: The circuitry implementation of chaotic oscillator for case 𝐴. andcommonlyused[28,29]. Assume canonical (in the sense of minimum circuit components) network shown in Figure 10 which represents a parallel connection of the thirdorder linear admittance and two-segment piecewise-linear resistor. The straightforward analysis gives us admittance function in Laplace transform, namely, 𝑌(𝑠)=𝑠3+𝜑1𝑠2𝜑2𝑠 ≈𝑐1𝑐2𝑐3𝑟1𝑟2𝑠3+𝑐1𝑐2(𝑟1+𝑟2)𝑠2+(𝑐1+𝑐2)𝑠, (15) whereresistorsandcapacitorshavenormalizedvaluesso far. To obtain the real passive element values let consider impedance normalization factor 𝜉𝐼and frequency norm 𝜉𝐹. By comparing the individual coefficients (15) with (1) we get following simple relationships 𝐶1=𝐶2=𝜑2 2𝜉𝐼𝜉𝐹,𝐶 3=𝜑2 2 𝜑2 1𝜉𝐼𝜉𝐹,𝑅 1=𝑅2=2𝜑1𝜉𝐼 𝜑2 2. (16) The amper-voltage characteristics of nonlinear resistor can be defined by two equations depending on the sign of input voltage, in detail: 𝐼=𝑓(𝑉)=−(1𝑟𝑎+1𝑟𝑐)𝑉+𝑉Sat 𝑟𝑎sign (𝑉),(17) where 𝑉sat ≈13Vrepresentsasaturationvoltageofthe used operational amplifier TL082and orientation of current