Dynamical Tangles in Third-Order Oscillator with Single Jump Function

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.


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 not restricted to stable equilibrium or limit cycles but there is 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 the same equilibrium is discovered, the so-called Shilnikov's chaos can be observed. As will be clarified later, additional information must be obtained before the start of the searching 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] can be found in publication [6]. Many associated problems like vector field geometry of the so-called double-hook or dual double-scroll attractors [7] are solved in the interesting book [8]. Also chaos evolution principles for simple driven systems can be found here.
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 2 The 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 chaotic attractors [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.

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 + 2 sign ( ) , 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: Note that there is a significant degree of vector field symmetry. Both systems have simply defined boundary plane BP : { = 0, ∈ R 1 , ∈ R 1 } 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 − . There are just three equilibriums, one per each vector field segment, located at It is evident that the characteristic polynomial and set of the eigenvalues is the same for both fixed points. To be more specific it is cubic polynomial: For the system case A, after substituting (2) in (1) we get a pair of the complex conjugated and a single real eigenvalue, namely, configuration R 3 ∈ 1 ⊕ 2 : Thus fixed point is saddle-focus with stability index one. The vector field geometry resembles double scroll attractor generated by well-known Chua's equation with suppressed inner segment [5].
The system case for the substitution (3) in (1) has the reverse stability index two, in detail R 3 ∈ 1 ⊕ 2 : 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 0 is 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.

Linear Transform of Coordinates
Suppose original dynamical system (1)ẋ and associated system after linear transformation of the coordinatesẋ written in the compact matrix form [14] x = Ax + b sign (w x) , where T is regular square matrix (3 × 3) of the real numbers. 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 first example of LTC is the generation of state matrix A in Jordan form [15]. This is fundamental conversion with transformation matrix T with columns composed of the real and imaginary part of the complex eigenvector and real eigenvector : The Scientific World Journal where state matrix A directly represents system (1) and upper index of and denotes the corresponding element of the eigenvector. The knowledge of such system can be 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 A and vector w can be prescribed and transformation between desired system and equivalent system in the normal form can be established accordingly to [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 Analogically for the second equivalent system the matrix T becomes 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 1 and 2, respectively. For these simulations final time equals max = 1000 with time step step = 0.01. The initial conditions were set: x 0 = (0.1, 0, 0) and x 0 = (0.1, 0, 0) . Note that LTC operation is demonstrated by means of Figure 3

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 and (6) prompted that both essential conditions desired by ST are satisfied simultaneously for system of classes and : 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.

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 − . Since we have unstable 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 be stopped as soon as trajectory leaves segment + . This corresponds to the unique mapping which should be satisfied: 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,  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 . Thus this approach also cannot provide a time required for trajectory to leave affine segment under inspection.

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. 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 where is an imaginary part of the complex eigenvector and represents its real part and auxiliary constants. Although the analytic solution can be derived the time necessary for this 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, direction of unstable manifold towards the opposite equilibrium. 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. Despite being 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].

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 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 used to prove that the region for chaos is wide enough 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 for the purpose of detecting chaos in the general class of 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 with discontinuity in the forced system can be found in [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 400 points. Due to the symmetry of the vector field two mirrored strange attractors are highly 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. For the system case two chaotic attractors have been found; see Figure 8.
The graphical illustration of the sensitivity to initial conditions can be seen in Figure 9. The total number of random generations is 10 ⋅ 10 3 with standard deviation = 0.01 around initial point 0 = (0, 0, 0) , total time end = 100, and integration step ℎ = 0.01.

Circuitry Implementation
In order to evaluate geometrical structural stability of equations a new circuit (not presented so far for (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 ⋅ 10 3 with standard deviation = 0.01 around initial point 0 = (0, 0, 0) , total time end = 100, and integration step ℎ = 0.01. On the left there is case A and on the right there is case . The green points represent initial conditions, red color is a final point, and gray is original attractor. and commonly used [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, where resistors and capacitors have normalized values so 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 The amper-voltage characteristics of nonlinear resistor can be defined by two equations depending on the sign of input voltage, in detail: where sat ≈ 13V represents a saturation voltage of the used operational amplifier TL082 and orientation of current  Figure 11: Several examples of the laboratory measurements (for case ): double-scroll attractor (a,b), single-scroll (c), AV curve of nonlinear resistor including saturation (d), limit cycle (e), and period doubling (f).
is outwards the resistor. In fact note that due to saturation in practice the function has four segments. By adopting the numerical constants (2) and norms = 10 4 and = 10 4 we get the final values of the passive elements to generate doublescroll attractor: = 65 kΩ, = 9.6 kΩ.
The experimental verification (measurements) of this chaotic oscillator is provided by means of Figure 11. Note that proposed circuitry represents only case dynamical system (2). Variant can be modeled by using slightly modified    Figure 12 give an idea if uncertainties in values of the passive components are capable to destruct chaotic nature of circuit. It is obvious that 1% tolerances cannot cause such changes, but 10% tolerances can lead to death of doublescroll in about 8% of cases. The total number of turns for each histogram is 10000 (for case ) with initial conditions i = (0.1, 0, 0) .

Conclusion
This contribution addresses problems associated with complete analysis of the third-order dynamical system with single jump-type nonlinearity. The dynamical system under inspection belongs to the generalized class of the Chuas equations published in [26]. The existence of the Shilnikov chaos has been proved by using optimization and numerical integration. The specialty of given dynamical system is in changing the sign of nonlinear function; the homoclinic orbits or heteroclinic orbits appear. That has not been demonstrated yet. However the electronic circuit capable to model the nonlinear behavior near the strategic orbits has not been presented so far.
There are many unsolved topics involving dynamical systems with possible chaotic solution. A fraction of them are mentioned in the paper [30]; but personal experience and attempts to solve real specific system necessarily unfold many others. The universality of chaotic behavior is uncovered in