A global study of nonlinear dynamical systems by a combined numerical-analytical approach

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Request accessible transcript

Transcript (if available)

Content A GLOBAL STUDY OF NONLINEAR DYNAMICAL SYSTEMS BY A COMBINED NUMERICAL-ANALYTICAL APPROACH by Michael C. Golat A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (MECHANICAL ENGINEERING) December 2002 Copyright 2002 Michael C. Golat R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. UMI Number: 3093955 Golat, Michael Charles All rights reserved. ® UMI UMI Microform 3093955 Copyright 2003 by ProQuest Information and Learning Company. All rights reserved. This microform edition is protected against unauthorized copying under Title 17, United States Code. ProQuest Information and Learning Company 300 North Zeeb Road P.O. Box 1346 Ann Arbor, M l 48106-1346 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. UNIVERSITY OF SOUTHERN CALIFORNIA THE GRADUATE SCHOOL UNIVERSITY PARK LOS ANGELES, CALIFORNIA 90089-1695 This dissertation, written by under the direction o f h dissertation committee, and approved by all its members, has been presented to and accepted by the Director o f Graduate and Professional Programs, in partial fulfillment o f the requirements fo r the degree of DOCTOR OF PHILOSOPHY Director Date '1&2X- Dissertation Committee Chair R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. D ed ication To Victoria, Trevor, and Taylor. R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. A cknow ledgem ents I wish, to express my deep gratitude to my research advisor, Professor Henryk Flashner. His patience, interest, and unwavering support were greatly appreciated and have made this work possible. I also thank Professor Marijan Dravinski and Professor Gary Rosen for serving on my thesis committee. Above all I would like to thank my wife Victoria for her constant encouragement and support during this extended period of time. R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. C ontents D edication ii A cknowledgem ents iii List Of Tables vi List Of Figures vii Abstract ix 1 Introduction 1 2 Background 6 2.1 Point Mapping A n a ly s is ........................................................................... 6 2.1.1 Periodic Solutions and Local Stability Properties ................. 8 2.1.2 Bifurcation C onditions.................................................................. 9 2.2 Truncated Point M appings....................................................................... 10 2.2.1 Jets and Truncation Operation .................................................. 11 2.2.2 Computation of Point M appings.................................................. 13 2.2.3 Dependence on Param eters............................................................ 15 2.3 Cell Mapping A n a ly sis.............................................................................. 16 2.3.1 Periodic M o tio n s ............................................................................ 18 2.3.2 Creating a Cell M ap p in g............................................................... 20 2.3.3 Global A nalysis............................................................................... 21 3 Expanded Point M apping 25 3.1 Introduction................................................................................................. 25 3.2 EPM Form ulation....................................................................................... 28 3.2.1 State Space D iscretization............................................................ 28 3.2.2 Cell Space A nalysis......................................................................... 28 3.2.3 Analytical Study of Trajectory E volution.................................. 29 3.2.4 Local Stability ............................................................................... 33 iv R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. 3.2.5 Bifurcation Analysis ...................................................................... 34 3.2.5.1 Bifurcation of P-K solutions into P-K solutions . . 36 3.2.5.2 Bifurcation of P-K solutions into P-2K solutions . 38 3.2.5.3 Bifurcation of P-K solutions into P-MK solutions . 39 3.3 EPM P rocedure........................................................................................... 40 3.4 EPM S u m m a ry ........................................................................................... 42 4 Pendulum w ith Vertical Support M otion 44 4.1 System Dynamic Equations .................................................................... 44 4.2 Bifurcation Prelim inaries........................................................................... 45 4.3 P-l to P-2 B ifu rc a tio n ............................................................................... 49 4.4 P -l to P-3 B ifu rc a tio n ............................................................................... 56 4.5 EPM Analytical V erification.................................................................... 60 4.6 Global A n aly sis............................................................................................ 67 4.6.1 Basins of Attraction .................................................................... 74 4.6.2 Analytical V erification................................................................. 76 5 Pendulum w ith Periodically E xcited Support in the Plane 78 5.1 System Dynamic Equations ..................................................................... 79 5.2 Global A n aly sis............................................................................................ 83 5.2.1 Periodic Solutions and Basins of Attraction .......................... 84 5.2.2 Stability A nalysis........................................................................... 87 5.2.3 Trajectory Evolution Accuracy ................................................. 96 5.2.4 Computational S tu d ie s................................................................. 99 5.3 Bifurcation Analysis ....................................................................................100 5.3.1 Mapping P relim inaries....................................................................101 5.3.2 Analytical Mapping F o rm u lae.......................................................104 5.3.3 Bifurcation Analytical R e su lts.......................................................106 5.3.4 Bifurcation P-K Solutions and Basins of A ttra c tio n ................108 6 Concluding Rem arks 121 Reference List 123 A ppendix A Expanded Point Mapping Source C o d e .............................................................128 A.l Declarations and In itializ a tio n s................................................................ 128 A.2 Preliminary EPM Com putations................................................................ 129 A.3 Trajectory Com putations............................................................................. 130 A ppendix B Floquet T h e o ry ....................................................................................................... 132 v R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. List O f Tables 5.1 EPM P-K points vs. Direct Integration P-l point............................... 8 6 5.2 CPU time requirements for P-K solutions and basins of attraction. 100 5.3 Bifurcation EPM P-K Solutions and Exact P-K Solutions.................. 114 vi R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. List O f Figures 2.1 Two-dimensional cell state space.............................................................. 19 2.2 Example of a simple cell mapping with periodic (shaded) and tran sient cells. z 2 is a P-K solution of period 1 . z3, z4, and z 5 constitute a P-K solution of period 3.......................................................................... 22 3.1 Simple cell mapping spatial discretization error.................................... 27 3.2 Two EPM mapping steps........................................................................... 31 3.3 EPM mapping for cell Zi............................................................................ 43 4.1 Pendulum with a vertically moving support........................................... 46 4.2 Bifurcation to P-2 solution, £ = 0 and a\ = 0.1.................................... 51 4.3 A stable P-2 trajectory associated with the solution (A i,A 2) for ( = 0, a\ = 0.1, and t* 2 = 2.00102............................................ 53 4.4 A stable P -2 trajectory associated with the solution (B1,B 2) for f = 0, a\ = 0.1, and a 2 = 2.12602............................................................. 54 4.5 A stable P-2 trajectory associated with the solution (6 j , C2) for £ = 0, af — 0.1, and a 2 — 2.12602............................................................. 55 4.6 Bifurcation to stable P-3 solutions B and P , £ = 0 and a\ — 0.1. . . 57 4.7 Bifurcation to unstable P-3 solutions A and D, £ = 0 and a\ = 0.1. 58 4.8 The stable P-3 trajectory associated with the solution ( P i,P 2 ,P 3) for £ = 0, a\ = 0.1, and a 2 = 1.62158....................................................... 61 4.9 The stable P-3 trajectory associated with the solution (Pi, C2, C ?,) for f = 0, al = 0.1, and a 2 = 1.62158....................................................... 62 4.10 The unstable P-3 trajectory associated with the solution (Ai, A2, A3) for ( = 0, a\ = 0.1, and a 2 = 1.62158....................................................... 63 4.11 The unstable P-3 trajectory associated with the solution (Th, D2, D3) for ^ = 0, af = 0.1, and a 2 = 1.62158....................................................... 64 4.12 The stable P - 6 trajectory associated with the solution (Ri ,..., Re) for £ = 0, af = 0.1, and a 2 = 1.64008...................................................... 65 4.13 The stable P - 6 trajectory associated with the solution (Si, ..., Se) for £ = 0, oi\ = 0.1, and a 2 = 1.64008....................................................... 6 6 4.14 P-3 solutions A, B , P , and I) for £ = 0.01 and eh ( = 0.1..........................69 4.15 The stable P-3 trajectory associated with the solution (Pi, P 2, P 3) for £ = 0.01, a\ = 0.1, and a2 = 1.59158................................................. 70 vii R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. 4.16 The stable P-3 trajectory associated with the solution (6 j , C 2 , C3 ) for £ = 0.01, a\ — 0.1, and a 2 — 1.59158................................................. 71 4.17 The unstable P-3 trajectory associated with the solution (Ai, A2, A3 ) for £ = 0.01, aj = 0.1, and a 2 = 1.59158................................................. 72 4.18 The unstable P-3 trajectory associated with the solution (P i, P 2, D 3 ) for £ = 0.01, c? x = 0.1, and cr2 = 1-59158................................................. 73 4.19 P-3 EPM Basins of Attraction.................................................................. 75 5.1 Parametrically forced pendulum with viscous damping....................... 81 5.2 EPM P-K Points.......................................................................................... 8 8 5.3 EPM Basins of Attraction.......................................................................... 89 5.4 Direct Integration Periodic Solutions....................................................... 90 5.5 Direct Integration Basins of A ttraction................................................... 91 5.6 SCM P-K Cells............................................................................................. 92 5.7 SCM Basins of Attraction........................................................................... 93 5.8 Trajectory Comparison. SCM versus Exact Integration...................... 97 5.9 Trajectory Comparison. EPM versus Exact Integration...................... 98 5.10 Bifurcation-Frequency Response Solutions................................................ 108 5.11 A stable P-l trajectory associated with the solution Pi for // = 0.1, X q = 0.5, and r = 0.88................................................................................... 109 5.12 A stable P-2 trajectory associated with the solution (Ai ,A 2) for fj, = 0.1, Xo = 0.5, and r = 1.00....................................................................110 5.13 A stable PA trajectory associated with the solution (P i, P 2, P 3 , P 4 ) for p = 0.1, X q — 0.5, and r = 1.075........................................................... I l l 5.14 A stable P - 8 trajectory associated with the solution (Cj,..., C$) for H = 0.1, X Q = 0.5, and r = 1.080............................................................... 112 5.15 A stable P-16 trajectory associated with the solution (P i,..., D ^) for fi = 0.1, Xo = 0.5, and r = 1.083......................................................... 113 5.16 P-l EPM Basins of Attraction.................................................................... 116 5.17 P-2 EPM Basins of Attraction.................................................................... 117 5.18 P-4 EPM Basins of Attraction.................................................................... 118 5.19 P - 8 EPM Basins of Attraction.................................................................... 119 5.20 P-16 EPM Basins of Attraction................................................................... 120 A.l Flow chart for the EPM program................................................................ 131 viii R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. A bstract A new numerical-analytical method for the combined global-local analysis of non linear periodic systems referred to as an Expanded Point Mapping (EPM) is pre sented. This methodology combines the cell to cell mapping and point mapping methods to investigate the basins of attraction and stability characteristics of equi librium points and periodic solutions of nonlinear periodic systems. The proposed method is applicable to multi-degree of freedom systems, multi-parameter systems, and allows analytical studies of local stability characteristics of steady state solu tions. In addition, the EPM approach allows the study of stability characteristics as a function of system parameters to obtain analytical conditions for bifurcation. In this dissertation, the theoretical basis for the EPM method is formulated and a procedure for the analysis of nonlinear dynamical systems is presented. Two ana lytical studies are performed. The first analysis consists of a study of a pendulum with a vertical moving support. Next a more general investigation of a pendulum with a periodically excited support in the plane is used to illustrate the method. The results demonstrate the efficiency and accuracy of the proposed approach in analyzing nonlinear periodic systems. ix R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. C hapter 1 Introduction A wide class of dynamical systems are governed by periodic differential equations. This class includes systems with periodically varying parameters as well as systems with periodic forcing. These problems are of great mathematical challenge and have attracted a great deal of theoretical research (see Cesari [4], Hale [20], Urabe [52], and Yakubovich and Starzhinskii [55]). Likewise, they have great practical importance because of their many applications in various fields of science and engineering (see Bolotin [3], Kauderer [30], Lichtenberg and Lieberman [31], and Hayashi [22]). The majority of analytical methods for this class of problems can be divided into two categories: perturbation analysis using asymptotic expansion in terms of small parameters and numerical simulation (see Bogoliubov and Mitropolsky [2], Nayfeh and Mook [36], Nayfeh [35], Moon [34], and Thompson and Stewart [46]). Perturbation methods allow for an analytical treatment of the problem but, 1 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. being local in nature, are valid only for small parameter variations and in a small neighborhood of an equilibrium or periodic solution. For instance, the existence of large damping for various systems usually restricts the scope of application of perturbation methods. Moreover, perturbation methods are in general difficult to implement for systems of order greater than two. In addition, global domains of attraction can only be qualitatively extrapolated from the localized behavior around fixed points and periodic solutions. Numerical simulation studies are well suited for investigating a system’s behav ior in the entire state space when one considers the extremely high computational capabilities of present day computers. It is straightforward to perform digital sim ulations on a nonlinear problem by implementing direct integration methods in order to gain insight into the behavior of a dynamical system. However, direct numerical computations are limited in their overall effectiveness. Even if a digital simulation represents exactly the dynamics of the nonlinear system being studied, it represents only a finite number of solutions for a given set of initial conditions and fixed parameter values. Consequently, such important topics as stability analysis and global domains of attraction cannot be addressed in a comprehensive theoret ical fashion. Furthermore, repeated simulation runs over a broad range of initial conditions and parameters require prohibitively large computational resources, es pecially in the case of parametric studies such as bifurcation analysis. 2 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. It is from these aforementioned deficiencies of direct numerical techniques that a number of combined numerical-analytical methods have been developed. An at tractive way to analyze periodic systems is to formulate their discrete-time dynam ics by defining a point mapping (see Bernussou [1]). By doing this, the dynamics of the original system is given in terms of difference (algebraic) equations, rather than in terms of time varying differential equations. Point mapping techniques are widely used today not only to study theoretical aspects of discrete-time systems but also to provide a computational basis for understanding their global dynamics (see Guckenheimer and Holmes [13] and Poincare [39]). Similarly, the method of cell to cell mapping has been extensively employed for studying the global behavior of nonlinear systems. A significant amount of research involving cell to cell map ping techniques has been performed by Hsu and his coworkers [25, 27, 28], Guttalu, et al [16]. Another cell mapping method called interpolated cell mapping has been developed by Tongue [47]. However, its main line of development differs in spirit from the cell to cell mapping method (see Hsu [26])). For recent modifications to the cell mapping method the reader is referred to Dellnitz, Guder, Kreuzer, van der Spek, van Campen, et al [14, 15, 53]. In this thesis, a new numerical-analytical method referred to as an Expanded Point Mapping (EPM) is introduced. The Expanded Point Mapping approach combines key elements of the cell to cell mapping and point mapping methods to form a comprehensive methodology to study the global behavior of nonlinear 3 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. periodic systems in conjunction with analytical investigations of stability char acteristics of equilibrium points and periodic solutions. The proposed method is applicable to multi-degree of freedom systems and the analytical procedure is formulated for any order of approximation. Consequently, the approach al lows convenient analysis of multi-dimensional and multi-parameter systems. The global analysis performed by the EPM formulation consists of obtaining all possi ble equilibrium and periodic solutions and determining their domains of attraction. Moreover, the EPM technique generates an analytical expression representing each periodic solution allowing local stability analysis of the equilibrium and periodic solutions computed in the global analysis part of the approach. In addition, the EPM approach provides the tools to express the solutions as a function of system parameters which allows the study of stability characteristics as function of these parameters and subsequently obtaining bifurcation conditions. The thesis is organized as follows. In Chapter 2 the point mapping method and an algorithm to obtain its approximation are summarized along with a review of the cell to cell mapping method. Fundamental definitions, concepts, etc. will be described in detail in order to provide the requisite theoretical background for the development of the EPM methodology. In Chapter 3 the Expanded Point Mapping method is formulated and a procedure for its application to analysis of nonlinear periodic systems is described. The proposed approach is demonstrated first on a pendulum with a vertically moving support in Chapter 4. Then in Chapter 5 the 4 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. EPM method is used to investigate the global behavior for the generalized case of an inverted pendulum with a support periodically excited in the plane. Finally, concluding remarks and suggestions for future research are provided in Chapter 6 . 5 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. C hapter 2 Background 2.1 P oin t M apping A nalysis Consider a dynamical system described by a set of N ordinary differential equations x(i) = f(*,x(i),s) (2.1) where t ( E R + denotes time, x(i) € R ^ is the state vector, s £ is the parameter vector, and f : R + x R ^ x RL — y R N is analytic for every x for a given value of the parameter vector s. Moreover, f is assumed to be periodic in t with period T, i.e. f(t,x , s) = f (t + T, x, s). To analyze the equations in (2.1), we observe the state of the system at discrete instances every T units of time. A point mapping is the dynamic relationship between the state of the system at t = nT and the state at t = (n + 1 )T which can be expressed as a set of difference (algebraic) equations 6 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. x(n + 1) = G (x(n), s) n = 0,1,2,... (2.2) where x(n) and x(n + 1) are the states of the system at t = nT and t = (n + 1)T, respectively. Note that since f is analytic, the differential equations (2.1) satisfy Lipschitz conditions and therefore G is one-to-one for a fixed value of s. The point mapping formulation (2.2) allows for convenient determination of periodic solutions and their stability characteristics. A periodic solution of period K T , or in short a P-K solution of the point mapping (2.2) for some s = s*, consists of K distinct points x*(_f), j = 1,2,..., K such that x*(m + l) = G m(x*(l), s*) m = 1,2,..., if — 1 (2.3) x*(l) = G A(x*(l), s *) where G m is the point mapping G applied m times. From this definition, an equilibrium state of (2.2) can be viewed as a P-l solution. Also, an equilibrium point x* of (2 . 1 ) which satisfies f(i,x*,s*) = 0 Vt is clearly a solution of the corresponding point mapping (2.2). The terms fixed point, critical point, and stationary point are also used in various literature to describe x*. R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. D efinition 2.1 The domain or basin of attraction J5X » of an asymptotically stable P-K point x* is defined as Bx. = {x € i N | G jK(x) -> x*, j -> 0 0 } (2.4) 2.1.1 Periodic Solutions and Local Stability Properties Finding a P-K solution of the point mapping G is equivalent to finding a periodic solution of the original continuous time system. Having found a P-K solution by solving the algebraic system of equations (2 .2 ), its local stability properties (in the sense of Liapunov) are determined by the eigenvalues of the m atrix H given by K BG H = J J H*(j) where H*(j) = — i = i (2.5) x = x * (i),s = s * The eigenvalues of H are denoted by A;(H) i = 1,2,..., N. Using the theory of linear difference equations (see Miller [33]), the local stability of a P-K solution can be summarized as follows. 1. A P-K solution of (2.2) is locally asymptotically stable if and only if |A t -(H)| < 1 = 1,2,...,A T . 2. A P-K solution of (2.2) is unstable if |A*(H) | > 1 for some i = 1,2,..., N. R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. 3. If |A t -(H)| = 1 for some i = j and |A t -(H)| < 1 Vi = 1,2,..., N i ^ j, then the linearization process is inconclusive. Moreover, if this condition occurs while s is varied, bifurcation conditions are established. 2.1.2 Bifurcation Conditions The bifurcation conditions for a given solution are stated in the following. If, for some s G S, there exists an integer M such that then a bifurcation of a P-K solution to a P-MK solution may occur. Alternately the bifurcation conditions (see Hsu [24]) may be expressed as det(I - H M(s)) = 0 (2.6) 1. Bifurcation of P-K — > P-K : det(I - H(s)) = 0 (2.7) 2. Bifurcation of P-K — > ■ P-2 A: det(I + H(s)) = 0 (2.8) 9 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. 3. Bifurcation of P-K — > P-MK: det(I — H M/,2(s)) = 0, M even, M /2 even; det(I - H(s) + H 2( s ) ------- + = 0, M even, M/2 odd- det(I + H 2(s) + H 3(s) + • • • + H m_1(s)) = 0, M odd. (2.9) It should be noted that for two-dimensional Hamiltonian systems, equation (2.6) is satisfied for all s. Therefore, a P-K solution cannot lose stability via bifurcations to P-MK solutions with M > 2. The mechanism of stability loss for Hamiltonian systems is only via bifurcations to P-K solutions or to P-2K solutions. 2.2 Truncated P oint M appings The algorithm for obtaining an approximation, in analytic form, for the point mapping G (x(n),s) of (2.2) was first introduced by Flashner [8] and Flashner and Hsu [10]. For a modification of this procedure see Flashner and Guttalu [9]. The underlying philosophy and the main steps of the algorithm are presented in the following. For more details about this procedure see Guttalu and Flashner [17, 19]. 10 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. 2.2.1 Jets and Truncation Operation Since the function f(px(t),s) in equation (2.1) is analytic, it can be expressed by a series of the form f(i,x,s) = p (t,x ,s) (2.10) p r where p(t,x, s) denotes a vector homogeneous multinomial of degree P in the state variables x,{ i — 1,2,..., N . The symbol denotes summation over all sequences j i , .. ■ ,jr such that XwtLi jk — r an<l r — 1 , 2 ,.. ., P, i.e. a vector homogeneous multinomial of degree P contains terms where the sum of the state variable exponents is less than or equal to P. To abbreviate notation denote the set of indices ji , . . . , j r by Er. Then we denote the coefficients of the vector homogeneous multinomial in (2.10) by the symbol b ^ ( t,s ) . These coefficients, in general, can be functions of the parameter vector s. The algorithm for determining point mappings is based on multinomial func tions by truncating at some degree k. Let p (i,x ,s ) be a multinomial of degree greater than k. Define the truncation operation at degree k by k p {k)(t, X , s) = p(t, X , s) k = bErf (^ S) X1 ■ ■ ■ XN (2-11) r li R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Note that in our definition of the multinomial p(t,x,s), we assume no terms of order zero in the components of x. More specifically, there is no hard forcing in the system (2.1). Let p(t, x, s) and q (t, x, s) be two vector homogeneous multinomials, then the truncation operation obeys the following rules for addition, multiplication, and composition (see Poston and Stewart [40]): For the case of P — > oo in equation (2.10), f(f, x, s) is approximated by a k-jet p + q* = p k + q k (2.12) p • q * = p k ■ q k (2.13) p(*,q,s)* = p k(t, q ki s) (2.14) denoted by j kf , i.e. by a Taylor series expansion. It can be shown (see Poston and Stewart [40]) that the k-jets obey the rules in equations (2.12-2.14). Namely j*(p + q) = i ‘p + i ‘q (2.15) / ( p - q ) = j ‘p j ‘q (2.16) (2.17) 12 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. 2.2.2 Com putation of Point Mappings The computational algorithm for evaluating a point mapping is based on the fact that the operations of multinomial addition, multiplication, and telescoping (com position) can be interchanged with the truncation operation. Moreover, the Runge- K utta method of integration can be expressed as a sequence of multinomial tele scoping and truncation operations which finally results in a truncated multinomial expression for the point mapping G (x(n),s) of (2.2). Using the Runge-Kutta method, the state of the system at time t = tp + h can be formulated as follows (see Henrici [23]): M x(tp + h) = x(tp) + h dmk m(t, x, s) (2.18) m— 1 where M is the order of the Runge-Kutta method, h is a fixed time step, and dm are certain constants determined by the Runge-Kutta method. The vectors km are given by kT O (f,x,s) = i(tp + ham,x(tp) + /icT O km -^s) m = l,2, . . . , M (2.19) where the constants am and cm are defined by the particular Runge-Kutta method implemented, with a\ = 0 and c\ — 0 . 13 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. To derive the algorithm, the period T is divided into Nt subintervals of length h = T /N t. We are interested in computing a A r-jet of the point mapping G in (2.2). The point mapping algorithm takes advantage of the structure of the dynamic equation for the system in (2.18). It can be shown that each function km can be expressed by telescoping (composition) and addition of truncated multinomials, resulting in j fckm(i; x, s) = P(£)(«; x, s) m = l , 2 , ...,M (2 .2 0 ) Similarly, equation (2.18) consists of additions of truncated multinomials yielding M j kx(tp + ih) = x + h dmp^k)(i; x, s) = Q(fc)(i; x, s) (2 .2 1 ) 771 — 1 Repeating the procedure Nt — 1 times yields x((n+ 1)T) = Q ik)(Nt- , (Q(fc)(W - 1; (• • • (Q(fc)(l; Q(fc)(0 ;x, s),s)), s) k = j kG(x(nT), s) = G w (x) (2.22) where Q (fc)(0 ; x, s) = x. It is extremely important to note that the truncated point mapping for G(x(n), s) is exact up to the k-th order multinomial term. Note that 14 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. because of tbe characteristics of the truncation operation (see equations 2.15-2.17), G (fc) can be computed by manipulating multinomials of degree not larger than k. Hence, for example, for linear approximations one needs to manipulate multino mials of degree one which reduces to m atrix multiplication. 2.2.3 Dependence on Param eters While obtaining the mapping in (2.22), a truncation of the parameter vector s can be performed simultaneously by a scheme similar to the one used for the truncation of the state vector x. The algorithm devised here separately keeps track of the powers of the components of the parameter vector s in addition to the state vector x. Assume that the coefficients of the elements of the parameter vector s of f given in equation (2 .1 0 ) are analytic, i.e. they consist of powers of the components of s. Therefore, we can take the p-jet of the coefficients with respect to the powers of Sj j = 1 ,2,..., L as follows S) = a ^E r,E q{t) 5 i* ' ‘' sl (2.23) 9 where the symbol denotes summation over all sequences kt . ..., kq such that J2m=i — Q and q — 1 , 2 ,p. Define the coefficients of f(t, •, •) as We can write: 15 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. JxJs f (*> X ,S) = VEr&S1) ■■■XN Si 1 sLL (2 -24) r q In order to keep track of the powers of the components of x and perform appropri ate truncation with respect to s, the multinomial telescoping routine is modified. Telescoping in the algorithm is performed by simple multiplication, where an in ternal multiplication scheme operating on the powers of s is implemented and appropriate truncation is performed. Since truncation on parameters is typically carried out to low order, separate truncation of parameter powers significantly speeds up computation. The corresponding truncated point mapping is denoted by x((n + l)T) = G(x(nT),s)fcp (2.25) = G ( M ( x (n T )’s) 2.3 C ell M apping A nalysis We now briefly discuss a general computational technique based on discretization of the state space which can be employed for analyzing the global behavior of a given dynamical system of the type (2.1). We use here the derivation and characteristic developments provided in Hsu [26]. In this formulation, the state space is thought 16 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. of not as a continuum but instead as a collection of cells. Each cell occupies a prescribed region of state space as defined by (2.26). Every cell in the discretized state space or cell state space is to be treated as a state entity. Discretization of the state space in terms of cells is motivated by the fact that in practice state variables can be measured only to a certain accuracy and that round off errors incur while carrying out computation of the state variables using digital computers. These two facts force one to abandon the idea of state space as being a continuum when studying the global properties of dynamical systems from a computational point of view (see Hsu [26]). The cell state space is constructed by simply dividing each state variable com ponent Xi i = 1,2,..., A T into a large number of intervals of uniform size hi. The new cell variables Zi are defined to contain the state variables Xi such that (z* - ^ hi < Xi < (^Zi + 0 hi (2.26) By definition, Z { in (2.26) is an integer. A cell vector z is defined to be an iV-tuple Zi z — 1,2,..., A ^. Clearly, a point x G with components a :,- belongs to a cell z G S with components Zi if and only if x8 - and Z { satisfy (2.26) for all i. The space S consisting of elements which are Ar-tuples of integers is referred to as an A-dimensional cell space. If NC i denotes the number of cells along the x 8-axis, then the cell space S contains a total of NC l X NC 2 x • • • x N C n cells. Figure 2.1 depicts an example of a two-dimensional cell space S and its cellular structure. 17 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. In the cell state space S, one can define a cell-to-cell mapping dynamical system in the form z(n + l) = C(z(n)), C : S - * S , n £ Z+ (2.27) where C is referred to as a simple cell mapping (SCM) and is perceived as a mapping from a set of integers to a set of integers. Equation (2.27) then describes the evolution of a cell dynamical system in TV-dimensional cell state space S. For a detailed treatment concerning the properties of the map C and its refinements, the reader is referred to the research monograph by Hsu [26]. 2.3.1 Periodic M otions The dynamics of the cell mapping system (2.27) is characterized by singular cells consisting of equilibrium cells and periodic cells. The terminology used to de scribe periodic solutions and basins of attraction for a point mapping G(x(n),s) in Section 2.1 also applies to the cell mapping. An equilibrium cell z* is given by z* = C(z*). A periodic cell of period K £ Z + (or a P-K cell) is a set of K distinct cells {z*(k) | k = 1, 2,..., K } such that z*(m + 1) = Cm(z*(l)) m = 1,2,..., iV — 1 (2.28) z*(l) = C *(z*(l)) 18 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. ^2 ..; ...■ ■ ■ / • — < r $ * M M I ' ' , ' ' ' ", ' ' . . ^ X * > ' f* * A ' • * , ' m m i i iM m IM llllli w t # § § i i f » f / s' |l|p: - • • '* - • / J * ^ , 4 / 0 ' 0 < 0 f 0 ^ y 0 S v M J t W H K B ttM |V 1 ' ■ ■ ^ s s ■ 7 : ^ , ' ' - / ■ ' ■ * ' t s , " ' • * ' ^ ' ' } | t | l p j l l p j l l l i i |iiilii|l * ' ' V ; v I s ' , Sint ||M NNM iHP ■ '■ ■ > .'/ s S pi i C 0 $ & 0 & g , ' % l \ ;■ ip ' i ' / ' 1 1 <r ' * ' '' ' ^ V * m im m m ilip ll! , / / / " , , t / , l l l l l l M M M i ' '" , A / t t Figure 2.1: Two-dimensional cell state space. 19 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. where Cm means the mapping C applied m-times. This set is said to constitute a periodic motion of period K or simply a P-K motion. An equilibrium cell in this context is a P-1 cell. A regular cell is one which is not singular. 2.3.2 Creating a Cell Mapping To obtain a cell mapping associated with the periodic dynamical system (2.1), the following procedure may be implemented. First, construct a cell space structure in the state space region with cell size hi in the aq-direction. Let x.^(tp) denote the center point of the cell z(tp) so that x f \ t v) — hiZi(tp) i = 1 ,2 ,..., N. Integrate the state equations (2.1) over one period of time from time t = tp to time t = tp-\-T. Suppose that the trajectory starting from x ^ ( i p) terminates at x ^ (£ p + T ). The cell in which x ^ ( t p + T ) lies is taken to be z(tp + T), the image cell of z(tp). Specifically, Zi(tp + T) — Ci(z(tp)) - Int Xjd)(tP + T) 1 hi 2 i = 1 ,2 ,. .. , A T (2.29) where Int(y), for any real number y, represents the largest integer which is less than or equal to y, i.e. Int(y) < y. This process of finding image cells is repeated for every cell in the cell state space S. The mapping C so determined is a cell mapping associated with the dynamical system (2.1). In this case, the center point 20 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. method (see Hsu [26] and Hsu and Guttalu [28]) was used to compute the image cells. It should be noted that the associated cell mapping C of (2.1) is evaluated by applying the same total integration time duration T uniformly to all the cells in the cell state space S. From a numerical viewpoint, it suffices to apply the rule given in (2.29) to determine an associated cell mapping. Figure 2.2 depicts a two- dimensional simple cell mapping given by C (zi) — > Z 2 , C (z2) — » ■ z2, C(zs) -> z4, etc. 2.3.3 Global Analysis Cells in the cell state space S are of two categories. 1. Singular cell: Either a equilibrium or periodic cell. 2. Regular cell: A cell which is not singular (meaning it is a transient cell). A trajectory of the cell-to-cell dynamical system (2.27) starting from an initial cell state z(0) is referred to as a cell sequence of (2.27) which is the set of integers given by {z(&), k = 1 ,2 ,... }. In practical applications, one is usually interested in a bounded region of the state space which contains a finite number of cells. Once the state variable goes beyond a certain positive or negative magnitude, we are no longer interested in further evolution of the system (2.1). Taking advantage of this, all the cells outside this finite region can be lumped together into a single cell, a sink cell, which will be assumed to map into itself in the mapping scheme. In 21 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. /z. Figure 2.2: Example of a simple cell mapping with periodic (shaded) and transient cells. z 2 is a P-K solution of period 1. z 3 , z 4 , and Z 5 constitute a P-K solution of period 3. R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. practical terms this means that once a cell evolves into the sink cell its trajectory remains locked into the sink cell for all future iterations. ■zsink{n + 1) = C ( z Sink{n)) (2.30) Due to the finite number of cells in state space, all cell sequences must terminate with a finite number of cell mappings into one of the steady cell states, namely an equilibrium cell, group of periodic cells, or the sink cell. This property is of crucial importance in developing the global cell mapping algorithm by Hsu and Guttalu [28]. The algorithm that unravels the global behavior of the system (2.1), provides the following characteristics of the system: 1. Location of the equilibrium states and periodic solutions in a given region of the state space. 2. Domains of attraction associated with the asymptotically stable equilibrium states and periodic solutions. 3. Step-by-step evolution of the global behavior of the system starting from any initial state within the cell state space. Figure 2.2 illustrates several trajectories for a hypothetical two-dimensional simple cell mapping. For this case z 2 is a P-1 cell and Z 3 — > z 4 — > z 5 are P-3 cells. Cells z 1 ; z 6 , and z 7 are regular (transient) cells lying within the domains of 23 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. attraction for the respective P-K solutions. The degree of approximation can be improved by simply reducing the cell size ht (equivalently, increasing the number of cells in S). Depending on the cell resolution unstable equilibrium points and unstable periodic solutions of (2 .1 ) may or may not be recovered. The associated cell mapping C given by (2.27) is viewed as an approximation of the original dynamical system (2.1). The cell mapping C as defined by (2.27), in general, replaces the stable equilibrium points and stable periodic solutions of (2 .1 ) in the state space M .N with sets of periodic cells in the cell state space S. We note that the locations of singular cells of (2.27) provide initial estimates for determining the equilibrium states and periodic solutions of the system (2.1). The cell mapping algorithm avoids repetitive time consuming calculation of trajectories of dynamical systems and is an efficient computational tool (see Hsu [26]). 2 4 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. C hapter 3 E xpanded P oint M apping 3.1 In trodu ction The simple cell mapping (SCM) method has a number of limitations with respect to the global analysis of a nonlinear dynamical system. We have demonstrated in Section 2.3 that under the SCM procedure each cell z , maps to another image cell z j = C ( z , - ) based upon the termination of the trajectory emanating from the center point of the cell z On e immediate consequence is the fact that all points belonging to cell z , are considered to map to the image cell z j = C ( z j ) . We are in effect approximating the continuous time solutions of (2 .1 ) by a sequence of iterations mapping finite regions of state space, in the form of cells, to other finite regions in state space. This leads us to a primary source of error within the SCM approach, namely SCM state space discretization error. We can readily illustrate this particular type of SCM error in Figure 3.1. First we note that cell Zi maps 25 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. to the image cell z 2 , cell z 2 maps to the image cell z 3 , and cell z 3 maps to the image cell z 4 . It is clear that at each mapping step C ( z j ) i — 1,2,3 the trajectory being approximated does not stop precisely at the center point of the image cell zj j = 2,3,4 being mapped to. Hence, a finite error j = 2,3,4 exists at each SCM mapping step. This inaccuracy can lead to erroneous results depending on the manner in which the error propagates for the actual nonlinear system being studied. For certain systems it may force one to utilize small cell sizes which in turn can impose computational difficulties with respect to memory requirements. The fundamental issue that needs to be addressed is that all points belonging to a particular cell do not necessarily map to the same image cell. In addition to SCM spatial discretization error there are also a number of secondary SCM shortcomings as well. For example, the SCM technique does not possess the analytical capability to determine the local stability properties of a particular P-K motion. Furthermore, the SCM method has difficulties studying the behavior of chaotic systems. Under chaotic conditions the SCM method cannot reconstruct arbitrarily long non-periodic trajectories. The best it can do under these circumstances is to indicate a P-K motion of relatively high periodicity. These observations with regard to simple cell mapping provided the initial impetus to develop the Expanded Point Mapping approach. In this chapter, the Expanded Point Mapping (EPM) methodology is introduced. Chapter 3 is divided into two parts. First, the Expanded Point Mapping method is formulated and its 26 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Figure 3.1: Simple cell mapping spatial discretization error. R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. theoretical justification discussed. Then a procedure for EPM-based global analysis and local analytical study of periodic solutions of nonlinear periodic systems is described. 3.2 E P M Form ulation 3.2.1 State Space D iscretization The region of interest in the state space is discretized into a finite number of cells Nc + 1 in the manner described in Section 2.3. Thus forming a cell state space of dimension Nc — NC l x NC 2 x • ■ • x NC n and a corresponding sink cell. 3.2.2 Cell Space Analysis Define 4>*(t) € to be a trajectory starting at the center point of each cell z* • for i — 1, 2,, Nc- As in the simple cell mapping approach the center point of each cell is integrated from t = tp to t — tp + T with tp — 0 for simplicity. Each cell center point is represented by <£*( 0 ) and the corresponding cell center point trajectory termination is < f> *(T). 28 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. 3.2.3 A nalytical Study of Trajectory Evolution The truncated point mapping approach is used to precisely locate trajectories inside each individual cell. The perturbation (error) about is denoted as £i(t) (see Figure 3.2) where x(*) = # (* ) + &(*) (3-1) Substitution of x(t) into equation (2.1) yields #i (*) + &(*) = f 0 > * 0 0 + &00>s) (3-2) Assuming that f is analytic, we apply a Taylor series expansion about to the right hand side of (3.2). Recalling that ) is a solution of equation (2.1), the perturbed trajectory satisfies ki(t) = P;(m)(*,&(*),s) + R i{m+i)(t,^{t),s) (3.3) where p ,-(m )(t, £s -(f), s) is a vector homogeneous multinomial of degree m in £,-(t). Notice that (3.3) is written as a perturbation about a trajectory. Consequently, it does not possess terms of order zero, and satisfies the hard forcing criteria of Section 2.2.1. The remainder term Rqm+i)(f, £-(t), s) is a higher order multinomial 29 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. of 0 (||& (*)ir+1)- In some neighborhood of for which |[£;(f)|| < p, the truncated system &(*) = i = l,2 ,...,N c (3.4) provides a good approximation to the dynamics of the perturbed system. A trun cated point mapping for each individual cell z; is obtained by integrating (3.4) from t = tp to t = tp + T where again tp — 0 for simplicity. U T ) = G i(m)(&(0), s) i = 1 , 2 ,..., Nc (3.5) Hence, the analytical m-th order approximation of the trajectory is given by: x(T) = cf>*(T) + G t -(m)(&(0), s) * = 1 , 2 ,..., Nc (3.6) Note that equation (3.5) defines the mapping of all points within the cell z, • and is accurate to degree m. This is the fundamental tenet of the Expanded Point Mapping approach. Defining the norm || • || to be the infinity norm, i.e. = Hiax|&rf-(*)l j = 1,2,..., A (3.7) with an appropriate scaling, then ||£,(0)||oo < o * defines a rectangular region. This rectangular region is considered by Hsu [26] as a definition of a cell occupying a 30 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. x ( k = l ) Cellm ^,*j(0) 1 m^k_2> m ' Cell i Cellp x ( k = 0 ) Figure 3.2: Two EPM mapping steps. 31 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. region of state space in the same spirit as previously developed in Section 2.3. Therefore, £,-(T) defines a map of any point within the entire cell z,. It is critical to reiterate the fact the map £t -(T) can be evaluated to any order of multinomial approximation m and can be developed for systems of arbitrary dimension N. It should be noted that the mapping defined in equation (3.6) is an alternative to the Interpolated Cell Mapping (ICM) procedure (see Tongue [48] and Tongue and Gu [49, 50, 51]) for determining the mappings of points within each cell. Furthermore, the mapping defined in (3.6) has the following advantages. 1. It is defined for high order systems while the extension of the Interpolated Cell Mapping for systems of order greater than two is not straightforward. 2 . The degree of approximation (interpolation) is arbitrary as opposed to Inter polated Cell Mapping for which higher order approximation, especially for high order systems is difficult. 3. There is no need to have special tests for defining the interior of a cell which is a major handicap of the ICM method. It seems that for higher order systems this determination is virtually impossible. 4. The mapping within the cell is parameterized analytically by a design pa rameter. No such feature exists in the ICM approach. 32 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. 3.2.4 Local Stability The analytical expression for the truncated point mapping (3.5) permits us to study the stability properties of the Expanded Point Mapping. The local stability char acteristics of periodic solutions can be obtained by using the point mapping theory developed in Section 2.1. By definition, a periodic solution (x(0), x(T),..., x(fcT)) of period K T satisfies: x(0) -► x(T) (3.8) x(T) x(2T) x((/T-l)T) -> x(KT) x(isTT) = x(0) Consider Gqm) (£;(&), s), it will consist of a linear part and subsequent higher order terms up to degree m. G i W « j(fc),s) = H i( s ) & (*) + ••■ (3.9) Applying Liapunov’s indirect method and its extension (see Jury [29] and Vidyasagar [54]) to the P-K solution x(j) j = 1,2,..., K yields 33 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. K H ( s ) = n H»w p - 10) i = i where the local stability of x(j) is determined by the eigenvalues of H(s). There exists a general relationship between linearized flows and Poincare maps (see Guckenheimer and Holmes [13]) using Floquet theory. For the nonlinear system x = f(f, x, s) the characteristic multipliers for the linearized flow are identical to the eigenvalues of the linearized Poincare map. Consequently, the characteristic multipliers of the linearized system, « •■ (< ) (3-11) X = 0 * (t) will be identical to the eigenvalues of the linear part of the cell truncated point mapping £(T) = G i(1)& (0 ), S) = H,-(s) &(0) (3.12) 3.2.5 Bifurcation Analysis The truncated point mapping analytical expression (3.5) also contains functional relations between system parameters. Therefore, the Expanded Point Mapping approach can be used to perform bifurcation studies by writing the functional 34 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. A * [ H ( s ) ] i = l,2,...,iV (3.13) Bifurcation produces new solutions from the existing solutions when the system parameter vector s changes its value. If for some s = s0 | A , - [ H ( s ) ] | = l ( 3 . 1 4 ) the periodic solution at hand may lose stability and new solutions may come into existence. We are interested in determining the value of the system parameter s at which a bifurcation from a periodic solution may take place and the conditions under which such a bifurcation is possible. Let us suppose an existing P-K solution is stable and represented by x* | * = 1,2,...,K : s = s0 (3.15) As the parameter s is varied the stable P-K solution (3.15) may become un stable producing new periodic solutions whose periodicity may be different from the original periodic solution x*. The new solutions, initially stable, may again become unstable as s changes giving birth to yet another set of periodic solutions. This process of bifurcation may continue without an end. In general, a given P-K solution may bifurcate into other P-MK solutions for M = 1,2, ■ • ■ , and the 35 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. original stable solution may become unstable. These new solutions can be analyt ically studied by writing a m -th order approximation of the periodic solution for s = s0 + p where ||p|| is small. 3.2.5.1 Bifurcation of P -K solutions into P -K solutions One convenient means of obtaining the necessary conditions for bifurcation of periodic solutions is to examine the eigenvalues of the stability m atrix H (p) locally at the P-K solution {x* | * = 1,2,..., K : s = So}. First we create a mapping G,-(m ig)(zi(n); p) of degree m in Z; and degree q in p where x,-(i) = x*(f) + z;(f) * = 1,2,..., if (3.16) s = S o + p (3.17) using the truncated point mapping procedure previously developed in Section 3.2.3. For EPM bifurcation analysis we expand about the P-K solution x*. Additionally, we now include parameter dependency by perturbing the parameter vector s about s0. It is important to note that the mapping G,-(T O i 3)(z;(n); p) is obtained via a Taylor series expansion about x* only. The parameter dependency p is obtained by direct substitution of (3.17) into (2.1). G,-(m > g )(z;(n); p) is a m -th order mapping about x* and a g-th order mapping about So, 36 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Zi(ra + 1) = G < ( m i , ) ( z i ( n ) ; p ) i = 1,2,..., h f (3.18) where z,-(n) and z,(n + 1 ) are perturbations of the system at t = nT and t = (n + 1 )T, respectively. By definition a bifurcation to another P-K solution is possible when z i(n + K) — z i(n) (3.19) By applying the mapping in equation (3.18) K times and combining this result with equation (3.19) and rearranging we obtain the following relationship. Z i ( " ) - G f m,,)(zi(n);p) = 0 * = 1,2 (3.20) A bifurcation may take place when the system of equations (3.20) has multiple solutions. Application of the Implicit Function Theorem (see Marsden [32]) to the set of equations (3.20) results in this condition for nonunique solutions. det [I - J (G f mi?)(z4 -(n); p ))] = 0 i = l ,2 , .. ., K (3.21) where J is defined as the Jacobian. Using the composite function rule for the Jacobian J we can express (3.21) as 37 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. det [ I - H ( p ) ] = 0 (3.22) where K K h (p ) = n J h - , . ) ! 2'!");? )] = t3-23) i~ 1 i= 1 Examining equation (3.22) it is obvious the condition for bifurcation of a P-K solution into other P-K solutions is to require one of the eigenvalues of H ( p ) to equal +1. 3.2.5.2 Bifurcation of P -K solutions into P-2K solutions Next, we consider the possibility of bifurcation of a P-K solution into other P-2K solutions. Using the same line of reasoning developed in equations (3.18-3.23) we can obtain Zi(n) - G f miq)(zt(n)-p) = 0 i = l ,2 , .. ., K (3.24) As before, using the Implicit Function Theorem the condition for multiple so lutions is det [ I - H 2 ( p ) ] = det [ I - H ( p ) ] det [ I + H (p)\ = 0 (3.25) 38 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Since det [I — H(p)] = 0 leads to a bifurcation to another P-K solution the condi tion for a bifurcation to a pure P-2K solution will be det [I + H (p)] = 0 (3.26) Clearly, this condition requires one of the eigenvalues of H(p) to be equal to — 1. Hence, a bifurcation from a P-K solution into other P-2K solutions is possible when one of the eigenvalues of H(p) moves across the value of — 1 as p is varied. 3.2.5.3 Bifurcation of P -K solutions into P -M K solutions In general, we can study the possibility of bifurcation into a P-MK solution for M = 1,2, •• • where equation (3.20) is generalized to Zi(n) - G ^ g)(zi(n); p) = 0 i = l ,2 , .. ., K (3.27) Again the conditions for bifurcation are the existence of multiple solutions. Application of the Implicit Function Theorem to (3.27) results in d e t [ I , ( *( ») ; ?) ) ] =0 i = l,2,...,K ( 3 . 2 8 ) By the composite function rule (3.28) may be reduced to N det [I - HM(p)] = JJ(1 - A f) = 0 (3.29) i=i 39 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. where A j j = 1,2,..., N are the eigenvalues of the m atrix H (p) associated with the P-K solution. It follows from (3.29) that the condition for bifurcation from a P-K solution into other P-MK solutions requires one of the eigenvalues of H(p) to satisfy A f = 1 \j = exp2 ™ p/M, * = v /=T, 1 ,2 ,...,M - 1 (3.30) In (3.30), the case p = M is deleted since it corresponds to A j = 1 which is the condition for the bifurcation of a P-K solution into other P-K solutions already considered in Section 3.2.5.1. Also, the periodicity of the solution requires that the ratio p/M be an irreducible number. Therefore, all the values of p for which p/M is not an irreducible number should be excluded from (3.30). 3.3 E P M P rocedure Based on Section 3.2 a trajectory x(fc) with t = kT is computed for any initial condition of (2.1) in the region of interest as follows: 1. Discretize the state space region of interest into Nc cells. 2. Find the mapping of each cell center point <f)*(T). 40 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . 3. Compute the m-th order approximation (£,•(&), s) of the mapping inside each cell. 4. Set tp — 0 — y k = 0. 5. Calculate the perturbation £*(&) = X(A) ~ 0*(O). 6. Apply the perturbation £;(&) to the truncated point mapping (3.5) and ob tain €i(k + 1) = Gi{m)(^(k),s). 7. Calculate where x(fc) maps to by using x(fc + 1) = + £,(& + 1). 8. Increment k — t k + 1 and repeat steps 5 thru 8 until the trajectory x(fc) satisfies one of three possible outcomes. (a) The trajectory x(fc) maps to the sink cell. (b) The trajectory x(fc) evolves into a P-K solution. (c) The trajectory x(fc) remains within the region under study and exhibits chaotic behavior, i.e. it is aperiodic. 9. Trajectory evaluation is performed for initial conditions of interest using steps 4-8 above. 10. Local stability analysis is performed by evaluating H(s) along the trajectory for s = S o- 41 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . 11. Bifurcation analysis is performed by studying the eigenvalues of H(s) as a function of s. Figure 3.2 depicts how the EPM algorithm operates for a two-dimensional system. For this case we perform the first two mapping steps k = 0,1. Details for the main features of the actual source code used in the implementation of the EPM method are contained in Appendix A. 3.4 E P M Sum m ary The end result of the EPM procedure is a comprehensive evaluation of global (cell mapping) and local (point mapping) analytical expressions. Namely • Each cell center point trajectory and the cell it maps to. • The analytical m-th order mapping Gqm) (£,-(()), s) (equation (3.5)) of all points within each cell, which permits accurate trajectory evaluation while performing global analysis on the system. • The coefficients of the truncated point mapping G t -(m)(£t -(0), s) for each cell, which are used to perform local stability and bifurcation condition studies. Figure 3.3 illustrates a cell truncated point mapping G 1(m)(^1(0), s) resulting from the EPM procedure. All of the points within cell Zi (dark shaded region) map to the lightly shaded EPM mapping region. Note that the EPM map partially 42 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . covers several cells, Z 2 ,...,Zg, and is a curvilinear region. For comparison the simple cell mapping is shown, where the SCM map is the crosshatched area within cell Z3. X9 EPM Mapping: Cell Zj SCM Figure 3.3: EPM mapping for cell Zi. 43 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . C hapter 4 P endulum w ith V ertical Support M otion An approach for analyzing the behavior of nonlinear periodic systems has been developed in Chapter 3. For our first example, we will demonstrate the application of the Expanded Point Mapping method on a pendulum subjected to an excitation in the vertical direction only. 4.1 S ystem D ynam ic E quations Consider a pendulum of mass m, length L , and damping c at its hinge as shown in Figure 4.1. The vertical support applies a periodic excitation of the form u(t) = Asm.{u>ji) (4-1) 44 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . where A is the amplitude of the vertical input and tVf is the forcing frequency. The angular displacement and angular velocity are denoted as the state variables Xi = 6 and x2 = 0, respectively. The state equations of the pendulum are given by xi = x2 (4.2) x2 = —2(aiX2 — (aq — a 2 sin 21) sin Xi where aq = 2 — a2 = ~ L o 2 n = j = 2& n 21 = uf t (4.3) Ljf L L mLz and L o n is the system natural frequency, g is the acceleration of gravity, £ is the rotational damping coefficient, and the argument of the periodic excitation u(t) has been transformed to 21. 4.2 B ifurcation P relim inaries The first step for this EPM bifurcation analysis is to determine the mapping about the perturbed solution x = (f)* + z. From the preceding chapter zi(t) can be 45 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . u(t) 8 W t. m m m Figure 4.1: Pendulum with a vertically moving support, approximated by a vector homogeneous multinomial of degree m in z, and degree r in p, namely — P j ( m ,r ) ( A p ) I — 1 , 2 , . . . , A ( 4 . 4 ) using the process developed in Section 3.2.5. Note that K is the order of the P-K solution under study. Recalling that 4 6 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . = 4> *0) + Z i(X ) S = Po+P (4.5) with 4>*(t) representing the solution to (4.2). The vector homogeneous multinomial Pi(m,r)(T z,-(<); p) is the following Taylor series expansion of f(t, x(f), s). K t t ) = L / i ( * > < £ * ( / ) ; Po + p ) ~ /i(*5 < t > * ( < ) ; Po)\ + Z2,i{t) = IM^ # ( * ) ; Po + p) - # ( * ) ; Po)] + -9/i d / 2 * = * !■ (* ) Z 2 ,- s = p 0+ p d f2 dx-; x= 4?(t) s = P 0+ P ^2 ,« + 1 d2/ 2 2! 5a;? dx\ 1 d3f 2 X = ^ f(t) S = P 0 + P 1 d4/ 2 4! da;} x=0*(t) ~ 1,' S = P 0+ P + + x = 0 f(t) S = P o + P 1 9n /2 3! da;? (4.6) zi,i + (4-7) x = 0 t «(t) z i , i + S = P 0 + P 11! da;}1 z1 1 f(t) *i,z S = P o + P Evaluation of the individual terms in equations (4.6-4.7) results in /i(*>*W;po + p) /i(T </*(*); Po) dfi dxi dfi dx2 djfi dx 0 0 1 0 (4.8) for j > 2 4 7 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . f 2 {t,(f>*(t);p0 + p) - / 2(f,$■(*); Po) = sm 2tsm xlP (4.9) 9 h , cbi <^2 dx 2 (a 2 sin2t — c t ^ ) c o s - 2 £ a i 0 for j > 2, k > 0, and / > 1 dx\dx l 2 r- = — (ct2 sm 2r — a x) sm xi »Uf2 dxl1 -(a 2 sin2f — a j) cosxi Substitution of (4.8-4.9) into (4.6-4.7) leaves us with a vector homogeneous multinomial of degree m — 11 in z* and degree r in p. z\,i(t) = ^2,i (4-10) = [ ( /= » o sin 2 t - a \ ) cos )] 2i,t - + [sin 2 t sin ( f > * hi(t)] p + [sin 21 cos <^,-(2)] zi,i p — 2£cri z2,i + ^ [ - ( / jo sin 2t - al) sin# > t -(*)] z\y i + ^ [ - sin2t sin #,,•(*)] 4,i P + [-(/)o sin 2 t - al) cos ^ -(f)] sJJ + [~ sin 2t cos <£*,.(*)] z\)i P 48 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . For this particular problem the parameter dependency is determined for a 2 by the relationship s = a.2 = Po + p (4-11) The TPM algorithm is applied to (4.10) in order to evaluate the mapping of interest, namely Zi(T) = G i(m,r)(Zi(0); p) i = 1 ,2 ,..., K (4.12) where T — n is the period of (4.2). It is vitally important to note that in this study Q !2 was varied between 0 to 2. This corresponds to a large motion of the support input. Therefore, any attem pt to find a solution by a perturbation method is not possible. 4.3 P - l to P-2 B ifurcation We proceed to investigate the behavior of the system starting with the trivial solution £i(0) = £ 2 (0 ) = 0 while the parameter a 2 is varied. The Expanded Point Mapping analysis was carried out for these fixed values of system parameters in (4.2). £ = 0 ot\ = 0.1 (4.13) 49 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm ission . First, we take the Jacobian J of equation (4.12) and determine the value of a 2 where a P-l to P-2 bifurcation may occur by examining the relationship det [I + H(p)] = 0 p - p - a 2 — po (4.14) The parameter order was set to r = 10 for p in equation 4.14. The trivial solution bifurcated into a single stable P-2 solution (A i,A 2) at a 2 = 1.64712 as depicted in Figure 4.2. As for the trivial solution, it became unstable. For this case po was set equal to 1.0 initially, thereby demonstrating the capability of the EPM method to accurately locate a bifurcation point for large values of the parameter p. The branches A\ and A2 were determined by solving the analytical expression z; - G qm,r)(z;; p) = 0 * = 1 , 2 (4.15) where G^mrj ( z p ) is the mapping (4.12) composed once. The mapping G?^m r^ ( z p ) was initially created by expanding about the trivial solution with p0 = 1.64712, a state order m = 11 for z a n d a parameter order r — 12 for p. Matlab (Version 6.0 Release 12) was the solver of choice used to determine various solutions of (4.15). The parameter p was varied between 0 to 0.08 and z* was determined from (4.15) and substituted back into x 4 - = ( f> * + z i i = 1,2 (4-16) 50 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Bifurcation to P - 2 solution: £ = 0, » 0.1 0.75 0.5 0.25 -0.25 -0.5 A2 -0.75 C2 2.2 2.3 2.1 1.7 1.8 1.5 a. “ 2 Expanded Point M apping Method Figure 4.2: Bifurcation to P-2 solution, £ = 0 and a\ = 0.1. to determine the newly created P-2 solutions as shown in Figure 4.2. At p = 0.08 the mappings G qT O ) r )(z*;p) i = 1,2 were created by expanding about the non trivial solution x t - i = 1,2 previously calculated at a 2 = 1.72712. For this step and all subsequent iterations the state order used was m = 7 and the parameter order used was r = 8 for Gqmjr)(z;; p) * = 1,2. The bifurcation solution branches were constructed by varying o:2 in increments of 0.08. At a 2 = 2.08302 the P-2 solution (Ai, A2) bifurcated into two stable P - 2 solu tions (Bi,B2) and (Ci,C2) while the P-2 solution (Ai, A2) itself became unstable. The location of this P-2 to P-2 bifurcation was determined precisely by solving the relation 51 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm ission . det [I - H (p)] = 0 (4.17) where K K H (p) = n J P)] = I I K = 2 (4-18) t = l j = 1 As a 2 was further increased to 2.12722 the two stable P-2 solutions , B2) and (Ci,C 2 ) became unstable. Figures 4.3-4.5 represent the continuous trajectories associated with the three different P-2 solutions at various values of a2. The solutions (J3i, B2) and (Ci, C2) were symmetrical with respect to each other. The determination of the stability character for the P-2 solutions was addressed by examining the eigenvalues of the stability matrix H (p) using equation 4.18. Finally, with regard to bifurcation solution accuracy the EPM method was extremely effective. The solution trajectories were accurate to six significant digits which is quite remarkable. Clearly, this demonstrates that the EPM bifurcation analysis was very robust with respect to solution accuracy over a large parameter variation. 52 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . 1 0.5 -0.5 -1 Figure 4.3: A stable P-2 trajectory associated with the solution (Ai, A 2 ) for £ = 0, af = 0.1, and « 2 = 2.00102. Stable P-2 trajectory: % = 0, oq = 0.1, a2 = 2.001 T ------------------------------1 ------------------------------ 1 ------------------------------ 1 -------------- A2 J _______________ L _ .... ..... I. I _______________ L 53 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Stable P-2 trajectory: % = 0, = 0.1, a 2 = 2.126 1.5 1 B2 0.5 0 -0.5 1 -1.5 L -1.5 1 1.5 0 0.5 ■ 1 -0.5 * 1 Figure 4.4: A stable P-2 trajectory associated with the solution (P i, P 2) for £ = 0, aj = 0.1, and a 2 = 2.12602. 54 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . 1 0.5 > ? 0 -0.5 -1 ’1'-U5 -1 -0.5 0 0.5 1 1.5 * 1 Figure 4.5: A stable P - 2 trajectory associated with the solution (Ch, C2) for £ = 0, af = 0.1, and a 2 — 2.12602. Stable P-2 trajectory: \ = 0, oq = 0.1, a2 = 2.126 -------------- 1 -------------------- 1 -------------------- 1 -------------------- r C2 : : J_______________ I _______________ L 55 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . 4.4 P - l to P -3 B ifurcation Again we begin the study by ascertaining where the trivial solution aq (0) = * 2 (0 ) = 0 bifurcates to new solutions. In this case the governing relationship for a P-l to P-3 bifurcation is det [I + H (p) + H 2(p)] = 0 (4.19) The trivial solution bifurcates into four P-3 solutions at a 2 = 1-35158. These solutions are denoted as A — (A*), B — (Pi), C = (Ci), D = (D i), * = 1,2,3 and are shown in Figures 4.6-4.7. For this case po was initially set equal to 1.0, again demonstrating that one does not have to be particularly close to a bifurcation point to locate it accurately. Again the parameter order was r = 10 for p in equation 4.19. The four branches were determined by solving the analytical expression zi - G \ mtr)(zi-p) = 0 t = 1,2,3 (4.20) The mapping G ^m r^ ( z p) was initially created by expanding about the trivial solution with po = 1.35158, a state order m = 11 for z a n d a parameter order r = 12 for p. At p = 0.03 the mappings GqT O jr)(z,-;p) i = 1,2,3 were created by expanding about the non-trivial solution x,- i = 1,2,3 previously calculated at a 2 = 1.38158. For this step and all subsequent iterations the state order used was m — 7 and the parameter order used was r = 8 for the mappings GqT O jr)(z;;p) 56 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Bifurcation to stable P-3 solutions: % - 0, = 0.1 0.75 R4 0.5 S4 0.25 S6 S3 C3 R5 R2 B2 -0.25 B3 R 6 -0.5 R3 S2 C2 -0.75 S5 1.2 1.3 1.4 1.5 1.6 1.7 1.8 Expanded Point Mapping M ethod Figure 4.6: Bifurcation to stable P-3 solutions B and C, £ = 0 and af = 0. R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Bifurcation to unstable P-3 solutions: % = 0, = 0.1 0.75 0.5 0.25 D2 A3 -0.25 A2 -0.5 -0.75 D3 1.8 1.9 a, 2 Expanded Point Mapping Method Figure 4.7: Bifurcation to unstable P-3 solutions A and D, £ = 0 and = 0.1. 58 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . i = 1,2,3. The parameter p was varied between 0 to 2.0 and Z; was determined from (4.20) and substituted back into Xi = 4>*i+Zi * = 1,2,3 (4.21) Once the various P-3 solutions were found for each branch the continuous tra jectories associated with them can be readily obtained. For illustrative purposes, Figures 4.8-4.11 depict these trajectories at a 2 = 1.62158. Trajectory stability was determined for each branch A = (Ai), B = (Bi), C = (Ci), D = (Di), i — 1,2,3 by formulating the stability m atrix H (p) K K H(p) = JJ J [Gi(m,r)(zi;p)] = J ] H ^ ) K = 3 (4-22) 4 = 1 « = 1 where J is the Jacobian. The solutions B and C were determined to be stable while solutions A and D were unstable. It is important to note that although the new solutions emanate from the trivial solution, the trivial solution itself remains stable. As oi 2 was increased the two stable P-3 solutions B and C became unstable at < * 2 = 1.63598 and bifurcated into two stable P-6 solutions R = (Ri), S = (Si), i = 1,..., 6. The two P-6 solutions R and S remained stable until « 2 = 1-64778 where they lost stability and bifurcated into two P-12 solutions. The continuous trajectories for the two P-6 solutions R and S a t ct2 = 1.64008 are shown in 59 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Figures 4.12-4.13, respectively. The precise locations for both the P-3 to P - 6 and P - 6 to P-12 bifurcations were determined by solving the relation det [I + H ( p ) ] = 0 (4.23) For this particular analysis the solution trajectories were accurate to five sig nificant digits over a large parameter range, i.e. a 2 was varied between 0 to 2 .0 . Again this demonstrates the robustness of the Expanded Point Mapping method with respect to large parameter variations. 4.5 E P M A n alytical V erification We recall for this particular system under study, i.e. equation (4.2), the rotational damping coefficient £ was set equal to zero. Therefore, in this case equation (4.2) is a Hamiltonian system since it satisfies the condition = 0 (4.24) An important property of Hamiltonian systems is they are volume preserving in state space. Consequently, this condition leads to the equality = det H (p) = 1 (4.25) 60 det 8 G < 9 x tr d t dx. R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Stable P-3 trajectory: £ = 0, = 0.1, a2 = 1.6216 0.5 B2 -0.5 -1.5 0.5 -1.5 -0.5 Figure 4.8: The stable P-3 trajectory associated with the solution (Bx, B 2 , B 3 ) for £ = 0, a\ = 0.1, and c * 2 = 1.62158. 61 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Stable P-3 trajectory: \ = 0, = 0.1, a 2 = 1.6216 C2 0.5 C3 -0.5 -1.5 L -1.5 -0.5 0 0.5 1 1.5 1 Figure 4.9: The stable P-3 trajectory associated with the solution (Ci, 6 * 2 , C3 ) for £ = 0, a± = 0.1, and 0 : 2 = 1.62158. 62 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm ission . U nstable P-3 trajectory: \ = 0, ot^ = 0.1, a g = 1.6216 0.5 A3 A2 -0.5 -1.5 0.5 -0.5 -1.5 Figure 4.10: The unstable P-3 trajectory associated with the solution (Ai, A2, A3) for £ = 0, af — 0.1, and 0 : 2 = 1.62158. 63 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Unstable P-3 trajectory: 5=0, = 0.1, a g = 1,6216 D3 0.5 D2 -0.5 -1.5 -0.5 0.5 -1.5 Figure 4.11: The unstable P-3 trajectory associated with the solution (Hi, H2, 0 3) for £ = 0, af = 0.1, and a 2 = 1.62158. 64 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Stable P-6 trajectory: % = 0, ct^ = 0.1, a 2 = 1.6401 1.5 R3 R6 0.5 R5 R2 M / / -0.5 -1.5 -1.5 -0.5 0.5 Figure 4.12: The stable P - 6 trajectory associated with the solution (Z?i,... ,Rq) for £ = 0, a\ — 0.1, and a 2 = 1.64008. 65 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Stable P-6 trajectory: £ = 0, = 0.1, c * 2 = 1.6401 S5 0.5 S2 S6 S3 -0.5 S4 -1.5 0.5 -1.5 -0.5 Figure 4.13: The stable P-6 trajectory associated with the solution (Si,..., 56) for £ = 0, af — 0.1, and « 2 = 1.64008. 66 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm ission . Hence, this relation can be used to assess the accuracy of the proposed EPM method. The determinant of H (p) was calculated at a 2 = 2.00102 for the stable P-2 solution (Ai, A2) to be detH (p) = 0.9999999966 £ = 0.0, a? = 0.1, a 2 = 2.00102 (4.26) The results in (4.26) have an error of 3.4E— 09 which for all practical purposes can be considered exact. 4.6 G lobal A nalysis To demonstrate the global analysis capability of the EPM procedure we introduce a finite value for the rotational damping coefficient £. When £ > 0 the behavior of the system (4.2) changes qualitatively. A bifurcation into a P-3 solution is impossible, because the eigenvalues of H (p) cannot reach the unit circle in the complex plane. However, this does not preclude the existence of P-3 solutions. As in Section 4.4 we created a mapping zi - G ? ,rajr)(zi;p) = 0 . = 1,2,3 (4.27) 67 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . and found that P-3 solutions suddenly appear when a 2 reaches a certain value as shown in Figure 4.14. Damping was set equal to £ = 0.01 (4.28) In this case, four P-3 solutions, two stable B = (A), C — (Ci) i = 1,2,3, and two unstable, A = (A'), D = (A), * = 1,2,3, appeared at a 2 = 1-46758. As a 2 increases, the various P-3 solutions behaved in a similar manner as in the undamped case in Section 4.4. The two stable P-3 solutions B = (A ), C — (Ci) remained stable until o 2 = 1.63718 and then bifurcated into two stable P-6 solutions while B = (P»), C = (Ci) themselves became unstable. The continuous trajectories of the four P-3 solutions A = (A), & — (A), C = (Ci), D = (Di) at Q '2 = 1.59158 are shown in Figures 4.15-4.18. For this analysis the solution trajectories were accurate to six significant digits over a large param eter range, i.e. a 2 was varied from 1.46758 to 2.0. This again demonstrates the ability of the Expanded Point Mapping method to accurately locate solutions over a wide parameter span. 68 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm ission . Bifurcation to stable P-3 solutions: % = 0.01, = 0.1 0.75 C1 0.5 0.25 D2 A3 -0.25 A2 -0.5 -0.75 D3 1.8 a, Expanded Point M apping Method Figure 4.14: P-3 solutions A, B, C, and D for £ = 0.01 and af = 0.1. 69 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Stable P-3 trajectory: \ = 0.01, o c^ = 0.1, a g = 1.5916 B3 0.5 B2 -0.5 -1.5 -0.5 0.5 -1.5 Figure 4.15: The stable P-3 trajectory associated with the solution (Bi, B2, B3) for £ = 0.01, a\ = 0.1, and a 2 — 1.59158. 70 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Stable P-3 trajectory: % = 0.01, ct^ = 0.1, a 2 = 1.5916 C2 0.5 C3 -0.5 -1.5 0.5 -0.5 Figure 4.16: The stable P-3 trajectory associated with the solution (Ci, C2, C3 ) for ( = 0.01, af = 0.1, and a 2 = 1.59158. 71 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Unstable P-3 trajectory: | = 0.01, = 0.1, o c2 = 1.5916 0.5 A3 A2 -0.5 -1.5 0.5 -0.5 -1.5 Figure 4.17: The unstable P-3 trajectory associated with the solution (A\, A2, A3) for £ = 0.01, af = 0.1, and a 2 = 1.59158. 72 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm ission . Unstable P-3 trajectory: \ = 0.01, = 0.1, a g = 1.5916 D3 0.5 D2 -0.5 -1.5 0.5 -1.5 -0.5 Figure 4.18: The unstable P-3 trajectory associated with the solution (P i, D 2 , D 3 ) for £ = 0.01, af — 0.1, and « 2 = 1-59158. 73 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . 4.6.1 Basins of A ttraction The next step will be a study of the system global behavior at a specific parameter value a 2, namely < x 2 = 1-59158 (4.29) For this EPM analysis the region of state space studied in R 2 was -7T < X\ < 7 T -7T < X2 < 7T (4.30) Using the procedure developed in Section 3.3 we performed a fifth-order EPM analysis of the system (4.2). The cellular resolution was 101 x 101 in the phase plane (equivalent to 10,201 cells in cell state space S ). A truncated point mapping, G,-(5)(£4 (0), s ) i = 1,2,..., 10201, was created for each individual cell. EPM trajec tory evolution was then performed for an initial condition resolution of 451 X 451. Figure 4.19 shows the EPM domains of attraction. The initial conditions in the phase space were divided into three categories. Specifically, trajectories that map to the two stable P-3 solutions B = (Pi), C = (Ci) represented by the symbols (o, + , *, and x), trajectories that map to unrelated P-K solutions (lightly shaded region), and trajectories that map to the sink cell (blank regions). 74 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Expanded Point Mapping Method Cell TPM order is 5. Cell resolution is 101 by 101. Cell TPM integration stepsize is 3.1416 /100. Trajectory initial condition resolution is 451 by 451. : ■ /> 9* ' - c ° c 0.8 0.6 0.4 0.2 - 0.2 -0.4 - 0.6 - 0.8 * > ■ ■ o « : ■ . : W 9 el * 0 S • ’ O O ° « eg ® o g A ■ 1 k W> £ * < 4? 9 ^ 0 $ . » o <o V o ;« o • ° c o . . v w . . * » ® ■ oft* % ■ ? O ^ .-.v.v.v.v.v. ^ < * > < * * > • o * "ftT ° s * » x ; O ° V*» 0 < * o # ■ ® ° c > £ . . -O .o . 1 _l_ _L -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 X, / it Basins of attraction for the P-K solutions 0.6 0.8 Figure 4.19: P-3 EPM Basins of Attraction. R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . 4.6.2 A nalytical Verification For periodic systems with linear viscous damping it can be shown that det dG dx — det H = exp(— 2£ct1Tl) (4.31) where T = K ■ 7 T is the periodicity of a P-K solution for the system (4.2). Substi tuting « i, and T into (4.31) for a P-3 solution yields Xi ■ A 2 = exp(— 2 • 0.01 • VO1 • 3tr) (4.32) Ai • A 2 = 0.94213422 to eight significant digits. Since Ai ■ A 2 is determined analytically, it is possible to verify the multiplication of the calculated eigenvalues using the stability matrix of a P-K solution. K H = n H i £ = 0.01 a\ = 0.1 K = 3 (4.33) i = 1 For the stable P-3 solution represented by B = (Bi) at « 2 = 1.59158 in Fig ure 4.15 we obtain 76 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. -1.2244748 -1.3005225 1.6894402 1.0249463 (4.34) where the eigenvalues of H<Bd were A 1 )2 = -0.099764215 ± 0.96549537 j with j — a/— 1. In this case Ai • A 2 = 0.94213421, which is identical to (4.32) to seven significant digits. These are very accurate results and can be considered as exact for all practical purposes. 77 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . C hapter 5 P end u lu m w ith P eriod ically E xcited Support in th e P lan e Pendulums with an oscillating point of support have been studied extensively by a number of researchers. Typically most pendulum studies have considered only vertical support excitations (see Phelps and Hunter [38], Struble [43, 44, 45]). However, there are numerous other types of support motion and pendulum type. Dadfar and Geer [6] studied horizontal excitations, a flexible pendulum was ana lyzed by Ryland and Meirovitch [41], and a pendulum with a support performing circular motions was examined by Schmidt [42] to name just a few examples. For all of these cases the resulting equation(s) of motion included some combination of periodic forcing and parametric excitations. Chapter 5 is organized as follows. The first section will provide the governing equation of motion for the pendulum under study. Then the details of the actual EPM example problem setup will be reviewed. Next EPM computational studies 78 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . analyzing both the global behavior and bifurcation characteristics of this example will be presented and discussed. 5.1 S ystem D yn am ic E quations Consider a pendulum subjected to simultaneous horizontal and vertical planar excitations as depicted in Figure 5.1. The system possesses a viscous damping torque proportional to the angular velocity, Tviscous = cm ,i ■ 6 = m ■ L ■ c ■ 9. The governing differential equation of motion for the pendulum is cP 6 dO g . 1 d2X 1 d2Y . , . i ? + c * + i s m 9 + i d p cos# + i ^ s m ', = 0 (5 J) where m = mass of the pendulum bob, L = length of the pendulum, c = rotational damping coefficient, g = acceleration of gravity, 6 = angular displacement, X (t) = x-axis excitation, and Y (t ) = y-axis excitation. Consider a periodic support mo tion of the form X{t) = LXocosut (5.2) Y (t) = LYq c o s (ut + < f> ) 79 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Let r = ut, l Jq ~ £, c = /L tw # , then the equation of motion can be written as follows d?B u , d6 1 . . „ , , . / , / /k —— 4 -: — I — - sm 6 — X 0 cos 6 cos r — lb sin 0 cos(r + < p ) = 0 (5.3) dr 2 r dr H where r = ^ is the ratio of the forcing frequency u) to the natural frequency L o g of the pendulum. Define the state variables aq = 0, x 2 = jp = 0, and the parameter variable x 3 = r, and get the state representation of the system as follows. & i = x 2 (5.4) x 2 = —~ x 2 -----~ sin c tq + X0 cos x\ cos r -f Y 0 sin*i cos(r + < p ) x3 x\ The following values of system parameters were selected. H = 0.1, X 0 = 0.5, Y 0 = 0.5, x3 - 0.75 (5.5) £ Equation (5.4) was used for the cell center point integration process, < ^ > * (0 ) — » 4>*(T) i = 1 ,2 ,..., Nc, with T = 2tt. Using the procedure previously developed in Section 3.2.3 we can approximate £,-(r) by a vector homogeneous multinomial. Recalling that x (r) = ) + £ ,(t) with representing the solution to (5.4) 80 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . t Jr j * jp jr Jf J1 j\ Figure 5.1: Parametrically forced pendulum with viscous damping, and the fact that P t ( m ) ( r , £ ,-( t ) , s ) is nothing more than the Taylor series expansion of f(T ,x(r), s), we obtain for m = 4 M r ) 6,i(r) d h dx\ d h dxi £ l , t + d h dx 2 d h h,' (5.6) x=0*(T ) dx -2 1 d 3 f 2 3! dx\ & + 1 d2h ’ * ^ 21 dx 2 x = t f ( r ) 1 d*h 4! dx\ P ■ S l , j 81 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . where the following right hand terms of equation (5.6) are evaluated to be d h dxi x=0*(r) d h dx 2 x=0*(t) & h dx x=0i(T ) d h dxi x=^i(T ) (5.7) = 0 for j > 2 = -X o sin [<f>li{T) + t ] - ~ 2 c o s xi, d h dx 2 x = 4 > * (r) d3h dx\dx 2 x = < f> i (T ) d2/ 2 < 9 :r2 x = < f> i(T) d3/ 2 dx\ (r) dA h dx\ X =0?M X 3 — 0 for j > 2, k > 0, and I > 1 = -X o cos [ < f > * lti(r) + r] + sin [^.-(t)] J/O 1 = X0 sin [ # ti(r) + r] + cos [^.-(r)] /y » * itq sm f e w ] (5.8) Substitution of (5.7-5.8) into (5.6) results in a ra-jet approximation of system dynamic equations (see equation (3.4)) for m = 4 given by: 82 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. d . M = &,i (5.9) 6w(T) = {-*osin[<fi-(T) + T]-~cos[^yT)]){i|i— + I x3 J X3 | - ^ o cos [< f> *j(r) + r] + ^ sin [ ^ -( t)] | ^ + |X 0 sin \4>li(T) + r} + ^ cos [^ -( t)] | £h + ^ j x 0 cos [ ^ { t) + r] - ^ sin [<^(r)] | ^ The truncated point mapping algorithm is applied directly to (5.9) in order to calculate each mapping of points within the cell, i.e. the cell’s truncated point mapping. This completes the preliminary EPM problem definition, computational parameter specifications, etc. 5.2 G lobal A nalysis To demonstrate the EPM procedure the global behavior of (5.4) will be presented. For the EPM analysis the region of state space studied in R 2 was — 7 r < x\ < rr and — 7 r < x 2 < n. This system was previously analyzed using simple cell mapping (see Guttalu and Flashner [18]). 83 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. 5.2.1 Periodic Solutions and Basins of A ttraction We proceed with a fourth-order EPM analysis of the system by taking 51 x 51 inter vals in the phase plane (equivalent to 2,601 cells in cell state space S) and determine each individual cell’s truncated point mapping, G,(4)(^(0), s) i = 1 ,2 ,..., 2601. Equation (5.10) illustrates how the EPM method forms an analytical expression for an individual cell truncated point mapping with s = x 3 = r = 0.75. s) = 0.36 6 , + 0.75 6 , , - + 0.50 + 0.32 - (5.10) 0.01 d i + 0.10 d i - 0.34 d i t h - 0.28 - 0.29 & - 0.20 ^ - 0.58 Z i4 li - 0-87 ~ 0.57 - 0.27 £ 2 4 , G2 ,i(4 )(&, s) = -0.35 6 , + 0.48 + 0.66 + 1-08 h t f v + 0-91 + 0.51 d i + 0.81 6 ,^ 2 , + 1-28 + 0-42 f 2,i + 0.36 + 0-27 C i4li + 0-60 + 0-59 + ° - 0 9 & EPM trajectory evolution was performed for an initial condition resolution of 451 X 451. The results indicated there were three periodic orbits (near x\ ~ — 0.67T, x\ « — 0.27r, x\ & 0.27r) as shown in Figure 5.2 where the location of the P-K points found by the EPM method were indicated by the symbol “+ ” enclosed in small boxes. The periodic orbits were obtained by direct integration of (5.4) with the initial states being the location of the P-K points. Two of the orbits 84 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. were found to be asymptotically stable (indicated by solid curves) while the third is unstable (indicated by a dashed curve). In Section 5.2.2 we will demonstrate how the EPM P-K solution stability character was determined. Figure 5.3 shows the corresponding EPM domains of attraction. The initial conditions in the phase space were divided into three categories. Specifically, trajectories that map to the two asymptotically stable orbits with the darker region corresponding to the smaller Group A orbit and the lighter region to the larger Group C orbit and those that map to the sink cell (blank regions). In order to assess the accuracy of the EPM method an “exact” periodic solution response and corresponding domains of attraction was constructed by continuously integrating each initial condition until the response either went to a periodic so lution or left the region of state space S under study. The Direct Integration periodic solutions and basins of attraction are depicted in Figures 5.4-5.5 respec tively. First, we observe that the locations of the P-K points determined by the EPM method were remarkably accurate to reproduce the “exact” orbits without requiring additional numerical calculations. Table 5.1 documents a representa tive sample of EPM P-K points. Groups A and C were virtually identical to the “exact” results while Group C differs by no more than 0.5%. Furthermore, it is important to note a fundamental difference between the EPM and Direct Integra tion periodic solutions displayed in Figures 5.2 and 5.4 respectively. The EPM method accurately located the unstable solution while Direct Integration did not. 85 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . An additional numerical technique would be required to ascertain the existence of unstable solutions if one utilizes Direct Integration. Table 5.1: EPM P-K points vs. Direct Integration P-1 point. Percent Deviation from Group A P -l points Direct Integration 1 (0.56325,0.23582) (+0.01,+0.08) 2 (0.56335,0.23573) (+0.03,+0.04) Direct Integration P -1 point (0.56319,0.23564) Group B P-1 points 1 (-0.57128,2.0179) (+0.30,+0.030) 2 (-0.57243,2.0174) (+0.50,+0.005) Direct Integration P-l point (-0.5695,2.0173) Group C P-l points 1 (-1.8631,0.81491) (+0.04,-0.08) 2 (-1.8618,0.81707) (-0.03,+0.20) Direct Integration P -l point (-1.8624,0.81553) W ith regard to the EPM and Direct Integration basins of attraction for the two asymptotically stable solutions (Groups A and C), they match extremely well as depicted in Figures 5.3 and 5.5. Within the limits of the initial condition dis cretization no discernable differences exist between the two basins of attraction. In addition, the EPM technique determines the unstable manifold boundary be tween the two asymptotically stable solutions in the form of initial conditions that map to sink cell whereas Direct Integration does not. Again, additional alternative methods would be required if one was relying solely on Direct Integration. 86 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . Finally, Figures 5.6-5.7 display the results of applying a simple cell mapping analysis for the same level of state space discretization utilized for the EPM study, i.e. a cell array of dimension 51 x 51. For this case the simple cell mapping method did not detect the unstable solution and the two asymptotically stable P-K cell locations (Groups A and C) were grossly inaccurate. Moreover, only the most coarse features of their corresponding basins of attraction were uncovered with accompanying distortions in the basin boundary positions. Clearly, Figures 5.6- 5.7 illustrate the SCM state space discretization error introduced when utilizing a relatively large cell size hi for the system (5.4). It is worthwhile repeating the fact that EPM results for the 51 x 51 cellular space very closely approximate the “exact” results obtained by direct numerical integration. 5.2.2 Stability Analysis It is not difficult to demonstrate a general relationship between linearized flows and point mappings (see Guckenheimer and Holmes [13]) using Floquet theory. A brief review of some basic concepts and definitions reqarding Floquet theory for a system of differential equations with periodic coefficients is contained in Appendix B. The interested reader is referred to [5, 21] for an in-depth discussion of Floquet theory. For a nonlinear system x = f(f, x, s) the characteristic multipliers for the linearized flow are identical to the eigenvalues of the linearized point mapping, 87 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Expanded Point Mapping Method Cell TPM order is 4. Cell resolution is 51 by 51. Cell TPM integration stepsize is 2 i c / 100. Trajectory initial condition resolution is 451 by 451. 0.8 0.6 : Group B Outer Stable 0.4 Group C 0.2 Group A inner Stable -0.2 -0.4 - 0.6 - t-'Outer Unstable - 0.8 0.8 1 0 0.2 0.4 0.6 -0.4 - 0.2 - 0.8 - 0.6 ■ 1 X ,/ j u P-K solution locations Figure 5.2: EPM P-K Points. R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . X/w Expanded Point Mapping Method Cell TPM order is 4. Cell resolution is 51 by 51. Cell TPM integration stepsize is 2 « /100. Trajectory initial condition resolution is 451 by 451. -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 X,/lC Basins of attraction for the P-K solutions Figure 5.3: EPM Basins of Attraction. 89 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Direct Integration Method Initial condition resolution is 451 by 451. Cell initial integration stepsize is SOk ! 1000. 0.8 0.6 Group B Outer Stable 0.4 Group C 0.2 Group A -0.2 -0.4 -0.6 - •'(-’ 'Outer Unstable - 0.8 - 0.8 -0.4 0.2 0.4 0.6 ■ 1 - 0.6 - 0.2 0 0.8 1 X , / T C Periodic solution locations Figure 5.4: Direct Integration Periodic Solutions. R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. */ X Direct Integration Method Initial condition resolution is 451 by 451. Cell initial integration stepsize is 8 0 re /1000. -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Basins of attraction for the periodic solutions Figure 5.5: Direct Integration Basins of Attraction. 91 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Simple Cell Mapping Method Cell resolution is 51 by 51. Cell initial integration stepsize is 2% /100. 0.8 0.6 Group B Outer Stable 0.4 Group C 0.2 Vr Group A -0.2 -0.4 - 0.6 - -{-''Outer Unstable - 0.8 0.2 0.4 0.6 0.8 1 -0.8 - 0.6 -0.4 - 0.2 0 ■ 1 X,/7C P-K cell locations Figure 5.6: SCM P-K Cells. R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Simple Cell Mapping Method Cell resolution is 51 by 51. Cell initial integration stepsize is 2 it /100. I 1 i ' » ■ V ;V V ■ *.• V 1 V » « 1 > " • ------------- 1 - v v v i 1 i '■ i t 1 V W > * • * : • * * 1 ........................ ! a * a e V ; ' . ' ' o ’ 1 , • : . » X X 0 » .< a © a, x « -. X . © . a » * ;« A X a y a a A * A a X A / A a X A a a X A ’ V X A X A V X " X A W X ' . * X A v a 7 x a * ia ■ X a a a> x X A a a x a a * A a * X A « a! a a x a i a a x a A A ’ A A A a V v x ■ » •. X A a a x a a x a a a X a a a x a‘ a a a a' a a * a a a i • » • V V V • « a ® > V v v ; : © • " a" a y '» . . . . . ... . . . . . v v ' © r • i - a v v ; v v ©• ■a 7 “ ■ « © • a a © *• 0 • . « . i . < . © ^ A • • • « [ « » * c - X a 0 »; a a © I I * ;* A ! a • .• * " x > a K ? * » a x a a a > ; A a a x a - •. a x *• X A © 7 I 1 « ~ •, > * « ' a V i v a ' a ;. • A' © V V a ' a 'a ’x’ 'a' V « ' x a ‘ a ‘a- V v a 1 i ‘ < • ■ >: ' v ' x v ! v V X ■ • • > V .• • X A x a a a x a a a x a a »: A a a a | A a a a ; •. •. . ■ ; * * !* > * « d © a «; 0 a a « ; © a a 0 ; ■ • a a » « c » « e 0. A A » . a A 0 :• • *•; • S .-1 © - X • '«• •© F * t •© -jf •* •<:■©• i, * - i -x- . * © .• * « * a s © X A C 0 X A C © . • * A .• •. V X a a a a a a x A a a x a a a x X - v x ■ a ' a x x a ' . » V . x [x V X * a v x a a x a a a i * v vi x a a . * a v x ! - * A % « a .• •• V V f a a X A a a x a a a x a a a x o- a • ■ A © a * a a a > • a 7 ; > * « * > » c © > » * © x © © k « 4 1 t x; 0 a « a; 0 .• x « ! x a > x a a > * a a x a a a a a| a x w a « a ; « x x n a a ;» x 7 a S . A a a a A a a x A < a x ' «: a a T T > a I a x a a X A a e a a a. a a a a . a a a A . o » •* e a a fe a 0 a A e e* a c © a ' a a e a ' ; * - « o » .» « © > .< 1 C X > * © © » * < « © » A « * > *• 0 » » «: 0 © - 0 . . a * * * 7:7 a » a A 7 X a > * a * i“ X A s © X *; V / . t ; 7 7 * . • A V a \ •» A X a ’* a x ■ » a X A a a w x A a V X A V a- X A u x • i * A a • - a > a » f » ,> a a a > x a a if A a a x a a x a > a x a! a a a a . a a a a ] • % a A A a • * * -f * x * a X A a a w x a a a- a a a a - A a a A • - .! * © * .x » . * a . » . > . e . .a. -a p . a * . k . c .9. K. e 0 .4 .«• - ■ ei . a . « .a ,a .i.a . jf. © . . a i — . A k If o a • - « o » •» « c o a c 0 e © * a c a > a. ! * * a x c a m « - > © x Y « * a A © x a c © x ■ c © » « * © ' x a © ©| x a a - a » < e » » > » .« X C A X c © » X « © X A c © A a « a « X. 0 > * a . 0 a x ! ' g a » x ;« a a '* a A y a X A a © x A a a 4 a a X A V * ' X A w x ; _ . if. ? t .*. .< •-*. t .•i A .*. > . * . < . .a. .a / . ? * . < ( . < .». > . y . t 9 . \ .x. . .V > . * . . a — a « > a x ’ A a X A V a a a a a X A a a X a a a* e a a a ' a a a A ' ! * e a a * !■ - < o .» © a x e a e a a a c 0 a a c O X A X ©! a a 0 > a r 0 • S X » * •« « « : ■ © * < a a a c 0 a * e e a « X » a a- e a a a' a a a 0 • < •- a x « .'© a » « » i> o > .* a x x i © X S C © X A C 0 i A * o| * A O X A 0 * - ? * . .V ? ?. 5 .* . ' . .• f 5 .*». .*. .c . ? . ; - .« .» .- . • *t . ? .V .° . * . f . A A a iv a x a .a a v x f. a X A V a X A a X A V a j « a a «. x a v a i x -. v * ; • X X a v a a a X x -a x a v x a a X A a a x a a a x a a a x X - a a x a - a a .* ! ' * a A a ;* * A a x a a X A a a x a a a x a a a x a a a] n A a a ] - a a • » « a * * • •« a v f • • •> a » a f a a a a X A a a x © a a © 0- a « « > . a a « ' » X « o » a ’< a e o a c a e 0 a a c » a a‘ ■ » a * ' c » a 1 1 I I 1 1 1 1 I I i -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 X,/it Basins of attraction for the P-K cells Figure 5.7: SCM Basins of Attraction. R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. x(n-f 1) = G (x(n), s). Consequently the characteristic multipliers of the linearized system (5.4) i ( r) <9x x=4>Ht ) C M (5.11) will be identical to the eigenvalues of the linear part of the cell truncated point mapping (5.12) For this specific example problem the linearized flow of (5.4) is kib { ~Xo s in [# (f ( t ) + r] - ^ c o s ^ ^ r ) ] } U r ) (5-13) Applying Liouville’s Formula (see equation B.4) directly to (5.13) we obtain 2 r f J J A , = exp / i=i L Jo i* T — 2tr ds (5.14) 9 4 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Substituting fi and * 3 from (5.5) into (5.14) yields Ai • A 2 Ai • A 2 exp 2 ? r - 0.1 075 (5.15) 0.432679 to six significant digits. Since A * • A 2 is determined analytically, it is possible to verify the multiplication of the calculated eigenvalues. As explained in Section 3.2.4 the local P-K solution stability is determined by the eigenvalues of the matrix H = IljL i H?(s). For the P-1 solutions represented by Group A in Figure 5.2 j j ( G roup A) -0.165032 0.465830 -0.803010 -0.355164 (5.16) where the eigenvalues of H were A i j 2 = — 0.260098 ± 0.604176 j with j — V — 4- The corresponding magnitudes were |Ai| = |A2| = 0.657784, implying asymptotic stability. Note that Ai • A 2 = 0.432680, which is identical to (5.15) to five signif icant digits. Turning our attention to the P-1 points represented by Group C in Figure 5.2 we have 95 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. j|(G roup C) _ 0.855133 -0.822730 1.298395 -0.743216 (5.17) with eigenvalues Aij2 = 0.0559584± 0.655399 j, magnitudes |A x | = |A 2 1 = 0.657784, and A x • A 2 = 0.432680; implying asymptotic stability of Group C. Again ■ A 2 is identical to (5.15) to five significant digits. For Group B in Figure 5.2 we have j j( G r o u p B) _ 3.54822 -3.85470 1.17039 -1.14955 (5.18) with eigenvalues Ai = 2.20219, A 2 = 0.196477, and again Ai-A2 = 0.432680. Hence, the Group B P-1 solution in Figure 5.2 represents a saddle point. The white space in Figure 5.3 shows the location of points mapping to the sink cell (corresponding to the unstable manifold of the saddle point); this is in agreement with the result that the stable and unstable manifolds of a saddle point separate various domains of attraction. 5.2.3 Trajectory Evolution Accuracy Using the previous results from Section 5.2.1 where the cell resolution was 51 x 51 we generated a SCM and EPM trajectory from the same initial condition and 96 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. 3 1.5 *r o -1.5 -3 Figure 5.8: Trajectory Comparison. SCM versus Exact Integration, compared each trajectory with the “exact” integration solution. The results are shown in Figures 5.8-5.9. In this instance, the SCM trajectory diverged rapidly from the “exact” integration and ultimately went to the incorrect periodic response as shown in Figure 5.8. This illustrates the potential pitfalls of a coarse SCM cellular resolution and the numerical error introduced due to cell discretization effects. In contrast, the EPM response closely matches the “exact” response as depicted in Figure 5.9. W ithin the limits of the plot resolution the EPM method has reproduced a trajectory virtually identical to continuous integration. 97 Trajectory Comparison Simple Cell Mapping versus Exact Integration i................ i.............. -"I .— f i— i i ..........:..... i................"i.... ......i".............................. Exact + SCM J_______ I _______ I _______ L ___ .. I _______ I _______ 1 _______ I _______ I _______ I ----- ------J L 0 10 20 30 40 50 60 70 80 90 100 110 Time % R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Trajectory Comparison Expanded Point Mapping versus Exact integration n r~ I .........1 ..... ■ 1 .... Exact EPM + 1.5 o -1.5 -3 - !- I ! li ill A A ' " T f t H-1 N l-ltr vl!-1II phi j i n H c i r i it s f I I I I i ! a : f t , T i t H _ l______ I ______ I ______ L . J_ _ _ _ _ _ _ _ _[ _ 0 10 20 30 40 50 60 70 80 90 100 110 Time t Figure 5.9: Trajectory Comparison. EPM versus Exact Integration. R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. 5.2.4 Com putational Studies The computational time requirements for the EPM, SCM, and Direct Integration methods were also compared. All computations were performed on a Motorola StarMax 3000. In Section 5.2.1 we examined 451 x 451 initial conditions in the region of interest and determined both an “exact” basin of attraction and EPM basin of attraction. To make a meaningful comparison a SCM basin of attraction plot was constructed with a cellular resolution of 451 X 451 and each cell center point matching the 451 x 451 initial conditions used for both the “exact” and EPM basins of attraction. The cpu times are listed in Table 5.2. The cost for Direct Integration was approximately 78 times as high as needed by the EPM technique. Similarly, the SCM cpu cost was approximately 4 times greater than the EPM cpu cost. Clearly, the EPM method possesses a significant computational advantage over both the “exact” and SCM methods. Futhermore, it was interesting to note the SCM P-K cell locations were still not as accurate as the EPM P-K points. Again, this points out SCM cell discretization error, albeit small, exists even at fine cellular resolutions. Another advantage of the EPM method worth mentioning was the 78 fold reduction in the total number of cells required in comparison with SCM. In order to obtain the same level of EPM P-K solution accuracy it would be necessary to increase the SCM cellular resolution to 451 x 451 (equivalent to 203,401 cells in cell state space S). Moreover, it is worthwhile repeating the fact that EPM results 99 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Table 5.2: CPU time requirements for P-K solutions and basins of attraction. Method CPU Seconds Initial Condition Resolution Cellular Resolution EPM 3,649 451 x 451 51 x 51 SCM 14,193 451 x 451 451 x 451 Direct Integration 283,424 451 x 451 N/A for the 51 x 51 cellular space very closely approximate the “exact” results obtained by direct numerical integration. Finally, an additional EPM study was performed whereby the cellular resolution of the example problem was systematically lowered and the effects on both basins of attraction and P-K solution accuracy were examined. It was found for a fourth- order Expanded Point Mapping the cellular resolution could be reduced to 17 x 17 and still accurately locate the two stable and one unstable P-K solutions. The only inaccuracies introduced were in the boundaries separating the various basins of attraction. Thus, the EPM approach was extremely robust with respect to large cell sizes. 5.3 B ifurcation A nalysis We now proceed to vary the ratio of the excitation frequency u to the natural frequency l jq of the system and search for potential solution bifurcations along with their corresponding domains of attraction. For this analysis we will examine how the outer stable P-1 solution previously determined in Section 5.2.1 evolves as r is varied. 100 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. 5.3.1 M apping Preliminaries To perform a bifurcation study the parameter variable x 3 is formulated as x3 = r0 + r. We begin by determining the mapping about the perturbed outer stable solution, x = < f)* + z. x ;(t) = 0*(r) + Zi(r) i = 1 ,2 ,..., It (5.19) s = x 3 = r 0 + r where K is the order of the P-K solution. Note that represents the solution to (5.4). Again we use the procedure developed in Section 3.2.5 and approximate Z;(r) by a vector homogeneous multinomial of degree m in z, and degree q in r, namely z;(t) = Pi(m,9)(r,z<(r);f) i = 1,2,..., A A (5.20) The vector homogeneous multinomial p,-(m > ? )(r, z,(r); r) is the following Taylor se ries expansion of the right hand side of equation (5.4), i.e. f(r, x (r ),s). 101 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. * m ( t ) = [ M T i<l>KT Y , r o + r ) ~ f i { T , ( p * ( T ) ) r 0)\ + h , i { T ) = [ / 2 ( r > * ( 7' ) ; r o + r ) - / 2 ( r , ^ * ( r ) ; r 0 ) ] + d h dx2 d f2 x = < £ * M s = r 0 + r ^ 2 ,1 (5.21) d f2 dx2 1 d 4 f 2 , 1 d^f 2 x = * * ( r ) ^ ,; + 2! dx2 s = r o i- r l 4! dxf x = * ? ( t ) 2 l , i + s = r 0 + r + x =^?(t ) s = r 0+ f 1 dr/ 2 + 1 d3f2 x = * ? ( r ) 2 l , i + s = r o + r 3! a®? X=*f(T)Zl,t ^ s = r 0+ r 71 dx\ x = * ? ( t ) s= r o + r Evaluation of the individual terms in equation (5.21) results in ro + r) - /i(r , r0) d h dxi d f i dx2 d j h d x h { r , < A * (r ); ro + r) - / 2(r, </>*(r); r0) d h dx i ^ 2 dx2 djf 2 dx\dxl 2 0 0 0 for j > 2 — s m x i fJ-x 2 1 1 _(r0 + r)2 r( OJ 1 1 (5.22) (5.23) (r0 + r) r0 + f r0 - X 0 sin(xi + r) - V r0 + r 0 for j > 2, k > 0, and I > 1 102 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. d2h dx\ = - X 0 cos(o;i + t ) + 1 (r0 + r ): sin x\ v . , , , , 1 = Xq sm(xi + T ) + 7---- 7 - 3 7 7 c o s x i dx\ (r0 + r)2 Substitution of (5.22-5.23) into (5.21) leaves us with a vector homogeneous multinomial ps(m,g)('r, zt -(r); r) of degree m — 7 in z, • and degree c / in r. *2,i(r) .(^o + r)2 r0 j V&A'r) 1 1 cos r0 + r r0 f c ( r )] t ZU + (5.24) r0 + r X osin[^*.(r) + r] - ^ + -)2 1 ^ -X0 cos [< ^(t) + r] + * ,2 sin [ # fi(r)] [ > + ■ • • + Z2,i + 2! 1 7! (r0 + r)2 i ( x 0 sin [ t f / r ) + r] + ^ ^ _)2 cos [ # j t .(r)] \ 4 ,i In equation (5.24) the quantities g(r0 + r) h(rQ - + - r) 1 r0 + r 1 (r0 + f)2 (5.25) 103 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. were expanded separately to include terms up to order q = 10 in r. The expanded versions of g(r 0 •+ r) and h(rQ + f) were then substituted back into (5.24) resulting in a vector homogeneous multinomial Pt(m,9 )(7 5 z;(T)i r) °f degree rn = 7 in z4 and degree q — 10 in r. The TPM algorithm was then applied to (5.24) and we obtain the mapping of interest Zi(T) = G i M (Zi(0); f) i = 1,2,..., K (5.26) where T = 2n is the period of (5.4). 5.3.2 A nalytical M apping Formulae For r0 = 0.75 equation (5.27) represents the analytical mapping 1 G q7 ,io )(zd r) for the cell containing the outer stable P-l solution. G j -( 7 iio)(z,; r) is a function of the two state variables zi and z 2 and a single parameter variable r. Truncation for the state variables z was performed to order 7 while truncation with respect to the parameter variable r was of order 10. The functional relationship consisted of 395 terms for each component G1 ),(7il0)(zi; r) and (72^(7,io)(z,-; r) respectively. Equation (5.27) illustrates how the EPM method forms an analytical expression for this particular bifurcation analysis. Computations were performed using double precision. However, for purpose of illustration, the mapping was displayed to 2 decimal place accuracy. 104 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Gm(7,io)(z*-; f) = 0.80 ^ - 0.78 *2ii + 2.52 z\± - • ■ • (5.27) 0.28 zitiz2,i - 187 z f / 2 + 48 z ^ z ^ r 2 4----- 323 zu zl / 2 ~ 246 zfjZ2/r2 - 29 z \jr 2 + • • • ^ 2,1(7,10) (z* 5 7 " ) = 1-25 Ziti — 0.68 Z2ti + 2.02 z2^ — 0.12 z\^Z2ti — ■ • • 59zl-f2 + 273z1}iz3 2/ 2 + 893z \ ^ z \ f - • • • v m z \ - z 2^ 2 - 112 4 , -f2 + • • • The linear part of the analytical mapping in equation (5.27) for the outer stable P-1 cell location was H ,-(f) 0.80 + ... + 4.8E+08 f 1 0 -0.78 + ... + 4.6E+08 r 1 0 1.25 + ... - 5.9E+09 rw -0.68 + ... + 2.8E+09 f 1 0 (5.28) The bifurcation condition for a P-1 — » ■ P-2 solution can be expressed by the relation det[I + H,-(r)] = 0. For this cell location we have 3.3E+09 f 1 0 - 6.2E+08 ry + -------34.6 f + 1.55 = 0 (5.29) 1 0 5 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Equation (5.29) is an explicit formula which, is satisfied when a bifurcation to a P-2 solution occurs. To reiterate the EPM method forms an analytical mapping ex pression for each cell z , - i = 1,2,..., K under study. These functional relationships consist of both state and/or parameter variables. 5.3.3 Bifurcation A nalytical Results Again, we initiate the bifurcation analysis of the system starting with the outer stable solution 0*(O) = — 1.86240, ^(O ) = 0.815534, and r0 = 0.75. The Expanded Point Mapping study was carried out for these fixed values of system parameters. p = 0.1, Ah = 0.5, To = 0.5 (5.30) The bifurcation (frequency response) diagram in Figure 5.10 was created by solving the analytical expression - G $;,10)(zq r) = 0 * = 1,2,..., K (5.31) where K is the order of the P-K solution. As in the previous example from Chap ter 4, Matlab (Version 6.0 Release 12) was the solver of choice used to determine the solutions of equation (5.31). The parameter f was initially varied between 0 to 0.08 and z,- was determined from (5.31) and then substituted back into 106 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. X i = < t> * i + z i i = l , 2 , . . . , K (5.32) to determine the newly created P -l solutions. At f — 0.08 another mapping Gj(7 )io)(Zj-; r) was created by expanding about the solution previously calculated at r = 0.83. The bifurcation diagram in Figure 5.10 was constructed by varying r from 0.72 to 1.10. A period-doubling sequence of the P-1 solution takes place as r is increased. At r — 0.9284 the P-1 solution (Pi) bifurcated into a stable P-2 solution (Aj, A2) while the P-1 solution (Pi) itself became unstable. At r = 1.0580 the P-2 solu tion (Ai,A2) bifurcates into the P-4 solution (B \, P 2, B :i, BA ) which subsequently bifurcates into the P-8 solution (C i, ..., Cs) at r = 1.0785. The P-8 solution (C i, ..., Cs) bifurcates into a P-16 solution (Pi,..., Pie) at r = 1.0826. The exact locations of the various P-K solution into P-2K solution bifurcations were determined precisely by solving the relation det [I + H(r)] = 0 (5.33) where K H (r)= n - i [GifT.iojfziip] - n H r (s -34) 107 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Bifurcation / Frequency Response Diagram: p = 0.1 B4 B3 -0.25 B2 “V -0.5 A 1 -0.75 1.2 0.7 0.8 0.6 r Expanded Point Mapping Method Figure 5.10: Bifurcation-Frequency Response Solutions. 5.3.4 Bifurcation P-K Solutions and Basins of Attraction Once the bifurcation sequence was determined by employing the EPM method the next logical step was to perform a study of the system global behavior. Five different values of r were selected corresponding to a stable P-l solution, stable P-2 solution, stable P-4 solution, stable P-8 solution, and stable P -l6 solution. The stability character for the various P-K solutions was verified by examining the eigenvalues of the stability matrix H (r) using equation (5.34). Figures 5.11- 5.15 represent the continuous trajectories associated with the five different P-K solutions at their respective values of r, i.e. r = 0.88 (P-l), r — 1.00 (P-2), r = 1.075 (P-4), r - 1.080 (P-8), r = 1.083 (P-16). 108 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Stable P-1 trajectory: r= 0.88, ji = 0.1, XQ = 0.5 0.8 0.4 0.2 - 0.2 -0.4 - 0.6 - 0.8 0.6 0.8 0.2 0.4 - 0.6 -0.4 - 0.2 - 0.8 X,/Jt Figure 5.11: A stable P-l trajectory associated with the solution Pi for /i = X 0 = 0.5, and r = 0.88. R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Stable P-2 trajectory: r= 1, n = 0.1, XQ = 0.5 1 0.8 0.6 0.4 0.2 A2 0 - 0.2 -0.4 - 0.6 - 0.8 ■ 1 1 0.4 0.6 0.8 -0.4 - 0.2 0 0.2 - 0.8 - 0.6 ■ 1 X^/n Figure 5.12: A stable P-2 trajectory associated with the solution (Al5 A2) for (i = 0.1, Ao = 0.5, and r = 1.00. 110 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Stable P-4 trajectory: r= 1.075, |x = 0.1, XQ = 0.5 0.8 0.6 0.4 0.2 B2 B4 B3 - 0.2 -0.4 - 0.6 - 0.8 0.6 0.8 0.2 0.4 - 0.8 - 0.6 -0.4 - 0.2 X, /it Figure 5.13: A stable P-4 trajectory associated with the solution (Bi, B 2 , B 3 , for ji = 0.1, Xq = 0.5, and r = 1.075. R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. 1 0.8 0.6 0.4 0.2 - 0.2 -0.4 - 0.6 - 0.8 - 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 X ,/7 c Figure 5.14: A stable P-8 trajectory associated with the solution (C \,. .. , Cs) for /r = 0.1, X 0 = 0.5, and r = 1.080. Stable P-8 trajectory: r= 1.08, n = 0.1, X = 0.5 112 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Stable P-16 trajectory: r= 1.083, n = 0.1, XQ = 0.5 0.8 0.6 0.4 0.2 D3.D4 / , Djl,D2 - 0.2 D13,D14,D15,Dt6 D5,D6,D7,D8 -0.4 -O.l - 0.8 - 0.8 - 0.6 -0.4 - 0.2 0.4 0.6 0.8 0.2 XJ k Figure 5.15: A stable P-16 trajectory associated with the solution (D i,..., Di6) for y, — 0.1, X q = 0.5, and r = 1.083. 113 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Table 5.3 compares the P-K solutions determined from the relationships es tablished in equations (5.31-5.32) with the exact solutions derived by direct inte gration. It was important to note that the P-K solutions found using the EPM method were identical to the exact solutions by five to eight significant digits. Again this demonstrates both the accuracy and robustness of the EPM method over a wide range of values for the parameter r. Table 5.3: Bifurcation EPM P-K Solutions and Exact P -K Solutions. EPM P -l point at r = 0.89 Exact P -l point at r = 0.89 ~~~(-l .6188690,0.22992797) (-1.6188694,0.22992808) EPM P-2 points at r = 0.99 Exact P-2 points at r = 0.99 (-1.6117005,0.15185878) (-1.6117005,0.15185878) (-1.0119130,0.0089349646) (-1.0119131,0.0089349638) EPM P-4 points at r = 1.075 Exact P-4 points at r = 1.075 (-0.68496529,-0.079007936) (-0.68496529,-0.079007936) (-1.5948360,0.060010654) (-1.5948360,0.060010952) (-0.39016597,-0.15009216) (-0.39016389,-0.15009340) (-1.5321871,0.19355734)_______ (-1.5321882,0.19355602) Finally, the basins of attraction for the five P-K solutions were constructed using a cellular resolution of 51 x 51 and an initial condition resolution of 451 x 451 as was done in Section 5.2.1. Different grey scales were used for the domains of attraction of the P-2 solution, P-4 solution, P-8 solution, and P-16 solution to differentiate between trajectories mapping to the various P-K points for a given P-K solution. Blank regions correspond to trajectories that map to the sink cell. The basins of attraction provide both qualitative and quantitative information for 114 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. the complex manner in which the system global behavior changes as the parameter r varies. This bifurcation analysis demonstrates the versatility of the EPM method by using the alternative formulation from equations (5.31-5.32) to find P-K solutions as a function of the parameter r. The analytical mapping capability of the EPM methodology is a powerful feature that allows one to analyze a wide variety of nonlinear periodic systems exhibiting extremely complex behavior that would be difficult to analyze using traditional methods of analysis. Moreover, the validity of the EPM analytical results for large variation in parameter values has been clearly demonstrated in this example. 115 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Expanded Point Mapping Method Ceil TPM order is 4. Cel! resolution is 51 by 51. Ceil TPM integration stepsize is 2 it / 100. Trajectory initial condition resolution is 451 by 451. -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 hc Basins of attraction for the P-K solutions Figure 5.16: P-l EPM Basins of Attraction. 116 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Expanded Point Mapping Method Ceil TPM order is 4. Cell resolution is 51 by 51. Cell TPM integration stepsize is 2 it /100. Trajectory initial condition resolution is 451 by 451. -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 X, / it Basins of attraction for the P-K solutions Figure 5.17: P-2 EPM Basins of Attraction. 117 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. X J % 1 - 0,8 - 0.6 - 0,4 - 0 .2 - , 0 - -0.2 - -0.4 - - 0.6 - - 0.8 - -1 - Expanded Point Mapping Method Cell TPM order is 4. Cell resolution is 51 by 51. Cell TPM integration stepsize is 2 rt/100. Trajectory initial condition resolution is 451 by 451. ( 1 -------- 1 -------- 1 -------- 1 -------- 1 -------- 1 -------- 1 -------- 1 ----- ; —r < • ., ^ ~ * * s . ^ ^ ^ ,, , f < , ,v<\ L''v" ’ '/ s ' ," -'' 'Af-tv* ' m Ptiiist l i i i * _i_ _L. -0,8 -0.6 -0.4 -0.2 0 0.2 0.4 X ,/ K Basins of attraction for the P-K solutions 0.6 0.8 Figure 5.18: P-4 EPM Basins of Attraction. 118 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. Expanded Point Mapping Method Cell TPM order is 4. Ceil resolution is 51 by 51. Cell TPM integration stepsize is 2%! 100, Trajectory initial condition resolution is 451 by 451. -1 -0.8 -0 .6 -0.4 -0 .2 0 0.2 0.4 0.6 0.8 1 X^i % Basins of attraction for the P-K solutions Figure 5.19: P-8 EPM Basins of Attraction. 119 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. X n!n Expanded Point Mapping Method Cell TPM order is 4. Cell resolution is 51 by 51. Cell TPM integration stepsize is 2 i r / 100, Trajectory initial condition resolution is 451 by 451. -1 -0.8 -0 .6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 X , / t c Basins of attraction for the P-K solutions Figure 5.20: P-16 EPM Basins of Attraction. 120 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. C hapter 6 C oncluding R em arks In this dissertation a new numerical-analytical method called Expanded Point Mapping (EPM) was formulated. The Expanded Point Mapping approach com bines the cell to cell mapping and point mapping methods to form a comprehensive methodology for the study of the global behavior of nonlinear periodic systems and analysis of stability characteristics of equilibrium points and periodic solutions. The proposed method has the following features. • The EPM method has a well defined theoretical basis stemming from cell mapping and point mapping theories. The analytic expression for the map ping of points within a cell is exact up to any desired degree of approximation m. The only limitation on m is the amount of computational power required for the system under study. ® The EPM method is formulated for multi-dimensional systems. This feature is an important advantage over other analytic approaches for which analysis 121 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. of systems of order greater than two is very cumbersome, and sometimes even impossible. • Local stability studies of periodic solutions can be performed analytically on the m-th order analytic approximation of the point mapping within a given cell. • The trajectory evolution for EPM is based upon a continuum in state space. Hence, the EPM approach has the ability to generate any number of trajec tories. This feature allows an accurate global analysis of the system and is an advantage over cell mapping-based global analysis. • The EPM method imposes no explicit limitation on the number of consecu tive mappings. Therefore, it is well suited for analysis of systems exhibiting chaotic behavior. • The ability to obtain an analytical expression for the functional relationship between system parameters makes the EPM approach suitable for bifurcation analysis. For future research a number of areas of interest come to mind. The analysis of systems of order greater than two should be explored. Another possibility is the investigation of periodic and aperiodic solutions. The focus of this dissertation was finding and analyzing periodic solutions for a given system. However, the system may also possess solutions which exhibit either aperiodic or chaotic behavior. The 122 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. application of the EPM method to engineering problems in control systems should be considered as well. Extensive numerical studies have been performed in the area of helicopter rotor blade structural dynamics, e.g. Peters and Hohenemser [37], Friedmann and Tong [12], et al. A numerical method was employed to determine the state transition matrix for a periodic system of period T, i.e. 3>(T, 0). System stability charac teristics were evaluated using Floquet theory in conjunction with the numerically calculated state transition matrix. The numerical technique developed by Fried mann, Hammond, and Woo [11] was restricted to the analysis of linear periodic systems. The application of the EPM technique for global studies of nonlinear helicopter rotor blade dynamics is open for further study. Finally, an improve ment on the Expanded Point Mapping computational procedure with regard to the location of periodic solutions is in progress. 123 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. R eference List [1] J. Bernussou. Point Mapping Stability. Oxford, Pergamon, 1977. [2] N. N. Bogoliubov and Y. A. Mitropolsky. Asymptotic Methods in the Theory of Nonlinear Oscillations. Hindustan Publishing Company, Delhi, 1961. [3] V. V. Bolotin. The Dynamic Stability of Elastic Systems. Holden-Day, San Francisco, 1964. [4] L. Cesari. Asymptotic Behavior and Stability Problems in Ordinary Differen tial Equations. Academic Press, New York, 3rd edition, 1971. [5] E. A. Coddington and N. Levinson. Theory of Ordinary Differential Equa tions. McGraw-Hill, New York, 1955. [6] M. B. Dadfar and J. F. Geer. Power series solution to a simple pendulum with oscillating support. SIAM Journal of Applied Mathematics, 47:737-750, 1987. [7] J. R. Dormand. Numerical Methods for Differential Equations A Computa tional Approach. CRC Press, Boca Raton, 1996. [8] H. Flashner. A Point Mapping Study of Dynamical Systems. Ph.D. disserta tion, University of California, Berkeley, 1979. [9] H. Flashner and R. S. Guttalu. Analysis of nonlinear nonautonomous systems by truncated point mappings. International Journal of Nonlinear Mechanics, 24:327-344, 1989. [10] H. Flashner and C. S. Hsu. A study of nonlinear periodic systems via the point mapping method. International Journal for Numerical Methods in En gineering, 19:185-215, 1983. [11] P. Friedmann, C. E. Hammond, and T. Woo. Efficient numerical treatment of periodic systems with application to stability problems. International Journal for Numerical Methods in Engineering, 11:1117-1136, 1977. 124 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. [12] P. Friedmann and P. Tong. Nonlinear flap-lag dynamics of hingeless helicopter blades in hover and in forward flight. Journal of Sound and Vibration, 30:9- 31, 1973. [13] J. Guckenheimer and P. Holmes. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. Springer-Verlag, New York, 1983. [14] R. Guder, M. Dellnitz, and E. Kreuzer. An adaptive method for the ap proximation of the generalized cell mapping. Chaos, Solitons and Fractals, 8(4):525-534, 1997. [15] R. Guder and E. Kreuzer. Control of an adaptive refinement technique of generalized cell mapping by system dynamics. Nonlinear Dynamics, 20(1):21- 32, 1999. [16] R. S. Guttalu. On Point Mapping Methods for Studying Nonlinear Dynamical Systems. Ph.D. dissertation, University of California, Berkeley, 1981. [17] R. S. Guttalu and H. Flashner. Periodic solutions of nonlinear autonomous systems by approximate point mappings. Journal of Sound and Vibration, 129:291-311, 1989. [18] R. S. Guttalu and H. Flashner. A numerical study of the dynamics of a pendulum with a moving support. ASM E Journal of Applied Mechanics, 192:293-311, 1994. [19] R. S. Guttalu and H. Flashner. Analysis of bifurcation and stability of periodic systems using truncated point mappings. In Nonlinear Dynamics: New The oretical and Applied Results, pages 230-255. Akademie Verlag, Berlin, 1995. [20] J. K. Hale. Oscillations in Non-linear Systems. McGraw-Hill, New York, 1963. [21] P. Hartman. Ordinary Differential Equations. Wiley, New York, 1964. [22] C. Hayashi. Nonlinear Oscillations in Physical Systems. McGraw-Hill, New York, 1964. [23] P. Henrici. Discrete Variable Methods in Ordinary Differential Equations. Wiley, New York, 1962. [24] C. S. Hsu. Nonlinear behavior of multibody systems under impulsive para metric excitation. In K. Magnus, editor, Dynamics of Multibody Systems Symposium, Berlin, 1978. Springer-Verlag. [25] C. S. Hsu. A theory of cell-to-cell mapping dynamical systems. ASM E Journal of Applied Mechanics, 47:931-939, 1980. 125 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. [26] C. S. Hsu. Cell-to-Cell Mapping: A method of global analysis for nonlinear systems. Springer-Verlag, New York, 1987. [27] C. S. Hsu and H. M. Chiu. A cell mapping method for nonlinear deterministic and stochastic systems, Part I. The method of analysis. ASM E Journal of Applied Mechanics, 53:695-701, 1986. [28] C. S. Hsu and R. S. Guttalu. An unravelling algorithm for global analysis of dynamical systems: An application of cell-to-cell mappings. ASM E Journal of Applied Mechanics, 47:940-948, 1986. [29] E. I. Jury. Inners and Stability of Dynamic Systems. Wiley, New York, 1974. [30] H. Kauderer. Nichtlineare Mechanik. Springer-Verlag, Berlin, 1958. [31] A. J. Lichtenberg and M. A. Lieberman. Regular and Stochastic Motions. Springer-Verlag, New York, 1983. [32] J. E. Marsden and M. McCracken. The Hopf Bifurcation and Its Applications. Springer-Verlag, New York, 1976. [33] K. S. Miller. Linear Difference Equations. Benjamin, New York, 1968. [34] F. C. Moon. Chaotic Vibrations. John Wiley and Sons, New York, 1987. [35] A. H. Nayfeh. Perturbation Methods. Wiley, New York, 1973. [36] A. H. Nayfeh and D. T. Mook. Nonlinear Oscillations. John Wiley, New York, 1979. [37] D. A. Peters and K. H. Hohenemser. Application of the Floquet transition matrix to problems of lifting rotor stability. Journal of American Helicopter Society, 16:25-33, 1971. [38] F. M. Phelps, III and J. M. Hunter, Jr. An analytical solution of the inverted pendulum. American Journal of Physics, 33:285-295, 1965. [39] H. Poincare. Les Methodes Nouvelles de la Mecanique Celeste, 3 Volumes. Gauthier-Villar, Paris, 1899. [40] T. Poston and I. Stewart. Catastrophe Theory and its Applications. Pitman Publishing, London, 1978. [41] G. Ryland, II and L. Meirovitch. Stability boundaries of a swinging spring with oscillatory support. Journal of Sound and Vibration, 51(4):547-560, 1977. 126 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. [42] B. A. Schmidt. Pendulum with a rotational vibration. ASM E Journal of Applied Mechanics, 48:200-203, 1981. [43] R. A. Struble. Oscillations of a pendulum under parametric excitation. Quar terly of Applied Mathematics, 21:121-131, 1963a. [44] R. A. Struble. On the subharmonic oscillations of a pendulum. ASM E Journal of Applied Mechanics, 30:301-302, 1963b. [45] R. A. Struble. On the oscillations of a pendulum under parametric excitations. Quarterly of Applied Mathematics, 22:157-159, 1964. [46] J. M. T. Thompson and H. B. Stewart. Nonlinear Dynamics and Chaos. Wiley, Chichester, 1986. [47] B. H. Tongue. On obtaining global nonlinear system characteristics through interpolated cell mapping. Physica, 28D(3):401-408, 1987. [48] B. H. Tongue. A multiple map strategy for interpolated mapping. Interna tional Journal of Nonlinear Mechanics, 25(2/3):177-186, 1990. [49] B. H. Tongue and K. Gu. A higher order method of interpolated cell mapping. Journal of Sound and Vibration, 125(1):169-179, 1988. [50] B. H. Tongue and K. Gu. Interpolated cell mapping of dynamical systems. ASM E Journal of Applied Mechanics, 55(2):461-466, 1988. [51] B. H. Tongue and K. Gu. A theoretical basis for interpolated cell mapping. SIAM Journal of Applied Mathematics, 48(5):1206-1214, 1988. [52] M. Urabe. Non-linear Autonomous Oscillations. Academic Press, New York, 1967. [53] J. A. W. van der Spek, C. A. L. de Hoon, A. de Kraker, and D. H. van Campen. Parameter variation methods for cell mapping. Nonlinear Dynam ics, 7(3):273-284, 1995. [54] M. Vidyasagar. Nonlinear Systems Analysis. Prentice-Hall, Englewood Cliffs, New Jersey, 1978. [55] V. A. Yakubovich and V. M. Starzhinskii. Linear Differential Equations with Periodic Coefficients. Vols 1 and 2. John Wiley, New York, 1975. 127 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. A p p en d ix A E xpanded P oint M apping Source C ode The EPM algorithm was written in ANSI C style with the objective of producing a highly structured and portable program. There were 66 separate functions in total that comprised the EPM program1. W hat follows is a summary of the pseudocode representing the EPM program. The EPM program consists of three major parts. A .l D eclarations and Initializations 1. We first discretize the region of state space under study into a cellular structure. System order iV, the number of parameter variables, cell size hi z — 1,2,..., location of the cell state space origin z = 0, and the range of cell vectors z, i = 1,2,..., N c are user defined input. 1 The original source code was developed using THINK C/Symantec C + + Version 7.05 on a Motorola StarMax 3000 128 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. 2. Set up the cell center point integration and TPM procedures. The primary user defined input is (a) Total integration time T for both cell center point integration and cell TPM processes. (b) Order of Runge-Kutta method used for individual cell TPM computa tions. (c) Order of truncation to be performed individually on each state and parameter variable exponent. Order of truncation to be performed on the sum of the state variable exponents. Order of truncation to be performed on the sum of the parameter variable exponents. 3. The trajectory initial condition spacing in each state variable direction. A .2 P relim inary E P M C om putations 1. Each individual cell center point is integrated from tp = 0 to tp = T using a 7-stage FSAL fifth-order Runge-Kutta adaptive time step algorithm per Dormand [7]. 2. An individual cell TPM of order m in state variables and order r in parameter variables is calculated for all cells, namely G i(ro,r)(t-(0),s) * = 1,2 (A.l) 129 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. A .3 T rajectory C om putations A trajectory is computed for each initial condition analyzed using the EPM pro gram. Each trajectory will eventually evolve into either the sink cell, a P-K so lution, or an aperiodic motion. Two primary user defined inputs are required for P-K and aperiodic motions. 1. A P-K solution occurs when the condition |x(fc + * - ) - x ( * ) L , < e (A2) \x(k + K is satisfied. 2. A trajectory is considered to be aperiodic if k > k M ax occurs. A flowchart for the EPM program is given in Figure A.I. 130 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. No No No Is x(k+l) a P-K solution Is x(k+l) a chaotic solution Does x(k+l) map ■ to the sink cell , Yes Yes Yes No Yes Record the results Define and initialize the variables needed to analyze a region of state space for zj i = 1 ,2 Nc Cell Center Point Integration: <|>*j(T) Cell Truncated Point Mapping: Gi(- m r )(|j(0), s) For each initial condition, I .C.( j ) j = 1,2,... , Ni f compute the trajectory x(k) k = 0,1,2,. 1) Determine perturbation §j(k) 2) Apply cell TPM resulting in |;(k+l) 3) Determine x(k+l) = < j> * j(T ) + |j(k+l) ■Max Figure A.l: Flow chart for the EPM program. R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. A p p en d ix B F loquet T heory We consider the linear system described by x(i) = A{t) x(t) x(t) £ R n t € R + A(t + T) = A{t) (B.l) Theorem B .l Let he the fundamental solution matrix of (B .l). Then for all t, corresponding to every <&(t), there exists a nonsingular continuous matrix P (t) of period T and a constant matrix R such that < f> {t) = P(Y) etR P (t + T) = P(t) (B.2) 132 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. One can also show there exists a constant nonsingular matrix C such that C = eTR. The eigenvalues of C, pi, p 2, ■. ■, Pn, are called the characteristic (Flo- quet) multipliers. Similarly the eigenvalues of R, pi,p 2 , • • •, PN, are referred to as the characteristic exponents. The characteristic exponents are related to the characteristic multipliers by pk — eTpk k = 1,2,..., iV. Definition B .l A non-trivial solution (fft) of (B .l) is called a n o r m a l ( m u l t i p l i c a t i v e ) s o l u t i o n if there is a number pi such that < f > i(t + T) = pi < j > i(t) (B.3) Theorem B.2 (Stability of solutions of (B .l)) 1. The normal solutions of (B .l) are asymptotically stable if and only if the characteristic multipliers of A , \pk\ < I, or alternatively Re(pk) < 0 k = 1.2...., JV. 2. The normal solutions of (B .l) are bounded if and only if \pk\ < 1 k = 1.2...., N. This corresponds to stable normal solutions. 133 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission. 3. The normal solutions of (B .l) are bounded if and only if \p,k\ = 1 appears as a l x l matrix in the Jordan canonical form of C. In general there are no universal methods available for the exact determination of the characteristic multipliers pk- The Floquet theory determines only the form of the solution rather than the solution m atrix $(f) itself. Two useful theorems for our purposes are stated without proof. T h eo rem B .3 (L iouville’s Form ula) det 3>(t) = det 3>(to) exp | tr [A(s)] dsj- (B-4) If we assume is the state transition matrix, where 3? (to) = I? we can prove T h eo rem B .4 For the system x = A(t) x, where A(t) has the period T, let the characteristic multipliers of the system be ■ ■ ■ , Pn - Then p.i = exp | y tr [A(s)j d s j (B.5) with a repeated characteristic multiplier being counted according to its multiplicity. 134 R eproduced with perm ission of the copyright owner. Further reproduction prohibited without perm ission.

Asset Metadata

Creator Golat, Michael Charles (author)

Core Title A global study of nonlinear dynamical systems by a combined numerical-analytical approach

Contributor Digitized by ProQuest (provenance)

School Graduate School

Degree Doctor of Philosophy

Degree Program Mechanical Engineering

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag engineering, mechanical,OAI-PMH Harvest

Language English

Advisor Flashner, Henryk (committee chair), Dravinski, Marijan (committee member), Rosen, Gary (committee member)

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c16-268732

Unique identifier UC11339362

Identifier 3093955.pdf (filename),usctheses-c16-268732 (legacy record id)

Legacy Identifier 3093955.pdf

Dmrecord 268732

Document Type Dissertation

Rights Golat, Michael Charles

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the au...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus, Los Angeles, California 90089, USA