EIRATA TECPNICAL REPORT NO. 102 A STUDY OF FORCED OSCILLATIONS IN A CLASS OF nth-OReER NONLINEAR FEEDBACK SYSTEMS by A.. Nylor xix Change polynominal to read polynomial 121 Eq. (5.36) change SO(S) to read SO(s) 159 Change A(jc) to read A(Jcs) and A(-Jm) to read A(-j^,) STT^TJ 5BJ7 B(-, B3-Jco 173 Eq. (6.8) change p(t, 1) - to read p(t, 1) " 173 Eq. (6.10) change ai(l) - to read ai(l) m 188 Eq. (6.32) changi- +(t)3 O to read + c3(t)' 189 Eq. (6,3t) chmge 2f06 to read 2 and 1.03 to read 1 189 Eq. (6.35) change 2,06 to read 2 and 03 to read 1 189 Eq. (6.36) change 8.2 to read 8 00 in both places, 190 Eq. (6.h1) should read: t -(20 x 10 d6)o Pl 191 Eq. (6.42) shoald read: l o (3.75x lo0)l 1o0033.3(100e)- ) +j 33.3+.6( 1Ot )- a} = m5 -500o40(100o) + o10o + (20o)'(0oo0^) - 2~ P

Page 191 Eq. (6.h3) should read: p(t) + hpl(t) 4 co [l01.8t + 0.156n + {0.92 x 104osc [lOI.8t + 0.395nJ + 9.5 x 1040cos [305.4t + 0478n1} 191 Eq. (6.14) should read: p(t) + 7pl(t) a 7. cos [lOlt + 0.156 ] + ( {150 x 104cos [lOlt + 0.305n] + 64.5 x. 0Ioco [o303 - 0.566n} 192 Eq. (6.h5) should read: p(t) + I pl(t), 12.75 coo [lOOt] + L{-313 x 10cos [1oot - 0.64113 + 324 x lO cos [300t +.lOOn) 192 Eq. (6.46) should read: p2(t);19.o48 x 108cos [0ll.t + 0.268n] + 143.0 x 10-8cos.[30t - 0.6o6n] + 2.83 x 10cos [509.0t + 0*687n] 192 Eq. (6.47) should read: p2(t) a 22.3 x 1O4eco [0lOt + 0.865] + 1.,9 x 104co [s303t - 0.59 + 00o56 x oo10 [o 5t 0.9 ] 192 Eq. (6.48) should read: p (t)~ 1.76 x 104cos [ioot + o081i] 0.79 x 104cos [300t] + 0.84 x 10U cos [0oot + 0.l~7u]

Page 19h Point 3: should read: a10 + all J 100 + At1.079 a + a21 j100 + L2.7l1 a30 +'La31, -100 + AL4.dJ 200 hth line from bottom: change principle to read principal

THE UNIVERSITY OF MICHIGAN RESEARCH INSTITUTE ANN ARBOR A STUDY OF FORCED OSCILLATIONS IN A CLASS OF nth-ORDER NONLINEAR FEEDBACK SYSTEMS Technical Report No. 102 Electronic Defense Group Department of Electrical Engineering By: A. W. Naylor Approved by: A. B. Macnee Project 2899 TASK ORDER NO. EDG-1 CONTRACT NO. DA-36-039 sc-78283 SIGNAL CORPS, DEPARTMENT OF THE ARMY DEPARTMENT OF ARMY PROJECT NO. 3-99-04-106 Submitted in partial fulfillment of the requirements for the Degree of Doctor of Philosophy in the University of Michigan January 1960

TABLE OF CONTENTS Page LIST OF ILLUSTRATIONS v LIST OF SYMBOLS vii ACKNOWLEDGEMENTS xvii LIST OF APPENDICES xviii ABSTRACT xix CHAPTER I. INTRODUCTION 1.1 Statement of the Problem 1 1.2 Review of the Literature Pertinent to the Problem 7 1.2.1 Approximation Techniques 7 1.2.2 Results from the Theory of Perturbations 15 1,3 An Outline of the Development of the Technique Developed for Determining Periodic Responses 21 1.4 An Outline of the Technique Developed for Handling the Stability Question 25 CHAPTER II. DEVELOPMENT OF CANONICAL FORMS 2.1 Introduction 31 2.2 Normal Perturbation Problem 31 2.3 Inverse Perturbation Problem 36 2.4 First Canonical Form 42 2.5 Second Canonical Form 52 CHAPTER III. FIRST APPROXIMATION AND EXISTENCE CONDITION 3.1 Introduction 56 3.2 Determination of the Generating Solutions for the First Canonical Form 60 3.3 Determination of the Generating Solutions for the Second Canonical Form 70 3.4 Relations Between the Jacobian Condition and Behavior of Generating Solutions as Functions of the Forcing-Function Frequency 73 CHAPTER IV. HIGHER-ORDER APPROXIMATIONS 4 1 Introduction 78 4.2 Development of the Recursion Procedure 80 CHAPTER V. STABILITY CRITERIA 5.1 Introduction 93 5.2 Development of a Suitable Variational Equation 100 5.3 Transformation of the Variational Equation 106 ii

TABLE OF CONTENTS —Continued Page CHAPTER V —Continued 5.4 Determination of the First Terms, ao (i =,2,...,n), in the Power Series Expansion for the Characteristic Exponents, a. (i) (i = 1,2,...n) 108 5.5 Determination of the Coefficients, ai (i =,.n) of,u in the Power Series Expansion for the Characteristic Exponents, a,(u) (i = 1,2,...,n) 115 5.6 Determination of the Coefficients, a 2 (i =,.n) of u2 in the Power Series Expansion for the Characteristic Exponents, cai(O) (i = l,2,...,n) 127 5.7 Determination of the Coefficients, ai (i = 1,2,...,n; k = 3,4,...), of the Higher-Degree Terms in the Power Series Expansion for the Characteristic Exponents, aci() (i = 1,2,...,n) 139 CHAPTER VI. ALGORITHM FOR THE APPLICATION OF THE TECHNIQUES DEVELOPED IN THIS STUDY 6.1 Introduction 140 6.2 Selection of a Canonical Form 142 6.3 First Approximation to the Periodic Response: First Canonical Form 146 6.4 First Approximations to the Stability Criteria: First Canonical Form 149 6.5 First Approximation to the Periodic Response: Second Canonical Form 153 6.6 First Approximations to the Stability Criteria: Second Canonical Form 155 6.7 Higher-Order Approximations to the Periodic Response: First Canonical Form 156 6.8 Higher-Order Approximations to the Stability Criteria: First Canonical Form 162 6.9 Higher-Order Approximations to the Periodic Response: Second Canonical Form 166 6.10 Higher-Order Approximations to the Stability Criteria: Second Canonical Form 166 6.11 An Example of the Application of the Algorithm Presented in this Chapter 167 6.11.1 Selection of the Pertinent Canonical Form 169 6.11.2 First Approximation to the Periodic Response 173 6.11.3 First Approximations to the Stability Criteria 185 6.11.4 Higher-Order Approximations to the Periodic Response 188 6.11.5 Higher-Order Approximations to the Stability Criteria 193 CHAPTER VII. SUMMARY AND CONCLUSIONS 7.1 Introduction 198 7.2 Conclusions 198 7*3 Suggestions for Future Research 206 iii

TABLE OF CONTENTS —Continued Page APPENDICES 213 LIST OF REFERENCES 230 DISTRIBUTION LIST 232 iv

LIST OF ILLUSTRATIONS Figure Page 1.1 Block Diagram of the Nonlinear System Considered in this Study 4 1.2 Block Diagram of the Nonlinear System Considered by Tucker 8 2.1 Example of a Continuable Periodic Solution 34 2.2 Graphical Representation of the Allowable A Intervals for a "Normal Perturbation Problem" 35 2.3 Graphical Representation of the Allowable A Intervals for an "Inverse Perturbation Problem" 37 6.1 Location of Poles and Zeros of a Possible B() 144 6.2 Root Locus Associated with the B(z) Portrayed in Fig. 6.1 144 6.3 Location of ~JN and ~+jw in the Complex Plane 145 1 e 6.4 Example of a Possible Graphical Method for Solving Eq. (3.16) 148 6.5 Block Diagram for the Third-Order Nonlinear Feedback System Considered in this Example 167 6.6 Diagram for a Physical System Which Might Be Characterized by an Equation in the Form of Eq. (6.4) 168 6.7 Location of the Poles and Zeros of (z) 170 6.8 Root Locus Pattern Associated with Fig. 6.7 171 A(jw ) 6.9 Locus of s. Os 17 5 2 6.10 n, Zeroth-Order Describing Function, vs. a 177 6.11 Illustration of the Graphical Technique Employed in the Solution of the Zeroth-Order Describing Function Relation 179 6.12 Amplitude, "a," of the First Approximation to the Periodic Response vs. w: Locus of the Points at Which the Jacobian Expression Vanishes: Locus of the Points at Which Either or Both a11 and a12 Vanish 180 6.13 Phase Angle, e, of the First Approximation to the Periodic Response vs. as: E = 2 x 10 181 v

LIST OF ILLUSTRATIONS — Continued Figure Page 6.14 Phase Angle, 0, of the First approximation to the Periodic Response vs. w E = x 10 182 6.15 Phase Angle, 6, of the First Approximation to the Periodic Response vs. c E = 0.7 x 10 183 6.16 Phase Angle, 0, of the First Approximation to the Periodic Response vs. u E = 0.6 x 10- 184 vi

LIST OF SYMBOLS1 First Used Symbol Description on page A Amplitude of Kryloff-Bogoliuboff Approximation 25 A(z) nth degree Hurwitz polynomial in z with real coefficients 4 A Initial conditions for the solution [pj(t,A,)] 24 A(4) Initial conditions for the solution [cj(tl)] 56 B(z) mth degree polynomial in z with real coefficients (m < n) 4 C(z) C(z) = A(z) - zn 46 C Complex amplitude of Tucker's and Smirnova's approximations 9 C1 Right-hand side of Eq. (5.63) 131 C2 Right-hand side of Eq. (5.64) 132 C3 Right-hand side of Eq. (5.73) 135 C4 Right-hand side of Eq. (5.74) 135 D(t,L) nxn matrix, periodic t with period T, analytic in i, part of the solution of the variational equation 28 E Half the amplitude of Tucker's or Smirnova's forcing function 9 E Identity matrix 14 E Fourier coefficient in the Fourier series for the forcing function, e(t) 67 ^* s E Fourier coefficients for the terms with frequency 6 in the Fourier series for e(t) 67 F(x,4) n component column matrix 18 F(t,x,4) n component column matrix 16 F(x) Transfer characteristic of the nonlinear element, an entire function of x 4 1 For additional information regarding symbols see Table 1,1, p. 26* vii

LIST OF SYMBOLS — Continued First Used Symbol Description on Page k F Coefficient of in the power series expansion k of F{p(t, )} 86 FI Fourier coefficient in the Fourier series for 0 F 67 0 ~ F Fourier coefficient in the Fourier series for k F7 ~s s F Fourier coefficients for the terms with frequency 2 in the Fourier series for F 67 ~+s F Fourier coefficients for the terms with frequency - in the Fourier series for F 89 f. 0 F Matrix made up of -} (j,k = l...,n) 16 x oxk H Condition associated with Theorem I 17 I Initial conditions for the solution (p(t,I,k) or y(t, I, ) 24 I(A) Initial conditions for the periodic solution p(t,4) 61 Ii (u) One of n linearly independent initial condition vectors 29 J(A,4) Jacobian, Eq. (3.9) 66 K Positive constant 103 L nxn constant matrix 14 M(t),M(t,p) Fundamental solution matrices for variational equations 105 N Real constant used in the approximation of F(x) by N 44 N. Those values of N for which Eq. (2.5) has a pair of purely imaginary complex conjugate characteristic roots 44 P(t) nxn periodic matrix with period T and mean value zero 14 P(t,u) nxn periodic matrix with period T and analytic in p. for j11 sufficiently small 28 viii

LIST OF SYMBOLS — Continued First Used Symbol Description on Page &Q Constant term in the proposed generating solution for servomechanism systems 209 R(p) nxn matrix, analytic in 4 for sufficiently small 1'1. The characteristic roots of this matrix are the characteristic exponents of the variational equation. 28 dF So(t) dx for x = p(t) = p(t,O) 110 Sk(t) Coefficiee of 1' in the power series representation for x 110 x=p(t, y) S Coefficient in the Fourier series for S 67 0 O *2s W S Coefficient of the terms with frequency -in the ~0 Fourier series for So(t) 69 2 s s S Coefficient of the terms with frequency - in the Fourier series for S 69 S0 Mean value of S (t) 69 0 O S1 Mean value of Sl(t) 130 T Least period of (t), the periodic response of the nonlinear feedback system:'(t) = p(til) 3 T Least period of undisturbed oscillation 1 T Ti >4-dependent least periodic of oscillation, T = T(4) 19 T Least period of the periodic forcing function, e e(t) 1 T Least period of the dominant term in the generating S solution p(t). Te = (p/q) Ts, qTe = T pT = T. 2 Ui(t,4) n component column matrix, analytic in j for 141 sufficiently small, periodic in t with period T, forms the periodic part of y(t,Ij.,) 29 U n component column matrix, Eq. (1.25) 29 V1 Vandermonde matrix 49 ix

LIST OF SYMBOLS — Continued First Used Symbol Description on Page V Domain in (t,x)-space 16 X.t Y nxn diagonal matrix made up of J -elements 49 Z(tz) Higher order terms in the expansion of F(t,x) about d(t) 214 a Amplitude of the dominant term in the generating a _jE solution: c10 = e 74 a. (j = l,2,...,n) component of A 24 aj(G) (j = 1,2,...,n) component of A(4) 62 a Coefficient of. in the power series expansion of aj(p); i.e., aj( a) = a. + ajl + 4 aj2 +.. 81 a Part of y = ax - bx 1 b Imaginary part of Smirnova's n 12 0 b Part of y = ax- bx 2 blb2 Amplitudes for the sine and cosine term in the Van der Pol approximation 25 [c.(t,u)] Solution of the canonical form corresponding to the periodic solution, p(t,.) 56 cj(t, ) (j = 1,2,...,n) component of [cj(t,4)] 25 [cjo] Generating solution in the r -coordinate system 56 c. (j = 1,2,...,n) component of [cjo] 57 k Cj Coefficient of A in the power series for cj(t,) 25 C., {cJC. } 81 jk toyCjiT' jk lk' c2k (- < I < a) Fourier coefficients in the Fourier series representations for c (t) and c2k(t); if the superscript I is replace-kby "+s," -s etc., the coefficient of t -j Jfst -JoSt etc. is intended 89 x

LIST OF SYMBOLS — Continued First Used Symbol Description on Page c Complex coefficient in the generating solution; i.e., -j t jcu t -j00t p(t) = c1O0 + clO 24 clkc2k Mean values of ck(t), c 2(t) 89' A(jw3)Uo A(-c)p d(IclolOs) d(IclOIws) 5 {S - _ ~ } +2s -2s - S S, [see Eq. (3.22)] 69 e(t) Periodic forcing function with fundamental angular frequency we 4 f. (j = 1,2,..,n), components of the column matrix i^ F(t,x) or F(t,x,i) 17 g Real part of Smirnova's n 12 g1 Term in the first canonical form 51 gl Coefficient of i in the power series expansion for either gl or g2 86 g2 Term in the second canonical form 54 g Either gl or g2 from the first and second canonical forms 58 k 9g Coefficient of k in the power series.expansion for either gl or g 81 h(t;w ) Part of the generating solution for the second canonical form 53 i Integer (i = 1,2,...,n), designates a characteristic exponent, ca (i), and terms associated with xi (G) 29 j As a coefficient: j = 9 j As a subscript: (j = 1,...,n) specifies a component of a column matrix 24 k As a subscript, k = 0,1,2,..., specifies the coefficient of AJk in a power series 25 xi

LIST OF SYMBOLS — Continued First Used Symbol Description on Page k As an exponent; k = 0,1,2,...; e.g., k 25 (k) As a superscript: indicate kth-derivative 6 I An integer, -oo < I < co which as a superscript designates the Ith coefficient in a Fourier series 67 m A positive finite integer (m < n) 5 n A positive finite integer 5 n Zeroth-order describing function 67 n Smirnova's or Tucker's "describing function" 9 0 nk kth-order describing function 80 p Positive integer; e.g., U) = pcX 2 p(t) Generating solution of the "a = 0" differential equation 41 p(t) n component column matrix, periodic with period T, which is the periodic solution of a matrix differential equation 16 p(tpk) Continuable periodic solution of the i-dependent family of differential equations 56 p(t,4) n component column matrix, periodic in t with period T, which is the continuable periodic solution of a matrix differential equation 16 Pl(t),p2(t) Linearly independent solutions of Eq. (1.15) 20 pk(t) (k = 1,2,...) coefficient of }4 in the power series for p(t,l), pk(t) is periodic with period T 173 pk (-a. < A < oo; k = 1,2,...) coefficient in the Fourier series for Pk(t); ps is the coefficient of Pk k of jwst e~, etc. 188 q Positive integer 2 xii

LIST OF SYMBOLS-~- Continued First Used Symbol Description on Page q (Ap) One of n expressions from Eq. (3-3) involving A and A 62 C.(A,,) (j = 1,2) one of two expressions from Eq. (3.5) involving A and p 64 r Positive integer, 0 < r < k 86 r. (j = l,2,...,n), dependent variables of the canonical forms 24 r n component column matrix, [rl,r2,...,rn] 49 s A complex variable 5 s Variable of integration 58 t Independent variable: time in seconds 5 u.(t,p) (i = 1,2,...,n; j = l,..,n) components of the column matrix, Ui(t,k) 97 u. (t) Coefficient of k in the power series represen30 tation of u i(t,4) (i,j = 1,2,...,n; k = 1,2,...) 97 iJ9 ujk Coefficient in the Fourier series for u (t) 19 U j t ik u. Coefficient of c in the Fourier series for u. (t): -s, 2s, -2s, etc. relate to -jO', jaS,.jk S S -jw s, etc. 195 j2W t j2wst U121(t) U121(t) = U121(t) - "121(0)' 30 211(t) 211(t) = U211(t) U211(0)E 135 u(t) Periodic function of t 27 v. Element of the matrix V, j = 1,2,...,n (see AppenJf dix E) 51 d w w= — 46 x Dependent variable in Eq. (1.1): input to the nonlinear element 4 x (t) The approximation of Smirnova or Tucker 9 xiii

LIST OF SYMBOLS — Continued First Used Symbol Description on Page xt nt x = c E +.. + c E ( <r < k) 86 r r ir nr -- x. n component column matrix, [~tj~d((n-1)x 49 dt x n component column matrix, [xlX2,...,xn] 16 y Output of the nonlinear element 4 y n component column matrix; dependent variable in the variational equation 28 y(t, I) General solution of Eq. (1.20) 28 z z -., or a complex variable 4 0 ag -. t )(t:l) rI)(t4,) (j = 1,2,...,n) 107 j. (t) Coefficient of u in the power series expansion for r.(t,), (k = 0,1,...). r. (t) = B(z + xj)So(t) 3 o0 3 0 -A(j) (j = 1,2,...,n) [Eq. (5.20)] jk(t) = B(z + Xj)Sk(t) (j = 1,2,...,n; k = 1,2,3,...) [Eq. (5.21)3 109 A An arbitrary determinant 74 [A. (t)] Transformed ((t) 26 Q:) Summation point 4 ~(t) Periodic solution of Eq. (1.1) with period T 6 CI(t,c, ) n component colurmn rmatrix; a solution of Eq. (1.11) 17 ca IEPart of the initial condition for the solution V(t a, ) 17 a(M) Initial condition for which (t,a,u) is periodic in t wi h pericd T 18 ca General charactertstic exponent 27 xiv

LIST OF SYMBrOLS — Continued First Used Symbol Description on Page ca.(u) u-dependent characteristic exponent (i = l,...,n) 28 ik Coefficient of k in the power series for ai(u) 97'(1) (2) 1) c'1) The two possibilities for a1 117 3B Normalized frequency variable: = s 45 N. 7 y =' A(O), A = 0 66 5 Small, positive constant 6 c Small, positive constant 6 Base of the natural logarithm system: e = 2.71828... 6 a je o Argument of c: c 10 7 6 Phase angle of the Kryloff-Bogolioboff approximation 25 X Complex variable 14 X. (j = 1,2,...,n) characteristic roots of the "I = 0" 3d differential equation: k = jW, 2 = -Jio 24. ((j = 1,2,,..,n) characteristic roots of Eq. (2.6) J 3 and Xj 1,2,...,n) 45 p A parameter 16 GO Upper bound on \l\ 21 )1 ~ Sample value for the parameter. 35 7r 3.14159265... 3 P12 P2 Arbitrary complex constants 20 [p. ] Solution of a canonical form corresponding to 0 cp 1 (tI) (j =,2,...,n) 58 p. Component of [pjo] 64 [p.(t,A,pu)] A solution of the canonical satisfying the initial conditions [pj(0,A,)] = A 56 J~~~C

LIST OF SYMBOLS — Continued First Used Symbol Description on Page Auxiliary independent variable: t = T 46 cp (t,I) Periodic solution of the ", = 0" differential equation with initial conditions I 48 cp(tI,?) Solutions of the ui-dependent family of differential equations. The initial conditions of c(t,I,,) are denoted I. 24 e ) = 2c-, where W = 10 e e T e CP e 2i Ws Ws =T S =s =pT = 2Li N) e (C = 3-N for the applicable N.) 44 xvi

ACKNOWLEDGEMENTS The author wishes to thank the members of his committee for their suggestions and helpful criticisms. He particularly wishes to acknowledge the guidance and encouragement of Professors Wilfred Kaplan and Alan B. Macnee. Further, the author wishes to acknowledge the hospitality and assistance given him by the faculty of the Technische Hogeschool of Delft, Netherlands, where the initial phases of this research were carried out. Finally, the author wishes to thank Mr. Robert E. Graham, Mrs. Frances P. Holt, and Mrs. Marilyn J. Miriani for their preparation of this paper for publication. The major part of the research reported in this paper was supported by the U.S. Army Signal Corps. xvii

LIST OF APPENDICES Page APPENDIX A. Stability Criteria for Periodic Solutions of Ordinary Differential Equations: Variational Equations 213 APPENDIX B. Identities for the Jacobian Condition in Chapter III 216 APPENDIX C. Identities Involving Describing Functions 220 APPENDIX D. Commutation of the Operators and B(z) 224 APPENDIX E. Inversion of the Vandermonde 1Iatrix 227 xviii

ABSTRACT The purpose of this study is to develop a new method for employing perturbation techniques for investigating the stability of and approximating the periodic responses of a class of nth-order nonlinear feedback systems under the influence of a periodic forcing function. The ordinary differential equation which characterizes the class of the system considered is assumed to be of the forn: A(z)x-B(z) [F(x) + e(t)] = o (z where A(z) is an nth-degree Hurwitz polynominal in z with real coefficients, B(z) is an mth-degree polynomial in z with real coefficients (m < n), B(z)/A(z) can be identified with the transfer function of a physically realizable network, x is a real variable, F(s) is an entire function (s a complex variable), B(O) = 0, A(O) f 0, and e(t) is a bounded, continuous, periodic function of the independent variable t with least period T (T > 0). Nontrivial solutions, x = ~(t) (there e e may be more than one), periodic in t with least period T, T = q Te (q a positive integer), are sought and their stability investigated. If 4(t) is dominated by a term of least period T, T = T/p, (p a posis s tive integer) then the following designations have been used to describe the types of responses: (1) if p = q, then synchronization, lock-in, or frequency entrainment have been used; (2) if p = 1 and q > I, then subharmonic response has been used; (3) if p > 1 and q = 1, then ultraharmonic response has been used; (4) if p > l, q > I, and p p q, then ultrasubharmonic response has been used. T is related to T by a rational s e factor q/p; i.e., Ts = (q/p)T, and T = pTs = qTe. xix

The method developed in this study for finding the desired periodic solutions, ((t), is based upon a restatement of the problem, which allows certain ideas from perturbation theory to be applied. This restatement, called the "inverse-perturbation problem" in this study, involves the regrouping of the terms of the differential equation and the introduction of a parameter L. Two canonical forms exist for this restatement. The form which is applicable in a particular case is determined by whether or not the frequency spectrum of e(t) contains a term of fre1 quency -. The result of this restatement is an imbedding of the oris ginal differential equation into a family of differential equations which is continuously dependent on p so that a separate member of the family corresponds to each value of u. The member which corresponds to. = 1 is the original differential equation. The member which corresponds to M = 0 is a linear differential equation with constant coefficients. This linear differential equation is homogeneous or inhomogeneous, depending on which canonical form for the restatement is applicable. In either case, this linear differential equation has a two-parameter family of solutions periodic in t with period T. Selected members of this last family, p(t), act as generating solutions or first approximations for the periodic solutions, x = ~(t), of the original or u. = 1 differential equation. A method for accomplishing this selection is presented here. Further, an existence condition, based on the work of Coddington and Levinson, is developed which guarantees the existence of solutions, x = p(t,p.), for a family of differential equations periodic in t with period T, continuous in p. for arbitrarily small 1pj, and such that each lim p(t,p.) = p(t), a generating -0 solution. These p(t,u) are developed recursively as power series in p. xx

As an aid to this recursion procedure a set of describing functions and relations is defined. In the first canonical form the member o~ this set associated with the first approximation or generating solution is shown to be the usual describing function from control theory. Assuming that the radius of convergence of these power series in u includes = 1, the solutions x = p(t, 1) are of the desired type, x = 4(t), for the original or 1 = 1 differential equation. They are not, however, necessarily all solutions for the original differential equation of this type. This is one limitation of the method developed in this study. The stability question for the periodic solutions, x = p(t,) is answered by use of a method based upon some of the work of Moulton in celestial mechanics. For each p(t,p,) a system of variational equations is obtained which has coefficients dependent on p. as well as being periodic in t. Associated with each of these variational systems is a set of n characteristic exponents dependent on 1. The signs of the real parts of these characteristic exponents determine whether or not the solutions, x = p(t,.), are stable. These characteristic exponents, along with certain particular solutions of the variational equations, are developed recursively in a power series of,.* These recursive developments are accomplished through transforming, by means of n linear timedependent transformations, each variational system [i.e., one for each p(t,p.)] into n different systems of n first-order linear differential equations. Because of the choice of transformation, each new system has a solution periodic in t with period T. Each of these linear, timedependent transformations is also dependent on one of the unknown characteristic exponents which in turn depends on A; consequently, each transformation depends on p. Therefore, these transformations are not known at the outset, but are obtained recursively along with the characteristic xxi

exponents. The needed recursion procedures are developed by employing the periodicity conditions on the periodic solutions of the linear systems which result after transformation. Each linear system is suitable for determining only one characteristic exponent. Once the series for the characteristic exponents are known, the stability question is answered. The final part of this study consists of a presentation of an algorithm based upon the results obtained in this study and a working of an example. xxii

CHAPTER I INTRODUCTION 1.1 Statement of the Problem It is well known that if a periodic forcing function acts upon an oscillating nonlinear feedback system it may cause synchronous and resonant-like responses. Probably the most widely known example of such an interaction is the van der Pol oscillator under the influence of a sinusoidal forcing function [1]. This sinusoidal forcing function is assumed to have a least period Te, which is within some small neighborhood of the undisturbed oscillator's least period of oscillation, Te. If T is imagined to vary slowly through T, then response can be observed in the oscillator which is typical of those of interest in this study. In particular, as T comes sufficiently close to T, the period of oscillation will shift from T to T, and the amplitude of oscillation will O e exhibit a resonance-like maximum as T passes through T. The shift in the period of oscillation is known by various names: frequency entrainment, synchronization, and lock-in. The changes in amplitude of oscillation are usually referred to as nonlinear resonance. When both phenomena are considered simultaneously the term "harmonic response" has been used [2]. The preceding van der Pol oscillator case is a comparatively simple example of the systems of interest in this study. In general, both the nonlinear feedback system and the forcing function are considerably more complex. The order of nonlinear feedback system may be higher than two; that is, the differential equation which characterizes the system may be of order higher than two. In addition, the nonlinear element in the feedback loop mayr be described by a transfer charecteristic which I

2 is more complicated than the y = ax - bx3 of the van der Pol oscillator. Further, the forcing function need not be limited to a simple sinusoid. It may be a more general periodic function of time. The least period 1 may be near some rational multiple, p/q, of To instead of simply T i.e., T = (p/q)T, where T is within some small neighborhood of T e s s o As in the case of the simple van der Pol oscillator example, many of these more general systems exhibit some special response as T passes slowly through T. As before, it is possible that as T approaches o s sufficiently close to T the period of oscillation will shift from T to T and a resonance-like reaction will occur in the amplitude of oscillation. It is not expected that a forcing function with a least period T, associated with every rational multiple of T, will cause such phee s nomena. In many physical systems only those forcing functions associated with the simpler rational numbers cause significant response: for example, it is more likely that in a given system a significant response will arise from a forcing function with least period (1/2)Ts rather than one with least period (249/323)T5. Moreover, it may be that even some of the simpler rational numbers (e.g., 1/2, 1/6, etc.) will not have significant responses associated with them. Whether or not this happens in a particular system can be determined only by a detailed examination of that system. On the other hand, it is possible that a given system will exhibit a periodic response with a least period which has no obvious relation to T. Moreover, it may be that the undisturbed system does not execute a nontrivial periodic motion; therefore, it has no T associated 1 Throughout this study each rational factor is assumed to have a numerator and denominator with no common divisor other than one.

3 with it. In any event, the periodic motions executed will preserve the feature of having their least periods, T, be integer multiples of the least period, T, of the forcing function. As an example of a system in which this might occur, consider an oscillator in which the linear network has one resonant peak at 1 and another at 2. Assume that a, and 2 are not related by a rational factor. Further, assume that the entire system is such that the undisturbed oscillator executes a periodic motion with least period T1, near 2i/uy. It is then conceivable that as a forcing function with least period T2, near 2t/w2, is applied to the oscillator, a periodic motion with least period T2 will be executed. There may be even more abstruse responses, but of all possible responses only those that are periodic with a least period T equal to some integer multiple of T are considered in this study. If, in addition to being periodic with period T = qTe, the response is dominated by a term of least period T = T/p (p a positive s integer), then the rational number p/q provides a convenient method for classifying the response. The simplest case is 1/q = 1, for which the designation is, as in the above van der Pol oscillator example, harmonic response. If p/q = p, then the response is called ultraharmonic. If p/q = l/q, then the response is called subharmonic. And, if p/q is any other rational number, then the response is called ultra-subharmonic. The purpose of the present study is to develop a technique for determining periodic responses, such as discussed in the preceding paragraphs, and ascertaining whether or not the determined solutions are 1 A periodic solution, 4(t), is said to be dominated by a term of period Ts if, in the frequency spectrum associated with 4)(t), a term with frequency 1 is significantly larger than a combination Ths of all the other frequency components of the spectrum.

stable. This is done for a special class of nonlinear feedback systems and a special class of periodic forcing functions. A block diagram which characterizes this class of systems is snown in Fig. 11. e(t) + 18, B(z). A(z) F(x) FIG. 1.1 BLOCK DIAGRAM OF THE NONLINEAR SYSTEM CONSIDERED IN THIS STUDY The nonlinear ordinary differential equation which characterizes this class of systems mathematically is as follows: A(z)x - B(z) [ F(x) + e(t)] =0 ( (1,1) The terms in this equation have the following definitions and restrictions. i) x is a real variable.

5 ii) t is a real variable: time. iii) n is a positive integer. iv) m is a positive integer, m < n. v) A(z) is an nth-degree Hurwitz polynominal in z with real coefficients (coefficient of zn = 1). vi) B(z) is an mth degree polynomial in z with real coefficients (m < n). vii) B(z)/A(z) can be identified as a transfer function of a physically-realizable linear passive network. viii) F(s) is an entire function (s a complex variable). ix) B(O) = 0, A(0) o 0. x) e(t) is a real, continuous, sufficiently-differentiable, periodic function of t with least period T e xi) y is the output of the nonlinear element. It is readily seen that this differential equation is suitable for characterizing a large class of nonlinear feedback systems. The linear network is limited only by the requirements that it be physically realizable and have a transmission zero at z = 0 (considering z, for the moment, as the complex variable in a transfer function). This latter requirement guarantees that dynamic shifts of the operating point of the nonlinear element will not be transmitted from its output back through the linear network to its input. The methods which are developed in this study are, therefore, not immediately applicable to low-pass servomechanisms. However, it is shown in Chapter VII that this limitation can be rather easily removed. The only restriction on the transfer characteristic of the nonlinear element, F(x), is that it be an entire function when considered as a function of a complex variable. This is a

6 mathematical convenience which does not limit seriously the types of nonlinear physical elements which can be characterized. However, F(x) is a function of x only and of none of the time derivatives of x; therefore, it is not suitable for characterizing nonlinear elements which exhibit a hysteresis effect. The forcing function is suitable for almost all periodic forcing functions which are encountered physically. The solutions 4 (t) sought of Eq. (1.1) are, if they exist, nontrivial, periodic in t with least period T = qT, and often dominated by a term of least period T = (q/p)T. In addition to finding the solutions CD (t), it is also desired to know whether or not they are asymptotically stable. Tihoughciut this study asymptotic stability is defined as follows: A periodic solution 4(t) of Eq. (1.1) with period T which is defined for t > 0 is said to be asymptotically stable if, given any e > 0, there exists a 6 > 0 such that any solution x of Eq. (1.1) satisfying x (k(o) -C (0k ) < 5 [k = o, 1, 2..., (n-l)] where k k x() dx and - d dtk dtk satisfies ax(k)(t) - (k)(t) < e (t > O) [k,= 1, 2,..., (n-l)] and x(k)(t) (k) (t) (t co) [k = O, 1, 2,.., (n-l) A similar definition applies to differential equations in matrix form. The usual situation of interest in this study is the one in which the system has a stable periodic motion when e(t) = O. This periodic

7 motion is then disturbed as a nontrivial e(t) is applied. On the other hand, the technique developed here is not limited to this situation. The undisturbed nonlinear system could have, among other possibilities, x - O as its only asymptotically stable motion. The method is, however, limited to a consideration of only periodic responses to periodic forcing functions. Transient behavior is considered only in connection with the stability question for the periodic responses. 1.2 Review of the Literature Pertinent to the Problem The literature which pertains to the problem of interest in this study can be divided into two groups. The first group consists of papers in which various approximation techniques are utilized in the analysis of physical systems similar to the one considered in this study. The second group consists of papers which consider various aspects of perturbation theory. This latter group is important to this study because perturbation theory is the central core of the techniques which are developed in it. A discussion of each of these groups of publications is given below. 1.2.1 Approximation Techniques: Many authors have treated the problem of forced oscillations in nonlinear feedback systems. References to their work is given in the bibliography under items [31, [4], [5], [6], [7], and [8]. The authors from this first group whose work is the closest to the present problem are Tucker [3] and Smirnova [4]. Each of these authors considers a nonlinear feedback system under the influence of a sinusoidal forcing function. The forcing function is, therefore, less general than that considered in this study. Their problem is further limited in that they consider only harmonic response and

8 not the subharmonic, ultra-harmonic, or ultra-subharmonic response. The refore, the type of response, as well as the type of forcing function, considered by Tucker and Smirnova is less general than that considered in this study. The nonlinear feedback system investigated by Tucker is a slight variation of the one shown in Fig. 1.1. A block diagram which characterizes Tucker's system is shown in Fig. 1.2. B(z) A(z) G 1 F(x+e) + e(t FIG. 1.2 BLOCK DIAGRAM OF THE NONLINEAR SYSTEM CONSIDERED BY TUCKER A comparison of this block diagram with the one shown in Fig. 1.1 shows that, except for the injection point of e(t) and the limitation of the forcing function to a sinusoid, Tucker's system is the same as that considered here.

9 The technique which Tucker employs for finding, or rather approximating, the periodic responses of the above system is based upon a set of plausible assumptions which are typical of linearization techniques such as the describing function method. The key assumptions are that a jeZt * -jwt simple sinusoidal function, x (t) = Ce + C (C = complex conjugate of C), is a satisfactory approximation for the periodic response and that the linear network, which is mathematically characterized by B(z)/A(z), has a pass band which includes )s and stop bands which include all higher order harmonics of w, i.e,, 2w, 3sn..., etc. Granting these two assumptions, the output of the nonlinear element, y(t), becomes a function periodic in t with least period 2n/w; i.e., F(x + e) = jW) t -j( t. F (E + C)c + (E + C)E ]. [The forcing function, e(t), is +jW t -j t equal to EC s + Es (E real).] This periodic function is not usually a simple sinusoidal function, but has an associated Fourier series with a fundamental period equal to 2/wos; consequently, the output of the nonlinear element, y(t), has frequency components at integer multiples of X. Since these higher order harmonics fall within the assumed stop bands s of the linear network, they will not be transmitted through this network. The fundamental, or wa, component is, therefore, the only component that s must be considered. Tucker defines a "gain" for the nonlinear element as follows: The magnitude of the fundamental component of the output of the nonlinear element (1.2) n = ~..(1.2) The magnitude of the sinusoidal input to the nonlinear element This gain is a function of 2 jE + Cj, the magnitude of the sinusoidal input to the non ar mnt an an be thought of as efinin a linear element element, characterized by n0, which replaces the nonlinear element,

10 characterized by F(x+e); i.e., F(x+e) ->n x + n e 0 0 After the preceding replacement has been made, Tucker derives a system equation in a manner similar to that which would be used if the system were actually a linear feedback system. The resulting system equation is A(jo3 ) A(jos) - E ] B(j) A B(jwn ) S The solutions of this complex equation determine which values of C are admissible for a given amplitude, E, and angular frequency, ae (recall that a = w D ), of the forcing function. It is possible that more s q e than one C will satisfy the above system equation. For each allowed C jW ts -js t there is a periodic response x (t) = CE + C e which may or may not be asymptotically stable. In addition to deriving Eq. (1.3), Tucker presents a graphical technique for solving it which involves the separation of Eq. (1.3) into real and imaginary parts. He also discusses the stability question. Unfortunately, his discussion of the stability question is not complete. As he himself points out, it was not possible for him to determine all the necessary stability criteria. On the other hand, Smirnova claims to have been able, using essentially the same technique as Tucker, to obtain stability criteria. These criteria are only approximate, in the same sense as the solutions; nevertheless, they may be useful in many systems. Unfortunately, some of the steps in her argument are difficult to appreciate. 1 Neither Tucker's nor Smirnova's symbols are adhered to here.

11 The system considered by Smirnova is, except for the limitation of e(t) to a simple sinusoid, the same as that shown in Fig. 1.1. The only difference between it and Tucker's problem is, therefore, the injection point for e(t). The nonlinear ordinary differential equation which characterizes this system is A(z)x-B(z)[F(x)+2Esin t] =0 z d ) (1.4) e dt The symbols in this equation are, except for a real constant E, defined above for Eq. (1.1). The equation in Smirnova's method which plays the role played by Eq. (1.3) in Tucker's method is A(jow ) E n - - 0 (1.5) B(jw) o C The difference between Eqs. (1.3) and (1.5) is due to the different injection point for e(t). Smirnova solves this equation graphically to obtain allowed values for C (magnitude and angle). It should be remarked here that Smirnova's graphical technique appears to be different from Tucker's, although both techniques are basically the same. Smirnova solves Eq. (1.5) directly in the complex plane with vector quantities, whereas Tucker considers the real and imaginary parts of this complex equation separately. Another difference between the methods of Tucker and Smirnova is the way in which each defines n, the equivalent linear gain for the 1 Neither Tucker's nor Smirnova's symbols are adhered to here.

12 nonlinear element. By Tucker's definition, Eq. (1.2), n is real, whereas 0 Smirnova's definition is intended to allow a complex n. Since Smirnova 0 appears to believe that she has generalized n0, it is of interest to inspect her definition. Assuming that n = g + jb (g and b real), the linear replacement for the nonlinear element in Smirnova's work is characterized as follows: F(x )~ 2gICisin(w t + L^) + 2blClcos(w t + LC) (1.6) S 5 and x (t) = 21Cjsin() t + LC) (1.7) O s The significant aspect of Eq. (1.6) is that it allows the linear replacement for the nonlinear element to introduce a phase shift between its input and output. However, according to Smirnova's definition, b must be identically zero. This can be seen by considering her definition which is as follows: 2,/i C o b ~ fF [21Clsin(wt + LC)] cos(wt + L C) dt (1.8) e a e i gral hs a periic i ra w pri n The above integral has a periodic integrand with period- and an interval of integration that is also C-; therefore, this integral is indepens dent of the placement on the t-axis of the interval of integration. If the variable t is transformed by translation such that the sine function in the above equation becomes a cosine and the cosine becomes a sine, then the function F becomes an even function with respect to the new origin and is multiplied by a sine function; therefore, the new integrand is an odd function with respect to the new origin. Since the value of the integral is independent of the placement of the interval of integration, this interval can be centered about the new origin. The result is the

13 integral of an odd function over an interval centered about the origin: the value of this integral is zero. It is apparent from the preceding argument that Smirnova has not generalized n as she thought she had. Moreover, it is clear that any nonlinear element which is characterized by a function of x only always has a real n associated with it. This is true whether Tucker's 0 or Smirnova's definition is considered. A complex n0 is legitimate if the nonlinear element concerned is a function of the time derivatives of x as well as of x itself. An example would be a nonlinear element that exhibited a hysteresis effect. This is the point that Smirnova missed. The approximate method used by Smirnova to handle the stability problem is an extension of the equivalent linearization method which she and Tucker each use to obtain the approximate periodic solutions, x (t). A discussion of the rigorous approach to a stability problem of this nature is given in Appendix A. Contrasted with this rigorous approach, Smirnova's method seems to depend on two implied assumptions. The first is that a variational equation (see Appendix A) based upon x (t), the approximate solution, answers the "yes or no" part of the stability question correctly even though the values of the characteristic exponents may be shifted away from those of a variational equation based upon the actual periodic solution. Unless the actual periodic solution is known exactly, this is an assumption which must always be made. The second assumption is that satisfactory estimates for the characteristic exponents of a linear differential equation with periodic coefficients 1 It is also shown in Chapter III of this study that a complex n may arise in problems involving subharmonic, ultraharmonic, and ulrasubharmonic response.

14 can be obtained if the periodic coefficients are replaced by constant coefficients equal to their mean values. For example, consider the system, = Lx + P(t) x. (1.9) L is an nxn constant matrix. P(t) is a periodic matrix with mean value zero; i.e., each element of the matrix has mean value zero. This system has n characteristic exponents associated with it. The second assumption is that each root of the characteristic equation, det[L - E] = 0, is close to one of the characteristic exponents of Eq. (1.9). Above all, it is assumed that a characteristic root having a given sign for its real part implies the same sign for the real part of the neighboring characteristic exponent. In order to appreciate the danger attending this assumption, recall the Mathieu equation [9] which is a classic example of this assumption not being justified. Unfortunately, Smirnova's steps in the implementation of these assumptions are not clear. As has been pointed out previously, Tucker's and Smirnova's methods for finding the approximate periodic solutions are equivalent to the describing function method from control theory. The graphical method of Smirnova is exactly that which would be used if the describing function method were applied to the type of problem she considers. It is only with the stability questions that differences occur. In addition to Tucker and Smirnova, Adler [5] has investigated a problem which is closely related to the one of interest in this study. The nonlinear system that he considers is not so general as that of Tucker or Smirnova; however, he does attempt to handle responses which are not limited to being periodic. The technique which Adler develops is based upon instantaneous frequency and phase concepts. His basic

15 result is a differential equation describing the behavior of the instantaneous frequency of oscillation. Adler's method is not related in an obvious manner to those of Smirnova or Tucker. Nevertheless, for very simple systems, similar results should be obtained from all three methods. Adler's method does not supply information regarding the amplitude behavior of the oscillation. Huntoon and Weiss [6] have expanded Adler's method so that it predicts the behavior of the amplitude of oscillation as well as the frequency. Further, their method is suitable for treating more complicated oscillators than those treated by Adler. 1.2.2 Results from the Theory of Perturbations. As was stated previously, the second group of authors whose work pertains to the present problem is composed of writers who have discussed aspects of the Theory of Perturbations. These authors do not, in contradistinction to the foregoing group, consider directly the physical problem of interest here or even a closely related one. They are concerned with the mathematical tools of perturbation theory. It is just these tools, however, that are used in the creation of the analysis technique presented in this study. It will be shown in detail in Chapter II that the original problem can be restated as a problem which can be treated with the techniques from perturbation theory. The results given below from the work of these authors are not intended to be a general survey. Only those results which are needed later in this study are selected and presented. The idea behind perturbation theory is the assumption that as a differential equation is continuously transformed into a neighboring differential equation a periodic solution of the first differential

16 equation is continuously transformed into a periodic solution of the second, or neighboring, differential equation. The application of perturbation theory can be considered as a process of continuation of a known periodic solution to an unknown periodic solution. The type of problem from perturbation theory to be considered here is as follows. An ordinary differential equation that depends on a parameter, say A, is assumed given. Let it have the form, I = F(t, x, 4), (1.11) dt in which d and F(t, x, p) are n component column matrices. As x appears dt0 in (t, x, p), it is an abbreviation for xl, x2, x3,..., x; where the xl, x2,..., and x are the components of the n-component column matrix x. The column matrix F(t, x,.) is periodic in t with period T (not necessarily a least period) and satisfies continuity and Lipschitz conditions. It is assumed that a periodic solution, p(t), with period T (also not necessarily a least period) is known for this differential equation for some value of., say p = 0. The problem is then to determine whether or not a solution p(t, p.) exists for sufficiently small 1ip1 which is continuous in p, periodic in t with period T, and such that lim p(t, p) = p(t). Further, it is useful to know the conditions under.110 which p(t,.) may be expanded in a convergent power series in p. An important theorem by Poincare regarding the preceding problem in perturbation theory is as follows: Theorem I. If F and F are real and continuous in (t, x, I) when'(t, x) is in some domain V of (t, x) space containing the curve 1 Unless stated otherwise, all periodic solutions are nontrivial in this study.

17 [t, p(t)j, and when \tl is small, and if the first variation of Eq. (1.11) at A = 0 with respect to the solution p(t) has no solution of period T, then for small 141 the Eq. (1.11) has a solution p(t, p.), periodic in t with period T, continuous in (t,.), and with p(t, o) = p(t). There is only one such solution for each.. f. F is the nxn matrix made up of f (j, k = l,..., n) where the f. are X Xk the components of the matrix F. Let the condition that the first variation of Eq. (1.11) with respect to p(t) has no solution periodic in t with period T be denoted by H. If the conditions of the above theorem are satisfied, then it is at least justifiable to consider periodic solutions, p(t, p.), for sufficiently small 1I. In this case, p(t, p) is the continuation of p(t) for j1l sufficiently small. Note that this theorem says nothing about the maximum range of p over which p(t, p) exists. A complete proof of this theorem is given on page 349 of [10]; however, since this proof follows a line of reasoning which is typical of theorems in perturbation theory, it is outlined here. The solution of Eq. (1.11) which assumes the initial value p(o) + a at t = 0 is denoted. = C(t, a, lp). It follows from the uniqueness of i( and the periodicity of F that for this solution to be periodic of period T it is necessary and sufficient that CD(T, a, p.) = C(o, a,.) or that Z(T, a, p.) - p(o) - a = 0. This last system of equations is an implicit relation between the components of a and.. For p = O, it has a solution a = 0. In order to see this, recall that p(t) is a periodic solution for p = 0 and that D(T, o, o) = p(o). If the Jacobian of 1 See Appendix A for a discussion of first variations.

this system taken with respect to the components of ac is nonvanishing at pA = 0, ca = 0, then it follows from the implicit function theorem that the system has a unique solution, ca = ca(), in the neighborhood of p = 0, a = 0, and a(p) is continuous in, with a(0) = O. The remainder of the proof shows, by some manipulation of the equations involved, that the nonvanishing Jacobian condition is equivalent to the condition that the first variation with respect to p(t) have no periodic solution with period T. The key idea in this proof is the employment of a periodicity condition in conjunction with the implicit function theorem. Note that the Jacobian which arises from the application of the implicit function theorem depends only upon the known solution p(t). Unfortunately, in many important cases the condition H is not satisfied. A classic example is an autonomous system of differential equations. Such a system would characterize an undisturbed oscillator. Consider the following system of differential equations: dx = F(x, A) (1.12) The symbols in the above equation are similar to those in Eq. (1.11). The first variation of this equation with respect to a known periodic solution p(t), where p(t) is defined as before, is I F x[(t), oj y. (1.13) dt x It is easily shown that dp/dt is a periodic solution with period T of the variational equation, Eq. (1.13). This immediately demonstrates that the condition H cannot be satisfied by systems of this type; therefore, Theorem I is not applicable to this important class of differential equations.

i9 It is shown on page 352 of [10] that if there exists only one periodic solution of the variational equation, then there exists for \41 small a unique solution p(t, p) of Eq. (1.12) continuous in A and periodic in t with period T, where T = T(), a continuous function of A. It is, on second thought, not surprising that the period of oscillation of an autonomous system should depend on A. Another example of an important system that does not satisfy condition H is the following: x= Lx + A F(t, x, A). (1.14) dt dx L is an nxn constant matrix; x, d, and F are n-component column matrices; the components of F satisfy the same conditions as the F in Eq. (1.11); x in (t, x, A) is shorthand for x1, x2,..., x; and A is a real parameter. If. = 0, then Eq. (1.14) reduces to the linear system, I Lx. (115) dt It is assumed, as before, that this reduced, or p = O0 differential equation has a periodic solution p(t), with period T. Since the first variation of Eq. (1.14) with respect to p(t) is the same as Eq. (1.15), p(t) must also be a periodic solution with period T of the variational equation. Therefore, the condition H is not satisfied, and Theorem I does not apply. Therefore, the difficulty which arises in the consideration of systems which are characterized by Eq. (1.14) is the determination of which solutions, p(t), if any, can be continued. As before, the word "continued" is intended to mean here the relation between p(t) and p(t,p),

20 The problem is that the A = 0 differential equation is linear; therefore, it has a family of periodic solutions. For example, if Pl(t) and p2(t) are two linearly independent periodic solutions with period T (not necessarily a least period) of the A = 0 differential equation, Eq. (1.15), then plpl(t) + p2P2(t) (P1 and P2 arbitrary constants) is also a periodic solution. If the linear differential equation has k linearly-independent periodic solutions with period T (not necessarily a least period), then it is said to have a k-parameter family of periodic solutions. It should not be surprising that in most cases only selected members of this kparameter family of solutions may be continued for A ~ 0. The linear system, Eq. (1.15), is the limit, as L - 0, of a nonlinear system, Eq. (1.14). The limit cycles associated with each p(t, A') are expected to approach limit cycles of the A = 0 differential equation as A - 0. It is not usually the case for nonlinear, 4 f 0 differential equations that these limit cycles of the p = 0 differential equation which are approached will comprise the entire k-parameter family of solutions; however, such a situation is not impossible. Systems of the type discussed in the preceding two paragraphs are a form of degenerate system. The order of their degeneracy is said to be equal to the number, k, of linearly independent solutions with period T of the variational equation. For a general discussion of degenerate systems see [11], [121, [131, and [14]. The particular degenerate problem of interest in this study, i.e., Eq. (1.14), has been treated by Coddington and Levinson [10]. They present a method for selecting the members of the k-parameter family of periodic solutions of the p. = 0 differential equation which can be continued. Further, they show that if a certain Jacobian is nonvanishing, then the existence and uniqueness of the continuations, p(t, A), are guaranteed for p. suffi

21 ciently small. The essential features of the Coddington and Levinson method are, as in the proof of Theorem I, the utilization of a periodicity condition to obtain an implicit relation between the initial conditions and u and the application of the implicit function theorem. Coddington and Levinson show also that if F is analytic in x and A, then the desired periodic solution, p(t, p), can be obtained recursively as a power series in i. Lewis [15] has considered a problem which includes the problem considered by Coddington and Levinson. The important difference between Lewis's results and those of Coddington and Levinson is that Lewis is able to estimate a range of p over which the solution can be continued, whereas Coddington and Levinson develop a solution for which the existence and uniqueness is guaranteed only for l1 sufficiently small. Unfortunately, Lewis's method depends upon making an a priori estimate of the maximum value of x to be expected. This is a disadvantage, because the guaranteed range of 4 may depend on this initial estimate. Moulton [16] has considered another problem in perturbation theory which pertains to the stability of the periodic responses of the nonlinear feedback systems considered in this study. The problem considered by Moulton is the finding of the solutions and characteristic exponents of a system of n first-order linear differential equations that have coefficients periodic in t with period T and analytic in pi for 0 < p < O. The importance of this problem to the present study is discussed in Section 4 of this chapter. 1.3 An Outline of the Development of the Technique Developed for Determining Periodic Responses The development of the method for determining the periodic responses of a class of nonlinear feedback systems to periodic forcing func

22 tions proceeds in three steps. Each of these steps is accomplished in one of the three following chapters. The first, and crucial, step is a restatement of the differential equation which characterizes the nonlinear feedback system under the influence of a periodic forcing function. This restatement results in a new differential equation amenable to treatment by the techniques of perturbation theory. This step is accomplished in Chapter II of this study. The restatement has two parts: (1) a regrouping of the terms of the original differential equation, and (2) the introduction of a parameter i. If, after restatement, p is set equal to zero, a linear system with constant coefficients is obtained which has two linearly-independent periodic solutions of period T (not necessarily a least period) and, as before, T = qT. If, also after restatement, p is set equal to one, the orige inal differential equation is obtained. The regrouping is accomplished by dividing the original differential equation into two sections. The first section contains all those terms which comprise the "p = 0" linear differential equation. The second section contains those terms which make up the difference between the original differential equation and the "L = 0" linear differential equation. This latter section contains all the nonlinear terms. The parameter p is introduced as a multiplying factor for this second section. The above restatement of the original differential equation can be thought of as an imbedding of the original differential equation into a family of differential equations. A separate member of this family corresponds to each value of,. The member of this family that corresponds to p. = i is the original differential equation. The member which corresponds to p = 0 is the aforementioned linear differential equation

23 with constant coefficients. Recalling the discussion in Section 2 of this chapter, it can be seen that this u-dependent family of differential equations is related to the type of problem treated by Coddington and Levinson. The t" = 0" differential equation has a known two-parameter family of periodic solutions with period T (not necessarily a least period). The group of terms multiplied by p is either independent of t or periodic in t with period T (T = qT ). The group is independent of t if the forcing function, e(t), is not a part of it. If the forcing function is a part of it, then the group is explicitly and periodically dependent on t. The period is, of course, T. Concomitantly,the ". = 0" linear differential equation is homogeneous if e(t) is multiplied by A, and inhomogeneous if it is not. The existence of the possibility of these two locations for e(t) motivates the use of two canonical forms in the restatement of the original differential equation. The choice between these two canonical forms for a given problem depends upon the frequency spectrum associated with e(t). If this spectrum contains a term of the same frequency as that of the periodic solution, p(t), for the "I = 0" linear differential equation, then the canonical form is chosen in which the forcing function is multiplied by p.. The reason for such a choice is simply to avoid the resonance response that would occur in such a inhomogeneous system, i.e., t sin t. This type of response would destroy the periodicity of the solutions of the "}J = 0" differential equation. The reason for choosing the other canonical form is more subtle, because it is related to the satisfaction of the nonvanishing Jacobian condition. The choice of canonical forms is discussed in Chapter III. In addition to the above restatement, a transformation of variables is carried out in Chapter II. The transformation used is similar

24 to that employed in the variation-of-constants technique. The result is the transformation of an nth-order differential equation obtained after restatement into a system of n first-order differential equations of the following form: dr. rl = f (rl, r2,''" rn) t) r2 = 4f2 (ri, r2''*0" rn t) rn = nf n (r1, r2,.. rn' t) The form of the solution to the original, untransformed differential equation with initial conditions I (I I1, I2,..., I ) is denoted, x = P(t, I, A), and the corresponding solutions of Eq. (1.16) are denoted, r. = p.(t, A ), (j = 1, 2,..., n), where lt X2t n t cp(t, I, ) = pl(t, A, L)e + p2(t, A, p)e +... + pn(t, A, A)e n (1.17) and A> (al, a2,..., an) is the transformation of I and the A's are the characteristic roots associated with the "p = 0" linear differential equation. This form of the solution is shown later to have a corresponding generating solution or first approximation of the form, jwt -jw t p(t) = c0 s + 0 e, (1.18) where c10 is a constant (c10 is the complex conjugate of c10) and(rl = c10, r2 = c10, r3 = 0,..., r = O) is a solution of Eq. (1.16) at p = 0.

This complex form for the generating solution is to be compared with that of van der Pol, p(t) = b1 cos wst + b2 sin w t, and that of Kryloff and Bogoliuboff, p(t) = A cos (a t + Q). In Chapter III, the method of Coddington and Levinson is adapted so that it can be applied to Eq. (1.16). The results of this chapter are the development of a graphical technique for the finding of the generating solutions (i.e., c10) and determining whether the Jacobian condition is satisfied or not. In Chapter IV, a recursion technique is presented for obtaining the pj[t, A(g), p.1 = cj(t, p.)'s [A(u) are those initial conditions which correspond to the selected solutions of Eq. (1.16)] which correspond to periodic solutions, p(t, p.), of the untransformed, u-dependent family of differential equations as a power series in'; i.e., oo c.(t,.) = cjk(t) 11 (j = 1, 2...., n) (1.19) k= 0 The actual periodic solution is obtained by substituting the above expression into Eq. (1.16). As part of this recursion technique, additional functions and relations are designated which are similar in nature to the describing function. It is, in fact, shown that the member of this group associated with the determination of the c. (t)'s, the generating Jo solutions, is just the ordinary describing function. Table 1.1 lists, defines, and relates some of the more important differential equation solutions and variables that appear in this study. 1.4 An Outline of the Technique Developed for Handling the Stability Question The general approach to the question of asymptotic stability for a periodic solution of a nonlinear differential equation is presented

26 TABLE 1.1 All elements of the left-hand column refer to the variable x, and all elements in the right-hand column refer to the rj variables: both columns are related by the coordinate transformations which relate the r.s to x. J Dependent variable in Eq. (1.1), and [rj] An n-component column ma input to the nonlinear element (also, which results from a tra: x |used as a column matrix in other (j = 1, 2,..., n) formation of x and its d. parts of this study) vatives. Dependent varis in both canonical forms. Periodic solutions of the ". = O" Solution of the canonica. diff. eq. with period T and initial [ Pjo3 form at 0 = correspond p(t I) | conditions I. Includes all such= 1 2 n) to Pp(t, ) solutions for the " = 0" diff. eq. Continuable (t, I); i.e., |jc Solution of canonical fo: p(t) p(t) = qp(t (O)) = p(t, O), (j, 2, n) corresponding to p(t) Generating solutions Solution of the I-dependent family [p.(t, A, 2)] Solution of canonical fo: p(t, I,p.) of differential eq's. with initial with initial conditions J conditions I. (j = 1,..., n) A is the transformation ~ Continuable periodic solutions of [c.(t, i)] Solutions of canonical f the e-dependent family of diff. eq's., corresponding to p(t, p) p(t2,) jp(t, O) = p(t) and p(t,4) = cp(t, I(), j =... n cj(t,O) = Cjo and cj(t, 4) pj(t, A(.),.) I Arbitrary initial conditions A Arbitrary initial condit I( L) Initial conditions for p(t, ) A() Initial conditions for [cj(t,.)] Desired periodic solutions of original (] Transformation of C (t) )(t)| differential eq., Eq. (1.1). It is [j It is assumed that assumed that (J = 1. n) | j(t) = cj(t, 1) (D(t) = p(t, 1)'( i, 2, ~. n)

27 in Appendix A. This approach leads immediately to a so-called first variation of the differential equation with respect to the periodic solution or, more concisely, the variational equation. This variational equation is, as has been mentioned before, a linear differential equation with periodic coefficients. The common period of these periodic coefficients is the same as that of the periodic solution for which the stability is being tested. According to Floquet [17], the preceding linear system has solutions made up of combinations of functions of the form E o u(t), where u(t) is a periodic function with a period equal to that of the coefficients, and a is a constant. The a is referred to as a characteristic exponent, and there are n of them associated with an nth-order variational equation. These characteristic exponents are comparable in nature with the characteristic roots associated with a linear differential equation with constant coefficients. If all the characteristic exponents have negative real parts, then the periodic solution of the nonlinear differential equation is asymptotically stable (see Appendix A). Unfortunately, it is usually difficult and laborious to determine these characteristic exponents. One approximate method is to replace the periodic coefficients in the linear differential equation by their mean values and treat the linear system as if it had constant coefficients. Markus [18] has shown that in certain cases this method gives satisfactory results. In more general problems, methods from perturbation theory, quite similar to those utilized here to develop recursively the periodic solution of the nonlinear differential equation itself, have been applied. However, in the present situation there exists the added difficulty that the solution about which the stability question is being asked

28 is known only approximately. This lack of knowledge causes difficulty because the periodic coefficients of the variational equation are functions of this solution; therefore, so are the characteristic exponents. This means that if the solution is known only approximately, then the "true variational equation" is known only approximately. In such a case, even if this "approximate variational equation" were solved exactly, there would still be some question about the stability. It is, therefore, clear that something more must be done than just finding a method for determining the characteristic exponents of a linear differential equation with periodic coefficients. The key to this problem is supplied by some of the work of Moulton in celestial mechanics [16]. The type of equation considered by Moulton is a linear differential equation that has coefficients dependent on a parameter ( as well as being periodic in t. The form of this type of problem is, = P(t, ()y ( y ) = (1.20) dt where y and y are n-component column matrices, and P(t, ~) is an nxn matrix analytic in 1 for 0 < p < IO and periodic in t with period T. P(t, 0) is a constant matrix. The general solution for this differential equation is of the form, y(t, I, ) = D(t, )cR(I. (1.21) D(t, A) is an nxn matrix analytic in p for some finite p interval and periodic in t with period T. R(p) is an nxn matrix analytic for p. in some finite p. interval. I is an n-component column matrix of initial values. The n characteristic exponents of Eq. (1.20) are denoted aki(,) (i = 1, 2,..., n) and are solutions of the equation,

29 det [R(a) - E] = 0, (1.22) where E is the identity matrix. The basic idea behind Moulton's approach to the problem of determining the aCi() is the recognition that there exist n linearly independent choices for the initial value vector, I, such that the solution associated with each of these choices is of the form,:i(I)t y(t, Ii(), ) = e Ui(t, ), (1.23) where I.(A) is the ith choice of initial condition vector I, and U.(t, L) 1^~~~~~~ ~~cia.(c)t l is an n-component periodic column matrix with period T. e is a scalar multiplier. A different characteristic exponent, ai, appears in this multiplier for each choice of an initial value vector from the abovementioned linearly independent group. The immediate consequence of this selection of initial conditions is that a complete set of linearly independent solutions to Eq. (1.20) is obtained, each member of which contains one and only one characteristic exponent. At this point Moulton introduces n linear time-dependent transformations of the form, a t y = ei U (1.24) After substituting Eq. (1.24) into Eq. (1.20), the following linear differential equation for U is obtained; U [= [-E + P(t, A)] U, (1.25) where E is the identity matrix. The characteristic exponent that appears in this equation is as yet unknown. It is known that one solution of the above differential equation is periodic in t with period T. This follows

30 from the discussion regarding Eq. (1.23). It also follows from the same discussion that the periodic solution is Ui(t). All this assumes, of course, that the proper value of ai has somehow been determined and been substituted into Eq. (1.25). The technique which is presented by Moulton for treating Eq. (1.25) is a recursion procedure somewhat similar to that used in Chapter IV to determine the periodic solutions of the nonlinear differential equation characterizing the nonlinear feedback system. The significant difference is that in addition to developing the periodic solution as a power series in p the characteristic exponent ai () is also developed as a power series in. The above method is modified and applied in Chapter V so that it constitutes a technique with which the stability question can be answered.

CHAPTER II DEVELOPMENT OF CANONICAL FORMS 2.1 Introduction The purpose of this chapter is the development of two canonical forms for the restatement of the original differential equation, Eq. (1.1), as a problem amenable to treatment by perturbation techniques. Assuming that the desired periodic solutions, D(t), of Eq. (1.1) are dominated by S terms of frequency 2, the first canonical form is intended for those systems in which the forcing function, e(t), contains a term of frequency S, and the second canonical form is intended for those systems in which e(t) does not contain a term of frequency 2s 2n A discussion of certain aspects of perturbation theory which is meant as a background for the actual development of the canonical forms is given in the next two sections of this chapter. The first of these two sections deals with what is called in this study the "normal perturbation problem;" the second deals with what is called in this study the "inverse perturbation problem." As these names imply, some type of inverse relation exists between these two problems. As is demonstrated in the next sections, this inverse relation arises from the fact that the "normal perturbation problem" is essentially an analytic problem, whereas the "inverse perturbation problem" is essentially a synthetic problem. The actual development of the two canonical forms is carried out in the last two sections of this chapter. 2.2 Normal Perturbation Problem Although, so far as application of perturbation techniques to physical problems is concerned, the "normal perturbation problem" is the 31

32 second stage and the "inverse perturbation problem" is the first, these two "problems" are discussed here in reverse order of appearance because a prior understanding of the "normal perturbation problem" facilitates comprehension of the "inverse perturbation problem." The word "normal" has been used in the title of the "normal perturbation problem" merely because it is "normally" the problem considered in the mathematical literature. Conversely, the "inverse perturbation problem" is usually given scant attention, if treated at all. In fact, the motivation for making a definite distinction between these two "problems" is the desire to emphasize their existence and their fundamentally different character, i.e., analysis versus synthesis. At the outset of a "normal perturbation problem," a >-dependent (u a parameter) family of differential equations is given. One member of this family, say the "i = 0" member, is assumed to have one or more known periodic solutions. Presumably, other members of this family also have periodic solutions, and it is reasonable to suspect that members of this family which are neighbors of one another may have, in some sense, neighboring periodic solutions. In particular, members neighboring the "L = 0" differential equation may have periodic solutions which neighbor, in some sense, certain known periodic solutions of the " O = 0" differential equation. However, it is easy to construct counterexamples which demonstrate that this suspicion is not always justified; i.e., x + px + x = 0 has a periodic solution for only. = O. On the other hand, there do exist L-dependent families of differential equations that do possess solutions periodic in t and continuous in u over some u interval. This type of periodic solution for the i-dependent family of differential equations is referred to in this study as "continuable." The "normal perturbation problem" is to discover the existence of and determine

33 expressions for continuable periodic solutions which approach distinct periodic solutions of the ".l = 0" differential equation as i -e O. An example of such a continuable periodic solution for a time-dependent system is shown in Fig. 2.1. Note that this figure is not the usual phase-space representation of the solution trajectories of a differential dx equation. It merely portrays the functional relation between x, d, and u for one continuable periodic solution, and it does not portray the more general or aperiodic behavior as would be the case in the usual phasespace representation. Also, the disappearance of periodic solutions for u > 3 should not be considered the only manner in which a solution can fail to be continuable. It can now be seen that the "normal perturbation problem" actually consists of three problems. The first of these is the problem of determining which, if any, of the known periodic solution of the "i' = 0" differential equation can be continued away from 4 = 0. These periodic solutions are referred to in this study as generating solutions and can be considered as first approximations to the continued periodic solutions for values of j other than zero. The second of these three problems is to determine the composition of these solutions for, f O. This is usually accomplished by a recursion procedure which develops the solution as a power series in Ai. The third problem is to determine the 4 intervals over which the periodic solutions can be continued. Whereas the two preceding difficulties can be handled more or less easily, this last problem still awaits complete solution. In this study a Jacobian existence condition is developed which guarantees the A interval for \1\ sufficiently small. As was mentioned in Chapter I, Lewis [15] has published results regarding this

34 [L2 / I i//dx Ay l dt LOCUS OF INITIAL CONDITIONS FIG. 2.1 EXAMPLE OF A CONTINUABLE PERIODIC SOLUTION

35 problem. A convenient method for graphically presenting the information regarding allowable i intervals is shown in Fig. 2.2. Each horizontal line in this figure is associated with a particular periodic solution and has a projection-on the a-axis which corresponds to the [i interval over which the solution in question can be continued. For example, solution number 1 can be continued over an interval including 4 = 0 and =.1' Solution number 2 can be continued away from 4 = 0, but not all 3 /____ 2 X /~=o ~= FIG. 2.2 GRAPHICAL REPRESENTATION OF THE ALLOWABLE i INTERVALS FOR A "NORMAL PERTURBATION PROBLEM" the way to u1. Solution number 3 can not be continued at all. Of these three examples, numbers 1 and 2 would satisfy the aforementioned Jacobian condition at u = 0, and solution number 3 would not. It is assumed here that the periodic solutions are denumerable, but if they are not, the

36 insight gained with this figure is still valid. Finally it can be seen, from the above statement, that the "normal perturbation problem" is the analysis of (i.e., the finding of the periodic solutions of) a given differential equation. 2.3 Inverse Perturbation Problem In this case a differential equation is also given at the outset of the problem, but it does not depend upon the parameter p.. It is desired to find the periodic solutions of this original differential equation with the aid of perturbation techniques, and in order to do this the original differential equation must be imbedded into some "normal perturbation problem." The selection or synthesis of this "normal perturbation problem" constitutes the "inverse perturbation problem." This selection can be divided into two separate parts: (1) the selection of a ". = 0" differential equation, and (2) the selection of a uI-dependent family of differential equations which contains both the "p = 0" and the original differential equations as members. It is, of course, not to be expected that every imbedding will be a satisfactory solution to the "inverse perturbation problem." To begin with, there is the unsolved problem of allowable p intervals. The sizes of the various intervals depend on the selections of both the "i = 0" differential equation and the i-dependent family of differential equations. A diagram, similar to Fig. 2.2, which presents a possible configuration for the continuable periodic solutions of some v-dependent family of differential equations is shown in Fig. 2.3. In this figure p = Kl is assumed to correspond to the original differential equation. Solution number I is an example of the desired type of continuable solution. Its

37 7 I a V 6 FIG. 2.3 GRAPHICAL REPRESENTATION OF THE ALLOWABLE U IUTERVALS FO.R ANI "TNVERSE PER-rURBAPTON PR9OBLEM.I" allowable 1 interval includes both 0 = O and 1 = 1L Considering it from left to right the continuable periodic solution of the = differential equation, which satisfies the Jacobian condition, can be thought of as predicting a periodic solution of the "a = o1" or original differential equation. The word "predicting" is used here to emphasize the danger which exists as long as no actual guarantee of the i intervals is employed. For instance, solution number 2 satisfies the Jacobian condition, but predicts a non-existent periodic solution of the it = I1" differential equation. Viewing solution number 1 from right to left, it can be seen that a periodic solution exists at A = Ct1 which is predicted by a periodic solution-at AI = 0. Solutions numbers 4 and 5 exist at L = >1' but are not predicted by periodic solutions at L = Ok Solution number 6 is continuable over a satisfactory I1 interval except at one point, "a". This point could be the transition 4 FIG. 2.3 GRAPHICAL REPRESENTATION OF THE ALLOWABLE ~ INTERVALS FOR AN "INVERSE PER'XURBATION PROBLEM" allowable ~ interval includes both \ = 0 and ~ = Kl' Considering it from left to right, the continuable periodic solution of the " = 0" differential equation, which satisfies the Jacobian condition, can be thought of as predicting a periodic solution of the "g = u1" or original differential equation. The word "predicting" is used here to emphasize the danger which instance, solution number 2 satisfies the Jacobian condition, but predicts a non-existent periodic solution of the "g = ~1" differential equation. Viewing solution number 1 from right to left, it can be seen that a periodic solution exists at ~ = ~l which is predicted by a periodic solution at g = 0. solutions at g = 0. Solution number 6 is continuabie over a satisfactory { interval except at one point, "a". This point could be the transition between stable and unstable periodic solutions. Solution number 7 does

38 not exist at either A = 0 or A = p1; therefore, it would have small importance as far as the original differential equation is concerned. It is now apparent that the relation, in the above context, between the "a~ = 0" and the original differential equations is a completely mutual relation; i.e., either of these equations may be viewed as a perturbed differential equation whose periodic solutions are continued until, assuming continuability, they become periodic solutions of the other. The only difference between them is that one happens to have periodic solutions which are relatively easy to find. As was previously mentioned, the a-interval problem is only partially solved in this study, and a Jacobian existence condition is relied upon; i.e., every periodic solution of either the "p. = 0" or the'u = l'" differential equations is assumed to be part of continuable periodic solution similar to solution number 1 in Fig. 2.3. Now along with those factors which relate to the.-interval problem, there are several desirable characteristics which the selected "i = 0" differential equation should exhibit. Naturally, it should possess periodic solutions; however, this is not enough. Unless these solutions are easily found, there is little point in introducing perturbation techniques. Some of these periodic solutions should satisfy sufficient conditions for their being continuable away from p. = 0, and the ascertaining of whether or not these conditions are satisfied should not be too laborious. These continuable periodic solutions are, of course, the generating solutions. Finally, these generating solutions should be good approximations to the corresponding periodic solutions of the original differential equation. With the above desired characteristics as a background, consider the possible choices for a "p = 0" differential equation. Nonlinear

39 differential equations can be immediately ruled out simply because their periodic solutions are not usually easy to find. The fact that nonlinear differential equations exist which have easily determined periodic solutions does not help a great deal. The main interest of this study is the creation of an analytic tool which can be used upon a wide range of nonlinear feedback systems. Therefore, the "t = 0" differential equation must be drawn from a class of differential equations, any member of which would be a satisfactory "l = 0" differential equation. The class of nonlinear differential equations is notoriously not such a class. However, it does seem possible that certain subclasses of nonlinear equations might be constructed —perhaps all nth-order differential equations with xj(t) = a. cos (ut + pj) as solutions —but this is outside the scope of this study. The range of possible choices for the "' = 0" differential equation is, therefore, rather quickly reduced to the class of linear differential equations. Linear differential equations with time-dependent coefficients are eliminated just as quickly using roughly the same argument. Consequently, the " = 0" differential equation becomes restricted to the class of linear differential equations with constant coefficients. The problem, then, becomes one of determining which linear differential equationwith constant coefficients is satisfactory for a given original differential equation. Needless to say, as in any synthesis problem, there is no single answer. The first decision to be made is whether to choose a homogeneous or an inhomogeneous differential equation. Some of the statements which follow regarding this decision are necessarily somewhat vague. Their clarification occurs in Chapter III, where the details of finding the generating

solutions are developed, especially the details of the Jacobian condition. It is shown in Chapter III that an inhomogeneous "I = 0" differential equation has continuable periodic solutions in cases where a related homogeneous equation would not have them; i.e., in the latter situation the Jacobian is not satisfied. Obviously, in such cases the inhomogeneous differential equation should be chosen. On the other hand, if the forcing function of the inhomogeneous equation contains a term equal in frequency to a natural mode of undamped oscillation for the reduced homogeneous equation, a resonant response of the form, t sin cwt, results which implies that the "i = 0" differential equation cannot have periodic solutions. There are other, less important, differences between homogeneous and inhomogeneous ", = O" differential equations. For instance, a givenorder inhomogeneous differential equation can have more involved generating solutions, whereas an nth-order homogeneous equation is limited to generating solutions which contain, at most, 2 different frequencies. This could be important in relation to the first approximation problem. Differences of similar importance between homogeneous and inhomogeneous "' = 0" differential equations will become apparent in the course of this study. Whether the selected equation is homogeneous or inhomogeneous, the problem of choosing the location of the characteristic roots exists. Those characteristic roots whose associated characteristic solutions are periodic have, of course, real parts equal to zero. Further, since the only solutions of interest in this study are those that are periodic with a period T (T = qTe), periodic characteristic solutions which form part of a generating solution must have period T (not necessarily a least period). Add to this latter consideration the assumption that the periodic solution

41 of the original differential equation is dominated by a term of frequency, and the first choice for the location of the characteristic roots becomes +~j and the remaining (n-2) with nonzero real parts. For a homogeneous "i = 0" differential equation of this type, the generating solution would be jwt -j t 5jst * 5jst p(t) = c10 e + c10 Assuming that the order of the "a = 0" differential equation were four or higher, this approach could be carried another step, and two characteristic roots could be placed at, for example, t j3ws (i.e., some harmonic suspected of contributing a significant term to the solution) in addition to those at +jos. Presumably, this would lead to a generating solution which was a better first approximation, but, because the coeffi+j3c t -jwst * cients of E and e as well as c10 and c10 would have to be determined, it would also lead to a generating solution whose determination was more laborious. This procedure could be continued until all the characteristic roots were at various, separate, integer multiples of j2 If singular perturbations [19] were allowed, a "4 = 0" differential equation of higher order than the original differential equation could be used; therefore, a generating solution of arbitrary complexity could be considered. However, at this point the determination of the generating solution becomes as formidable as the determination of the periodic solution of the original differential equation. After the " 0 = 0" differential equation has been chosen, a Udependent family of differential equations must be selected. Clearly, 1 Recall that Cas = 2 and that Ts = (p)T where q is a rational number. o~~~~~

42 there are an infinite number of these families which contain both the "4 = 0" and original ditf'erential equation. The somewhat arbitrary choice made in this study was one of the simplest: the difference between the original and "i = 0" differential equations was multiplied by 4 and added to the "I = 0" differential equation. Therefore, the original differential equation corresponds to the "p = 1" member. The implications of other choices for the i-dependent family have not been investigated by the author. Note that in the above discussion the method for choosing T s is not actually given. The reason for this is that there does not exist a straightforward a priori procedure for ascertaining T. The usual procedure employed is simply to make an intuitive guess. In relatively simple nonlinear systems this guess is usually satisfactory. 2.4 First Canonical Form With the preceding three sections as a foundation, it is now possible to appreciate the steps taken in the following development of the first canonical form. This canonical form can be considered a partial solution of the particular "inverse perturbation problem" associated with the nonlinear feedback system characterized by Eq. (1.1). The part of this "inverse perturbation problem" solved in this section is the selection of a "p. = 0" differential equation and a 4-dependent family of differential equations. The part left unsolved is the determination of the allowable ranges for pL. Equation (1.1) is repeated below. A(z)x - B(z)[F(x) + e(t)] = 0 (1.1)

43 The symbols in this equation are defined on pages 4 and 5 of Chapter I. It can be seen from the discussion in the preceding sections of this chapter that the first move in the solution of this "inverse perturbation problem" is the selection of a "0 = 0" differential equation. As was previously stated, the type of "p = 0" differential equation selected in this study for the first canonical form is a homogeneous linear differential equation with constant coefficients. However, the selection of a satisfactory "I = 0" differential equation from this class of differential equations is not straightforward. This is especially true when Eq. (1.1) is of large order (i.e., n large). Even assuming that two of the n characteristic roots of such a "I = 0" differential equation are known (e.g., + = +j and X2 = -jis), the selection is still complicated by the existence of a seemingly-infinite number of possibilities for the location of the other (n-2) characteristic roots, and, in the sense that there is an infinite number of linear differential equations with constant coefficients "close to" the original or "p = 1" differential equation, an infinite number of these possible "4 = 0" differential equations is satisfactory. Moreover, there exists another consideration which affects the selection of the "0 = 0" differential equation for both the first and second canonical forms. In addition to having a "L = 0" differential equation that satisfies all the requirements connected with the perturbation problem, it is desirable for the selection of a "I = 0" differential equation to result in recognizable and tractable functions appearing in the later stages of the analysis. For instance, it is obviously desirable for any selection of a "p = 0" differential equation to keep the functions A(z), B(z), F(x), and e(t) somehow intact and recognizable

44 as building blocks. It is this last consideration, together with the considerations discussed in the previous sections, which motivates the particular choice of the first canonical form which follows. The main idea behind this choice is the development of a formal procedure for selecting suitable locations for all the characteristic roots associated with the "i = 0" differential equation. As a first step in this procedure for the first canonical form, introduce a real constant, N, into Eq. (1.1). This constant, N, serves in the formation of a linear approximation for F(x): Nx m F(x). It is not, however, equivalent to a describing function which appears later and is denoted n. The following differential equation is obtained from Eq. (1.1) by adding the quantity, -B(z)Nx, to both sides of the equation and rearranging some terms. A(z)x - B(z)Nx = B(z)[F(x) - Nx] + B(z)e(t) (2.4) Now consider the linear differential equation with constant coefficients which is obtained if the left side of Eq. (2.4) is equated to zero, A(z)x - B(z)Nx = 0. (2.5) The characteristic roots of the above differential equation are a function of the real constant, N. It is assumed that for one or more values of N, two characteristic roots are purely imaginary and the complex conjugates of one another. Let these values of N be designated N.; denote the corresponding pair of imaginary characteristic roots by + jo. Thus, Ii

45 the differential equation, A(z)x - B(z)N.x = 0, (2.6) has ~ jwNN as two of its characteristic roots. The other (n-2) characJ teristic roots are assumed to have nonzero real parts and to be distinct. Therefore, Eq. (2.6) has the following periodic solution: +jLN t ".jt x(t) = clo E + c10 (2.7) (Since throughout this study only real solutions are of interest, complex conjugate coefficients, clO and clO, are used.) Equation (2.6) has, therefore, some of the features desired for a "a = 0" differential equation. However, the angular frequency, aN, of the periodic solution Eq. J (2.7) does not necessarily have the proper relation to e, the fundamental angular frequency of the forcing function. Recall that it was desired that the least period of the periodic solutions of "a = 0" differential equation be related to the least period of the forcing function by a rational factor; i.e., Ts = (p)Te. Therefore, the characteristic roots, + jWN,must be shifted somehow to + jw. J The shift of these characteristic roots is accomplished by subtracting selected terms from Eq. (2.6) so that if the characteristic roots of the equation were located at l, X2,..., n, then the characteristic roots of the resulting differential equation would be located at Al = 51' = 2 = ~x2,-. An = An (3 is a real constant). In particular, X1 = jws and 2 -jo. The selection of the terms which are to be subtracted S can be most easily understood with the aid of a transformation of the time scale.

n A(z) z + C(z). (2.8) Equation (2,5) becomes z x + C(z)x - B(z)N.x 0. (2.9) Transfo rm t to a new independent variable, T, with the linear transformation, t = T. (2.10) This transformation implies that k k dk 1 dk k 1 k d i d. k lk k'k dk lk dk or =., w, (2.11) dt CT k where w. Equation (2.9) becomes d-c wnx + 3nC( )x - 3n ( )N.x 0. (2.12) 3 3 This is, of course, still a linear differential equation with constant coefficients that has a known periodic solution. This periodic solution is +J T. * 1N.3 X(J) - c e +* c (2.13) The angular frequency of this periodic function of T is 1ON. If 8 is chosen such that -3 (2.14) 3

then the angular frequency of the function x(T) becomres ws. Correspondingly, the characteristic roots of Eq. (2.12) will be at Xi p \ = 3..., k = pk n n If in Eq. (2.12), t is substituted for'r and z for w, iq. (2.12) becomes a suitable " = 0" differential equation. The terms which rmust be subtracted from Eq. (2.9) in order to obtain Eq. (2.12) are [C(z)x - j(B(z)x. - ) [c( )x - NJ3()x (r1 The choice for the!' = 0" differential equation is then x + [c() - N ()]x = 0. (2.16) The choice for the 4-dependent family of differential equations is znx + nC()x - NB3()x] = [n{C()x - N.B()x} - (2.17) {C(z)x - NjB(z)x} + B(z) {F(x) - Njx + B(z)e(t)] Irmediately notice that the terms in Eq. (2.17) multiplied by, are nothing more than the difference between the original differential equation, Eq. (1.1) and the'" = 0" differential equation, Eq. (2.16). The member of this family corresponding to i = 1 is, therefore, the original differential equation. Also, note that most of the desired functions have been kept intact. Even A(z) is present and intact, although it doesn't appear that this is the case in Eq. (2.17). The "g = 0" differential equation has periodic solutions of the form,

48 +j&2t 88 PP + C-, from which the generating solutions, p(t) = c1 e. + clO e.,must be selected. Since Eq. (2.17) includes both the "1. 0' and the'e J differential ecuations, it is a possible a-dependent family of different.ial equations. However, it is not in a convenient form for further analysis. A representation of this equation as a systemrr of n first-order equations is a form which leads to easier manipulation. The coordinate transforna:tion from the variation-of'-constants technique is utilized below to transform Eq. (2.17) into such a form. This coordinate transformation is as follows: XIt n X..e.. E r1 n x ^ ~Xle lt.... \ e r. (2.19) (n-l (n-) t A A pre5sed in x f. ss " = 0" differential equation, and the r.'s are the new coordinates. ExP~res nmti om-E.( eoeJ pressed in matrix form, Eq. (2.19) becomes 1 Throughout this study symbols are used for matrix equations that are similar to symbols used in scalar equations. In case a possibility of confusion exists, the type of equation intended will be specifically stated.

49 where x r, r2 X = X r= rz 1. L * J,(n-i(n-1(n-) (n ) x2 I i.... and AXt e 0..... *...0 Y. 0 0 e. 0.... (2.22)................ 00 x2t..~~...o....0 6

50 -1 Note that V is the Vandermonde matrix which can be relatively easily inverted (see Appendix E). A direct consequence of Eq. (2.19) is that the following relation is satisfied for all t. (2.23) l1t X2t Xnt E e e r 2 1 2 %2t nt dr. X1e 2 n 3 j - dt (n-2). t (n-2)X2t (n-2) nt 1 2 ^n n This relation can be verified by a consideration of the time derivative of Eq. (2.19). Note that this relation involves n unknowns and (n-l) equations. Now substitute the appropriate components of Eq. (2.19) into Eq. (2.17). The resulting equation is (2.24) ( ) { (n-2) t + x n(n-1) Xnt It Ein th It n t j=z ( j=l j Note that the subscript "j" of N. in Eq. (2.24) is not a summation index.

51 functions of time they can be considered to be constants as far as the linear operators B(z) and C(z) are concerned. This is true because these operators do not involve derivatives of order higher than (n-l). Referring to Eq. (2.19), it can be seen that derivatives of x of order less than n do not involve derivatives of the r.'s. Therefore, the J operation by B(z) or C(z) and multiplication by the r.'s can be commuted as follows: \%t t it n t B(z)[ r1 +... + r n = rnB(2)c +... +rB(z) (2.25) For the same reason, B(z)F contains no derivative of an r.. J Further, from Eq. (2.16), n Ci n [c(s) - NjB()] e) (2.26) x.t This is true because each e is by definition a characteristic function for this "p = 0" differential equation. Substituting Eqs. (2.26) and (2.8) into Eq. (2.24) and combining the resulting equation with Eq. (2.23) gives the following system equation: -hlt ~^t 2 = iVd E Lgl(t'rl'''n) (2.27) -k t n = ~Vnn n () rn = 4"vnn e gl(t,rl',...,'' rn)' where the v. are elements of the matrix V [Eq. (2.21)] and on

52 %It %nt gl(tr r) = B(z)F{ r1E +... + rn (2.28) n X.t + B(z)e(t) - Z r. A(z)e j=l The above two equations constitute the first canonical form. Note that periodic solutions of Eq. (1.1) do not transform into periodic solutions of Eq. (2.27). The needed "periodicity" condition is developed in Chapter III and is of the form X.T Cj(Tp) C - cj(O,4) = 0 ( = 1,2,...,n) (2.29) where rj = cj(t,.) (j = 1,2,...,n) is the desired solution of the canonical form in question. 2.5 Second Canonical Form The second canonical form is developed in a manner that is similar to the development of the first canonical form. The basic difference is that the term B(z)e(t) is not included in the group of terms which is multiplied by p; therefore, the "p = 0" differential equation becomes an inhomogeneous differential equation. The main reason for introducing this second form is that the first form does not lead to a satisfaction of the Jacobian condition if the frequency spectrum of e(t) does not contain a term of frequency w. Therefore, the second canonical form pre5 supposes that the frequency spectrum of e(t) does not contain a term with frequency 2-. Moreover, as was mentioned previously, if this assumption were not made, the "I = 0" differential equation for the second canonical form would not have a periodic solution. The "I = 0" differential equation for the second canonical form

53 is the following inhomogeneous equation: znx + P [1c() - NjB(z)] = B(z)e(t). (2.30) The "t = 0" differential equation has periodic solutions of the following form: jw t -jw t Cp(t,I) = P1o c + P20 + h(t;Cs) (2.31) where h(t;w ) is the particular solution corresponding to the forced oscillation of this linear system caused by B(z)e(t). Here, h(t;w ) is considered to be a function of w (W = P a) as well as t in order to s s q e simplify later consideration of the Jacobian condition. e(t) can be expanded in a Fourier series as follows: ~ J kTt e(t) = E E iT c (2.32) =-oo (Note that WT is used instead of We. Since e = qT only every qth E (Note that mT e e q OIIt can be nonzero.) Therefore, h(t;ws) has the following form: S j WTt h(t;w ) = E 6~~ (2.33) (=-_ _ [(jLT)n + n{C( T) - (is = Pr) Since the e(t)'s which can be used with the second canonical form do not contain terms of frequency s, h(t;w ) does not contain a term of so frequency S~. The generating solutions, p(t), are, of course, a set of periodic solutions of the " 0 = 0" differential equation, Eq. (2.30), which must be selected from among the qp (t,I)'s in Eq. (2.31). The a-dependent family of differential equations similar to

54 Eq. (2.17) of the previous section is: znx+ n[C()x - NjB(z)x] = B(z)e(t) + 4 [n{c()x - NjB()x} (2.34) -C(z)x - NjB(z)x} + B(z) {F(x) - Nx}]. Now, similarly to Eq. (2.19), introduce the following coordinate transformation: Xt %nt x En*e n........ r1 h(t;w ) lt Xt. * * * * * * * * * * * * * * * *. * * * * * * * * n 2 s + (2.35) ^n^ (n-i) klt (n-l) knt (n-l) x n c......X E n h (t;) It can be shown that Eq. (2.23) is satisfied by the rj's defined in Eq. (2.35). Substituting Eqs. (2.33) and (2.35) into Eq. (2.34), gives (n-l) lt. (n-l) t n X.t r1X\ e +. + rnn E [-A(z)h(t; s) - rjiA(z) j=l (2.36) xkt X t + B(z)F{re + + rne + h(t;s)}] Combining Eqs. (2.35) and (2.36), gives the second canonical form: -klt rl = v ln g2(t,rl,.,rn)................ (2.37) -k t r = v nn g(t,rl,...,rn ),

55 where kit k t g2(t,rl,...,r) = B(z)F rle +... + r e + h(t;s)} 2 V'.'rn 1 n s (2.36) n X.t - r.A(z) J -A(z)h(t;w ). j=l s Therefore, as far as the r. variables are concerned, formally the two canonical forms are quite similar. In particular, both forms reduce to the same equation at i = O. However, it should be kept in mind that the untransformed "I = 0" differential equations, Eqs. (2.16) and (2.30), are, respectively, homogeneous and inhomogeneous. Finally, it should be realized that there does exist some latitude regarding the two canonical forms which are presented in this chapter. It is possible that a mixture of the two might be used. A part of the forcing function e(t) could be included among those terms multiplied by A, and the remaining part could be independent of L.

CHAPTER III FIRST APPROXIMATION AND EXISTENCE CONDITIONS 3.1 Introduction Two canonical forms, Eqs. (2.27) and (2.37), for the restatement of the original differential equation, Eq. (1.1), as a form which is amenable to treatment by perturbation techniques, are developed at the end of Chapter II. Each of these two canonical forms is a system of n first order l-dependent differential equations in n variables (rl, r2,.., rn). The solution of this system, which takes on the arbitrary initial conditions A=(al, a2,..., a ), is denoted r. = pj(t,A,u) (j = 1, 2,..., n). The desired solutions of this system, denoted [cj(t,L)]j = [pj(t,A(I),,)j [where A(,) is a special set of 1-dependent initial conditions], are those solutions which are analytic in. for Jl\ < 1 and which satisfy the following "periodicity" condition for 11j < 1. X.T 3 3 c.(T,)e J - c (0,=) =, 2,..., n) This requirement is called a "periodicity" condition because it is the necessary and sufficient condition for the periodicity in t with period n X.t T of p(t,p), where p(t,4) = Z c (t,)e J [see Eqs. (2.19) and (2.35)3. j=-l As i - O, a desired solution, [cj(t,pL)], approaches a solution, [c jo, of the transformed. " = O" differential equation. This solution, [c o], 1 Note that there are actually two "S = 0" differential equations for a given problem; (1) the untransformed "1 = O" differential equation [e.g., Eq. (2.16)], and (2) its transformed representation in the rj variables [e.g., Eq. (2.27) with A set equal to zero]. Similarly, the generating solutions for a given problem occur in pairs; (1) the continuable periodic solution, p(t), of the untransformed "n = 0" differential equation, and (2) its transformed representation, [cjo], which is a solution of the transformed representation of the "G=0" differential equation. 56

57 is referred to as the generating solution for the desired solution which approaches it as A 0. In this chapter, perturbation techniques are employed for the determination of the generating solutions. In the following chapter, a recursion procedure is developed for the determination of the desired cj(t,l)'s as power series in i; i.e., cj(t,1) = c. + i 2 cjl(t) + i cj (t).... The method used throughout both these chapters is an adaptation of the one employed by Coddington and Levinson ([10], pages 363-4) for the treatment of a perturbed system which is closely related to the class of systems encompassed by both the canonical forms in this study. Before becoming enmeshed in the details and proliferation of symbols which follows, it is worthwhile to consider a simple outline of the argument presented in this and the following chapter. Referring to Eqs. (2.27) and (2.37), it can be seen that both canonical forms have the following " 0 = 0" differential equation. 0, -:1 =0 dr. r2 0 ( dt r = 0 n Therefore, the general solution of the "p = 0" differential equation is an arbitrary set of constants, (al,..., an). The purpose of this chapter is to develop a method for choosing those solutions (al, a2,..., a ), if any, of the above " = 0" differential equation which are generating solutions [Cjo]. The approach to this problem has two parts, (1) the development of two necessary conditions for the existence of a [c.(t,)] which

58 select what might be called tentative generating solutions, and (2) the development of a set of sufficient conditions for the existence of [cj(t,.)] for 1 1 sufficiently small. These sufficient conditions are suitable for determining whether or not a tentative generating solution is a true generating solution. At the outset of the development of the aforementioned two necessary conditions, a solution, [c.(t,)l] of the desired character is J assumed to exist. The two necessary conditions are then developed by considering the behavior, as. -> O, of the equation which is obtained by requiring the following integral equation expressions for the c.(t,4)'s, J which are obtained from a straightforward integration of either canonical form, cj(t,) = c (0O ) + vjn j e g(s,cl,...,cn) ds, (j = l,...n) to satisfy the "periodicity" condition; i.e., %.T A.T T -X.t (e J -1) c(O,) + vjn e J f c g(t,cl,...,cn) dt = (j = l,...,n). 0o The first necessary condition is simply that the generating solutions, [c. j, must be selected from among those solutions of the "i = 0" differential equation which satisfy the "periodicity" condition. These latter solutions are denoted [pjo] and are a special subset of all possible solutions of the " = 0" differential equation. The generating solutions, [c jo, are, in turn, a subset of the [pjo]. This first necessary condition is arrived at by merely noting that as 4 - 0 the above integral expression for the "periodicity" condition reduces to A T k.T (E J -1) j(0,0) = (e j -1) cj = 0, (j = l,...,n) 3 JO~~

59 The second necessary condition, which in effect selects the [cjo]'s from among the [Pjo]'s, is arrived at by noting from the above integral expression for the "periodicity" condition that for j = 1 and j = 2 the following integrals must vanish for all [ within the domain of existence of [cj(t,) j. T -A.t e g(tc1,..., cn) dt ( = (j =, 2) o X.T That this condition is justified follows from the fact that (e -1) = 0 for j = 1 or j = 2. The above condition must be satisfied at A = 0, and, as is shown later in this chapter, it is suitable for selecting tentative [c. j's from among the [p. ]'s. Jo JO The second part of the approach employed in this chapter is the determination of a sufficient set of conditions for the existence of [c (t,i)]. If the existence of [cj(t,4)] is not assumed a priori, then the previous two necessary conditions can be thought of as selecting tentative generating solutions, [c. j, which may or may not be part of jo a continuable solution, [cj(t,P)]. It is shown in this chapter that by viewing the integral expression for the "periodicity" condition as an implicit relation between A and the initial conditions of the desired solution, [cj(t, )], the implicit function theorem can be employed to develop a Jacobian condition which, if satisfied, guarantees the existence of [cj(t,{.)] for 1lj sufficiently small. This Jacobian condition can be thought of as selecting the actual generating solutions [cjo] from among the tentative generating solutions. Throughout this chapter the expressions which determine the first approximations or generating functions are reduced to describingfunction-type relations. Furthermore, the last section of this chapter

60 presents some useful relations between the Jacobian condition and the behavior of the generating solutions as the forcing function frequency, co 2-, is varied. In Chapter IV, a recursion technique is developed for the determination of the cj(t,p)'s as a power series in p; i.e., cj(t,) = Cjo + C jl(t) + p cj2(t) + * The key tool in the development of this recursion procedure is the requirement that each c. (t) must satisfy the "periodicity" condition. jk 3.2 Determination of the Generating Solutions for the First Canonical Form Consider the first canonical form, which is repeated below: -1lt rl= - vln gl(t,rl'"^,rn ) -2 t r2 =[V2n~ 9l(trl, *,rn) r2 = vW 2 ^^r............... (2.27) -A t n rn= 4VnnC gl (trl''',rn), where %kt k t 1 n n gl(trl,...,r) = B(z)F {r1e +...+rnE n } (2.28) n + B(z)e(t) - Z r. A(z)c j=l o Suppose that Eq. (2.27) has a unique solution [rj = pj(t,A,)] with initial conditions, A> (al...,a ), which exists for t in some interval containing

61 the interval 0 < t < T, and is continuous in A for. sufficiently near = 0. The solutions of Eq. (2.27) which are sought are those that correspond to continuable periodic solutions of Eq. (2.17) with period T, p(t,4) = C[t,I(4),]j [where I(4) denotes a special set of i-dependent initial conditions]. The corresponding solutions of Eq. (2.27) are denoted [cj(t,)j = {pj(t,A(i), k)}, and, just as p(t,) - p(t) as 4 - 0, [cj(t,I)] - [Cj j as 4 -~ 0. Owing to the nature of the transformation, Eq. (2.19), the solutions, cj(t,4), are not themselves necessarily periodic with period T, although the corresponding p(t,4L) is. However, the periodicity of the p(t,4)'s can be utilized to develop a corresponding "periodicity" condition for the cj(t,))'s. Since Eq. (2.17) is, in its explicit dependence on t, periodic with period T, and since the solutions of Eq. (2.17) are unique, the following is a necessary and sufficient condition that a solution of Eq. (2.17) be periodic in t with period T. pk)(T,I, ) - c(k)(oI,) = 0 (k = 0,l,2,...,n-1), (3.1) where cp(k) _. From Eq. (2.19), (2.20), and (2.22), Eq. (3.1) becomes 1 T NT P2(T,A,4)e 1 - p2(0,A,4) V-2 = 0 (3.1A) A T pn(TA,)e n - p (O,A, ) Since the rmi t is n ero the determinant of V is nonzero (as are distinct), the only solution of Eq. (3.1A) is the trivial solution; therefore, the "periodicity"

62 condition on the individual pj(t,A, )'s is X T p (T,A,)e t - p.(O,A,) = 0 (j = 1, 2,.., n). (3.2) Integrating Eq. (2.27) and substituting the result into Eq. (3.2), this "periodicity" condition becomes A T A.T T -.t qj(A,)= (e -)pj(O,Ai) + vvj J E J g(tl,..,pn)dt = 0 (3.3) (j = 1, 2, 3,..., n) [q.(A,4) is introduced for subsequent convenience.] Therefore, the above equation is the necessary and sufficient condition that a solution of Eq. (2.27), [pj(t,A,.)], corresponds to a solution of Eq. (2.17), p(t,.), with the desired periodicity. Considering this condition, Eq. (3-3), as an implicit relation between A and A, the desired type of continuable solutions of Eq. (2.27) corresponds to solutions, A(4), of this implicit relation which are continuous and unique for jI< < 1; i.e., [c.(t,)j = {p.(t,A(),~)}. In order to ascertain the necessary conditions for the existence of A(A) for 1j1 sufficiently small, and, thereby, the necessary conditions for the existence of continuable periodic solutions of Eq. (2.17), p(t,4) for 11j sufficiently small, assume that such a function, A(A)> [al(i),...,an(u)], exists. As -. 0, Eq. (3.3) reduces to.T (e -l)pj[O,A(O),O] = 0 (j = 1, 2,..., n), (3.4) where p [O,A(O),O] = a.(0) (j = 1, 2,..., n). 3 3[,()O jO

63 Since at ~ = 0 Eq. (2.27) reduces to r2 = O, (3.4A) n it is apparent that pj[t,A(O), O] = pj[O,A(O),O] = aj(O); therefore, c. = a.(0). Now referring back to the transformation, Eq. (2.19), that relates dx d^ n x x, it (nil) to (rl, r2,..., r), it can be seen that Eq. (3.4) is equivalent to requiring the solution [Cjo] = {pj[t,A(0),O]} to correspond, through the transformation, Eq. (2.19), to a periodic solution of the untransformed = 0" differential equation, Eq. (2.16). This necessary condition was, of course, expected. Now let [pjo] denote all solutions of Eq. (3.4A) which satisfy the necessary condition, Eq. (3.4). Not every solution [p. j] is a generating solution, [Cjo]: the [c j's are a selected subset of the [pj ]'s. The [Pj ]'s can be further delineated by taking into account the assumed distribution of the'Xts and the requirement that eventual expression of any solution in d (n-l) dx, d ( 1) x Lx,,..e, ~''', -variables must be a real function of t. Since all Xj's except k1 and X2 have nonzero real parts, the condition, Eq. (3.4), can be satisfied only if pjo= 0 for j = 3, 4,..., n. On the other hand, AlT and k2T are both integer

64 XT \T multiples of j2:; therefore, (E -1) = 0 and (E -1) = 0. Consequently, the condition, Eq. (3.4), is satisfied by [pjo] = [P' P2 0, 0,...,0], 2t where p10 and p20 are arbitrary complex numbers. Finally, because e lt is the complex conjugate of e (l = js, X2 = -jOs), and because of 1 S 2 S the reality condition on the x, it is necessary that p20 = p10 In order to determine which p10 corresponds to c10, Eq. (3.3) muT b csedgn2T must be considered again. Because (e -1) = 0 and (e -1) = 0, independently of p., the condition, Eq. (3.3), reduces to -\.t q.(A,) = vj/ e gl(t,,...,'' p)dt = 0 (j = 1, 2) (3.5) 0 for j = l, 2 and \jl sufficiently small. As it appears in Eq. (3.3), this integral is multiplied by p.; therefore, it might be thought that Eq. (3.5) does not hold at p. = 0. However, the integral is continuous in the p.'s, which, by assumption, are continuous in y; therefore, the integral is continuous in p., and the condition, Eq. (3.5), must be satisfied at p. = 0. Therefore, employing the conclusions reached regarding the first necessary condition, Eq. (3.4), at p = 0, Eq. (3.5) becomes q.(A,0) V jnv E g(t, P10' P10' 0,..., O)dt = 0 (j = 1, 2) (3.6) 0 This appears to be an over-determined system with two equations, j = 1 and j = 2, and one unknown, Plo0 However, referring to the definition of gl, Eq. (2.28), it can be seen that gl = gl for all values of [Pjo] = [P.o PlO' 0,..., 0]; therefore, since = the two equations of Eq. (3.6) form a consistent set. The satisfaction of either of the above two equations is the second necessary condition for the existence of a generating solution. Furthermore, all p101s which satisfy Eq. (3.6) can be thought

of as tentative generating solutions. So far, only necessary conditions have been given for the existence of a continuable solution which terminates in [cjo] at p = O. Now the sufficient conditions for the existence of such a solution will be presented, i.e., sufficient conditions that a tentative generating solution from Eq. (3.6) be a genuine generating solution. First, of course, the necessary conditions in Eqs. (3.4) and (3.6) must be satisfied by the tentative generating solution in question. Then the standard theorems regarding the dependence of solutions of differential equations on parameters and initial conditions (see Ref. [10], pages 32-37) guarantee the existence of unique solutions, [pj(t,A,4,)], for Eq. (2.27) which are continuous in the initial conditions and p for initial conditions sufficiently close to those of the tentative generating solution and 4i\ sufficiently small. These latter solutions do not all satisfy Eq. (3.2). The desired solution, [cj(t,)], if it exists, is that one whose initial conditions form a unique, continuous function, A(.L),which satisfies Eqs.(3.3) and (3.5), and is such that A()^- [al(O),a 2(),..., a (O)J = (c1O, C10, 0,..., 0), as p - O, where the generating solution (c10, c10,...,0) is equal to the tentative generating solution under consideration. The implicit relation between A and p, which is obtained through considering Eq.(3.3) for j = 3,4,...,n and Eq.(3.5) for j = 1,2 is, because of its known solution at p = O, in just the form needed for the application of the implicit function theorem. According to this theorem ([20], p. 191) if the following Jacobian condition is satisfied, then for 1j1 sufficiently small a unique, continuous function, A(>), exists such that A(>) [al(O), a2(O),..., an(0)] = (pl, pl, 0,...,0), as-.

66 n 7 q21 6q2 121 E Y 2 y n y J(A,4) = det... o, (3.9) aqn I a6n I L17 27 n y where X implies A = A(O) and A = O. qj (j = 3, 4,..., n) is defined in Eq. (3.3) and qj (j = l, 2) is defined in Eq. (3.5). There can, of course, be more than one tentative generating solution which satisfies the necessary conditions given in Eqs. (3.4) and (3.6). However, each such solution corresponds to a continuable periodic solution of Eq. (2.17) only if the above Jacobian condition is also satisfied. Those which do are denoted (clO, cO1 0,..., O) in order to distinguish them from the tentative encrratin. so-ut:ns, ic, te (p, which satisfy Eqs. (3.4) and (3.6). At this point it is possible with the aid of the relations, Eqs. (3.4), (3.6), and (3.9), to determine which solutions of the " = 0" differential equation can be continued away from A = O. However, these relations are not expressed in a convenient form. A simplification of these expressions is carried out in the remainder of this section. Zeroth-Order Describing Function Relation: Assuming for the moment that a generating solution, [Cjo], is known to exist, substitute [cjo] into Eq. (2.28), the definition for gl, and then substitute Eq. (2.28) into Eq. (3.6). The necessary condition, Eq. (3.6), thereby becomes,since v. (j = 1,..., n) 4 0 [see Appendix Ej, jn

67 T fT t jr * -jWt { f e [B(z)F {c1 + c1O s + B(z)e(t) - 0 (3.10) jo t -j) t cloA(js)e coA(-sj )E jdt = 0 (j = i, 2), where 1 = jcU and 2 =-jo Assume that e(t) Z E e (recall that a = pT), (3.11) 2= ja * jt -ja t oo F {10Ce + c = Fe (3.12) 0 and d f +jm)t -j') t c } I j 1't 1 -x10 { -10 o ) SC 1 C. (3.13) = —oo (This function is needed in the following discussion of the Jacobian condition.) Then, substituting Eqs. (3.11) and (3.12) into Eq. (3.10) and integrating gives the following two consistent, equations: B(jws)FOs + B(jw )E+s - c A(jw') = 0, (3.14) J s 0 - c ams =1, and B(-j )F-S + B(-jw )E - c* A(-jw) = 0, (3.15) s o S 10 ~s - where E and F are the Fourier coefficients of the terms of frequency 0 3s +s + 2 i.e., such that 2eo = W s That E is nonzero, is a prerequisite for the use of the first canonical form. The necessary condition for an allowable c10 is thus A(j3O) +s n -^~^T — + B- =0 (3.16) s 10^C..

68 F+S F-S where n = = _. (a real function of the c10 c magnitude of lc10, see Appendix C) If the existence of [cjo] is not assumed, then the solutions of the above equation are the tentative generating solutions which must be tested with the Jacobian condition. Note that n is equivalent to the normal describing function [21]. Further, note that Eq. (3.16) is formally the same as Eq. (1.5), which is Smirnova's system equation. The difference between Eqs. (1.5) and (3.16) is that the solution of Eq. (3.16) is intended as a first approximation to the periodic response of the nonlinear system to the general forcing function, e(t), rather than to just a simple sinusoidal forcing function, as is the case for Eq. (1.5). The implication here is that for a first approximation to the response, only the term in the frequency spectrum of e(t) equal in frequency to the dominant term of the solution need be considered. The other terms of this frequency spectrum enter into higher-order approximations, which are developed in Chapter IV. Simplification of the Jacobian Condition: Now consider the simplification of the Jacobian condition, Eq. (3.9). Recalling that pj(O,A,p) = aj, it can be seen from Eq. (3.3) that qJJ _J = 0 for j 4 k; j, k 4 1, 2 (3-17) ~k O and Jq. A.T q = (6 J -1) for j + 1, 2. (3.18) ba -

69 Therefore, the Jacobian condition, Eq. (3.9), reduces to ~ _ l q2 a l a ( -1)4 0. (3 19) Lai a 2 a2 1a'i j=3 But all j.'s for j $ 1, 2 have nonzero real parts; therefore, the product in Eq. (3.19) is always nonzero. The Jacobian condition is, therefore, further reduced to -qI 6q 2 6ql 6q2 0o. (3.20) L l aa2 6a2 aly It is shown in Appendix B that following a procedure similar to that used in the derivation of Eq. (3.16), the Jacobian condition, Eq. (3.20), becomes VlnVnB(jws)SO - A(j)} {B(-jw)S - A(-J)} (3.21) in 2n { s } { ( S 0 s} (3.21) V V2 B(j )B(- )s+2s S-2 + o -n 2n s o o +2s where S0 and S are defined in Eq. (3.13). Assuming vln2n B(jw )B(-Jco) ( Eq. (3.21) becomes A(jw )j A(-j ) -2s_ I B(jcD) B(-jw )I If the above expression is considered to be a function of as and Ic10l (Appendix C shows that it is independent ofLcl0 ) and denoted by d, then the equation, d = d(os, I clO), describes a surface above (or below) the Cs, Ic 10-plane. Those combinations of w and IC101 which lie on the intersection of the d-surface with the w, j c0l-plane are the points at which the Jacobian condition is not satisfied. This point is discussed 1 In the remainder of this study, the expressions in (3.21) and (3.22) are both referred to as the "Jacobian expression". In all cases, the intended reference is apparent.

70 again in the last section of this chapter. The determination of generating solutions for the first canonical form has now been reduced to a simple two-step procedure: (1) tentative generating solutions are selected by the condition, Eq. (3.16); (2) these tentative generating solutions are tested with the condition, Eq. (3.22), and those that pass are up-graded to genuine generating solu*^~ [" ~ cdx d(n-)X tions, (cc0, 0,..., 0). Transformed into the x, dt,..., ~0 dt coordinate system these generating solutions become jW t -j* t P(t) 10 c( 3 (3.23) 3.3 Determination of the Generating Solutions for the Second Canonical Form Since the argument which must be made regarding necessary and sufficient conditions for the existence of continuable solutions of the second canonical form follows closely the argument used for the first canonical form, the development of the conditions associated with the second form is only sketched in this section. Recall that the second canonical form is 1 = -vlne g2(t, r,..., rn) -A t -2t r2 = V2nE g2(t, rl,.'.,'n) (2.37) =n nn g2(t, rl.., rn) = n^W " ^(^^ ~.2L where t n t g2(t, r,., n)= B(z)F{rle +... + rn } (2.38) Cn n n k.t - L r. A(z)e J - A(z)h(t; X ) j=l J

71 The "periodicity" condition for the desired solutions of Eq. (2.37) is once again %.T p.(T,A,)e J - p.(O,A,) = 0 (j = 1, 2,.., n). 3 J And, similarly to Eq. (3.3), this condition becomes %.T qj(A,) = (E J -l)pj(OA,A) A T.t (3.24) +Vjn / Pi.? g2(t*, P12..., n)dt = (j = 2.. n). 0 The first necessary condition for the pj(O,A(O), O)'s is i T (e j -1) pj(O,A(O), O) = 0 (j = 1, 2,..., n), (3.25) 3 which has, once again, solutions of the form (Plo1 P1O',..., 0). The tentative generating solutions are then selected from among the solutions of Eq. (3.25) by the following necessary condition which is developed following an argument similar to that followed for the condition, Eq. (3.6). T I c^ g2(t, P0l P10, O,...O)dt = ( = 1, 2) (3.26) 0 The solutions of these two, again consistent, equations are the tentative generating solutions. In order to simplify Eq. (3.26), let jcWt -jw t ~ rt" a==oo F cClO s+ c s +h(t;)} = F e (3.27) 2h where a = T' and F+S F-s o * 0 n = ~ and n = -. (3.28) O C10 O C1

72 +S As before, Fo is the Fourier coefficient of the term in Eq. (3.27) with angular frequency + w. Substituting Eqs. (2.38), (3.27), and (3.28) into Eq. (3.26) and then integrating gives the following simplified expression for Eq. (3.26). A(jw ) n - =o. (3.29) o B(jn) s Note that n is no longer necessarily a real number. The presence of 0 h(t; a ) in the argument of F introduces the possibility of complex n's, s 0 even though F(x) is not a function of x. Moreover, n may now be a func0 tion of w as well as the amplitude and argument of c10l s The Jacobian condition is again Fql a2 ql a q2 ^1l.2] 40. (3.30) a, aa2 aa2 aj Referring to Appendix B, this Jacobian condition reduces to 0 A A(-js) 1 +2s s-2s + 0. (3-31) o B(jw ) 0 B(-j) 0 0 Here it is assumed that jWo t -jU t 00 j &t dF {c o + * E = ~EiLEr SO (3.32) dx + 1010 -3}- 1^0 ~ +2s where C = T, and So is the Fourier coefficient of the term in Eq. (3.32) with angular frequency +20c. — " s

73 3.4 Relations Between the Jacobian Condition and Behavior of Generating Solutions as Functions of the Forcing-Function Frequency A. First Canonical Form It is shown in this section that for either of the canonical forms a simple relation exists for the determination of whether or not the Jacobian condition is satisfied. First Canonical Form: Consider Eq. (3.16) again: A(jw ) E+s n - + c- = 0 (3.16) s 10 If E+S were zero, this equation would reduce to A(jo ) n -, j = o (3.34) o B(j) - s which is the usual describing function for an undisturbed oscillator. However, it is possible for E+S to be zero and e(t) f O. In such a situation the first approximation, given by a solution of Eq. (3.34), would not reflect the presence of a forcing function and would not, in most cases, even result in a generating solution whose frequency was commensurate with that of the forcing function. Moreover, substituting Eq. (3.34) into the Jacobian condition, Eq. (3.22), gives { S0 n} {S- n S2 -2s (335) But from Eqs. (C5) and (C6) in Appendix C, + S = - no} {S- n} (336

74 Therefore, if E+ = 0, the Jacobian condition in Eq. (3.35) cannot be satisfied. Thus, there are two reasons for not using the first canonical form if E+ = O. First, since the n, as defined for the first canonical form, is real, Eq. (3.34) does not usually have a nontrivial solution which is a continuous function of. Secondly, even if Eq. (3.34) is S satisfied, the Jacobian condition cannot be satisfied. Now the conditions under which the Jacobian condition is not satisfied with E+ $ 0 will be investigated. a jQ ii Let c = e, where a" and 0 are real. Equation (3.16) 10 2 becomes A(jw) 2E+s jQ n (a)E - + E = 0. (3.37) B(jo) ) Assume that this equation has a solution of the form: w = (a) and 0 = S s G(a). The following system of equations is obtained by differentiating Eq. (3.37) and its complex conjugate with respect to "a." d A(Js) S+dw d dn 2E+S d ( s 2E -jQ dQ o E -j ds B (jwu) da a da da 2 s L s a (3.38) A(~jwjdw d dn -s d A(-jus) l s. {R2E +jQ dQ 0 o 2E +ji ds BC(-j a da da 2 s LB s, a Assume that the determinant of this system is not zero; i.e., 2E -s +jO d s A(j 2E jQ d (339s j B(jw) j a e' s L 4 I 3' Then the solutions of the system, Eq. (338B(- aJs Then the solutions of the system, Eq. (3.38), are

75 E s dn -s dn s dca a __- a k Z* C* + Ba c ) 40) S 1 -10 10 1o E-S10 da a A(-jw A( ) aa d Es d ^^ ^\ ^ ((^\^- +,, -A(j cJL d. B(-j(.) + dw+ / 1) and o0 d o _ E-s d S doO C10 |d \ a ds 1 From Appendix C and Eq. (3.16) this expression is equivalent to ( ~ _A(j) ( A(-j \) ) -2s ( de _ - 1 dn 0 S 10 S ^o B~O)/ ^o+ B(-j — )/ )"o (o ~3.3 ~~~~~~~da a Eq. (3S 2). Therefore, a sufficient condition for the Jacobian condition not to be satisfied is that the slope of versus for constant E, be zero. Further, except for those special points at which both the numerator and denominator of Eq. (340) vanish, this is also a necessary condition. Since wc is determined by we(cg == P oe), the testing of whether or not the Jaconsiderondition is numerator of the reuces to investigating the slope ofq. (3.40) the generating solution amplitude versus the forcing function frequency, we / dn -s dn E 0 da,0 E E _ C* C I - a 00 (3.42) _ C do da E 2< * BFm Aend C an. S) this e sson i Ceonical For The Jacobian condition for the second(3.). canonical form is Therefore, a sufficient condition for the Jacobian condition not to be satisfied is that the slope of ms versus "a," for constant E+s, be zero. Further, except for those special points at which both the numerator and denominator of Eq. (3.40) vanish, this is also a necessary condition. Since ms is determined by oe(ws = P me), the testing of whether or not the Jacobian condition is satisfied reduces to investigating the slope of the generating solution amplitude versus the forcing function frequency, We 27t B. Second Canonical Form'

76 s B(A ) s (-j s) +2s 2 (3s31) jIo o B(-jcs) o o and the equation which determines the tentative generating solutions is A(jwD) n (a,Q,) - A = 0. (3.29) B(j%) S Assume that this equation has a solution of the form: w = a (a) and S s 9 = Q(a). Differentiating this equation and its complex conjugate with respect to "a" gives A(jw ) an dw ~n 6 an d. _ s_ o S d* _ = ds B(jw ) CU da 69 da 6a (3.44) A(-js ) ) n du n an d [s' o % o ds Q dL B(-jw) ) W da a9 da 6a Assuming that the determinant, o d s o o do S,~ -d A(-jW) 6n Fd A(jW n] a = aT do) B(-J~) ) a- d- B(j%) - + o (3. the solutions of Eq. (3.44) are ds 1 n an an 1 n d 0 0 0 0 (3.4) and 1 O / d A(jw) an A(-jw) an* dQ i1 0_ dn S o 0 d A s (' 7) oda Aa da \dn) B(J%) o oJT" adws B(-ju) ) S s S/S S Noting from Appendix C that j(S~-0) S+2s (C10)

77 1 o +2s-j2@ o= 1 (S~ + S+ E - n ) (C12) it can be seen that the expression in brackets on the right-hand side of Eq. (3.46) reduces to 2j (So -n ( -n - +2s s. a o 0 0 0 0 0 But, from Eqs. (3.29), A(jws) * A(-jws) n = ~ and n =. 0 B(jwo) 0 B(-jos) Therefore, Eq. (3.46) becomes S 2j S (jW A (-j s) -2s ^ ^s ^ S - Ss S (348') |da a B(jOW) 0 B(-jcs) 0 2j This numerator is, except for the coefficient a, the same as the expression which appears in the Jacobian condition. Therefore, as in the first canonical form, a sufficient condition that the Jacobian condition not be satisfied is dws/da = 0. Also similarly to the first canonical form, as long as the determinant, Eq. (3.45), is nonzero, this condition is necessary.

CHAPTER IV HIGHER-ORDER APPROXIMATIONS 4.1 Introduction In the preceding chapter a procedure was developed for determining the generating solutions, [cj ]. In this chapter the next step jo is taken: a procedure is developed for determining continuable solutions [cj(t,)]. This procedure is a recursion procedure which determines the continuable cj(t,u.)'s in power series of A. The corresponding p(t,)'s are developed in the following form: First Canonical Form it co k X.nt 0 k p(t,) = ck(t)k +.+. e + c (t). (4.1) k=O k=O Second Canonical Form ptt ) k Xnt k p(t,) = e clk(t) +... + Z cnk(t)i + h(t;Os), (4.2) k=O k=O where h(t;s ) and the Xj's are as previously defined in Chapter II, and each cJk(t) [not the same functions in Eq. (4.1) and Eq. (4.2)] is a function of t satisfying the following "periodicity" condition: X T j = l,2,...,n cjk(T)e I - cjk(O) = 0, (3) k = 0,1,2,... This will be recognized as a sufficient condition for the periodicity of p(t,4) in t with period T. There is a separate series such as Eq. (4.1) or Eq. (4.2) for each continuable periodic solution p(t,i). At - = 0, these series reduce to the generating solutions as follows: 78

79 First Canonical Form Just * -jCD t p(t,O) = p(t) = cloe + c1e (4.4) Second Canonical Form Juost *-j st p(t,O) = p(t) = clOe + c10 +h(tw). (4.4A) At l = I, the series, Eqs. (4.1) and (4.2), correspond to the desired periodic solutions,<P(t), of the original differential equation, Eq. (1.1). Because of the formal similarity between the canonical forms, Eqs. (2*27) and (2.37), it is not necessary to develop completely separate recursion procedures for each canonical form. Those minor differences that exist between the first and second canonical-form cases are pointed out toward the latter part of the development of the recursion procedure which follows. Finally, in the recursion procedure, a hierarchy of equations arises which can be viewed as a set of describing function relations. Each succeeding equation of this hierarchy is associated with a succeeding step in the recursion procedure. It is shown that the first relation of this hierarchy is —not too surprisingly —exactly the same as that obtained in Chapter III as a necessary condition for continuability of a generating solution. Depending upon the canonical form considered, the first equation of the hierarchy is, therefore, either First Canonical Form, A(jw ) Es no - B(jo = 0 (3.16) s C0 or Second Canonical Form, A(js ) B(j-n) -0. (3 29) The succeeding equations have similar forms but involve more complicated

80 describing functions, which are labeled nk. Again, the initial steps of this development are variations on the work of Coddington and Levinson ([10], pages 356-364). 4.2 Development of the Recursion Procedure Assuming that a generating solution, [cjo], is known, its continuation, [cj(t,L)], exists for \i\ sufficiently small by definition of the term "generating solution." Integration of either of the canonical forms, Eqs. (2.27) or (2.37), gives the following integral equation expression for [cj(t,)]. t k\ls c (t,)) = a1 () + l f g(s,cl...C)ds 0.......................... (4.5) t -A s cn(t, ) = a () + nn n g(scl,.n)ds, where the function g is either gl or g2. A(A) is the set of initial conditions for the desired solution. The first question which must be answered is whether or not the cj(t,))'s are analytic in A for l1l sufficiently small. That this is indeed the case follows easily from a consideration of the usual theorems regarding the dependence of differential equations upon parameters (Ref. [10], pages 32-37) and the implicit function theorem for analytic functions (Ref. [20], page 191). The former theorem guarantees that all the solutions of Eq. (2.27) or Eq. (2.37), with initial conditions within a sufficiently small neighborhood of the initial conditions of the generating solution and with I[r sufficiently small, are analytic in both A and the initial conditions. The latter theorem guarantees that the initial conditions of the desired solutions are analytic in u for |jf sufficiently

81 small. The cj(t,4l)'s can, therefore, be represented as follows: c,(t, =) =. c(t) ( = 1,2,...n), (4.6) k=O for 11J sufficiently small. The initial conditions of such a solution are expressible as follows: aj(L) = ajO + ajl + + 2 + (j = 1,2,...,n). (4.7) The periodicity condition on the cJ(t,i)'s is, from Chapter III, X.T cj (T,)eJ - cj(0,) - 0, (3.2) for IJ{ sufficiently small. Substituting Eq. (4.6) into Eq. (3.2) gives C [c j(T)e - -c() 0 (j 1,2,...n), (48) k=O for |J} sufficiently small. Therefore, each coefficient in this power series of u must be equal to zero, and the condition, Eq. (4.3), is justified. In order to obtain a recursive procedure for developing Eq. (4.6), substitute Eqs. (4.6) and (4.7) into Eq. (4.5), and equate equal powers of i. The resulting recursion formulas are c. (t) = ajo, (a constant) (j = 1,2,.,n) (4.9) and t -js Cj(k+l)(t) = j(k+l)Vjn e gk(S'C lk''nk)ds ( where 00 g(s,cl,...,Cn).= Z gk(SClk, A'Cnk)ll (4.11) k=O

82 and A Cjk (Cjo'CjllCj2.'..cjk) Each cj(k+l)(t) depends only upon the preceding cjk(t)'s, which have presumably been determined in earlier stages of the recursion procedure. It remains only to choose the aj(k+l) s so that: (1) the cj(k+l)(t)'s satisfy the condition, Eq. (4.3), and (2) the cj(k+2)(t)'s can be chosen to satisfy Eq. (4.3). The a.'s are the first set of aj s which must be chosen. Jo jk Since p(t,4) must be real and = jws and k2 = -jO and all other Xj's 2L s sj have nonzero real parts, it is apparent, as well as expected from Chapter III, that the a's which lead to c.'s that satisfy Eq. (4.3) are Jo a10 = a10 (arbitrary) * a20 = a10 a - 0 a30 an 0. Thus, (n-2) of the ajo's are determined by the requirement that the c.'s jo jo satisfy the condition, Eq. (4.3). The remaining two a.'s, aO and a20 are determined by requiring that it be possible to select a set of cl(t)'s which satisfy Eq. (4.3). From Eq. (4.10), the expression for cjl(t) is t -js cjl(t) = ajl + Vjn f e gO[s,ao0]ds (j = 1,2,..,n). (4.13) As all cj's and a.'s, except a10, are determined before this stage, go is considered to be a function of only t and a10 in Eq. (4.13). Applying

83 the condition, Eq. (4.3), to Eq. (4.13) gives k.T k T T -Xt (E J- l)aJl +Vine e gO[t,alO]dt = 0 (j = 1,2,..,n) (4.14) 0 klT k2T Since (e - 1) = 0 and (e - 1) = 0, the above condition for j = 1 or j = 2 reduces to T -k.t e g[t,a10]dt = (j = 1,2). (4.15) 0 These two equations are consistent for the same reason that the two equations in Eq. (3.6) are consistent, and either one of them determines the real and the imaginary parts of a10. This completes the determination of [Co]. Jo It is shown below that gO[t,a10] is a periodic function of t with period T for any a10; therefore, Eq. (4.15) is equivalent to requiring that g (ta0) have no components in its frequency spectrum of frequency + -. This last statement can be justified by merely noting that Eq. (4.15) is the definition for the pertinent Fourier coefficients. Moreover, referring back to Chapter III demonstrates that Eq. (4.15) is the same as Eq. (3.6) or Eq. (3.26). At this point the function gO in Eq. (4.13) is completely determined; consequently, the cjl(t)'s, except for their initial conditions ajl are, after integration, also completely determined. The stage in the recursion procedure has therefore been reached at which the second set of initial conditions, the a's, must be selected. Since the X.'s for j's other than 1 or 2 have nonzero real parts, (n-2) of the ajl's can be easily obtained from Eq. (4.14) as follows: jJl

84 ~.T v. E T -X.t ai - X. T I e gO[t,a10]dt, ( = 3,4,...,n). (4.16) (i-E ) 0 It is important to note that Eq. (4.16) is equivalent to requiring that cjl(t) have mean value zero in the sense of Cesari [22]. That is, after integration the integral in Eq. (4.13) results in a constant plus a -\.t periodic function of t multiplied by e. If a.l is chosen such that the sum of this constant and al is zero, then cjl(t) has a zero mean value in the sense of Cesari. For example, let gO[saO] = ~ go e O^ 10 " g0 E 0 _=-0 then t -s - g o )t fe g0 [Sa10ds = j [ (-j+JT) Therefore, the mean value of cjl(t) [see Eq. (4.13)] in the sense of Cesari is zero, if X go ajl Vjn ~ (% + T jT which can be shown to be equivalent to Eq. (4.16). The choice of the al's for j = 3,4,...,n is, therefore, quite straightforward. It is necesajl sary merely to choose them such that the mean values of the cjl(t)'s (j = 3,4,...,n) in the sense of Cesari are zero. This requirement can be seen to be related to the requirement in the preceding step of the recursion procedure that the a. = 0 (j = 3,4,5,...,n). Now there remains only the determination of all and a21 to complete the determination of the cll and c21. Once again, due to the reality condition on p(t,), a21 = all. The initial condition, all, is determined

85 analogously to a1 in the previous step, by requiring that it be possible to select a set of cj2(t)'s which satisfy Eq. (4.3). Now consider the general case. Assume that, except for the initial condition alk (a2k = a), the c are known completely. Substituting Eq. (4.10) into the condition, Eq. (4.3), gives, as in the case of Eq. (4.14), X.T X.T T -X.t (~ j - l)aj(k+l) + vjne f gk[t,alk]dt = 0 (4.17) 0 (j = l,2,...,n). Again, gk is a function of only t and alk. It will be shown later in this discussion that gk[t,alk] is a periodic function of t with period T for any alk. Following an argument similar to that which led to Eq. (4.15), alk is determined from either of the following two equations. T -.t f e gk[t,alk] dt = (= 1,2). (4.18) o0 This alk is, in turn, substituted into the integral expression, Eq. (4.10), which determines Cj(k+l). Integration of the resulting integral in Eq. (4.10) determines cj(k+l)(t) up to its initial condition, a.(k+l). As in the case of Eq. (4.16), the aj(k+l) for j = 3,4,...,n, are determined by X T v. e T -X.s j(k+l) - k.T e gk[s,alk]ds (j = 3,4,...n). (4.19) (1-c ) 0 The above expression is again a requirement that the mean values in the sense of Cesari of the cj(k+l)(t)'s be zero. The remaining initial condition, al(k+l) (again, a2(k+l) = al(k+l)), is determined by requiring that it be possible to select a set of c (k+2) which satisfy Eq. (4.3). This requirement results in an equation similar to Eq. (4.18), wherein

86 k is replaced by (k+l). The proof for the uniqueness of the above recursion procedure follows the argument used by Coddington and Levinson ([10], pages 363364). An adequate recursion procedure has now, in principle, been developed. However, for the purposes of this study, the structure of the terms which arise in this procedure must be investigated in more detail so that simpler and more familiar expressions can be obtained. First consider the above recursion procedure applied to the first canonical form. In this case, g = gl. [Note that there is a possibility of confusion between the gl of the first canonical form and the gl of Eq. (4.11).] The definition of gl is repeated below: Xt xt gl(t,cl...,cn) = B(z)Fc 1 +... + n } (2.28) n X.t + B(z)e(t) - Z cjA(z) J. j=i t X t By substituting Eq. (4.6) into F{c 1 +. + cne n} and recognizing nt} that F{c1e +... + ce } is an analytic function of i for IjP sufficiently small, the following expression is obtained: { t t co F E +... + c = z Fk. (4.20) k=0 A A A As Fk appears in this equation it is a function of t,cl,c2k...cnk; however, at the stage in the recursion procedure where Fk is to be employed, k the cjk's are presumably known completely except for the initial condition of Clk i.e, ai.e. Therefore, Fk can be considered as a function of t and a. Morever sinc tt int alk. Moreover, since the Cj. s in each x =C l +... + Cre, r lr

8"7 where 0 < r < k, are determined by the steps in the recursion procedure so that each xr is periodic in t with period T and xk itself is periodic in t with period T for all values of alk, Fk(talk) must be periodic in t with period T for all alk. Therefore, Fk(t,alk) can be expressed as follows: X jX~PTt Fk(talk) = Fk(alk) l (4.21) Q=-oo where F(alk) is considered as a function of only ak. Recalling Eq. (3.11), e(t) can be expressed as e(t)= E (3.11) =-oo Substituting Eqs. (4.21), (3.11) and (4.6) into Eq. (2.28) gives the following expression for the gk of Eq. (4.11): For k = 0 gO(t,a0l) - B(z) FO(al)0 ~=-oo (4.22) o jlPTt 2.t +B(z) E -e c A(Z). =-o j=l In the above equation use is made of the fact that cJo = 0 for J 3,4,...,n. For k > 0 (tja) = B(z) ^ iFk, (alk) T c j(t,alk)A(z)e n \.t - Z cjk(t)A(z)e j=3 It can be seen by inspection of the above two equations that gk(t, alk) (k = 0, 2,...) is periodic in t with period T for all a l.

88 Note that cjk is a function of alk for j = 1,2 and independent of alk for j = 3,4,...,n. This follows directly from the definition of Cjk see Eq. (4.10). Applying the condition, Eq. (4.15), to Eq. (4.22) gives A(j)s) +S nO - B ~JWs) + - = O4 no - B(jw )s 0E (4.24) which determines al10 This equation is, of course, exactly the same as Eq. (3.16). As was previously shown, the initial condition alk is chosen in the course of the recursion procedure such that gk(t,alk) contains no terms of frequency s; therefore, alk is chosen such that Eq. (4.23) contains no terms of frequency 2. However, the making of this choice is not so simple as it was for gO in Eq. (4.22). Since the cjk's which appear in Eq. (4.23) are functions of t rather than constants as are the c.'s, it is necessary to investigate the composition of these ck's in jo k order to discover which of their components contribute to terms of frequency 2 in Eq. (4.23). From Eq. (4.10), the equation for cj is t -Xis jk(t) = ak +Vjnf gkl[Sal(kl)s] s (j = 1,2,...,n). (4.25) Assume that all cjk(t)'s have been completely determined except for the initial condition alk. This means, among other things, that the initial condition ak.1 has been chosen such that gk l[t,a k1] contains no terms 0) of frequency i -~. This in turn means, since X = jo and 2 = -jw and 791, 3W t S ~- t all other Xj's have nonzero real parts, that only Clke and c2ke 1 X.t 8 of the Cjk(t)e's in Eq. (4.23) contain terms of frequency * 2. Further, clk and c2k must be periodic in t with period T. Consequently,..........t 1 Recall that all cjk(t)e J's are periodic with period T. jk

89 o j 1WTt Clk(t) = elk c (4.26) a=-oo and I I J 1pTt C2k(t) 2k (4.27) =-*co Again, only real solutions are of interest for all l; i.e., c2k = Clk; consequently, (clk) = c2 Substituting Eqs. (4.26) and (4.27) into Eq. (4.23) and applying Eq. (4.16) to the result gives the following two conditions which, when 0) satisfied, guarantee that gk contains no terms of frequency + A(Jrs) o ^2si [nk (jO Clk Bw [A(-jws)2k (4.29) s s and [ ^ (-ns)]Ok = [( sj c (4.30) nk B(-Jw)] o _'2k B)c ] l s s F+S F-S where nk = nk = k o k o clk c2k [Note that clk is a simple linear function of alk; see Eq. (4.25)]. How-2s = 0 and +2s ever, re-examination of Eq. (4.25) shows that c 2 = and c = because 1k 2k 0s gk1 does not contain terms with frequency i* -. Therefore, the final system equation reduces to the following familiar form: A(jjws) nk(alk) - B(J) o 0 (4.31) This equation determines alk and, thereby, completes the determination of all the cjk(t)'s and gk. If it should happen that Eq. (4.31) did not have

90 0 0 a solution, then the correct clk would be ck = 0. It is known from the existence of [cj(t,p)] and the uniqueness of this recursion procedure that Eqs. (4.29) and (4.30) have one and only one solution. The cj(k+l)(t)'s are determined by Eq. (4.10) which is repeated below. t -X.s Cj(k+l) = aj(k+l) + Vjn f gk [s alk]ds 0 (4.10) (j = 1,2,...,n), where aj(k+l) for j = 3,4;...,n are determined by X.T Vinc J T -X.s = jn e TJ gk [sa lkIds. (4.19) aj(k+l) =.T (l- 3 ) o That is, for j 1,2 the cj(k+l) s have zero mean value in the sense of Cesari. The initial condition, al(k+l), is determined by A(j s) nk+l(al(k+l)) B( = 0, (432) and a2(k+l) = al(k+l)' At this point g is completely determined, and a further step of the recursion process can be executed. The entire recursion procedure has thus been reduced to a relatively simple process, and it can now be seen that the ordinary describing function relation is just the first in a hierarchy of such relations. Equation (4.32) would be the (k+l)th-member of this hierarchy. Now consider the application of this recursion procedure to the second canonical form. As is shown below, there is little difference

91 here between the first and second canonical forms. In this case, the g in Eq. (4.5) is g2 as defined in Eq. (2.38) and repeated below. klt? t g2(tcl, 1**Cn) = B(z)F cl 1 + + cne n + h(t;wS)} n..t-~ ^-3^(2.38) c.A(z) - A(z)h(t;U ) j=l Again, let klt \nt k F {c +. +cnE + h(t;w)} = E Fk (4.33) k —O and cxa 00 F I~a JITt kt,ak) = (alk) (4.34) =-oo Substituting Eqs. (4.6), (4.33), and (4.34) into Eq. (2.28) gives the following terms in the series expansion, Eq. (4.11), for g2. For k = 0 1 jICUTt 2 A.t g(ta10) = B(z) E F(a)c t-. cjA(z)e - A(z)h(t;w ), (4.35) L=-oOa j =1 where use is made of the fact that c. = 0 for j = 3,4,...,n. Once again, Jo a20 = a10O C10 = a10 and c20 = a20 For k > 0 co BI jl Tt 2 t gk(t,alk) = B(z) ~ Fk(alk)e - e c k(t, )A(z )e l-=oo J=l J'alk n Ajt (4.36) - E cjk(t)A(z)E J=3 It can be seen by inspection of the above two equations that gf(t, alk) (k = 0, 2...) is periodic in t with period T for all alk.

92 Again, cjk is independent of alk for j = 3,4,...,n. Applying the condition, Eq. (4.15), to Eq. (4.35) gives, as expected, A(jw ) no(a10) - B(= (437) F+S where n = -. This equation determines a1o, which in turn completes o clO the determination of the co. s. If o is considered as a variable, n is also a function of w: recall that this was the point of view in Chaps ter III. In a general step of the recursion procedure, the cj(k+l)'s are determined by Eq. (4.10). The initial condition alk is chosen such that gk(t,alk) contains no terms of frequency 2. From Eqs. (4.36) and (4.18) and an argument similar to that which led to Eq. (4.31), the equation which determines alk, is A(jws) nk(alk) B( 0. (4.38) The cj(k+l)'s are, after integration in Eq. (4.10), determined, except for their initial conditions aj(k+l). For j = 3,4,...,n, the aj(k+l)'s are determined by the requirement that c (k+l)i j = 3,4,...,n, have a mean value in the sense of Cesari equal to zero. The reality of p(t,4) requires that a2(k+l) = al(k+l) and al(k+l) is determined in the following step of the recursion procedure by the equation A(jo ) nk+l(a l(k+l)) - B 0. (439) S

CHAPTER V STABILITY CRITERIA 5.1 Introduction At this point the [cj(t,)]'s which correspond to periodic solutions, p(t,), of either Eq. (2.17) or Eq. (2.34) can be determined recursively as a power series in I; i.e., 00 c.(t,) = E ck(t)k, (j= 12.,n) k=o j The relation between a [c.(t,p)] and its corresponding periodic solution, p(t,4), is X1t 2t t p(tp) = cl(t)e + c2(t,4)e +... + cn(t,) for the first canonical form, and Xt xt t P(tp) = Cl(t,)e + c~(t,)e +... + cn(t)e n + h(t;s) for the second canonical form, where, in both cases, lX, X2,..., and Xn are the characteristic roots of the "1 = 0" differential equation. The remaining problem is to determine whether or not this periodic solution, p(t,4), is asymptotically stable as defined in Chapter I. This stability question is treated in this chapter. The first step of the treatment of this question, carried out in the next section of this chapter, is the development of a suitable variational equation, where a variational equation is considered to be suitable if its trivial solution possesses stability characteristics which are similar to those of the periodic solution, p(t,4). This variational 93

94 equation is, as expected (see Section 4 of Chapter I and Appendix A), a,dependent system of n simultaneous first-order linear differential equations with coefficients which are periodic in t with period T and analytic in u for luj sufficiently small. This variational equation can, recalling the outline in Section 4 of Chapter I regarding the stability question, be written in the following form: y= P(t,))y ( - 1), (1020) where y and y are n-component column matrices and P(t,p) is an nxn matrix which is periodic in t with period T and analytic in u for IJ1 sufficiently small. The relation between the stability characteristics of the periodic solution, p(t,4), and those of the trivial solution of the above variational equation are shown to be as follows: (1) if the trivial solution of the variation equation is asymptotically stable, then p(t,~) is asymptotically stable; and (2) if the trivial solution of the variational equation is not stable, then p(t,4) is not stable. Consequently, the key to the stability question is a knowledge of the behavior of the solutions of the variational equation. It is known ([10], pages 78-81) that the general solution of the above variational equation must be of the following form: y(t,I,u) = D(t,4)c R)I, (1.21) where y(t,I,u) is an n-component column matrix, D(t,4) is an nxn matrix which is periodic in t with period T and analytic in 1 for CIl sufficiently small, R(ja) is an nxn matrix which is analytic in ju for }4I{ sufficiently small, and I is an n-component column matrix of initial conditionso The

95 characteristic roots of the matrix R(u) are referred to as the characteristic exponents of the variational equation and are denoted ai(u) (i = l,2,....,n) in this study. If all these characteristic exponents have negative real parts, then the trivial solution of the variational equation is asymptotically stable; and if one or more of these characteristic exponents has a positive real part, then the trivial solution of the variational equation is not stable. After the variational equation has been developed in Section 2, the remainder of this chapter is devoted to determining the characteristic exponents, a.i(), recursively as power series in A. The technique employed for this determination is an adaptation of one utilized by Moulton ([16], pages 331-348) in celestial mechanics. The important feature of this technique is the introduction of n new coordinate transformations, each of which depends upon a separate, as-yet-unknown, characteristic exponent, ai (I). These transformations are of the form a. ()t y = U, (i = 1,2,...,n) (1.24) i (A)t where e is a scalar and U is an n-component column matrix. These transformations are carried out in Section 3 of this chapter. As is shown subsequently in this chapter, when the variational equation is expressed in terms of any one of the n new possible sets of variables, U [one set for each ai(u)], it is known to possess a solution which is periodic in t with period T. As an illustration consider the following first-order linear system with periodic coefficients: y = [cz() + P(t,p)]y, where P(t,4), a scalar here, is periodic in t with period T, of mean

96 value zero, and analytic in A for \I1 sufficiently small; and a(4) is analytic in, for 1j1 sufficiently small. The general solution of this equation is t oa(4)t f P(s,4)ds y(t,y(O),4) = y(O)e e t f P(sp)ds Note that e is periodic in t with period T. Now, introducing the coordinate transformation y = u()tu, the differential equation becomes u {[a() + P(t,)] - a(c)}u. This new equation has the following i-dependent periodic solution: t f P(s,))ds u(t,u(O),) = u(O)e In the case of higher-order systems, a similar argument can be followed. It is necessary to note only that for n selected, linearly independent, initial condition vectors, Ii (i = 1,2,...,n), the solutions of the variational equation reduce to a periodic column matrix multiplied by an a.(Wrt appropriate scalar, E, where each ai(4) (i = l,2,..,n) is associated with one of the n selected initial condition vectors, Ii. Therefore, in place of the one variational equation, the n coordinate transformations result in n new variational equations of the form U = [P(t,) - a()EJU (i = 1,2,...,n) (1.25) dU where U 5 d E is the identity matrix, and each of these n systems is

97 known to possess a periodic solution. Moulton's technique then involves expanding all the quantities in Eq. (1.25) in power series in 4, carrying out the indicated multiplications of power series, and equating the coefficients of equal powers of, on either side of the equality sign to one another. The result of this equating is a hierarchy of differential equation systems which can be solved recursively for both the characteristic exponent, Ci(~) = ^io +l + i0 L +.. (i = and the known-to-exist periodic solution, 2 U.(t4) = ui.o(t) + iJl(t) +. u(t) +... (ij = The recursion procedure employed is similar to that already used in Chapter IV to determine [cj(t,)]. As in Chapter IV, the key to carrying out each step in this recursion procedure is the requirement that each term, u. k(t) (i,j = l,...,n; k = 0,1,2,...), satisfy a periodicity condition. A typical step in the recursion procedure employed in this chapter, say, for the determination of aik. follows the ensuing pattern: (a) At the outset of this step a system of n inhomogeneous linear differential equations involving the uijk(t)'s (j = l,2,...,n) is known. This system is the appropriate member of the above-mentioned hierarchy. The forcing functions of this system are made up of the, as-yet-undetermined, aik and the uijo(t),..o,uij(k. )(t) (j = l,2,...,n) which have, with the possible exception of one of the initial conditions for the

98 uij(kl)(t)ts, been cetermined in the previous steps of the recursion procedure. (b) In the first part of this step of the recursion procedure, a. and the undetermined initial condition, I. remaining from the previous step, if it is present, are selected so that the solutions of the system of inhomogeneous differential equations can have a periodic solution with period T. This selection determines the a. and completes, with the determination of the initial condition, the determination of the u. (t)'s. ij(k-l) (c) In the second part of this step of the recursion procedure, the initial conditions, uijk(O) (j = 1,2,...,n),are selected so that the uijk(t)'s (j = l,2,...,n) are periodic in t with period T. At this point it can happen (for i = 1 or i = 2) that one initial condition can be chosen arbitrarily without destroying the periodicity of the u jk(t)'s. This undetermined initial condition must then be determined in the following step of the recursion procedure. The foregoing recursive development is carried out in Sections 4, 5, 6, and 7 of this chapter. The organization of these sections is as follows:

99 (a) In Section 4, all the first approximations,.io (i = l,2,...,n), to the n characteristic exponents, cz.(l), and all the u. i(t)'s (i,j = l,2,...,n), ex1 lJo cept for ul20(0) and u21(0)), are determined. Note that these determinations are the first steps in the application of the recursion procedure to n different systems of differential equations; i.e., there is a separate system associated with the determination of each a,(,). (b) In Section 5, all the second terms, ail (i = 12,...,n), in the series for the n characteristic exponents are determined. At the same time, the undetermined initial conditions, u120(O) and u 210() from the preceding step in the recursion procedure and all the ui j(t)'s (i,j = 1,2,...,n), except for the initial conditions, ul2l(0) and u21l(0), are determined. Note that these determinations are the second steps in the application of the present recursion procedure to n different systems of differential equations. (c) In Section 6, all the third terms, ai2 (i = 12,...,n) in the series for the n characteristic exponents are determined. At the same time, the undetermined initial conditions, u121(0) and u211(0), from the preceding step in the recursion procedure and all the ui. (t)'s (i,j = 1,2,...,n), except for the initial conditions, u122(0) and u212(0), are determined.

100 (d) In Section 7, the procedure to be followed in the determination of the cik'S and uijk(t)'s for k > 2 is presented. At each step in the recursion procedure relations which simplify the determination of approximate stability criteria are pointed out. These relations involve previously defined variables such as *2s n, S, SO A(jw ), and B(jos). Throughout the following discussion of the stability question, special or degenerate cases —such as those with repeated roots —are ignored. Although such cases have mathematical importance and even, sometimes, physical importance, for the purposes of this study they can be neglected. This can be appreciated by noting that in such cases a slight change of a component or parameter in the physical system often removes the special or degenerate condition. 5.2 Development of a Suitable Variational Equation The solutions [c (t,p~)] which are of interest are known to possess certain properties. First, [cj(t,p)] is a solution of the following differential equation dr1 dr -Xt dt Vln 2 n) dr2'-2t d = V2n g(t;rl, 2n.,r dr -t t dn _ ^e n g(t;rlr}-. ), where g is gl [Eq. (2.28)] in the first canonical form and g2 [Eq. (2.38)]

101 in the second canonical form. Secondly, the components, cj(t,) (j = l,...,n), of a given solution [cj(t,)] have been developed by the procedure of Chapter IV in the following form: co c.(t,) - c k(t)4 ( = 1,2,...n) (4.6) k=O k Thirdly, each cjk(t) satisfies the following "periodicity" condition: X.T J - 1,2,*...n cjk(T) e c jk(0) = { (43) k = 02k..= This condition on the cjk(t)'s implies that [cj(t,)J also satisfies this condition. Fourthly, the [cj(t,i)J's are related to the continuable periodic solutions p(t,4) by one of the following two transformations. First Canonical Form: |,(ttA S.. e n 1 p~l~t' q (2.19) P(t,4) C.. e Cl(tp) L(n-l,(n-l)i ) (l)n lt n-l) %nt Second Canonical Form: p(tA) h(t' ) hiP(^-^)' X L j 0 g Lh(nt;w ) p~~~tCI) I I II~~~~~~~~~~~

102 where [h(n-1)(t (n' l)-] and h(t;w), which is defined in Eq. (2.33), and its first (n-l) derivatives with respect to t are continuous and periodic with period T. The function p(t,,) is, of course, a periodic solution of either Eq. (2.17) or Eq. (2.34). Now, the actual stability problem is related to the continuable periodic solutions, p(t,p), and not to the corresponding [cj(t,4)]. It can be seen that these two points of view can give different answers to a stability investigation by considering only the following relation: Xlt X t p(tl) = cl(t,) 1 +... + cn(tP,)e. (5.1) Let X. be one of the A's in the above expression. Further, assume that J X. is real and negative. Then in the process of determining whether or not the periodic solution p(t,I) is asymptotically stable, expressions of the following form arise: qp(t,IjA) - p(t,4) -.0 as t - o, (5.1A) where cp(t,I,I) is a solution whose initial conditions neighbor those of p(t,p). But, from Eq. (5.1), this implies that X.t e J tpj(t,Ap) - cj(tC)1 - 0 as t -, (5.2) where pj(t,Ap) is the jth component of a solution whose initial conditions neighbor those of cj(ta). It is immediately apparent that the above condition can be satisfied without even having Jpj(tAp) - c.(t,p){ bounded as t -, oo. It is necessary only that

103 (lM.l-5)t lpj(t,A,) - cj(t,)l < K (3) where K and 5 are real positive constants. Therefore, asymptotic stability for a solution [cj(t,4)], assuming the obvious extension of the definition of asymptotic stability to aperiodic solution such as the [cj(t,)]3, is not a necessary condition for the asymptotic stability of the corresponding p(t,4). Moreover, it will not be a sufficient condition if some of the j's have positive real parts. In order to investigate the stability of the continuable periodic solution, p(t,)), the canonical forms must be transformed to new forms whose solutions exhibit the same stability qualities as the p(t,4)'s. One obvious choice for this transformation is just to re-express the canonical forms in the original (xx,...,x(n-l)) system; however, this choice is not so convenient as the system which results from the following coordinate transformation: x1 ~ 0 ~ 0 r X1 | F 1.. 0 r^ X2t ~~2 2~ (5.4) X t x 0 0... rn n n This new set of variables (x,..,xn) is related to the (x,,...,x(nl)) variables by the following transformation:

104 (5.5) x 1 1... 1 x ~x \l X2... n X *~ 2 2 2 1 2 n 3 0 9 9 0 00 0 0 0 0 0 0 I3 * 0 n1)(n-l (-) (n-l) (nli) 2 n f X Since the XI.s are distinct, the above transformation is nonsingular. 3 The new variables are a linear, time-independent, combination of the old ones; consequently, the stability qualities of the periodic solutions are invariant under this transformation, and the solutions p(t,4) remain periodic in t with period T if they are transformed by Eq. (55). Introducing the coordinate transformation, Eq. (5.4), into either of the canonical forms results in the following system of the differential equations for the characterization of the nonlinear feedback system. (5.6) -\lt -X t -Xt - t L X0 0...0 X1 vnng(t,xl,...XE n ) where the g is gl in the first canonical form and g2 in the second. Note that g is now a function of (t,xe,.., ) instead of (t,rl,1.,rX)n that g is now a function of (t,x1,...,x) instead of (~~.~)

105 The stability problem now becomes the problem of determining whether or not the periodic solutions of Eq. (5.6) are asymptotically stable. This is just the problem which is considered in Appendix A; therefore, the methods given in that appendix can immediately be applied. The variational equation is (5.7) ag t g nt 1 Yo 0... 0.. o y1 I r=c n r=c Y2 2 2 0... O.. Xt g tIn ~ Yn~r Vnn - V)r1''nn-t L Yn 1r=c n r=c where r = c implies {r1 = c(t,),...,r = c (tl)}. Note that in the time-cependent coefficient matrix the following identity has been employed: x m 8g ~EA (j = 1,2,...,n). (5.6) j P(t,4) j r=c This has been done so that the final stability criteria can be expressed in terms of the same variables, c.j, as the solutions [cj(t,u)]. From the definition of either gl, in Eq. (2.28), or g2, in Eq. (2.38), it can be seen that the s e J's are periodic in t with period T and analytic r=c in p for 141 sufficiently small. Floquet [17] has shown that a homogeneous linear differential equation with periodic coefficients has a fundamental solution matrix of the form, M(t,^) - D(t,^cl)E (5~9)

106 where D(t,p) is periodic in t with period T, and R(u) is independent of t. The characteristic exponents are the roots of the following nthdegree polynomial in a: det [R(4) - aE] = 0, (5.10) where E is the identity matrix and the characteristic exponents are denoted ao(n),..an (A). The so-called characteristic multipliers are 1" n(c )T an(G)T obtained as follows: e,)..., 5.3 Transformation of the Variational Equation In this section the variational equation, Eq. (5.7), will be transformed with the aid of n time-dependent transformations into n separate systems of linear differential equations with periodic coefficients. Each of these new linear systems is known to have a solution, denoted [uij(t,)l3 (i,j = l,.o.,n), which is periodic in t with period T and analytic in ~. In Section 4 of this chapter a recursion procedure similar to that used in Chapter IV for obtaining [c.(t,4)] is developed which gives each of these latest periodic solutions, [u. (t,4)], in a power series in A, and concomitantly results in developing each of the ai(O) in a power series in i. Owing to the composition of the fundamental solution matrix M(t,0), as expressed in Eq. (5.9), it is possible, after assuming the characteristic multipliers to be distinct, to select n linearly independent solutions of the variational equation which have the following form:

107 yi(t,) ui(t,) Yi2(t ()t i2(t) = ~ (i = 1,2,...,n), (5.11) Yin(t,ll) Uin(t, ) where a.i() is one of the characteristic exponents, and each u ij(t,) (i,j = l,...,n) is a periodic function of t with period T for |14 sufficiently small. If the n coordinate transformations, 1 U1 Y2 ca (i)t U2 = e (i= 1,2,...n), (5.12) y U n n are introduced into the variational equation, the following n new variational systems are obtained (one for each value of the integer index i, i = 1,2,...,n): 1 + [oxi() - kl]u = uVln[ rl(t,4)ul+ o..+ r(t,4)un] U2 + [C.(G) - k]u2 = LV2n[ rF(t,9)ul+...+n (t,4)u] 2 z1 22 2n 1 1n (5.13) \* *0 0 O * * 0 * * 0 0 0 0 * O * 0 * * * * 0 *nn Un [ai(4) - An]un = vn[ rl(t,~)ul+...+ n(t,)Un],

108 where r.(t,) = | (j =,2,...,n).. It is known from Eq. aj r=c (5.11) that each of the above systems [i.e., a different system for each of the n a.i()'s] has a periodic solution with period T, [uij(t,4)]. Moreover, because the coefficients in the variational equation, Eq. (5-7), are analytic in i in the same domain in which the solutions [cj(t,4)] are analytic in i, the characteristic exponents, ai(j), are analytic in 4 within this domain with the possible exception of a finite number of branch points. Similarly, the general solution, [uj(t,u.(0),4)],is analytic in A within the same domain. In particular, the periodic solution, [u,.(t,>)], 13 of Eq. (5.13) with period T, which is known to exist, is analytic in u for 4 within the above domain. For a further discussion of these last two points see Cesari {122], page 59} and Moulton {[161, pages 331-352}. 5.4 Determination of the First Terms, aio (i = 1,2,...,n), in the Power Series Expansion for the Characteristic Exponents, Ci()) (i = 1,2,...,n) Since the characteristic exponents, ai(4) (i = l,...,n), are analytic within the same domain as the coefficients of the variational equation, Eq. (5.7), with the possible exception of a finite number of branch points, they can be expanded in power series as follows: (4) =.ikl.n (5.1^i) k=O where 4 is within the domain of analyticity for cai(u). It will be assumed throughout this discussion that, with the exception of one case discussed later, all characteristic multipliers are distinct for Ilu < 1. This assumption removes the branch points as a source of difficulty. Similarly, the solutions of Eq. (5.13) which are periodic in t with period T can be represented as follows:

109 o0log uij(t,4) = Ui j(t) (i,j = 1,2,..n), (515) k=O ijk where the index i indicates that periodic solution associated with a.i(u)' Since the uij(t,4) are periodic in t with period T over some ~interval of nonzero length, the individual u. (t)'s are also periodic in t with period T. The rj(t,p)'s in Eq. (5.13) can also be represented in series form; i.e., co r.((t) = r (j =,2,. n). (5.16) k=O jk From Eq. (5.13), ag - t r.(t,) = e (j = 1,2...n) 3;j r=c From Eqs. (2.28) and (2.38), g is either { t \nt l(t,r,...,rn) = B(z)F{rl 1 +... +rn } (2.28) n A.t + B(z)e(t) - Z r.A(z)e 0 j=-l or kit Xnt g2(t,rl,...n ) = B(z)F{rle +...+rnE + h(t,js)} n X.t (2.38) - r.A(z)c -A(z)h(t(e )A j=en, since (see Appendix ) -lThen, since (see Appendix D)

110 k it knt [ dF3 t. )B(z)F{rC +.+r C - B(z)[ld (D1) aW n ]r=c dxx=p(t,4) 6 t X t Xt r B(z){re +...+rn + h(t,W) = B(z)[, (D13) r=c x=p(t,) where (j = 1,2,...,n), and B(z) I XBt(z + x.)[ B(z )[dx|e ]= iB(z x)[d1 ]. (5.17) x==p(t),) x=p(t,ui) the expression for r.(t,4) in either of the canonical forms is r(t^,) -= B(z + )[ dx A(.) (5.18) 3 3Xj x=p(tp) Let dF 2 dx = S0(t) + SSl(t) + S2(t)+..., (5.19) x —p(t, i1) where the S (t),S(t),... are periodic functions of t with period T. The expressions for the r.k(t) are, therefore, as follows: Jk For k = 0, rIo(t) = B(z + j)So(t) - A(=,2,...,n); (5.20) For k > 0, rk(t) = B(z + Xj)Sk(t) { (5.21) 3 k = 1,2,...,n Now, the object of this chapter is to develop a procedure for determining the aik's [see Eq. (5.14)] recursively. Not too surprisingly, this involves determining the uiJk(t)'s simultaneously with the aiks.

ill The initial step of this development is to substitute Eqs. (5.14), (5.15), and (5.16) into Eq. (5.13). If the indicated multiplications of power series are carried out and the coefficients of like powers of A on each side of the equation are equated to one another, a sequence of linear inhomogeneous systems of equations is obtained. The system of equations which arises from the i. terms is il0+ (+ io - l) ilO ai20 + (aio 2)Ui20 - (5-22) a. +(c^ -X)u = 0,!ino (Rio n )uino where (i = 1,2...,n). The system of equations associated with the i terms is Uill + ( io - )uill = - a110 +v ln{ (t)uil0O+.0+ no(t)uino} i21+ (io 2)ui21 - il u20 + v2n{ rl(il no ()ino} (5.23) inl + (Cio - A) = ilino nn rO)ilo+ + rno (t)ino} where (i = 1,2,...,n). Note that the right-hand side of Eq. (5.23) is, assuming that a first step in the recursion has already determined the appropriate u.. (t)'s, a known forcing function. Similar systems of 34 equations are associated with i, ~...~ terms.

112 The first step in the recursion procedure is to select the a.'s and the initial conditions for the u.. (t)'s in Eq. (5.22). The key to this selection is the requirement that all u..k(t)'s be periodic with period T. Therefore, one possible choice for the a.'s is 10'l =Xo a2 =o =.X s(5.24) 1 x2,',n0 n 2 As far as the periodicity requirement is concerned, the choice of a. ios is not unique. Any integer multiple of joT can be added to any of the ai's in Eq. (5.24) without destroying the required periodicity of the particular solutions [u~. (t)]. This ambiguity is just what is expected from Floquet Theory, which states that the characteristic exponents can be determined only up to integer multiples of JOT' On the other hand, the choice in Eq. (5.24) does have a simple relation with the "u = 0" differential equation of Eq. (5.7). That is, as A approaches zero the characteristic exponent, ai(i), approaches the characteristic root, \i. After a particular aio has been substituted into Eq. (5.22) the initial conditions must be chosen so that the u.. (t)'s are all periodic with period T. Since Eq. (5.13) is a linear system, each set of initial conditions corresponding to a periodic solution [one set for each 0i(1)3 can be determined only up to an arbitrary factor. Therefore, one component of each set of initial conditions, say uii(O,), can always be chosen arbitrarily and each of the others adjusted so that the ratio between it and uii(O,.) is correct. Furthermore, uii(O,) can be chosen to be independent of A so that uii(O,) = uii (0) and uiik(O) = 0 for k = 1,2,... For the remainder of this chapter, the initial conditions are selected according to the preceding pattern. Therefore, the periodic solutions [u.io(t)] are as follows:

113 For C10 = l ullO(t) =- u1(O) (an arbitrary constant) -j2 t u20 (t) u20 () e (u"120(0), undetermined) (5.25) u jo(t) -0 (j = 3,4,,n); Uljo For a20 2, +j2u t 210(t) = u210(O ) (u210(0), undetermined) u220(t) u20(0) (an arbitrary constant) (5.26) u2jo(t) ( 3, n) 2jo For ai = i. (i = 3,4,.,n), Uilo(t) 0 u20(t) - (5.27) ui0i(t) ui o(0) (an arbitrary constant) Ui (t) 0.

114 It is immediately apparent from Eqs. (5.25) and (5.26) that, although Eq. (5.22) is sufficient for the determination of the a.'s it is not sufficient for the complete determination of the [uio (t)]. The difficulty is that ul20 and u210(0) are undetermined. Reminiscent of the recursion procedure in Chapter IV, u12 a(0) and u21(0) must be determined in the next step of the recursion procedure. Equations (5.25) and (5.26) demonstrate also that oal and a20 are special cases and require slightly-more-involved treatment than the other aiots (i = 3,~..,n). The essential difference between these two groups of io s is that (a10 - a20) is an integer multiple of joT. Since ca () -, o and a2(j) -* a20 as A ->0, and the characteristic multipliers associated with these two characteristic exponents are e and 2(g)T *jw T e, as p - 0 the two characteristic multipliers approach 1 (e = 1). Therefore, the significance of (alo' a20) being an integer multiple of jWT is that the two associated characteristic multipliers coalesce as u - O0, even though the characteristic exponents remain distinctl. Recall that it was assumed previously that the characteristic multipliers were distinct. It is now apparent that this assumption must be modified to allow the above coalescence. Therefore, it will be assumed that all the characteristic mulaol(1)T a2(4)T tipliers are distinct with the exception of e and e at u = 0. Moulton ([161, page 335) shows that in this situation al(G) and a2(l) are still analytic in the neighborhood of ~ = 0 but that the recursion 1 The word "distinct" here is somewhat ambiguous. It is known from Floquet Theory that the characteristic exponents can be determined only up to integer multiples of jow; therefore, in that sense, ah(u) and a (u) are not distinct at u = O because jo = -jow (mod j ~ On the other hand, in the ordinary meaning of the word, a (,u) ana a2(j) are distinct at, = 0. It is just this ambiguity that 7ery often makes it easier to discuss the characteristic multipliers, which are unchanged no matter how many integer multiples of joC are added to the characteristic exponents, than the characteristic exponents.

115 procedure based upon a periodicity condition for the [uij(t,)]'s leads to two possible analytic functions cl(4) such that al(>) -,io o as -, O and, similarly, there are two possible analytic functions a2(G) such that a2 () - a20 as 4 - O. However, these four analytic functions are not unrelated. The difference between either of the two Oal()'s and a corresponding one of the two a 2(4)'s is just an integer multiple of joT. Therefore, these four functions characterize only two i-dependent characteristic multipliers. One of the al(4)'s is associated with one of the two 4-dependent characteristic multipliers, and the other al(4) is associated with the remaining i-dependent characteristic multiplier. Consequently, in a recursion procedure it is necessary to guarantee that an al(4) and an a 2() are developed which correspond to distinct characteristic multipliers. This guarantee exists as long as the ca1() and a2(4) are selected from the four functions so that they do not differ by an integer multiple of jWT' 5.5 Determination of the Coefficients, a i(i = l,2,...,n), of A in the Power Series Expansion for the Characteristic Exponents, ai(4) (i = 1,2,...n) The uij(t)'s (i,j = 1,2,...,n) satisfy Eq. (5.23) and are periodic in t with period T. The terms which appear on the right-hand side of Eq. (5.23) are all, except for ail (i = 1,2,...,n), U120(O), and u21(0)(, known from the first step in the recursion procedure. These known terms are given in Eqs. (5.20), (5.25), (5'.26), and (5.27). For each of the n systems of equations encompassed by Eq. (5.23) (one system for each choice of index i) the appropriate, undetermined quantities on the righthand side must be chosen so that the uijl(t)'s (i,j = 1,2,...,n) can be periodic in t with period T. Thus, for i = 1, o11 and ul12(0) must be

116 chosen so that the u1jl(t)'s can be periodic with period T. For i = 2, a21 and u210(0) must be chosen so that the u2jl(t)'s can be periodic with period T. And, for i = 3,4,...,n, the il must be chosen so that the ui (t)'s can be periodic with period T. For a1: Consider Eq. (5.23) for i = 1 and with Eq. (5.25) substituted for the uljo(t)'s. Recall that c0 = 1, Xl = jo. (5.28) Ulll = -o11 U110() + Vln{ 10(t)ull0() -j2w t + F20(t)u120O(o) } -j2w t U121+ j 2w U121 = -' lU120 (0) + 2n(r10(t)ul10(0 ) -j2w t +~20 (120 (0 } + j 20(t)1+ 0(O)E S Ui31+ (Js i 3) 1313n rFl(t)110(0 -j2w t +r20(t)U 120o( Ulni + (OWs n)Ulnl + vn{ p10(t)u110(O) It can be seen from the first euation of the above system that (t) It can be seen from the first equation of the above system that u111(t) is periodic with period T if the mean value of the right-hand side of this first equation is zero. From the second equation, u121(t) is periodic with period T if the right-hand side of this second equation contains no term

117 of angular frequency -2w. These two conditions are satisfied if a11 and s u12(0) are chosen so that they satisfy the following pair of simultaneous equations vl T -nv T -j2w st { n I r (t)dt - u110o() + T r20()}dt 120 O O O (5.29) v2n T +2t 2n T { f T F J(t) dt}U110 () +- T o r20(t)dt - a11}u120(0) 0. o 0 This system has a nontrivial solution if the following determinant is zero. {l- a}{ T v2n T {T rlO (t)d- -1 { T I 20(t)dt -al } 0 0 (5 30).v T -j2s t dt v V2n T +j2w = t -{ f 1720(t)e d T r {10( e dt} O O The above equation can be reducedwith the aid of Eq. (5o20), to the following form: 2! -. ~nl{^v [B(Js ) o - ^(j)] + v B( ) A(-J ^ )]} + VlnV2n[B(jws)So - A(jws)][B(-Jcs)~ - A( -jIs)] (5.31) Vl nV2n B(j)B(-j s s)So.S 2 = o The two roots of the preceding equation, say, o1) and ()' are an evidence of the previously-discussed fact that two ae(i) functions are possible. Moreover, in the determination of a21 it is shown that a21 also

118 satisfies Eq. (5.31); therefore, there appears to be a certain redundancy in the determination of all and a21' However, recall that if the same root of Eq. (5.31) is used for both oll and a21, the first two terms of the series for ol(0) and a2(G) are jws + iall and'J)s + ll, respectively; but in such a case the characteristic exponents, as presaged by these first two terms, are separated by an integer multiple of juT (os = pT}, p an integer), which indicates that the two series are associated with the same characteristic multiplier. The apparent redundancy can, therefore, be resolved by choosing oll equal to either one of the roots of Eq. (5.21), and choosing a21 equal to the remaining root of Eq. (5.31). Consequently, 11C21 = VlnV2n[B(jw s)S -A(j' )][B(-Jjw)S - A(-j )] (5.32) Vv B(J )B(-j )s+2s -2s VlnV2nB( js)B(-js)so So and L11 +C21 = vn[B(jws)SO - A(js)] + v2n[B(-s)So - A(-js] (5 3) A comparison of Eq. (5.32) with either Eq. (3.21) or Eqo (3.31) reveals that the right-hand side of Eq. (5.32) is the expression which appears in the Jacobian condition. Therefore, the Jacobian condition yields information regarding the roots of Eq. (5.31). Since in many nonlinear feedback systems the answer to the yes-or-no part of the stability question can be obtained from only the first two terms of the power series representation for the characteristic exponents, ci(u) (i = l,...,n), this relation between the roots of Eq. (5.31) and the Jacobian condition is important. In particular, from the discussion regarding Eqs. (3.40) and

119 (3.46), Eq. (5.32) leads to the following observations: (a) If the Jacobian expression is less than zero, then one of the pair, c1, and a21, is real and negative and the other is real and positive. (b) If the Jacobian expression is greater than zero, then both all1 and c21 are real numbers with the same sign or a complex-conjugate pair located anywhere in the complex plane. (c) If the Jacobian expression is zero, then one (repeated roots are ruled out) of the pair, a1 and a21, is zero. Assuming that io + aCil (i = l,..,n) is a satisfacdw tory approximation to a.i(u) and recalling that da = 0 if the Jacobian expression vanishes, a simple stability do criterion results; that is, d- = 0 at the points in da (a,O, W )-space which form a boundary between the region in which the product of the roots of Eq. (5.31) is positive (all121 > 0) and the region in which their product is negative (all21 < 0). Note that the roots of Eq. (5.31) can be in the right-hand plane even though the Jacobian expression is positive. (d) Unfortunately, it is not true that the sign of the dw Jacobian expression is determined by the sign of dda The van der Pol oscillator under the influence of a sinusoidal forcing function is a classic counterexample. It has stable solutions with both positive

120 dw and negative da (e) If the Jacobian expression is positive, then the other stability criterion associated with Eq. (5.31) is Re { vl[B(jws)S - A( is)]} < 0 (534) This criterion follows from the fact that each of the two terms which make up the coefficients of Ocll in Eq. (5.31) is the complex conjugate of the other. It is shown later that this criterion is related to the stability criteria developed for the ci's (i = 3,4,.*,n). Note that the coefficients in Eq. (5.31) depend only on the generating solution, [Cjo], and on none of the higher-order approximations, [Cjk] (k = 1,2,...). Furthermore, it is not necessary to complete the determination of the uljl(t)'s (j = l,...,n) in order to find a11 Therefore, the first two zik's of each series can be obtained with relative ease. On the other hand, if the first two terms of the power series representations for the aG(u)'s (i = l,...,n) do not provide sufficient information regarding stability and it is desired to determine more terms in these power series, the determination of the periodic solutions, Uljl(t) (j = l,...,n), of Eq. (5.28) must be completed. This determination is carried out in the remainder of this section. Assuming that all has been chosen as one of the roots of Eq. (5.31), it can be substituted into Eqo (5.29) to determine u120(0). Utilizing, in Eq. (5.29), the definitions for rj(t) andr'20(t) from Eq.

121 (5.20) gives the following two equivalent relations for the determination of u120(O): ua120() {11 - s[B(j)) A(j) u120(~) ) A ) ] }+ 110.. v B(. )S... in s o (5.35) U120(o) v2nS20B(-jS ) u110(O) 11 2n[B(Jw5)So - A(ij )]} This relation completes the determination of all terms on the right-hand side of Eq. (5.28). The functions, uljl(t) (j = l,...,n), can then be determined by straightforward integration of the n first-order, linear, inhomogeneous differential equations which comprise Eq. (5.28). It is now necessary only to choose the initial conditions of the uljl(t)'s so that they are periodic in t with period T and so that u111(O) = 0. The results of this integration are as follows: t u111(t) = U110(o)f [-11 + v ln [B(z + j )S0(S) - A(jws)] "O~~~~~ (^~~(5.36) u120 (0) -ja s + u-O i s [B(z - js )S0(s) - A(-j S) }]ds, -j20t t -j2 (t-s) (120(0) - s Ul21(t) = U121(0)e + Ullo c L 1 UlO() E O 110 + V2n[B( + js )S0(s) - A(jws)] (537) Ul20(0) -j2ws s + 2o C [B(z - j )So.(S) - -Js)}] ds, and

122 (%.-jw )t t (.j-jo )(t-s) Uljl(t) = Uljl(O)E s + U11o(O) j s 0 [vjn {[B(z + jws)S (s) - A(jo )] (5.38) -j~ ss u1(0) + U - u[B(0 - ( s)S0(s) - A(-j s)}] ds; where (j = 3,4,...,n) in Eq. (5.38). Since Uljl(t) in Eq. (5.38) must be a periodic function of t with period T, the initial condition, uljl(O), is chosen as follows: Un11(O)v jn T (jws-%i )s uijiUllO) VIn ('j)s [B(z + j s)So(S ) - A(s)] ~(e ~ J~ -.^ ^1) ~~ (5.39) -j2a s u20 (0) + Ec 10) [B(z - js)S (s) - A(-j)]} ds. 110O Except for the initial condition of u121(t), the above equations completely determine the u ll(t)'s (j = l,2,...,n). The u121(t) is periodic in t with period T for any u121(0); therefore, a periodicity condition on u121(t) is not sufficient for the determination of u121(0). As in the case of the initial condition u120(0), the initial condition u121(0) must be determined in the next step of the recursion procedure for cz(4). For a2: Consider Eq. (5.23) for i = 2 and with Eq. (5.26) substituted for the u2jo(t)'s (j =,2,...,n). Recall that a20 =,2' k2 = -j.s See Eq. (5.40) on the following page.

123 +j2a t +ju) t U211 - j2 s U211 = - o21U210(0)E S { r10(t)210() + r0)20(t)u0()} +j^w S 2=21 = -21u220(0) + Vn2 r10(t)u210() S + r2o(t)u220(0)} +j2m t 231 + ( - 3)U231= + V3n r10(t)u210(o)E (5.40) + r2(t)u220(o)} + 20 (t)U220(o)} +j2 All the terms on the right-hand sides of the equations in Eq. (5.40) are known, except for a21 and u210(O). As in the case of the determination u210(0)] must be chosen so that the right-hand side of the first equation in Eq. (5.40) contains no term of angular frequency +2w and the rights hand side of the second equation has mean value equal to zero. These two conditions lead to the following pair of equations: {vln B( s) o S220(0) +{Vln[B(iws)0 A(js)]- a21}210(O) = 0; (5.41) {V2n[B( jws)s - A(-jws)] - a23}U220(0) +{V2nB( -js)So2 }u210(0) =.

124 This system has a nontrivial solution only if its determinant is zero. The equation resulting is, as has already been mentioned, the same as the quadratic equation, Eq. (5.31), for Ca1l. Therefore, the previous statement that all and a21 satisfy the same quadratic equation is substantiated As soon as the c21 is determined, the undetermined initial condition, u210(0), can be determined from Eq. (5.41) as follows: U210 () vlnB(jos )S u220 0 U220 {1 - lvl [B(j )S0 - A(w } )] (5.42) or 210(~) { 21 v2n[B(-j s)SO - A(-jw )]} U220 () v2nB(-jw )So Equation (5.42) completes the determination of all the terms on the righthand sides of the equations in Eq. (5.40). The functions u2jl(t)'s (j = l,...,n) can, therefore, be obtained from a straightforward integration of the differential equations in Eq. (5.40). The results of these integations are given below. j2o t t j2w (t-s) U210() +ji s s u210(O) +jSusS u211(t) = u211(O)E + u ) u~2200220( + v{ 210 s [B(z + jW )S (s) - A(jw5) (5.43) + [B(z - js)S(s) - A(-js)] ds [Recall that uii(O) (i = l,2,...,n) has been chosen equal to zero.]

125 u22(t) = u20() f[- a21 + v{ U2 0() [B(z + jW )S(s)- A( )] O 220 (5.44) + [B(z - jws)So(s) - A(-jws)]}] ds (kj+jas)t u2jl(t) = 2jl(0)E t (kj+jWs) (t-s) u210(O) j'2sS + 22(0) e [ v 20(0) C [B(z + jws)S (s) (5.45) 0 220 - A(j)] + [B(z - j s)So(s) - A(-j s) }] ds where j = 3,4,...,n. Since u2jl(t) must be periodic in t with period T, the initial condition, u2jl(O) (j = 3,...,n), is chosen as follows: U2jl (0): on ~O 0ms s u2(0) v U220(o) T -(jw5+X)s {U210(0)e [B(z +J s u 1(0)".T Uo(O)) - [B(z + j- )S (3) j 220 (5.46) A(js)] + [B(z - jWs)So(s) - A(-jws)]}ds. The above equation completes, except for the initial condition u11(0), the determination of the u2jl(t)'s (j = l,...,n). The u211(0) must be determined in the next step of the recursion procedure for'2(6)* For il (i = 34,...,n): Consider Eq. (5.23) for (i = 3,4,...,n) and with Eq. (5.27) substituted for the uijo(t)'s (j = l,...,n).

126 Uill + (xi - l)Uill = Vlnuiio(O)rio(t) i21 + (hi 2)ui2 = 2 nuiio ()ro(t) (5.47) Uiil = {Cil + inio(t)}uiio(0) ~Ii~ ~ ~ o ~ tl ~ ~ I ii I (o ~ s ~}\ lo ~ * e e * * o0 * 0 * S 0 5 0 * 0 * 0 * * v v * * a 0 Uini + (xi - n)UinI - VnnUi () o(t) All terms on the right-hand sides of the equations in Eq. (5.47) except for al are known. Cil and the initial conditions for the u i(t)'s (j = l,2,...,n; i = 3,4,...,n) must be chosen so that the ui j(t)'s are periodic with period T. It can be seen that uiil(t) is periodic with period T if the right-hand side of the equation in Eq. (5.47) involving uiil(t) and ail has mean value zero. This requirement is satisfied if v. T in l = T fio(t)dt (i = 3,4,..,n). (5.48) 0 Using Eq. (5.20), this equation reduces to il = in{B(Xi)S A(i) (i = 3,...,n). (5 49) Note the similarity between this relation and Eq. (5.34). Utilizing the value for a.l determined by Eq. (5.49), integration of the appropriate equation in Eq. (5.47) gives the following expression for uiil(t) (i = 3,4,...,n):

127 U. (t) Ui (~) -.+ viB (z + X.)So(s) - A(i)}]ds lii i1o Lnil i I0 ( o3 (550o) (i = 3,4,..n). Again recall that uiil(O) = O. The expression for the other uijl(t)'s, obtained from integrating the other equations in Eq. (5.47), is (Aj -i)t Uijl(t) = Uijl(O)E (5.51) t (kj-i)(t-s) jn~iio i = 3,4,...,n J = 1,2,..b(i-l), (i+l),..,n i L j. The initial condition, Uijl(O), must be chosen so that Uij(t) is periodic with period T; therefore, v. u (0) T -(X.-i.)s U l() = ( T f e { B(z + Xi)So(s) - A(i)}ds (5.52) i = 3,4,...n J = 1,2,..,,(i-l),(i+l),...,n j L iT Equation (5.52) completes the determination of all the uijl(t)'s with the exception of the initial conditions for u121(t) and u211(t). 5.6 Determination of the Coefficients, ai2 (i = l,2,...,n), of \i in the Power Series Expansion for the Characteristic Exponents, ai(G) (i = 1,2,,...,n) Equations (5.22) and (5.23), from which aio and ail (i = 1,2,...,n)

128 are determined in the preceding two sections of this chapter, are a consequence of expanding both sides of Eq. (5.13) in power series, and equating coefficients of the like powers of p on each side of the equality sign to one another. Equation (5.22) is the result of equating the coefficients of the u terms, and Eq. (5.23) is the result of equating the coefficients of the i terms. The following equation results from 2 equating the coefficients of the, terms. 1il2 + (io i - ) 22 = -i2 il (t) - i. il(t) + v [rl(t)u (t)nl + + (t)Uin(t)] + ln O(t)ill (t) +... + n(t)Uin() (5.53) in2+ (czio - X)u in2 -a i2ino - lu(t) + v [rl(t)ui (t) +. + nl(t)uino(t )] v [FlO(t)ui (t) +.. +T n(t)Uinl(t)] nn 10 ill no in The determination of the ai2's and the uij2(t)'s follows the same pattern as that followed in the previous section for the determination of the il's and the u j(t)'s. For a12 Consider Eq. (5.53) for i = 1 and with Eq. (5.25) substituted for the uljo(t)'s (j = l,...,n) and the results of the preceding section substituted for the uljl(t)'s (j = l,...,n). Recall that clO = l, Ai == S

129 112 = (12 110(0) - U11ll(t) -j2w t + Vln[rl(t)ullo( O) +F21(t)u120(O)c + vn[l(t) (t) +... +o u (t)(t)] ini ii il no i-nl -ja~2 t u +22+ j2wsU122 - 12U120) - 11U121(t) -j20 t + V2n[11 (t) ( 0) 1 + 0 (t)uu120( O) c S (5.54) +V2n[[rlO(t)U ll(t) +... n +fnt (t)lnl] Uln2 + (jWs n)Uln2 =- a11unl(t) -j2o t +v nntli(t)uO1 () +2+r1(t)u o 20(0) S + Vnn[(t)U (t) +... +F (t)u n(t)] nn 10 il ~ no Inl The terms on the right-hand sides of the equations in Eq. (5.54) are all, except for c12 and the initial condition u121(0), known. These two undetermined constants must be chosen so that the solutions of Eq. (5.54) can be periodic in t with period T. The function u 12(t) is periodic if these constants are chosen so that the mean value of the right-hand side of the first equation in Eq. (5.54) is zero. The function u 2(t) is periodic if the right-hand side of the second equation in Eq. (5.54) contains no term of angular frequency -aw. These two requirements lead to the followin pair of simulaneous equaions ing pair of simultaneous equations.

130 Vln T {-jst u121(0) 11 Tull,(t) 12' T{ J F r20(t)e dt -().f U (0dt o Uo 11 TVn u 120(0) -j2a st [n r ( t) + r2(t) U ) ]dt 0 111 V T ull(t) 1(t) Ulnl(t) Vln Till(o +1v lTl ( ujt) + f [Fr0(t) u(0) + + r n(t) Ul ( dt T 0I~r~o~~ u (t) 20 u1' (t) uu a 2n T 121(0) ll T 121 (t) +nj2o t 12whee u 1( r= u (t) - U (). I is poile o ilif ~, T 10(0) +j2Wost vT T 2ll (t) an f Iril(t) ~'120(0' ) ~ + F (t)]dt + [ (t) u121(t) (t) +jws = ) + rF(t) u2O + + t (t)U {(5 5t -j2w t where ul2(t) = Ul2l(t) - U121(0)e It is possible to simplify Eqs. (5~55) and (5-56) by use of the relations which have been developed in the foregoing sections of this chapter. It follows from Eq. (5.20) that Vln T'Jaw t 1n f oWt) S dt vlnB(jwO)S +2s k5* o and v2n T of 17 2o(t)dt - a11 = V 2n B(-js)S - A(-Js)}- ai 1 (5.58) From Eq. (5.21) Vl T 120(0) -j. t (5s9) f [Fl(t) + r21(t) U1(~)' jlt viB(c){s 1 U120(0) s'^^T^V 1

131 and v2_ T u 1(0) +j2t (t)dt t= ) (56 ~0'120(o) (5.~ r Uz2o0O) E(0) U110(0) -2s Equations (5.59) and (5.60) can be even further reduced with the aid of Eq. (5.35), which implies that i U12o0(~) s - B( o 11JO~ ~0 v1(){s +O U110() 2s} v -2s (5.62) so Utilizing the foregoing reductions, Eqs. (5.55) and (5.56) become a12 -{ - } Ujlo( O) o +2s aln;p^ 1 1 (0 U121(.6) nB~i )1S1 + - Vin[B(jw )SO - A()]} 2(5.63) 0 +an f [{'(t) } T-+ P7 (t) uii(O +P (t) Uo]dt, 0 151 10n 110 110 and.

132 2[V2n{B(-ijs)SO - A(-js}'- ] ]o+ v-Co) ^12u1)-[ \ 121) - 2n ^s -2s v11 V2n[B(-jw)S - A( jw)]} s1s V2n +T[ i lllp(t) 6 + 11 - 2n so s o s)5 -2s +T f IFO(t) U (5.64) S{ o 120 +{rP2t) - l 120(t) + 30(t(t ) + o+ r (t) i dt This system has a unique solution if its determinant does not vanish; i.e., sV2n{B(-j-)SO - A(-JCs)} - 1 [VB(j) S+2 A am_ [ 2nt s ~ ( ) s } 11] + Uno(o)~ ], 0. (5.65) 1 20~~) u' () ^ Utilizing Eq. (5.35), this requirement becomes A = { 2alz - v B (js) S~ - A(j( ) (5 66) -v2n[B(-j )sO - A( iw)]} 0 But, introducing Eq. (5.33), Eq. (5.66) reduces to all - a21 (5.67) = U120(O): o. (5.67) Therefore, it is necessary only that all and a21 be distinct, but this has already been stated as a precondition for this analysis. Consequently, Eq. (5.67) is always satisfied for the problems considered. The solution of the two simultaneous equations, Eq. (5.63) and Eq. (5.64), is

133 1. 2 [ 12[ - V2n[B(-jws)S~ A(j )w]} ~ c| a2 - a21 -'V2n so -' s} (5.68) + {a11- Vn[B(Js) s - A(Js)]} ] u120(0) u121(o) = (2 ( C1) (5.69) where C = The right-hand side of Eq. (5.63) C2 = The right-hand side of Eq. (5.64). The functions ulj2(t) (j = 1,2,...,n) can now be determined by substituting the values for a12 and u121(0) from Eqs. (5.68) and (5.69) into Eq. (5.54) and integrating the differential equations involved. The initial conditions, ulj2(0) (j = 3,...,n), are chosen so that the ulj2(t)'s are periodic with period T. The initial condition ull2(0) is, by assumption, zero, and, as in the case of Eq. (5.37), U 22(t) is periodic for any choice of initial condition, ul22(0). Therefore, the undetermined constant, ul22(0), must be determined in the next step of the recursion procedure. For a22: Consider Eq. (5.53) for i = 2, with Eq. (5.26) substituted for the u2jo(t) (j = l,2,...,n), and the results of the foregoing section of this chapter substituted for the u2jl(t)'s (j = l,2,..,n). Recall that a20 = 2 2 = -J..

134 +jj t t + lnI rl(t)u210(0O) S + r21(t)u220(0)] + V 10(t)U211(t) +r20(t)u22(t) +... +rno(t)u2nl(t)] U222 =- a22u220(0) - 21u221(t) + V2n[ rll(t)U210(O)e +21l(t)u220(o)] (5*70) + V2n[ rO(t)u211(t) +.. +rno(t)u2nl(t)] u2n2' (Js + kn)u2n2 = 212nl(t) +j2S t + Vnn[ r(t)u210(0) + r21(t) 220(o( v r[o10(t) 2(t) +.. + rn(t)Unl(t)] 1nn l 211 no'2nl As in the case of Eq. (5.54), all the terms on the right-hand side of the equations in Eq. (5.70), except for a22 and u211(0), are known. As before, the a22 and u211(0) must be chosen so that the u2j2(t)'s (j = l,2,...,n) can be periodic in t with period T. The function u212(t) is periodic with period T if the right-hand side of the first equation in Eq. (5.70) contains no terms of angular frequency +20s, and u222(t) is periodic with period T if the right-hand side of the second equation in Eq. (5.70) has mean value zero. These two requirements lead to the following two equations.

135 2n T U210(0) T u211 (t) -jw (t a22 { + a21 - T J Flo (t)} - _ T f - e dt v+T T u220(0) -j2w t T+ 2 [ rO(t). +~ 220() rO (t)s ]dt (5.71) 0 ou210() 21 +n T r 211 (t) -w rsdt + F- - [J O (t) 211) +In C t -dt o ou210 (0) no (0 ]st T 0 10 U 210 o) no u210 2nT +j2t u2(0) a2 Tu22(t) a22 -{ T O rlO(t) dt} u ( T I O dt o 220 0 220 V2n T U210(0) +jst (72) 1 T u2(t) (t) +2n f F 11,, ()"2ni + f J[ so(t) + *** +Rno(t) ]idt 0 1 220 u220 ja2w t here u211(t) 2(t) - 211(n )e Then, utilizing relations analogous to those given in Eqs. (5.57), (5.58), (5.59), (5.60), (5.61) and (5.62), Eqs. (5.71) and (5.72) can be reduced to a21 ln [B(jw )SO - A(j )] a2 { ~ U~210(2) } u211(o ) +2s vlB(j )S +{ 21 - Vln[B(jws)SO - A(w)]} S (573) + T — [(FlO(t) %_A)+(o0T - no() l Ta u211(t) r (t) _2_ ]t ( -(t) U 20(0) +. (t) (Ut dt o....... no u210kO) ] ~ vnv B(-jw )S 2 4 - {2,u 0(o,... o A211o} 1 220 + V nB(-jc )S s +{a., -V[B(-jws)SO- i(-~.~} 1 s. V2n T8 1 -2n^ s - s) ro,qo^^ ^1^,^ ^ -21 21.' ^~~~~~~~~~~ ~oo ^^^^!t

136 The determinant of this pair of simultaneous equations is A nB(-js )sr - )2S Vln[B(j s)SO - (j).2. 20 2100 (575) and, with the aid of Eqs. (5.42) and (5.33), this becomes u (0) [ll - C21] (5. 76) This determinant is nonzero for the same reason that Eq. (5.67) is nonzero. Therefore, the above pair of simultaneous equations has a solution which is as follows: |l - 21l[{21 - V2n[B(ijWs)s _ A(-j ))]} C3 (5 77) +"21 - vln[B(jws)S - A(j- )]} C4] and u210(O) u21(0) = 1 c [C4 - C,3, (5.78) where C = The right-hand side of Eq. (5.73) C4 = The right-hand side of Eq. (5.74). The functions u2j2(t) (j = l,...,n) can now be determined by substituting the values for a22 and u211(0) from Eqs. (5.77) and (5.78) into Eq. (5.70) and integrating the differential equations involved. The initial conditions, u2j2(0) (j A 2; j = 1,3,4,...,n), are chosen so that the u2j2(t)'s are periodic with period T. As in the case of Eq. (5.43), u212(t) is periodic for any choice of initial condition, u212(0). This undetermined constant must be determined in the next step of the

137 recursion procedure. The initial condition, u222(0), is chosen equal to zero (i.e., all uiik(0) = 0, k = 1,2,...). For ai2 (i = 3,4..,n): Consider Eq. (5.53) for i = 3,4,...,n with Eq. (5.27) substituted for the uiJo(t)'s (j = 1,2,...,n) and the results of the foregoing section of this chapter substituted for the uij(t)'s (i,j = l,..,n). Recall that ao = Xi. Uil2 + (hi - 1)Uil2 = - cilUill(t) + Vln ril(t)Uiio(O)] + Vln[ lO(t)Uill(t) +. + no(t)inl(t) ii2 - i2Uiio() - ailUiil(t) + Vin[ ril(t)uiio(O)] (5.79) +vi[ lO(t)Uill(t) +... + rno(t)uin(t)] in O ill) no inl in2 + (i - n)in2 = - CilUinl(t) + Vn [ril(t)Uiio()] v Inr (t(t) +*' ++ (t) nl(t)] nn )UO^ ill no inl (i = 3,4,...,n) Once again, all the terms on the right-hand side of the above system of equations, except for ail, are known. In this case there is only one undetermined constant because ( - Xk) (j A k; J,k = 3,4,...,n)is not an integer multiple of joT. The constant ai2 must be chosen so that the

138 uij2(t)'s can be periodic with period T. Referring to the ith equation in Eq. (5o79) it can be seen that uii2(t) can be periodic with period T only if the mean value of the right-hand side of this equation is zero. Therefore, ai2 is chosen as follows:'cil T uiil(t) Vin T i2 u ii( r dt + T i ri(t)dt 0 UO 0 (5.80) Vin T Uill (t) Uinl(t) Vin/ If[F (t) Ui + + no(t) u(0 ]dt (i = 3,4,..,n), or (t) i2 = Vin B(i)S iZ I T ~vi~Tr~ uill(t) ~U~ ~ ~(t)(5.81) +in - [rl (t) Uill +. +I'no (t) (i = 3,4,...,n). T 0_Uiin) Y now be (i = The functions ui2(t) (j = l,2,...,n; i = 3,...,n) can now be determined by substituting the value of ai2 from Eq. (5.81) into Eq. (5.79) and integrating the differential equations involved. The initial conditions, uij2(), are chosen so that the u. j(t)'s are periodic with period T. Since uii (0) is arbitrarily chosen equal to zero, none of the initial conditions remains undetermined; therefore, this situation is somewhat different from that for i = 1 or 2.

139 5.7 Determination of the Coefficients, aik i = 1,2,...,n; k = 3,4,...), of the Higher-Degree Terms in the Power Series Expansion for the Characteristic Exponents, a.(C) (i = 1,2,...n) The determination of the aik's for k > 2 follows exactly the same pattern used to determine the ail's and the ai2s. The system of differential equations which involves the uijk(t)'s is obtained by equating the coefficients of the k terms from each side of Eq. (5.13) to one another. This system of differential equations is similar to Eqs. (5.54), (5.70), or (5.79), depending on which integer i is considered. For i = 1 or 2, the right-hand side of this system is known except for cik and one initial condition from the previous recursion step, and, just as in the previous cases, these two undetermined constants are chosen so that the uijk(t)'s can be periodic with period T. The system differential equations are then integrated, and the initial conditions uijk(0) chosen so that the uik(t)'s are periodic with period T. One initial condition remains undetermined and must be determined in the next step of the recursion procedure, and u i(O) is chosen equal to zero. For (i = 314,...,n), the procedure is, as was just demonstrated for ai2 (i = 3,4,...,n), simpler. The right-hand side of the appropriate system of differential equations is known except for Cik. This constant is chosen so that the uijk(t)'s can be periodic with period T. The system is then integrated, and the initial conditions, uiJk(O), chosen so that the uijk(t)'s are periodic with period T. Again, Uiik(O) is chosen equal to zero. In this case, there is no undetermined initial condition which must be determined in the next step of the recursion procedure.

CHAPTER VI ALGORITHM FOR THE APPLICATION OF THE TECHNIQUES DEVELOPED IN THIS STUDY 6.1 Introduction In this chapter the key results from the foregoing chapters are restated so that they form an algorithm for the application of the techniques developed in this study. An attempt has been made to make this chapter self-contained except for a few references to other chapters. It is hoped, thereby, to obviate a detailed examination of the first five chapters in this study by the reader who is interested only in applying the techniques which are developed here to a given nonlinear feedback system. This algorithm is divided into two main parts. The first part presents the procedures to be followed for the determination of the first approximation to the periodic response and the first approximations to the stability criteria. As is demonstrated by the example worked out at the end of this chapter, these first approximations can be quite satisfactory. The second part of the algorithm presents the procedures to be followed for the determination of higher-order approximations to both the periodic response and the stability criteria. Owing to the existence of two canonical forms, each of the two main parts of the algorithm is further divided into two parallel groups of sections. One group is applicable to the first canonical form, and the other group is applicable to the second canonical form. In addition, one common section, which ac acts as an initial step, is presented in the first part of the algorithm. This section discusses the selection of a canonical form. 140

Consequently, this chapter is arranged in the following pattern: a) In Section 6.2, the selection of the pertinent canonical form is discussed. b) In Section 6.3, the method is presented for determining the first approximation to the periodic responses of a nonlinear feedback system which is characterized by the first canonical form. c) In Section 6.4, the method is presented for determining the first approximations to the stability criteria for a nonlinear feedback system which is characterized by the first canonical form. d) In Section 6.5, the method is presented for determining the first approximation to the periodic response of a nonlinear feedback system which is characterized by the second canonical form. e) In Section 6.6, the method is presented for determining the first approximations to stability criteria for a nonlinear feedback system which is characterized by the second canonical form. The above five sections constitute the first part of the algorithm which is presented in this chapter. The second part of the algorithm is composed of the following sections. f) In Section 6.7, the method is presented for determining higher-order approximations to the periodic responses of a nonlinear feedback system which is characterized by the first canonical form. g) In Section 6.8, the method is presented for determining

142 higher-order approximations to the stability criteria for a nonlinear feedback system which is characterized by the first canonical form. h) In Section 6.9, the method is presented for determining higher-order approximations to the periodic responses of a nonlinear feedback system which is characterized by the second canonical form. i) In Section 6.10, the method is presented for determining higher-order approximations to the stability criteria for a nonlinear feedback system which is characterized by the second canonical form. The final section of this chapter, Section 6.11, presents an example of the application of this algorithm to a third-order system which is characterized by the first canonical form. 6.2 Selection of a Canonical Form If it is assumed at the outset that the periodic response in question is dominated by a component in its frequency spectrum of frequency ~, then the selection of the pertinent canonical form is straightforward. The first canonical form applies if the frequency spectrum of the forcing function, e(t), contains a term of frequency, and the second canonical form applies if the frequency spectrum of the forcing function, e(t), does not contain a term of frequency. Unfortunately, the above assumption regarding a knowledge of the dominant term in the response hides the basic problem in the selection of the pertinent canonical form; namely, the a priori determination of the angular frequency, co5, of the dominant term. As has been discussed

143 previously (see the last paragraph of Section 2.3), there is no general method for making an a priori selection of T (T = -). The reason that..... s s 0 s there is no such method is intimately connected with the lack of a complete solution of the "inverse perturbation problem" as discussed in Chapter II. In any event, it is sufficient for the purposes of this chapter to say that it is necessary in every case to guess o. Conses quently, it is worthwhile to consider briefly how one might make such a guess. The development of both the first and second canonical forms in Chapter II begins with a consideration of the following linear system: A(z)x - B(z)Nx =. (2.5) The characteristic roots of this system are functions of N. It is assumed that for each of certain values of N this system has a pair of complexconjugate, purely-imaginary roots which are denoted + Jo. These special values for N are denoted Nj, where j is an integer index. It is shown in Chapter II that by the use of suitable transformations the characteristic roots + Jo. can be shifted so that their magnitude becomes a rational multiple of the fundamental angular frequency of the forcing function, e(t). The question of "guessing an o" is really the question of deciding which pair -jcN. is to be shifted and to where it is to be shifted [i.e., which rational multiple of the fundamental angular frequency of e(t)j. In many cases the decisions are obvious. For example, consider a system whose linear part, A, has poles and zeros located in the complex plane as shown in Fig. 6.1.

144 Im(z) POLE ZERO Re(z) POLE B(z) FIG. 6.1 LOCATION OF POLES AND ZEROS OF A POSSIBLE A(z The root locus pattern for such a system is sketched in Fig. 6.2. i lm(z) LOCUS OF POLE LOCATION AS N VARIES Re(z) FIG. 6.2 ROOT LOCUS ASSOCIATED WITH THE AB PORTRAYED IN FIG. 6.1 It is apparent that for the above system there exists only one N.; consequently, there is only one pair + jN.. Let the forcing function under

145 consideration be a simple sinusoid of angular frequency we, where a) is 6 ei in the neighborhood of 20l, as illustrated in Fig. 6.3. 1 Im (z) +jWe +jWN1 Re(z) -jW FIG. 6.3 IOCATION OF ~j and ~jow IN THE COMPLEX:HANE 1 a) Clearly, the poles + JoN should be shifted to ~ j e 1 2 Fortunately, this type of rather-crude reasoning suffices for many nonlinear feedback systems. This is especially true of those systems for which the characteristic roots of Eq. (2.5) near the imaginary axis are clustered in a narrow frequency interval. In more complicated systems the problem becomes more difficult, and it may be necessary to try several possible )'s before a suitable one is discovered. Just as the "inverse perturbation problem" has more than one solution, there may be more than one suitable a. However, once an ta is selected the selection of the s s pertinent canonical form follows immediately, as is discussed at the beginning of this section. Note that the same canonical form is not

146 necessarily the pertinent one for each of the possible w's. For example, consider a simple system under the influence of a purely-sinusoidal forcing function of angular frequency ). Let this system have two possible oN's; e the first in the neighborhood of w, and the second in the neighborhood e of 2. If the first aN. is chosen, the first canonical form is applicable; and if the second N is chosen, the second canonical form is applicable. It is quite possible that either will give the correct solution. However, one may lead to a better first approximation or a more rapidly converging series representation than the other. 6.3 First Approximation to the Periodic Response: First Canonical Form The step presented in this section is quite simple and comprises nothing more than the solving of the following equation: A(jw) ) +s n -. + E =, (3.16) 0 B(isw) 10 0 where: (1) The solution of this equation is of the form c10 = cO(E+s, 1 ao), a complex function of two real variables. Letting a jQ 10 (a, @ real), the solution of Eq. (3.16) can be written a = a(E+s, X ) e = e(E+S, ) 1 E can always be assumed to be real without a loss in generality.

14.7 Note that the solutions "a" and 0 may be considered to be functions of system parameters other than E+ and ws; e.g., some component value in the linear network characterized A(jwo ) by (2) The first approximation to the periodic response of the nonlinear feedback system is jW t -jo t 10 +c10 = a cos [c t + e] s (3) n, the zeroth-order describing function, is a function of "a" and is defined as follows: 1 T/2 ju t * -Jut -jw t O C1O ~r~,C1o* jt coto cc T | { 10 E }d dt, -T/2 where (see Appendix C) the function F(x) characterizes the nonlinear element in the feedback system and T, in this case, equals an integer multiple of 2 A(jo ) ) (4) B~y is the reciprocal of the transfer function for the linear part of the nonlinear feedback system. (5) )s is the angular frequency of the first approximation. The determination of this quantity is discussed in the section preceding this one. +jw t (6) E+S is the Fourier coefficient of the e term in the Fourier series for the forcing function, e(t). The solution of Eq. (3.16) can be accomplished by many standard techniques. Probably the most useful one is a graphical technique. In this technique

148 A(jw ) the locus of n as "a" varies and the locus of - as w varies are o B(jwS) s B(Jis plotted on the same complex plane. For a given amplitude, "a", and a A(j) ) given angular frequency, w, the vector whose head is at B(J and whose s +S S) tail is at n (a) gives the magnitude and argument ofE. An example o c 10 is shown in Fig. 6.4. IMAGINARY AXIS A(jwo) - LOCUS OF B() \B(Jo) LOCUS OF no(i) 0I0 REAL AXIS a= a INCREASING w INCREASING FIG. 6.4 EXAMPLE OF A POSSIBLE GRAPHICAL METHOD FOR SOLVING EQ. (3.16). Needless to say, there are many other ways for solving Eq. (3.16). A convenient and significant way of presenting the solution of Eq. (3.16) is to plot "a" vs. ws for constant E, and e vs. as for constant E. Among s s other things, such a presentation allows those points at which the Jacobian condition is not satisfied to be easily determined. These will be the points at which a~ becomes infinite. SE

149 6.4 First Approximations to the Stability Criteria: First Canonical Form As is discussed in detail in Chapter V, there are n so-called characteristic exponents associated with a periodic solution of an nthorder nonlinear system. The signs of the real parts of these characteristic exponents determine whether or not the periodic solution in question is asymptotically stable: if all the real parts are negative, it is; if at least one of the real parts is positive, it is not stable. In this study, the characteristic exponents are denoted, a1(p), xg2(),., a n(A), and are functions of the perturbing parameter A which is introduced in Chapter II. The true stability criteria —that is, the ones corresponding to the physical system —are obtained by considering al(>), a2(i),..., n (AL) evaluated at A = 1. A procedure which develops the ai(u)'s recursively as power series in A is presented in Chapter V. With this procedure, the characteristic exponents, ai(n) (i = 1, 2,..., n), for a given problem can be determined in the following form: 2 o^i(~) = i0o + pCil + L Xa.i2 + and the "true" characteristic exponent will be ai(1) = CiO + ail + 2 + (i =, 2,..., n) The approximate stability criteria presented in this section are based upon the assumption that the first two terms in the power series for a.i(p) form a satisfactory approximation to ci (p.); that is, that ai(l) o + a (i = 1, 2,... n). In particular, it is assumed that the real parts of ai(l) and (ai0 + ail) (i = 1, 2,..., n) have the same sign, although they may differ in magnitude.

150 Moreover, it may be assumed that a (1) and a (1) are the significant or dominant characteristic exponents and that the remaining (n-2) characteristic roots can be ignored as far as a first approximation to the stability criteria is concerned. This last assumption must, however, be employed judiciously. For example, if the linear network has several pairs of poles in its transfer function which are near + jws, it is a questionable practice to consider a (1) and a 2(1) only. On the other hand, if the linear network has one pair of poles in its transfer function near the + ji and its remaining (n-2) poles far in the left-hand plane, the assumption that a (1) and a2(1) dominate is a reasonable one. In what follows, approximate stability criteria based upon the assumption that ao(l) and a2(1) dominate are presented first. These criteria are then modified so that they take all the n characteristic exponents into consideration. Now if it is assumed that a (l) and a2(1) are the dominant charac teristic exponents and that (alO + all) and (a20 + a2 ) are satisfactory approximations to a0(1) and a2(1), respectively, it is necessary to consider al and a21 only because 10 = +jaw and a0 = -jw [see Eq. (5.24)]. II s i20 IU~ S 20 S It is known from Eq. (5.31) and indirectly from Eq. (5.41) that all and a21 are equal to the two roots of the following quadratic equation: a2 Vln [B(s )So - A(jw)] + 2n [B(_j)S - A(J)]} + Vln2n [B(Js )so - A(jw)] [- )S - A(-s)] (531) ^ln~n s s S,] - +2S -2s - V B(jws )B(-jw s)SO S = 0 In 2n S 0

151 T/2 jWt -jwt o - Tpx dt=10-] -T/2 x=p(t) and T/2 S2s= 1 / dF I J2stF "o TJSx d~t. -T/2 x=p(t) SO is a function of "a" only, and S+ and S2 are functions of both 9 0 o o and "a". Since the constant term in Eq. (5.31) is equal to the product of cll and C21, it is apparent that one stability criterion is [ A(jw )) 0 A(-jw s +25 -2s - B jo) s o B(-jo J) o o In addition to being an approximate stability criterion, the above expression is also the Jacobian condition [see Eq. (3.22)]. This connection between the Jacobian condition and the approximate stability criterion is pointed out in Chapter V. It is particularly important to recall that along the boundary, ____- - s +2S - 2s o A B( js) BA(-jS) +2s -2s )J Lo wj 0 0 (S +s) the slope(i~ + ) becomes infinite. A second approximate stability criterion is, considering Eq. (5-31) again, clearly Re vln B(jws) [S -( I < (5.34 In order to employ this criterion, vln must be evaluated. However, since vln is a function of all the characteristic roots,; (j = 1, 2,..., n),

152 of the "p. = 0" differential equation (see Section 2.4 and Appendix E), it is necessary to calculate the.j's before this second stability criterion can be utilized. As is shown in Chapter II, these.X's (j = 1, J 2,..., n) are related to the characteristic roots, j ( = 1, 2,..., n), of the linear differential equation, A(z)x - N. B(z)x = 0, (2.6) by the relation x. = A. (j = 1, 2,..., n), 0 ) s where P =. Note that the parameter P is a convenient measure for 10Nj angular frequency deviation from some fixed angular frequency, c. It, I, is employed extensively in the working of an example in the last seo tion of this chapter. In particular, the v.'s (j = 1,..., n) can be considered to be the following functions of P (see Appendix E): Vn (P = 1) vJ,() j= (E3) where v (=l) = (-) (E2) (1 I)" (j-) I) ((j+l) i n j and j = 1, 2,..., n. The term v ln() can be determined from the above two expressions and, then, substituted into the stability criterion, Eq. (5.34). Note that for this stability criterion it is necessary to know only the argument of the complex term, vln(P), and it can be seen from Eq. (E3) that the argument of v n(P) (j = 1,..., n) is independent of 8 (a real positive number). If it is not assumed that the two characteristic exponents,

153 al(1) and a2(1), dominate, then the first approximations, aio + ail (i = 3, 4,... n), to the remaining (n-2) characteristic exponents are determined by the following expressions: aio (i = (i 3, 4,..., n), (5.24) and aY: Vin B(%i)( S n i(ki) i1 -i\n i { -o B- } (i = 3, 4,... n) (5 B(o i) ail = vin() B( i){ (i = 3, 4..., n). 6.5 First Approximation to the Periodic Response: Second Canonical Form The determination of the first approximation, jw t. -jwo t p(t) = c1O + c1O + h(t; o ), 10 10 s to a periodic response is slightly more complicated in the case of the second canonical form. The complication is caused by the presence of the term h(t; w ) in the first approximation. It is shown in Section 2.5 that h(t; w ) is equivalent to the s steady-state response of a linear system which is characterized by the transfer function, z(z) = B(z) z + [ C() - N B( to the forcing function, e(t). Therefore, if e(t) is expanded in a Fourier series,

154 e(t) = Ee e, (2.32) = -oo where = qe (q an integer), the expression for h(t; w ) is I-o (.2)n + n [ C( ) - NiB ] where s = pa (p an integer) and, again, = -. Note that h(t; cs) is a function of a) (equivalently, A). S After h(t; a) ) is determined, c10 is determined from the requirement that it satisfy the following describing-function-like relation. A(jwo ) n - (jcS = 0, (3.29) o B(^TT = 0, wnere T/2 T/2 1 jWj t. -jw t -jw t o c T F C s + c* e + h(t; )} dt. (C8) 10 -T/2 It is also true that c 1 satisfies the complex conjugate of Eq. (3.29). The significant difference between the zeroth-order describing function defined above and the one defined for the first canonical form is that the above zeroth-order describing function is not necessarily real and is a function of 9 and os (or 3) as well as "a." (Recall that c1 = a JQ 2~.) Once the appropriate expressions are substituted into Eq. (3.29) A(jw ) for n and B(j-' its solution can be accomplished by a variety of methods. 0 FC2WT The best method will depend upon the special characteristics of the particular problem under consideration. Graphical techniques are, of course, often useful.

155 6.6 First Approximations to the Stability Criteria: Second Canonical Form The determination of the first approximations to the stability criteria for the second canonical form resembles closely the same determination for the first canonical form (see Section 6.4). If it is assumed that the two characteristic exponents, a(l) and a (1), dominate, then the first approximate stability criterion is, again, A(jw A(-j s- S S > 0o [ B j 0Lo o B-j o o where T/2 o F 1 dt "o T dx -T/2 x-p(t) and T/2 s2- 1 dF 2st _s 1 / d~t o T J/ 0x c -T/2 x=p(t) The difference between the above case and the first canonical form is that here p(t) is of the form, j3 t -jo t p(t) = c10E + ce + h(t; u) whereas, in the case of the first canonical form, it is of the form, jc t -jo t p(t) = c10C + cc S The second approximate stability criterion is, similarly to the first canonical form, A(j ) Re vlB(jw <(j ) 0. (5.34) Re{Vln B(Js) [SO B ] The term vl is evaluated just as it was for the first canonical form.

156 If it is not assumed that the two characteristic exponents, a(l) and a2(1), dominate, then the first approximations to the remaining (n-2) characteristic exponents are determined, again similarly to the first canonical form, by the following expressions: ci= i= ii (i= 3, 4..., n) (5.24) and ail.in() B(Bi) { SO - () } (i = 3 4,..., n). (5.49) il in i~ Io B(_iJ 6.7 Higher-Order Approximations to the Periodic Response: First Canonical Form In Chapter II, a perturbation parameter is artificially introduced into the differential equation which characterizes the nonlinear feedback system of interest. The periodic solutions of the resulting i-dependent differential equation are denoted p(t, i). At u = 1, these v-dependent periodic solutions, p(t,;L), correspond to periodic responses of the original nonlinear feedback system; i.e., the desired responses are P(t, 1). A coordinate transformation is also introduced in Chapter II. This transformation converts the single nth-order differential equation, Eq. (2.17), which characterizes the problem into a system of n firstorder differential equations, Eq. (2,27). The desired solutions of this system of equationsaredenoted [cj(t, p)], where [cj(t, ))] is an ncomponent column matrix. The desirable property which [cj(t, 0p)] possesses is that it satisfies the following "periodicity condition": X T [c (T, )e T - c,(O, pi)] - ( = 1,..., n), (3.2)

157 where the j.'s are the characteristic roots of the "I = 0" differential equation. This condition is referred to as a "periodicity condition" because it guarantees that the solution, p(t, p), of the nth-order differential equation, which is related to the cj(t,.)'s (j = 1, 2,..., n) by A i t2t t A t p( ) (t, ) = c l(t, ) + c(t, ) +.. + c(t, )e n is periodic in t with period T. Therefore, determination of the c (t, i)'s is equivalent to a determination of the periodic solution, p(t, A). A recursion procedure is developed in Chapter IV for the determination of the desired cj(t, p)'s as power series in A; i.e., cJ(t, ) = cjO + (t) +.c.+... (j= 1, 2,..., n). The first terms in these power series, cj0 (j = 1, 2,..., n), determine the first approximation, +jw t -J t p(t) = c10 + c to the periodic response. It is shown in Chapter IV that cj = 0 for * j = 3, 4,..., n and c20 = c1. The term c1O is determined by the methods of Section 6.3. The higher-order terms in the power series expansions for the c(t, )'s are determined by a rather straightforward recursion procedure. A typical step in this recursion procedure follows the ensuing pattern: a) At the outset of this step all the cj(t)'s, except the jk mean values of clk(t) and c2k(t) (equivalently: the initial conditions alk and a2k), are completely known. Furthermore, all the terms in the power series for the cj(t, pI)'s which

158 precede the c jk(t)-terms are completely known. b) It is first necessary to determine the mean values for lk (t) and c2k(t). [Recall that, because k = jws and 2 = -jWo, clk(t) and c2k(t), in contradistinction to ck(t) (j = 3, 4,..., n), are periodic with period T. jk [See Eq. (4.3).] Owing to the requirement that the solution p(t,.) be real, the mean value of ck(t) must equal the complex conjugate of the mean value of c k(t). This is a first relation which can be employed in the determination of these unknown mean values. The second and key relation is the kth-order describing function relation, A(jw ) nk - -0 (4.31) k B(jw) where nk is the kth-order describing function. This describing function is derived in the following manner. First, co 00. t t 1v C c k n k x = e E clk~ +... +C nk k=O k is formally substituted into F(x), and the resulting function is expanded in a power series in C; i.e., x\t Xt 00o FcC1 +.. + c n } = k Fk 4. (4.20) Each 1, 2,.) is a f n of k=O Each Fk (k = 0, l, 2,...) is a function of clO,..., Clk, c20, *.'. C2k'., CnO'..0 Cnk' and t. Now, assuming that all the steps in the recursion procedure occurring before the present step has been carried out properly, all the ck's, except for the mean value of clk(t), are known functions of time. Consequently, Fk can be considered, in

159 the present step in the recursion procedure, to be a function of t and mean value of clk(t) only. Moreover, Fk is periodic in t with period T for any value of the mean value of clk(t) [see the discussion following Eq. (4.20]. Expanding F(t c) [ck = the mean value of ck(t)] in a Fourier di F k(t, Clk) lk k series gives Co o o F(t, clk)= Fk (Clk)E (4.21) a=-oo 2m where = T-. The kth-order describing function, nk, is defined as follows: F+s k n, = clk jik +S jO t where F is the coefficient of the e -term in the Fourier k series, Eq. (4.21). It should be noted at this point that the Jacobian expression arises in the solution of the kth-order describing-function relation, Eq. (4.31). If Eq. (4.31) and its complex conjugate are considered to be a pair of simultaneous equations in clk and (clk)*, the resulting system of equations, 0~ A(jS) 0 nk - B(J ) Clk - s a^, o * A(-js) o ) lk - ^ ~J (c7^ = 0, reduces to a linear inhomogeneous system of equations whose system determinant is just the Jacobian expression. This statement is not proved here, but it is shown to be the case

160 for the example which is worked out in Section 6.11. A rigorous proof of this statement is given by Coddington and Levinson ([10], pages 363-364) in their proof of the uniqueness of the solutions obtained by a recursion procedure which is closely related to the recursion procedure of Chapter IV. In any event, this is another demonstration of the importance of the Jacobian expression. c) After the mean value of clk(t) is determined, the next terms, in the power series Cj(k+l)(t) (j = 1, 2,..., n), can be completely determined, except for the mean value of cl(k+l) (t), through the following expression: t - \.s Cj(k+l)(t) = aj(k+l) + Vjn g(s)ds (j =, 2,., n) (6.1) 0 (k = 1, 2,... ) [See Eq. (4.10)], where aj(k+l)is the initial condition for Cj(k+l) (t) and c(k+l) n= B(z) t gk(t) = B(z)Fk(t) -. cjk(t)A(z) n, (z- t) j=l where (k = 1, 2,... ). This expression for gk is obtained by substituting the c (t)'s (j = 1, 2,..., n) which have jk been determined into Eq. (4.23). A significant characteristic of this substitution is that it results in a periodic function, gk(t), whose Fourier series contains no terms of angular frequency X. Therefore, it is known that the cojw ts -jw t efficients of e and e in Eq. (6.2) always add up to zero if a suitable set of c.k(t)'s has been substituted into

161 Eq. (4.23). Consequently, at this point in the recursion procedure, if the mean value of clk(t) has been determined ~Jst by Eq. (4.31), the e terms in gk can be neglected. It is shown in Section 6.11, in the example, that this characteristic of gk(t) simplifies the operation characterized by Eq. (6.1). After the integration in Eq. (6.1) has been carried out, the initial conditions, aj(k+l) must be determined. This is easily done for j = 3, 4,..., n. It is only necessary to select the aj(k+l s (j = 3, 4,., n) so that the corresponding Cj(k+l)(t) (j = 3, 4,..., n) satisfies the "periodicity condition", A.T Cj(k+l)(T)c - C(k+l)(O) = 0. (43) Since the definite integral in Eq. (6.1) is, after integration, made up of a constant plus an exponentially-damped periodic function, it is merely necessary to choose aj(k+l) so that it cancels the constant term which results from the integration. Moreover, the procedure need not really be carried out: the constant term which results from the integration can just be neglected. In the case of j = 1 and j = 2, the determination of aj(k+l) is slightly more complicated than for j = 3, 4,..., n. Because j = Jos and 2 = -j(, the "periodicity condition", X T cj(k+l) (T)E c(k+l)() = 0 ( =, 2), is satisfied for all al(k+l) and a2(k+l) Therefore, the

162 "periodicity condition" does not determine al(k) and a2(k+l)' In other words, both Cl(k+l)(t) and 2(k+l)(t) are periodic with period T for all values of their initial conditions. However, since these initial conditions determine the mean values of cl(k+l)(t) and c2(k+l)(t), the fact that al(k+l) and a2(k+l) are undetermined is equivalent to Cl(k+l) and C2(k+) being undetermined. However, just as the present step in the recursion proO O cedure was begun with clk and ck undetermined, the next 0 step in the recursion procedure begins with cl(k+l) and c2(k+l) undetermined. Therefore, the foregoing constitutes one step in the recursion procedure. 6.8 Higher-Order Approximations to the Stability Criteria: First Canonical Form As is discussed in Section 6.4, whether or not a periodic response is asymptotically stable is determined by its associated characteristic exponents, aci(l) (i = 1, 2,..., n). A recursion procedure is developed in Chapter V for the determination of these characteristic exponents as power series in.; i.e., 2 ai(,) = ao + i +... (i = 1, 2,..., n). The approximate stability criteria which result from considering a i and ail only are given in Section 6.4. The procedure to be followed for the determination of aik (i = 1, 2,..., n; k = 2, 3,....) is outlined in this section. Since in many practical situations a knowledge of aiO, ail, and a.i2 only will suffice and since each of the steps in the recursion procedure developed in Chapter V follows the same pattern, only the step which determines a2i (i = 1, 2,..., n) is outlined here.

163 The determination of i2(i = 1, 2,..., n) follows, then,the ensuing pattern: a) At the outset of the determination of ai2 (i = 1, 2,..., n), both ai0 and cil (i = 1, 2,..., n) are assumed to have been determined by the methods of Section 6.4. b) The functions r1F(t) (j = 1, 2,..., n) are evaluated from the expression, rjo(t) = B(z+Aj)So(t) - A(Xj) (z - ) (5.20) where j = 1, 2,..., n and S (t) d dF Ed jcn t -jW t x = CeOe + C10E 12 0 210 u (o) u (o) c) The quantities -IY and 20) are evaluated from 110) u220( u120(0) ll - vln [B(s - A(s)] } (5 110 (0) +2s vln B(jws) So and u (0) { - V [B(-jW5)S- A(-w)] 210 21 2n S o, 220 V2n B(-jw ) so..42) 0 -2s n -2s d) The quantities S, S1, and S1 are determined from the following expressions: o 1 T/2 1= - T s(t) dt, -T/2

164 T/2 j2 t S~ = t)j2 t 1 I f 1t)e S, -T/2 where Sl(t) is the second term in the power series expansion for; i.e., d x = IF So^ ^2S + ^ +....L dx ] 1 So +ls~l + >2S2 + *4 * (5.19) x=p(t, y) The next two steps pertain only to the determination of a12' e) The appropriate quantities evaluated in b) and c) are substituted into Eqs. (5.36), (5.37), and (5.38). These equations determine Ulll(t), U121(t),.*, Ulnl(t). The term involving acl in the integrand of Eq. (5.36) has been chosen so that the mean value of this integrand will be zero; therefore, the constant terms in this integrand can be ignored. The term involving all in the integrand of Eq. (5.37) has been -j2w t chosen so that this integrand will have no e -term in -j2w t its Fourier series; therefore, the e s -terms can be ignored in the integrand of Eq. (5.38). The initial conditions, U111( ) 21(O), lnl(O) must be chosen so that all the functions, Uljl(t) (j = 1, 2,..., n), are periodic with period T. Therefore, for j = 3, 4,..., n, the initial conditions, Uljl(O), are chosen so (Xj-Jw )t Uljl(O)E (J = 3, 4, tem, n),

165 (A.- j )t cancels the -term which results from the integration in Eq. (5.38). In the case of u12(0), the function ul21(t) is periodic with period T for all u121(0); therefore, a simple periodicity condition is not sufficient for the determination of ul2l(0), and it must be determined later in the recursion procedure. The term ull (0) is merely a scale factor which is never evaluated, and ull1(0) has been, as it may be, set equal to zero arbitrarily. f) 12 is then determined by substituting the quantities and functions which have been determined in the foregoing steps into Eq. (5.68). The next two steps pertain only to the determination of a22. g) The functions u21l(t), u221(t),..., unl(t) are determined through Eqs. (5.43), (5.44), and (5.45). The procedure to be followed here is analogous to that followed in the determination of Ulll(t), Ul2l(t),. and ulnl(t) [see e). h) a22 is determined by substituting the appropriate quantities and functions into Eq. (5.77). All the necessary quantities and functions have been determined in the foregoing steps. The next two steps pertain only to the determination of ai2 (i = 3, 4,...,n). i) The functions uil(t) (i = 3, 4,..., n; j = 1, 2,..., n) are determined through Eqs. (5.50) and (5.51). The term ail in the integrand of the integral in Eq. (5.50) has been chosen so that this integrand has mean value zero; therefore, the constant terms in this integrand can be ignored.

166 j) 4i2 (i = 3, 4,..., n) is determined by substituting the appropriate quantities and functions into Eq. (5.81). All the necessary quantities and functions have been determined in the foregoing steps. The determination of the aik's for k > 2 follows a pattern which is analogous to the foregoing pattern. This determination is discussed in Section 5.7. 6.9 Higher-Order Approximations to the Periodic Response: Second Canonical Form The determination of the higher-order approximations to the response of those nonlinear feedback systems which are characterized by the second canonical form is analogous to the same determination for nonlinear feedback systems which are characterized by the first canonical form. The only difference between the method presented in Section 6.7 for the first canonical form and the method needed here for the second canonical form is that here the functions, Fk, are derived from the series expansion, K.it %nt ~~00 F{c 1e.* + cne + h(t; =s) Z k Ffc^t +...+c ch(t, w5)} = F^, (4.33) k=O instead of Eq. (4.20). 6.10 Higher-Order Approximations to the Stability Criteria: Second Canonical Form The determination of the i2's (i = 1, 2,..., n) for those nonlinear feedback systems which are characterized by the second canonical form is analogous to this determination for nonlinear feedback systems which are characterized by the first canonical form. The only difference between the method presented in Section 6.8 for the first canonical form and the

167 method needed here for the second canonical form is that here the functions, Sk, are derived from the series expansion, dF 2 I = So + aS1 + L 2S2 +. Alt A t x c= C +... +ce n + h(t;w) ~i-11 n s instead of dF 2 d So + LS1 + + 2. 0x + 1 S2+2 Xt X t x= ce 1 +.. +c n l n 6.11 An Example of the Application of the Algorithm Presented in this Chapter In this section, the algorithm which is presented in this chapter is used to determine approximations to the periodic responses of a thirdorder nonlinear feedback system and, further, approximations to the characteristic exponents associated with these periodic responses. The block diagram of the nonlinear feedback system which is considered here is shown in Fig. 6.5. e(t) 2Es cos c t + e + (0.8 x lO z ~W(^ 3 2 2 6 + |^ z + lOOz + 10,800z + 10 ~~ (1y00xl6x 66x106 y = (-1500 x 10-6)x + (6.66 x 10-6)x3 FIG. 6.5 BLOCK DIAGRAM FOR THE THIRD-ORDER NONLINEAR FEEDBACK SYSTEM CONSIDERED IN THIS EXAMPLE

168 Consequently, A(z) = z3 + lOOz2 + 10,800z + 106 B(z) = - (0.8 x 106)z, (6.3) F(x) = (-1500 x 106 )x + (6.66 x 10 )x3 jo t -jt t e(t) =E + e (E, real), and the differential equation which characterizes this nonlinear feedback system is d3x d2x dx + 100 + 10800 + 10 x dt dt + (0.8 x 10 ) t [(-1500 x 10 )x + (6.66 x 10 )x3 (6.4) + 2E cos wet] = O. e It is assumed throughout this example that w is within some small neighbore hood of 100. The diagram for a physical system which might be characterized by an equation in the form of Eq. (6.4) is shown in Fig. 6.6. 1) IIX(t) Ideal Current Generator: 2E cos W t e FIG. 6.6 DIAGRAM FOR A PHYSICAL SYSTEM WHICH MIGHT BE CHARACTERIZED BY AN EQUATION IN THE FORM OF EQ. (6.4).

169 It is assumed here that the circuit diagram in Fig. 6.6 characterizes an electron tube oscillator under the influence of an external signal. In order for this physical system to be characterized by an equation in the form of Eq. (6.4), it must be assumed that grid current can be ignored and that plate current is independent of plate voltage. Granting these assumptions, the grid voltage in Fig. 6.6 will correspond to the independent variable,x, in Eq. (6.4). The current from the ideal current generator in Fig. 6.6 will be, of course, equal to 2E cos c t. (Note e that, contrary to normal practice, E represents a current.) 6.11.1 Selection of the Pertinent Canonical Form. The first problem which arises in the use of the algorithm presented in this chapter is the determination of which of the two possible canonical forms is applicable. This determination is, in turn, as is discussed in Section 6.2, intimately related to the selection of w, the angular frequency of the dominant term in the first approximation to the response. The selection of cU is accomplished by investigating the characteristic roots of the s following linear differential equation: A(z)x - NB(z)x = 0. There are, of course, several methods for finding the characteristic roots of the above equation. For example, the root locus ([21], pages 559-612) technique from control system theory is one possibility. The poles and zeros of the transfer function B) are located as shown in Fig. 6.7. The root locus pattern associated with Fig. 6.7 would be of the form shown in Fig. 6.8. -6 It is clear from Fig. 6.8 that the critical N, N1, is -1000x10O. The corresponding characteristic roots for this value of N are:

170 [m(z) X POLE; (-2 + j102.04) POLE: (-96 + jO)Re(z) ZERO: (0 + jO) POLE: (-2 -.j102.04) X FIG. 6.7 LOCATION OF THE POLES AND ZEROS OF A(z)

171 Im (z) AT N=-1000 X 10-6 THIS POLE HAS SHIFTED TO z=+jIOO N —-00 N DECREASING TOWARD -00 ZERO -\ \~ ~ Re(z) C AT N = -1000 X 106 THIS POLE HAS SHIFTED TO z-I100 N-0 -CO AT N -1000 X 10'6 THIS POLE HAS SHIFTED TO z jlOO' — FIG. 6.8 ROOT LOCUS PATTERN ASSOCIATED WITH FIG. 6.7

172 a= jiO (= j( ) X2 = -jiOO (- -jN ) (6-5) 3 = -100. Since o is assumed to be in some small neighborhood of 100, the e best choice for w is obviously wa = m. Therefore, this example is a case S S e of harmonic response as opposed to subharmonic, ultraharmonic, or ultrasubharmonic response. Furthermore, because aw = at, the forcing function s e e(t) contains a term of angular frequency cu; consequently, the first canonical form is the pertinent canonical form. The characteristic roots of the " " = O" differential equation, A1, X2, % 3, are, from Eq. (2.14) and the paragraph following it, A = jl003, k2 -jlOO0, (6.6) %3 = -100o, where s e where 1 = It is important to appreciate that P can be considered 100 to be a normalized frequency variable. If oe is equal to cNl' P equals 1. If w is greater than ON' then P is greater than 1; and if w is less than e e oN1' then P is less than 1. This point is emphasized because P is used as the frequency variable in much of that which follows. At this point, then, the pertinent canonical form has been selected and certain preliminary calculations have been completed. The remainder of this section is devoted to determining an approximation to the periodic

173 response and approximations to the associated characteristic exponents, In each case the coefficients of the related power series are determined 2 up through the coefficient of the. -term. Thus, if the periodic solutions, p(t, p), are expressed in a power series in., p(t, p) = p(t) + pPl(t) + 2p2(t) +... (6.7) the coefficients determined in this section are p(t), pl(t), and p2(t), and the approximation to the response is p(t, 1) = p(t) + Pl(t) + p2(t). (6.8) Similarly, if the characteristic exponents are expressed in a power series in p., ai() ~ aiO + ail + + (i = 1, 2,..., n), (6.9) the coefficients determined in this section are CaiO, ail, and ai2 (i = I, 2,..., n), and the approximation to the characteristic exponent associated with the periodic response is i(l) = ai + ail + ai2 (i = 1, 2,..., n). (6.10) 6.11.2 First Approximation to the Periodic Response. The first approximation to the periodic response of the nonlinear feedback system which is considered in this example is determined by the methods of Section 6.3. As is discussed in that section, these methods are principally concerned with the solution of the zeroth-order describing function relation, -.... -- lt A t 1 Note that pk(t) = clk(t)e +... + cnk(t) for k = 1, 2,..., and it is periodic with period T.

174 A(jo ) +s n - B( + i- =0 (3.16) 5;;\;s 10 where the solution of the equation, c10, forms the first approximation to the periodic response as follows: jW t -just p(t) = clOs + ce s t=c10 + c10e +j(o t s The term E+ in Eq. (3.16) is the coefficient of the -term in the Fourier series for the forcing function e(t). In the present case, e(t) = 2E cos [e t] (E real) and ca =o s e therefore, E+ = E and E = E. A(jo ) The term B7 in Eq. (3.16) is given by A(ja) -jw3 - 1002 + j 10,800 + 106 s s S S (s sB(j) j(.8 x 106 (6.): -j(o 8 x 1o )W The locus of this term as w varies is sketched in Fig. 6.9. Further, A(jo) for the purposes of this example, B(j~) is sufficiently-well approximated by the following linear expression: A(jos) -6 -6 B(j% ) [-1000 + 250(10Q)]J x 10 - j[250(100z )] x 10, (6.12) where Zb = P - 1 is a measure of the deviation of ws from X = 100; i.e., = (1 + z@)100. The term n0, the zeroth-order describing function, is determined by the following relation:

~~ ~ ~~ 25_ ~25 -20 50 15 60 < XA 100 70 80 5 0 i p03Re A(jws) B(jws) 100 NOTE: FROM EQ (610 A(jws) -ujs5- IOOws2+j 10800cus+ 106 B( js) -j(0.8 x 106)c5 FIG. 6.9 LOCUS OF B(^+) vs. s 2175\~20

176 T/2 j1 t -jwst -jow t O c / F{co dt -T/2 which in this case is T/2 1 -6 s n0 -o T [(-1500 x 10) { C10S + CE -T/2 (6.13) (j t * -jw t 3 -jw t + (6.66 x 10) {cO s +C } ] s t, no = (-1500 x 10'6) + (20 x 10-6) clcl. (6.14) a 1j Letting cl = 2 e, the expression for n becomes no =(-1500 x 10 6) + (5 x 10 6) a. (6.15) 2 A plot of n vs. a is presented in Fig. 6.10. All the functions which appear in Eq. (3.16) are now known explicity and the determination of the desired c10 is now a matter of substituting these functions into Eq. (3.16) and finding the solutions of the resulting equation. The solutions of this equation are determined here through use of the graphical technique illustrated in Fig. 6.4. A typical step in this graphical technique is illustrated in Fig. 6.11 and follows the ensuing pattern: A(jo ) a) The locus of B- vs. s is plotted. b) The locus of n vs. "a" is plotted on the same complex plane 0 A(jw ) that the locus of 7B vs. )s is plotted on. In this case B(j(o) s 5 -6 the locus of no starts at -1500 x 10 for a = 0 and moves toward the right along the real axis as "a" increases. Fig. 2) which corresponds 6.10 gives the value of "a" (actually, a ) which corresponds

a 0 50 1I(0o 150 200 250 LETTING CIO= yE no BECOM -5 00 i FIG. 6.10 no, ZEROTH-ORDER DESCRIBING FUNCTION, vs. (5) ^ n0=(- 1500x10Y)+(20 x10)C10C,0* LETTING Cio= -l-E, no BECOMES -1000 ~~~ ~ ^^^ ~~~ no(a)=(-1500x1l6)+(5 xI0'6) a2 -1500~ FIG. 6.10 no, ZEROTH-ORDER DESCRIBING FUNCTION, vs. a2

178 to a given point on the locus of n. c) A value is assumed for E, say E = 2 x 103. d) A value is assumed for "a", say a = 12. E+S e) The magnitude of — is calculated: clO =333 x 10 c10 a f) The point on the locus of n vs. "a" which corresponds to a = 12 is determined: from Fig. 6.10, n (12) = -780 x 106. A(jw ) g) The intersections with the B(j -locus of a vector which S-6 E+S originates at n (12) = -780 x 10 and is of length E = -6 10 333 x 10 determine the possible values for we (recall that ) = a ). It can be seen from Fig. 6.11 that these values S e are w = 99.61 and c = 101.28. e e h) The vectors in g) also determine the angle 9. This can be seen by noting that the vectors are determined by 2E - jQ A(jw ) a B(jWOs) o therefore, the angle of each vector is equal to the negative of its associated 9. Then, from Fig. 6.11 it can be seen that the angle associated with we = 99.61 is 9 = 200, and the angle associated with ao = 101.28 is 63. e The above procedure is repeated for new values of "a" so that the behavior -3 of "a" vs. au and 9 vs. ) for constant E, in this case E = 2 x 10, can e e be determined. A new E is then selected, and the procedure is repeated. The results of the use of this graphical technique are presented in Figs. 6.12, 6.13, 6.14, 6.15, and 6.16. Additional information which pertains to later parts of this example is presented in Fig. 6.12.

WOs 98 \^~~~~~~~~~ A A(jWs) LOCUS OF 106 A(jC.) AS Wos VARIES Ws=99 \ 250 \ 106 x- 0 Rn WS 01 -250 102~~~~~~16R ^ ^\\ \/ /,06 p, A(J~~~~~~~~~s)0 ^\ \^ -OO^^^^ 10^ Re n ^ LOCUS OF IO6noAS "a" VARIES^ 630 ^ -Q> OJ^IOI \^ \ ~~~~~10-20 FIG. 6.11 ILLUSTRATION OF THE GRAPHICAL TECHNIQUE EMPLOYED IN THE SOLUTION OF THE ZEROTH-ORDER DESCRIBING FUNCTION RELATION

I-~~..~ 1~111-.-I..~... -i..- —:- NOTE: POINTS DESIGNATED - KI~"^ ~~~~-1 -~-~~ t -" t~~~11 t 1tXt. —N ——... CORRESPOND TO THE - 7h[Tt~~z-:1GENERATING SOLUTIONS FOR WHICH.T +1XWX77.44~ in - 7 HIGHER-ORDER APPROXIMATIONS ARE-T - 7~~ — ^ ~~ — ~. T~ ~- 1 — 7~11'~:~ EVALUATED. -- 7|-t ~I~~0 ------— ~-~.l-lLFIr ~ I { - E I - - ~ I2 -—;A1M S ~1 1- - j%7: —:-K~7~.! I I ~ E =2xl0"^ AMPS —-F -~ ~^~^ — -- 7-~ —-t -1-41~ -~1 -jt i'"^"'1 T- t- =_....... 7 aI 11 OR a2, EQUALS ZERO t - 71- 4_I E I xIO AMPS - I t j j v: -.:.F-:'..... ~,~i.1i=_.!j {~-t~t - JACOBIAN EXPRESSION.7f-_:_ * I I t 1:,LESS THAN ZERO 7 LOCUS OF POINTS AT WHICH BO or -TH-....\.. Re all AND Re a-, EQUAL ZERO R I - - - - ~ --— +-.... -.. —.7' — JACOBIAN EXPRESSION GREATER THAN ZERO -4 7' —"-' E = 0.7 x 0 AMPS!::-' - LOCUS OF THE POINTS ATI~}-"~ RESPONSE vs w.: LOCUS OF THE POINTS AT WHICH~ LWHICH THE JACOBIAN EXP O0S: L S OF TE P S AT W H E R OR BH EXPRESSION VANISHES:LOCUS OF POINTS AT W.ICH BOTH Rea,, AND Rea2 VANISH 777~ —-L-.-~j~1.-^- - -LOCUS OF POINTS AT WHICH BOTH Re On ~ o..~~~.......-~-~ —~- ~~-._.AND Re 021 EQUAL ZERO,__ x lO.~...^ A M P -..... " -'..- -...._ —........:......... 9 10 1 102 103' —FI.6.12:AMPLITUDE,: ":'". THE:''............... E= 0.6x IO- AMPS a -7, - v............'"." FIIG. 6.12 AMPLITUDE, ".,OF THE FIRST APPROX IMAT ION TO THE PERIODIC RESPONSE vs ws: LOCUS OF THE POINTS AT WHICH THE JACOBIAN EXPRESSION VANISHES: LOCUS OF THE POINTS AT WHICH EITHER OR BOTH all AND ac, VANISH: LOCUS OF POINTS AT WHICH BOTH Real, AND Rea( VANISH

360~ w LU) w ~~ 270~ 30-^~~~~~~~~~~~~~~~ ~18090 90~ 00 97 98 99 100 101 102 103 We (RADIANS/SECOND) FIG. 6.13 PHASE ANGLE, e, OF THE FIRST APPROXIMATION TO THE PERIODIC RESPONSE vs. w,:E2xlO3

3600 w w O3 LC E=lx 10-3 1800 ro 90" ~~~~~~L Xt I~~~~o00 97 98 99 100 101 102 103 we (RADIANS/SECOND) - FIG. 6.14 PHASE ANGLE, 8,OFTHE FIRST APPROXIMATION TO THE PERIODIC RESPONSE vs. w:E=lxlO3

1800 E = 0.7 x 10-3 (DO00 ~ ~ ~ ~ ~ a9~- 9 90~ 97 98 99 100 101 102 1 we (RADIANS/SECOND) - ~~_____________ 615 PHAI -450 __________________ ___________ FIG. 6.15 PHASE ANGLE, 8, OF THE FIRST APPROXIMATION TO THE PERIODIC RESPONSE vs. ws: E=0.7 xlO3

3600 w L^^^~~~ \4~~ ~LUL 270~ E = 0.6 x 10-3 180 900 37 98 99 100 101 102 103 We (RADIANS/SECOND) FIG. 6.16 PHASE ANGLE, 8.OF THE FIRST APPROXIMATION TO THE PERIODIC RESPONSE vs. w E=0.6x 13

185 The foregoing completes the determination of the first approximation to the periodic response. The next problem is the determination of the first approximations to the stability criteria. 6.11.3 First Approximations to the Stability Criteria. The first approximations to the stability criteria are obtained by employing the methods discussed in Section 6.4. According to Section 6.4, the first approximate stability criterion is that the Jacobian expression should be positive; that is, A(jo)) A(-jsw ) +2s 25 (6.16)'o" B'jw Lo B(-Jw ) o So >0 where SO and S 2s are defined in Section 6.4. Referring to Eq. (6.3), o o it can be seen that dF (-1500 x 106) + 3(6.66 x 106) x dx Consequently, S = (-1500 x 106) + (40 x 106 )clc-6- = (-1500 x 110 6)a2 -6 - 2 j26 (6.7) S+ =(20 x 10)c0 =(5 x 10 )a e, and (617) S-2s: (20 x 10-6)C0 (5 x lO-6)a c,- j2 o 10 So = (20 x 106)cl = (5 x 10 )ae2 Substituting Eqs. (6.17) and (6.12) into the inequality (6.16), gives the following relation: 1012{[-500 + lOa2 - 250(100YL) ] + 250(100t)2 - 25a}> 0 (6.18) where (100O) = (W -100). The locus of those points at which this Jacobian s expression vanishes is plotted in Fig. 6.12.. Note that this locus passes through the points on the "a" vs. )s curves at which d~ becomes infinite. s

186 This is just as expected. from the discussion in Chapter III of the Jacobian condition. Further, note that the condition (6.18) is not satisfied inside the closed curve formed by this locus. Finally, it can be seen from Eq. (5.31) that points on the locus are points at which all and/or Ca21 equals zero. The second approximate stability criterion presented in Section 6.4 is A(jw ) Re{v13 B(j) [SO j ]} In this example, condition (5.34) reduces to x 10 (500 - lOa) < 0, (6.19) where vl3 is, from Eqs. (E2) and (E3), 1 -4. (6.20) v13() = (1 + j) x104 (6.20) 1~342 Similarly, 2 () = -— ( ) x 10l; (6.21) v33 () 1 x 10 (6.22) 33 2 22 It is readily seen that condition (6.19) further reduces to a > +/0 7.0 7. (6.23) The boundary between regions in which condition (6.23) is satisfied and regions in which it is not satisfied is shown in Fig. 6.12 as a horizontal line through a = /~0. It can be seen from Eq. (5.31) that points on this line at which the Jacobian expression is positive are points at which Re all and Re a21 both equal zero. The first approximation to the third characteristic exponent,

187 a 3(1), associated with a periodic response of the system considered in this example is determined from Eqs. (5.24) and (5.49) (see Section 6.4). Thus, from Eq. (5.24), a 30= -100, (6.24) and from Eqs. (5.49) and (6.22), 40 -4 2 (6.25) a31 x 10 {-500 4-+ la2 + 250(100)}. (6.25) The corresponding approximate stability criterion is Re (C30 + 31) < 0, (6.26 and, clearly, the boundary between stable and unstable regions is Re (a30 + 31) = 0. (6.27) Substituting Eqs. (6.24) and (6.25) into Eq. (6.27) gives -lO0 + - x 104 {-500 + lOa2 + 250(100t,)} 0. (6.28) For small B3 and P 1, the solution to Eq. 6.28 is well-approximated by a = 2550 + 2500A. (6.29) It can be seen from Eq. (6.29) that "a" must be quite large relative to the range of "a" considered in Fig. 6.12 before the first approximation to the third characteristic exponent passes into the right half-plane. At this point the entire first approximation is completed; that is, first approximations have been obtained to the periodic response and all three stability criteria (characteristic exponents). The next problem is to determine higher-order approximations to the periodic response.

188 6.11.4 Higher-Order Approximations to the Periodic Response. The higher-order approximations to the periodic response are determined by employing the methods of Section 6.7. In order to facilitate the employment of these methods, a set of new symbols is introduced. Recall that the periodic solution which is being determined in a power series in A is denoted p(t, p). Further, the power series expansion is represented as 2 p(t, p) = p(t) + p+ pl2(t) +) +... (6.7) The first approximation and each of the pk(t)'s are a periodic function of t with period T. Expanding each Pk(t) in a Fourier series gives 00 j Po() I I (k = 1, 2,...). (6.30) Pk(t) = Pk o=-o The Fourier coefficients, Pk' are employed extensively in what follows. The symbol pk designates the coefficient of the term in Eq. (6.30) for -S 3s which j~oy = ju. The symbols k, pk, etc. follow the same pattern. The first term in the power series shown in Eq. (6.7) is, of course, jWst -js t p(t) = C + c 10 + c10 The next few paragraphs are devoted to the determination of Pl(t). Determination of pl(t): The function Pl(t) is related to the cjl(t)'s (j = 1, 2, 3) as follows: jW t -jo t ^ t Pl(t) = Cll(t)e + (t)(t)e 3, (6.31) or p1(t) = cll(tjlOO t + c2l(t)3jlOO1t + 1(t) -10t. (6.32)

189 The first step in the determination of the functions cll(t), c21(t), and c31(t) is the carrying-out of the integration indicated in Eq. (6.1). As is pointed out in Section 6.7, the proper selection of c10 implies that the integrand in Eq. (6.1) contains no terms of angular frequency.s Therefore, Eq. (6.1) reduces to p-K.s jw s -jCw s Cjl(t) = ajl + 3 j C B(z) Fo(s) S- -F_ dis. jl J3 o 0 0 (6.33) (j = 1, 2, 3) For this example, ~ J) s -j) s ^ j3 S *j33 s [F(s) - F e F e ] = 6.66 x 106 [c 0 s + cO vj (P) (j = 1, 2, 3) are given in Eqs. (6.20), (6.21), and (6.22); and J3 B(z) is given in Eq. (6.3). Substituting these expressions into Eq. (6.33) and carrying out the indicated integration gives the following expressions for cll(t), c21(t), and c31(t) -11o) + ) x 4 ) x l0 s j 0o 2.06 4 3 j4o st c11(t) c=C11 + (1+j) x 10l c 0 (6.34).4 3g - j~st + Ll (1-j) x 10 c e c 10 2 64 3 4 j4~ t c 1(t) 825 x 10 -5(3+j)c0 -3 j3 1 2 2 x (6.36) ~8.a5 -3-(%3+ j3ws)t c3 8(t) 2 x 10 (3-j) cO E p 1

190 In Eqs. (6.34) and (6.35), the mean values of cll(t) and c21(t), cll and c21, are, as yet, undertermined, and must be determined by the nl-describingfunction relation. (Recall that these mean values being undetermined is equivalent to the initial conditions, all and a21, being undetermined.) The initial condition, a31, has been chosen for c 31(t) so that c31(t) contains no constant term. The nl-describing-function relations are as follows: A(j(w ) n - B s 0 (6.37) s where F+S n F1 (6.38) 01l In this example, F1 = [(-1500 x 10-6) + (20 x 10 6) p(t)] pl(t); (6.39) consequently, nl = [(-1500 x 106) + (40 x 106)c0c] + (20 x 10-6) c 11 C11 -6 2 3s (6.40) 3$ + (20 x 106) co - ~ Note that a knowledge of p1 is implicit by Eqs. (6.34), (6.35), (6.36), 5 0 -s 0* 0 and (6.31). It should also be realized that pi = cll and pi =(C11) c21. Substituting Eqs. (6.40) and(6.12) into Eq. (6.37) and simultaneously considering the complex conjugate of the resulting equation gives the following linear inhomogeneous system of equations: [{-500-250(100lo )+ 40 cAcO + j20(10ooB)} x 10o6 ]c. + [(20x 0o )c0.] (c.1) - (20 x 10-6) cp3 (6.41) (Second equation of this system appears on the following page.)

191 [(20 x 10 6)c0 ] C11 + [ -500-250(100P) + 40coco0 - J250(100nsP) xlO(c] (c)*.6 2 -2 3 - (20 x 106) co p1-3 A comparison with Eq. (6.18) shows that the determinant of the above system of equations is equivalent to the Jacobian expression. 3s -3s Substituting the expressions for p3 and pP3s which are obtained from Eqs. (6.34), (6.35), and (6.36) into Eq. (6.41) and solving the resulting system for cl gives, ~ c, (3.84xo10)a4[ {100+33.5(100P)-.a} +j { 33.5+66.6(100B)-a}] (6.42) 11 10 22 2 2 4 [-5oo-20(100oo) + lOa ] + (250) (10ooB) - 25a a JQ where, again, cl O = A. The above equation completes the determination 2hraa 10 2 of pl(t). It is interesting to note that cl exists as long as, once again, the Jacobian condition is satisfied. Moreover, there is clearly only one 0 possible cll. Evaluation of Eqs. (6.34), (6.35), (6.36), and (6.42) for the three sample points indicated in Fig. 6.12 gives the following expressions for p(t) + ~Pl(t): Point 1, p(t) + uPl(t) ~ 4 cos [101.8t + 0.156l ] (6.43) + {o0.93 x 10-4 cos[101.8t + 0.395n] + 9.8 x 10-4cos[305.4t + 0.4783]}; Point 2, p(t) + Pli(t) ~ 7.5 cos[lOlt + 0.156i ] (6.44) + i {462 x 10 4cos[101t + 0.305c] + 66.4 x 104cos [303t - o.566t] };

192 Point 3, p(t) + 4iPl(t) 12.75 cos[100t] (6.45) + u {-320 x 10-4cos[l00t - 0.64t] + 334 x 10-4cos[300t + 0.100]}. The next few paragraphs are devoted to the determination of p2(t). Determination of p (t): The determination of p2(t) follows exactly the same pattern followed in the determination of Pl(t). Equation (6.33) is replaced by the corresponding equation for the cj2(t)'s. The mean value of cl2(t), cl2 is determined by the n -describing-function relation. The system which results from this describing-function relation is, just as Eq. (6.41), a linear inhomogeneous system whose determinant is the Jacobian expression, Carrying out the indicated computations, the expressions for p2(t) corresponding to the three sample points are as follows: Point 1, p2(t) " 19.48 x 10 cos[101.8t + 0.268r] -8 + 147.0 x 10 cos[305.4t - 0.606n] (6.46) + 2.92 x 10cos[509.0t + 0.687i]; Point 2, P2(t) = 22.3 x 10 4cos[101t + 0.8659] + 1.54 x 10 4cos[303t - 0.596n3 (6.47) + 0.058 x 10 cos[505t + 0.944n3; Point 3, p2(t) 1.76 x 10 cos[100t + 0.815 ] - 0.812 x 10 cos [300t (6.48) + 0.864 x 104 cos[500t + 0.1717]

193 The expressions for p(t) + tpl(t) + A, p2(t) can be obtained through a combination of Eqs. (6.43), (6.44), and (6.45) with Eqs. (6.46), (6.47), and (6.48). The one remaining problem is the determination of the higher-order approximations to the stability criteria (characteristic exponents). 6.11.5 Higher-Order Approximations to the Stability Criteria. As was stated previously, the characteristic exponents, ai(l) (i = 1, 2, 3), 2 are evaluated here up through the 2-term in their power series representation. The first two terms in these series, aio and ail (i = 1, 2, 3), can be determined by the methods of Section 6.4. In particular, for a. a10 = +jP100, 20 = -j100, (6.49) a 3= -P100; and, for ail + 2 x 104 l 500-10Oa2)2 -2 {[-500-250(100P)+l0a2] +2502(100B)225a},} {20 1 04 (500 - 102) 65 r20 -4 21 2 (6.51) 20 -4 ~l00 _10a2 2 2,2 2 4 20l104 /500-10a2)2 2 {[-500-250(100B)+10a2] +2502(100B) -25a}, and 31'-x 10 a31 x 10 { -500 + lOa + 250(100)}. (62) The above ai0's and ail's (i = 1, 2, 3) are evaluated for the three sample solutions, 1, 2, and 3 for which the higher-order approximations have been evaluated. The results of this evaluation are as follows: io Ji1 sapesolutio n s 1.~ 2~ and 3.h for wi c t h e nv hn-ige-ore apoimat ions hav

194 Point 1, ca10 + Cla 0 j0l1.8 + 4[0.67 + j2.43] a20 + 4aC21 -jlO1.8 + [0.67 - j2.43] (6.52) a30 + 4a31 o -101.8 + 4[0.43] Point 2, a10 + all j101 + 4[-0.12 + j0.37] a20 + a21 j-jlO1 + 4[-0.12 - j0.37] (6.53) a30 + a 31 -101 + 4[1.23] Point 3, a10 + all a j100 + 4[-1.72] a20 + a 21 -j100 + 4[-2.72] (6.54) a30 + a31 -100 + [4.44] The next terms in the power series representations for the a.(l)'s (i = 1, 2, 3) can be determined by the methods of Section 6.8. For example, the determination of a1 follows the ensuing outline. 120(0) a) The term Ul0( is determined from Eq. (5.35). In this case, u110(0) o120(o) 4(i+j)x104 20 "4(202io6 20(0) - ): { all — xlO'4 [o0a (5 00+j(-500 -500(100) + lOa2)]}. (6.55) U110 800 co~ 2 sO b) The terms r10(t) and r2(t) are determined by Eq. (5.20). In this example, these terms are

195 rJo(t) ({2 B + j [1.92 - 2.0 - 80 So ]}x 10o -{ ji240 S+2S 6 X10} e s (6.56) W-J2w t +{ j80 So2 x16 }e s } r, ) r +2s i2wo t -j2wt t r (t 0 +r C2s E S + I'-2s C s (6-57) 10 10 10 10 and r20(t) { 2 - j [1.98P - 2.0 - 80o SO]} x 106 -{j80o 2 x 106 } (6.58) + {j24 S2 x 10}e 6 -s j2wst -2s -jt t r (t r r2 +r (6.59) r(^ c~ - ^'+ r Co 20 (6 59) 20 20 20 c) The foregoing terms are substituted into Eqs. (5.36), (5.37), and (5.38), which partially determine ul2(t), and completely determine ulll(t) and u131(t). These substitutions result in functions of the form j2cu t -j2o t -j4 t / 0^\t o 2s s -2s s -4s s Ull(t) = u'll + U lle lsl + u + (6.60) ^111^^ "^111 ^11 11 1 11 11 j2w t -j2w t 0, o 2s s -2s s -4s' - j t u121(t) =U121 +121 e + u 21 + U11 E (6.61) and j2w t 2 t t o 2s s -2s s - 4s -j4 ~t U131(t): U131 + u131 ~ + Ul31 ~ + Ul31 ~ (6.62) where ul21 is undetermined. An example of a typical expression for one of the above Fourier coefficients is

196 -2s u20 (0) o +2s F 120......1 (0) vlo + u10(0) r20 - urO u20() r2e 112.u ~~~~~,+ =io() -. (6.63) uill ullo(O)v-3 j0 3~ P - S o T+2s -2s d) The terms S, S, and S1 are evaluated as in Paragraph d) of Section 6.8. For this example, = 40 x 10 [clO(C)* + c0 c1] (6.64) S+2s = 40 x 106 c1O ll + cO 1 (6.65) S1 [ o1 10 3 and S = 40 x 10 [c(c)* + c ] (6.66) e) a12 is determined by substituting the terms evaluated above into Eq. (5.68). A similar procedure is carried out for a21 and a31' The results for the sample points 1, 2, and 3 are as follows: Point 1, a2 0.014 + jO.032 a22 0.014 - j0.032 (6.67) a32 -0.014' 32 Point 2 al2 -0.027 - j0.074 022 2 -0.027 + j0.074 -0.011(6.68) a32X -0.011

197 Point 3 ca2 0.16 a22 O 0.23 (6.69) a32~ 0.076 In each case, the approximation to the characteristic exponent, ai(l) (i = 1, 2, 3) is ai(l) = aiO + ail + i2

CHAPTER VII SUMMARY AND CONCLUSIONS 7.1 Introduction A new method for utilizing the techniques of perturbation theory for the analysis of a class of nth-order nonlinear feedback systems under the influence of a general periodic forcing function has been developed. This method is suitable for the determination of periodic responses and the investigation of whether or not these periodic responses. are asymptotically stable. In addition, this method can be extended so that it is useful for the analysis of nonlinear systems which are even more general than those considered in this study. 7.2 Conclusions A. It is pointed out in Chapter I that there are roughly two bodies of literature which pertain to the problem which is of interest in this study. The first contains writings by those authors (Tucker, Smirnova, Adler, Hunton and Weiss) who analyze certain physical systems by methods based on "averaging" or equivalent linearization and plausibility arguments (see the outline of Tucker's argument in Chapter I). The systems studied by these authors form a subclass of the class of nonlinear systems considered in this study. The difficulties with these methods are several: their development does not clarify the sense in which the solutions obtained are approximations, nor do they admit to a refinement of the approximation if this is desired; the stability criteria are limited to an approximation of just two characteristic exponents out of the n which are needed for a complete description of the stability or 198

199 instability of an nth-order system; moreover, there is no obvious manner in which to refine the approximation of these characteristic exponents; finally, these methods cannot be extended to more general nonlinear systems with much consistency or confidence. The second body of literature contains writings by those authors (Poincare, Coddington and Levinson, Cesari, and many others) who discuss the application of perturbation techniques to the analysis of differential equations. As is pointed out in Chapter II, these authors usually stress the "normal perturbation problem" at the expense of the "inverse perturbation problem", which is a satisfactory restatement of the differential equation which characterizes the nonlinear feedback system as a problem amenable to solution by perturbation techniques. They are more interested in developing the mathematical tools of perturbation theory than in the application of these techniques to the analysis of nonlinear feedback systems. B. It is shown in Chapter II that the "inverse perturbation problem" is the central issue in any application of perturbation techniques to the analysis of nonlinear systems. The essentially synthetic character of this problem is pointed out, and it is emphasized that there is no unique solution to the "inverse perturbation problem." Moreover, the mutual relation between the original and the "I = 0" differential equations is stressed. C. A new method for obtaining a partial solution of the "inverse perturbation problem" associated with the application of perturbation techniques to the class of nonlinear systems considered in this study is presented in the last two sections of Chapter II. This method

200 involves a division of the original differential equation into two parts: (1) a part which if equated to zero would constitute a linear differential equation that possesses a periodic solution with period commensurate with that of the forcing function, e(t); and (2) a part which comprises the difference between the preceding part and the original differential equation. The second part of this division is multiplied by the artificially introduced parameter, i, so that at, = 0 the equation reduces to a linear differential equation, and at u = 1 the equation becomes the original differential equation. The above division of the original differential equation is carried out so that it allows the entire analysis which follows to be conducted in terms of easily identifiable expressions from the nonlinear feedback system. For example, the numerator and denominator of the transfer function for the linear network, B(z) and A(z), are, except for new arguments in certain cases, retained intact throughout the analysis. This is true also for the transfer function of the nonlinear element, F(x), and the forcing function, e(t). The principle advantage of this restatement of the original differential equation is that the steps in the recursion procedures for the solutions and the characteristic exponents reduce to relatively simple operations involving describing-function-type relations. In fact, it is reasonable to consider the work in this study partly as a generalization of the describing function technique. Another advantage of this method of restatement is that by introducing two canonical forms it allows, as is clarified in Chapter III, a wider range of forcing functions, e(t), to be considered.

201 D. In Chapter III, equations which determine the first approximations to the true periodic responses of the nonlinear system are developed by application of perturbation techniques to the two canonical forms. Because of the special choice in Chapter II for the canonical forms, the equations which determine the first approximations are in a form which is similar to the normal describing function relation. This is especially true for the first canonical form. However, these equations can be applied to a wider range of nonlinear feedback systems than can the normal describing function. Further, since these equations are the first step in a recursion procedure which determines a higher order approximation with each step, the approximation determined by these equations can be refined at will. This is not the case for the describing-function-type relation as it is usually developed; e.g., Tucker [3]. E. In the latter part of Chapter III, the usual Jacobian existence condition from perturbation theory is shown to be equivalent, except for the special case in which a certain determinant vanishes [Eq. (3.39) for the first canonical form, and Eq. (3.45) for the second canonical form], to requiring that the slope of the frequency of the first approximation versus its amplitude be nonzero. At those points where this slope does become zero the first approximation is not continuable. Furthermore, it is shown in Chapter V that there is a relation between this Jacobian existence condition and one of the approximate stability criteria. Therefore, the relation between the aforementioned slope and the Jacobian existence condition is doubly important: (1) because it allows the Jacobian condition to be readily applied, and (2) because it

202 simplifies one of the approximate stability criteria. F. In Chapter IV, the application of perturbation techniques to the analysis of the class of nth-order nonlinear systems which is of interest in this study is continued. A recursion procedure is developed whereby the periodic responses, p(t,p,), of the nonlinear system are obtained as a power series in the parameter,, as follows: First Canonical Form, n co t p(tp,) = C z c (t)cl t (7.1) J=1 k=O Second Canonical Form, n Xo k h t p(tF) - ~ cj(t) ke + h(t; ) (7.2) j=1 k=O where h(t;s ) is a periodic function with period T and the cjk(t)'s satisfy the following "periodicity" conditions XT = cjk(t+T)e - Cj(t) { } (7.3) k = O,1,... The actual responses of the nonlinear feedback system correspond to the p(t,,u)'s with u. set equal to one, and the first approximations to the responses obtained from the equations developed in Chapter III correspond to the p(t,n)'s with p. set equal to zero. As a part of this recursion procedure, a hierarchy of describingfunction-type relations are defined. The first or zeroth order relation of this hierarchy for either canonical form is just the appropriate equation from Chapter III. The succeeding members of this hierarchy are each

203 associated with a particular one of the steps of the recursion procedure. The advantage of this hierarchy is that it simplifies the recursion procedure and places the usual equivalent linearization approaches (i.e., zeroth order describing function relation) in a context instead of allowing the equivalent linearization approximations to be an end in themselves. G. The method developed in Chapters II, III, and IV allows the periodic responses of the nonlinear system to be determined by perturbation techniques as a power series in the artificially introduced parameter, A. In Chapter V, perturbation techniques are utilized to discover whether or not these periodic responses are asymptotically stable. An important feature of the perturbation techniques employed in' Chapter V is that they harmonize with those used in Chapters II, III, and IV. It is, therefore, possible to maintain a consistent point of view throughout the entire analysis. As each step in the recursion procedure for the determination of the periodic solution, p(t,,u), is carried out, the corresponding step in the recursion procedure for the determination of the characteristic exponents, CaP(L) (i = l,...,n) can be carried out. This allows the estimate of the characteristic exponents to be refined at the same time that the approximate solutions themselves are being refined Again, because of the choice of canonical forms in Chapter II, the expressions which occur in the determination of the characteristic exponents are simple, recognizable, and easily related to the nonlinear feedback system. This is particularly true for the expressions which determine the o *'s and the ails (i = l,...,n). Since these two terms are often sufficient to answer the yes-or-no part of the stability

204 question, this is an important simplification. Another important advantage that this method of investigating the stability question has is that it resolves the two usual difficulties present in stability criteria which are based on equivalent linearization arguments. Since a periodic solution whose stability is being investigated is known only approximately, its associated variational equation is known only approximately. Therefore, even if the characteristic exponents of the resulting variational equation could be determined exactly, they would still be approximate in relation to the true periodic solution. But, even if the variational equation were known exactly, there would still be the problem of carrying out some type of approximate analysis for the determination of its characteristic exponents. The perturbation procedure used in Chapter V carries out both of the necessary approximations —variational equation and the characteristic exponents thereof —simultaneously, and allows the estimate of the characteristic exponents to be refined at will. The usual stability criteria developed from equivalent linearization arguments do not have this flexibility. It is pointed out in the course of Chapter V that there is a relation between the Jacobian condition, and,1. and a21. In particular, it is shown that ala21 > 0 if the Jacobian expression is positive and a a ll 0 if the Jacobian expression is negative. Moreover, assuming that the equation, (3.16) or (3.29), which determines the first approximation, p(t), has a solution of the form [u) = W (a),e = e(a)], al a2 = do! 0 if - = 0. Since the a.'s and the a.'s are often sufficient for the da io il formulation of suitable stability criteria, and, further, since cl(#) and a2(p) are often (but not always) the most likely characteristic exponents to cross into the right-hand plane [recall that cal() and c2(G)

205 are both on the imaginary axis at L = 0 and that the other ai(p.)'s (i = 3,...,n) start in the left-hand plane], the relation between all and a21 and the Jacobian condition is especially useful and important. H. There are several disadvantages and limitations associated with the methods employed and developed in this study. First, there are two related and unsolved problems: (1) the determination of the allowable 4-intervals, and (2) an estimation of the error resulting from approximating the a.i(L)'s with a finite number of terms from their power series expansion in u. It is always possible to circumvent these two problems by considering only "weakly-nonlinear" systems —that is, systems which are almost linear. On the other hand, there are many "strongly-nonlinear" systems which can be treated by the methods developed in this study. Consequently, a solution to these two problems which would designate which "strongly-nonlinear" systems could be considered would be of great importance. Another, less serious difficulty is that, in Chapter V, the characteristic multipliers are assumed to be distinct for 0 < P < 1, q )1()T (:x2(1)T except for e and e which coalesce, as L -*0, at 0. In those rare cases in which it is absolutely necessary to consider multiple characteristic multipliers, special methods are available (see Moulton [16], pages 331-348). A limitation of the methods developed in this study is that they do not, as presented, apply to servomechanisms. This comes about because the linear network in the nonlinear feedback system has been restricted thus: B = O. However, this limitation can be overcome by an extension of the methods of this study, which is outlined in the following section of this chapter.

206 As is mentioned in Chapters II and VI, there is still the problem of determining the period, Ts, of the dominant term in the periodic response of the nonlinear system. Finally, it is true that the determination of the higher-order approximations for both the solution and the characteristic exponents rapidly becomes laborious as k becomes larger than two. 7.3 Suggestionsfor Future Research Since the methods for analyzing nonlinear systems developed in this study are explicitly based upon the theory of perturbations and not upon approximations which arise from merely "averaging over a cycle," it is possible to extend these methods in a consistent manner in order that they can be applied to a class of nonlinear systems which is even more general than that considered in this study. Some of the possible extensions are as follows: A. It was pointed out in Chapter I that the limitation of the transfer function, F(x), of the nonlinear element to being just a function of x ruled out the consideration of some important nonlinear physical phenomena; e.g., hysteresis effects, back-lash, etc. Therefore, the extension of the methods of this study to nonlinear feedback systems which contain nonlinear elements which are characterized by a function of the form, F(xx,...,x(k)) (x()_ dx), is of considerable practical dtk importance. Fortunately, this extension does not appear to offer serious difficulties as long as the term, B(z)F(x,,...,x(k)), contains no derivative of order higher than (n-l). In cases in which this is not true, the extension could become difficult, especially if the derivatives of order n appeared in terms of degree higher than one. In any event, it

207 is to be expected that the describing functions associated with such nonlinearities may be complex and frequency-dependent. This is in contrast to the present study, in which only the describing function (zeroth order) associated with the second canonical form is complex and frequencydependent. B. In addition to extending the methods of this study to include nonlinearities which depend on the time-derivatives of x as well as on x itself, it should be possible to extend these methods to other, more complicated systems. A system containing more than one nonlinearity could be treated. In a vacuum-tube oscillator the grid-circuit nonlinearity as well as that of the plate circuit could be considered. In many problems of this type it should be possible, by following the pattern of Chapter II, to recast the differential equations which characterize these nonlinear systems as problems amenable to treatment by perturbation techniques. Moreover, it should be possible to carry out the operations of the perturbation analysis in terms of quantities which are easily related to the actual nonlinear system. In other words, it should be possible to continue the generalization of the describing-function approach which has been commenced in this study. C. If e(t) - 0, the resulting nonlinear system cannot be analyzed with either the methods developed for the first canonical form or those developed for the second. In both these cases the Jacobian condition is not satisfied. As is well known (Truxal [21], pages 601-611; Coddington and Levinson [10], pages 364-369), in such systems not only the solutions but the period of undisturbed oscillation, To, must be expanded in a power series of ~ and obtained by a recursion procedure.

208 It would be worthwhile if a recursion procedure possessing a hierarchy of "describing function relations" analogous to that developed in Chapter IV were developed for these homogeneous systems. A classic and important example of such homogeneous systems is the self-excited oscillator. D. Since the methods of this study apply to any periodic forcing function, it may be that by considering a random forcing function to be a limit of a periodic forcing function as its period approaches infinity the responses of nonlinear feedback systems under the influence of noise can be analyzed with the methods of this study. E. In this study only a special class of responses has been considered, i.e., periodic. It should be possible to combine various perturbation techniques to handle, among others, almost-periodic responses. For example, the solutions could be made up of two parts: (1) a idependent part periodic in t with a period which is commensurate with that of the forcing function, and (2) a 4-dependent part periodic in t with a 4-dependent period. F. As has already been pointed out, the methods of this study cannot be applied, without modification, to a low-pass servomechanism. This limitation of the methods is a consequence of the restriction on the linear network that = 0. This restriction was required in order that the necessity of considering constant terms (e.g., dynamic shifts in the operating point) in the approximate solutions could be eliminated. However, the possibility of such constant terms must be admitted in low-pass servomechanisms and other nonlinear feedback systems for which AA 0.

209 Therefore, it is apparent that the type of generating solution or first approximation which should be sought is a constant term plus a sinusoidal term. For example, such a generating solution for a canonical form similar to the first canonical form of this study would be Jo t -Jwo t p(t) = c10C + cOe +e (7.4) where Qo is an, as yet, undetermined constant. Just as the allowable clo's would have to be determined as part of the analysis, so would the proper Qo have to be determined in the course of the analysis. A possible " = 0" differential equation corresponding to a generating solution such as Eq. (7.4) would be [see Eq. (2.16)] znx+ n [C(Z) - NB(l)]x - (7.5) n[ c() - NB(O)], and a corresponding L-dependent family of differential equations could be [see Eq. (2.17)] nx + n[C(^) - NjB()]x = i[C(O) - NjB()] Q + 4[n {c(l) - NjB(I)}x -{C(z) - NB(z}x (7.6) + B(z){F(x) - N.x} + B(z)e(t) - n{c(0) - NjB(o)}Q]. Then, using the following coordinate transformation:

210 % t Xnt x e rQ 1 Q o Xt t 1e n r 2 0 _ + (7.7) x(-i) (n-1) e.. (n-) l)e r 1 n n the following canonical form could be obtained: -k t rl = Vln g(trl, r1Q ).* * *.. *. * * * * * * * (7.8) -X t Kn = n ne n g(t,rl,.''rnQo) where \t t g(t,rl,...,rnQ) = B(z)F{r,e +... + rnn+ Q0} n t -_ r.A(z)e -A(O)QO +B(z)e(t) j=l Then, following a line of reasoning similar to that used to justify Eq. (3.6), the following necessary condition can be obtained for generating solutions: T -jwt t f e g(t,Cl0,clo00,..,O,Qo)dt = 0. (7.9) 0 This necessary condition is a single relation between the two unknowns, c10 and Qo, and is consequently, not sufficient for the determination of c^ and Q. Clearly, some other relation must be found between c10 and Q~

211 But, purely on the basis of requiring that the generating solution be continuable, any generating solution, Eq. (7.4), satisfying Eq. (7.9) and a Jacobian condition is satisfactory. Therefore, continuability is not enough for the determination of c10 and Qo, and, because the inverse perturbation problem does not have a unique solution, the choice of a second relation between c10 and Q is somewhat arbitrary. However, the most straightforward choice is T f g(t,cc10o,O,c 10..,O,Qo)dt = 0. (7.10) However, this is certainly not the only choice. Given the preceding equations as a start it should be possible to construct a set of relations for this new canonical form analogous to those constructed in Chapters III, IV, and V for the two canonical forms considered in this study. It should be noted that Eq. (7.10) indirectly raises an important point in regard to approximations to periodic solutions of differential equations. It is generally not necessarily true that that approximate solution which results in the smallest remainder when substituted into the differential equation is the best approximation. For example, consider x = F(x), which has a periodic solution, )(t). Let p(t) be a periodic approximation to (t) such that IP - F(p(t))l < e {for t over an interval equal} to one period of 4(t) f

212 That p(t) for which e is smallest is not necessarily the best approximation. This assumes, of course, that the class of functions from which p(t) can be selected does not contain~((t). For a general statement about so-called c-approximate solutions, see Coddington and Levinson {[101, page 8}. In any event, it can be seen that choosing Q on the basis of minimizing the remainder of the original differential equation when the first approximation has been substituted into it may or may not minimize the error between the actual periodic solution and the first approximation. G. There are still the previously-mentioned unresolved problems associated with the method developed in this study. These are: the a-interval problem, the estimation of the error involved in considering only a finite number of terms in the series for the ai(0)'s; and determination of T for cases in which it is not obvious. s H. It might be worthwhile to allow singular perturbations as well as those considered in this study. Among other possibilities, it should be possible to have "I = 0" and "~ = I" differential equations of different order. As was pointed out in Chapter II, this would make it possible to refine the first approximation at will. I. Note that throughout the development of higher-order approximations, no use is made of any assumed behavior of A( versus z. On the other hand, one of the cornerstones of the usual argument used in the B(jws) development of the describing function is that A(Js ) does not pass higherorder harmonics. Consequently, it should be possible for a special class of nonlinear systems to introduce this concept into the recursion procedure of Chapter IV.

APPENDIX A STABILITY CRITERIA FOR PERIODIC SOLUTIONS OF ORDINARY DIFFERENTIAL EQUATIONS: VARIATIONAL EQUATIONS Asymptotic stability is defined in Chapter I. The purpose of this appendix is to outline the development of the fundamental criteria with which it can be determined whether or not a given periodic solution of a differential equation is asymptotically stable. Consider the following differential equation: F(t, x), (Al) dt = dx dx1 dx where t is the column matrix (dt dts"), x=> (xl,..., X ) and Fdt x f (X...e xn)b and F(t,x) is the column matrix (fl(t,x),..., f(t,x)). Let p(t) be a given solution of Eq. (Al) which is periodic in t with period T. Further, let F(t,x) be analytic in x. (j = 1, 2,..., n) in the region of (t,x)space which contains the solution curve [t, p(t)j, 0 < t < oo, and let F(t,x)be continuous and periodic in t with period T (not necessarily a least period). Let Q be another solution of Eq. (Al), whose initial conditions are within some small neighborhood of those of p(t). Then, letting y = C(- p(t), y =F(t, y + p(t)) - F(t, p(t)). (A2) dt Expanding F(t, y + p(t)) in a power series about x = p(t), Eq. (A2) becomes dy = Fx(t, p(t))y + Z(t,y), (A3) 213

214 6f where F (t, p(t)) is the nxn matrix made up of (j, k = 1, 2, x "xk x-p(t)..., n). Since F(t, x) is analytic in each x. and continuous in t, Z(t, y) = o(lyl) (A4) for small Iyl uniformly in t over any finite t interval. If Z(t, y) is omitted from Eq. (A3), the resulting linear system with periodic coefficients, = Fx(t, p(t))y, (A5) is referred to as the first variation of Eq. (Al) with respect to p(t) or, more succinctly, the variational equation. It is shown by Coddington and Levinson ([10], pages 321-2), among others, that if the trivial solution of Eq. (A5) is asymptotically stable, then p(t) is asymptotically stable. Conversely, if the trivial solution of Eq. (A5) is not stable, then p(t) is not stable. Therefore, the question of the asymptotic stability of p(t) reduces to the question of the asymptotic stability of the trivial solu - tion of the linear system with periodic coefficients in Eq. (A5). Fortunately, the theory for this type of linear system is well developed (see: [10], pages 78-81] and often referred to as Floquet Theory [17]. This theory shows that for every fundamental matrixl, M(t), there exists a periodic nonsingular matrix D(t) with period T, and a constant matrix 1 A fundamental matrix of Eq. (A5) has the following properties: M = F (t, p(t))M (i) det M(t) 0O for -oo < t < oo (ii) In addition, it is convenient to require that M(O) = E, the identity matrix. This allows the solution of Eq. (A5) which satisfies the initial condition I to be written y(t, I) = M(t)I.

215 R such that M(t) = D(t)etR (A6) TR The characteristic roots of the matrix,e, mi(i =,..., n), are referred to as the characteristic multipliers. If the magnitudes of all the characteristic multipliers, mi., are less than one, the trivial solution of Eq. (A5) is asymptotically stable. If one or more of the characteristic multipliers, mi, are greater than one in magnitude, the trivial solution of Eq. (A5) is not stable. The characteristic roots of the matrix R, a.(i = 1,..., n), are referred to as the characteristic ex1 ponents. If all the characteristic exponents have negative real parts, the trivial solution of Eq. (A5) is asymptotically stable. If one or more of the characteristic exponents have positive real parts, the trivial solution of Eq. (A5) is not stable. The relation between the characteristic multipliers and the characteristic exponents is Ta. M.i = 1, 2,...,n (A7) Note that any integer multiple of j —- can be added to a without changT can be added to ai ing the value of m.. It is, in fact, shown in Floquet Theory that the characteristic exponents can be determined only up to integer multiples of j~. This indeterminacy can be appreciated by inspection of Eq. t2g (A6). Clearly, terms of the form e T can be assigned either to D(t) tR tR or e; therefore, the matrices D(t) and e are not unique. On the TR other hand, the matrix e is, in this sense, unique; consequently the characteristic multipliers, mi, are uniquely determined.

APPENDIX B IDENTITIES NECESSARY FOR THE JACOBIAN CONDITION IN CHAPTER III 1. First Canonical Form The Jacobian condition for the first canonical form which appears in Chapter III is as follows: 6a 6a 6a 6a, ~'~q -- - - | ~ 0 7>{A = A(O), i = O}, (3.20) 1 3 2 2 1 where [ from Eq. (3. 5)] T -jWt T +jW t q2(A, 1) = Vln I gl(t, P1,..., pn)dt, (B2) o and %lt A t gl(t, P,..., p) = B(z)F P C1 +. + P (2,28) n n.t + B(z)e(t)- p.A(z)c J j=l ~ Then r'js t jt ~gn iPj ~ v I e E s~ ~] dt) (B3) T +jo gl t n, =a = tV2no c p a2 L=1 aPj a2 ]it (B4) 216

217 T 6ql'- -j t gn ag p ap | 1 vK — s - dt, (B5) a2 0o J= aJ &2 and — _ +ja% t n 6g, 6 q2_ 2n + tn 1 -- 1 dt. (B6) 6 a 2 n 0 bf pj 6az] Note that from Appendix D, kit Xt Therefore, from Eqs. (2.28) and (B7), p (tA ) = aj + =v0 j E g j(s, 2,.. p')d n (j = i, 2..., n),(B9) from which it follows that i | = 1, ~ = O (k f j; = i, 2,..., n) (k = 1, 2,...,n),(B10) Equations (B3), (B4), (B5), and (B6) then become 1 o t BAssumig t (z)[ dX at n) dt. (B7) Equation (2.27 gives Assufrom which ngt follows that 6al ~ ~ ~ ~ ~ ~~I o

218 +jW t -jw \ t c j4t d P s + POE S () = w p) (B12) and assuming that this Fourier series can be differentiated term-by-term through the highest-order derivative contained in B(z), Eq. (Bll) becomes. = Vln T {B(Jcs)So - A(js) }. (B13) Similarly,... = V2n {T { B(-js)SO - A(-jws) } (B14) _a_ = Vl TB(j s)s+2S, (B15) = in - s o and 2 (Ely = v2 TB(-jw )SO2s (B16) 2. Second Canonical Form The Jacobian condition for the second canonical form which appears in Chapter III is as follows: L Sl, -2 O ] #, >{A A(O),= 0 }, (3.30) aal a2 a2 a1 I and, similarly to Eq. (3.5), l('~ ) = Vln f ~ g2(t, Pl'..' Pn)dt, (B7) 0

219 T T +jw t 2 (A,9) = v2n g2(t P..., p)dt (B18) 0 where [from Eq. (2.38)] g (t, p,., p ) B(z)F {p +... + p n + h(t; )} (B19) n At - z p.A(z)e - A(z)h(t; os) j=l p Assuming that jw t, -jcut " a oscyct 2 (B20) J + * E + h(t; )} Z Se ( ) (B20) S=-o0 the necessary identities are, again,: VlnT {B(jws)So - A(jcs)}, (B21)' 7 a2 2n T {B(-j)s - A(- )} (B22 6q2 ^~0~~~~2... = V TB(j )ss, and (B23) where is the coefficient of the term of Eq. (20) with frequency + (B24) 0o

APPENDIX C IDENTITIES INVOLVING DESCRIBING FUNCTIONS 1. First Canonical Form By definition on page 68 F+S T/2 j t -jw t -jw t 0c o s} S (C1) n =. F c + ce e (CO) C10 C10T -T/2 Letting c10 = e (a, 9 real), Eq. (C1) becomes 2 ja s 9 a ) -j(aWt+@) -j (st+) n (a,@) = 2 F t+) -j ) -j( t+ -T/2 Since the above integral has a periodic integrand and an interval of integration equal in length to the period of the integrand, the integral is independent of the location of the interval of integration on the t-axis; correspondingly, it is independent of @. n0 can also be shown to be real for the same reason. Therefore, an 0 -0 bn T/2 j(ot+@) -j(( t+) -j(W t+@) o 2 dFa a s *a a FT J + j dt. (C3) -T/2 r a-+T/ 2 -T/2 Recalling that 7 => A _ A(O) and, = 0, and 220

221 T/2 +2s T/22w t f 1 r dF j d+s 1 dF s j2S it can be seen that Eq. (C3) reduces to on 0o +2s -j29 o o o Similarly, So - -2s +j2@ - (c5) an. Now consider -. Differentiating Eq. (Cl) with respect to "a" gives the following identity: 6n a: + n = S+ -j2@ S+2s (C6) bao 0 0 Similarly, o +- +3 -2s a + ^ n So + +j2@ 2s (C7) 0a 0 0 0 Combining Eqs. (C5) and (C6) gives O a (7A) o 2 iar 2. Second Canonical Form By the definition on page 71 Fo 1 T / jSw t -jw t -jw st n = f / Fc 10 +h(t; C%)}E dt. (C8) 1 -T/2

222 Again, letting c1 =a e (a, 9 real) Eq. (C8) becomes T/2 j(W t=+) (a -J(t+) -j( t+@) (a,,s) = 2__ F{ e + E e + h(t; W) }e dt. (C9 -T/2 Because of the presence of h(t; ws), no(a,Q,wo) is not necessarily independent of 9 and oa or need it be a real number. s Let j(W t+a) a j( t+Q) 0 j Jy t dF {a J( t+h( )} = St (C9A) -— 00 where T = T- Then the following identities can be obtained by differentiating Eq. (C9) with respect to @: = j[so _ S+2s -j2Q - ^S - S - n ](c-o) o 0 0 bn = -j[Is s2s +j2@ _ n (Cll) o o o +2s where S- is the coefficient of the term in the Fourier series, Eq. (C9A), with angular frequency + 20. And, differentiating with respect to "a" s gives the following identities: 6n a S + S+2s -j2 - n, (C12) a a o 0 0 6n* a + s-2s ~+j2 n. (C13) 0 0

223 Combining Eqs. (ClO) and (C12) gives 6n an S a o +n (Cl 4) o 2 I n. (Ci4 0 & - C' 0

APPENDIX D COMMUTATION OF THE OPERATORS - AND B(z) 1. First Canonical Form The purpose of this section is to show that along a solution of the first canonical form, [rj] = [pj(tA,,)], the following equation is valid: CXt nt*X1 B FaF t j^rLB(z)F{rl + *. + rne = B(z)[ (D1) where (j = 1,2,...,n) and z dt. That this is true follows immediately from the nature of the transformation, Eq. (2.19), which relates the (x,x,..,x ))-variables to the rj-variables. First note that the term xt n t B(z)F{r 1 + + rne (D2) implies the following sequence of operations: First, B(z) operates on F(x), B(z)F(x) = G{x,x,...,x(n)}; (D3) then, by the transformation given in Eq. (2.19), {( } Eq. (2.19) G~x~x,...,x^~"^~ G^,r^,..2 r^}, (D4) where G{r,r2,..,rnt} = B(z)F{r 1+... + r }. (D4A) Therefore, 224

225 [ nt} [ 8 { } (D5) [ar B(z)F 1e +... +n = J, r2...trn (D) - J ~J for j = 1,2,..,n. But it follows from Eq. (2.19) that k- (rle+ + r +r ) [ ]i = (krl 1 +... + Xkrne ) (D6) I Li —'n n t "[rj] = r] [ P] for k =,l,...,(n-1). This implies that k lt t -- [(rle +... + ren ) J (D7) -.k X t X t(re 1+... + ren Lk 1 1 n [rj] = [pj for k = 0,1,2,...,(n-1). Since F(x) is an analytic function of x, this, in turn, implies that rd {r %It + nt B()[F{rle + l + rne Ldt n [r. [pj j J' (D8) B(~)F{re 1+.. + r ent } ] = [ Consequently, substituting Eq. (D8) into Eq. (D4A) and noting that the resulting equation is valid for arbitrary [pj], the following relation is obtained. \ t k t G r -rn}= B(F{rle 1 + ** + rne }(D9) It then follows that t xnt a Grl'2 n = F B(- 1 Fr E +... + " } (D10) j j After an interchange of the order of partial differentiation, Eq. (D1O) becomes

226 G rlr2,...,rt} = B(s) -. F rle +. + r n 1) n (Dll) dFJ B(~) j 3 xt - B ) 1 iJ x = re +...+ rne where j = 1,2,...,n. Then, letting [rj] = [pj] and employing essentially the same argument which led to Eq. (D8), the following equation is obtained: XB ~dFI k.t knt t ] {B( F)[dF j kJ n x = rle +... + rne [rj] = [pj] (D12) dr;X.t FdF t ]d -B( — ) 1 (zLx = y(t,I)J dt where j = 1,2,...,n. Therefore, Eq. (D1) is justified. 2. Second Canonical Form In this case it is desired to show that xt x t kn { [(B(z)F{r 1 +... + r n h(t;w)}} 3r + " rn [rj]= [pj] (D13) dF t ] ( d B(z) [x J (z - a x = q(t,Il) for j 1,2,...,n. That Eq. (D13) is valid, follows from an argument which is analogous to the one employed for the first canonical form.

APPENDIX E INVERSION OF THE VANDERMONDE MATRIX The Vandermonde matrix, which occurs in this study, has the following form: I i.. 1 = 1 2..n 2 2 2 ~~~~2 n * 0 0 0 0 0I* 0 (nml) (n-l) (n-i) \ -2 N* where the k.'s (j = 1,2,...,n) are the characteristic roots of the "u = 0" J differential equation. The X.'s (j = 1,2,...,n) are assumed to be distinct; therefore, the above matrix is nonsingular. The inverse of this matrix is, by definition, v21 v22 v2 V - = * 0 0 0........... 0. 0 0 0 Vnl Vn2. vnn The elements of the last column of this inverse matrix (VlnV2nVn) occur repeatedly throughout this study. Expressions for these elements 227

228 are developed in this appendix. From the standard procedure for matrix inversion, [V]n jn det V where [V]n is the (n,j)th cofactor of V; that is nj 1 1... 1 1... 1 \ 2 A,.. \j-i) A(j+l) ** * An IV = ('1)n+Jdet | A2 *' ~ j -1) j+l)''' A (n-2) (n-2) (n-2) (n-2) (n-2) 1 A2 (- 1) j) ( * n Since det V =( 2-A) ( )%32) (43-kl) (4- 43) ( 24- 2) (A-1) )* *-**X (-Xni) ('n (n-2))'' (Xn'1) and, correspondingly, [Vnj = (.1) n+j )(3-l ) (X3-)l)... (A(j.1)-(j. )). (.1)'-) ^nj 3-.-(k)(j+l) -k(i.1 ) )((j+lj2)- (i)*'' ((+l)') n (n'(n-1)).. ( n-n(j+l) (Xn-(jl)). * (n- 1), 1 Kaplan, W., "Ordinary Differential Equations," p. 128, Addison-Wesley Publishing Company, Inc., Reading, Massachusetts, 1958.

229 it follows that n+l (l)n+l (El) Furthermore, it is usually desirable to consider vjn to be a function of f. Recalling that and j = j a (see Sections 4 and 5 of Chapter 2), it can be seen that J (.1)n+l v (). _; (E2) iVn(B) = (-X1 )(2n) ((j therefore, vN(t) t t( *) (E3) Note that the argument of the complex function Vjn(t) is independent of 3. This fact is useful in certain of the stability criteria.

LIST OF REFERENCES 1. J. J. Stoker, Nonlinear Vibrations in Mechanical and. Electrical Systems (New York: Interscience Publishers, Inc., 1950), pp. 147-163. 2. Ibid., p. 8* 3. D. G. Tucker, "Forced Oscillations in Oscillator Circuits, and the Synchronization of Oscillators," J. Am. Inst. Elec, Engrs., XCII (1945), 226-234. 4. I. M. Smirnova, "On Stability of Approximately Determined Periodic Regimes of Automatic Regulator Systems," Transactions of the Second All-Union Congress on the Theory of Automatic Control, I (MoscowLeningrad: Izdat. Akad. Nauk SSSR, 1955), 193-203. 5. R. Adler, "A Study of Locking Phenomena in Oscillators," Proc. I.R.E., XXXVI (1948), 351-357. 6. R. D. Huntoon and A. Weiss, "Synchronization of Oscillators," Proc. I.R.E., XXXV (1947), 1415-1423. 7. E. V. Appleton, "The Automatic Synchronization of Triode Oscillators," Proc. Cambridge Phil. Soc., 1922-23, XXI, 231. 8, B. van der Pol, "The Nonlinear Theory of Electric Oscillations," Proc. I.R.E., XXII (1934), 1051. 9. N. W. McLachlan, Theory and Application of Mathieu Functions (Oxford: Clarendon Press, 1947) 10. E. A. Coddington and N. Levinson, Theory of Ordinary Differential Equations (New York: McGraw-Hill Book Company, Inc., 1955). 11. S. P. Diliberto and G. Hufford, "Perturbation Theorems for Nonlinear Ordinary Differential Equations," Contributions to the Theory of Nonlinear Oscillations III (Princeton: Princeton University Press, 1956), 207-236. 12. M. D. Marcus, "Repeating Solutions for a Degenerate System," Contributions to the Theory of Nonlinear Oscillations, III (Princeton: Princeton University Press, 1956), 261-265. 13. G. Hufford, "Banach Spaces and Perturbation of Ordinary Differential Equations," Contributions to the Theory of Nonlinear Oscillations, III (Princeton: Princeton University Press, 1956), 173-195. 14. K. O. Friedrichs, "Fundamentals of Poincar4's Theory," Proc. of the Symposium on Nonlinear Circuit Analysis April 23, 24, 1953, Vol. II (New York, Polytechnic Institute of Brooklyn). 230

231 LIST OF REFERENCES —Continued 15. D. C. Lewis, "On the Perturbation of a Periodic Solution when the Variational System has Nontrivial Periodic Solutions," J. Rat. Mechanics and Analysis, IV (1955), 795-815. 16. F. R. Moulton, Differential Equations (New York: Dover Publications, 1958). 17. G. Floquet, "Sur les Equations Differentialles Lineaires a Coefficients Periodiques," Ann. 9cole Norm., XI2 (1883), 47-89* 18. L. Marcus, "Continuous Matrices and the Stability of Differential Systems," Math. Zeitschr., LXII (1955), 310-319. 19. W. R. Wasow, "Singular Perturbation Methods for Nonlinear Oscillations," Proc. of the Symposium on Nonlinear Circuit Analysis April 23, 24, 1953, Vol. II (New York: Polytechnic Institute of Brooklyn). 20. L. Bieberbach, Lehrbuch der Funktiontheorie I (Leipzig: B. G. Teubner), 190. 21. J. G. Truxal, Automatic Feedback Control System Synthesis (New York: McGraw-Hill Book Company, Inc,, 1955). 22. L. Cesari, Asymptotic Behavior and Stability Problems in Ordinary Differential Equations (Berlin: Springer-Verlag, 1959).

DISTRIBUTION LIST Copy No. Copy No. 1-2 Commanding Officer, U. S. Army Signal 27 Commander, Air Proving Ground Center, Research and Development Laboratory, ATTN: Adj/Technical Report Branch, Fort Monmouth, New Jersey, ATTN: Senior Eglin Air Force Base, Florida Scientist, Countermeasures Division 28 Commander, Special Weapons Center, Kirt3 Commanding General, U. S. Army Electronic land Air Force Base, Albuquerque, New Proving Ground, Fort Huachuca, Arizona, Mexico ATTN: Director, Electronic Warfare Department 29 Chief, Bureau of Ordnance, Code ReO-l, Department of the Navy, Washington 25, 4 Chief, Research and Development Division, D. C. Office of the Chief Signal Officer, Department of the Army, Washington 25, 30 Chief of Naval Operations, EW Systems D. C., ATTN: SIGEB Branch, OP-347, Department of the Navy, Washington 25, D. C. 5 Chief, Plans and Operations Division, Office of the Chief Signal Officer, 31 Chief, Bureau of Ships, Code 840, DeWashington 25, D. C., ATTN: SIGEW partment of the Navy, Washington 25, D. C. 6 Commanding Officer, Signal Corps Elec- 32 Chief, Bureau of Ships, Code 843, Detronics Research Unit, 9560th USASRU, partment of the Navy, Washington 25, D. C. P. O. Box 205, Mountain View, California 33 Chief, Bureau of Aeronautics, Code EL-8, 7 U. S. Atomic Energy Commission, 1901 Con- Department of the Navy, Washington 25, stitution Avenue, N.W., Washington 25, D. C. D. C., ATTN: Chief Librarian 34 Commander, Naval Ordnance Test Station, 8 Director, Central Intelligence Agency, Inyokern, China Lake, California, ATTN: 2430 E Street, N.W., Washington 25, Test Director-Code 30 D. C., ATTN: OCD 35 Commander, Naval Air Missile Test Center, 9 Signal Corps Liaison Officer, Lincoln Point Mugu, California, ATTN: Code 366 Laboratory, Box 73, Lexington 73, Massachusetts, ATTN: Col. Clinton W. Janes 36 Director, Naval Research Laboratory, Countermeasures Branch, Code 5430, Wash10-19 Commander, Armed Services Technical In- ington 25, D. C. formation Agency, Arlington Hall Station, Arlington 12, Virginia 37 Director, Naval Research Laboratory, Washington 25, D. C., ATTN: Code 2021 20 Commander, Air Research and Development Command, Andrews Air Force Base, Wash- 38 Director, Air University Library, Maxwell ington 25, D. C., ATTN: RDTC Air Force Base, Alabama, ATTN: CR-4987 21 Directorate of Research and Development, 39 Commanding Officer-Director, U. S. Naval USAF, Washington 25, D. C., ATTN: Chief, Electronic Laboratory, San Diego 52, Electronic Division California 22-23 Commander, Wright Air Development Center, 40 Office of the Chief of Ordnance, DepartWright-Patterson Air Force Base, Ohio, ment of the Army, Washington 25, D. C., ATTN: WCOSI-3 ATTN: ORDTU 24 Commander, Wright Air Development Center, 41 Chief, West Coast Office, U. S. Army Wright-Patterson Air Force Base, Ohio, Signal Research and Development LaboraATTN: WCLGL-7 tory, Bldg. 6, 75 S. Grand Avenue, Pasadena 2, California 25 Commander, Air Force Cambridge Research Center, L. G. Hanscom Field, Bedford, 42 Commanding Officer, U. S. Naval Ordnance Massachusetts, ATTN: CROTLR-2 Laboratory, Silver Springs 19, Maryland 26 Commander, Rome Air Development Center, 43-44 Chief, U. S. Army Security Agency, Griffiss Air Force Base, New York, ATTN: Arlington Hall Station, Arlington 12, RCSSLD Virginia, ATTN: GAS-24L 233

OF* MICHIGAN 3 9015 03483 4252 DISTRIBUTION LIST (Cont'd) Copy No. Copy No. 45 President, U. S. Army Defense Board, 61-62 Commanding Officer, U. S. Army Signal Headquarters, Fort Bliss, Texas Missile Support Agency, White Sands Missile Range, New Mexico, ATTN: SIGWS-EW 46 President, U. S. Army Airborne and Elec- and SIGWS-FC tronics Board, Fort Bragg, North Carolina 63 Commanding Officer, U. S. Naval Air 47 U. S. Army Antiaircraft Artillery and Development Center, Johnsville, Guided Missile School, Fort Bliss, Texas, Pennsylvania, ATTN: Naval Air DevelopATTN: E & E Department ment Center Library USAF Security Service, San. 64 Commanding Officer, U. S. Army Signal 48 Commander, USAF Security Service, San. Antonio, Texas, ATTN: CLR Research and Development Laboratory, Fort Monmouth, New Jersey, ATTN: U. S. 49 Chief of Naval Research, Department of Marine Corps Liaison Office, Code AO-4C the Navy, Washington 25, D. C. ATTN: Cthe 91Navy, Washington 25, D. C. A: 65 President U. S. Army Signal Board, Fort Monmouth, New Jersey 50 Commanding Officer, U. S. Army Security Agency, Operations Center, Fort Huachuca, 66-76 Commanding Officer, U. S. Army Signal ReAgency, Operations Center, Fort Huachuca, search and Development Laboratory, Fort Arizona^.~~~~~~~ search and Development Laboratory, Fort ~~~~~Ar~i~zona ~Monmouth, New Jersey 51 President, U. S. Army Security Agency Board, Arlington Hall Station, Arlington, 1 Copy - Deccal Docuens eer ~ie2~~~~~~ V\r'^~~i i Copy - Technical Documents Center 12, Virginia ADT/E 1 Copy - Chief, Ctms Systems Branch, 52 Operations Research Office, Johns Hopkins - C tms stems Branch, University, 6935 Arlington Road, Bethesda Countemeasures Division 14, Maryland, ATTN: U. S. Army Liaison 1 o - i Detection & LoOfficer cation Branch, CounterOfficer measures Division 53 The Johns Hopkins University, Radiation 1 Copy - Chie Jaming & DeLaboratory, 1315 St. Paul Street, Balti- ception Branch, Countermeasures Division more 2, Maryland, ATTN: Librarian 1 Cop ie i o ai 1 Copy - File Unit No. 4, Mail & Records, Countermeasures 54 Stanford Electronics Laboratories, Stan- Records, Countemeasures ford University, Stanford, California, Division ATTN: Applied Electronics Laboratory 1 Cp - Ce Vnerability Br., Electromagnetic EnvironDocument Library ment Division 1 Copy - Reports Distribution Unit, 55 HRB-Singer, Inc., Science Park, State 1 - ots Distribution Unit, Countermeasures Division College, Penna., ATTN: R. A. Evans, te Manager, Technical Information Center Fie 3 Cpys - Chief, Security Division (for retransmittal to BJSM) 56 ITT Laboratories, 500 Washington Avenue, retransmittal to BJSM) Nutley 10, New Jersey, ATTN: Mr. L. A. Nutley 10, New Jersey, ATTN: Mr. L. A. 77 Director, National Security Agency, Ft. DeRosa, Div. R-15 Lab. DeRosa, Div. R-5 L. George G. Meade, Maryland, ATTN: TEC 57 The Rand Corporation, 1700 Main Street, r r Santa Monica, California, ATTN: Dr. J. L. Defense Group, niversity of Michign H~~~~~u-\~^~l~ ~Defense Group, University of Michigan Research Institute, Ann Arbor, Michigan 58 Stanford Electronics Laboratories, Stanford University, Stanford, California, 79-99 Electronic Defense Group Project File, ATTN: Dr. R. C. Cumming University of Michigan Research Institute, Ann Arbor, Michigan 59 Willow Run Laboratories, The University of Michigan, P. O. Box 2008, Ann Arbor, 100 Project File, University of Michigan Michigan, ATTN: Dr. Boyd Research Institute, Ann Arbor, Michigan 60 Stanford Research Institute, Menlo Park, California, ATTN: Dr. Cohn Above distribution is effected by Countermeasures Division, Surveillance Dept., USASRDL, Evans Area, Belmar, New Jersey. For further information contact Mr. I. O. Myers, Senior Scientist, telephone PRospect 5-3000, Ext. 61252. 23)4