Nonsmooth Chaos

Photocopying permitted by license only the Gordon and Breach Science Publishers imprint.

Printed in Malaysia.

Using Chaos Synchronization to Estimate the Largest Lyapunov Exponent of Nonsmooth Systems

ANDRZEJ STEFANSKIandTOMASZKAPITANIAK*

DivisionofDynamics, Technical Universityof^Lodz,Stefanowskiego 1/15,90-924Lodz, Poland (Received15 April1999)

We describe the method of estimation ofthe largest Lyapunov exponent ofnonsmooth dynamicalsystemsusingthepropertiesof chaossynchronization.The method is based on the coupling of twoidentical dynamical systems and is tested on two examples of Duffing oscillator:(i)withaddeddry friction,(ii)withimpacts.

Keywords." Chaos synchronization,Lyapunovexponents, Nonsmooth systems

1. INTRODUCTION

The estimation of

Lyapunov

exponents is one of the fundamental tasks inthe studies of dynamical systems.

Lyapunov

exponents measure the expo- nential ratesof divergenceorconvergenceofnearby trajectories in state space. Theoretical studies of Oseledec [1] and numerical algorithm of Benettin etal. [2]andWolfetal. [3]allowaneasyestimation ofthe spectrum of

Lyapunov

exponentsfor smooth systems described by the known motion equation.

If such equations forthesystemare notknown^orit is nonsmooth, the estimation of

Lyapunov

expo- nents is not straightforward. In recent years methods ofthe calculationof

Lyapunov

exponents for dynamical systems with discontinuities have beenproposedin somepapers[3,4].

Manyrealengineering systemscan beconsidered as discontinuous, i.e. themechanical systems with dryfriction orwithimpacts.Thedryfriction inany engineering systems often leads tothe appearance ofthephenomenon of stick-sliposcillations which also haveadiscontinuous nature.

Inthispaperweproposeamethod ofestimation of the largest

Lyapunov

exponent using chaos synchronization[5,6]. Ourmethodis motivatedby Fujisaka and Yamada’stheoretical studies

[5].

They have found thelineardependencebetween the syn- chronizationvalue of the couplingcoefficientoftwo identical continuous dynamical systems and the largestvalueof

Lyapunov

exponent of suchcoupled systems. Asanexample they have used the Lorenz system. Wehave shown thatthis approach canbe applied for nonsmooth dynamical systemsaswell.

Correspondingauthor.

207

2. THEORETICAL BACKGROUND FOR ESTIMATION PROCEDURE

describes the norm (q R

>_ 0)

ofthe state

differ-

encevector Consideradynamicalsystemwhich consistsoftwo

identical n-dimensionalsubsystems coupled byone- to-onenegative feedbackmechanism with apairof coupling coefficients

kx

^and

ky.

^{The first order}

differentialequations describing suchasystemcan bewritten as

k

-f(x) +

^kx(y-

^x),

: ^{f(y) + ky(x y),}

"k ^w where x,

yER"

and

kx,y [kx,y; kx,y;..., x,y]

R’_>

0 is a coupling vector. Let us assume that for

kx-ky-O

each system

dx/dt=f(x)

and dy/

dt=f(y) evolves on an asymptotically stable chaotic attractorA.

Inthefirstpartofourconsiderationsweassume that thevalues of thecouplingcoefficientsareequal to zero

(kx-ky-O).

^{In this case we obtain two} separate, identical dynamical systems given by equations

k

--f(x),

)?

=f(y). (2)

The solutions of

Eq. (2)

starting from different initial conditionsrepresenttwoindependent trajec- tories onthe attractorA.Ifinitial conditions are the samethe evolutionofboth subsystems is identical andwehaveideal synchronization

(x

y).

Nowweintroduce a new variablezrepresenting thestate

difference

between both subsystems during the time evolution. This variable is defined by the expression

z=x-y,

(3)

where z R

.

^In further analysis to simplify the notation and allow for the visualization we took n 3 and assumed that the evolution on the attr- actorA is characterizedby onepositive

Lyapunov

exponent.The absolutevalueofstate

difference

q-[I

2’ x-y

(4)

e-2. (5)

The image of the vector equation (Eq.

(5))

ⁱⁿ phase spaceisshown inFig. 1.Thestate

difference

vector is fixed to theendof vector

07

^{anditsend is in}

touchwiththeendof vector

(see

Fig.

1).

The time evolutionofthestate

difference

^{vector5’isdefined}

inz; systemof coordinatesandcanbe determined as a sum of its 5’; contributions in all directions of phase space

(see

Fig.

2):

*-- +i. ⁽⁶⁾

i=1

Theorigin^{of this}system is fixed to theorigin ^of vector2’.The axes of

zi-directions

^{areparalleltothe} appropriateaxes ofphase space.

The time variations of the norm of state

difference

^{vector are determined by first time}

derivative ofEq.

(4):

q-k- forx>_y,

O--k

^{for x}

^<

^y,

(7)

andafter substitutingEq.

(2)

ⁱⁿ

(7)

^weobtain gl

f ^(x) f ^Y)

forx_>y,

O f ^Y) f (x)

^forx<y.

xi,Yi

FIGURE The image ofthe vector equation (Eq. (5)) in phasespace.

FIGURE2 Thespatial orientationof thesystem of coordi- nates associated with: (1) A-distance vector (dotted line), (2)statedifference^vector(continuous line).

For small state

difference

^{vector, i.e. q}

<< ]A],

where

IAI

^R

^>

0 is the attractor size(maximum distance between two points on the

attractor)

ⁱⁿ phase space, ^we can assume that the distance between trajectories of the subsystems under con- sideration is given by^linearizedequation resulting from definitionofLyapunovexponent

8i- lSi0le x’t, (9)

where

A;

is a

Lyapunov

exponent,

6

^{is the norm}

of the contribution

6i

^of A-distance vector 6 in direction associated with /-number of Lyapunov exponent and 6;0is an initial distance in the same direction. The total A-distance vector is a sum of contributions

(see

Fig.

2)"

i=1

Thenormofthisvector isgivenby therelation

e lll (6i0ea")

²

(11)

i=1

The

Lyapunov

exponents ^are related to the expandingorcontractingnatureofdifferent direc- tionsinphase space.The spatialorientation of the

coordinates system ofdirections associated with a given exponentvaries inacomplicatedwayduring the evolution on the attractor. Forthat reasonthe axesof z,-directions

(state difference ^vector)

^{are not}

coveredwith the axes associated with A-directions (k-distance

vector).

Both systems ofaxes have the same origin and the angular difference between themis defined in eachmomentoftimeby angles

O,(t)

^whichdependon theposition of trajectory^on theattractorduringthe time evolution

(see

Fig.

2).

The direction corresponding to exponent A;=0 is always tangentialtophase trajectory.

From the above considerations the equality of state

difference

^vector

’

^{andA-distancevector8}

^(for

nearbyorbits)^arises:

5’- 6.

(12)

To make further considerations easier we have introducedthe followingnotations:

(1)

/max maximum

Lyapunov

exponent,

(2) 6o

initial distance in kma-direction,

(3) eXAm.-

^{difference between /max and other}

Lyapunov

exponent

(AAj >

0),

(4)

,)kj--/max-

/’j

the rest of

Lyapunov

exponents,

(5) @=m/g0

initial distance in

A.-direction, m..

is aconstantvalue.

Putting the above substitutions in Eq.

(11)

we obtain

2 -2AA.

6

(5e

^2Amaxt

+ m/e

^.t

i=l

(13)

During the time evolution the sum in

Eq. (13)

decreases to zero and the distance associated with

Amax-direction

^{becomesdominant.}

^Hence,

^{the norm}

ofthe A-distance vector can be finallywritten in a reducedform:

5

1601exmaxt. ⁽¹⁴⁾

Fromthe equality given by

Eq. (12)

the similar equality of the norms of both vectors results.

Itallowsusto makeasubstitution:

q-

I01exmaxt. ⁽¹⁵⁾

The first time derivativeof

Eq. (15)

is givenby the relation

0 ^/max]01 eAmaXt" (1 6)

The comparison of Eq.

(8)

and

(16)

gives the equalities:

x>_y"

x<y"

f ^(x) f ^y) /maxl0le

^Amaxt

f(y)-f(x)- ,maxl0le

^A

t. (17)

In the next step we introduce the nonzero coupling parameters

kx

^and

ky

ⁱⁿ

Eq. (1).

In this casethe equation describing thefirst time derivative of the norm state

difference

vector assumes the following forms:

x>_y:

x<y:

il f ^(x) f ^{y) (kx} + ^ky)q

O -f(Y) -f(x) (kx + ^{ky)q. (8)}

Using

Eq. (17)

ⁱⁿ

(18)

weobtain:

dq

/maxl50}e

^Amxt

(kx

-I-

ky)q. (19)

dt

Toseparatevariables we haveto divideEq.

(19)

by expression q/dt andnextusethesubstitution q exp(-/max

1)

arisingfrom

Eq. (15).

^Now

Eq. (19)

isgivenintheform:

dq

[,max (kx -+- ky)]

^dr,

(20)

q

andhas thefollowingsolution"

q qo

e[Amx-(kx+k)] (21

where q0 is a constant defined by the initial con- ditions.The synchronization betweencoupled sub- systems ofEq.

(1)

takesplace when thefixedpoint ofstate

difference

^{becomes a} stable attractor, i.e.

z

Zs

^{=0. From Eq.}

(21)

one finds that fulfilling theinequality

/max %

kx

^q-

ky, (22)

guarantees stability and then the norm of state

difference

^{vectordecreasesto zero.}

According to Eq.

(22)

the synchronization of two identicalsystems ^is possible when the sum of coupling parametersis largerthanthe value of the largest

Lyapunov

exponent of these systems. If

"max <

⁰then both systems synchronizeinastable fixed point. For /max--0 such systems have a periodic solution and the synchronization takes place just after a small coupling is introduced or

evep_.withoutcoupling, whenthe excitations inboth systems have the same phase. If

"max >

^{0 we can} observe a linear covering of two independent effects:

(1)

exponential divergence of nearby trajectories associated withpositive

Lyapunov

exponent,

(2)

exponential convergence due to the introduced

coupling.

The first of these phenomena occurs only in direct neighbourhood of the synchronized state wherelineareffects aredominant. The second one acts inthe entirephase space. For that reasonthe describedcovering of both effectstakesplace only nearbythe synchronizedstate.

The fulfilling ofthe inequality (Eq.

(22))

causes that synchronization manifold x=y becomes an attractor,i.e.theevolutionof6-dimensionaldynam- ical system (Eq.

(1))

^is reduced to the attractor A located on x=y manifold. In opposite case (,)kma

> kx--ky)

the manifold x=y is a repeller and synchronizationisimpossible.Inthe remaining parts of this paper we give evidence that this property can be very useful in the estimation of

Lyapunov

exponents ofthe systems for which the straightforwardapproach[2,3]^isimpossible.

3. PRACTICAL ESTIMATION PROCEDURE

Based onthe abovementionedproperties of chaos synchronization we have proposed a method of estimationofthe largest

Lyapunov

exponent.This method is independent of the type of dynamical

system because it is not necessary to assume con- tinuous character ofexpressions

f(x)

and f(y) to obtain synchronization of

dx/dt=f(x)

and dy/

dt--f(y) subsystems. In this paper we have con- centrated on the application of our method for finding thelargest

Lyapunov

exponentofdynami- cal systemswith discontinuities.

To simplify the procedure we have assumed unidirectionalcouplinginEq.

(1).

In this caseone of the coupling coefficients is equal to zero (say, kx=O,

ky=k)

^and

Eq. (1)

^is reduced to the followingform:

k

-f(x),

j

f(y) + k(x y). (23)

Nowthe condition ofsynchronization(Eq.

(22))

isgivenby the relation

,max <

^k.

(24)

The smallest value ofthe coupling coefficientk, for which the synchronization takesplace

ks

^isequal

tothe maximum

Lyapunov

exponent

Amax.

Toapplyourmethod for any dynamical systemit isnecessaryto buildadouble systemwithcoupling accordingto

Eq. (23).

^{The next} stepis a numerical research of the synchronization parameter

ks

for such augmented system. If thetested coupling parameterk reachestheboundary value

ks

^{then the}

largest

Lyapunov

exponent of investigated system amounts to

ks.

The simplestway to findthevalueof

ks

^ismaking

a bifurcation diagram of the state

difference

^z

between coupled subsystems of

Eq. (23)

in function of tested coupling parameter k.

An

example of suchadiagramispresentedinFig.3.Wecanobtain thesearchedvalueof

ks

^{asapointonthehorizontal}

axis (coefficient

k)

wherestate

difference

^zreaches

zerovalue.

The main problem which has been observed during practical applications ofourmethodisalong time of transient motion before the system under consideration achievesthe synchronizedstate. This problem appears near the boundary value

ks.

0.000 0.125 k ^0.250

FIGURE 3 The bifurcation diagram showing state dij]br- ence _z2 versus coupling coefficientk forcoupled Duffing os- cillators (Eq. (27)) and estimated value of the largest Lyapunovexponent;d 0.100.

To make faster achievement of searched

ks

^value

possible we have divided the phase space of the researchedsysteminto tworegions. The first one is a direct neighbourhood of the synchronization subspace x-y. This region is bounded by a coefficient which is a radius of n-dimensional sphere with a centre in the fixed point

Zs-0

(Fig.

4(a))

or it can be considered as a radius of n-dimensional "corridor"surrounding thesynchro- nization manifold (Fig.4(b)).^Thesecond regionis the remaining part of phase space. The way to achieve faster synchronization is to eliminate the long periodsof time whenphasetrajectories evolute farawayfromsynchronizedstate. Forthispurpose wehave applied the method calledelastic coupling whichforcescoupled systemstothe evolution inthe neighbourhoodofx-ymanifold.The main idea of this method is the sudden jump of a coupling coefficientwhenstate

difference

^{crossesthebound-}

ary e value. Inside the spherez

<_

e (Fig.

4(a))

the coupling parameter in

Eq. (23)

has testedkvalue.

Beyond thesynchronization region

(z >

e) Eq.

(23)

assumesthefollowing form:

f(y) + K(x (25)

where Kis the strong coupling parameter

(K >>

k) which does not allow the system evolving beyond e-radiussphere.

(a)

zi

(b)

Yi

Zi+l

Strong

coupling Tested coupling

xi

FIGURE 4 The idea of elastic coupling (see description in Section 3).

Thepropervalues ofcandKcanbeelaboratedin awayof individual studies ofthe given dynamical system.

4. EXAMPLES

In this section we present two types ofdynamical systemswith discontinuities.Asanexampleofsuch systemwehave used theclassicalDuffingoscillator:

(i)^withaddeddry friction, (ii)^withimpacts.

4.1. Duffing Oscillator with

Dry

Friction

The equation describing the first example is as follows"

1

^-2,

.

^-p

^{o(t) + , ( )}

d

+

^T]I

X21 ⁺

^#0-[-

T]2X22

sign(x2).

(26)

Thesystem(Eq.

(26))

has beencreatedby adding to the well known Duffing equation a term describing dry friction according to Popp-Stelter formula [7]. Its discontinuous nature makes any direct calculationof

Lyapunov

exponentsdifficult.

After putting the system under consideration (Eq.

(26))

in

(23)

^{we obtain}the augmented system inthe following form:

-1

^X2,

22

^-p

cos(cot) +

^cx,

(1 x) hx2

d

+

r/1

x2] +

#0

+ zx

^sign(x2),

1

2^@

(X1 --1), (27)

2:

^p

cos() + , ( )

d

1

₁

y2l +

0

+ + (x2

In ^our numerical simulations we have assumed the following values of parameters: ^c=1.00, h 0.25, p 0.30, 1.00, ₁ 1.42, ₂ 0.005, 0 0.25, 0.05.

Thenthelargestvalue of

Lyapunov

exponentof system(Eq.

(26))

has been determinedaccordingto the way described in Section 3. The results of calculations aregivenbythegraph showninFig. 5 which presents the changes ofthis exponent as a function of d parameter representing amplifying dryfrictionforce.

The bifurcation diagram of system (Eq. (26)), showninFig. 6presents thedependenceofvelocity x2ondparameter value. The dparameter rangeis identical to thatinFig. 5. Wecan seebands when chaotic motion occurs

(see

Fig.

6)

whenLyapunov

0.16

--ks

0.12

0.08

0.04

0.00

0.00

0’.04

^{0.08 0.12} _d ^0.16

FIGURE 5 The largest Lyapunov exponent XmLx ^{versus d}

parameter value for the investigated Duffing oscillator with dryfriction(Eq.(26)).

impacts. The equations describing such a system canbe written informs:

-p

co ( t) +

Ixl

²

^Xo:

^{X2a =--Rx2b,}

⁽²⁸⁾

for two oppositebuffers,and

x <Xo: 2

=x2

X1

X0:

X2a----Rx2b,

(29)

FIGURE6 The bifurcation diagram of velocityx^{versus d}

parameter forDuffingoscillator withdryfriction(Eq. (26)).

exponent takespositivevalue(Fig.

5).

Also periodic motion occurs when zero value of Lyapunov exponent is foundin Fig. 5. These results support our hypothesis that the values of the largest

Lyapunov

exponents obtained as above for sys- tems with dry friction can agree with real ones (impossible to determine analytically) with high approximation.

4.2. Duffing Oscillator withImpacts

As

the second example of discontinuous system we present again Duffing oscillator but now with

forsingle buffer.

In the above equations

X0

^{is a position of the} buffer, Ris the coefficient ofrestitution, x2v is the velocity in a moment before impact and _x2 is the velocity just after impact.

Impact

motion has a bit more complicated discontinuous naturethanthe motionofthesystem in the previous example. In phase space of the systemwith impacts the simultaneous existence of different types of attractors is possible. For that reason to determine the largest

Lyapunov

expo- nents of the systems (Eqs.

(28)

and

(29))

^we have applied theabovementionedmethod calledelastic coupling.

We have estimated the largest Lyapunov expo- nents of systems under consideration for chosen values ofthe parameters. The results ofnumerical calculations are presented in Figs. 7 and 8 which show thephase portraitsand associated with them the largest

Lyapunov

exponents

(on

the top of picture) which have been obtained using the describedsynchronization method.

Wecan seethatchaotic motion occursby positive

Lyapunov

exponents. It has also been confirmed (Figs.

7(a)

and

8(a))

that synchronization of periodic systems can appear without coupling

(ks "max- ⁰⁾

between testedsystemswhenexcita- tionhas the samephase.

(a)

1.50

max

^{k =0.000}

X2

0.00

-1.50

-0.80 0.00 0.80

Xl

,max--ks^=0.178 (b)_1.50

0.00

-1.50

-0.50 0.00 0.50

Xl

FIGURE 7 Thephase portraits presenting the solutions of Duffing oscillator with impacts (twobuffers Eq. (28)) and the largest Lyapunov exponent Ama associated with them:

(a) periodic solution c=1.00, h=0.80, p=1.00,

=

^1.00,

X0=0.80, R=0.65; (b) chaotic solution c=1.00, h=0.10, p 1.00,co=1.00,X0^{=0.50, R=0.65.}

(a) Lmax ks^=0.000

1.20

0.00

(b)

x,\

.60 0.00

,max k =0.083

0.50

2.50

X

0.00

-2.00

-2.00 0.00 0.50

Xl FIGURE 8 The phase portraits showing the solutions of Duffing oscillator with impacts (single buffer Eq. (29)) and the largest Lyapunov exponent Ama ^{associated with them:}

(a) periodic solution c=1.00, h=0.80, p=1.00, co=1.00, X0=0.50, R=0.65; (b) chaotic solution c=1.00, h=0.10, p 1.00,c=1.00,X0^{=0.50, R=0.65.}

5. CONCLUSIONS

Thepresentmethod allowsusto estimatethe value of the largest

Lyapunov

exponent of nonsmooth systems based onthe propertiesof chaos synchro- nization. The developed method will be useful in quantifying, predictingandunderstanding chaosin nonsmooth discontinuous systems for which the straightforwardcalculationofthe

Lyapunov

expo- nents is notpossible.Theapproach presentedin this papercanbe generalizedtothehigherdimensional

(n > 3)

systems. This method can be applied for bothnumericalandexperimentalestimationofthe largest

Lyapunov

exponent.Inthe firstcaseonehas

to know the equation of motion. In the experi- mental case twoexamples ofthe system haveto be created andcoupled together.

References

[1] V.I.Oseledec:Amultiplicative ergodic theorem:Lyapunov characteristic numbers for dynamical systems, Trans.

MoscowMath. Soc. 19, 197-231 (1968).

[2] G. Benettin, L. Galgani, A. Giorgilli andJ.-M. Strelcyn:

Lyapunov exponents for smooth dynamical systems and Hamiltonian systems;amethod for computing all ofthem, Part I: Theory; Part II:Numericalapplication,Meccanica15, 9-20;21-30(1980).

[3] J.B. Wolf, H.L.Swift, Swinney andJ.A. Vastano:Determin- ing Lyapunovexponents fromatime series,Physica D 16, 285-317(1985).

[4] P.M/iller:CalculationofLyapunovexponents fordynamical systems with discontinuities, Chaos, Solitons and Fractals 5(9),1671-1681(1995).

[5] H.Fujisaka andT.Yamada:Stability theoryof synchronized motion incoupled-oscillator systems,ProgressofTheoretical Physics69(1),32-47(1983).

[6] L.M. Pecoraand T.L. Carroll: Synchronization ofchaos, PhysicalReviewLetters 64,221-224(1990).

[7] K.PoppandP. Stelter: Nonlinear oscillations of structures inducedby dry friction,in: Non-linear DynamicsinEngineer- ingSystems,Ed.W. Schiehlen, Springer, NewYork(1990).