Anton A Lipovka
email: aal@cifus.uson.mx

Physics on the adiabatically changed Finslerian manifold and cosmology

In present paper we confirm our previous result [4] that Planck constant is adiabatic invariant of electromagnetic field propagating on the adiabatically changed Finslerian manifold. Direct calculation from cosmological parameters gives value h=6x10(-27) (erg s). We also confirm that Planck constant (and hence other fundamental constants which depend on h) is varied on time due to changing of geometry.

As an example the variation of the fine structure constant is calculated. Its relative variation ((da/dt)/a) consist 1.0x10(-18) (1/s).

We show that on the Finsler manifold characterized by adiabatically changed geometry, classical free electromagnetic field is quantized geometrically, from the properties of the manifold in such manner that adiabatic invariant of field is ET=6x10(-27)=h.

Electrodynamic equations on the Finslerian manifold are suggested. It is stressed that quantization naturally appears from these equations and is provoked by adiabatically changed geometry of manifold. We consider in details two direct consequences of the equations: i) cosmological redshift of photons and ii) effects of Aharonov -Bohm that immediately follow from equations. It is shown that quantization of system consists of electromagnetic field and baryonic components (like atoms) is obvious and has clear explanation.

Introduction

The problem of formulation of general theory which could naturally unify General Relativity (GR) and Quantum Theory (QT) is of most fundamental and actual problem of modern physics. But to resolve this problem we should first to know what is the Planck constant. This is the question which open the door and allow us to find unified theory for all branches of physics. To obtain this key we have first to remember some important facts.

1) Planck constant in quantum mechanics always appears together with derivatives and in equations them appears in the same power. This fact clearly points out on the possible relation between Planck constant and geometry.

2) Einstein [1] and later Debye [2] have shown from thermodynamical approach that electromagnetic field is quantized alone, without any assumption on the nature of oscillators. So this is another hint that it should be quantized from geometry.

3) Recently the cornerstone result was announced [3] in respect to the quantization at very small scales. It was found that at the small scales the behavior (movements) of quantum system becomes classical one (see paper [3] for details). This result also argue that QT is pure geometrical phenomenon and it disappears at small scales when geometry becomes to be approximately Euclidean.

4) Recently we have publish result which clearly point out into geometrical origin of the Planck constant [4]. In this paper we have shown that the Planck constant is adiabatic invariant of electromagnetic field on the adiabatically changed Finslerian manifold. From this fact the variation of the Planck constant on the time directly follows (and hence variations of fundamental constants, for example variation of the fine structure constant, due to adiabatically changed geometry [4,5]).

On the one hand we have these serious arguments for the geometric origin of quantization, but on the other hand we have also serious problems with Riemannian geometry. General relativity was created as theory on (pseudo-) Riemannian manifold (we will use farther in this paper "Riemannian" instead of "pseudo-Remannian"). Such choice was not favored by some serious arguments but only by the fact the real geometry is curved. Metrical function L(x) in this case depends only on the coordinates and this fact applies some restrictions on the theory, and lead to the serious problem with singularity. However, as it follows from experiments on measurement of the cosmological constant, our Universe is expanding with acceleration and for this reason complete geometry of Finsler should be favoured. On the Finslerian manifold the metric function depends not only of coordinates, but also on velocities L = L(x,

• x) and (as it will be shown in this paper) this fact allows us to introduce the cosmological constant by natural way, calculate from geometry correct value of Planck constant and unify gravity, electrodynamics and QT. From the Finsler geometry in natural way follows Hubble constant, cosmological constant, quantization (as we will show in this paper) and much more and we can conclude that Finsler geometry naturally comply with all observational data.

In this paper we obtain classical equations of motion for a system on the manyfold with adiabatically changed geometry and supplied by Finslerian metrical function L(x,

•

x). We show how the Planck's constant naturally appears from geometry and within (3+1) formalism write exact relation between h on the one hand and scalar curvature and cosmological constant on the other hand.

We show how the classical electromagnetic field is naturally quantized due to existence of adiabatic invariant of the field on adiabatically changed Finslerian manifold.

Finally we write equations of electrodynamics, which classically (geometrically) describe quantization of electromagnetic field. We clearly show two important particular cases for these equations: i) case of free electromagnetic field, when quantization appears from geometry and leads to losses of energy by propagating photon (so called cosmological redshift), and ii) the Aharonov -Bohm effects, which immediately follows from equations.

All-around in this paper we suppose that Latin indexes i, j, k, l, m = 1, 2, 3 and greek α, β, ...µ, ν... = 0, 1, 2, 3 . Signature of metric is (1, -1, -1, -1) .

2 Adiabatic invariant and general formalism let M be an 3-dimensional, class C 3 manifold characterized by scalar curvature R =2/R 2 , where R is the curvature radius in point x k , where x k is a local coordinate on an open subset U ⊂ M . Let suppose M be supplied by Finsler metric L(x,

•

x) and write a 1-parameter family of hypersurfaces on the X n defined by equation

S m (x k , • x k ) = S M (x k , • x k ) (1)
Here S m stay for matter action and S M corresponds to a 1 -parameter family of hypersurfaces on the M .

Our aim is to write the geodesic equations (Hamilton or Lagrange -Euler equations) for this general case. In classical physics right hand part is constant or zero and we get usual classical Hamilton (or Lagrange -Euler) equations. In the case of General Relativity we put there the only invariant we have in pseudo-Riemann geometry -scalar curvature of manifold R. Let us consider what happens in general case, when the right-hand term consist of the adiabatically changed parameters of the Finslerian manifold. In this case by varying ( 1)

δS m = δL m (x k , • x k )dt = δ(p k dx k -Hdt) = δS M = δL M (x k , • x k )dt (2)
we immediately obtain Hamilton -like equations

dp k dt = - ∂H ∂x k - ∂L M ∂x k (3)
and

dx k dt = ∂H ∂p k + ∂L M ∂p k . ( 4
)
One can see there appears an additional force in right -hand term of equation (3) and an additional velocity in (4), which we naturally can attribute to cosmological constant (acceleration) and to the Hubble constant (v = Hx). It actually corresponds to the fact that there no exist absolutely closed systems and these additional terms appear due to adiabatic changes of geometry (tensor of metric) because of expansion of the Universe.

Absolutely the same way we can write equations of Lagrange -Euler varying (1):

∂L m ∂x k - d dt ∂L m ∂ • x k = ∂L M ∂x k - d dt ∂L M ∂ • x k . ( 5
)
As one can see in right part of this equation again appear two additional terms due to expansion of the Universe (due to changing of geometry of the manifold, as the system under consideration is moving).

Exact Planck's constant value and quantization of electromagnetic field

Let us calculate value of the Planck constant from the parameters which characterize the Finslerian manifold. Consider a generalized system distributed over volume. Let T p (M ) and T * p (M ) be respectively tangent and cotangent bundles on M , where p α ∈ T p (M ) and p α ∈ T * p (M ) are covariant and contravariant components of corresponding impulse.

Total momentum density of our system, in unit volume summed over all directions is given by expression

p = c 3 8πG R = c 3 4πG 1 R 2 . ( 6
)
Here, as usually, R =2/R 2 is the scalar curvature and R is the radius of curvature.

But on the Finsler manifold the tensor of metric depends on x and

•

x for this reason the momentum p 2 = g µν (x, • x)p µ p ν also depends on x and • x. Let's fix the spatial part of the coordinate system x 0 . In this case we have p(t,

• x) and, hence R(t,

• x) too. So we can write

δp = c 3 2πG 1 R 3 δR(t, • x) = c 3 2πG 1 R 3 δR(t) + ∂R(H) ∂H δH , (7)
here we note that R(t,

•

x) = R(t, Hx 0 ) where H is Hubble constant. But by taking into account that

∂R ∂t = ∂ ∂t c 2H = - c 2H 2 ∂H ∂t (8)
we can write

∂R ∂H = ∂R ∂t ∂t ∂H = - c 2H 2 . ( 9
)
Let us consider variation of Hubble constant.

From the relation

• x = Hx we can find δ • x = xδH + Hδx (10)
or

xδH = •• xδt -Hδx (11)
But the only cosmological acceleration we have experimentally measured is associated with the cosmological constant Λ , so we can write for this variation

δH = c 2 Λ -H 2 δt ( 12
)
Substituting these expressions into (7) we find

δp = 2c 3 H πG 2H 2 c 2 -Λ δt , (13)
or taking into account that R =2/R 2 and by changing volume from the spherical coordinates to the euclidean ones, we obtain

δp = 3c 3 H 8π 2 G (R -4Λ) δt . ( 14
)
This is variation of momentum (in unit volume in 3 directions) of our generalized system located on the Finslerian manifold, due to adiabatically changed geometry.

Now we are ready to write complete adiabatic invariant for a free propagating electromagnetic field. The components of the 4 -momentum p of free electromagnetic field propagating on the Finslerian manifold with adiabatically changed geometry are varied on time. This variation proceeds adiabatically and can be considered as lineal function i.e. for energy ε of the field, for example, we have

δε ε = - δt t (15)
so, the adiabatic invariant we are interested in is

εt = - δε δt t 2 . ( 16
) But δε = cδp k . ( 17
)
By substituting (14) into (16) we can write finally (we divide δp by factor 3 because we are interesting only in one direction of the momentum)

εt = -ct 2 ∂p ∂t = - c 4 H 8π 2 G (R -4Λ) t 2 = η 0 (18)
so, for one second and measured values of H = 73 kms -1 M pc -1 = 2.4 • 10 -18 s -1 and Λ = 1.7 • 10 -56 cm -2 [6] we have for this adiabatic invariant η 0 = h = 6 • 10 -27 (ergs • s.) as it should be. In the same way we can obtain similar relations for other components of 4-momentum:

p γ x γ = η γ ( 19
)
(there is no summation over γ in this relation and for the photon propagating in direction x 3 the components p 1 = p 2 = 0). Introduced here 4-vector η γ has components in unit volume:

η γ = (h, h, h, h) . ( 20
)
Here the adiabatic invariant for electromagnetic field (Planck constant), which depends clearly on the parameters of the manifold R and Λ (and consequently depends on time) is:

h = - c 4 H 8π 2 G (R -4Λ) t 2 (21)
that give for unit time at present epoch h = 6•10 -27 (ergs • s.) as it was mentioned above. From this relation it is easy to see the Planck constant depends on time as h ∼ 1/t.

Hilbert integral

Now let's consider integral of Hilbert for particular case of the free electromagnetic field propagating along geodesic on the adiabatically changed Finslerian manifold:

p α dq α = ∆S M , (22)
where p α is 4-momentum of the field, q α is generalized coordinate and right hand term δS M , as before, corresponds to the changing of the 1-parameter family of hypersurfaces due to adiabatic variation of geometry as system under consideration is moving on M . For electromagnetic field with Lagrangian L m = F µν F µν /16π we have from ( 22)

1 c ∂L m ∂A µ,ν A µ,σ dx σ = ∆S M (23)
By taking into account that tensor of energy -momentum is

T ν σ = ∂L m ∂A µ,ν A µ,σ -δ ν σ L m (24)
we can write for propagating classical electromagnetic field

1 c (T ν σ + δ ν σ L m ) dx σ = ∆S M (25)
If the field propagate in the direction x 3 , than electric field E = E 1 and magnetic field H = H 2 . So the only non-zero components of the field tensor F µν are F 01 = -F 10 = E 1 ; F 13 = -F 31 = -H 2 and hence T 00 = T 33 = T 30 = (E 2 + H 2 )/8π .

In this case for the 00 -component, for example, by taking into account (18) we have

E 2 + H 2 8π = h T = hν (26)
and similar relations one can write for other components. So, as one can see the classical electromagnetic field is quantized due to adiabatic variation of the Finslerian manifold and we do not need some artificial methods to quantize it.

Variation of the fine structure constant

As it was shown before [5] even on the Riemannian manifold the value of the fine structure constant is changed adiabatically on time. In the case of Finslerian manifold this variation is smaller by factor 2/3 due to presence of the cosmological constant. To show this let's start from (14). For one direction (divided by factor 3) we have from ( 14):

δp = c 3 H 8π 2 G (R -4Λ) δt . ( 27
)
But the fine structure constant is α = V /c where V is electron velocity at the first Bohr orbit. Momentum in this case we can write as

P = mαc √ 1 -α 2 (28) so δP = mc (1 -α 2 ) 3/2 δα (29)
and

δα = 1 -α 2 3/2 Hc 3 mc8π 2 G (R -4Λ) δt (30)
that give us value [5]), H = 73 kms -1 M pc -1 = 2.4 • 10 -18 s -1 and Λ = 1.7 • 10 -56 cm -2 ).

• α/α = -1.03 • 10 -18 (for 1 second), R =2/R 2 , R = c/2H (see

Electrodynamics on the Finslerian manifold

In Riemannian geometry the first pair of equations of electrodynamics is follows directly from the properties of the field tensor.

F µν = A ν;µ -A µ;ν = A ν,µ -A µ,ν (31)
where

A µ;ν = A µ,ν -Γ σ µν A σ ( 32
)
And for this reason (Γ σ µν = Γ σ νµ ) the first pair of equations on Riemannian manifold with constant scalar curvature is

∂ σ F µν + ∂ µ F νσ + ∂ ν F σµ = 0 (33)
On the Finsler manifold we can obtain the first pair in the same way, but in this case the field tensor is

∽ F µν = A ν;µ -A µ;ν (34)
where covariant differentials DA µ include now terms with the Cartan connections

C µνσ = 1 2 ∂gµν ∂ • x σ like this C σ µ ν A σ d • x ν
and also Γ σ µ ν A σ dx ν . In the most important case, we are interested in, when the scalar curvature is small and the tensor of metric has spatial structure described by the Robertson-Walker metric, the additional terms C σ µν A σ d

• x ν and Γ σ µν A σ dx ν in covariant derivative can be evaluated easily as ≈ A σ dx ν /R (here 1/R = 2H/c = 5 • 10 -29 is the inverse radius of scalar curvature) so we have

∽ F µν = A νµ -A µν -t µν ( 35
)
where our estimation for the small components consist

t µν ≈ A µ dx ν /R ≈ A µ • 10 -29 .
As one can see these components probably will be significant only in the vicinity of black holes and can be omitted in our present consideration.

For this reason the first pair of the equations can be written as

∼ F µν;σ + ∼ F νσ;µ + ∼ F σµ;ν = 0 (36)
or, by taking account our estimations

∂ σ F µν + ∂ µ F νσ + ∂ ν F σµ = O(t µν,σ ) (37)
So one can see that the first pair of electrodynamic equations remains to be the first pair of the Maxwell equations with high precision.

The second pair of equations of electrodynamics follows directly from variation of functional (1) if we consider a charge characterized by 4-current j α , and the electromagnetic field on the Finslerian manifold:

S m = S M (38)
Here S M as before corresponds to family of the hypersurfaces on the expanded manifold. By varying S m we have

δS m = - 1 c Ω 1 c j α δA α + 1 16π δ(F µν F µν ) dΩ (39)
here we put F µν instead of ∽ F µν because, as we have seen, small additional terms, corresponding to small components of t µν are insignificant in the case of Robertson-Walker metric. Integrating the second term by parts, we obtain

δS m = - 1 c Ω 1 c j µ + 1 4π ∂F µν ∂x ν A µ,σ δx σ dΩ (40)
By varying S M we have (see eq. (16 -20))

δS M = Ω η σ (x σ ) 2 δx σ dΩ (41)
where, as it was shown before, η σ = (h, h, h, h) in unit volume (here h is the Planck constant). The equations under discussion one can write as follows

1 c Ω 1 c j µ + 1 4π ∂F µν ∂x ν A µ,σ δx σ dΩ = - Ω η σ (x σ ) 2 δx σ dΩ + O(η 2 σ ) (42)
or finally 1 c

1 c j µ + 1 4π ∂F µν ∂x ν A µ,σ = - η σ (x σ ) 2 + O(η 2 σ ) (43)
(there is no summation over σ here). This is the second pair of equations of electrodynamics on the adiabatically changed Finslerian manifold. The bounded electromagnetic field (second term) in this case is explicitly included into consideration, as it take place in the case of Bohmian formalism when this field appears in QT as quantum potential (see [7] for details and also results of paper [8]). Here we consider two important cases which immediately follow from these equations.

Cosmological redshift

It is well known as the photon propagate through expanding universe its frequency (or wave length) is changed. This loss of energy by free electromagnetic field, named as cosmological redshift, appears in our equations by natural way as losses of the energy by photon due to adiabatically changed geometry of manifold.

4πc

∂F µν ∂x ν A µ,σ = 1 8πc ∂(F µν F µν ) ∂x σ = - η σ (x σ ) 2 + O(η 2 σ ) (44)
(there is no summation over σ).

The Aharonov -Bohm effects

Another important case that follows directly from the second pair of equations is the Aharonov -Bohm effects. As is known, a necessary condition for the existence of the Aharonov-Bohm (AB) effects is the presence, in the overall structure of the equations, of the "zero field" potentials which cannot be removed by gauge transformations and do not create electromagnetic fields [9,10]. These "zero -potentials" are the result of "non-trivial topology" of the area on which the particle moves [9,[START_REF] Bagrov | Trudy FIAN[END_REF]12]. Such a situation arises in electrodynamics of anisotropic media where the structure of Maxwell's equations eliminates the possibility of satisfying the boundary conditions. To satisfy regularly the boundary conditions in anisotropic media, usually introduce the zero -potential, which do not create electromagnetic fields (see [9] and references therein). In the case of adiabatically expanding Finslerian manifold, the anisotropy of space occurs for any moving body automatically as right part of eq. ( 43).Therefore, it is safe to say that in the case of AB effects we are dealing directly with the anisotropy of space, due to adiabatically changed Finslerian manifold as the particle moves along its trajectory. In this case, the role of the zero potentials (which do not generate electromagnetic fields) performs variable geometry as it follows directly from (43).

In absence of electric and magnetic fields on the path of propagation of the particle under consideration, the second term disappears (but it still take place inside of the solenoid and affects our particle: "In spite of the fact that the magnetic field vanishes out of the solenoid, the phase shift in the wave functions is proportional to the corresponding magnetic flux inside of the solenoid" [13] ) and we obtain (we neglect here by the small term O(η

2 σ )) 1 c 2 j µ A µ,σ = - η σ (x σ ) 2 (45)
but j µ = (ρc, ρV k ) (here ρ is charge density and V k is 3-velocity) so if we put ρ = e and remember that δj µ = 0 (for this reason A µ ∂ σ j µ = 0 and A µ,σ j µ = ∂ σ (A µ j µ )) we obtain by using the Gauss theorem

e c A 0 = - η 0 (x 0 ) , e c A k V k = - η k (x k ) . ( 48
)
These equations describe the electric and magnetic effects of Aharonov -Bohm (here x 0 and x k are fixed). Namely for µ = 0 we have for the phase variation ∆Φ

e t ϕdt = -h ∆t 1 sec . = -h∆Φ (49) mc du µ ds - e c F µν u ν = - η µ (x µ ) 2 . ( 55
)
By taking into account that classical period for orbital movement of electron is

T = πe 2 m 2 |E| 3 , (56)
in classical limit v << c the straightforward calculations give the energy for first Bohr orbit, obtained from classical electrodynamic on Finslerian manifold:

E 1 = me 4 2 2 (57)
that coincide with quantum calculations. Relativistic corrections are obvious.

Conclusions

In this paper we confirm our previous result [4] that Planck constant is adiabatic invariant of electromagnetic field propagating on the adiabatically changed Finslerian manifold. Direct calculation of the Planck constant value made from cosmological parameters gives h = 6 • 10 -27 (ergs • s.) that is in great agreement with measured value. We also confirm that Planck constant (and hence other fundamental constants which depend on h) is varied on time due to changing of geometry of the manifold.

As an example we suggest calculation of the fine structure constant variation. Obtained value consist • α/α = -1.03 • 10 -18 (for 1 second) and this variation is expected to be measured in nearest future.

We show that on the Finslerian manifold characterized by adiabatically changed geometry, classical free electromagnetic field is quantized geometrically, from the properties of the manifold.

Equations for electrodynamics on the Finslerian manifold are suggested. It is shown that quantization naturally appears from these equations and is provoked by adiabatically changed geometry of manifold. We consider in details two direct consequences of the equations -cosmological redshift of photons and effects of Aharonov -Bohm.

Finally we show how quantization of systems which consist of electromagnetic field and baryonic components (like atoms and molecules) appears.

electric effect of Aharonov -Bohm, and when µ = k (here k = 1, 2, 3) we have relation

describes magnetic effect of Aharonov -Bohm.

To conclude this part we would like to stress again that whereas the bounded field (second term in (43)) do not appears in these relations, it actually affects the moving particle through potentials A µ [13] and this bounded field corresponds to quantum potential in the Bohmian formalism [7,8].

Complete theory

In previous part we have obtained electrodynamic equations. They are applied in the case when the movement of charge (or j µ ) is defined. To construct self-consistent theory we should treat j µ as variable from the beginning.

We consider now a charge characterized by 4-current j α , and the electromagnetic field on the Finslerian manifold:

Here S M as before corresponds to family of the hypersurfaces on the expanded manifold. In complete form, when δj µ = 0 , we have the action

which describes quantum properties of our system. It is clear there are a lot of different systems and applications. For this reason it is impossible to write here a general theory, but as an example let us consider the hydrogen atom. In order to coincide with quantum mechanical calculations we should neglect by third term in (52) which corresponds to quantum potential [7,8] and gives the zero energy correction (for example in the case of harmonic oscillator it gives term 1/2 in expression for energy [4,7]). In this case by varying (52) we have

and we can write


G W Luo

Hopf-flip bifurcations of vibratory systems with impacts

Keywords: Vibration, Clearance, Impact map, Periodic motion, Hopf-flip bifurcation

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Impact oscillators arise whenever the components of a vibrating system collide with rigid obstacles or with each other. Such systems with repeated impacts exist in a wide variety of engineering applications, particularly in mechanisms and machines with clearances or gaps. The principle of operation of vibration hammers, impact dampers, shakers, pile drivers, machinery for compacting, milling and forming, etc., is based on the impact action for moving bodies. With other equipment, e.g., mechanisms with clearances, heat exchangers, fuel elements of nuclear reactors, gears, piping systems, wheel-rail interaction of high speed railway coaches, etc., impacts also occur, but they are undesirable as they bring about failures, strain, shorter service life and increased noise levels. Researches into repeated impact dynamics have important significance in optimization design of machinery with rigid obstacles or clearances, noise suppression and reliability analyses, etc. The physical process during impacts is of strongly non-linear and discontinuous characteristics. The presence of the non-linearity and discontinuity complicates the dynamic analysis of repeated impact systems considerably, but it can be described theoretically and numerically by discontinuities in good agreement with reality. Compared with single impact, vibro-impact dynamics is more complicated, and hence, has received great 1 attention. The large interest in analyzing and understanding the performance of such systems is reflected by vast and ever increasing amount of research effort devoted in this area in recent years, a small sample of which is reported in Refs. [START_REF] Aidanpää | Periodic and chaotic behaviour of a threshold-limited two-degree-of-freedom system[END_REF][START_REF] Balachandran | Dynamics of an elastic structure excited by harmonic and a harmonic impactor motions[END_REF][START_REF] Budd | The effect of frequency and clearance variations on single-degree-of-freedom impact oscillators[END_REF][START_REF] Souza | Calculation of Lyapunov exponents in systems with impacts[END_REF][START_REF] Feigin | The increasingly complex structure of the bifurcation tree of a piecewise-smooth system[END_REF][START_REF] Feng | Modeling of the mean Poincare map on a class of random impact oscillators[END_REF][START_REF] Huang | Stationary response of multi-degree-of-freedom vibro-impact systems under white noise excitations[END_REF][START_REF] Ivanov | Stabilization of an impact oscillator near grazing incidence owing to resonance[END_REF][START_REF] Knudsen | Dynamic stability of weakly damped oscillators with elastic impacts and wear[END_REF][START_REF] Leine | Nonlinear dynamics and modeling of various wooden toys with impact and friction[END_REF][START_REF] Luo | An unsymmetrical motion in a horizontal impact oscillator[END_REF][START_REF] Luo | Dynamic stability of weakly damped oscillators with elastic impacts and wear[END_REF][START_REF] Mikhlin | Direct and inverse problems encountered in vibro-impact oscillations of a discrete system[END_REF][START_REF] Murphy | Grazing instabilities and post-bifurcation behavior in an impacting string[END_REF][START_REF] Nordmark | Non-periodic motion caused by grazing incidence in an impact oscillator[END_REF][START_REF] Peterka | Transition to chaotic motion in mechanical systems with impacts[END_REF][START_REF] Wagg | Rising phenomena and the multi-sliding bifurcation in a two-degree of freedom impact oscillator[END_REF][START_REF] Yu | The periodic impact responses and stability of a human body in a vehicle traveling on rough terrain[END_REF]. However, these studies focused mainly attention on codim 1 bifurcations of the vibro-impact systems. The purpose of the present study is to focus attention on codim 2 bifurcation of single-impact periodic motion. The results from qualitative analysis and numerical simulation shows that the vibro-impact systems, near the points of codim 2 bifurcations, can exhibit richer and more complicated dynamical behavior.

Mechanical model and Poincaré map

The mechanical model for a vibro-impact system with masses M 1 and M 2 is shown in Fig. 1. The vibro-impact system is similar to the model of the impact-forming machine with double masses (see p. [START_REF] Luo | Dynamic stability of weakly damped oscillators with elastic impacts and wear[END_REF]. Displacements of the masses M 1 and M 2 are represented by X 1 and X 2 , respectively. The masses are connected to linear springs with stiffnesses K 1 , K 2 and K 3 , and linear viscous dashpots with damping constants C 1 , C 2 and C 3 . The excitations on both masses are harmonic with amplitudes P 1 and P 2 . Here is the excitation frequency, and the phase angle. The mass M 1 impact mutually with the mass M 2 when the difference of their displacements equals the clearance . The impact is described by the conservation law of momentum and a coefficient of restitution R, and it is assumed that the duration of impact is negligible compared to the period of the force.

The motion processes of the system, between any two consecutive impacts, are considered. Between any two consecutive impacts, the time T is always set to zero directly at the instant when the former impact is over, and the phase angle is used only to make a suitable choice for the origin of time in the calculation. The state of the vibro-impact system, immediately after impact, has become initial conditions in the subsequent process of the motion. Between consecutive impacts, the non-dimensional differential equations of motion are given by m 0 0 1

ẍ1 ẍ2 + 2 (1 + 1 ) -2 -2 2 (1 + 3 ) ẋ1 ẋ2 + k 1 + 1 -1 -1 k 3 + 1 x 1 x 2 = 1 -f 20 f 20 sin( t + ), (1)
where the non-dimensional quantities

m = M 1 M 2 , k i = K i K 2 , i = C i C 2 , i = 1, 3,
f 20 = P 2 P 1 + P 2 , = M 2 K 2 , t = T K 2 M 2 , = C 2 2 √ K 2 M 2 , = • K 2 P 1 + P 2 , x i = X i K 2 P 1 + P 2 (2)
have been introduced. In Eq. (1) a dot (•) denotes differentiation with respect to the non-dimensional time t. When the impact occurs, for x 1 (t)x 2 (t) = , the velocities of the masses M 1 and M 2 are changed according to the conservation law of momentum, and the impact equation of both masses and the coefficient of restitution r are given by

K 3 K 2 K 1 M 1 M 2 C 3 C 2 C 1 X 1 X 2 P 2 sin ( T + ) P 1 sin ( T + )
m ẋ1+ + ẋ2+ = m ẋ1-+ ẋ2-, R = ( ẋ2+ -ẋ1+ )/( ẋ1--ẋ2-), (3)
where the subscript minus sign denotes the states just before impact and the subscript plus sign denotes the states just after impact, ẋiand ẋi+ (i = 1, 2) represent the velocities of immediately before and after the impact of the masses M i , respectively. Impacting systems are conveniently studied by use of a mapping derived from the equations of motion. Each iterate of this map corresponds to the mass M 1 striking the mass M 2 once. Under suitable system parameter conditions, the system shown in Fig. 1 can exhibit periodic-impact behavior. We can characterize periodic motions of the system by a symbol q = p/n, where p is the number of impacts and n is the number of the forcing cycles. The symbol q, originated by Peterka, is of important significance in analyzing the periodic-impact motions and bifurcation characteristics, and quantity q has rational and irrational value for periodic and chaotic motions, respectively; see Ref. [START_REF] Peterka | Transition to chaotic motion in mechanical systems with impacts[END_REF]. Let = t, we can choose a Poincaré section = {(x 1 , ẋ1 , x 2 , ẋ2 , ) ∈ R 4 × S, x 1x 2 = , ẋ1 = ẋ1+ , ẋ2 = ẋ2+ } to establish the Poincaré map of the vibro-impact system. The q = 1 1 motion and its disturbed map are derived exactly in Ref. [START_REF] Luo | Periodic motions and bifurcations of vibro-impact systems[END_REF]. The disturbed map of q = 1 1 motion is expressed briefly by

X = f (v, X), ( 4
)
where

X ∈ R 4 , v ∈ R 1 or R 2 ; X = X * + X, X = X * + X , X * = ( ẋ1+ , x 20 , ẋ2+ , 0 )
T is a fixed point in the hyperplane , X and X are the disturbed vectors of X * . Linearizing the Poincaré map at the fixed point X * results in the matrix

Df (v, X * ) = jf (v, X) jX (v,X * )
.

The stability and local bifurcation of q = 1 1 motion can be analyzed by computing eigenvalues of the linearized matrix Df (v, X * ) of Poincaré map. If the moduli of all eigenvalues of Df (v, X * ) are less than one (i.e. | i (v)| < 1), then the corresponding q = 1 1 motion is stable. In contrast, if the modulus of one of the eigenvalues of Df (v, X * ) is greater than unity, the fixed point is unstable. In general, bifurcations occur in various ways according to the numbers of the eigenvalues on the unit circle and their positions on the unit circle. In the following sections, Hopf-flip bifurcation of the repeated impact system will be analyzed.

Center manifold and normal form map

We continue to consider the Poincaré map X = f (v, X). X * is a fixed point for the map for v in some neighborhood of a critical value v = v c at which Df (v, X * ) satisfies the following assumptions:

H1. Jacobian matrix Df (v, X * ) has the eigenvalues 1 (v c ) = -1, 2 (v c ) = ¯ 3 (v c ) and | 2,3 (v c )| = 1;
H2. the remainder of the spectrum of Df (v, X * ) are strictly inside the unit circle. For all v in some neighborhood of v c , the map (4), under the change of variables

1 = v 1 -v 1c , 2 = v 2 -v 2c , = ( 1 , 2 ) T and X = X * + P Ỹ , becomes Ỹ = F ( ; Ỹ ), ( 5
)
where P is the eigenmatrix [START_REF] Hassard | Theory and applications of Hopf bifurcation[END_REF], 5) is expressed by

Ỹ = (y 1 , y 2 , y 3 , y 4 ) T , F = (F 1 , F 2 , F 3 , F 4 ) T . Let z 1 = y 1 , z 2 = y 2 + iy 3 , z2 = y 2 -iy 3 , z = (z 1 , z 2 , z2 ) T , W = y 4 , G 1 = F 1 -1 z 1 , G 2 = F 2 + iF 3 -2 z 2 and H = F 4 -4 W , the map (
z 1 = 1 z 1 + G 1 (z 1 , z 2 , z2 , W ; ), z 2 = 2 z 2 + G 2 (z 1 , z 2 , z2 , W ; ), W = 4 W + H (z 1 , z 2 , z2 , W ; ), ( 6
) in which 1 = ˜ 1 ( ) = 1 (v c + ), 2 = ˜ 2 ( ) = 2 (v c + ), 1 (0) = -1, 2,3 (0) = ± i , | 2,3 (0)| = 1.
For the map (6), there exists a local center manifold W (z 1 , z 2 , z2 ; ), which can be determined by the equation

W (z 1 , z 2 , z 2 ; ) = 4 W (z 1 , z 2 , z2 ; ) + H (z 1 , z 2 , z2 , W (z 1 , z 2 , z2 ; ); ). ( 7
)
On the center manifold the local behavior of the map ( 7) can be reduced to a three-dimensional map ˜ (z; ), which is now

z 1 = 1 ( )z 1 + G 1 (z 1 , z 2 , z2 , W (z 1 , z 2 , z2 ; ); ), z 2 = 2 ( )z 2 + G 2 (z 1 , z 2 , z2 , W (z 1 , z 2 , z2 ; ); ), ( 8
)
where the Taylor series expansion of W (z 1 , z 2 , z2 ; ) and G i (z 1 , z 2 , z2 , W (z 1 , z 2 , z2 ; ); ) about (0, 0, 0, ) can be determined by a similar method used in Ref. [START_REF] Luo | Hopf bifurcation of a two-degree-of-freedom vibro-impact system[END_REF]. By using the center manifold theorem technique and normal form method of maps, we can reduce the map (8) to the normal form map (z; ), which is given by

z 1 = -z 1 + 1 z 1 + az 3 1 + bz 1 |z 2 | 2 + O(|z| 5 ), z 2 = 2 (0)z 2 + ˜ 2 z 2 + cz 2 1 z 2 + dz 2 |z 2 | 2 + O(|z| 5 ). ( 9
)
The normal form map (9) (Y ; ), in the real form, is expressed by

y 1 = -y 1 + 1 y 1 + ay 3 1 + by 1 (y 2 2 + y 2 3 ) + h.o.t y 2 = ( + 2 )y 2 -( + 3 )y 3 + (cy 2 -ey 3 )y 2 1 + (dy 2 -fy 3 )(y 2 2 + y 2 3 ) + h.o.t y 3 = ( + 3 )y 2 + ( + 2 )y 3 + (cy 3 + ey 2 )y 2 1 + (dy 3 + fy 2 )(y 2 2 + y 2 3 ) + h.o.t (10) in which, = ( 1 , 2 , 3 ) T , i = i ( ), i (0) = 0.

Local codim 2 bifurcation of the normal form map

Let us assume that there exist period two fixed points for the map (Y ; ). In view of the normal form map, the period two points Y (2) 1 and Y

(2)

2 satisfy the equation 2 (Y (2) i ; ) = Y (2) i . ( 11
)
Ignoring the terms of high order, the solutions of Eq. [START_REF] Leine | Nonlinear dynamics and modeling of various wooden toys with impact and friction[END_REF] are

Y (2) 1 = -1 /a, 0, 0 T , Y (2) 2 = --1 /a, 0, 0 T . ( 12
)
If 1 /a < 0, there exist the period two points Y (2)

= (Y (2) 1 , Y (2)
2 ) T , and they are symmetrical about the origin. The linearized maps of (Y ; ) at the fixed point and 2 (Y ; ) at the period two points, respectively, are given by

Q 1 = j (Y ; ) jY (0, ) , Q 2 = j 2 (Y ; ) jY (Y (2) i , ) . ( 13
)
The partial bifurcation sets for the normal form map can be determined by computing and analyzing the eigenvalues of Jacobian matrices [START_REF] Luo | Hopf bifurcation of a two-degree-of-freedom vibro-impact system[END_REF]. However, a full understanding of the normal form map [START_REF] Knudsen | Dynamic stability of weakly damped oscillators with elastic impacts and wear[END_REF] requires more than Jacobian matrices [START_REF] Luo | Hopf bifurcation of a two-degree-of-freedom vibro-impact system[END_REF]. So it is necessary to change the normal form map to the polar coordinate form (r, ; 0 ) ∈ R 2 × S, (x, r, ) → (x , r , ).

Let

˜ 2 = ˜ 2 (0)˜ 20 , c = ˜ 2 (0) c0 , d = ˜ 2 (0) d0 , the normal form map (9) becomes z 1 = -z 1 (1 -1 -az 2 1 -b|z 2 | 2 ) + O((|z 1 | + |z 2 |) 5 ), z 2 = ˜ 2 (0)z 2 (1 + ˜ 20 + c0 z 2 1 + d0 |z 2 | 2 ) + O((|z 1 | + |z 2 |) 5 ). ( 14
)
In polar coordinates, the map (z; 0 ) ∈ R 3 is changed to (x, r, ; 0 ) ∈ R 2 × S, which is given by

x = -x(1 -1 -ax 2 -br 2 ) + h.o.t, r = r(1 + 20 + c 0 x 2 + d 0 r 2 ) + h.o.t, = + 0 + 30 + e 0 x 2 + f 0 r 2 + h.o.t. ( 15
)
in which 20

= 2 + 3 , 30 = 3 -2 , c 0 = c + e, d 0 = d + f , e 0 = e -c, f 0 = f -d.
We ignore the influence of phase angle to the map [START_REF] Luo | Dynamic stability of weakly damped oscillators with elastic impacts and wear[END_REF] temporarily so that a simplified map is obtained, which is

x = -x(1 -1 -ax 2 -br 2 ) + h.o.t, r = r(1 + 20 + c 0 x 2 + d 0 r 2 ) + h.o.t. ( 16
)
Let the iterated map ˆ 2 (x, r; 0 ) = (x (2) , r (2) ; 0 ) T , then we can write

x (2) = x + 2(-x 1 -ax 3 -bxr 2 ) + h.o.t, r (2) = r + 2(r 20 + c 0 x 2 r + d 0 r 3 ) + h.o.t. ( 17
)
The map [START_REF] Mikhlin | Direct and inverse problems encountered in vibro-impact oscillations of a discrete system[END_REF] has two fixed points Y 0 = (0, 0) T and Y * 1 = (0, √ -20 /d 0 ) T , and the iterated map [START_REF] Murphy | Grazing instabilities and post-bifurcation behavior in an impacting string[END_REF] has two pairs of fixed points

Y * 21 = ( √ -1 /a, 0) T , Y * 22 = (- √ -1 /a, 0) T ; Y * 31 = b 20 -d 0 1 ad 0 -bc 0 , c 0 1 -a 20 ad 0 -bc 0 T , Y * 32 = - b 20 -d 0 1 ad 0 -bc 0 , c 0 1 -a 20 ad 0 -bc 0 T .
Three main cases associated with scheme of coefficients of high order terms, i.e.,

Case (I): a > 0, b > 0, c 0 < 0, d 0 < 0, ad 0 -bc 0 < 0, a + c 0 > 0 and b + d 0 > 0; Case (II): a > 0, b > 0, c 0 > 0, d 0 > 0, ad 0 -bc 0 < 0; Case (III): a < 0, b < 0, c 0 < 0, d 0 < 0, ad 0 -bc 0 < 0,
L 1 : 1 = 0, 20 < 0; L 2 : 20 = c 0 a 1 , 1 < 0 ; L 3 : 20 = d 0 (a+c 0 ) a(d 0 +b) 1 , 1 < 0; L 4 : 20 = 0, 1 > 0; L 5 : 20 = d 0 b 1 , 1 < 0.
The bounds of the regions, shown in Fig. 2(c), can be listed as follows

L 1 : 1 = 0, 20 < 0; L 2 : 20 = 0, 1 > 0; L 3 : 20 = d 0 b 1 , 1 > 0; L 4 : 20 = d 0 (a+c 0 ) a(d 0 +b) 1 , 1 > 0, L 5 : 20 = c 0 a 1 , 1 > 0.
Now we consider the influence of phase angle to the map [START_REF] Luo | Dynamic stability of weakly damped oscillators with elastic impacts and wear[END_REF]. By making a comparison between the map (16) and normal form map, we can unfold qualitative analyses for the normal form map. The fixed point Y * 1 of the map (16) corresponds to the invariant circle of the normal form map, which is associated with the fixed point Y (1) 0 = (0, 0, 0) T . The fixed points Y * 2i of the map ( 16) correspond to period two fixed points of the map [START_REF] Knudsen | Dynamic stability of weakly damped oscillators with elastic impacts and wear[END_REF]. The fixed points Y * 3i of the map ( 16) correspond to the invariant circles associated with the period two fixed points for the map [START_REF] Knudsen | Dynamic stability of weakly damped oscillators with elastic impacts and wear[END_REF]. For the normal form map in the simple Cases (I) and (II), Hopf and period-doubling bifurcations associated with the fixed point Y (1) 0 occur as parameters crossing the lines L 4 and L 1 from the region R 1 , respectively; crossing the line L 2 Hopf bifurcation associated with period two fixed points Y (2) i occurs. The Hopf bifurcation associated with the fixed point Y [START_REF] Aidanpää | Periodic and chaotic behaviour of a threshold-limited two-degree-of-freedom system[END_REF] 0 is supercritical for the Case (I), subcritical for the Case (II). For the normal form map in the complex Case (III), crossing the lines L 2 and L 1 Hopf and period-doubling bifurcations associated with the fixed point Y (1) 0 occur, respectively; as the parameters crossing the line L 3 from the region R 2 , the closed circle becomes non-attracting and the torus bifurcation occurs, which causes that the quasi-periodic attractor represented by two attracting closed circles is born. The quasi-periodic attractor represented by two closed circles is attracting in region R 3 , and non-attracting in region R 4 . It is to be noted that the period-doubling bifurcation associated with the fixed point Y [START_REF] Aidanpää | Periodic and chaotic behaviour of a threshold-limited two-degree-of-freedom system[END_REF] 0 is subcritical for the complex case.

In view of Jacobian matrixes (13), we can find that there exist the period two fixed points Y (2) = (Y ( 2)

1 , Y (2)
2 ) T as 1 /a < 0. According to the unfolding of the simple Case (I), the bifurcation set of the normal form map [START_REF] Knudsen | Dynamic stability of weakly damped oscillators with elastic impacts and wear[END_REF], near = (0, 0, 0) T , can be further illustrated by Fig. 2(d) in which < 0. The bounds of the regions shown in Fig. 2(d) can be listed as follow

L 11 : 1 = 0, 2 > 0; L 12 : 2 = 0, 1 > 0; L 13 : 2 = e + c a 1 , 1 < 0, L 14 : 2 = d + f b 1 , 1 < 0; L 15 : 2 = ( d + f )(a + c + e) a(b + d + f ) 1 , 1 < 0, (18)
The bifurcation sets of the normal form map [START_REF] Knudsen | Dynamic stability of weakly damped oscillators with elastic impacts and wear[END_REF], associated with the Cases (II) and (III), can be obtained by using a similar method.

According to the center manifold theory, local behavior of map f (v, X), near the bifurcation point v c , is equivalent to that of (Y ; ) for near = (0, 0, 0) T . By virtue of the analyses of local bifurcation of the normal form map [START_REF] Knudsen | Dynamic stability of weakly damped oscillators with elastic impacts and wear[END_REF], we can find out dynamical behavior of the vibro-impact system in the case of Hopf-flip bifurcation. Local behavior of map f (v, X), near the bifurcation point v c , conforms to the numerical results below.

Numerical analyses

In this section the analysis developed in the former section is verified by the presentation of results for the vibroimpact system. The existence and stability of q = 1 1 motion are analyzed explicitly. Also, local bifurcations at the points of change in stability, discussed in the previous section, are considered, thus giving some information on dynamical behavior near the point of Hopf-flip bifurcation. The system, with parameters: m = 2.0, k 1 In the second region of forcing frequency associated with q = 1 1 motion, the motion will undergo period doubling bifurcation with increase in , and q = 2 2 motion stabilizes. And then Hopf bifurcation of q = 2 2 motion occurs so that the system exhibits the quasi-periodic impact motion associated with q = 2 2 motion. The transition process of q = 1 1 motion is plotted locally in an amplified form in Fig. 4(a1). The results from simulation show that no period doubling cascade of q = 1 1 motion occurs, near the point of Hopf-flip bifurcation, due to occurrence of Hopf bifurcation of q = 2 2 motion. Moreover, subcritical Hopf bifurcation of q = 1 1 motion, in the second window, occurs with decrease in as seen in Figs. 4(a),(b) and (c), and the type of q = 1 1 fixed point, from stable focus to unstable focus, is changed. Some projected Poincaré sections are plotted in Figs. 567. A representative analysis is given by choosing the damping ratio = 0.01225 and changing the forcing frequency in numerical analysis. The numerical results show that the system exhibits stable q = 1 1 motion with ∈ (5.073489, 5.074559) as seen Fig. 3 and Fig. 4(a1). The q = 1 1 motion undergoes Hopf bifurcation as passes through c1 = 5.073489 decreasingly. The eigenvalues of Df ( , X * ) are given for c1 as follows

1,2 ( c1 ) = -0.3773860 ± 0.9260564i, | 1,2 ( c1 )| = 1.0000001, 3 ( c1 ) = -0.9876892, 4 ( c1 ) = -0.4073028.
We can conclude, by the result from simulation, that subcritical Hopf bifurcation of q = 1 1 motion occurs for < 5.073489, and the type of q = 1 1 fixed point, from stable focus to unstable focus, is changed; see Fig. 5 which shows transient points as well as stable q = 9 9 fixed points, starting from the initial condition near the unstable fixed point of q = 1 1 motion for = 5.0724. The Period doubling bifurcation of q = 1 1 motion occurs as is increased gradually and passes through c2 = 5.074559, and the system exhibits stable q = 2 2 motion; see Fig. 4(a1). The q = 2 2 motion changes its stability for > 5.07585, and Hopf bifurcation of q = 2 2 motion occurs so that the system exhibits quasi-periodic impact motion associated with q = 2 2 points as seen in Figs. 6(a) and (b). With further increase in , the closed circle becomes quasi-attracting. The quasi-attracting invariant circle is attracting for the map point inside The q = 1 1 motion, associated with the first window, undergoes Hopf bifurcation as passes through c3 = 4.74017 increasingly. With increase in the forcing frequency , the quasi-periodic attractor transits to chaos via phase locking; see Fig. 7.

The dynamic behavior near the point of Hopf-flip bifurcation, occurring in the system of Fig. 1, corresponds with the unfolding of the simple Case (II) shown in Fig. 2(b).

In order to observe Hopf-flip bifurcations and transition phenomena of periodic-impact motions associated the complex case, we consider a vibro-bouncing system shown in Fig. 8. A body with mass m bounces on the flat horizontal surface of a vibro-bench with masses M 1 and M 2 . Displacements of these masses m, M 1 and M 2 are represented by Y 0 , X 1 and X 2 , respectively. The masses M 1 and M 2 are connected to linear springs with stiffnesses K 1 and K 2 , and linear viscous dashpots with damping constants C 1 and C 2 . The excitations on both masses of the vibro-bench are harmonic with amplitudes P 1 and P 2 . The mass M 1 impacts mutually with the bouncing mass m when they are on the same height, so the mass m exhibits the bouncing motion.

The motion processes of the system, between any two consecutive impacts, are considered. Between consecutive impacts, the non-dimensional differential equations of motion are given by 1 0 0

m ẍ1 ẍ2 + 2 -2 -2 2 (1 + c ) ẋ1 ẋ2 + 1 -1 -1 1 + k x 1 x 2 = 1 -f 20 f 20 sin( t + ), ( 19
) ÿ = -, (20)
where the non-dimensional quantities are given by

m = M 2 M 1 , k = K 2 K 1 , c = C 2 C 1 , f 20 = P 2 P 1 + P 2 , = M 1 K 1 , t = T K 1 M 1 , Y 0 M 1 M 2 C 2 C 1 P 1 sin ( T + ) P 2 sin ( T + ) K 1 K 2 m X 1 X 2
Fig. 8. Schematic of the vibro-impact system.

= C 1 2 √ K 1 M 1 , x i = X i K 1 P 1 + P 2 , = M 1 g P 1 + P 2 , = m M 1 , y = Y 0 K 1 P 1 + P 2 . ( 21
)
When the impact occurs, for x 1 (t)=y(t), the velocities of masses m and M 1 are changed according to the conservation law of momentum, and the impact equation of both masses m and M 1 and the coefficient of restitution R are given by ẋ1-+ ẏ-= ẋ1+ + ẏ+ , ẋ1+ -ẏ+ = -R( ẋ1--ẏ-).

(

) 22
where the subscript minus sign denotes the states just before impact and the subscript plus sign denotes the states just after impact, the velocities of two masses M 1 and m, just immediately before and after the impact, are represented respectively by ẋ1-, ẏ-, ẋ1+ and ẏ+ . A Poincaré section = {(x 1 , ẋ1 , x 2 , ẋ2 , y, ẏ, ) ∈ R 6 × S, x 1 = y, ẋ1 = ẋ1+ , ẏ = ẏ+ } is chosen to establish Poincaré map of the impact system, which can be derived analytically from the equations of motion. The disturbed map of period one single-impact motion is expressed by

X = f (v, X), ( 23
)
where

X ∈ R 6 , v ∈ R 1 or R 2 is varying parameter; X=X * + X, X =X * + X , X=( ẋ1+ , x 1 , ẋ2 , x 2 , ẏ+ , ) T and X = ( ẋ 1+ , x 1 , ẋ 2 , x 2 , ẏ + ,
) T are the disturbed vectors of X * , X * = ( ẋ1+ , x 10 , ẋ2 , x 20 , ẏ+ , 0 ) T is a fixed point in the hyperplane .

An interest torus doubling bifurcation is found to exist in the vibro-impact system near the point of Hopf-flip bifurcation associated with the complex case. The torus doubling bifurcation makes the quasi-periodic attractor associated with q = 1 1 motion transit to the other quasi-periodic attractor represented by two attracting closed circles. The phenomenon concerning torus doubling bifurcation is demonstrated by the example in the following text.

The vibro-bounce system, with parameters: = 0.3, k = 1, c = 1, = 0.5, f 20 = 0, = 0.05 and R = 0.6, has been chosen for analysis. The forcing frequency and distribution of masses m are taken as the control parameters, i.e. v = ( m , ) T . The eigenvalues of Jacobian matrix D f (v, X * ) are computed. All eigenvalues of D f (v, X * ) stay inside the unit circle for v = (0.72, 0.5275) T . By gradually decreasing m and from the point v = (0.72, 0.5275) T to change the control parameters v, we obtain the critical value v c = (0.7084388, 0.5275136) T , at which all eigenvalues

Conclusions

Dynamical behavior of the vibro-impact systems, near the point of Hopf-flip bifurcation, is investigated by qualitative analyses and numerical simulation. The vibro-impact systems, under the condition of codim 2 bifurcation, can exhibit more complicated quasi-periodic impact motions than those which occur in non-resonance and weak resonance cases [START_REF] Luo | Hopf bifurcation of a two-degree-of-freedom vibro-impact system[END_REF]. Near the point of Hopf-flip bifurcation there exists not only Hopf bifurcation of q = 1 1 motion, but also Hopf bifurcation of q = 2 2 motion or torus bifurcation. The period doubling bifurcation of q = 1 1 motion is commonly existent for the simple cases near the point of Hopf-flip bifurcation. However, no period doubling cascade of the motion emerges due to occurrence of Hopf bifurcation of q = 2 2 motion. The q = 2 2 fixed points are symmetrical about the corresponding q = 1 1 point. The results from simulation show also that the attracting invariant circles associated with q = 2 2 points are smooth in nature and symmetrical about the corresponding q = 1 1 point near the value of Hopf bifurcation of q = 2 2 fixed points. As the value of moves further away from the value of Hopf bifurcation, two attracting invariant circles expand, and the smoothness and symmetry of the quasi-periodic attractor are changed by degrees until they are destroyed.

An interest torus doubling bifurcation is found to exist near the value of Hopf-flip bifurcation. The torus doubling bifurcation makes the quasi-periodic attractor associated with period one single-impact motion transit to the other quasi-periodic attractor represented by two attracting closed circles. The torus bifurcation is qualitatively different from the typical torus doubling bifurcation occurring in the vibro-impact systems.

The strict condition of codimension two bifurcation is not easy to encounter in practical application of engineering. However, there exist the possibilities that actual nonlinear dynamical systems, with two varying parameters or more, work near the critical value of codim 2 bifurcation due to change of parameters. The impact-forming machinery is a typical example. Besides the forcing frequency , the clearance varies also with different thickness of the formed workpieces [START_REF] Wen | Theory and Application of Vibratory Mechanism[END_REF]. Another representative example is the inertial shaker, of which the distribution of masses is generally metabolic with the casts with different masses, and the forcing frequency is also important parameter changed [START_REF] Wen | Theory and Application of Nonlinear Oscillation[END_REF]. The change of multi-parameters possibly leads to the results that the vibro-impact systems work near the critical parameters of codim 2 bifurcation. It is necessary to study the bifurcations caused by change of multi-parameters and dynamical behavior of nonlinear systems near the points of bifurcations.

Fig. 1 .

 1 Fig. 1. Schematic of the impact oscillator with double masses and a clearance.

  are considered, respectively. The Cases (I) and (II) are called the simple cases; the Case (III) the complex case.The bifurcation set for the map (16), near = ( 1 , 20 ) T = (0, 0) T , can be illustrated by Figs. 2(a)-(c). Only the positive (x, r) quadrant is shown in Figs. 2(a)-(c). Since the portraits are symmetric under reflection about r axis. The bounds of the regions, shown in Figs. 2(a) and (b), can be listed as follows:

Fig. 2 .

 2 Fig. 2. The partial bifurcation sets of the simplified map and normal form map: (a) Case (I), the simplified map; (b) Case (II), the simplified map; (c) Case (III), the simplified map; (d) Case (I), the normal form map.

  = 4.0, k 3 = 3.0, = 0.02, f 20 = 0 and R = 0.7, has been chosen for analysis. The forcing frequency and damping ratio are taken as the control parameters, i.e. v = ( , ) T . The eigenvalues of Df (v, X * ) are computed with ∈ [5.025, 5.2] and ∈ [0.01, 0.02]. All eigenvalues of Df (v, X * ) stay inside the unit circle for v = (0.02, 5.025) T . By gradually increasing and decreasing from the point v = (0.02, 5.025) T to change the control parameter v, we can obtain the eigenvalues 1 (v c ) = -1.0000002 and 2,3 (v c ) = -0.3768488 ± 0.9262742i (| 2,3 (v c )| = 0.9999994) which are very close to the unit circle, and the fourth eigenvalue ( 4 (v c ) = -0.4031064) still stays inside the unit circle as v equals v c = (0.0121275, 5.074519) T . The eigenvalues 1 (v) and 2,3 (v) have scaped the unit circle as (increasingly) and (decreasingly) pass through v = (0.0121274, 5.07452) T . The eigenvalues 1 (v)and 2,3 (v) almost escape the unit circle simultaneously, so v c = (0.0121275, 5.074519) T may be approximately taken as the value of Hopf-flip bifurcation.Local behavior of the vibro-impact system, near the point of Hopf-flip bifurcation, is obtained by numerical simulation. The partial bifurcation set near the critical value is plotted in Fig.3. The bifurcation diagrams are shown for a series of values of in Fig.4. There exist two windows of q = 1 1 motion in Figs. 4(a), (b) and (c), and only the second one occurs near the value of Hopf-flip bifurcation. It should be noted that the second region of forcing frequency associated with q = 1 1 motion is narrow, and vanishes for < 0.0121275 (see Figs. 4(d) and (e)). With increase in damping ratio ( > 0.0121275), the second window of q = 1 1 motion becomes wide gradually; see Figs. 4(a)-(c). The q = 1 1 motion, associated with the first window, will undergo Hopf bifurcation with increase in ; see Fig. 4.

Fig. 3 .

 3 Fig. 3. The partial bifurcation set near the point of codim 2 bifurcation.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Bifurcation diagrams near the point of codim 2 bifurcation: (a) = 0.01225; (a1) = 0.01225; (b) = 0.0135; (c) = 0.0175; (d) = 0.012; (e) = 0.007.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Projected Poincaré sections: (a) attracting invariant circles of q = 2 2 points, = 5.076, = 0.01225; (b) attracting invariant circles of q = 2 2

Acknowledgements

The author gratefully acknowledge the support by National Natural Science Foundation (10572055, 50475109), Natural Science Foundation of Gansu Province Government of China (ZS-031-A25-007-Z (key item)) and 'Qing Lan' Talent Engineering by Lanzhou Jiaotong University.

The partial bifurcation set near the critical value is plotted in Fig. 9. Numerical results are shown for m =0.70845 and a series of values of in Fig. 10. The process of torus bifurcation can be observed clearly from the projected portraits of Poincaré map. As the forcing frequency passes through the critical value c1 = 0.527511 decreasingly, a complex conjugate pair eigenvalues of D f ( , X * ) escape the unit circle, and supercritical Hopf bifurcation associated with q = 1 1 motion occurs. The quasi-periodic attractor is shown for = 0.5275 in Fig. 10(a). With decrease in , instability of the closed circle associated with q = 1 1 point occurs, and the torus bifurcation leads to two attracting invariant circles; see Fig. 10(b). With further decrease in , the two invariant circles become non-attracting and tori doubling are born, which subsequently becomes unstable and chaotic; see Figs. 10(c) and(d). As the forcing frequency passes through c2 = 0.5275321 increasingly, a real eigenvalue of Jacobian matrix D f ( , X * ) escapes the unit circle from the point (-1, 0), and period doubling bifurcation associated with q = 1 1 motion occurs. However, the period doubling bifurcation is subcritical. Stable q = 1 1 motion and unstable q = 2 2 motion coexist in the region ∈ (0.527511, 0.5275321). After passes through c2 = 0.5275321 increasingly, the q = 1 1 motion becomes unstable and there exists no q = 2 2 motion. The dynamic behavior near the point of Hopf-flip bifurcation, occurring in the example of the system shown in Fig. 8, corresponds with the unfolding of the complex Case (III) shown in Fig. 2(c).


Roman Iakymchuk
email: riakymch@pdc.kth.se

Stepan Shakhno
email: sshakhno@lnu.edu.ua

Halyna Yarmola
email: hyarmola@lnu.edu.ua

Combined Newton-Kurchatov method for solving nonlinear operator equations

We investigate local and semi-local convergence of the combined Newton-Kurchatov method under the classical and generalized Lipschitz conditions for solving nonlinear equations. The convergence order of the method is examined and the uniqueness ball for the solution of the nonlinear equation is proved. Numerical experiments are conducted on test problems.

Introduction

Let us consider the nonlinear equation

H(x) ≡ F (x) + G(x) = 0, (1)
where F and G are nonlinear operators defined on a convex subset D of a Banach space X with their values in a Banach space Y . F is a Fréchet-differentiable operator and G is a continuous operator; the differentiability of G is not required.

For solving the problem (1), we proposed in [2] the two-point iterative process

x n+1 = x n -[F (x n ) + G(2x n -x n-1 ; x n-1 )] -1 H(x n ), n = 0, 1, . . . , (2)
which is built on the Newton method and the Kurchatov method of linear interpolation [1]. We studied the semi-local convergence of the method (2) and showed that its convergence order is (1 + √ 5)/2 ≈ 1.618 . . . , which is lower than the convergence order of the underlying methods. In this work, we extend our studies to the local convergence of the method (2) under the classical and generalized Lipschitz conditions for the first-order derivatives of the operator F and for the first-and second-order divided differences of the operator G. Moreover, we investigate the semi-local convergence of the method (2) under the classical Lipschitz conditions and determine its radius of the convergence ball and the quadratic convergence order.

Local convergence analysis under the generalized Lipschitz conditions

Let us denote B(x, R) = {x : x -x < R} an open ball of radius R with a center at the point x ∈ D.

Theorem 2.1 Let F, G : D ⊆ X → Y be continuous nonlinear operators, where X and Y are Banach spaces. Suppose, that: 1) H(x) = 0 has a solution x * ∈ D and there exists a Fréchet derivative H (x * ) that is invertible; 2) F has the derivative of the first-order and G has divided differences of the first-and second-order on B(x * , 3r) ⊂ D, satisfying the generalized Lipschitz conditions on B(x * , 3r)

H (x * ) -1 (F (x) -F (x τ )) ≤ ρ(x) τ ρ(x) L 1 (u)du, H (x * ) -1 (G(x; y) -G(u; v)) ≤ x-u + y-v 0 L 2 (z)dz, H (x * ) -1 (G(u; x; y) -G(v; x; y)) ≤ u-v 0 N (z)dz, where 0 ≤ τ ≤ 1; x, y, u, v ∈ B(x * , 3r); x τ = x * + τ (x -x * ); (x) = x -x * ; L 1 , L 2 , N are positive nondecreasing functions; 4) r > 0 satisfies the equation 1 r r 0 L 1 (u)udu+ r 0 L 2 (u)du+2r 2r 0 N (u)du / 1- r 0 L 1 (u)du+ 2r 0 L(u)du+2r 2r 0 N (u)du = 1.
Then for all x -1 , x 0 ∈ B(x * , r) the iterative process (2) is correctly defined and the generated sequence {x n } n≥0 , which belongs to B(x * , r), converges to x * and satisfies the inequality

x n+1 -x * ≤ 1 ρ(xn) ρ(xn) 0 L 1 (u)udu + ρ(xn) 0 L 2 (u)du + xn-xn-1 0 N (u)du x n -x n-1 1 - ρ(xn) 0 L 1 (u)du + 2ρ(xn) 0 L 2 (u)du + xn-xn-1 0 N (u)du x n -x n-1 x n -x * .
Corollary 2.2 The Newton-Kurchatov method (2) has the quadratic convergence order.

3 Semi-local convergence analysis under the classical Lipschitz conditions Theorem 3.1 Let F, G : D ⊆ X → Y be continuous nonlinear operators, where X and Y are Banach spaces. Suppose that: 1) A 0 = F (x 0 ) + G(2x 0 -x -1 ; x -1 ), where x -1 , x 0 ∈ U 0 = {x : x -x 0 ≤ r 0 } ⊂ D, is invertible operator; 2) F has the derivative of the first-order and G has divided differences of the first-and second-order on V 0 = {x : x -x 0 ≤ 3r 0 } ⊂ D, satisfying the Lipschitz conditions on V 0

A -1 0 (F (x) -F (y)) ≤ 2l 0 x -y , A -1 0 (G(x; y) -G(u; v)) ≤ p 0 ( x -u + y -v ), A -1 0 (G(x; y; z) -G(u; y; z)) ≤ q 0 x -u ;
3) a, c, r 0 are non-negative numbers such that

x 0 -x -1 ≤ a, A -1 0 H(x 0 ) ≤ c, c < a, r 0 < 1 -2q 0 a 2 -(l 0 + p 0 )c 2(l 0 + p 0 ) , r 0 ≥ c 1 -γ , γ = (l 0 + p 0 )c + q 0 a 2 1 -q 0 a 2 -2(l 0 + p 0 )r 0 , 0 ≤ γ < 1.
Then the following inequalities hold for all n ≥ 0

x n -x n+1 ≤ t n -t n+1 , x n -x * ≤ t n -t * , where t -1 = r 0 + a, t 0 = r 0 , t 1 = r 0 -c, t n+1 -t n+2 = (l 0 + p 0 )(t n -t n+1 ) + q 0 (t n-1 -t n ) 2 1 -q 0 a 2 -2(l 0 + p 0 )(t 0 -t n+1 ) (t n -t n+1 ), n ≥ 0,
{t n } n≥0 is a non-negative, nonincreasing sequence that converges to a certain t * such that r 0 -c/(1 -γ) ≤ t * < t 0 ; the sequence {x n } n≥0 , generated by the iterative process (2), is well defined, remains in U 0 and converges to the solution x * . Theorem 3.2 Let F (x * ) + G(x * ) = 0. Suppose that: 1)

A 0 = F (x 0 ) + G(2x 0 -x -1 ; x -1 )
, where x -1 , x 0 ∈ U 0 , is invertible operator; 2) F has the derivative of the first-order and G has divided differences of the first-order in D, satisfying the Lipschitz conditions

A -1 0 (F (x) -F (y)) ≤ 2l 0 x -y , A -1 0 (G(x; y) -G(u; v)) ≤ p 0 ( x -u + y -v ),
where x, y, u, v ∈ D; 3) r 0 satisfies the inequality r 0 < (1 -2p 0 a)/(2(l 0 + p 0 )).

Then the problem (1) has the unique solution x * ∈ U 0 . If there exists r 1 > r 0 , that U 1 = {x : x -x 0 ≤ r 1 } ⊂ D and 2p 0 a + (l 0 + p 0 )(r 0 + r 1 ) < 1, then the problem (1) has the unique solution x * ∈ U 1 .

Numerical experiments

We carried out a set of experiments and present in Tab. 1 the results for the following system of nonlinear equations:

3x 2 y -y 2 -1 + |x -1| = 0, x 4 + xy 3 -1 + |y| = 0.
We used initial approximations x 0 = (d, 0), x -1 = (5d, 5d) and stopping conditions x n+1 -x n ∞ ≤ , H(x n+1 ) ∞ ≤ . The obtained results confirm that the Newton-Kurchatov method (2) has the highest convergence speed.

Table 1 :

 1 The amount of iterations required to find the approximation to the solution of the problem (1) with the accuracy = 10 -15 . d Newton-Kurchatov method (2) Kurchatov method[1] Newton-type method [


Clifford A Reiter
email: reiterc@lafcol.lafayette.edu

CHAOTIC ATTRACTORS WITH THE SYMMETRY OF A TETRAHEDRON

Functions that are equivariant with respect to the symmetries of a tetrahedron are determined. Linear combinations of these functions that give rise to chaotic attractors are used to create images in 3space of attractors with the symmetry of the tetrahedron. These a1 tractors are visually appealing because of the tension between the pattern forced by their symmetry and the randomness arising from their chaotic behavior.

INTRODUCTION

Chaotic attractors arising from the iterative solution to systems of differential equations have been the subject of study since Lorenz noticed strange behavior in a model related to weather prediction [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF]. Even the iteration of simple functions such as the logistic function can be used to illustrate dynamics including chaotic behavior; see, for example Ref. [START_REF] Peitgen | Chaos and Fractals[END_REF]. More recently, attractors arising from the iteration of functions in the plane have been studied and equivariant functions have been used to create attractors that appear chaotic while having rotational and/or reflectional symmetries [START_REF] Field | Symmetric chaos[END_REF][START_REF] Fi | Symmetry in Chaos[END_REF][START_REF] Fi~1d | Symmetric chaos: how and why[END_REF] . Examples of attractors with the symmetry of the cube in 3space and n-cube in n-space have also been created [START_REF] Brisson | Symmetric attractors in three-dimensional space[END_REF] [START_REF] Reiter | Attractors with the symmetry of then-cube[END_REF]; in addition, examples of attractors in 3-space with dueling planar symmetries have been studied [8f In this paper we create attractors in 3-space which empirically appear to be chaotic but which have the symmetry of the tetrahedron.

Discussions of the symmetry groups of various geometric shapes can be found in a variety of places. In particular, Coxeter [START_REF] Ccxeter | Regular Polytopes[END_REF] describes the symmetry groups of the cube and tetrahedron. There are twelve rotations that preserve the four vertices of a regular tetrahedron. These include two third-turns that fix any one vertex of the tetrahedron yielding a total of eight third-turns. See Fig. 1. Also, there are three half-turns about the line connecting midpoints of opposing edges that switch two pairs of vertices. Lastly, the identity function also preserves the tetrahedron giving the twelve desired rotations. The results of composing these rotations can readily be observed geometrically. If reflections are also allowed, then the size of the symmetry group is twice as much: 24. While Fig. 1 shows the tetrahedron oriented with a vertex at a peak, it will be more convenient for us to consider the tetrahedron embedded in the cube by selecting alternate vertices of the cube. See Fig. 2. This representation makes it clear that the rotational symmetry group of the tetrahedron. is a subgroup of the rotational symmetry group of the cube. Notice the rotational symmetry group of the cube has quarter-turns (and some halfturns) that do not preserve the tetrahedron.

The next section develops the mathematics we need in order to construct our chaotic attractors with the symmetry of the tetrahedron. In particular, functions equivariant with respect to those rotations are determined. The last section gives illustrations of such attractors selected for visual appeal and diversity and describes our computations.

RNCTIONS WITH THE SYMMETRY OF A TETRAHEDRON

Let T denote the tetrahedron with vertices (1, I , I), (1,-1,-1), (-1,1,-1) and (-1,-1,1). One can check by direct computation that all twelve of the rotational symmetries of the tetrahedron are generated by composition from the two rotations rr(x,y , z)(--x , -y, z) and r(x, y, .:) = (z, x , y) . Notice that r is a third-turn which fixes the first vertex of T. The convenience of our choice of T is given by the fact that we can describe this third turn on T without adding coordinates or introducing a factor of J3 into the description of the rotations. The rotation () is a half-turn that interchanges two pairs of vertices of T. Moreover, all the reflectional and rotational symmetries of the tetrahedron are generated by () and r, along with the reflection v where v(x, y , .:) = (y , x , z) as can be verified by checking these generate the required 24 functions.

In order to generate attractors with the symmetry of the tetrahedron, we need functions from R 3 to R 3 that preserve those rotations and reflections in a certain manner. In particular, a function f is said to be equivar iant with respect to the group of symmetries of T if for all a in that group of symmetries and point are the same as the rotation of the iterates of the point. This means that the attractor associated with f tends to have the desired symmetries; however, it is possible that the attractor has only the symmetry given by an admissible subgroup of symmetries [START_REF] Ashwin | Symmetry groups of attractors[END_REF]II] or it may have the desired symmetry only in a trivial manner. Thus, even after identifying the appropriate equivariant functions one still needs to do some work to find examples that highlight the chaos and symmetry.

X E R 3 ,f(c((X)) = a(f(X)). If oc is
It is straightforward to check that a function that is equivariant with respect to the generators of a group will be equivariant with respect to every element of the group. With that observation, we are ready to determine the polynomial functions equivariant with respect to the symmetries of T.

- ---r-------------------------~-----~------------ where j and k have the same parity (even or odd) which is different from the parity of i.

(ii) The polynomial functions that are equivariant with respect to the rotational and reflectional symmetries c{ the tetrahedron T are linear combinations of functions of the form TiJk(x,y, :) + T;kj(X,y, z) where j and k have the same parity which is different from the parity of i .

Proof. A general polynomial function, P, from R 3 to R 3 has the form where N is the maximal degree of any coordinate that appears. P is equivariant with respect to the rotational symmetries of the tetrahedron if it is equivariant with respect to a and r. In particular, we must have r(P < x,y,z >) = P(r < x,y,z > ). Now each coordinate (but the exponents will be rearranged). Thus, each polynomial function equivariant with respect to r is a linear combination of functions of the form \ """"'

r(P(x, y , z)) = / L CijkXi),
i+J i J k P(a(x,y,z)) = 6 (-!) a!;kX) z , O~ij,k ~N """"' ( l )i+jb ;v'• ~k 6 - !ikX . ~ , O ~iJ,k ~N
which implies that a nonzero term aiJkxiyi~ can occur only when the parity of i and j are different. Moreover, such a nonzero term corresponds to C Jkixi / ::i in the third coordinate and a-equivariance implies j and k have the same parity. Thus, the only way for a function to be r and a-equivariant is if it is a linear combination of TiJk's where j and k have the same parity which is different from the parity of i.

Now consider the condition implied by the third coordinate of the equality: v(P(x,y,z)) = P( v(x.y,.:) ). We see that = from which we see

. . ..k c iJk .v' :..-!"" i . k Cj ikX /z c!ik = CJik and hence CJki = CkJi
Now this means that any time there is a nonzero term in the third coordinate there are two such terms Gust one if j = k) and if we have equivariance with respect to the rotations then there are two corresponding terms in each other coordinate. In particular, we see that a function equivariant with respect to a, r, and v must be a linear combination of functions of the form T iJk (x , y , z) + T ikJ(x, y, ::) Moreover, we can directly check that such func-Fig. 5. A cube wrapping attractor that has only the symmetries of the tetrahedron.

tions are equivariant with respect to (J, r, and v giving the desired result. 0

For example, we see the function f(x,y,z) = JT324(x,y, z) + 5Ttoz (x, y,.:) is equivariant with respect to the rotational symmetries of T. The function

f( x, y , z) = ST324(x,y, z) + ST342(x,y, z)
is equivariant with respect to the rotational and reflectional symmetries of T. Also, anytime i and j have different parities and we take j = k , it follows that

T .. •(x ,, ' 7) -lT .. •(x )' -) + lT .. •(x v ..,) •U , ' .r, --2 I}} • ' ' "- 2 l)J ' . , ~
is equivariant with respect to the rotational and reflectional symmetries of T.

ILLUSTRATIONS AND COMPUTATIONS

We noted above that linear combinations of the functions TiJk for the appropriate choice of indices gives functions that are equivariant with respect to the symmetry group of T. We created linear combinations of those functions in the programming language .I which is available from http://www.jsoftware.com/. The indices and coefficients were selected at random and the resulting function was tested to determine whether it created a nontrivial attractor (however, we always include the indices 1 0 0 since we want the fixed point at 0 0 0 to be repelling). If the attractor was nontrivial, then an image was created and observed. Images that appeared promising were mutated (:parameters varied along a search direction [START_REF] Reiter | Attractors with dueling symmetry[END_REF]) so that variations on the promising attractors could be observed. Finally, higher quality images were created for several functions using parallel computations of the attractors [START_REF] Brisson | Parallel processing in[END_REF]. While the details of J code often cannot be read by the uninitiated, we offer. in Table l, three lines of the details so that readers interested in replicating our work have Fig. 6. An attractor with the symmetries of the tetrahedron and platelike components.

sufficient information to implement our functions in J for comparison or experimentation. The first line constructs the function builder (conjunction) T; the second line creates a sample function f from the previous section and the third line shows how to create some iterates of that function.

Figure 3 shows the result of 200 million iterations of the function defined by f = :

1. 66332962 T 1 0 0 + 0. 419143344 T 2 55 + 0. 5422 2 1824 T 3 0 0 + 1. 00 T 0 1 1

Notice that the underlying shape of a tetrahedron is apparent, there are triple crossovers near the vertices and there are some hot spots nearby that are connected to each other in a complicated fashion . There are also hot spots near the edges that are interesting and this attractor also has reflectional symmetry. The color scheme shows pixels that are visited a small number of times in red; as the frequency of visits increases the colors run through the hues of the spectrum with magenta corresponding to the highest frequencies. Figure 4 also shows an attractor where a tetrahedron is also apparent. However, in this case notice that opposite each vertex there is a swirling ear/neck which lack reflectional symmetry and the attractor seems to fill space much more than the attractor in Fig. 3. The function used for Fig. 4 is given by the following. which appears to roughly form three girdles about a cube. However, notice that the rotations around the central axes perpendicular to any one of the girdles require a half-turn to preserve the attractor. Quarterturns cause the oscillations to be misoriented. The point being that although there is a cubical underlying form, the symmetries of the attractor are the symmetries of the tetrahedron and not the symmetries of the cube. Figure 6 shows the attractor produced by the function This figure is qualitatively quite different from the others since it seems to be constructed as a composite of many plates. This attractor also has reflectional symmetry. Figure 7 shows the attractor produced by the function Here the vertices of the tetrahedron lie along an arc connecting two points and the connections are reminiscent of wings. While it is difficult to prove results about the chaotic behavior of attractors, we can investigate this empirically via Ljapunov exponents [START_REF] Peitgen | Chaos and Fractals[END_REF][START_REF] Parker | Practical Numerical Algorithms for Chaotic Systems[END_REF]. Table 2 shows the Ljapunov exponents for the attractors shown in our images. Each of these has at least one positive Ljapunov exponent: this is associated with chaotic behavior. It is interesting to note that Figs 4 and (i have two positive Ljapunov exponents and these appear to be the heaviest attractors. Of course, these computations are not definitive. For example, the function is closely related to the function used to produce Fig. 3; however, our experience with this function shows the iterates of points converging to a 126-cycle after tens of thousands of iterations. One cannot exclude the possibility that our examples do not likewise become trivial after some very large number of iterations. We have seen that by choosing a careful representation of a regular tetrahedron we could describe generators for the rotational and reflectional symmetries in terms of sign changes and coordinate rearrangements. This allowed us to describe the functions equivariant with respect these symmetries as linear combinations of some simple generating functions. Random linear combinations of these functions could then be investigated yielding examples of chaotic attractors with the desired symmetries. Often the tetrahedral symmetry leads to attractors in the form of a tetrahedron but the attractor can be in the shape of a cube with some oriented features that destroy the cubical symmetry and even more bizarre shapes can appear.

•

  Fig. I. A tetrahedron and axes of rotation for a third-turn and a half-turn.

Fig. 2 .

 2 Fig.2. A regular tetrahedron embedded in a cube.

Fig. 3 .

 3 Fig.3. An attractor with tetrahedral shape and symmetry.

Fig. 4 .

 4 Fig. 4. An attractor with tetrahedral symmetry and opposing ears.

Figure 5

 5 Figure 5 shows the attractor produced by the function f = : 1.6724 51176 T 10 0 + 1.55056093 2 T 14 4 + __ 0.340028192 T 14 0 + _ 0.416 2611 08 T 3 0 0

2 Fig. 7 .

 27 Fig. 7. An attractor with the symmetries of the tetrahedron and wings.

  Fig. 1 2 3 3 0.111 -0 . 030 -0.271 4 0 .447 0 . 256 -0.072 5 0.22 8 -0 . 297 -1.315 6 0.556 0 . 555 -0.003 7 0.238 0.114 -0. 227

  Jzk , L aukxiyi:Jc ,

	0.5	0.4	0.5
	0.62875	0.503	0.4024
	0.514	1. 00391	0.514
	0.70765	1. 38212	2.69946
			\o~ij.k~N		O~ij,k, ; ; . N
			O ,;;. i,j.k ~N
				"""	)	k	"""	i . k
			(	O~ij,k.;,N		O.;,ij.k.;,N
				""" . k)

while " .. k) L.,; c iik z' xl y = ~ akiJ:X. y' z , L....., bkiJX y' z , L.,. CkijX 1 y'Z O.;, ij,k .;,N where the second equality holds by renaming the indices. Equating like coefficients yields Thus any nonzero term in any coordinate will be associated with a term with the same coefficient in Table l. Details of J code T=:2: 'm.&*@ ( (*/@: •"1)&( (-i.3) 1 ."O_n.))' f=: 3 T 3 2 4+5 T 1 0 2 f": (i.4) 0.5 0.4 0.5

Table 2 .

 2 Lj apunov exponents for the symmetric attractors

I.