Chaotic dynamics of a three particle array under Lennard--Jones type forces and a fixed area constraint
aa r X i v : . [ m a t h . D S ] S e p Chaotic dynamics of a three particle array underLennard–Jones type forces and a fixed area constraint
Pablo V. Negr´on–Marrero ∗ Department of MathematicsUniversity of Puerto RicoHumacao, PR 00791-4300
Abstract
We consider the dynamical problem for a system of three particles in which theinter–particle forces are given as the gradient of a Lennard–Jones type potential.Furthermore we assume that the three particle array is subject to the constraintof fixed area. The corresponding mathematical problem is that of a conservativedynamical system over the manifold determined by the area constraint. We studynumerically the stability of this system. In particular, using the recently introducedmeasure of chaos by Hunt and Ott (2015), we study numerically the possibility ofchaotic behaviour for this system.
Keywords : dynamical system, chaos, expansion entropy, constrained optimization
AMS subject classifications:
In this paper we consider the dynamical problem for a system of three particles, wherethe interactions between the particles of the system are due to forces given as the gradientof a Lennard–Jones [12] type potential. Furthermore we assume that the three particlearray is subjected to the additional constraint that the area of the triangle generated bythem is fixed. The motivation to consider the area constraint comes from the followingphenomena observed both in laboratory experiments and molecular dynamics simulations(see e.g., [1, 2]). As the density of a fluid is progressively lowered (keeping the temperatureconstant), there is a certain “critical” density such that if the density of the fluid is lowerthan this critical value, then bubbles or regions with very low density appear within thefluid. This phenomenon is usually called “cavitation” and it has been extensively studiedas well in solids. (See for instance [4, 9] for discussions and further references.) ∗ email: [email protected] A c after which they become unstable. Furthermore, usingtechniques of equivariant bifurcation theory they showed that bifurcating from the pointcorresponding to A c , there are solutions curves corresponding to isosceles triangles thestability of which could only be determined numerically. Moreover for an instance ofthe Lennard–Jones potential, they showed numerically that stable equilibrium pointscorresponding to scalene triangles exist for a range of values of the area parameter.The study of chaotic behaviour in dynamical systems is a very important and activearea in this field. Verifying whether or not a particular system is chaotic can be a difficulttask. In Section 5 we review a new criteria for chaos introduced recently by Hunt and Ott[10]. This new characterization of chaos is based on the so called expansion entropy (cf.(5.5)) which does not require the identification of a compact invariant set, and naturallyleads itself to a computational method for detecting chaotic behaviour in a dynamicalsystem. In Section 5 we describe this numerical scheme and give several examples of theuse of this method applied to various dynamical systems of known chaotic or non–chaoticbehaviour. Finally in Section 6 we use this method to show that the dynamical system forthe three particle array under the area constraint, exhibits chaotic behaviour essentiallyfor all values of the area parameter in the area constraint. We consider a system of three equal particles interacting via an inter–particle potential φ . The potentials we are concerned are those of the form: φ ( r ) = c r δ − c r δ , (2.1)where c , c are positive constants and δ > δ >
2. These constants determine thephysical properties of the particles or molecules in the array. (The classical
Lennard–Jones [12] potential is obtained upon setting δ = 12 and δ = 6.) We study the planardynamics of such a system subject to the constraint of fixed area for the triangular array.2et r i ( t ), i = 1 , , r ij ( t ) = r i ( t ) − r j ( t )for i < j . Let r i = || r i || = ( r i · r i ) , r ij = || r ij || . If m is the particle mass, then the kinetic and potential energies of the system are givenrespectively by: K = m X i =1 || ˙ r i || , U = X i 116 ( a + b + c ) . (2.3)For any A > g ( r , r , r ) = Γ( r , r , r ) − A . (2.4)Thus if we set g = 0, then we would be specifying that the triangle defined by thepositions of the three particles is of fixed area A . We shall need the following formulasfor the partial derivatives of g : ∂g∂ r = 1 r ∂ Γ ∂a ( r − r ) + 1 r ∂ Γ ∂b ( r − r ) , (2.5a) ∂g∂ r = − r ∂ Γ ∂a ( r − r ) + 1 r ∂ Γ ∂c ( r − r ) , (2.5b) ∂g∂ r = − r ∂ Γ ∂b ( r − r ) − r ∂ Γ ∂c ( r − r ) , (2.5c)where the partial derivatives of Γ are evaluated at ( r , r , r ). Note that ∂g∂ r + ∂g∂ r + ∂g∂ r = . (2.6) We assume that the particle corresponding to i = 3 is fixed at the origin, while that for i = 2 is restricted to move along the x –axis. Using the notation r i = ( u i , w i ), i = 1 , , C functions r , r , r such that:min E ( r , r , r )subject to g ( r , r , r ) = 0 , r = , w = 0 . m ¨ r = − φ ′ ( r ) r ( r − r ) − φ ′ ( r ) r ( r − r ) − λ ∂g∂ r ( r , r , r ) ,m ¨ r = φ ′ ( r ) r ( r − r ) − φ ′ ( r ) r ( r − r ) − λ ∂g∂ r ( r , r , r ) − µ e , (2.7) m ¨ r = φ ′ ( r ) r ( r − r ) + φ ′ ( r ) r ( r − r ) − λ ∂g∂ r ( r , r , r ) − ν e − ν e , g ( r , r , r ) , r = , w = 0 . Here λ , ν , ν , µ (which are functions of t ) are the Lagrange multipliers correspondingto the constraints in the problem, and { e , e } is the standard basis for R . From theseequations and (2.6), it follows now that m X i =1 ¨ r i ( t ) = − ν e − ( ν + µ ) e , which describes the motion of the center of mass of the system. The system (2.7) isconservative as the energy K + U is conserved along solutions.Applying the constraints r = and w = 0, the equations above simplify to: m ¨ r = − φ ′ ( r ) r ( r − r ) − φ ′ ( r ) r r − λ ∂g∂ r ,m ¨ u = φ ′ ( r ) r ( u − u ) − φ ′ ( r ) r u − λ ∂g∂u , φ ′ ( r ) r w − λ ∂g∂w − µ, = φ ′ ( r ) r r + φ ′ ( r ) r r − λ ∂g∂ r − ν e − ν e , g. where g and all of its partial derivatives are evaluated at ( u , w , u , , , r = | u | , the first two equations together with the last one, can be solved indepen-dently of the third and fourth equations. Once r , u , and λ are determined, one can find ν , ν , µ using the third and fourth equations. Thus we are lead to consider the systemin component form: m ¨ u = − φ ′ ( r ) r ( u − u ) − φ ′ ( r ) r u − λ ∂γ∂u ( u , w , u ) ,m ¨ w = − (cid:20) φ ′ ( r ) r + φ ′ ( r ) r (cid:21) w − λ ∂γ∂w ( u , w , u ) , (2.8) m ¨ u = φ ′ ( r ) r ( u − u ) − φ ′ ( r ) r u − λ ∂γ∂u ( u , w , u ) , γ ( u , w , u ) , γ ( u , w , u ) = g ( u , w , u , , , r = q ( u − u ) + w , r = q u + w , r = | u | . (2.9)It is easy to check now that γ ( u , w , u ) = 14 w u − A . Thus the area constraint γ = 0 is simply the square of the familiar formula of “areaequals one half base times height”. We shall exploit this very shortly but before doingso, we pause to characterize the equilibrium points of the system (2.8). We now show that the equilibrium points of the system (2.8) are precisely those studiedin [14]. We remark that while λ in (2.8) is in general a function of t , in the followingdiscussion it is a constant corresponding to a steady state. Proposition 3.1. The equilibrium points of the system (2.8) correspond to position vec-tors r = ( u , w ) and r = ( u , for which a = r , b = r , and c = r (cf. (2.3) , (2.9) )satisfy: φ ′ ( a ) + λ Γ ,a = 0 ,φ ′ ( b ) + λ Γ ,b = 0 ,φ ′ ( c ) + λ Γ ,c = 0 , Γ( a, b, c ) = A , (3.1) where Γ ,a = ∂ Γ ∂a , etc.Proof : After setting the right hand side of (2.8) equal to zero, we are led to φ ′ ( a ) a ( u − u ) + φ ′ ( b ) b u + λ ∂γ∂u = 0 , (3.2a) φ ′ ( a ) a w + φ ′ ( b ) b w + λ ∂γ∂w = 0 , (3.2b) − φ ′ ( a ) a ( u − u ) + φ ′ ( c ) c u + λ ∂γ∂u = 0 , (3.2c) γ ( u , w , u ) = 0 , (3.2d)where we have set a = r , b = r , and c = r . The constraint γ ( u , w , u ) = 0 isequivalent to Γ( a, b, c ) = A . Since A > r = ( u , w ), r = ( u , w = 0 , u = 0 . w = 0 and r = together with the definitions of a, b, c , we get from(2.5a) and the first component of (2.5b) that: ∂γ∂u = 1 a Γ ,a ( u − u ) + 1 b Γ ,b u ,∂γ∂w = 1 a Γ ,a w + 1 b Γ ,b w ,∂γ∂u = − a Γ ,a ( u − u ) + 1 c Γ ,c u . If we now substitute these equations into (3.2), and add and subtract (3.2a) and (3.2c),then the first three equations of (3.2) are equivalent to: H b u + H c u = 0 , (2 H a + H b ) u − (2 H a + H c ) u = 0 , (3.3)( H a + H b ) w = 0 , where H a = 1 a ( φ ′ ( a ) + λ Γ ,a ) , etc.Since w = 0 we must have that H a = − H b . Since u = 0 the determinant of thecoefficient matrix of the first two equations in the system (3.3) must be zero. A calculationshows that this determinant is 2 H a . Hence H a = 0, and since H b = − H a , then H b = 0.It follows now from the first equation of (3.3) and using that u = 0, that H c = 0. Since H a = H b = H c = 0 is equivalent to the first three equations in (3.1), the result follows.For future reference we record here the basic result in [14] concerning the existenceand multiplicity of solutions of the system (3.1). Theorem 3.2. For any value of A > , the system (3.1) has a solution of with a = b = c = a A corresponding to an equilateral triangle with corresponding Lagrange–multiplier λ A where: a A = 2 √ A √ , λ A = − φ ′ ( a A ) a A . (3.4) This equilibrium point is stable, that is, a minimizer of the potential energy functional U in (2.2) subject to the fixed area constraint, if and only if ρ ( A ) ≡ φ ′′ ( a A ) + 3 a A φ ′ ( a A ) > . (3.5) Moreover, if A is a simple root of ρ ( · ) , then there exist three branches of solutionscorresponding to isosceles triangles, bifurcating from the branch of equilateral triangles atthe point where A = A . The stability of the bifurcating branches in Theorem 3.2 can only be determinednumerically. In [14] they give numerical examples in which these bifurcations are of thetrans–critical type and also that there can be secondary bifurcations into stable scalenetriangles. 6 The reduced problem The system (2.8) is an example of a dynamical system over a manifold. One can inprinciple solve this system directly using some of the numerical techniques for thesetypes of problems, that essentially are based on variations of a predictor–corrector (orprojection) method (cf. [5], [6]). However because of the simplification of the constraint γ = 0 in (2.8), the system (2.8) can be reduced further to one in the variables ( u , w ).This simplifies greatly the calculations in Section 6 when we apply to our problem thenew criteria for detecting chaos in [10].As we mentioned before, the constraint γ = 0 reduces to:14 w u − A = 0 , Since w = 0 and u = 0 and A > 0, we may assume that w > u > 0. It followsnow that ∂γ∂u = 0 , ∂γ∂w = 12 w u , ∂γ∂u = 12 w u . (4.1)Using these and the fact that the constraint γ = 0 is equivalent to w u = 2 A , we caneliminate λ and the equation of u from (2.8), to get the following reduced system for u , w : m ¨ u = − φ ′ ( r ) r ( u − u ) − φ ′ ( r ) r u , (4.2a) m (cid:2) w + u (cid:3) ¨ w = 2 A (cid:20) − φ ′ ( r ) r ( u − u ) + φ ′ ( u ) (cid:21) − (cid:20) φ ′ ( r ) r + φ ′ ( r ) r (cid:21) w + 8 mA ( ˙ w ) w . (4.2b)In these equations w > u should be replaced with 2 A/w .For the numerical calculations of the next section, we shall need to transform thesystem (4.2) into one of first order. If we let v = ˙ u and v = ˙ w , then the above systemis equivalent to: ˙ u = v , ˙ w = v ,m ˙ v = − φ ′ ( r ) r ( u − u ) − φ ′ ( r ) r u , (4.3a) m (cid:2) w + u (cid:3) ˙ v = 2 A (cid:20) − φ ′ ( r ) r ( u − u ) + φ ′ ( u ) (cid:21) − (cid:20) φ ′ ( r ) r + φ ′ ( r ) r (cid:21) w + 8 mA v w . (4.3b)Given any ( p, q, α, β ) with q > 0, it follows from the standard existence and uniquenesstheorem for ode’s that these equations have a unique solution satisfying u (0) = p, w (0) = q, v (0) = α, v (0) = β. 7f we let r = ( u , w , v , v ) and p = ( p, q, α, β ), then we shall write r ( t ; p ) to denote thedependence of the solution of (4.3) on the initial conditions p . It follows from standardresults on the dependence on initial values for the solutions of initial value problems (cf.[7]), that r ( t ; p ) is a differentiable function of p . The study of chaotic behaviour in dynamical systems is a very important and active areain this field. A commonly used definition of chaos, originally formulated by Robert L.Devaney , says that for a dynamical system to be classified as chaotic , it must have thefollowing properties [8]: • it must be sensitive to initial conditions • it must be topologically mixing • it must have dense periodic orbitsVerifying whether these properties hold or not for a particular system can be a difficulttask. Recently Hunt and Ott [10] developed a new criteria for chaos based on the so called expansion entropy (cf. (5.5)) which does not require the identification of a compactinvariant set. In this section we introduce some of the notions and definitions in [10]leading to a practical, from the numerical point of view, characterization of chaos for adynamical system.We consider an autonomous dynamical system of the form:˙ r ( t ) = f ( r ( t )) , (5.1)where f : Ω → R n is a smooth function, and Ω ⊂ R n is open with a smooth boundary.We denote by r ( · ; p ) the solution of (5.1) such that r (0) = p . r ( · ; p ) is called the flow or evolution operator of (5.1). Let u ( · ; p ) = D p r ( · ; p ) . (5.2)Note that u ( t ; p ) is an n × n matrix for each t . Now since r (0; p ) = p , it follows that u (0; p ) = I , where I is the n × n identity matrix. After differentiating in (5.1) withrespect to p , we get that ˙ u ( t ) = D r f ( r ( t )) u ( t ) . Thus, the flow r ( · ; p ) and the matrix function u ( · ; p ) are solutions of the initial valueproblem: ˙ u ( t ) = D r f ( r ( t )) u ( t ) , (5.3a)˙ r ( t ) = f ( r ( t )) , (5.3b) u (0) = I , r (0) = p . (5.3c)8he function u is central in the definition of chaos in [10]. In particular, let G ( u ) bethe product of the singular values of u that are greater than 1. (If none of the singularvalues are greater than 1, we set G ( u ) = 1.) For any subset S of Ω, called a restrainingregion , let S T = { p ∈ S : r ( t ; p ) ∈ S, t ∈ [0 , T ] } , and define E T ( r , S ) = 1 | S | Z S T G ( u ( T ; p )) d p , (5.4)where | S | is the volume of S . The expansion entropy H ( r , S ) of the flow r over the set S is defined as: H ( r , S ) = lim T →∞ ln E T ( r , S ) T . (5.5)According to [10] the system (5.1) is chaotic if for some restraining region S , we havethat H ( r , S ) > Remark 5.1. The expansion entropy function has several nice properties, one of thembeing that H ( r , S ′ ) ≤ H ( r , S ) whenever S ′ ⊂ S . Thus if chaos is detected for anyregion S ′ , it will also be detected for any other region containing S ′ . Remark 5.2. The definitions in this section can be extended to non–autonomous systemsand for any manifold Ω in R n (cf. [10]). Example 5.3. We consider the constant coefficient linear system˙ r ( t ) = A r ( t ) , (5.6)where A is an n × n matrix. The solution of this system is given by r ( t ) = e tA p , t ∈ R , where p ∈ R n is arbitrary. Here e tA = ∞ X k =0 t k k ! A k . It follows from the above representation of the solution r ( · ) that the linear dynamicalsystem above is non chaotic. We show now that the definition of chaos in [10] givesindeed a non chaotic behaviour in this case.For simplicity of exposition we assume that A is symmetric with eigenvalues λ ≤ λ ≤ · · · ≤ λ n and corresponding orthonormal eigenvectors { x , x , . . . , x n } . In thiscase, the matrix e tA is symmetric as well with eigenvalues e λ t , e λ t ,...,e λ n t with the samecorresponding eigenvectors as A . Since the eigenvalues of e tA are all positive, it followsthat the singular values of e tA coincide with its eigenvalues. Since the matrix u in (5.2)is given by e tA , we get that G ( u ( t ; p )) = e t P λi> λ i , 9r one if the sum in the exponent is empty. In particular G ( u ( t ; p )) is independent of p .Writing p = P ni =1 c i x i , we get that r ( t ) = n X i =1 c i e λ i t x i . For the restraining region, we consider the set S = ( p = n X i =1 c i x i : | c i | ≤ a ) , for any a > 0. If all λ i are nonpositive, then S T = S , G ( u ( T ; p )) = 1, and we get that E T ( r , S ) = 1, and thus that H ( r , S ) = 0.Assume now that there are p > λ i . In this case S T = ( n X i =1 c i x i : | c i | ≤ a, ≤ i ≤ n − p, | c i | ≤ a e − λ i T , n − p + 1 ≤ i ≤ n ) , and G ( u ( T ; p )) = e T P ni = n − p +1 λ i . As | S T | = (2 a ) n e − T P ni = n − p +1 λ i , we get that E T ( r , S ) = 1 | S | Z S T G ( u ( T ; p )) d p , = 1(2 a ) n G ( u ( T ; p )) | S T | = 1 . Once again we get that that H ( r , S ) = 0.For an arbitrary compact restraining region S ′ , we take a in the definition of S above sufficiently large such that S ′ ⊂ S . It follows from Remark 5.1 that H ( r , S ′ ) ≤ H ( r , S ) = 0, and thus that the linear dynamical system is non chaotic according to [10].The definition of H ( r , S ) given in this section leads naturally to a very practicalnumerical scheme for estimating this quantity and thus to numerically test for chaos fora given dynamical system. If we take a set { p , . . . , p N } of uniformly random vectors in S , then we can approximate (5.4) by ˆ E T ( r , S ) wereˆ E T ( r , S ) = 1 N N X k =1 p k ∈ ST G ( u ( T ; p k )) . (5.7)This sum is computed for different values of T ∈ [0 , T max] for some prescribed T max.The whole computation is repeated for different random samples of the { p , . . . , p N } . Wethen compute the average over the samples of ln[ ˆ E T ( r , S )], for each of the chosen T ’s in[0 , T max]. From a plot of these averages vs T , we can estimate H ( r , S ) as the asymptoticslope of this graph. (See [10].) We close this section with an application of this numericalscheme to several dynamical systems whose possible chaotic or non–chaotic behaviour isknown. The computations in this section and the rest of the paper were performed usingthe ode45 routine of MATLAB using an event subroutine to detect when an orbit mayexit S for the first time. 10 10 20 30 40 5000.10.20.30.40.50.60.70.8 Figure 1: A simulation of the Hunt and Ott algorithm on a linear system with constantcoefficients. Example 5.4. As a first example we consider the special case of (5.6) in which A = (cid:20) − . − . − . − . (cid:21) . We already shown that this system is non–chaotic according to the Hunt and Ott criteria.The coefficient matrix has eigenvalues − . 1. The result of using the numericalscheme described above with N = 5000 and 40 samples, and restraining region S = { ( x, y ) : | x | ≤ , | y | ≤ } , is shown in Figure 1. The slope of approximately 8 × − of the best line in this case isconsistent with a non chaotic system. Example 5.5. We consider the following dynamical system called a “Sprott system ofcase A”: ˙ x ( t ) = y, ˙ y ( t ) = − x + yz, ˙ z ( t ) = 1 − y . This system is an example of a chaotic dynamical system with no equilibrium points(cf. [11], [16]). We tested on this system the numerical scheme described above with N = 10000 and 40 samples, and restraining region S = { ( x, y, z ) : | x | ≤ , | y | ≤ , | z | ≤ } . We show in Figure 2 a plot of the average ln[ ˆ E T ( r , S )] vs T with their correspondingerror bars (in terms of their variances). The positive slope of the best linear fitting in thepicture (indicative of the chaotic behaviour of the system), is approximately 0 . S for all forward time [10]. 11 20 40 60 80 100-0.500.511.522.533.54 Figure 2: A simulation of the Hunt and Ott algorithm on a chaotic system with noequilibrium points. Example 5.6. We consider now the motion of a double pendulum consisting of twopoint masses m , m joined together by two massless shafts or bars of lengths l , l . (SeeFigure 3.) We assume that there is no friction on the joints. If we let θ i ( t ) to be theangle that the shaft l i makes with the vertical direction, i = 1 , 2, then the equations forthe motion of such a system are given after simplification by: l ¨ θ ( t ) = T m sin( θ ( t ) − θ ( t )) − g sin θ ( t ) ,l ¨ θ ( t ) = − T m sin( θ ( t ) − θ ( t )) , where g is the constant of gravity acceleration and T , T , the tensions in the shafts, aresolutions of the system (cid:20) m − m cos( θ ( t ) − θ ( t )) − m cos( θ ( t ) − θ ( t )) m + m (cid:21) (cid:20) T T (cid:21) = (cid:20) l ˙ θ ( t ) + g cos θ ( t ) l ˙ θ ( t ) (cid:21) . This system is well known to be chaotic (cf. [3], [15]). The Hunt and Ott algorithm wastested in this case with N = 1000 and 20 samples, and restraining region S = n ( θ , ˙ θ , θ , ˙ θ ) : | θ i | ≤ , (cid:12)(cid:12)(cid:12) ˙ θ i (cid:12)(cid:12)(cid:12) ≤ , i = 1 , o . In Figure 4 we show the results for this simulation. The slope of approximately 5 . In this section we present various numerical results for the behaviour of solutions of (4.3)when the potential φ is given by (2.1). In particular, using the Hunt and Ott method,we test the system (4.3) for chaotic behaviour.12igure 3: Frictionless double pendulum. Figure 4: Simulation using the Hunt and Ott Frictionless double pendulum.13 .5 0.55 0.6 0.65 0.7 0.7511.051.11.151.21.251.31.351.41.451.5 A a scalene trianglesisosceles triangles(a=b) equilateral triangles Figure 5: Bifurcation diagram for the equilibrium points of the system (4.3) for theLennard-Jones potential (2.1) with values (6.1), green representing stable equilibria whilered unstable ones.We consider the case of (2.1) in which the parameters are given by c = 1 , c = 2 , δ = 12 , δ = 6 . (6.1)The stability and multiplicity of the equilibrium states of (4.3) (cf. Proposition 3.1) wasfully analysed in [14] for this case. We briefly review those results that are more relevantto the present discussion. In Figure 5 (reproduced from [14]) we show a projection of thebifurcation diagram of equilibrium states for the particular case (6.1). (The notation hereis as in Proposition 3.1 where a, b, c represent the sides of the triangular configuration ofthe array.) In this figure there are three bifurcation points at the following approximatevalues of the area parameter A in (4.3): A = 0 . A =0 . A = 0 . A − = 0 . A . The stability, multiplicity andtype of the equilibrium point is summarized in Table 1. The “type” column refers as towhether the shape of triangular array is an equilateral, isosceles, or scalene triangle.In the first simulation we examine the dynamics of the system (4.3) near the equilib-rium point for the value A = 0 . 55. This equilibrium point corresponds to an equilateralconfiguration of the array with sides 1 . u = 0 . w = 0 . − on one of the sides of the triangular array and used the resulting values of u , w , withvelocities v , v both O (10 − ), as initial conditions for (4.3). In Figure 6 we show theprojection onto the u – w plane of the computed orbit of the dynamical system. Figure7 shows the evolution of u ( t ) and w ( t ) as functions of time. The initial point is markedin blue and the equilibrium point in green. Both figures are consistent with a stable (not14 type stability multiplicity A ∈ (0 , A − ) equilateral stable oneequilateral stable one A ∈ ( A − , A ) isosceles stable threeisosceles unstable threeequilateral unstable one A ∈ ( A , A ) isosceles stable threeisosceles unstable threeequilateral unstable one A ∈ ( A , A ) scalene stable sixisosceles unstable sixequilateral unstable one A ∈ ( A , ∞ ) isosceles stable threeisosceles unstable threeTable 1: Type, stability, and multiplicity for the equilibrium points of (4.3) for the valuesin (6.1).asymptotically stable) fixed point. Similar results are obtained for initial values close toother stable equilibrium corresponding to different values of A .We now consider the case in which A = 0 . 65. According to Table 1 we have oneequilateral configuration which is unstable, six isosceles unstable configurations, and sixstable equilibrium points corresponding to scalene triangles. For the scalene case, thetriangle sides are 1 . , . , . u – w plane of anorbit generated with an initial point (indicated in blue) not necessarily close to any of thesix stable equilibrium points (in green). Not that the orbit appears to visit “regularly”all the stable equilibrium points. In Figure 9 we show the evolution of u ( t ) and w ( t ) asfunctions of time, with the apparent random nature characteristic of a chaotic system.In Figure 10 we tested the orbit in Figure 9 for sensitivity to initial conditions. Thisfigure was generated introducing a perturbation of O (10 − ) into the initial conditionused for Figure 9. The resulting figure differs substantially from that in Figure 10, againcharacteristic of a chaotic system.To test for possible chaotic behaviour of the system (4.3), we used the numericalscheme described at the end of Section 5 to approximate the expansion entropy H forour system. Note that the entropy depends now on the area parameter A . As a restrainingregion we used: S = { ( u , w , v , v ) : − ≤ u , v , v ≤ , < w ≤ } . We computed (5.7) with N = 5000 and 50 data samples. In Figure 11 we show a plot Both orbits in Figures 9 and 10 were computed setting the absolute and relative error tolerances inode45 to 10 − . .562 0.5625 0.563 0.5635 0.564 0.5645 0.565 0.56550.97520.97540.97560.97580.9760.97620.97640.97660.97680.977 u w Figure 6: Sample orbit for initial point (blue) near the stable equilibrium point (green)corresponding to A = 0 . t u t w Figure 7: Graphs of u and w vs t for the components of the orbit in Figure 6 corre-sponding to A = 0 . 55. 16 u w Figure 8: Sample orbit visiting all six stable equilibrium point (green) corresponding to A = 0 . t u t w Figure 9: Graphs of u and w vs t for the components of the orbit in Figure 8 corre-sponding to A = 0 . 65. 17 50 100 150 200 250 30000.20.40.60.811.21.4 t u t w Figure 10: Graphs of u and w vs t for the components of the orbit with initial conditionsas in Figure 9 plus a perturbation O (10 − ) corresponding to A = 0 . E T ( r , S )] for a range of values of T corresponding to A = 0 . 65. The red bars in the graph give intervals of plus or minus one sample standarddeviation from each computed mean. As can be seen from the figure, the computedaverages lie approximately on a line for large values of T . The slope of this line gives anapproximation of H ( r , S ) for our system. This slope is approximately 2 . A = 0 . 65. We performeda similar calculation for A = 0 . 55. We show in Figure 12 the corresponding graph, againwith the computed averages lying approximately on a line for large values of T , withslope of approximately 3 . A as well. The criteria for chaos given by Hunt and Ott [10] leads itself to a practical numericalmethod for detecting chaos in dynamical systems. The method is applicable to con-tinuous as well as to discrete dynamical systems, even non–autonomous systems. Thecomputations can be very intensive as the method requires a large number of randompoints over the restraining region, and this must be repeated another number of times inorder to compute the required averages. However the method is easy to run in parallelwhich helps to reduce the computational time.The values of A = 0 . 55 and 0 . 65 are typical values for different ranges of this parame-ter. Thus we expect chaotic behaviour as well for values of A in certain intervals. For the18 20 40 60 80 100-100-50050100150200250 Figure 11: Graph of ln( ˆ E T ) vs T for the system (4.3) corresponding to A = 0 . Figure 12: Graph of ln( ˆ E T ) vs T for the system (4.3) corresponding to A = 0 . A = 0 . 55, the equilibrium point of the system (4.3) is stable.If this equilibrium point is unique, then the system must have a hidden strange attractor.Examples of this type of dynamical systems are rather limited (cf. [13]). Further analysisin this direction as well as for the case A = 0 . 65, and the implications of the results inthis paper to the study of “cavitation” in fluids mentioned in the introduction, shall bepursued elsewhere. References [1] Bazhirov, T.T., Norman, G.E., Stegailov, V.V., Cavitation in liquid metals undernegative pressures , Molecular dynamics modeling and simulation. J. Phys. Condens.Matter , , 114113.[2] Blander, M., Katz, J., Bubble Nucleation in Liquids . AlChE J. , , 833–848.[3] Deleanu, D., Concerning the behavior of the harmonically forced double pendulum ,Universitatii Maritime Constanta, Vol. 12, Issue 16, pp. 229–236, 2011.[4] Fond, C., Cavitation Criterion for Rubber Materials: A Review of Void-GrowthModels , J. Polym. Sci. Part B Polym. Phys. , , 2081–2096, 2001.[5] Hairer, E., Geometric Integration of Ordinary Differential Equations on Manifolds ,BIT Numerical Mathematics, 41, pp. 996–1007, 2001.[6] Hairer, E., Solving Differential Equations on Manifolds, Universit´e de Gen`eve, Lec-ture Notes, June 2011.[7] Hartman, P., Ordinary Differential Equations, Birkh¨auser, Boston, 1982.[8] Hasselblatt, B. and Katok, A., A First Course in Dynamics with a Panorama ofRecent Developments, Cambridge University Press, 2003.[9] Horgan, C.O. and Polignone, D.A., Cavitation in nonlinearly elastic solids: A review , Appl. Mech. Rev. , , 471–485, 1995.[10] Hunt, B. R. and Ott, E., Defining chaos , Chaos, 25, 097618, 2015.[11] Jafari, S., Sprott, J., and Nazarimehr, F., Recent new examples of hidden attractors ,The European Physical Journal Special Topics, 224(8):1469-1476, 2015.[12] Lennard–Jones, J.E., On the Determination of Molecular Fields , Proc. R. Soc. Lond.A 106 (738), 463–477, 1924.[13] Molaie, M., Jafari, S., Sprott, J., and Golpayegani, S., Simple chaotic flows with onestable equilibrium , International Journal of Bifurcation and Chaos, Vol. 23, No. 11,pages 1350188, 2013. 2014] Negr´on-Marrero, P. V. and L´opez-Serrano, M., Minimal energy configurations offinite molecular arrays , Symmetry, 11, 158, 2019.[15] Shinbrot, T., Grebogi, C., Wisdom, J., and Yorke, J., Chaos in a double pendulum .American Journal of Physics, American Association of Physics Teachers, 60, pp.491–499, 1992.[16] Sprott, J., Some simple chaotic flows , Physical Review E50, pp. R647–R650, 1994.