# Energy Conservation, Singular Orbits, and FPGA Implementation of Two New Hamiltonian Chaotic Systems

Authors: Enzeng Dong; Guanghan Liu; Zenghui Wang; Zengqiang Chen
Journal: Complexity
 (2020)

---

## Abstract

Since the conservative chaotic system (CCS) has no general attractors, conservative chaotic flows are more suitable for the chaos-based secure communication than the chaotic attractors. In this paper, two Hamiltonian conservative chaotic systems (HCCSs) are constructed based on the 4D Euler equations and a proposed construction method. These two new systems are investigated by equilibrium points, dynamical evolution map, Hamilton energy, and Casimir energy. They look similar, but it is found that one can be explained using Casimir power and another cannot be explained in terms of the mechanism of chaos. Furthermore, a pseudorandom signal generator is developed based on these proposed systems, which are tested based on NIST tests and implemented by using the field programmable gate array (FPGA) technique.


---

## Body

## 1. Introduction

Many dissipative chaotic systems have been widely studied in the past decades. Lorenz discovered the first chaos system named Lorenz system in 1963 [1]. In 1999, Chen and Ueta proposed a new chaotic attractor called Chen attractor, which is topologically more complex than Lorenz attractor [2]. Lü and Chen developed a new chaotic system named Lü system in 2002, which bridged the gap between the Lorenz system and the Chen system [3]. The complex dynamical behaviors of chaotic systems have been gradually excavated. Both conservative systems [4, 5] and dissipative systems [6–8] have some phenomena of transition from period to quasi-period and then into chaos. For the past few years, due to the potential applications of chaotic systems with complex dynamic behavior in speech signals and digital image/video encryption, analysis in this area has also become one of the main directions of chaos research [9, 10]. Dong et al. rigorously verified the chaotic properties of a 3D autonomous four-wing chaotic attractor using topological horseshoe theory and numerical calculations [11, 12]. As we know, the dissipative chaotic system (DCS) has strange attractors, which can be reconstructed by using the delay embedding method [13]. Compared with DCS, CCS has no general attractors, and the conservative chaotic flows have pseudorandomness and white-noise-like properties [14, 15]. Recently, there are some conservative chaotic systems that have been constructed and analyzed. Cang proposed several conservative chaotic systems with only linear and quadratic terms [16]. Soon afterwards, a new method for constructing high-dimensional conservative chaotic systems was proposed by Dong and Yuan, and a four-dimensional conservative chaotic system and a 5-dimensional conservative hyperchaotic system were specifically investigated [17].Existing conservative systems can be divided into two categories: Hamiltonian conservative chaotic system (HCCS) and non-Hamiltonian CCS (non-HCCS). The HCCSs contain two subclasses: traditional HCCS [18] and generalized HCCS. The typical examples of traditional HCCSs are He’non-Heiles system (the maximum Lyapunov exponent LE1=0.02) [19] and quantum chaotic systems, and the classic generalized HCCSs are Cang-case A system LE1=1.2 [20] and Qi system with large LEs LE1=78194.27 [21]. HCCS has zero-sum LEs and satisfies volume conservation ∇·f=0 and energy saving H˙=0. However, non-HCCSs only meet the zero-sum LEs and do not need to satisfy the conditions of ∇·f=0 or H˙=0. For example, Sprott A system LE1=0.014 has divergence ∇·f=x3x3∈R [22], and Vaidyanathan and Volos system LE1=0.0395 has ∇·f=2x3 [23] and the 6D conservative hyperchaotic system LE1=LE2=0.055 has ∇·f=2x5 [24]. As can be seen from the above classification and comparison, Hamiltonian conservative systems need to meet stricter condition than other conservative systems.For the mechanism of chaos generation, Arnor’d first proposed a Kolmogorov system of using ordinary differential equations to describe a dissipative-forced dynamic system, and this system includes three field components: inertial, dissipative, and external components [25]. Afterwards, Qi transformed Qi four-wing chaotic system into a Kolmogorov-type system and then analyzed the cycling of energy among potential energy, kinetic energy, dissipation, and external energy [26]. Qi and Zhang and Wiggins [27, 28] found that the conservative part of the 3-dimensional Euler equations is symplectic structure in the generalized Hamiltonian system and the Kolmogorov system. Qi proposed six complete 4D Euler equations by extending a system from 3D to 4D [21].However, how to construct a system conforming to the four-dimensional Euler equations has not been fully studied. Currently, some conservative chaotic systems are verified by analog circuits [29], while the digital circuit implementation of Hamiltonian conservative chaos is rarely reported. In this paper, two Hamiltonian conservative chaotic systems are constructed. We first select three linear functions and construct a four-dimensional system using the methods of Reference [17]. It accords with the conservative four-dimensional Euler equations of Hamilton and Casimir function, but the four-dimensional system does not produce chaos. Then, two linear terms are added based on the modified 4D Euler equations and a new constructive method. Choosing proper parameters and initial values, the constructed two systems show chaotic dynamics. For high-dimensional systems, the equilibrium point analysis [30, 31] and bifurcation analysis [32] are necessary. In general, the noncanonical Hamiltonian systems have Casimir functions [33], but the structure matrix of the new two systems are different: one is noncanonical and another is canonical. Hence, we use two methods to analyze the mechanism of chaos generation. Pseudorandomness is confirmed following NIST statistic tests. In addition, the field programmable gate array technique is used to design pseudorandom number generators.The rest of this article is organized as follows. In Section2, the construction of two new 4D HCCSs is given. In Section 3, these two systems are analyzed using different theories. Section 4 contains the verification of strong pseudorandomness of these systems and a field programmable gate array-based digital pseudorandom number generator using the proposed HCCSs. Finally, the conclusions are drawn in Section 5.


## 2. Modeling HCCS

### 2.1. Preliminaries

Euler equation is the most important fundamental equation in the nonviscous fluid dynamics, and it is also very effective for describing rigid bodies with symmetric rotation [34].Σ
24 are certain 4D Euler rigid body equations, which can be described as the form of a Hamiltonian vector field:(1)x˙=J24x×∇Hx,where(2)J24x=0−x40x2x40−x4x3−x10x40−x2−x2x1−x3x20,and ∇ is the gradient operator. Hx is the Hamiltonian energy function.(3)Hx=12Π1x12+Π2x22+Π3x32+Π4x42,where Πi=Ii−1, Ii is the moment of inertia of the rigid body, and xi is the angular momentum. Then, system Σ24 can be obtained as(4)x˙1=Π4−Π2x2x4,x˙2=Π1−Π4x1x4+Π4−Π3x3x4,x˙3=Π2−Π4x2x4,x˙4=Π2−Π1x1x2+Π3−Π2x2x3.Remark 1.
Because the matrixJ24 satisfies the skew-symmetric structure, system Σ24 is conservative both in the Hamiltonian and Casimir energy. The Casimir function, also known as energy momentum, is an important physical quantity in fluid dynamics. Its derivative is Casimir power, and in some cases, the Casimir power can be used as a criterion for determining whether a system can produce chaos. The system is conservative when C˙≡0 with C≠0. The Casimir function C is properly selected as Cx=1/2x12+x22+x32+x42, which satisfies the Lie–Poisson brackets [35], C,H=0,∀H∈C∞g∗.
According to the method from Reference [17], select functions f3y,z=hy,f4x,w=aw,f5x,z=ex and then obtain the constructed system(5)x˙=hyw,y˙=azw+exw,z˙=−ayw,w˙=−h+exy.
SetΠ4−Π2=h, Π4−Π3=a, Π1−Π4=e, Π2−Π1=−h+e, and Π3−Π2=0; equation (5) can be transformed into the form of Euler equations (4). Because the constructed system (5) lacks a constant function, it cannot produce chaos, which corresponds to the Hamilton energy and Casimir energy conservation of four-dimensional Euler equations.


### 2.2. Case 1

The proposed system, equation (5), satisfies the four-dimensional Euler equation form. However, from Remark 1, this type of 4D Euler equations is conservative for Hamiltonian and Casimir energy, and it cannot generate chaos. So, we need to break the conservation. In J24Hx∇Hx, Hx is Hamiltonian, and the vector field form is immutable. However, the matrix J24Hx is symplectic and can be transformed into skew-symmetric form. So, J24H is obtained as(6)J24Hx=0−x4cx2x40−x4x3−x1−cx40−x2−x2x1−x3x20,where c is constant. And system Σ24H is described as(7)Σ24Hx˙1=Π4−Π2x2x4+cΠ3x3,x˙2=Π1−Π4x1x4+Π4−Π3x3x4,x˙3=Π2−Π4x2x4−cΠ1x1,x˙4=Π2−Π1x1x2+Π3−Π2x2x3,,where Σ24H represents the system that breaks the symplectic structure in Σ24. System Σ24H has 5 adjustable parameters: Πi,i=1,2,3,4, and c, which can influence the orbital modes of system Σ24H, and they have physical meanings. The initial value, xi0, determines the total Hamiltonian; when the initial Hamiltonian is fixed, Πi is proportional to the frequency of the rigid body. For system Σ24H,(8)H˙x=∇HxTJ24Hx∇Hx=0,and(9)C˙=xT·x.=xT·J24Hx·∇Hx=cΠ3−Π1x1x3.So, the Hamiltonian energy of systemΣ24H is conservative and the Casimir energy is nonconservative. In order to facilitate subsequent analysis, the parameter c is reserved. Since the Euler equations satisfy that the sum of the coefficients of the coupling term is zero, it is calculated that h,a, and e have a fixed relationship, that is, a=h,e=−0.5h. And set Π3=d,Π1=b; then, the modified system Σ24H can be transformed into our proposed equation form:(10)x˙=hyw+c·dz,y˙=hzw−0.5hxw,z˙=−hyw−c·bx,w˙=−0.5hxy.Whenh,d,b,c=−1,2.5,2,0.2 and the initial states x0,y0,z0,w0=0.5,0.5,1,1, system (10) is chaotic since its LEs are 0.2405,0,0,−0.2405.


### 2.3. Case 2

Selecting functionf2=g, where g is a constant, the constructed system is(11)x˙=hyw+gz,y˙=azw+exw,z˙=−ayw−gx,w˙=−h+exy.Seth,g,a,e=0.5,10,6,5.5 and x0,y0,z0,w0=1,1,1,1; then, system (11) is chaotic since its LEs are 1.0739,0,0,−1.0739.Remark2 illustrates that the Case 1 and Case 2 systems are fundamentally different.Remark 2.
Case 2 system can only be expressed as(12)x˙1x˙2x˙3x˙4=0−x40x2x40−x4x3−x10x40−x2−x2x1−x3x20Π1x1Π2x2Π3x3Π4x4+gx10gx30.
This is becauseΠ3≠Π1; if we set g=c·Π3, it is clear that g cannot be equal to c·Π1. But this form cannot get the Hamiltonian energy of Case 2 system. It should be expressed as follows:(13)x˙1x˙2x˙3x˙4=00ghx200ax4ex1−g−ax400−hx2−ex100x1x2x3x4.
So, the Hamiltonian energy of Case 2 system isHx=1/2x12+x22+x32+x42.


## 2.1. Preliminaries

Euler equation is the most important fundamental equation in the nonviscous fluid dynamics, and it is also very effective for describing rigid bodies with symmetric rotation [34].Σ
24 are certain 4D Euler rigid body equations, which can be described as the form of a Hamiltonian vector field:(1)x˙=J24x×∇Hx,where(2)J24x=0−x40x2x40−x4x3−x10x40−x2−x2x1−x3x20,and ∇ is the gradient operator. Hx is the Hamiltonian energy function.(3)Hx=12Π1x12+Π2x22+Π3x32+Π4x42,where Πi=Ii−1, Ii is the moment of inertia of the rigid body, and xi is the angular momentum. Then, system Σ24 can be obtained as(4)x˙1=Π4−Π2x2x4,x˙2=Π1−Π4x1x4+Π4−Π3x3x4,x˙3=Π2−Π4x2x4,x˙4=Π2−Π1x1x2+Π3−Π2x2x3.Remark 1.
Because the matrixJ24 satisfies the skew-symmetric structure, system Σ24 is conservative both in the Hamiltonian and Casimir energy. The Casimir function, also known as energy momentum, is an important physical quantity in fluid dynamics. Its derivative is Casimir power, and in some cases, the Casimir power can be used as a criterion for determining whether a system can produce chaos. The system is conservative when C˙≡0 with C≠0. The Casimir function C is properly selected as Cx=1/2x12+x22+x32+x42, which satisfies the Lie–Poisson brackets [35], C,H=0,∀H∈C∞g∗.
According to the method from Reference [17], select functions f3y,z=hy,f4x,w=aw,f5x,z=ex and then obtain the constructed system(5)x˙=hyw,y˙=azw+exw,z˙=−ayw,w˙=−h+exy.
SetΠ4−Π2=h, Π4−Π3=a, Π1−Π4=e, Π2−Π1=−h+e, and Π3−Π2=0; equation (5) can be transformed into the form of Euler equations (4). Because the constructed system (5) lacks a constant function, it cannot produce chaos, which corresponds to the Hamilton energy and Casimir energy conservation of four-dimensional Euler equations.


## 2.2. Case 1

The proposed system, equation (5), satisfies the four-dimensional Euler equation form. However, from Remark 1, this type of 4D Euler equations is conservative for Hamiltonian and Casimir energy, and it cannot generate chaos. So, we need to break the conservation. In J24Hx∇Hx, Hx is Hamiltonian, and the vector field form is immutable. However, the matrix J24Hx is symplectic and can be transformed into skew-symmetric form. So, J24H is obtained as(6)J24Hx=0−x4cx2x40−x4x3−x1−cx40−x2−x2x1−x3x20,where c is constant. And system Σ24H is described as(7)Σ24Hx˙1=Π4−Π2x2x4+cΠ3x3,x˙2=Π1−Π4x1x4+Π4−Π3x3x4,x˙3=Π2−Π4x2x4−cΠ1x1,x˙4=Π2−Π1x1x2+Π3−Π2x2x3,,where Σ24H represents the system that breaks the symplectic structure in Σ24. System Σ24H has 5 adjustable parameters: Πi,i=1,2,3,4, and c, which can influence the orbital modes of system Σ24H, and they have physical meanings. The initial value, xi0, determines the total Hamiltonian; when the initial Hamiltonian is fixed, Πi is proportional to the frequency of the rigid body. For system Σ24H,(8)H˙x=∇HxTJ24Hx∇Hx=0,and(9)C˙=xT·x.=xT·J24Hx·∇Hx=cΠ3−Π1x1x3.So, the Hamiltonian energy of systemΣ24H is conservative and the Casimir energy is nonconservative. In order to facilitate subsequent analysis, the parameter c is reserved. Since the Euler equations satisfy that the sum of the coefficients of the coupling term is zero, it is calculated that h,a, and e have a fixed relationship, that is, a=h,e=−0.5h. And set Π3=d,Π1=b; then, the modified system Σ24H can be transformed into our proposed equation form:(10)x˙=hyw+c·dz,y˙=hzw−0.5hxw,z˙=−hyw−c·bx,w˙=−0.5hxy.Whenh,d,b,c=−1,2.5,2,0.2 and the initial states x0,y0,z0,w0=0.5,0.5,1,1, system (10) is chaotic since its LEs are 0.2405,0,0,−0.2405.


## 2.3. Case 2

Selecting functionf2=g, where g is a constant, the constructed system is(11)x˙=hyw+gz,y˙=azw+exw,z˙=−ayw−gx,w˙=−h+exy.Seth,g,a,e=0.5,10,6,5.5 and x0,y0,z0,w0=1,1,1,1; then, system (11) is chaotic since its LEs are 1.0739,0,0,−1.0739.Remark2 illustrates that the Case 1 and Case 2 systems are fundamentally different.Remark 2.
Case 2 system can only be expressed as(12)x˙1x˙2x˙3x˙4=0−x40x2x40−x4x3−x10x40−x2−x2x1−x3x20Π1x1Π2x2Π3x3Π4x4+gx10gx30.
This is becauseΠ3≠Π1; if we set g=c·Π3, it is clear that g cannot be equal to c·Π1. But this form cannot get the Hamiltonian energy of Case 2 system. It should be expressed as follows:(13)x˙1x˙2x˙3x˙4=00ghx200ax4ex1−g−ax400−hx2−ex100x1x2x3x4.
So, the Hamiltonian energy of Case 2 system isHx=1/2x12+x22+x32+x42.


## 3. Equilibrium Point and Energy Analysis

In this work, Matlab ODE45 algorithm is utilized to solve differential equations in the numerical simulations.


### 3.1. Equilibrium Point Analysis

Since the difference between these two systems is only the coefficients of the linear terms of the first and third equations, the coefficients do not affect the type and stability of the equilibrium point. Equilibrium points of Case 2 system are analyzed for simplicity. Letx˙=0,y˙=0,z˙=0, and w˙=0. The equilibrium point set of Case 2 is origin and two lines, and the two lines are the y-axis and the w-axis. Linearize (11) at S1=0,0,0,0, and the characteristic equation is(14)λ2λ2+100=0.Solveλ1=λ2=0,λ3=10i,λ4=−10i. So, equilibrium point S1 is center. Linearize (11) at S2=0,y∗,0,0, and the characteristic equation is(15)λλ3+100+3y∗2λ−360y∗2=0.The eigenvaluesλ1=0,λ2=σ1,andλ3,4=−σ2±jω1. Here, σ1, σ2, and ω1 are the functions of variable y∗. When y∗≠0, the eigenvalues are 0,σ1,−σ2±jω1. Here, σ1=s−m/s,σ2=m/2s−s/2,ω1=3/2m/2s+s/2i, m=y∗2+100/3,ands=m3+32400y∗41/2+180y∗21/3. So, equilibrium point S2 is saddle foci.Similarly, linearize (11) at S3=0,0,0,w∗, and the characteristic equation is(16)λλ3+100+33.25w∗2λ+360w∗2=0.The eigenvaluesλ1=0,λ2=−σ3andλ3,4=σ4±jω2, where σ3, σ4, and ω2 are the functions of variable w∗. When w∗≠0, the eigenvalues are 0,−σ3,σ4±jω2, where σ3=n−p/n,σ4=p/2n−n/2,ω2=3/2p/2n+n/2i, p=133w∗2/12+100/3,andn=32400w∗4+p31/2−180w∗21/3.Finally, the equilibrium points of Case 1 system and Case 2 system are summarized in Table1. The systems are Hamiltonian conservative; therefore, these two systems cannot arrive at an equilibrium point asymptotically. This implies the nodes and spirals cannot exist in the conservative system, while the elliptic (center) and hyperbolic (saddle) equilibrium points can occur.Table 1
Equilibrium point type of Case 1 and Case 2 systems.
System
System model
Equilibrium pointy∗,w∗∈R
Eigenvalueσ,ω∈R+
Equilibrium type
Case 1
x
˙
=
h
y
w
+
c
·
d
z
y
˙
=
h
z
w
−
0.5
h
x
w
z
˙
=
−
h
y
w
−
c
·
b
x
w
˙
=
−
0.5
h
x
y
0,0,0,0
or
0
,
y
∗
,
0,0
or
0,0,0
,
w
∗
0,0
,
±
0.5
j
or
0
,
σ
5
,
−
σ
6
±
j
ω
3
or
0
,
−
σ
7
,
σ
8
±
j
ω
4
Center
Saddle
−
foci
Saddle
−
foci
Case 2
x
˙
=
h
y
w
+
g
z
y
˙
=
a
z
w
+
e
x
w
z
˙
=
−
a
y
w
−
g
x
w
˙
=
−
h
+
e
x
y
0,0,0,0
or
0
,
y
∗
,
0,0
or
0,0,0
,
w
∗
0,0
,
±
10
j
or
0
,
σ
1
,
−
σ
2
±
j
ω
1
or
0
,
−
σ
3
,
σ
4
±
j
ω
2
Center
Saddle
−
foci
Saddle
−
foci


### 3.2. The Mechanism of Chaos and Dynamic Analysis

#### 3.2.1. Case 1

The mechanism analysis of the Casimir power is an effective method in the study of the generation process of chaos. Casimir energy is conservative when no external or dissipative forces are included in systems. For Case 1 system, the Casimir energy is not conservative, which means that there must be an external torque or a dissipative torque. Arnor’d [25] put forward a Kolmogorov system in a 3D form, which described the dissipative-forced dynamical systems or the hydrodynamic instability:(17)x˙=x,H−Λx+f,where x=x1,x2,x3T, {,} is the Lie–Poisson structure, and Λx and f represent the dissipative and external forces, respectively. The Kolmogorov form is useful for studying energy cycle in systems. In general, x,H are quadratic coupling terms, and Λx is a linear term. Use this idea to separate the conservative and dissipative parts of system Σ24H. To illustrate this more clearly, equation (4) is rewritten as(18)x˙1x˙2x˙3x˙4=Π4−Π2x2x4Π1−Π4x1x4+Π4−Π3x3x4Π2−Π4x2x4Π2−Π1x1x2+Π3−Π2x2x3+cΠ3x30−cΠ1x10.From Remark1, it can be found that the term J24x∇Hx is conservative in both Hamiltonian and Casimir functions. Therefore, this term is an inertial torque, which corresponds to the first vector term to the right of the equal sign in equation (18). It is the fictitious torque generated by a free rotational rigid body with no external torque or dissipative torque. However, the term(19)JHx∇H1x=cΠ3x30−cΠ1x10T,is a nonconservative torque, which corresponds to the second part of equation (18). And this term is the external torque, which induces the energy exchanges between dissipative and supplied energy, thereby bringing about the oscillation of the Casimir power.It is worth noting that the external torque in equation (18) actually includes the dissipative torque. Since the sign of cΠ3x3 and −cΠ1x1 is uncertain, it is impossible to determine whether it is dissipative or external torque. There are two terms in the total Casimir power C˙=cΠ3−Π1x1x3: (1) cΠ3x1x3 (the upper curve in Figure 1(a)), in which the power is supplied by the external torque cΠ3x3, and (2) −cΠ1x1x3 (the lower curve in Figure 1(a)), in which the power is supplied by the external torque −cΠ1x1. The term, which is greater (less) than zero, takes a supply (dissipation) role. It was concluded that both of them were responsible for dissipative or external torque.Figure 1
Casimir power, Hamilton energy, and Casimir energy of Case 1 system. (a) Two parts of Casimir powers. (b) Hamilton energy and Casimir energy.
(a)
(b)Equations (8) and (9) are also verified by the numerical simulation results. The green line in Figure 1(b) represents the Hamiltonian energy of Case 1, which does not change over time and remains at 2.625. The irregular oscillating blue curve in Figure 1(b) is the Casimir energy. Chaotic orbits displaying Hamiltonian energy and Casimir energy are shown in Figures 2(a) and 2(b), respectively.Figure 2
Chaotic orbits of Case 1 system. (a) Hamilton energy. (b) Casimir energy.
(a)
(b)The initial value is0.5,0.5,1,1 and parameters are h,d,b=−1,2.5,2. When c=0, Case 1 system is the 4D Euler equation Σ24; therefore, it is conservative in both the Hamiltonian with Ht=2.5625 and Casimir energy with Ct=1.25. The system produces periodic states (red curve in Figure 3(a)), the Casimir power is on the zero line (red curve in Figure 3(b)). If we set c=0.2, the Hamiltonian is still conservative with Ht=2.5625, but it oscillates irregularly in the Casimir power (blue curve in Figure 3(b)), and the original periodic orbit is perturbed and turns into chaotic (blue curve in Figure 3(a)). Although the excitation is small as far as the value of c is concerned, the response of the system has essentially changed because the system has completely changed its orbital mode from regular to chaotic, as verified by the Lyapunov exponents LE=0.2405,0,0,−0.2405. The red curves in Figures 4(a) and 4(b) are the periodic orbits in 3D space of x1−x2−x4 and the phase portrait of x2−x4 produced by Case 1 without the excitation and the blue curves are the chaotic orbits produced by Case 1 under excitation, respectively.Figure 3
Time series ofx3 and Casimir powers of Case 1 system, with h,d,b=−1,2.5,2 and x0,y0,z0,w0=0.5,0.5,1,1. (a) Time series of x3. (b) Casimir powers.
(a)
(b)Figure 4
Phase portraits of Case 1 system on the different plane, withh,d,b=−1,2.5,2 and x0,y0,z0,w0=0.5,0.5,1,1. (a) 3D space of x1−x2−x4. (b)x2−x4 plane
(a)
(b)Bifurcation diagram and Lyapunov exponent spectrum of Case 1 system are shown in Figures5(a) and 5(b). Normally, the bifurcation transition occurs from source to periodic orbit, double periodic orbits, and multiperiodic orbits until chaos or an inverse evolution of dynamical bifurcation [36]. However, Case 1 system does not undergo such a bifurcation transition. It enters the chaotic mode suddenly as long as parameter c>0, as shown in the bifurcation diagram, where c varies in the interval c∈0,0.2.Figure 5
Lyapunov exponent spectrum and bifurcation diagram of Case 1 system versusc. (a) Bifurcation diagram. (b) Lyapunov exponent spectrum.
(a)
(b)


#### 3.2.2. Case 2

But there is an interesting phenomenon during the mechanism analysis of Case 2 system.In addition to the invariants such as energy and momentum that are associated with explicit symmetries, the noncanonical Hamiltonian systems generally possess Casimir invariants or distinguished functions. The canonical form is represented as(20)J=0E−E0,where E is an n×n identity matrix and 0 is an n×n zero matrix. From Remark 2,(21)x˙=Jx⋅∇Hx=00ghx200ax4ex1−g−ax400−hx2−ex100x1x2x3x4.It is obvious that the structure matrix of Case 2 system is canonical, and a canonical system has no nontrivial Casimirsfunctions [33].Since Casimir power cannot explain the chaos generation mechanism of Case 2 system, we have to find another way. Kokubu and Roussarie verified the existence of one specific heteroclinic cycle, called the singularly degenerate heteroclinic cycle, which consists of an invariant set formed by a line of equilibria together with a heteroclinic orbit connecting two of the equilibria [37]. They also found that the rupture of this cycle would separate the chaotic attractors. Are there such orbits in Case 2 system? In this part, by choosing the appropriate parameters and performing a large number of Matlab numerical simulations, the cycle was discovered. According to Section 3.1, the system has nonisolated equilibrium points Ey and Ew. The characteristic equation of Ey is(22)λλ3+hh+ey∗2+g2λ−g·a·h+ey∗2=0,where τ1=0, τ2=hh+ey∗2+g2, and τ3=−g·a·h+ey∗2, and it is easy to know that each Ey is hyperbolic. The characteristic equation of Ew is(23)λλ3+a2−h·ew∗2+g2λ+g·a·h+ew∗2=0,where ρ1=0, ρ2=a2−h·ew∗2+g2, and ρ3=g·a·h+ew∗2, and it is easy to know that each Ew is hyperbolic. According to the Routh–Hurwitz criterion, Remark 3 is established.Remark 3.
If parametersa, h, e, g, and y∗w∗ satisfy τ1>0ρ1>0, τ2>0ρ2>0), τ3>0ρ3>0, and τ1τ2−τ3>0ρ1ρ2−ρ3>0, then EyEw is a stable node-like or focus-like equilibrium point. Otherwise, EyEw is unstable or more degraded.
TakeEy as an example. Selecting appropriate parameter values and initial values, a large number of Matlab numerical simulations were carried out, and the following conclusions are obtained. If parameters a, h, e, g, and y∗ do not satisfy τ1>0, τ2>0, τ3>0, and τ1τ2−τ3>0, then the 1D (or 2D) unstable manifold WuEy1 of each saddle point or unstable focus-like Ey1=0,y1,0,0 will tend to normal hyperbolic stable manifold WsEy2 of a stable focus Ey2=0,y2,0,0, as shown in Figure 6. Case 2 system has infinite singularly degenerate heteroclinic cycles, and Figure 7 shows that chaotic attractors are generated near the singularly degenerate heteroclinic cycle.
For Case 2 system, the phase trajectory and Poincare section are drawn in Figures8–11 to verify the dynamic behaviors under different parameters. The bifurcation diagram (Figure 12(a)) and Lyapunov exponent spectrum (Figure 12(b)) shows that when g∈0,16, Case 2 system generates a chaotic flow and when g∈16,25, Case 2 system enters into quasiperiodic motion. It can be seen from the phase trajectory (Figure 8) and the Poincare section (Figure 9) that when g=10, the system is chaotic. If g is 20, the system shows a quasiperiodic state (Figures 10 and 11). It is found that the parameters h and g have great impact on the dynamic behaviors of Case 2 system through many simulations. If we fix the initial values with x0,y0,z0,w0=1,1,1,1 but change the parameters, the maps of four Lyapunov exponents related to parameters h and g can be obtained as shown in Figures 12(a) and 12(b), and the LEs with different values are highlighted with different colors. It can be seen from Figure 13 that LE1 and LE4 are symmetrical and LE2 and LE3 are symmetrical, which is also consistent with the characteristics of the Lyapunov exponent of the conservative system [38, 39].Figure 6
Singularly degenerate heteroclinic cycles of Case 2 system: seth,a,e,g=−25,−2,1,0.7 and x0,y0,z0,w0=±0.0918,y∗,±0.00469,±0.00290, where y∗=0.05,0.08. (a)x−y−w. (b) x-y.
(a)
(b)Figure 7
Phase portraits of Case 2 system: seth,a,e,g=−25,−2,1,0.7, (a)x0,y0,z0,w0=0.0918,0.08,0.00469,0.00290, and (b)x0,y0,z0,w0=−0.0918,0.08,−0.00469,−0.00290.
(a)
(b)Figure 8
Phase portraits of Case 2 system on the different plane, withh,g,a,e=0.5,10,6,5.5 and x0,y0,z0,w0=1,1,1,1. (a) x-y plane. (b) x-z plane. (c) y-w plane.
(a)
(b)
(c)Figure 9
Phase portraits of Case 2 system on the different plane, withh,g,a,e=0.5,10,6,5.5 and x0,y0,z0,w0=1,1,1,1. (a) x-y plane. (b) x-z plane. (c) y-w plane.
(a)
(b)
(c)Figure 10
Phase portraits Poincare sections of Case 2 system on the different plane, withh,g,a,e=0.5,20,6,5.5 and x0,y0,z0,w0=1,1,1,1. (a) x-y plane. (b) x-z plane. (c) y-w plane.
(a)
(b)
(c)Figure 11
Poincare sections of Case 2 system on the different plane, with (h, g, a, e) = (0.5, 20, 6, 5.5) and (x0, y0, z0, w0) = (1, 1, 1, 1). (a) x-y plane. (b) x-z plane. (c) y-w plane.
(a)
(b)
(c)Figure 12
Lyapunov exponent spectrum and bifurcation diagram of Case 2 system versusg. (a) Bifurcation diagram. (b) Lyapunov exponent spectrum.
(a)
(b)Figure 13
Dynamical evolution map of Case 2 system withh∈−10,10,g∈−5,5. (a) LE1. (b) LE2. (c) LE3. (d) LE4.
(a)
(b)
(c)
(d)


## 3.1. Equilibrium Point Analysis

Since the difference between these two systems is only the coefficients of the linear terms of the first and third equations, the coefficients do not affect the type and stability of the equilibrium point. Equilibrium points of Case 2 system are analyzed for simplicity. Letx˙=0,y˙=0,z˙=0, and w˙=0. The equilibrium point set of Case 2 is origin and two lines, and the two lines are the y-axis and the w-axis. Linearize (11) at S1=0,0,0,0, and the characteristic equation is(14)λ2λ2+100=0.Solveλ1=λ2=0,λ3=10i,λ4=−10i. So, equilibrium point S1 is center. Linearize (11) at S2=0,y∗,0,0, and the characteristic equation is(15)λλ3+100+3y∗2λ−360y∗2=0.The eigenvaluesλ1=0,λ2=σ1,andλ3,4=−σ2±jω1. Here, σ1, σ2, and ω1 are the functions of variable y∗. When y∗≠0, the eigenvalues are 0,σ1,−σ2±jω1. Here, σ1=s−m/s,σ2=m/2s−s/2,ω1=3/2m/2s+s/2i, m=y∗2+100/3,ands=m3+32400y∗41/2+180y∗21/3. So, equilibrium point S2 is saddle foci.Similarly, linearize (11) at S3=0,0,0,w∗, and the characteristic equation is(16)λλ3+100+33.25w∗2λ+360w∗2=0.The eigenvaluesλ1=0,λ2=−σ3andλ3,4=σ4±jω2, where σ3, σ4, and ω2 are the functions of variable w∗. When w∗≠0, the eigenvalues are 0,−σ3,σ4±jω2, where σ3=n−p/n,σ4=p/2n−n/2,ω2=3/2p/2n+n/2i, p=133w∗2/12+100/3,andn=32400w∗4+p31/2−180w∗21/3.Finally, the equilibrium points of Case 1 system and Case 2 system are summarized in Table1. The systems are Hamiltonian conservative; therefore, these two systems cannot arrive at an equilibrium point asymptotically. This implies the nodes and spirals cannot exist in the conservative system, while the elliptic (center) and hyperbolic (saddle) equilibrium points can occur.Table 1
Equilibrium point type of Case 1 and Case 2 systems.
System
System model
Equilibrium pointy∗,w∗∈R
Eigenvalueσ,ω∈R+
Equilibrium type
Case 1
x
˙
=
h
y
w
+
c
·
d
z
y
˙
=
h
z
w
−
0.5
h
x
w
z
˙
=
−
h
y
w
−
c
·
b
x
w
˙
=
−
0.5
h
x
y
0,0,0,0
or
0
,
y
∗
,
0,0
or
0,0,0
,
w
∗
0,0
,
±
0.5
j
or
0
,
σ
5
,
−
σ
6
±
j
ω
3
or
0
,
−
σ
7
,
σ
8
±
j
ω
4
Center
Saddle
−
foci
Saddle
−
foci
Case 2
x
˙
=
h
y
w
+
g
z
y
˙
=
a
z
w
+
e
x
w
z
˙
=
−
a
y
w
−
g
x
w
˙
=
−
h
+
e
x
y
0,0,0,0
or
0
,
y
∗
,
0,0
or
0,0,0
,
w
∗
0,0
,
±
10
j
or
0
,
σ
1
,
−
σ
2
±
j
ω
1
or
0
,
−
σ
3
,
σ
4
±
j
ω
2
Center
Saddle
−
foci
Saddle
−
foci


## 3.2. The Mechanism of Chaos and Dynamic Analysis

### 3.2.1. Case 1

The mechanism analysis of the Casimir power is an effective method in the study of the generation process of chaos. Casimir energy is conservative when no external or dissipative forces are included in systems. For Case 1 system, the Casimir energy is not conservative, which means that there must be an external torque or a dissipative torque. Arnor’d [25] put forward a Kolmogorov system in a 3D form, which described the dissipative-forced dynamical systems or the hydrodynamic instability:(17)x˙=x,H−Λx+f,where x=x1,x2,x3T, {,} is the Lie–Poisson structure, and Λx and f represent the dissipative and external forces, respectively. The Kolmogorov form is useful for studying energy cycle in systems. In general, x,H are quadratic coupling terms, and Λx is a linear term. Use this idea to separate the conservative and dissipative parts of system Σ24H. To illustrate this more clearly, equation (4) is rewritten as(18)x˙1x˙2x˙3x˙4=Π4−Π2x2x4Π1−Π4x1x4+Π4−Π3x3x4Π2−Π4x2x4Π2−Π1x1x2+Π3−Π2x2x3+cΠ3x30−cΠ1x10.From Remark1, it can be found that the term J24x∇Hx is conservative in both Hamiltonian and Casimir functions. Therefore, this term is an inertial torque, which corresponds to the first vector term to the right of the equal sign in equation (18). It is the fictitious torque generated by a free rotational rigid body with no external torque or dissipative torque. However, the term(19)JHx∇H1x=cΠ3x30−cΠ1x10T,is a nonconservative torque, which corresponds to the second part of equation (18). And this term is the external torque, which induces the energy exchanges between dissipative and supplied energy, thereby bringing about the oscillation of the Casimir power.It is worth noting that the external torque in equation (18) actually includes the dissipative torque. Since the sign of cΠ3x3 and −cΠ1x1 is uncertain, it is impossible to determine whether it is dissipative or external torque. There are two terms in the total Casimir power C˙=cΠ3−Π1x1x3: (1) cΠ3x1x3 (the upper curve in Figure 1(a)), in which the power is supplied by the external torque cΠ3x3, and (2) −cΠ1x1x3 (the lower curve in Figure 1(a)), in which the power is supplied by the external torque −cΠ1x1. The term, which is greater (less) than zero, takes a supply (dissipation) role. It was concluded that both of them were responsible for dissipative or external torque.Figure 1
Casimir power, Hamilton energy, and Casimir energy of Case 1 system. (a) Two parts of Casimir powers. (b) Hamilton energy and Casimir energy.
(a)
(b)Equations (8) and (9) are also verified by the numerical simulation results. The green line in Figure 1(b) represents the Hamiltonian energy of Case 1, which does not change over time and remains at 2.625. The irregular oscillating blue curve in Figure 1(b) is the Casimir energy. Chaotic orbits displaying Hamiltonian energy and Casimir energy are shown in Figures 2(a) and 2(b), respectively.Figure 2
Chaotic orbits of Case 1 system. (a) Hamilton energy. (b) Casimir energy.
(a)
(b)The initial value is0.5,0.5,1,1 and parameters are h,d,b=−1,2.5,2. When c=0, Case 1 system is the 4D Euler equation Σ24; therefore, it is conservative in both the Hamiltonian with Ht=2.5625 and Casimir energy with Ct=1.25. The system produces periodic states (red curve in Figure 3(a)), the Casimir power is on the zero line (red curve in Figure 3(b)). If we set c=0.2, the Hamiltonian is still conservative with Ht=2.5625, but it oscillates irregularly in the Casimir power (blue curve in Figure 3(b)), and the original periodic orbit is perturbed and turns into chaotic (blue curve in Figure 3(a)). Although the excitation is small as far as the value of c is concerned, the response of the system has essentially changed because the system has completely changed its orbital mode from regular to chaotic, as verified by the Lyapunov exponents LE=0.2405,0,0,−0.2405. The red curves in Figures 4(a) and 4(b) are the periodic orbits in 3D space of x1−x2−x4 and the phase portrait of x2−x4 produced by Case 1 without the excitation and the blue curves are the chaotic orbits produced by Case 1 under excitation, respectively.Figure 3
Time series ofx3 and Casimir powers of Case 1 system, with h,d,b=−1,2.5,2 and x0,y0,z0,w0=0.5,0.5,1,1. (a) Time series of x3. (b) Casimir powers.
(a)
(b)Figure 4
Phase portraits of Case 1 system on the different plane, withh,d,b=−1,2.5,2 and x0,y0,z0,w0=0.5,0.5,1,1. (a) 3D space of x1−x2−x4. (b)x2−x4 plane
(a)
(b)Bifurcation diagram and Lyapunov exponent spectrum of Case 1 system are shown in Figures5(a) and 5(b). Normally, the bifurcation transition occurs from source to periodic orbit, double periodic orbits, and multiperiodic orbits until chaos or an inverse evolution of dynamical bifurcation [36]. However, Case 1 system does not undergo such a bifurcation transition. It enters the chaotic mode suddenly as long as parameter c>0, as shown in the bifurcation diagram, where c varies in the interval c∈0,0.2.Figure 5
Lyapunov exponent spectrum and bifurcation diagram of Case 1 system versusc. (a) Bifurcation diagram. (b) Lyapunov exponent spectrum.
(a)
(b)


### 3.2.2. Case 2

But there is an interesting phenomenon during the mechanism analysis of Case 2 system.In addition to the invariants such as energy and momentum that are associated with explicit symmetries, the noncanonical Hamiltonian systems generally possess Casimir invariants or distinguished functions. The canonical form is represented as(20)J=0E−E0,where E is an n×n identity matrix and 0 is an n×n zero matrix. From Remark 2,(21)x˙=Jx⋅∇Hx=00ghx200ax4ex1−g−ax400−hx2−ex100x1x2x3x4.It is obvious that the structure matrix of Case 2 system is canonical, and a canonical system has no nontrivial Casimirsfunctions [33].Since Casimir power cannot explain the chaos generation mechanism of Case 2 system, we have to find another way. Kokubu and Roussarie verified the existence of one specific heteroclinic cycle, called the singularly degenerate heteroclinic cycle, which consists of an invariant set formed by a line of equilibria together with a heteroclinic orbit connecting two of the equilibria [37]. They also found that the rupture of this cycle would separate the chaotic attractors. Are there such orbits in Case 2 system? In this part, by choosing the appropriate parameters and performing a large number of Matlab numerical simulations, the cycle was discovered. According to Section 3.1, the system has nonisolated equilibrium points Ey and Ew. The characteristic equation of Ey is(22)λλ3+hh+ey∗2+g2λ−g·a·h+ey∗2=0,where τ1=0, τ2=hh+ey∗2+g2, and τ3=−g·a·h+ey∗2, and it is easy to know that each Ey is hyperbolic. The characteristic equation of Ew is(23)λλ3+a2−h·ew∗2+g2λ+g·a·h+ew∗2=0,where ρ1=0, ρ2=a2−h·ew∗2+g2, and ρ3=g·a·h+ew∗2, and it is easy to know that each Ew is hyperbolic. According to the Routh–Hurwitz criterion, Remark 3 is established.Remark 3.
If parametersa, h, e, g, and y∗w∗ satisfy τ1>0ρ1>0, τ2>0ρ2>0), τ3>0ρ3>0, and τ1τ2−τ3>0ρ1ρ2−ρ3>0, then EyEw is a stable node-like or focus-like equilibrium point. Otherwise, EyEw is unstable or more degraded.
TakeEy as an example. Selecting appropriate parameter values and initial values, a large number of Matlab numerical simulations were carried out, and the following conclusions are obtained. If parameters a, h, e, g, and y∗ do not satisfy τ1>0, τ2>0, τ3>0, and τ1τ2−τ3>0, then the 1D (or 2D) unstable manifold WuEy1 of each saddle point or unstable focus-like Ey1=0,y1,0,0 will tend to normal hyperbolic stable manifold WsEy2 of a stable focus Ey2=0,y2,0,0, as shown in Figure 6. Case 2 system has infinite singularly degenerate heteroclinic cycles, and Figure 7 shows that chaotic attractors are generated near the singularly degenerate heteroclinic cycle.
For Case 2 system, the phase trajectory and Poincare section are drawn in Figures8–11 to verify the dynamic behaviors under different parameters. The bifurcation diagram (Figure 12(a)) and Lyapunov exponent spectrum (Figure 12(b)) shows that when g∈0,16, Case 2 system generates a chaotic flow and when g∈16,25, Case 2 system enters into quasiperiodic motion. It can be seen from the phase trajectory (Figure 8) and the Poincare section (Figure 9) that when g=10, the system is chaotic. If g is 20, the system shows a quasiperiodic state (Figures 10 and 11). It is found that the parameters h and g have great impact on the dynamic behaviors of Case 2 system through many simulations. If we fix the initial values with x0,y0,z0,w0=1,1,1,1 but change the parameters, the maps of four Lyapunov exponents related to parameters h and g can be obtained as shown in Figures 12(a) and 12(b), and the LEs with different values are highlighted with different colors. It can be seen from Figure 13 that LE1 and LE4 are symmetrical and LE2 and LE3 are symmetrical, which is also consistent with the characteristics of the Lyapunov exponent of the conservative system [38, 39].Figure 6
Singularly degenerate heteroclinic cycles of Case 2 system: seth,a,e,g=−25,−2,1,0.7 and x0,y0,z0,w0=±0.0918,y∗,±0.00469,±0.00290, where y∗=0.05,0.08. (a)x−y−w. (b) x-y.
(a)
(b)Figure 7
Phase portraits of Case 2 system: seth,a,e,g=−25,−2,1,0.7, (a)x0,y0,z0,w0=0.0918,0.08,0.00469,0.00290, and (b)x0,y0,z0,w0=−0.0918,0.08,−0.00469,−0.00290.
(a)
(b)Figure 8
Phase portraits of Case 2 system on the different plane, withh,g,a,e=0.5,10,6,5.5 and x0,y0,z0,w0=1,1,1,1. (a) x-y plane. (b) x-z plane. (c) y-w plane.
(a)
(b)
(c)Figure 9
Phase portraits of Case 2 system on the different plane, withh,g,a,e=0.5,10,6,5.5 and x0,y0,z0,w0=1,1,1,1. (a) x-y plane. (b) x-z plane. (c) y-w plane.
(a)
(b)
(c)Figure 10
Phase portraits Poincare sections of Case 2 system on the different plane, withh,g,a,e=0.5,20,6,5.5 and x0,y0,z0,w0=1,1,1,1. (a) x-y plane. (b) x-z plane. (c) y-w plane.
(a)
(b)
(c)Figure 11
Poincare sections of Case 2 system on the different plane, with (h, g, a, e) = (0.5, 20, 6, 5.5) and (x0, y0, z0, w0) = (1, 1, 1, 1). (a) x-y plane. (b) x-z plane. (c) y-w plane.
(a)
(b)
(c)Figure 12
Lyapunov exponent spectrum and bifurcation diagram of Case 2 system versusg. (a) Bifurcation diagram. (b) Lyapunov exponent spectrum.
(a)
(b)Figure 13
Dynamical evolution map of Case 2 system withh∈−10,10,g∈−5,5. (a) LE1. (b) LE2. (c) LE3. (d) LE4.
(a)
(b)
(c)
(d)


## 3.2.1. Case 1

The mechanism analysis of the Casimir power is an effective method in the study of the generation process of chaos. Casimir energy is conservative when no external or dissipative forces are included in systems. For Case 1 system, the Casimir energy is not conservative, which means that there must be an external torque or a dissipative torque. Arnor’d [25] put forward a Kolmogorov system in a 3D form, which described the dissipative-forced dynamical systems or the hydrodynamic instability:(17)x˙=x,H−Λx+f,where x=x1,x2,x3T, {,} is the Lie–Poisson structure, and Λx and f represent the dissipative and external forces, respectively. The Kolmogorov form is useful for studying energy cycle in systems. In general, x,H are quadratic coupling terms, and Λx is a linear term. Use this idea to separate the conservative and dissipative parts of system Σ24H. To illustrate this more clearly, equation (4) is rewritten as(18)x˙1x˙2x˙3x˙4=Π4−Π2x2x4Π1−Π4x1x4+Π4−Π3x3x4Π2−Π4x2x4Π2−Π1x1x2+Π3−Π2x2x3+cΠ3x30−cΠ1x10.From Remark1, it can be found that the term J24x∇Hx is conservative in both Hamiltonian and Casimir functions. Therefore, this term is an inertial torque, which corresponds to the first vector term to the right of the equal sign in equation (18). It is the fictitious torque generated by a free rotational rigid body with no external torque or dissipative torque. However, the term(19)JHx∇H1x=cΠ3x30−cΠ1x10T,is a nonconservative torque, which corresponds to the second part of equation (18). And this term is the external torque, which induces the energy exchanges between dissipative and supplied energy, thereby bringing about the oscillation of the Casimir power.It is worth noting that the external torque in equation (18) actually includes the dissipative torque. Since the sign of cΠ3x3 and −cΠ1x1 is uncertain, it is impossible to determine whether it is dissipative or external torque. There are two terms in the total Casimir power C˙=cΠ3−Π1x1x3: (1) cΠ3x1x3 (the upper curve in Figure 1(a)), in which the power is supplied by the external torque cΠ3x3, and (2) −cΠ1x1x3 (the lower curve in Figure 1(a)), in which the power is supplied by the external torque −cΠ1x1. The term, which is greater (less) than zero, takes a supply (dissipation) role. It was concluded that both of them were responsible for dissipative or external torque.Figure 1
Casimir power, Hamilton energy, and Casimir energy of Case 1 system. (a) Two parts of Casimir powers. (b) Hamilton energy and Casimir energy.
(a)
(b)Equations (8) and (9) are also verified by the numerical simulation results. The green line in Figure 1(b) represents the Hamiltonian energy of Case 1, which does not change over time and remains at 2.625. The irregular oscillating blue curve in Figure 1(b) is the Casimir energy. Chaotic orbits displaying Hamiltonian energy and Casimir energy are shown in Figures 2(a) and 2(b), respectively.Figure 2
Chaotic orbits of Case 1 system. (a) Hamilton energy. (b) Casimir energy.
(a)
(b)The initial value is0.5,0.5,1,1 and parameters are h,d,b=−1,2.5,2. When c=0, Case 1 system is the 4D Euler equation Σ24; therefore, it is conservative in both the Hamiltonian with Ht=2.5625 and Casimir energy with Ct=1.25. The system produces periodic states (red curve in Figure 3(a)), the Casimir power is on the zero line (red curve in Figure 3(b)). If we set c=0.2, the Hamiltonian is still conservative with Ht=2.5625, but it oscillates irregularly in the Casimir power (blue curve in Figure 3(b)), and the original periodic orbit is perturbed and turns into chaotic (blue curve in Figure 3(a)). Although the excitation is small as far as the value of c is concerned, the response of the system has essentially changed because the system has completely changed its orbital mode from regular to chaotic, as verified by the Lyapunov exponents LE=0.2405,0,0,−0.2405. The red curves in Figures 4(a) and 4(b) are the periodic orbits in 3D space of x1−x2−x4 and the phase portrait of x2−x4 produced by Case 1 without the excitation and the blue curves are the chaotic orbits produced by Case 1 under excitation, respectively.Figure 3
Time series ofx3 and Casimir powers of Case 1 system, with h,d,b=−1,2.5,2 and x0,y0,z0,w0=0.5,0.5,1,1. (a) Time series of x3. (b) Casimir powers.
(a)
(b)Figure 4
Phase portraits of Case 1 system on the different plane, withh,d,b=−1,2.5,2 and x0,y0,z0,w0=0.5,0.5,1,1. (a) 3D space of x1−x2−x4. (b)x2−x4 plane
(a)
(b)Bifurcation diagram and Lyapunov exponent spectrum of Case 1 system are shown in Figures5(a) and 5(b). Normally, the bifurcation transition occurs from source to periodic orbit, double periodic orbits, and multiperiodic orbits until chaos or an inverse evolution of dynamical bifurcation [36]. However, Case 1 system does not undergo such a bifurcation transition. It enters the chaotic mode suddenly as long as parameter c>0, as shown in the bifurcation diagram, where c varies in the interval c∈0,0.2.Figure 5
Lyapunov exponent spectrum and bifurcation diagram of Case 1 system versusc. (a) Bifurcation diagram. (b) Lyapunov exponent spectrum.
(a)
(b)


## 3.2.2. Case 2

But there is an interesting phenomenon during the mechanism analysis of Case 2 system.In addition to the invariants such as energy and momentum that are associated with explicit symmetries, the noncanonical Hamiltonian systems generally possess Casimir invariants or distinguished functions. The canonical form is represented as(20)J=0E−E0,where E is an n×n identity matrix and 0 is an n×n zero matrix. From Remark 2,(21)x˙=Jx⋅∇Hx=00ghx200ax4ex1−g−ax400−hx2−ex100x1x2x3x4.It is obvious that the structure matrix of Case 2 system is canonical, and a canonical system has no nontrivial Casimirsfunctions [33].Since Casimir power cannot explain the chaos generation mechanism of Case 2 system, we have to find another way. Kokubu and Roussarie verified the existence of one specific heteroclinic cycle, called the singularly degenerate heteroclinic cycle, which consists of an invariant set formed by a line of equilibria together with a heteroclinic orbit connecting two of the equilibria [37]. They also found that the rupture of this cycle would separate the chaotic attractors. Are there such orbits in Case 2 system? In this part, by choosing the appropriate parameters and performing a large number of Matlab numerical simulations, the cycle was discovered. According to Section 3.1, the system has nonisolated equilibrium points Ey and Ew. The characteristic equation of Ey is(22)λλ3+hh+ey∗2+g2λ−g·a·h+ey∗2=0,where τ1=0, τ2=hh+ey∗2+g2, and τ3=−g·a·h+ey∗2, and it is easy to know that each Ey is hyperbolic. The characteristic equation of Ew is(23)λλ3+a2−h·ew∗2+g2λ+g·a·h+ew∗2=0,where ρ1=0, ρ2=a2−h·ew∗2+g2, and ρ3=g·a·h+ew∗2, and it is easy to know that each Ew is hyperbolic. According to the Routh–Hurwitz criterion, Remark 3 is established.Remark 3.
If parametersa, h, e, g, and y∗w∗ satisfy τ1>0ρ1>0, τ2>0ρ2>0), τ3>0ρ3>0, and τ1τ2−τ3>0ρ1ρ2−ρ3>0, then EyEw is a stable node-like or focus-like equilibrium point. Otherwise, EyEw is unstable or more degraded.
TakeEy as an example. Selecting appropriate parameter values and initial values, a large number of Matlab numerical simulations were carried out, and the following conclusions are obtained. If parameters a, h, e, g, and y∗ do not satisfy τ1>0, τ2>0, τ3>0, and τ1τ2−τ3>0, then the 1D (or 2D) unstable manifold WuEy1 of each saddle point or unstable focus-like Ey1=0,y1,0,0 will tend to normal hyperbolic stable manifold WsEy2 of a stable focus Ey2=0,y2,0,0, as shown in Figure 6. Case 2 system has infinite singularly degenerate heteroclinic cycles, and Figure 7 shows that chaotic attractors are generated near the singularly degenerate heteroclinic cycle.
For Case 2 system, the phase trajectory and Poincare section are drawn in Figures8–11 to verify the dynamic behaviors under different parameters. The bifurcation diagram (Figure 12(a)) and Lyapunov exponent spectrum (Figure 12(b)) shows that when g∈0,16, Case 2 system generates a chaotic flow and when g∈16,25, Case 2 system enters into quasiperiodic motion. It can be seen from the phase trajectory (Figure 8) and the Poincare section (Figure 9) that when g=10, the system is chaotic. If g is 20, the system shows a quasiperiodic state (Figures 10 and 11). It is found that the parameters h and g have great impact on the dynamic behaviors of Case 2 system through many simulations. If we fix the initial values with x0,y0,z0,w0=1,1,1,1 but change the parameters, the maps of four Lyapunov exponents related to parameters h and g can be obtained as shown in Figures 12(a) and 12(b), and the LEs with different values are highlighted with different colors. It can be seen from Figure 13 that LE1 and LE4 are symmetrical and LE2 and LE3 are symmetrical, which is also consistent with the characteristics of the Lyapunov exponent of the conservative system [38, 39].Figure 6
Singularly degenerate heteroclinic cycles of Case 2 system: seth,a,e,g=−25,−2,1,0.7 and x0,y0,z0,w0=±0.0918,y∗,±0.00469,±0.00290, where y∗=0.05,0.08. (a)x−y−w. (b) x-y.
(a)
(b)Figure 7
Phase portraits of Case 2 system: seth,a,e,g=−25,−2,1,0.7, (a)x0,y0,z0,w0=0.0918,0.08,0.00469,0.00290, and (b)x0,y0,z0,w0=−0.0918,0.08,−0.00469,−0.00290.
(a)
(b)Figure 8
Phase portraits of Case 2 system on the different plane, withh,g,a,e=0.5,10,6,5.5 and x0,y0,z0,w0=1,1,1,1. (a) x-y plane. (b) x-z plane. (c) y-w plane.
(a)
(b)
(c)Figure 9
Phase portraits of Case 2 system on the different plane, withh,g,a,e=0.5,10,6,5.5 and x0,y0,z0,w0=1,1,1,1. (a) x-y plane. (b) x-z plane. (c) y-w plane.
(a)
(b)
(c)Figure 10
Phase portraits Poincare sections of Case 2 system on the different plane, withh,g,a,e=0.5,20,6,5.5 and x0,y0,z0,w0=1,1,1,1. (a) x-y plane. (b) x-z plane. (c) y-w plane.
(a)
(b)
(c)Figure 11
Poincare sections of Case 2 system on the different plane, with (h, g, a, e) = (0.5, 20, 6, 5.5) and (x0, y0, z0, w0) = (1, 1, 1, 1). (a) x-y plane. (b) x-z plane. (c) y-w plane.
(a)
(b)
(c)Figure 12
Lyapunov exponent spectrum and bifurcation diagram of Case 2 system versusg. (a) Bifurcation diagram. (b) Lyapunov exponent spectrum.
(a)
(b)Figure 13
Dynamical evolution map of Case 2 system withh∈−10,10,g∈−5,5. (a) LE1. (b) LE2. (c) LE3. (d) LE4.
(a)
(b)
(c)
(d)


## 4. National Institute of Standards and Technology (NIST) Test and Field Programmable Gate Array (FPGA) Implementation

### 4.1. NIST Test for Two Cases

The proposed systems can be used to make a random or pseudorandom number generator. National Institute of Standards and Technology [40] provides the SP800-22 standard, which is a very common method that can be used to evaluate pseudorandom sequences identified in the field of encryption. It contains 15 tests. The number of pseudorandom sequences tested in Case 1 and Case 2 is 95, and the length of each group is 100 million bits. The test results of these two systems are listed in Tables 2 and 3, which contain the P−value of each NIST test, the distribution of the P−values, and the pass proportion of the 95 groups. For these 15 tests, each P−value is the probability that a perfect random number generator would have produced a sequence less random than the tested sequence. A significance level“α” is chosen for the tests, which represents the probability of a Type-I error, and its range is from 0.001 to 0.01. If P−value≥α, then the null hypothesis is accepted, i.e., the sequence appears to be random. If P−value<α, the sequence is nonrandom due to the rejected null hypothesis.Table 2
NIST test result of Case 1 system.
C1
C2
C3
C4
C5
C6
C7
C8
C9
C10
P value
P value
Proportion
Statistical test
11
8
9
11
11
9
10
7
12
7
0.9440
0.9707
0.989
Frequency
12
12
9
8
7
3
10
9
7
18
0.0646
0.1225
0.989
Block frequency
10
10
14
10
9
14
7
4
9
8
0.104170
0.6130
0.979
Cumulative sums
9
11
6
16
7
11
6
8
12
9
0.3595
0.5526
1
Runs
7
13
4
9
12
7
11
13
10
9
0.4580
0.6528
1
Longest run
14
6
7
11
9
9
8
16
7
8
0.2912
0.4728
0.968
Rank
14
14
6
5
11
7
10
10
6
12
0.2464
0.4150
0.968
FFT
11
9
21
5
8
6
9
6
9
11
0.0137
0.0221
1
Nonoverlapping template
5
17
8
10
12
15
4
6
10
8
0.0339
0.0622
0.989
Overlapping template
13
11
13
10
12
4
5
7
11
9
0.3074
0.4925
0.989
Universal
9
7
11
11
7
14
7
8
11
10
0.7728
0.8849
0.989
Approximate entropy
3
4
8
6
2
12
2
8
6
4
0.0288
0.0734
1
Random excursions
2
11
4
2
6
7
5
6
3
9
0.0628
0.1552
0.982
Random excursion variant
10
10
7
16
10
8
11
7
3
13
0.1626
0.2931
1
Serial
8
11
6
8
13
9
12
5
13
10
0.5229
0.7106
1
Linear complexityTable 3
NIST test result of Case 2 system.
C1
C2
C3
C4
C5
C6
C7
C8
C9
C10
P
‐
value
P
‐
value
T
Proportion
Statistical test
16
12
10
8
9
5
5
10
12
8
0.2465
0.415
1
Frequency
10
9
13
12
9
9
7
8
11
7
0.888137
0.954
0.979
Block frequency
15
12
6
13
7
10
9
4
5
14
0.0742
0.1407
1
Cumulative sums
10
8
11
7
4
14
10
10
13
8
0.4580
0.6528
0.979
Runs
12
10
9
8
4
11
7
11
12
11
0.6602
0.8148
0.989
Longest run
12
9
6
8
12
10
7
8
9
1
0.6833
0.8302
1
Rank
5
11
13
13
6
14
8
7
11
7
0.2757
0.4533
0.989
FFT
2
8
12
4
11
12
16
10
14
6
0.01724
0.029
1
Nonoverlapping template
6
13
14
7
13
13
9
8
6
6
0.2328
0.3964
0.989
Overlapping template
7
11
15
7
8
9
11
5
14
8
0.3074
0.4925
0.989
Universal
8
11
8
11
17
7
14
7
5
7
0.1184
0.2204
1
Approximate entropy
13
6
6
3
3
3
3
6
5
6
0.0519
0.1122
0.963
Random excursions
5
2
8
4
7
7
3
7
1
10
0.0856
0.1808
0.981
Random excursion variant
8
11
13
7
9
11
11
15
2
8
0.1528
0.2775
0.989
Serial
5
16
10
9
11
9
11
8
6
10
0.3973
0.593
0.989
Linear complexityAccording to the SP800-22 standard [40], the test results have to satisfy three conditions:(a)
AllP−values must be greater than the significance level, α=0.01.(b)
A proportion must pass the confidence interval formulated as(24)p^±3p^1−p^m,wherep^=1−α,m=95, implying the proportion should be between 0.9601 and 1.0298.(c)
Since theP values of the straight random sequence obey uniform distribution, the randomness can be described according to the uniformity of the distribution of P−values. The P−value is a number between 0 and 1, so we divide 0 to 1 into 10 intervals: 0.0,0.1, 0.1,0.2, 0.2,0.3, 0.3,0.4 , 0.4,0.5, 0.5,0.6, 0.6,0.7, 0.7,0.8, 0.8,0.9, and 0.9,1.0. The value of C1–C10 is the number of these 10 P−values in the corresponding interval. Since the data are divided into 95 groups, the sum of C1 to C10 is 95. When the sequence is sufficiently large, the P−value should be averagely distributed in the 10 intervals, and the frequency of occurrence in each interval is counted as Fi. If the number of sequences is m, F1+F2+⋯+F10=m, then the statistics are calculated as(25)χ2=∑i=110Fi−m/102m/10.Then, calculate theP−valueT:(26)P−valueT=igamc92,χ22.IfP−valueT>0.001, then the P−values are subject to uniform distribution; in other words, the sequence is random; otherwise, it is nonrandom.The performance of the tests is successful for the proposed two systems. From Tables2 and 3, all P−values given in the statistical tests are larger than 0.01, and condition (a) is met. The relevant proportions are within the acceptance interval 0.9601,1.0298, so condition (b) is also satisfied. Furthermore, all P−valueTs are greater than 0.001. For simplicity, the Nonoverlapping template (test no. 9 in Tables 2 and 3) is considered as an example. The distribution of P−values for nonoverlapping template is uniform (as shown in Figures 14(a) and 14(b)). Other 14 tests display a similar uniformity. Hence, condition (c) is satisfied as well. It can be concluded that Case 1 and Case 2 systems are eligible as pseudorandom generators according to the NIST standard.Figure 14
Probability distribution ofP−value for the nonoverlapping template. (a) Case 1 system. (b) Case 2 system.
(a)
(b)


### 4.2. Field Programmable Gate Array Implementation for Two Cases

The parameters of analog circuit components are affected by temperature, aging, and other environmental factors, and chaos system is extremely sensitive to the initial value of the system, and using analog circuit to realize chaotic systems is quite limited. Hence, the digital circuits are better to realize chaotic systems. Generally, there are two implementation methods: one is based on Euler algorithm and Runge–Kutta discrete algorithm to write the underlying hardware code to realize the chaotic system and the other method is using MATLAB’s DSP Builder library to build the chaotic digital circuit and automatically generate Verilog Language. This paper adopts the latter method. Because the conservative chaotic system is highly sensitive to initial values, parameters, and computational errors, its hardware implementation is more difficult than common dissipative chaotic systems. In this section, the implementation of the proposed Case 1 and Case 2 systems in FPGA is discussed. To build DSP builder model (Figure15), these two systems are first discretized.Figure 15
DSP Builder model of Case 1 and Case 2 systems.For Case 1: settingh,d,b,c=−1,2.5,2,0.2, one obtains(27)xn+1=−ΔT·yn·wn+0.5·ΔT·zn+xn,yn+1=−ΔT·zn·wn+0.5·ΔT·xn·wn+yk,zn+1=ΔT·yn·wn−0.4·ΔT·xn+zn,wn+1=0.5·ΔT·xn·yn+wn.For Case 2: settingh,g,a,e=0.5,10,6,5.5, one obtains(28)xn+1=0.5ΔT·yn·wn+10·ΔT·zn+xn,yn+1=6ΔT·zn·wn+5.5·ΔT·xn·wn+yk,zn+1=−6ΔT·yn·wn−10·ΔT·xn+zn,wn+1=−6·ΔT·xn·yn+wn.Secondly, the chaotic digital circuits are built in Matlab/Simulink/DSP Builder, with the fixed step of10−2s, initial values 0.5,0.5,1,1, 1,1,1,1, and sampling time T=10−5s. The system models (Figure 15) are converted to the Verilog Hardware Description Language and finally synthesized and compiled to the register transfer level code. The process of compiling and downloading to FPGA core board is completed through QuartusII 15.1. Table 4 shows the consumption of system resources. The occupancy rate of total logical elements is 12%, and the occupancy rate of total memory bits is less than 1%. The field programmable gate array board is used to perform digital implementation, and a 14-bit dual-channel DA converter AD764 is used to convert time series. The converted data are displayed by a digital oscilloscope (Tektronix MSO4054B) to demonstrate phase portraits (Figure 16(a)). The hardware design of the circuit implementation is shown in Figure 16(b). In order to compare the computer simulation with the digital circuit results, the phase diagrams of Case 1 and Case 2 generated by the field programmable gate array are shown in Figures 17 and 18, respectively. It can be seen that the phase diagrams of these two systems are consistent with the results of computer simulation.Table 4
System resource consumption.
Logic elements
Multiplier elements
Registers
Memory bits
17260
48
11932
5690Figure 16
(a) The experimental circuit of the chaotic orbits. (b) Hardware design.
(a)
(b)Figure 17
Phase portraits of Case 1 chaotic system using FPGA. (a)x-y plane. (b) x-z plane. (c) y-w plane.
(a)
(b)
(c)Figure 18
Phase portraits of Case 2 chaotic system using FPGA. (a)x-y plane. (b) x-z plane. (c) y-w plane.
(a)
(b)
(c)


## 4.1. NIST Test for Two Cases

The proposed systems can be used to make a random or pseudorandom number generator. National Institute of Standards and Technology [40] provides the SP800-22 standard, which is a very common method that can be used to evaluate pseudorandom sequences identified in the field of encryption. It contains 15 tests. The number of pseudorandom sequences tested in Case 1 and Case 2 is 95, and the length of each group is 100 million bits. The test results of these two systems are listed in Tables 2 and 3, which contain the P−value of each NIST test, the distribution of the P−values, and the pass proportion of the 95 groups. For these 15 tests, each P−value is the probability that a perfect random number generator would have produced a sequence less random than the tested sequence. A significance level“α” is chosen for the tests, which represents the probability of a Type-I error, and its range is from 0.001 to 0.01. If P−value≥α, then the null hypothesis is accepted, i.e., the sequence appears to be random. If P−value<α, the sequence is nonrandom due to the rejected null hypothesis.Table 2
NIST test result of Case 1 system.
C1
C2
C3
C4
C5
C6
C7
C8
C9
C10
P value
P value
Proportion
Statistical test
11
8
9
11
11
9
10
7
12
7
0.9440
0.9707
0.989
Frequency
12
12
9
8
7
3
10
9
7
18
0.0646
0.1225
0.989
Block frequency
10
10
14
10
9
14
7
4
9
8
0.104170
0.6130
0.979
Cumulative sums
9
11
6
16
7
11
6
8
12
9
0.3595
0.5526
1
Runs
7
13
4
9
12
7
11
13
10
9
0.4580
0.6528
1
Longest run
14
6
7
11
9
9
8
16
7
8
0.2912
0.4728
0.968
Rank
14
14
6
5
11
7
10
10
6
12
0.2464
0.4150
0.968
FFT
11
9
21
5
8
6
9
6
9
11
0.0137
0.0221
1
Nonoverlapping template
5
17
8
10
12
15
4
6
10
8
0.0339
0.0622
0.989
Overlapping template
13
11
13
10
12
4
5
7
11
9
0.3074
0.4925
0.989
Universal
9
7
11
11
7
14
7
8
11
10
0.7728
0.8849
0.989
Approximate entropy
3
4
8
6
2
12
2
8
6
4
0.0288
0.0734
1
Random excursions
2
11
4
2
6
7
5
6
3
9
0.0628
0.1552
0.982
Random excursion variant
10
10
7
16
10
8
11
7
3
13
0.1626
0.2931
1
Serial
8
11
6
8
13
9
12
5
13
10
0.5229
0.7106
1
Linear complexityTable 3
NIST test result of Case 2 system.
C1
C2
C3
C4
C5
C6
C7
C8
C9
C10
P
‐
value
P
‐
value
T
Proportion
Statistical test
16
12
10
8
9
5
5
10
12
8
0.2465
0.415
1
Frequency
10
9
13
12
9
9
7
8
11
7
0.888137
0.954
0.979
Block frequency
15
12
6
13
7
10
9
4
5
14
0.0742
0.1407
1
Cumulative sums
10
8
11
7
4
14
10
10
13
8
0.4580
0.6528
0.979
Runs
12
10
9
8
4
11
7
11
12
11
0.6602
0.8148
0.989
Longest run
12
9
6
8
12
10
7
8
9
1
0.6833
0.8302
1
Rank
5
11
13
13
6
14
8
7
11
7
0.2757
0.4533
0.989
FFT
2
8
12
4
11
12
16
10
14
6
0.01724
0.029
1
Nonoverlapping template
6
13
14
7
13
13
9
8
6
6
0.2328
0.3964
0.989
Overlapping template
7
11
15
7
8
9
11
5
14
8
0.3074
0.4925
0.989
Universal
8
11
8
11
17
7
14
7
5
7
0.1184
0.2204
1
Approximate entropy
13
6
6
3
3
3
3
6
5
6
0.0519
0.1122
0.963
Random excursions
5
2
8
4
7
7
3
7
1
10
0.0856
0.1808
0.981
Random excursion variant
8
11
13
7
9
11
11
15
2
8
0.1528
0.2775
0.989
Serial
5
16
10
9
11
9
11
8
6
10
0.3973
0.593
0.989
Linear complexityAccording to the SP800-22 standard [40], the test results have to satisfy three conditions:(a)
AllP−values must be greater than the significance level, α=0.01.(b)
A proportion must pass the confidence interval formulated as(24)p^±3p^1−p^m,wherep^=1−α,m=95, implying the proportion should be between 0.9601 and 1.0298.(c)
Since theP values of the straight random sequence obey uniform distribution, the randomness can be described according to the uniformity of the distribution of P−values. The P−value is a number between 0 and 1, so we divide 0 to 1 into 10 intervals: 0.0,0.1, 0.1,0.2, 0.2,0.3, 0.3,0.4 , 0.4,0.5, 0.5,0.6, 0.6,0.7, 0.7,0.8, 0.8,0.9, and 0.9,1.0. The value of C1–C10 is the number of these 10 P−values in the corresponding interval. Since the data are divided into 95 groups, the sum of C1 to C10 is 95. When the sequence is sufficiently large, the P−value should be averagely distributed in the 10 intervals, and the frequency of occurrence in each interval is counted as Fi. If the number of sequences is m, F1+F2+⋯+F10=m, then the statistics are calculated as(25)χ2=∑i=110Fi−m/102m/10.Then, calculate theP−valueT:(26)P−valueT=igamc92,χ22.IfP−valueT>0.001, then the P−values are subject to uniform distribution; in other words, the sequence is random; otherwise, it is nonrandom.The performance of the tests is successful for the proposed two systems. From Tables2 and 3, all P−values given in the statistical tests are larger than 0.01, and condition (a) is met. The relevant proportions are within the acceptance interval 0.9601,1.0298, so condition (b) is also satisfied. Furthermore, all P−valueTs are greater than 0.001. For simplicity, the Nonoverlapping template (test no. 9 in Tables 2 and 3) is considered as an example. The distribution of P−values for nonoverlapping template is uniform (as shown in Figures 14(a) and 14(b)). Other 14 tests display a similar uniformity. Hence, condition (c) is satisfied as well. It can be concluded that Case 1 and Case 2 systems are eligible as pseudorandom generators according to the NIST standard.Figure 14
Probability distribution ofP−value for the nonoverlapping template. (a) Case 1 system. (b) Case 2 system.
(a)
(b)


## 4.2. Field Programmable Gate Array Implementation for Two Cases

The parameters of analog circuit components are affected by temperature, aging, and other environmental factors, and chaos system is extremely sensitive to the initial value of the system, and using analog circuit to realize chaotic systems is quite limited. Hence, the digital circuits are better to realize chaotic systems. Generally, there are two implementation methods: one is based on Euler algorithm and Runge–Kutta discrete algorithm to write the underlying hardware code to realize the chaotic system and the other method is using MATLAB’s DSP Builder library to build the chaotic digital circuit and automatically generate Verilog Language. This paper adopts the latter method. Because the conservative chaotic system is highly sensitive to initial values, parameters, and computational errors, its hardware implementation is more difficult than common dissipative chaotic systems. In this section, the implementation of the proposed Case 1 and Case 2 systems in FPGA is discussed. To build DSP builder model (Figure15), these two systems are first discretized.Figure 15
DSP Builder model of Case 1 and Case 2 systems.For Case 1: settingh,d,b,c=−1,2.5,2,0.2, one obtains(27)xn+1=−ΔT·yn·wn+0.5·ΔT·zn+xn,yn+1=−ΔT·zn·wn+0.5·ΔT·xn·wn+yk,zn+1=ΔT·yn·wn−0.4·ΔT·xn+zn,wn+1=0.5·ΔT·xn·yn+wn.For Case 2: settingh,g,a,e=0.5,10,6,5.5, one obtains(28)xn+1=0.5ΔT·yn·wn+10·ΔT·zn+xn,yn+1=6ΔT·zn·wn+5.5·ΔT·xn·wn+yk,zn+1=−6ΔT·yn·wn−10·ΔT·xn+zn,wn+1=−6·ΔT·xn·yn+wn.Secondly, the chaotic digital circuits are built in Matlab/Simulink/DSP Builder, with the fixed step of10−2s, initial values 0.5,0.5,1,1, 1,1,1,1, and sampling time T=10−5s. The system models (Figure 15) are converted to the Verilog Hardware Description Language and finally synthesized and compiled to the register transfer level code. The process of compiling and downloading to FPGA core board is completed through QuartusII 15.1. Table 4 shows the consumption of system resources. The occupancy rate of total logical elements is 12%, and the occupancy rate of total memory bits is less than 1%. The field programmable gate array board is used to perform digital implementation, and a 14-bit dual-channel DA converter AD764 is used to convert time series. The converted data are displayed by a digital oscilloscope (Tektronix MSO4054B) to demonstrate phase portraits (Figure 16(a)). The hardware design of the circuit implementation is shown in Figure 16(b). In order to compare the computer simulation with the digital circuit results, the phase diagrams of Case 1 and Case 2 generated by the field programmable gate array are shown in Figures 17 and 18, respectively. It can be seen that the phase diagrams of these two systems are consistent with the results of computer simulation.Table 4
System resource consumption.
Logic elements
Multiplier elements
Registers
Memory bits
17260
48
11932
5690Figure 16
(a) The experimental circuit of the chaotic orbits. (b) Hardware design.
(a)
(b)Figure 17
Phase portraits of Case 1 chaotic system using FPGA. (a)x-y plane. (b) x-z plane. (c) y-w plane.
(a)
(b)
(c)Figure 18
Phase portraits of Case 2 chaotic system using FPGA. (a)x-y plane. (b) x-z plane. (c) y-w plane.
(a)
(b)
(c)


## 5. Conclusion

This paper constructed two HCCSs by two methods. The first method is based on 4D Euler equations, which generate chaos by breaking the Casimir energy conservation. The second method is based on a proposed method of constructing a class of Hamiltonian conservative chaotic system and then choosing the appropriate functions to construct HCCS. The dynamic behaviors of the proposed two HCCSs have been investigated based on the equilibrium point analysis, singular orbits, and energy analysis. The NIST tests indicate that the proposed two new HCCSs have good application values in information encryption. Finally, a pseudorandom chaotic signal generator is implemented based on the field programmable gate array digital development platform, which is useful for potential engineering applications.



---

*Source: 8693157-2020-04-25.xml*