A global analysis of a non-linear system under parametric excitation

Journal of Sound and Vibration (1984) 97(3), 399-427

A GLOBAL ANALYSIS OF A NON-LINEAR SYSTEM UNDER PARAMETRIC EXCITATION R. s. GurrALut

AND

c.

s.

Hsu

Department of Mechanical Engineering, University of California, Berkeley, California 94720, U.S.A. (Received

17 June

1983, and in revisedform

10 February 1984)

A strongly non-linear mechanical system consisting of a rigid damped bar subjected to a periodic parametric excitation is treated here in an exact manner. The emphasis is on the global behavior of this system which is carried out by using the point mapping and the cell-to-cell mapping methods. Even though the mechanical system is a simple one, yet it has a very complex global behavior. It is shown here that the newly developed theory of cell-to-cell mappings offers a tremendous advantage in obtaining the global domains of attraction of strongly non-linear dynamical systems. 1. INTRODUCTION

The steady state response of a non-linear dynamical system subjected to impulsive periodic parametric excitation is studied here by using the point mapping and cell mapping techniques. Impulsive periodically excited systems belong to a class for which one can carry out an analysis of the equilibrium states and periodic states in a genera1 and systematic manner without approximation; sometimes this analysis results in a closed form analytical stability criteria [ 1,2]. Once all the equilibrium states and periodic states with their local stability character become available, then one can proceed to study the global domains of stability associated with each of the stable states. However, finding the locations of all the equilibrium states and periodic states of a non-linear system is usually a very difficult task and the determination of global stability regions can often be very tedious. In this regard, a recently developed method, known as cell-to-cell mapping [3,4] makes it possible to obtain the stable equilibrium states, periodic states and their global regions of attraction in a simple manner without any prior knowledge of the location and distribution of the equilibrium states and periodic states. In section 2, we explain some basic terminology regarding point mapping and cell mapping dynamical systems with due emphasis on the methods of determining the global behavior of nonlinear systems. The mechanical system to be studied here consists of a rigid bar subjected to two impulsive actions in a period. In section 3, an exact point mapping for this system is derived from the ordinary differential equation of motion. In section 4, equilibrium states and a few periodic states bifurcated from the trivial equilibrium states are obtained and their local stability determined. Finally, in sections 5 and 6 we present and compare the global domains of attraction obtained by point mapping and cell mapping methods. Even though the mechanical system studied here is simple enough, yet it has a rich variety of solutions and its global behavior is rather complicated. The new ideas developed in the cell mapping approach is seen to be very advantageous in analyzing non-linear dynamical systems. To conserve space only the most necessary references are given in the reference list of this paper. For other important papers on point mapping and global analysis of strongly non-linear systems the reader is referred to the reference lists given in references [l-9]. t Now at the Department of Mechanical Los Angeles, California 90089, USA.

Engineering,

University

of Southern

California, University Park,

399 0022-460X/84/230399+29

$03.00/O

0

1984 Academic

Press Inc. (London)

Limited

400

R.

2. PRELIMINARIES

S.

GUTI‘ALU

AND

<.

S. HSU

REGARDING POINT MAPPING DYNAMICAL SYSTEMS

AND

CELL

MAPPING

We briefly describe here some methods currently used to study local and global behavior of non-linear dynamical systems. We first consider a differential dynamical system described by a system of ordinary differential equations, F(t+~,x)=F(t,x),

k(t) = F(t, x(t)),

(1)

function assumed to be where x is an N-dimensional state vector and F is a real-valued periodic in time t of period T. It is of interest to study the equilibrium states and periodic states of equations (1) and their local stability character; if the system depends on a parameter, then we would like to know how these states change as a function of this parameter. Finally, we would like to study the global behavior of the system by determining the domains of attraction of asymptotically stable equilibrium states and periodic states. When F is non-linear, for a periodic system like equations (1) it is advantageous to consider its equivalent discrete-time formulation. Equations (I) can be integrated over one period of time to generate a functional relationship between the state of the system at the beginning of one period to its state at the end of the same period. This results in a difference dynamical system, or a point mapping G, described by a system of difference equations x(n+l)=G(x(n)),

(2)

where n is an integer. Once one obtains a solution of equation (2), one can always go back to system (1) to determine its continuous time history if needed. It may be difficult to obtain G analytically from equations (1); however, it can be obtained in principle by either numerical integration or by using difference approximations. For certain types of systems, notably impulsively excited systems [l], one can determine G exactly. The simplicity of expressing the systems in the discrete form (2), which is particularly suited to digital computers, helps one to study the local and global behavior of the system. We assume that G is continuous and its inverse W exists. Usually G is referred to as the forward point mapping and W as the backward point mapping: i.e., x(n) = W(x(n + I)). In order to study the point mapping system, one defines a discrete trajectory starting with an initial state x(0) as the sequence of states obtained by the mapping G: i.e., G’(x(O)) withj=O, 1,2,. . . . An equilibrium state, or a jixed point x* of system (2) is given by x* = G(x*) and a periodic state of period K (denoted a P-K distinct states x*(k), k = 1,2, . . _ , K such that x*(k+

l)=G”(x*(l)),

(3) solution)

k=1,2,...,K-1;

is given by a sequence

~*(l)=G~(x*(l)).

of K

(4)

An equilibrium state is also referred to as a P - 1 solution. If one defines the Jacobian of G(x) by J(G(x), x), then the stability of a P-K solution is ascertained by studying the eigenvalues of the stability matrix H: H = (aGK/ax)x=x*m

= [J(GK(xL

= [J(G(xL x)I,=,yK)[J(G(x),

x)lx=x~,, xlx=xv--lj

. . . [J(G(x), x)l,=.-c~).

(5)

The reader is directed to reference [1] for details regarding the local stability analysis of periodic solutions, their classification for second-order systems, and the conditions for local stability and bifurcation of periodic solutions. We merely state here that, for asymptotic stability, all the eigenvalues of H have their absolute value less than one; for

GLOBAL

ANALYSIS

the second order systems the conditions

IdetH( < 1,

OF

A NON-LINEAR

SYSTEM

401

for asymptotic stability become

l-traceH+detH>O,

l+traceH+detH>O.

(6)

Once the various P-K solutions of (2) are obtained and their local stability determined, then it is important to ascertain how a trajectory starting from x(0) is attracted by a stable P-K solution and if more than one P-K solutions are stable which one would attract it. It may also be important to study how this trajectory is repelled from an unstable P-K solution. The domains of attraction can then be determined by ascertaining the behavior of the total collection of trajectories in the state space. In what follows, we briefly identify different methods of determining the global stability regions of asymptotically stable periodic solutions of second-order systems: for details the reader is directed to references [l] and [2]. The direct method or iterative method consists of choosing a set of points in the phase plane and simply applying G consecutively to each of these points until either the trajectory is seen to approach a P-K solution or becomes unbounded. Evidently, only a limited number of points can be chosen and a large amount of computations is involved. The second method depends on the fact that for two-dimensional systems the separatrices of a non-linear system at a saddle point are tangent to the eigendirections of the locally linearized system. A point x(O) is chosen to be close to the saddle point lying on the straight line whose direction coincides with one of the eigenvectors. To ’ applied to x(O) if the separatrix generate a discrete trajectory, forward mapping GK IS leaves the saddle point or backward mapping WK is applied if the separatrix approaches the saddle point. Even though the unstable invariant manifold so obtained often acts as a boundary between two separate domains of attraction, the complicated windings of the separatrices (due to the presence of many unstable solutions or heteroclinic and homoclinic points) makes it very difficult to delineate the global stability regions. The third method consists of expanding a small region around an asymptotically stable P-K solution. One chooses a closed curve around a P-K point (a member of the P-K solution) such that when WK mapping applied once to this curve results in another curve which lies outside of the original curve. By repeated application of backward mapping, it is possible to expand the domain of attraction to a desired level. The fourth method of determination of global behavior depends on a new way of looking at system (2), which we now explain in more detail. A simple and practically realistic approach of solving equations ( 1) and (2) is as follows. From the considerations of round-off errors in numerical computations and limited resolution of physical measurements, it is appropriate to consider the state space not as a continuum but as a collection of intervals or cells. That is, the state variables are no longer continuous, but, instead, they consist of a discrete set of values. Such a description of dynamical systems is referred to as cell-to-cell mapping systems described by

Z(n+I)=C(Z(n)),

(7)

where Z(n) is the state defined by an iV-tuple of integers and C is the cell mapping which maps a set of integers to a set of integers. For a basic description and study of cell mapping the reader is referred to reference [3]. Here we consider cell mappings derived from point mapping (2). In order to construct a cell mapping from equations (2), the state space is first divided into a large number of cells, each cell having a size hi, i = N, in the direction of the xi axis such that the interval Zi is defined to be an 1,2,..., integer satisfying (zi-t)hi~xi<(Zi+~)hi,

i=1,2

,...,

N.

(8)

402

R. S. GUI-I-ALU

Let x4 be the center process:

point

AND

C. S. HSU

of a cell. Then the cell mapping

x:‘(n) = hiZi( n),

is generated

by the following

x,(n+l)=G,(~~(n)),

(9,lO)

Z~(~+l)~C,(Z(n))~Int[(l/hi)(x~(n+l))+~],

(II)

where Int [u], for any real number u, represents the largest integer such that Int [u]~ u. Then, in a manner analogous to P-K solutions, one can define P-K cells. The mapping C now contains essentially the information regarding the global behavior of the original dynamical system. In order to unravel this information, an algorithm has been devised [4] which effectively locates the equilibrium cells, periodic cells and the global domains of attraction all at once. In general, the cell mapping C as defined by equation (11) replaces the stable periodic solutions of equations (2) by periodic cells; the unstable solutions of equations (2) may not be recovered. Each set of periodic cells, in general, may consist of a group of “true periodic cells” which replace the original periodic points of the point mapping and groups of “pseudo-periodic cells” which surround the true periodic cells. These pseudo-periodic cells are the result of discretization process; the state space occupied by them has been shown [3] to decrease as the cell sizes are made smaller. In using equation (1 I), one views the associated cell mapping as an approximation of the point mapping; the degree of approximation can be improved by simply reducing hi.

3. AN IMPULSIVE

PERIODIC

PARAMETRICALLY

EXCITED

DYNAMICAL

SYSTEM

We briefly describe here how to obtain the point mapping (also known as the Poincare map) of a parametrically excited system as shown in Figure 1. The mechanical system

Figure

I. A rigid damped

bar with its support

subjected

to a periodic

displacement

f(7)

consists of a rigid bar in the plane with a moving support at A which has a rotational damper with damping coefficient c. The support is subjected to a periodic displacement f(7) of period T,,. Let C be the center of mass of the bar, IA its moment of inertia about A, and 8 its angular displacement measured from the reference direction xl. The gravitational field is not involved in this problem. The equation of motion of the bar is given by IA d20/dr2+

c de/dr-

JA(d2f/dr2)

sin 8 = 0,

JA = mlo,

(12)

GLOBAL

ANALYSIS

OF A NON-LINEAR

403

SYSTEM

where m is the mass of the bar, and I,, the distance of the center of mass from A. A saw-tooth form of the support motion is assumed,

;7:::2:.

(1-4r/rcl)p, fcT)={(-3+4r,To)p,

Introducing

the non-dimensional

(13)

quantities 2p = cro/ IA,

t = r/ 70,

(14)

and substituting (13) in (12), one has the following non-dimensional of motion of the bar:

form of the equation

2

S+2p$+

1

(Y ; G(t-m)-a In=-cc

f

In=-x

@t--$-m)

I

sin0=0,

(15)

with LY= 4p~~J*/ 1,. In the differential equation (15), the parametric excitation is now of period one and consists of two equal and opposite actionst within one period. The first impulsive action occurs at the instant t = 0 with strength LYsin 8 and the second at t = f with strength -_(y sin 8. Let

a=!zc 2P ’

b = eep.

(16)

Consider the first period from t = 0 to t = 1 as a typical period. The angular displacement xi of the bar being continuous at t = 0 and t = i, one obtains the following jump conditions satisfied by the angular velocity of the bar: x,(0+) -x,(OP) = --(Ysin x,(O),

x2($+) - x,(f-) = +a sin x,(i).

(17)

Between the impulses the non-linear system ( 15) reduces to a linear system whose fundamental matrix can be readily found. The continuous trajectory in the phase plane for this linear system is x*(t)

-x,(0)1.

=X2(0)-W%(t)

(18)

Then, upon starting with the initial condition x(0) corresponding to the instant t = Ome just before the action of the first impulse and by successive use of the fundamental matrix and jump conditions, x( 1) can then be determined and the general period-to-period mapping corresponding to equation (15) can be established. The reader is referred to reference [5] where a very general procedure is established to obtain the point mapping of a general dynamical system of the kind (15), or to reference [2] where the above procedure is followed. The forward mapping G is found to be x,(n+l)=G,(x(n))~x,(n)+a(l+b)X(n)+aasin[x,(n)+aX(n)], ~~(n+l)=G,(x(n))~b~X(n)+cybsin[x,(n)+cuX(n)], X(n)=-cusinx,(n)+x,(n). Here n is an integer. The above equations transformation Ly’= -Ly, t For the cases reference [2].

where

there

x,(n) = t=+ are two

arbitrary

(19)

can be seen to be invariant

1h+Y,(nL impulses

x2(n)

within

one

=

period,

under the

Yz(nL the

reader

is referred

to

404

R. S. GUT-I-ALU

AND

C. S. HSU

where k is an integer. Hence one can assume with no loss of generality Jacobian matrix J(G(x), x) of the mapping G is given by

J(GW, 4 =

1-~a(l+b)cosx,(n)+aaY(n)cosQ(n)

a(l+b)+CYa2cosQ(n)

-ab’cosx,(n)+abY(n)

[

that cy > 0. The

cos Q(n)

Y(n)=l-cuacosx,(n),

b*+aab

cos Q(n)

O(n) = x,(n)+aX(n), det J = b2.

traceJ=l+b*-czu(l+b)cosx,(n)+cya(b+Y(n))cosQ(n), For /L > 0, det J = e-2cI < 1 and one has an area contraction mapping W can be easily obtained from equations (19) as x,(n)=

x,(n)=

mapping.

a(l+b) ----.r,(n+l)+ysin bZ

W,(x)~x,(n+l)-

W2(x)P$x2(n+l)-fsin

(20)

The backward

T(n+l),

T(n+l) u(l+b) ----x,(n+l)+ysin b2

T(n+l)

,

(21)

T(n+l)=x,(n+l)-Zx,(n+1). If there is no damping,

1’

then p =O,

For this case, the forward

mapping

b= 1.

u =;,

G simplifies

(22)

to

x,(n+l)=G,(x)~x,(n)+X(n)+(cu/2)sin[x,(n)+~X(n)], X(n)=-a

x,(n+l)=G,(x)~X(n)+(~sin[x,(n)+~X(n)],

sinx,(n)+x,(n), (23)

and the Jacobian

matrix

J(W4,x) =

1 -a

1

J(G(x),

x) becomes

cos x,(n)+(a/2)Y(n)

cos Q(n)

-a cos x,(n)+aY(n)

Y(n) = 1 -(a/2)

cos Q(n)

cos x,(n),

cos Q(n) cos Q(n) 3 ’

O(n) =x,(n)+fWn), detJ=

traceJ=2-~~0s~,(n)+(~/2)(1+Y(n))cosQ(n),

4. PERIODIC

1 +(a(/4) lf(a/2)

1.

(24)

SOLUTIONS

only those In this section we briefly list various periodic solutions of (19), highlighting relevant for global analysis to be presented in sections 5 and 6. The reader is directed to reference [2] for an elaborate discussion on periodic solutions of equation (19). 4.1. P-l

SOLUTIONS

Since the mechanical system under consideration is elastically unrestrained, one can modify the definition of P-l solutions as follows to permit determination of physical equilibrium states which have their successive positions differing by modulus of 27~: x,(n+l)=x,(n)+2m?r,

x2(n+l)=x*(n),

(25)

GLOBAL

ANALYSIS

OF A NON-LINEAR

SYSTEM

405

where m is an integer. One notes that in equations (25), the case m = 0 corresponds to the regular P-l solutions. For m # 0, we will refer to the P-l solutions obtained by using equations (25) as advancing P-1 solutions. The meaning of an advancing P-l solution is that the angular displacement x,(n) of the bar changes by an amount equal to 2m7r in every period. In other words, physically the bar undergoes m complete rotations at the end of each period. Substituting equations (25) into equations (19) gives the following equations to determine the P- 1 solutions x*( 1): asinx~(1)-asin{x~(l)-[aasinx~(1)-2m~b]/(l+b)}+4m~~=O, xT(2)=xT(1)+2m7r,

x~(l)=[b/a(l+b)][cuasinx~(1)+2mlr].

(26)

First one can consider the regular P- 1 solutions for which m = 0. The stability conditions can be obtained directly from equation (6). The trivial solution x* = (/CT,0) is found to be unstable for a2 > 2( 1+ b2)/a2.

stable for (Y’< 2( 1 + b2)/ a’,

(27)

Besides x7( 1) = kr, there are two other types of solution which one can designate as (P- 1, I, n) and (P- 1, ZZ,n) and they are given as follows: (P-l,

z, n): (P-l,

x:(1)=2nrrb/a,

asinx:(l)=2nr(l+b)/a,

ZZ,n):

(28)

LYsin x7( 1) = ~[2x:(l)-(2n+l)rr], *f(l)=~[2x:(1)-(2n+l)*],

(29)

where n is an integer. The Type I solutions exist for LY> 2n7r( 1 + b)/ a and are found to be stable for a’<

unstable for a2> 16( 1 +4n27r2).

16( 1+4n27r2),

(30)

Similarly, the Type II solutions for CY > 0 are stableforO
cosxT(1)<2(1+b)/a,

unstableforacosx~(l)
or

crcosx~(l)>2(1+b)/a.

(31)

For the case ZL= 0, the solutions of equation (26) are designated as (P-1, III, m, n) and (P- 1, IV, m, n) and are as follows: (P-l,

ZZZ,m, n):

(Y sin x7( 1) = 4nT,

xT(1)=2(m+n)7r, (P-l,

IV, m, n):

x?(2)=xf(l)+2m7r,

m, n integers, (m + n) even;

asinxT(1)=8xT(l)-4n?~,

xT(1)=4x:(1)+2(m-n)7r,

(32)

x7(2) =x7( 1)+2mr, (m+n)odd.

(33)

The Type III solutions exist for a > 4nrr and are stable for a2 cos’ x7( 1) < 16 (or 4n7r < a < Al+

n2rr2),

unstable for (Y*cos’ xt( 1) > 16.

(34)

The Type IV solutions are found to be stable for 0 < a cos xf( 1) < 8, unstable for CYcos xf( 1) < 0 or a cos x7( 1) > 8. (35)

406

R. S SUI-I-ALU

AND

C. S. HSlJ

For the case /L Z 0, one has to solve equation (26) using as a guide the results obtained for the case p = 0. The solution curves for m = * 1 and p = 0.02~ are shown in Figure 2. Here the convention is used that a solid line represents the stable part of a solution curve and a dotted line represents the unstable part. The branch A,B, corresponds to m = +I while the branch C,D, corresponds to m = -1. As the value of p is descreased to zero, these two branches coalesce into one branch and the branches represented by E, and F, to x7( 1) = -7r, x7( 1) = 0, and x7( 1) = +n, and a segment of the xT( 1) axis. F,

-0

25

-1

Figure

4.2. P-2

2. Advancing

P- 1 solution

curves for the case p = O.O27r, m = *I

SOLUTIONS

P-2 solutions are obtained by solving two vector equations, x(2) = G(x( I)), x( 1) = G(x(2)), which consist of four equations in four variables x1( I), xq( l), x,(2) and x,(2). Eliminating x2( 1) and x,(2) gives two coupled non-linear algebraic equations to determine x7( 1) and x?(2). The resulting equations are -bx,(l)+x,(2)+

(1+b2)[-~,(1)+~,(2)]+cyaX,2+~a(1-b)sin

+x,(2)+x,(l)+ l-b

(1+b2)[-x,(2)+x,(l)]+aaX,,+aa(l-b)sin

~a

x

~a l-b2

x

l-b2

l-b

(36) x2(1)=-

b a(1 -b) 1

xI(l)-xI(2)-5X12

I

b

x2(2) = ~ a(l-b) X,* = b sin x,(l)

-sin

x,(2) -x,(l) x,(2),

-$-x2,

1

,

,

X,, = b sin x,(2) -sin

x,( 1).

GLOBAL

ANALYSIS

OF

A NON-LINEAR

The equations for the case p = 0 can be obtained (that is, a -+ ‘2,b + 1) and are

407

SYSTEM

from equations

(36) by letting

CL+ 0

-x,(l)+x,(2)+$r{sinx,(l)-sinx,(2)}=0, -x,(2)+x,(l)+~cy{sinx,(2)--sinx,(l)}=O, x,(1)=$

x,(2) = fa sin x,(2).

sinx,(l),

(37)

Since equations (36) are complicated to solve for xT( 1) and x;(2) and may have more than one solution, the first step is to investigate the possibility of P-2 solutions bifurcated from the trivia1 P-l solutions of xT( 1) = ka. The resulting solutions are as follows: (P-2,

n):

2(1+b2)t-(-l)“Lua(l+b)sinE (l+b)t-(-l)“c~asin[

+(-l)“aa(l-b)sin

x7(2) = -5-t

xT(l)=(+?Wr,

n?T,

x~(1)=-x~(2)=[b/a(l-b)][2~-(-l)“aa The corresponding results limiting case as p + 0: (Ysint=(-1)“4&

The stability

for the case p=

xT(l)=[t-nn,

matrix

H associated H =

The P-2 solutions

sint].

0 can be obtained

x7(2) = -c+

with the P-2 solutions

(38)

from equations

n7r,

of equations

(39) for p = 0 are found

(38) in the

xT( 1) = -x?(2)

[J(Wd, x)lx=,*,z,[J(G(x),(x)1x=x-c, ).

given by equations

= 0,

l-b

= 25. (39)

(38) is given by (40)

to be

stableforO4.

(41)

For p # 0, one can obtain the (P-2, n) solutions by solving the non-linear algebraic equation (38), using the case p = 0 as a guide. These solutions are shown in Figure 3(a) for p = 0.1 m and n = 0, f 1 with the continuous trajectories in the phase plane corresponding to (Y= (Y, shown in Figure 3(b). As shown in the figure, there are three P-2 solutions for a given value of (Y consisting of the groups AE, A,E, and A2Ez and together they contain six P-2 points. The solution branch AE is a (P-2,0) solution bifurcating from the trivia1 P- 1 solution of xf( 1) = 0. The solution branches A,E, and AzE2 are, respectively, (P-2,+1) and (P-2,-1) solutions bifurcating from xT(l)=+r and x:(1)=-z The (P-2, n) solutions start at (Y= 4.081 497 and are stable up to LY= 5.054 97 1 3 and become unstable thereafter. If the initial state for a given cy is a point on the branch A, A, or AZ, then it is mapped at the end of the first period to a corresponding point on the branch E, E, or EZ, respectively, which returns to its original state at the end of the second period (or after the second mapping). The same process repeates if, instead, one chooses a point on the branch E, E, or E2. The branches A, and E, for Jo f 0 merge into a single branch x7( 1) = +r for p = 0 while the branches A2 and E2 merge into x:( 1) = -TL It is interesting to note from Figure 3(b) that the continuous trajectories in the phase plane for the three sets of P-2 solutions are identical except for the fact that the periodic states of the system for the solution branches A,E, and A2E2 correspond to half-period states for the solution branch AE.

408

R. S. C’cl-l-ALU

AND (‘. S. HSU

__--_------_------

____________________~~ (-AlE, -=z-__-__----_______x

1.0

(a:

t-

--._ -.---_ ----____ -E ,,A2

-10

/ 2

0

4a,

6

a

10

12

14

(b) F

H

H1

61 Cl

G2

>A

4

I

1

-1.0

-050

1

0.0

1

0.50

J

10

x,/lT

Figure 3. (a) (P-2, n) solution curves in the (a, xi) planes for the case p = 0.1 T, n = 0, *l ; (b) trajectories in the phase plane for (Y, = 4.5.

4.3. P-3

continuous

SOLUTIONS

In order to determine P-3 solutions one must solve three vector equations, x(2) = G(x( 111,

x(3) = W(2)),

x( 1) = %(3)),

consisting of six scalar equations in unknowns x(l), x(2) and x(3). Again, eliminating x2(1), x2(2) and x2(3), one obtains the following three coupled non-linear algebraic equations to determine x,(l), x,(2) and x,(3): - d2Xlz3 + (~a sin (d3ZLz3- &X1& = 0,

x,(l)-x1(2)+d,Y,,,

yz3, - d2Xz3, + (YQsin (d3Zz3, - d4Xz3,) = 0,

x,(2)-x,(3)+d,

d, Ys,z -dzX,,2+aa

x,(3) -x,(l)+ xX(1) = 4 Ym + dJm,

sin (d3ZJj2-d4X3,J

x2(2) = 4 Y,,, + d.Jn,,

=O,

x2(3) = d, Y3,z+ &X3,2,

Xiz3 = sin x,( 1) + b* sin x,(2) - b sin x,(3), Y12~=(l-b)x,(l)+bx,(2)-x,(3), d

I

Z,,,=x,(l)+@x,(2)--x,(3),

=b(l+N* l+b3

+df-

1 + b”

* d

=

5

l+b

=41+w

d

’

l+b3

b(l+b) a(l+b’)’

’

d3 = 1 + b3’ ab de=-1+ b3’

(42)

GLOBAL

ANALYSIS

OF

A NON-LINEAR

409

SYSTEM

For the case p = 0, the above equations reduce to x,(l)+x,(2)-2~~(3)-~aX,~~+~~

sin{x,(l)+x,(2)-x,(3)-$aX,,,}=O,

x,(2)+~,(3)-2~,(1)-~(1X~~,+~cz

sin{x,(2)+x,(3)-x,(1)-f~X,,,}=0,

x1(3)+x,(1)-2x,(2)-~cyX,,,+~~

sin{x,(3)+x,(l)-~,(2)-~~X,,~}=O, x2(2)=2[x1(3)-x1(1)i+&X,23,

x2( 1) = 2[x,(2)-x,(3)l+~~Xm,

X,,,=sinx,(l)+sinx,(2)-sinx,(3).

x*(3)=2[x,(l)-x,(2)1+t~X*,I,

0 75

_zr:Z-

Al

__z:=

D5

-,: .;;: ,* ,*

/;-

05

,‘,a /‘,’ ,/.’

:

\k

z -

A6 -------________ ______________-____--------____---

o-o

*<

(43)

\\

A2

‘:\

I

‘C\

‘;\

-0-5

‘Q.

‘. .‘.;.

.-:z>_

-0.75

-1.0 0

4

2

I

1

a,

6

--::__ -_3_ Dz A4

~

/

6

10

12

14

a

(b) 1.10’

1 05

,,:y====;=+__ / I ,,‘,/’ Ez / / ,’ I : 86

.

k 2

100 Y

>,-

b 095f

0

,_-________-_______--__--_-_~ _____----__-fl ,,_* ,’

_I’

1

E3 04

,’

\

901 0

-.m 85

“:‘______~_~~__---------E1 2

4

a,

6

8

10

12

14

a Figure 4. (a) P-3 solution the line xt( 1) = T.

curves in the (a, XT) plane for the case /* = 0.00277; (b) P-3 solution

curves near

410

R.

S.

GUT-l-ALU

AND

<‘.

S.

HSU

For the case p = 0, investigation of the possibility of P-3 solutions bifurcated the trivial P-l solutions of xT( 1) = kr, results in the following solutions: (P-3, 2[-5,

z, n):

+&-&a[-sin

from

5, - & - $cz sin 5, +l,a sin [t, - :(Y sin r,] = 0, 5, +2 sin &]+$(Y sin (-5,

+2&--a(~[-sin

El+2 sin &I} = 0, (44)

xT(l)=5,+2n7r, xT(l)=iasin[,,

x7(2) =x7(3) x2*(2) = 2[-5,

X;(3) = 2[5, -&+&[-sin (P-3,

ZZ, n): xT(l)=[+n?T,

x;(l)=-2(+(-l)“(~sin&

asin.$+$

= &;+2n7r,

+ &]+:(Y sin [,,

5, +2 sin &I,

sin[-25+(-l)“$sint]=(-1)“3.$, x772) = -5 + n?T, x?(2) = -25,

x7(3) = m-r, x;(3)

=45-(-1)“a

sin 5.

(45)

One notices from equation (45) that only two cases of even and odd integer values of be considered. P-3 solutions for p # 0 are determined by solving the first three of the non-linear algebraic equations (42). Note that a bifurcation from trivial P-l solutions into P-3 solutions (or from P-3 to P-3M, M integer) is not possible as indicated in reference [ 11. The solutions obtained may again be designated as (P-3, Z, n) and (P-3, ZZ, n) as in the case of p = 0. Figure 4(a) shows P-3 solution curves near the xT( 1) = 0 axis in the ((Y, x7( 1)) plane. Figure 4(b) shows the P-3 solutions near xT( 1) = n to an enlarged scale. In this figure the mirror images with respect to the line x7( 1) = 7~ are not shown. Such solutions also exist near x7( 1) = --n and are not shown. Totally there are twelve P-3 solutions (thirty-six P-3 points), four each near xf( 1) = 0, xT( 1) = +n. and xT( 1) = -r in the ((Y, x7( 1)) plane. Six of these solutions are asymptotically stable from (Y= 3.640 839 to a = 4.022 089, two each near xT( 1) = 0 and xT( 1) = +7r, respectively. Solutions labeled D,D2D3, D4D5D6, and EIE2E3 are stable, while those labeled A1A2A3, &A,& and B4B5B6 are unstable. When p is decreased to zero, the branch A, merges with the branch A3, As with &, D, with DZ, DS with D,, and D, and D6 coincide with the xT( 1) = 0 axis. For other specific periodic solutions not discussed here the reader is directed to reference [2]. n need

5. DETERMINATION

OF DOMAINS OF ATTRACTION MAPPING TECHNIQUES

BY POINT

This section is concerned with the determination of the domains of attraction associated with asymptotically stable P-K solutions of the point mapping system (19) by using some of the methods described in section 2. Unless otherwise mentioned, we will obtain the global regions of attraction in the phase cylinder (instead of the phase plane) by taking xi modulo 27r. This is made possible because the dynamical system (19) is elastically unrestrained. First we. will consider a case for which only P-l solutions of equation (19) exist by taking the damping coefficient and the strength of the impulse to be p = 0.027r and (Y= 1 *O. Table 1 gives details on various P-l solutions that exist, their location in the phase plane and their local stability character. There are five stable and four unstable P-l solutions for the above set of parameter values. In the figures to follow, stable solutions are designated by the symbol “+” and the unstable solutions by “x”. The result of the method of tracing of separatrices applied to the unstable P- 1 solution x2 is shown

GLOBAL

ANALYSIS

OF

A NON-LINEAR

TABLE

Location in thephaseplane

P-K

P-l

and local stability characteter of P- 1 solutions of ( 19) for the case /A= o.o027r, (Y = 1.0 Eigenvalues (A,, Al)

(XI, xz)

P- I

XT = (0.0) xT=(+v,O)

P-l

x$ = (+0.539

x)=(+r,O) 973~)

P-I

x,* = (-0.129

1807r, +1,937

826~)

XT = -xh* Advancing

P- 1

xf = (-0.870 x9* = -xR*

819q

+I.937

826~r)

Stability character

(0.823 529kO.451

3451’)

Stable spiral

(0.823 52910.451

3451’)

Stable spiral

3451)

Stable spiral

(0,823 529kO.451 4707q +0.152

XT =-x4* Advancing

411

1

Location

solution

SYSTEM

( 1.548

390,0.469

57)

Saddle

(I.548

390,0.569

57)

Saddle

I I

(0.176609*0.9223291’)

Stable spiral

(0.176609iO.9223291’)

Stable spiral

13.559 1554,0.247

787)

Saddle

(3.559

787)

Saddle

i554,0,247

I I

in Figure 5(a). Here, the pair of separatrices leaving the saddle point is seen to approach the stable P- 1 solutions XT, xf (in the phase cylinder XT and XT are identical). The other pair of separatrices approaching the saddle point traverse all over the phase cylinder and do not interesect again, as they are expected not to in general. For the unstable P-l solution XT, which is the counterpart of x2, the result of tracing the separatrices is a symmetric reproduction with respect to the origin of Figure 5(a). For the unstable advancing P-l solution xg, the result is shown separately in Figure 5(b). Here, we emphasize that the set of two pairs of separatrices shown in Figure 5(a) is entirely different from the set shown in Figure 5(b). Those shown in Figure 5(a) are from x,* while those from Figure 5(b) are from xg. We also wish to point out that when interpreting the results shown in Figures 5(a) and (b), one should keep in mind that we have used mod 271. for the x,-co-ordinate and that xz is an advancing type P-l solution. In Figure 5(b), the pair of separatrices leaving the saddle point approaches the stable advancing P-l solution x2 and the stable P-l solutions XT and x z. When Figures 5(a) and (b) are combined together with their symmetric parts with respect to the origin, then the regions between appropriate separatrices give the global asymptotic stability regions for the P- 1 solutions XT, x:, XT, xt and XT. The plots of Figures 5(a) and (b) clearly indicate the complex nature of the fine divisions of domains of attraction far removed from the stable singular points. To illustrate these fine divisions of domains of attraction one can pick three points very close to each other in the phase cylinder given by x, = (-0.35~7, -2.4630~), x2 = (-0*35%-, -2.4635~) and x3 = (-0.357r, -2.464Orr). These states are indicated by the symbol “+” in Figures 6(a)-(c), respectively, where the evolution of the dynamical system starting from these states is also shown by dots. It is observed that the initial state x, approaches the stable advancing P- 1 solution XT, x2 approaches XT and xj approaches XT. We now consider a case for which P-l and P-3 solutions of system ( 19) exist for the parameter values p = 0.0027r and (Y= 3.75. Table 2 lists various P- 1 and P-3 solutions together with their local stability character. Many advancing P-l solutions exist and higher and higher angular velocities are required to reach them. Only a few of them have been indicated in Table 2. The method of expansion of region can be used, first in the phase plane and then in the phase cylinder, to obtain the domains of attractions associated with the various P-l and P-3 solutions. Figure 7 shows the domains of attraction of the P-l and P-3 solutions XT, xT8 and xTs, which are located near the origin. The procedure for obtaining such regions is described in reference [l]; to obtain the global stability regions associated with the P-3 solutions, the backward mapping W3 is needed. Figure 8

412

R. S.

-1.0

GUI-I-ALU

AND

C.

S.

HSU

I

-1.5 -2.0

I

-2.5 t -3.0 -3.5 -1.0

I -0.75

(b)

I -0.5

-0.25

0.0

0.25

0.5

0.75

1

x,/lr

Figure 5. Separatrices (Y= 1.0.

of the unstable P-l point (a) x2 and (b) x$ of the mapping (19) for the case p = 0.02~,

GLOBAL

ANALYSIS

OF A NON-LINEAR

SYSTEM

+

I

I

I

I

..,.

I

I

..

I

413

414

R. S. GUI-I-ALU TABLL

AND

(‘.S. HSU 2

Location in thephaseplane and local stability character of P- I and P-3 solutions of ( 19).for the case p = 0_0027r, (Y = 3.75 P-K solution

Location

P-l P-l

XT = (0,O) xT=(+r,O) xf=(-Tr,O)

P-l

xl = (+0.635 XT = -x4* x; = (-0.033 XT = -x6* x; = (-0.996 x9*= -xs* x& = (-0.992 XT, =-x,0 *

Advancing

P- 1

Advancing

P- 1

Advancing

P- 1

Advancing

P- 1

Advancing

P- 1

Advancing

P- 1

P-3

P-3

P-3

P-3

P-3

P-3

Eigenvalues

(x,, x2)

(A,, A,)

8281r, +0.541 606~) 5 1 I T, + I .993 723~) 649q

+ I ,937 826~)

388q

+4.014

199~r)

x& = (-0.02 I 03 1T, +3.960 694~) XT3= -x,2* xx = (-0.65 I 897w, +4.528 464n) x& = -x1** x&,=(-0.169463~,+3446429~) XT, = -x,6* x&( 1) = (+0.219 106q +0.426 5057~) x&(2)=(-0.010632q -0.128 3337~) x&(3) = (-0.200 3437r, -0.290 314~) x?Jj) = -xt(j),j = 1,2,3 x&( 1) = (+0.242 149q +0.417 029~) x;,(2) = (-0087 0357q -0.252 899~) x$,(3)=(-0,115593~,-0.125931~) G,(j) = -x$(j), j = 1,2,3 x;,( 1) = (+0.946 263q +0.329 700~) x;,(2)=(+1.0547297~, +0.087919~) x;,(3) = (+0.994 817~, -0.409 761 a)

Stability character

(-0.753 02 1 f 0.648 402i) (-0.753 021 kO.648 402i) (-0.753 021 10.648 402i) (5.490821,0.179848) (5.490821,0.179848) (-0.985959iO.124083i) (-0.985959*0.124083i) ( 12.856 59,0.007 681) (12.856 59,0.007 681) (-0.745 209 f 0.657 399i) (-0.745 209* 0.657 3991’) (-0.745 209 f 0.657 3991) (-0.745 209 f 0.657 399i) (5,472 803,0.180439) (5.472 803,0.180439) (5.472 803.0.180 439) (5.472 803.0.180 439) (0,750 274kO.632 5291’)

Stable Stable Stable Saddle Saddle Stable Stable Saddle Saddle Stable Stable Stable Stable Saddle Saddle Saddle Saddle Stable

(0.750 274kO.632 5291’) (1.879446,0.512386)

Stable spiral Saddle I

( I ,879 446,0.5

12 386)

(0.750 274kO.632

5291)

spiral spiral spiral I

I spiral spiral I I spiral spiral spiral spiral I

I I I spiral

Saddle I Stable spiral

(0.750 274jzO.632 5291’) (0.750 274* 0.632 5291’)

Stable spiral Stable spiral

x&(j) = -G(j), j = 1,2,3 x&(l) = (+0.947 571*, +0.068 995T) x&(2) = (+I.032 978~. +0.296 209~) x&(3)=(+1.039813~,-0.403403~)

(0.750 274iO.632 (I.879 446,0.512

5291) 386)

Stable spiral Saddle I

G(j) = -G(j),j = 1,2,3 x&(l)=(+l~O52429~, -0.0689957~) x&(2) = (+0.967 022q -0.296 209~) x&(3) = (+0.960 187q +0.403 4037~)

(1.879 446,0.512 (1.879 446,0.512

386) 386)

Saddle I Saddle I

G(/) = -x;k(J)J = 1,233

(1.879 446,0.512

386)

Saddle 1

x&(j) xf,( I) x&(2) x;J3)

= = = =

--XT,(j), (+I.053 (+0.945 (+I.005

j= I, 2,3 737q -0.329 700~) 271 T, -0.087 919~) 182q +0.409 761~)

shows the domains of attraction associated with the P-l and P-3 solutions x,*, x& and x& situated near the line x, = T. Figure 9 is the domain of attraction associated with the advancing P-l solution xz for m = - 1 for which we have used x, modulo 2n. In these figures, it can be noticed that the various domains of attraction wind about one another in a very intricate and complex manner. Further expansion of these regions indicate much more complex interactions of the various stability regions. In Figure 9 one can observe five regions around which fingers begin to appear and later wind in a complex manner.

GLOBAL

ANALYSIS

OF A NON-LINEAR

Figure 7. Domains of attraction associated with the Pnear the origin for the case p = 0$W2~, a = 3.75.

SYSTEM

1 and P-I solutions

of the mapping

415

(19) situated

---_-__

0.6, 05 c 0.4t 03 1 02’ 01;

I

k ;N

1

o-o

-0.1

-02,

t

t i

-03; -0.4 7 -051 -0.61.~ 0.92

0.94

0 96

0 98

1 00

1

02

1 04

106

'

x,/7?

Figure 8. Domains of attraction associated with the P- 1 and P-3 near the line x, = VTfor the case p = 0.002~, (I = 3.75.

solutions

of the mapping

(19) situated

416

R. S. GUT-l-ALU 3:

AND

C. S. HSU

,-

3.c l-

2.5 ,-

:

2.c

1.:

1.c )-

0.: i_ -( 1.f 5

Figure 9. Domains a =3.75,

p =0.002~,

of attraction m = -1.

,

I

-0.2

-0.4

associated

I

0.0

with the advancing

Figure 10. Domains of attraction associated here Figures 7-9 are combined into one graph

1

0.2

P-

with the various P-l in the phase cylinder.

,

0.4

06

1 solution of the mapping

and P-3 solutions

(19) for the case

of the mapping

(19);

GLOBAL

ANALYSIS

OF A NON-LINEAR

417

SYSTEM

This probably indicates the existence of other periodic solutions in those regions; later cell mapping will actually find a stable P-5 solution there. Figures 7-9 are combined into one plot in Figure 10, which is drawn in the phase cylinder, to indicate the relative sizes of the various domains of attraction. TABLE

3

Location in thephaseplane and local stability character of P- 1 and P-2 solutions of ( 19),for the casep=O.l7r, ff =4.5 P-K Solution

Location

x::= (0

P-l P-l

(Al, 4)

0)

(+L,

XT = 0) x:=(-T,,)

P-l

Stability character

Eigenvalues

(XI. x2)

-

(1.916 292,0.278 396) ( 1.916 292,0.278 396) (1.916292,0.278 396) (5.421 865,0.098 396) (5.421 8650.098 396) (0.730 391 i-0.004078i) (0.730391 *O.O04078i) ( 10.498 60,0.050 8 15) ( 10.498 60,0.050 8 15) (0.502 140*0~180 18Oi)

Saddle I Saddle I Saddle I Saddle I Saddle 1 Stable spiral Stable spiral Saddle I Saddle I Stable spiral

P-2

x2 = (+0.656 544q +0.532 958~) Xf = -x: x$ = (-0.144 543q +I.702 263n) XT = -x.$ xt = (-0.855 4577r, + I .702 263 P) x; =-x8* xTs( 1) = (+0.241920~, +0.381 479~)

P-2

x7,( I) = (+0.982

181q

-0.442

144~)

(0.502 140*0.180

18Oi)

Stable spiral

P-2

x?,(2) = (+1.017 x:,(1)=(-0.982 x:,(2) = (-1.017

SlSrr, +0442 1811r,+O442 818q -0.442

144rr) 144~) 144~)

(0.502 140*0.180 (0.502 140*0~180

18Oi) 18Oi)

Stable spiral Stable spiral

Advancing

P-l

Advancing

P-I

xX(2) = -x%(I)

Figure solutions

Il. Domains of attraction associated with the advancing of the mapping (19) for the case p = 0.1 rr, n = 4.5.

P-l

solution

(with

m = *I)

and the P-2

418

R. S. GUTTALU

AND

C. S. HSU

Finally, we consider a case when only P-l and P-2 solutions of (19) exist for the parameter values p =O*lr and (Y= 4.5. The various solutions and their local stability characters are indicated in Table 3. In this case XT, x; and XT are all unstable singular points. The global regions of P- 1 and P- 1 solutions x,*, x7, xTO, XT, and XT, in the phase cylinder are shown in Figure 11. Thicker curves in the figure indicate finely divided stability regions.

6. DETERMINATION OF DOMAINS OF ATTRACTION BY CELL MAPPING TECHNIQUES In this section the global regions of attraction of the point mapping G given by system (19) are determined by using the global cell mapping algorithm. The main purpose is to compare the domains of attraction obtained by the cell mapping method with those obtained by the point mapping methods of section 5. The cell mapping C(Z( n)) associated with the point mapping G(x( n)) is obtained by the following equations: Z(n+l)=C(Z(n)),

X(n)=-asin(h,Z,(n))+h,Z,(n).

(46)

Here hi and hz are the cell sizes in the xi and x2 directions, respectively. Since the domains of attraction will be obtained in the phase cylinder, the cell space covers the following range of xl values: -(l+l/N,,)~~x,<(l-l/N,,)n;

(47)

where N,, is the number of cells in the 2, direction. This also makes it possible for the center point of the cells (covering the lines xi = fr) to coincide with x, = frr. Since x* = (*n, 0) happen to be singular points of the mapping G in equation ( 19), the number of pseudo-singular cells of the mapping (46) corresponding to these singular points will be reduced by this cell division. It should be noted that in the phase cylinder, the point (+n; x2) is identified with (-rr, x2). of equation (19) exist, the first We first consider the case when only P- 1 solutions problem treated in section 5. For Figures 12(a)-(d) we have used N, = 101 00 regular cells covering the range of the phase plane -1.017r c xi s 0*997r, and -3.57r s x2 < 3.5 7~ with N,, = 100, NC, = 101, h, =O-O27r, hZ=0*069 3077r, and with the parameters of (46) given by p = 0*027r, LY= 1.0. Figure 12(a) shows the O-step domains of attraction or the location of various periodic cells. With reference to Table 1, the stable P- 1 point XT of equation (19) at the origin is replaced by five P-l cells represented by the symbol z; the by the symbol stable P- 1 points x& = (f r, 0) are replaced by five P- 1 cells represented +; the unstable P-l points x: and XT are each replaced by four P-l cells represented by T and s, respectively. The stable advancing P- 1 solutions x2 and XT are each replaced by one P-l cell and a group of P-4 cells represented by / and \, respectively. It should be noted that the periodic cells corresponding to the unstable advancing P-l solutions xg and xz are missing from this run of the cell mapping algorithm. Figure 12(b) shows the IO-step domains of attraction for the stable solutions XT, x&, respectively by symbols z, +, /, and \. The cells with symbols s x6* and XT, represented and T are cells which straddle over the stable manifolds (i.e., approaching separatrices) of the saddle points xt and XT. Figure 12(c) shows the 25-step domains of attraction

GLOBAL 3.5

ANALYSIS

OF A NON-LINEAR

419

SYSTEM

(a)

3.0 r----2.51 / 2.01

, ,:I

I 1.51 1.01

1 L

03

t *”

0.0

.. .

.“.

1

z

1

*

,,I.

. .

-0.5

1

-10

-1.51 -2.01 -2.5 -3.0

-3.51, -1.0

3.5

-0.75

-0.5

-0.25

0.0 x,/n

(b)

0.25

0.5

0.75

1

1.0

---I I

3.0 t 2.5 111 111111,111, ,,,,,,,,,, ,,;$“‘;;”

-1.0

-1.5 1 -2.0 1 -2.5 t

I::[,, , , / , , / , -1.0

-0.75

-0.5

-0.25

0.0

0.25

0.5

0.75

l-

Figure 12. (a) O-step, (b) IO-step, (c) 2S-step and (d) total (51~step) domains of attraction for the periodic cells of the mapping (46) for the case p = 0.02~, (2 = I .O, with TV,,= 100, NC, = 101, TV,= IO 100.

420

R.

2-5 2.0

S. GUTTALU

AND

C.

S.

HSU

I,, ii,,, “’ “’ I,, I 11, I, ,,,,,,,,,, “’ ,,,,,,, IN,,,,,,,,,,, ’ , y:,y ,,,:::::: ,,,,, :::;::::” ,,,,,, ““’ ’ , $,,, I, ,111, ,111, I. ,,,,,,,,,,,,,,,,, ,l, m,::,,,,:::: z::::::::::::, ,; ,, ,:::::::;:::‘;:::;:;::;::“,,,:, s 11, ,, ,111 lljjlllj” ‘::::::::“:::”

‘/ I

I 1

z “m

0.5 t *”

0.0 -0.5

r

.....I...........:;;........,.......... ..,II.... .....2 ..,.... ...I~llll~ll.l,ll~,ll,l~~,~~, .1, ...........I,,.......I ..........11,,...-~1,z ..,..............I.... at,. ..,I..... ,,,z,,, I z*‘.‘%::::.x “‘: T:‘:..? “:‘;:“‘::‘:i:‘:“‘“:;;:‘:::; ,.,....... ;.'..:::;$:::.

2,:

xx:*

-1.0

.:.;:;.:,...,:....::, . . .. . . I ,:

:yl'

. . .. . .

1

:,,,:z-::::.::: ;.

* :-

I

:.

.,1,1.....,.,

. . . . . . . ..I... II ::::':"...::::.""':::"'

;::"+$m_

22

. . . .. .. .

:i:"'._,..

.

~L~z..ir~:"" i:: ..:: ::f::r "Si:":.. .I, 1121.1*... .zz .., Xl . .. .. .. . .;,z_r.. ll.,l. .. . I> 7 . .. . I 1.11.1.

. . . ..I...

z ?

'

.

. . *:::, -zz ;y,::......,

"'1'*'=

-1-5 -2.0 -2.5

-1.0

-0-75

-0.5

-025

0.0

0.25

05

0.25

0 5

o-75

10

x,/n 3.5 3.0 2.5

-2.0 -2.5 -3.0 -3.5

t. -1.0

-0.75

-0.5

-0.25

0.0 X(/lr

Figure 12 (cont.).

0.75

1.c

GLOBAL

ANALYSIS

OF

A NON-LINEAR

421

SYSTEM

where the complicated interweaving pattern of the domains of attraction is beginning to emerge. Figure 12(d) is the total domains of attraction where the largest number of steps any cell is away from an attracting periodic cell is 5 1. Out of 10 100 regular cells, 1587 cells belong to the domain of attraction of xi* , 5779 to x:, 1296 to xz and 1296 to XT. To each of the unstable core of P-l cells (denoted by the symbols s and T) belong 71 cells on account of being on the stable manifolds of the saddle points x2 and XT. Here none of the regular cells is mapped to the sink cell. Figure 12(d) not only delineates very nicely the distribution of the four domains of attraction but also shows clearly the symmetry of the domains with respect to the origin. Figure 12(d) should be contrasted with Figures 5(a) and (b) obtained by point mapping methods. We now consider the case for which p = 0*0027r, CY= 3.75, when P- 1 and P-3 solutions of equation (19) exist. The cell mapping (46) will generate many pseudo-periodic cells for this case since p is very small. Whether or not these periodic cells are true periodic cells or pesudo-periodic cells may be determined by simply applying the cell mapping (46) with smaller cell sizes to a smaller region containing them. Of course, one needs to use the appropriate GK mapping, K 2 1, to generate the associated cell mapping. Another way is to iterate the center point of the periodic cells under the mapping GK to determine which of the P-K solutions attracts it. If all of the periodic cells in the core are attracted by the same P-K point, then one knows that the periodic cell in the core closest to the P-K point is the true periodic cell while the rest of them are pseudo-periodic cells. This method of testing, however, does not reveal the true nature of the periodic cells lying on the separatrices of the point mapping. Some of these considerations were taken into account in obtaining the total domains of attraction in Figures 13(a)-(c) in the phase o-5

0.4

0.3

0.2

0.1

5

0.0

-0.1

-0.2

-0-3

-0.5 -0 3

-0.2

-0.1

o-o x,/r

01

0.2

0.:

Figure 13. Total domains of attraction (a) 47-step (see Figure 7), (b) 118-step (see Figure 8) and (c) 84-step (see Figure 9) for the periodic cells of the mapping (46) for the case @ = O.O027r, a = 3.75, with NC, = N, = 101, N,=lO201.

422

R. S. GUT-I-ALU

AND

C. S. HSU

1.00

l-02

0.6 t(b)

-0.6

, 0.92

0.94

0.96

0.96

Figure

13 (cont.).

104

1.06

l-08

GLOBAL

ANALYSIS

OF A NON-LINEAR

SYSTEM

423

plane, which should be compared with Figures 7-9, respectively, produced by the point mapping method. In all these figures, the region of the phase plane indicated in the graphs is divided into a total of N, = 10 201 regular cells with N,, = N,, = 101. The blank spaces mean that the regular cells in these locations are mapped into the sink cell and hence their evolution is not determined by this run of the cell mapping algorithm. For Figure 13(a), the range is -0.37~~~~ %0*37r, -O.~TSX,SO.~~~, and h, = 0.005 9407~, h, = 0.009 901 rr. We have used the G3 mapping to generate the associated cell mapping which allows us to obtain the domains of attraction of the P-3 solutions separately. Hence, the number of steps in Figures 13(a)-(c) should be interpreted appropriately. The domains of attraction of the P-l solution are denoted by the symbol +, and the domains of attraction of the P-3 solutions by z, v, /, \, x and Y, while the P-3 cells lying on the stable manifolds of the saddle points are indicated by A, B, c, D, E and F. The largest number of steps is 47. Here one notices that the domains of attraction of the P-3 points denoted by / and \ each contain only two cells as determined by this run of the algorithm; a larger number of regular cells is needed to obtain more cells there. For Figure 13(b), the range is 0.92~~ x1 s 1.08~, -0.67~ s X?S 0*67~, and h, = 0.001 5847, hZ= 0.11 8817~ The symbols here have the same meaning as in Figure 13(a). Again, we have used here G3 mapping; the largest number of steps is 118. This figure excellently reproduces the domains of attraction of the P-l and P-3 solutions of Figure 8. Figure 13(c) is for -0*6?r~x,~O*6~, 0.57r4?cz<3.57r, with h,=O=Oll 88171; h,== 0.029 7037r. We have used the G5 mapping to generate the cell mapping for the reason that there exists a pair of advancing P-5 solutions: in the phase plane one is given by ~x~(~),~=1,2.3,4,5)={(+0~286018~,+1~665467x),~+I~5425347~,+1~778213~~ (+4.450

663~, 12.855

3177~), (+5.650 629~, +0.726933x),

(+8.018

5107r, +2.942 686n)}

and the other one is given by {x;(j) = -xt(j),j = 1,2,3,4,5}. These advancing P-5 solutions are such that G’(xT) = XT + 10~ and G’(xT) = XT - 10~. The domain of attraction of the advancing P-l solution in Figure 13(c) is indicated by the symbol + and the domains of attraction of the advancing P-5 solutions are indicated by A, 6, c, D and E, and we have taken x, modulo 27r. A P-5 cell in the region A is mapped in sequence to the regions 6, c, D and E before returning to A. There are only a few cells in the domains of attraction of P-5 points found by this computer run. A larger number of regular cells are required if one wants to see better the structure of the domains. It is interesting to note in Figure 9 that there are five regions around which the domain of attraction of the advancing P-l point winds, and they contain a stable advancing P-5 solution which is easily discovered by the cell mapping algorithm. Figure 14 is included to show how sensitive the evolution of the dynamical system (19) is to two nearby starting states given by x, = (-0.450~, +1.7327~) and x2 = (-0.4557r, +1.7327~). The initial states are indicated by the symbol + in Figure 14. These two starting points are very close and yet their evolution leads to an advancing P-l solution for x, (as in Figure 14(a) and to the advancing P-5 solution for x2 (as in Figure 14(b)). Moreover, Figure 14 shows that in the early stage the point trajectories of the two evolutions are close, but eventually they converge to different solutions. In Figure 14(a) symbols 1, 2, 3, . indicate the evolution from x,. In order to appreciate the existence of this P-5 solution, we have included the continuous trajectories of this solution in Figure 15 in the phase plane. In Figure 15. starting at state A, the system occupies the states B, C, D and E at the end of each period before turning to A at the end of the fifth period.

424

R. S. GUT-t-ALU AND C. S. HSU 3.57--(b)

(a)

3.0

I

E,,..

4

2

IC

2,5" ..,: ..... ... <' _.:. _.. .::. ,.,:;$Y ....

t *" 2.0l*

l3i

5

:

1.5-

n

l.O3

b.

0

1.6, -0.6

-0.4

I

I

I

i

0.0

0.2

0.4

0.6

I

I

-0.2

Figure 14. Evolution of the dynamical system (19) for the case p = 0.003q a = 3.75, starting from the initial state (a) x, = (-0.45Orr, +1.7327~), (b) x2= (-0.45&r, +1.7327a). x, converges to advancing P-l and x2 to advancing-P-5.

3 o-

-72

2 5-

k > 2 oA

A

I 5-

J, 05 00

2.0

40

60 X,/T

60

I

100

I

12 0

Figure 15. Continuous trajectories of the advancing P-5 solution in the phase plane for the case p = 0&)2~, (I = 3,75.

We now consider the case for which the system (19) has P-l and P-2 solutions (see Table 3). Figure 16(a)-(c) shows the domains of attraction obtained by using the cell mapping algorithm for the case p = O*lr, (Y= 4.5, with IV,, = 100, NC, = 101, NC = 101 00. Here, we have used the G* mapping to determine the associated cell mapping. The O-step

GLOBAL

-0.5

ANALYSIS

OF A NON-LINEAR

425

SYSTEM

c

-1-o i

-1.0

-0.75

-0.5

-0.25

0.0 x,/TT

O-25

0.5

0.75

1,’

-0.75

-0.5

-0.25

0 0 X,/H

O-25

05

o-75

1.1 0

3.5 3-o 2.5

-2.0 -2.5 -3.0 -3.5

L

.O

1

Figure 16. (a) O-step, (b) five-step and (c) total (ISstep) domains of attraction for the periodic mapping (46) for the case p =O.l T, (1=4.5, with NC, = 100, NC, = 101, IV, = 10 100.

cells of the

426

25

-2-5

domains of attraction are shown in Figure 16(a). The unstable P-l point XT at the origin is replaced by one P- 1 cell and one group of P-2 cells represented by s; unstable P- 1 points x& by one isolated P-l cell T; advancing P-l points x$ and XT by four P-l cells / and \, respectively; each of the two P-2 points xTO(l) and x:(,(2) by a core of five periodic cells, with the cells in one core designated by + and the other by X; each of the two P-2 points xT,( 1) and x?,(2) by a core of two periodic cells designated by o and z, P-l points respectively. The unstable P-l points x ,* and xc and the unstable advancing xt and x$ do not show up as periodic cells in this run. Figure 16(b) shows the five-step domains of attraction which should be compared with Figure 11 obtained by the point mapping technique. Figure 16(c) is the total domains of attraction for the six groups of attracting cells and corresponds to 15-step domains of attraction. Of the 10 100 regular cells, 2538 cells belong to x,*, 2538 to xf, 2 x 1249 to x&, and 2 x 1204 to XT,. On account of being on the stable manifolds of the saddle points of the point mapping, a total of 118 cells represented by the symbols s and T have appeared separately. Here one can of course also combine the two domains of attraction of a single P-2 solution into one domain of attraction for the stable P-2 motion.

7. CONCLUDING

REMARKS

In this paper we have presented an analysis of a strongly non-linear parametrically excited system. In particular, we have used this analysis to assess the virtues of various methods available to study the global behavior of non-linear systems. The problem studied is an instructive one, because, on the one hand, the formulation in terms of a point

GLOBAL

ANALYSIS OF A NON-LINEAR

SYSTEM

427

mapping is exacf, and yet, on the other hand, it demonstrates very well the extreme complexity of global behavior a seemingly simple mechanical system can have. A global analysis of a strongly non-linear system is always a very difficult one. In this paper we have used three different techniques to determine the global domains of attraction: the method of separatrices, the method of expanding domains, and the cell mapping algorithm. It turns out that so far as the computation1 effort is concerned, the cell mapping technique is several orders of magnitude more efficient than the other two. Thus, the cell mapping seems to offer a very attractive approach to the global analysis of non-linear systems. ACKNOWLEDGMENT

The material is based upon work supported by the National Science Foundation

under

grant no. MEA-8217471.

REFERENCES 1. C. S. Hsu 1977 Adoances in Applied Mechanics 17,245-301.

On non-linear parametric excitation

problems. 2. R. S. GU~TALU 3.

4. 5. 6. 7. 8. 9.

1981 Ph.D. dissertation, University of California, Berkeley. On point mapping methods for studying nonlinear dynamical systems. C. S. Hsu 1980 Journal of Applied Mechanics 47, 931-939. A theory of cell-to-cell mapping dynamical systems. C. S. Hsu and R. S. GUTTALU 1980 Journal of Applied Mechanics 47, 940-948. An unravelling algorithm for global analysis of dynamical systems: an application of cell-to-cell mappings. C. S Hsu, H. C. YEE and W. H. CHENG 1977 Journal of Sound and Vibration 50, 95-116. Steady-state response of a non-linear system under impulsive periodic parametric excitation. J. GUCKENHE~MER and P. HOLMES 1983 Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. Berlin: Springer-Verlag. V. I. ARNOLD 1983 Geometrical Methods in the Theory of Ordinary Diferential Equations. Berlin: Springer-Verlag. G. IOOSS 1979 Bifurcation of Maps and Applications. Amsterdam: North Holland. J. BERNOUSSOU 1977 Point Mapping Stability. Oxford; Pergamon Press,