Coupled chaotic elements, extended in space, often display complex time evolutions whose description cannot be captured by a low-dimensional dynamical model. Spatially distributed coupled chaotic elements was proposed as a simple model for high dimensional chaos involving spatial pattern dynamics [1]. The coupled map lattice (CML) approach can be viewed as a combination of two processes: collection of elements that locally exhibit chaotic dynamics on a lattice, and a coupling between these elements leading to a diffusive process. CML has predicted various classes of qualitative behavior. For example, spatial bifurcation  [1], frozen random chaos [2], pattern selection with suppression of chaos [3], spatiotemporal intermittency [4], soliton turbulence [2], and quasistationary supertransients [5 ]. These phenomena have found experimental verification in many areas, like Benard convection, chemical reactions, and biological media.

Umberger et al. [6] studied a 1-D chain of forced ideal Duffing oscillators, coupled via a linear dispersive term. They observed many classes of qualitative behavior found in CMLs. They argue that the dispersively coupled lattice ODE system might be expected to model typical characteristic phenomena to be expected in a forced, spatially extended nonlinear wave system with dispersion and damping.

Previous studies of the nonlinear dynamics of impact oscillators have generally involved, either literally or effectively, a 3-dimensional state space. In an attempt to get a qualitative understanding of the spatial complexity in the dynamics of solids, structures, and multimesh gear systems which are common in automotive, rotorcraft and industrial transmissions [7], we study the behavior of a 1-D array of impact coupled oscillators.

MODEL

In our study, the local element placed on each lattice site is a linear dissipative driven oscillator which by itself shows no chaos. The complex behavior in a 1-D array of these simple systems is caused by the impact coupling. The coupling between the oscillators occurs only in the event of a collision between the neighboring elements. The particular model studied consists of n oscillators, of equal mass, which are driven globally by a sinusoidal drive. Using the reduced coordinates, introduced earlier, the equations of motion of the i^th oscillator, between collisions with the neighboring oscillators, is given by

QUALITATIVE ANALYSIS

Numerical investigations of the above model have been done under two different boundary conditions. In the first case, the oscillators at the left and right end are free to oscillate whereas, in the second case, they are assumed to remain fixed at their equilibrium positions serving as fixed end walls.

In the first case of free boundary conditions, we found that the system eventually settled into simple periodic behavior. This is because of the tendency of the system to settle down to a pattern of oscillation, where the oscillators avoid collisions. In this case, the collisions, with dissipation, eventually reduce the coupling between the oscillators to zero, in which case they behave essentially like n independent oscillators.

In order to force the oscillators to collide more often, we changed the boundary conditions from free boundary conditions to fixed boundary conditions, in that the oscillators at the ends are held fixed at their equilibrium position thus, acting as rigid walls. As a result, the oscillators colliding with these walls now reverse their velocity. This velocity is reduced in magnitude depending upon the value of the coefficient of restitution. In contrast to the case of the oscillation patterns under free boundary conditions, simulations done under fixed boundary conditions show extremely complicated patterns of oscillations.

SPACE-AMPLITUDE PLOT

In order to investigate the dynamical response of linear chain of oscillators coupled through collisions between nearest neighbors, we graphically represent one of the phase variables, say, the deviation of the position of each oscillator from its equilibrium point observed at some fixed drive phase plotted versus the equilibrium position of the oscillators. This is repeated for several drive cycles and the points are overlayed on top of each other. This is like taking the snapshot of the system once every drive cycle and overlaying these snapshots for many successive drive cycles. We refer to this graphical representation as space-amplitude plots.

The plots above, show the space-amplitude plot representing different dynamical responses of a 16 oscillator system for different parameter settings. The position of each oscillator, when the drive phase goes through zero, is plotted on the Y-axis, and the equilibrium site of each oscillator is plotted on the X-axis. The points for each drive cycle are connected by a line. For every plot shown in these figures there are 16 oscillators and a wall at the left end and at the right end. The parameters are set at K = 0.7 and b = 0.05. The plots show the pattern of oscillations for four different values of A and w: (a)(A = 3.5, w = 1.5); (b)(A = 1.8, w = 1.5); (c)(A = 1.7, w = 1.5); and (d)(A = 1.8, w = 1.47). In Fig.1(a) all the oscillators are oscillating in a period-one state, Fig.1(b) shows when most of the oscillators are oscillating in a period-two state and the rest are in period-one state, Fig.1(c) shows where all the oscillators are in a chaotic mode, and Fig.1(d) shows a pattern of oscillation where the first half of the eight oscillators are in chaotic mode and the other half are in a periodic mode (period-one).

The plot below shows how the dynamic response of imapct-coupled oscillators change with slow change in the drive amplitude, A. They show the space-amplitude plot for set of 15 oscillators as a function of A for fixed values of w = 1.5, K = 0.7 and b = 0.05. The value of A from (a)-(f) changes from A = 2.16 to A = 2.21 in steps of 0.01. fig.(a) shows the oscillators in period 2 for A = 2.16 and at A = 2.17, shown in Fig.2(b) the oscillators show intermittent oscillations. As the amplitude is increased to 2.19 the system settles down to high period orbit with different oscillators in different periods. Most of the oscillators show period 4 behavior and at A = 2.18 the system shows intermittency. Until at A = 2.21 all the oscillators seem to oscillate chaotically.

The phenomenon of multistability is very commonly seen for linear chain of impact oscillators. Depending upon the initial conditions, the system settles down to one of the many accessible attractors. The figure above shows the space-amplitude plot for A = 2.2, w = 1.5, K = 0.7, and b = 0.05 for the two possible patterns of oscillations the system settles down to depending on the initial conditions.

SPACE-TIME-AMPLITUDE PLOT

To show how the complex oscillation patterns evolve in time, we represent the data obtained by numerical simulation, for coupled impact oscillators, using space-amplitude-time plots. Time evolution is displayed by displacing each successive space-amplitude plot (recorded when the phase of the drive passes through zero) along a time axis plotted in the third dimension. The positions of the oscillators at zero phase are plotted along Z-axis, the equilibrium sites of each oscillator are plotted along the X-axis and the time is along Y-axis. All plots are shown after initial transient behavior were discarded. Fig.4 and 5 show such a plot for the A = 1.7, K = 0.7, w = 1.5, and b = 0.05 where the system of 16 oscillators shows intermittent bursts of nearly periodic and aperiodic behavior. These plots show temporal evolution of oscillation pattern of the impact-coupled oscillators. Assuming that the start time for Fig.4(a) is zero, the figure shows the time evolution for the first 60 cycles. Figure4(b) shows the temporal evolution for the next 60 cycles. Similarly, Fig.4(c) and (d) show the temporal evolution for the next 120 cycles. Thus, each plot shows the temporal evolution for 60 cycles. Fig.4(a) shows the nearly periodic oscillation for all the oscillators in the chain. In Fig.4(b) the left half of the chain shows nearly periodic response and the right half shows aperiodic response. The system, as a whole, soon settles down again to nearly periodic response as seen in Fig.4(b) and Fig.4(c). In Fig.4(c) around the 180^th cycle the left half of the chain is beginning to respond

aperiodically and in fig.(d) nearly the whole chain begins to respond aperiodically. Similar patterns of burst of periodic and aperiodic response can be seen in the next four plots shown in Fig.5 for a further 240 cycles. It should be mentioned that the particular case shown in figs.(4) and fig.(5) has the chain split into two independent pieces since oscillator 8 and 9 never collide. Such intermittent bursts of nearly periodic and aperiodic response of the system were also seen for several other parameter settings.

SHOCK WAVE ANALYSIS

To understand the cause of such intermittent response and the transition of the system from one mode of oscillation to another, we analyze the response of the system in terms of the cascading effect of the collision of the first oscillator with the left wall down the chain and a similar cascading effect when the last oscillator collides with the right wall. We make this notion concrete in what we call a shock wave plot. The plot is made in the following way. Without loss of generality, we assume the following convention: if there is a collision between k and k+1 oscillator, we say k^th oscillator has collided. We also know the time at which the oscillator collided. To identify the shock wave generated from the left wall, we assume that the first oscillator has collided with the left wall, say, at time t = t₀. We look for the collision between the first and the second oscillator, and then second and the third oscillator and so on. We record the site index and the time at which the collision took place. In the event of a collision between the left wall and the first oscillator we say a new shock wave is generated at the left wall. Thus, at any given time there are, in general, more than one shock wave propagating down the chain from the walls. The collision that make the shock wave is the one that has its collision index one more than the last entry made in the construction of the shock wave. For the case of shock wave generated at the right wall we follow the same

Shock wave plot for set of 16 oscillators for A = 1.7, w = 1.5, K = 0.7 and b = 0.05. The plot explains the cause of intermittent bursts of periodic and aperiodic response shown in Fig. 4.	Shock wave plot for set of 16 oscillators for A = 1.7, w = 1.5, K = 0.7 and b = 0.05. The plot explains the cause of intermittent bursts of periodic and aperiodic response shown in Fig. 5.

procedure as above except that the collision that make the shock wave is the one that has its collision index one less than the last entry made in the construction of the shock wave. The time of collision of the oscillator versus the equilibrium position of the colliding oscillator shows the time progression of a `wave' of impacts and we refer to the graph as a shock wave plot.

Fig.6 and 7 show the shock-wave plot originating from both left and right wall for the oscillation patterns shown in Figs.4 and 5. The shock-waves originating from the left wall are shown in the lower half of the plots and the ones originating at the right wall are shown in the upper half of the plots. It should be noted that the corresponding space time amplitude plot in Fig.4 and shock wave plot in Fig.6 have the same time range. This is also the case for Fig.5 and Fig.7. This makes it easy to see the effect of new shock waves, generated at either wall, on the temporal evolution of oscillation patterns. The shock-wave plot shown in Fig.6 explains the reason for intermittent bursts of aperiodic response in the space-time-amplitude plot shown in Fig.4. Similarly, Fig.7 explains the reason for transition from periodic response to intermittent bursts shown in Fig.5.

Comparing Fig.4a and Fig.6a, it is clear that the response of the system would be regular for a regular shock wave originating at the left and the right wall. The aperiodic response in Fig.4b can be attributed to the creation of a new form of shock wave generated at the right wall (Fig.6b). This additional shock wave could be thought of as a perturbation to the system and dies out after a few drive cycles after which the system settles down to its periodic response as seen in the latter half of Fig.4b and the first 58 cycles in Fig.4c. The generation of new shock waves originating first at the left wall and then at the right wall and its more frequent creation (Fig.6d) accounts for the aperiodic response of the whole system for a relatively longer time, as shown in Fig.4d. It should be noted that the intermittent burst of aperiodic response and periodic response of the system was observed for 100,000 drive cycles and thus appears not to be a transient behavior of the system.

The figure below shows the Poincaré section, taken at f = 0, for the first oscillator from the left wall. It clearly shows the non-stationary nature of oscillations. The arrow indicates the movement of the periodic orbit along the unstable direction of the fixed point at the Poincaré section. The plot is made for 50 cycles after the transients have died out.

BIFURCATION PLOT

In this section, we investigate the dynamic response of the impact-coupled oscillator systems as a function of the system parameters. With the help of both qualitative and quantitative analysis, we establish the effect of changing dissipation in the system response. In particular, we show that by increasing the dissipation in the system (by decreasing K or increasing b) for fixed value of A (or w) the system tends to respond periodically for a wider range of w (or A) respectively. This is shown by plotting the bifurcation diagram and by calculating the Lyapunov spectra of the system as one of the system parameters is changed.

In Fig.10, the bifurcation diagram for set of 8 coupled oscillators are plotted. The parameters are set at K = 0.5, w = 1.9, and b = 0.05. Figure 10(a-h) shows the value of the position of each oscillator (leftmost oscillator to rightmost oscillator) measured when the drive phase goes through zero as a function of A. The value of A is decreased from 6 to 1, in steps of 0.001, after every 300 drive cycles. The first 250 drive cycles are discarded and the last 50 cycles are recorded for each value of A.

It is clear from the plot that for K = 0.5, for most of the parameter space in A (between zero and six) the system responds periodically. Between A = 1.0 to A = 2.1 all the oscillators are in the non-impacting mode. At around A = 2.1 there is a collision between the oscillators on each end of the chain and the walls leading to collisions between the neighboring oscillators. From Fig.10d and 10e - which is the bifurcation plot for the 4^th and the 5^th oscillator from the left wall, it is clear the shock waves originating at the left and the right wall do not make it to the middle two oscillators. Thus, for A between 2.1 and 2.6, three oscillators, from either walls, exhibit chaotic response and the middle two oscillators show non-impacting periodic response. It should be noted that the chaotic response of, say, first oscillator and the last oscillator are not the same, despite the symmetry of the system. The same is true of the 2^nd and the 7^th pair and the 3^rd and the 6^th oscillators. This implies multistability.

Figure10 shows sudden destruction of attractors, jump phenomena and various kinds of period 1 orbits with different sequence of collisions with the neighboring oscillators or the walls. As seen from Fig.10 and number of other parameter settings for which the simulation was done, the oscillators in the chain, for a wide range of A seem to settle down to a periodic response. This is surprising given the fact that all oscillators in the chain are colliding with their neighbors.

Next, we discuss the effect of dissipation and coupling on the dynamic response of the impact-coupled systems. We first discuss the effect of viscous damping on the system. Numerical experiments done on a set of 8 coupled impact oscillators indicate that as we decrease the viscous damping in the system by decreasing b, the volume in the parameter space for which the system responds chaotically increases. The sum of the positive Lyapunov exponents decreases with increasing b. In Fig.11 we show the bifurcation sketch for the first oscillator from the left wall as a function of w. The parameters are set to A = 2.0, K = 0.7, and for figures on the top row are for b = 0.05, 0.1, 0.15. The Lyapunov spectra for the corresponding parameter values are plotted below. The details of the derivation and method, used to calculate the Lyapunov spectra, is given in section 5. In Fig.11 only the largest 10 Lyapunov exponents are plotted.

It is clear from the plot that the parameter range for which the system responds chaotically is more for the case of low damping as compared to the case of higher damping. Similar qualitative result was also found by increasing the value of K. That is, for higher values of K (the measure of dissipation and coupling) the sum of the positive Lyapunov exponents was higher than for lower values of K.

Arranging the Lyapunov exponents of an attractor of a continuous-time dynamical system in the descending order, like l₁ ³ ¼ ³ l_n and if j be the largest integer such that l₁ + ¼+ l_j ³ 0. The Lyapunov dimension (D_l) as defined by Kaplan and Yorke [10] is

along the y axis is same and it is clear from the plot that with increase in damping the Lyapunov dimension decreases. Similarly decreasing K has an overall effect of lowering D_l or the range of parameter region for which the system exhibits chaotic or hyperchaotic behavior.

QUANTITATIVE ANALYSIS

In a dynamical system with N-dimensional phase space in Â^N, there exists N independent directions. Corresponding to these directions, there exist N Lyapunov exponents that characterize the rate of divergence (or contraction) of small deviations along these directions. In this section, based on Nordmark's work [11], we derive the Jacobian of the local Poincaré map at the collision surface.

DERIVATION

Let the flow in Â^N be given by [(x)\dot] º [(dx)/ dt] = F(x) and a smooth surface in Â^N be defined by scalar function H(x). Also, let x_s be a point on the surface S and

In the rest of this section, we derive the Jacobian of the local Poincaré map at the colliding surface when (a) two oscillators collide and (b) when an oscillator collides with a stationary wall. For the case of collision between two oscillators, we assume that the equilibrium position of the oscillator are separated by a distance a, the left wall is at zero position and the right wall is at a distance of (n+1)a from the left wall, where n is the number of oscillators in the chain. As a notational convention, the primed coordinates are measured in an absolute reference frame with the left wall as the origin and unprimed coordinates are measured with respect to the local coordinate system with the origin at the equilibrium position of each oscillator. Thus, x_k^¢ = x_k + ka and x_k+1^¢ = x_(k+1)+(k+1)a are the position of the k^th and (k+1)^th oscillator in the chain. To carry out the analysis, it is advantageous to express these variables in the center of mass and relative coordinates. Introducing the center of mass coordinate as R_k = (x_k^¢+x_k+1^¢)/2 and the relative coordinate as r_k = x_k^¢-x_k+1^¢ we can express these variables in terms of the local coordinate system. To keep the discussion and notation simple, we assume there are two oscillators in the system and we are considering the collision between them. Thus, the positions and the velocities in the center of mass and relative coordinates, for our impact system, can be expressed as,

To calculate the flow projection mapping from constant phase plane to the surface of collision(while approaching the impact surface), we first use Eq. 8 to calculate the Jacobian and then we take the projection. The surface of collision, when the two oscillators collide, is defined by H(r) = r = 0 and thus the gradient ÑH = (1,0,0,0,0) and ÑH ·F = V_r. Thus, using Eq. 8, we get the Jacobian of the map. To find the flow projection map, we pick the columns corresponding to variables (r, V_r, R, V_R) and the rows corresponding to variables (f, V_r, R, V_R) (because locally, about the colliding surface, the flow projection map is D_x P₂ : S_{[`(f)]</sub> ® S<sub>[`r]},), this gives

The map relating phase variables just before to phase variables just after the impact can be readily found by realizing that at the impact surface, the variables f and R go through no change and the velocity relations are given by V_r⁺ = -KV_r^- and V_R⁺ = V_R^-, where K is the coefficient of restitution. We pick columns corresponding to (f, V_r, R, V_R) and rows corresponding to (f, V_r, R, V_R) (because locally, about the colliding surface, the flow projection map is D_x P₃ : S_{[`r]</sub> ® S<sub>[`r]}). Thus,

To calculate the flow projection mapping from surface of collision (leaving the impact surface), to constant phase plane, we use Eq. 8 to calculate the Jacobian and then take the projection. The surface of collision, when the two oscillators collide, is defined by H(f) = f- f_s = 0 and thus the gradient ÑH = (0,0,1,0,0) and ÑH ·F = 1. Thus, using Eq. 8, we get the Jacobian of the map. To find the flow projection map, we pick the rows corresponding to variables (f, V_r, R, V_R) and the columns corresponding to variables (r, V_r, R, V_R) (because locally, about the colliding surface, the flow projection map is D_x P₄ : S_{[`r]</sub> ® S<sub>[`(f)]},), this gives

The derivation for the collision matrix for the case of oscillator colliding with the left wall can be derived using the same arguments as above and the expression, in local coordinates is

The determinant of the matrices, given by Eqs. 14 through  16, is K². This implies that the volume in phase space reduces by this factor for both collision between oscillators and between oscillator and the wall.

LYAPUNOV SPECTRA CALCULATION

In this section we present results of successful computation of Lyapunov spectra using the results derived in the previous section. We show the Lyapunov spectra calculated for a set of 16 oscillators at the Poincaré section taken when the global drive phase goes through zero. The phase space dimension of the system, at the Poincaré section, is 32 - position and velocity of each oscillator. Thus, there are 32 Lyapunov exponents to be calculated.

The computation of the Lyapunov exponents is made as follows. After the transients have died out, we follow the system from when the drive phase was zero the previous time. All the oscillators are integrated forward in time until a collision between oscillators are detected. The time between the last piercing of the Poincaré surface of section and the collision is is used to evolve the linearized set of equations forward. The linearized variables are evolved forward, between impacts, by a (32×32) matrix. The matrix has nonzero entries along the diagonal in the form of (2×2) matrix given by Eq. The matrix is given by

g = Ö{1-b²} and f = g(dt). dt is either the time elapsed between the last crossing of the Poincaré section and the first collision, or between two successive collisions. Assuming there was a collision between the i^th and the (i+1)^th oscillator, we use the following matrix to evolve the linearized variables from the state just before the collision to state just after the collision

All other nondiagonal elements of the J matrix are zero and the diagonal elements are 1. This ensures that we change the phase variables of only those oscillators that are involved in the collision. For the case of collision between the left wall and the oscillator we have J₁₁ = J₂₂ = -K, J₁₂ = 0 and J₂₁ = (K*A^-+A⁺)/V^-. All other nondiagonal terms are zero and the diagonal elements are 1. For the case of collision between the right wall and the oscillator we have J_2N-1,2N-1 = J_2N,2N = -K, J_2N-1,2N = 0 and J_2N,2N-1 = (K*A^-+A⁺)/V^-. All other nondiagonal terms are zero and the diagonal elements are 1. Having used the appropriate J matrix to evolve the linearized phase variables from state just before the collision to state just after the collision, we use the P matrix until the next collision. Thus, we will have series of P and J matrices in a given drive cycle. We compose the D matrix from these P and J matrices as follows

Lyapunov spectra calculation for a linear chain of 16 impact-coupled oscillators as a function of the drive amplitude A. The other parameters are set to

w = 1.9, b = 0.05, and K = 0.7. The scan parameter,A, is changed in steps of 0.005. Every time A is changed, we wait for 100 drive cycles for transients to die out and use the next 300 cycles to calculate the average Lyapunov exponents in the units of bits/cycle.

REFERENCES

[1]

K. Kaneko, Phys. Lett. 111A, 321 (1985).

[2]

K. Kaneko, Physica 23D, 436 (1986).

[3]

K. Kaneko, Physica 34D, 1 (1989).

[4]

K. Kaneko, Prog. Theor. Phys. 72, 480 (1984). J. Keeler and J. Framer, Physica 23D, 413 (1986). H. Chate and P. Manneville Physica 32D, 409 (1988).

[5]

K. Kaneko, Phys. Lett. A 149, 105 (1990).

[6]

D. K. Umberger, C. Grebogi, E. Ott, and B. Afreyan, Phys. Rev. A, 39 4835, (1989).

[7]

C. Padmanabhan and R. Singh, J. Sound Vib. 155, 209-230 (1992).

[8]

C. J. Begley and L. N. Virgin, Journal of Sound Vibration,211, 1998, 801;

[9]

W. H. Press, B. P. Flannery, S. A. Teukoslky, W. T. Vetterling, ``Numerical Recipes in C'' Cambridge University Press, pp. 60-71, 1990.

[10]

J. L. Kaplan and J. A. Yorke, Chaotic behavior of multidimensional difference equations, Lecture Notes in Mathematics, Springer-Verlag, New York, N.Y., pp. 228-237 (1979).

[11]

A. B. Nordmark, J. Sound Vib., 145, 279 (1991); A. B. Nordmark, J. Bifurcation Chaos, 2 597 (1992).

[12]

A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano, Physica D 16, 285 (1985).