TRANSIENT ANALYSIS OF NONLINEAR NON-AUTONOMOUS SECOND ORDER SYSTEMS USING JACOBIAN ELLIPTIC FUNCTIONS by PETER GEORGE DOUGLAS BARKHAM B . S c , Southampton U n i v e r s i t y , 1967 A THESIS SUBMITTED IN PARTIAL FULFILMENT OP THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE i n the Department of E l e c t r i c a l E n g i n eering We accept t h i s t h e s i s as conforming to the ..required standard Research Supervisor.. Members of Committee Head of Departnfeift Members of the Department of E l e c t r i c a l Engineering THE UNIVERSITY OF BRITISH COLUMBIA Ju l y , ' 1 9 6 9 ,--4:y!J'' In presenting this thesis in p a r t i a l fulfilment of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library shall make i t freely available for reference and Study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of ECacfcrCca-l. £.».c, C n.e e r i«.c, J The University of B r i t i s h Columbia Vancouver 8, Canada Date 3o tg, V K X -__7 ABSTRACT A method i s presented for determining approximate solutions to a class of grossly nonlinear, non-autonomous second order differential equations characterized by 2 2 3 d x +m (x + rx ) -t- u.f (x, dx, X ) — 0 , (-1 < r < oo) do? " d r with the restriction that resonance effects be negligible. Solutions are developed in terms of the Jacobian e l l i p t i c functions, and may be related directly to the degree, of non-linearity in the dif f e r e n t i a l equation. An integral error definition, which can be applied to any particular d i f f e r -ential equation, i s used to portray regions of va l i d i t y of the approximate solution i n terms of equation parameters. In practice the approximate solution i s shown to be of greater accuracy than would be expected from the error analysis, and use of the error diagram leads to a pessimistic estimate of solution accuracy. Two autonomous equations are considered to f a c i l i t a t e comparison between the e l l i p t i c function approximation and that obtained from the method of Kryloff and Bogoliuboff. The e l l i p t i c function solution i s shown to be accurate even for heavily damped nonlinear autonomous equations, when the quasi-linear approximation of Kyrloff-Bogoliuboff cannot with v a l i d i t y be applied. Four examples are chosen, from the fields of astrophysics, mechanics, ci r c u i t theory and control systems to ill u s t r a t e , some areas to which the general approximation method relates. CONTENTS Page 1. Introduction 1 1.1 Thesis Outline 5 1.2 Notation 6 2. • The E l l i p t i c Function Approximation 8 2.1 Introduction 8 2.2 Development of the approximation 8 2 .3 E r r o r analysis 18 2 .4 A p p l i c a t i o n of the method 23 2.5 Discussion 28 2.6 Conclusion 29 3. Refinement of the F i r s t Approximation . . . . . . . . . . . 31 3.1 Introduction 31 3.2 The r e f i n e d approximation 31 3*3 Error analysis 37 3- 4 A p p l i c a t i o n ' o f the refinement 43 3.5 Discussion 5 0 3.6 Conclusion 5 0 4. Saturating N o n l i n e a r i t i e s . . ^ 5 2 4.1 Introduction ' • 5 2 4.2 Development of the f i r s t approximation 5 2 4.3 E r r o r analysis 61 4- 4 A p p l i c a t i o n 66' 4.5 Discussion 70 4.6 Conclusion 71 5. Autonomous Nonlinear Systems • 73 5.1 Introduction 73 5.2 Linear d i s s i p a t i o n 74 5- 3 The modified Van der Pol equation 77 5.4 The damped Duffing equation 83 .5-5 Comparison with'the K-B method '86 5-6 Conclusion 88 6. Some Ap p l i c a t i o n s 89 6.1 Introduction 89 6.2 Environmental studies of. mechanical systems 89 6.3 S t e l l a r p u l s a t i o n 90 6.4 Nonlinear frequency modulation 92 6.5 Nonlinear c o n t r o l systems 93 . 6.6 Conclusion 95 7. Conclusion 96 8. References 99 8.1 Van der Pol equation 100 8.2 Duffing equation 101 8.3 C i r c u i t theory 101 8.4 Control systems 102 i i i Appendix 9-1 E l l i p t i c function relationships 103 9.2 Tabulation of E and 5 , 1 0 3 Table 1 . 104 Table 2 109 LIST OF ILLUSTRATIONS Figure . Page 2.0 Negative-resistance o s c i l l a t o r 14 2.1 Error i n t e g r a l as a function of time and (3 f o r ' f (x,x,t) = tx, with p = .2.0. . . • 20 2.2 Error i n t e g r a l as a fu n c t i o n of time and (3 f o r f ( x , i , t ) = tx, with p = 2.0. . . . 20 2.3 Error i n t e g r a l as a function of time, and 6 f o r f( x , x , t ) = -x cos t/m, with p = 2.0. 21 2 .4 E r r o r i n t e g r a l as a fu n c t i o n of time and p f o r • f ( x , i , t ) = tx, with 6 = 0.10 21 2 .5 Error i n t e g r a l as a function of time and p f o r f ( x , i , t ) = tx, with 6 = 0.06. 22 2.6 Error i n t e g r a l as a function' of time and p f o r f ( x , i , t ) = -x cos t/m, with 6 = 0.60. . . 22 2.7 Approximate and exact s o l u t i o n s . o f dfx + 5x + 10x 5 + 1 .118 X x = 0 26 d* 2 2.8 Approximate and exact solutions of d 2x 4 5x + 10x 5+ 0.2 t dx =0. . . . . . . . . . . . . . . 26 d* 2 d r 2.9 Approximate and exact solutions of d 2x -t 5x + 10x 5 + 0 . 5 t d x =0 27 d^ 2 d l r 2.10 Approximate and exact solutions of dfx + 5x + lOx 3 - 5x c o s t = 0. 27 d t 2 3-1 Error i n t e g r a l as a fu n c t i o n of time and 6 f o r f( x , x , t ) = tx, with p = 2.0. 39 3-2 Erro r i n t e g r a l as a fu n c t i o n of time and 6 f o r f(x , x , t ) = tx, with p = 2.0 39 3.3 Error i n t e g r a l as a fu n c t i o n of time and (3 f o r f(x , x , t ) = -x cos t/m, with p = 2.0 40 .3-4 Error i n t e g r a l as a function of time and p f o r f ( x , i , t ) = tx, with 8 =0.10. . 40 v 3-5 E r r o r i n t e g r a l as a f u n c t i o n of time and p f o r f ( x , i , t ) = t x , w i t h 6 = 0.10 41 3.6 E r r o r i n t e g r a l as a f u n c t i o n of time and p f o r f ( x , i , t ) = -x cos t/m, w i t h 8 = 0.60 41 3 . 7 Approximate and exact s o l u t i o n s of d 2x + 5x + 10x 5 + 1.118t x = 0 46 3.8 Approximate and exact s o l u t i o n s of d 2x + 5x t 10x 5 -+ 0.2X dx =0. 46 3 . 9 Approximate and exact s o l u t i o n s of 2 3 d x + 5 x + 10x^ - 5x c o s f = 0 . . 47 d* 2 4 . 1 E r r o r i n t e g r a l as a f u n c t i o n of time and B f o r f ( x , i , t ) = t x , w i t h p = 0 . 5 63 4 - 2 . E r r o r i n t e g r a l as a f u n c t i o n of time and 6 f o r f ( x , i , t ) = t x , w i t h p = 0 . 5 63 4 - 3 E r r o r i n t e g r a l as a f u n c t i o n of time and 6 f o r f ( x , x , t ) = -x cos t/m, w i t h p = 0 . 5 64 4 - 4 E r r o r i n t e g r a l as a f u n c t i o n of time and p f o r f ( x , i , t ) = t x , w i t h 6 = 0.018 64 4»5 E r r o r i n t e g r a l as a f u n c t i o n of time and p f o r f ( x , x , t ) = tic, w i t h B = 0 . 0 2 5 . 65 4 - 6 E r r o r i n t e g r a l as a f u n c t i o n of time and p f o r f ( x , i , t ) = -x cos t/m, w i t h 6 = 0.165 65 4 . 7 Approximate and exact s o l u t i o n s of dfx +- 5x - 2 . 5 x 5 + 0 . 2 t x = 0 . 67 4-8 Approximate and exact s o l u t i o n s of dfx 4- 5x - 2 . 5 x 5 4- 0 . 1 2 5 ? dx = 0 67 4 . 9 Approximate and exact s o l u t i o n s of dfx •»- 5x - 2.5x 5 - 0.825x c o s t = 0 68 dt 2 5-1 Approximate and exact s o l u t i o n s of x 4- x 4- 2 x 3 - 0.1 (1 - x 2 ) x = 0 82 v i 5.2 Phase p o r t r a i t of the s o l u t i o n to x +- x + 2 x 5 - 0.1 (1 - x 2 ) i = 0 82 5.3 Approximate and exact s o l u t i o n s of 'x -+ x + 2x^ + x = 0 85 6.1 L-C c i r c u i t w i t h n o n l i n e a r and time-varying c a p a c i t o r s . . . . 92 6.2 The Lur'e type c o n t r o l system 94 v i i ACKNOWLEDGEMENT I should l i k e to express my gratitude to a l l those i n the Department of E l e c t r i c a l Engineering who, by t h e i r advice, encouragement and i n t e r e s t , stimulated the completion of t h i s work. Sp e c i a l mention i s unquestionably deserved by Mr. A. MacKenzie f o r his . h e l p i n preparing the i l l u s t r a t i o n s , and by Mr. H.H. Black f o r h i s advice on matters photographic. I am e s p e c i a l l y g r a t e f u l to Dr. B.J. K a b r i e l f o r .reading the manuscript at.-such short notice and with such a l a c r i t y , and f o r h i s valuable suggestions. The f i n a l manuscript was typed by Miss D. Mackenzie, and was proof-read by Miss C. S t e i n and Mr. R. Riches. Their assistance i s deeply appreciated. This research would not have been possible without the advice of Mr. K. Morgan, who was the f i r s t to say "Go west, young man", and would equally have been impossible without the guidance of my supervisor, Dr. A.C. Soudack. His encouragement, enthusiasm and i n t e r e s t i n t h i s work has been a constant source of i n s p i r a t i o n . The research presented here was supported by National Research Council of Canadn Grant A-3138, and to i t g r a t e f u l acknowledgement i s given. v i i i 1 1. INTRODUCTION " ... p r a c t i c a l l y a l l of the problems i n mechanics simply are .nonlinear from the outset, and the l i n e a r i z a t i o n s commonly prac t i c e d are an approximating device which i s often simply a confession of defeat i n the face of the challenge presented by the nonlinear problems as such". J . J . Stoker [ 22 ] These words, which appear i n the i n t r o d u c t i o n of Stoker's book "NONLINEAR VIBRATIONS", may w e l l be taken as an i n d i c a t i o n of the complexity with which engineers and s c i e n t i s t s are confronted when endeavouring to describe the operation, of p h y s i c a l systems. In the a n a l y s i s of nonlinear systems the i n v e s t i g a t o r has three basic tools at h i s d i s p o s a l : a) graphical analysis, which, although i t may be applied to grossly nonlinear problems, u s u a l l y implies a lack of p r e c i s i o n . b) d i g i t a l or analog computer simulation techniques which, i n common with graphical a n a l y s i s , can only be applied to a p a r t i c u l a r problem, and c) techniques of t h e o r e t i c a l a n a l y s i s . Of these three approaches only the l a s t enables r e s u l t s to be obtained i n terms of the system parameters. The system may then be designed to conform with s p e c i f i c a t i o n s , and therein l i e s i t s value to the engineer. U n t i l recently, however, l i t t l e research has been conducted i n t o the behaviour of systems which are g r o s s l y nonlinear, and e x i s t i n g a n a l y t i c a l techniques are r e s t r i c t e d to i n v e s t i g a t i o n s of q u a s i - l i n e a r systems. As an example of a q u a s i - l i n e a r d i f f e r e n t i a l equation, B. Van der Pol [23] considered x + w x = j U f ( x , x ) i & dx dt which i s the equation of a simple harmonic o s c i l l a t o r perturbed by a small 2 function f(x,x), the degree of "smallness" being determined by the constant term |t (0 < jx 0) was determined by a method comparable to the K-B method. The importance of t h i s equation i n describing nonlinear o s c i l l a t o r y systems stems from the natural occurrence of odd, i n preference to even, n o n l i n e a r i t i e s . Often a cubic polynomial i n the dependent v a r i a b l e i s s u f f i c i e n t to model such a n o n l i n e a r i t y ; a higher order polynomial may s t i l l , however, be accommodated e i t h e r by"incorporating the high order terms i n the function f(x,x) or by approximating the polynomial by a cubic. An exhaustive analysis of t h i s l a t t e r technique i s to be found i n Soudack [ 2 0 ] . The o s c i l l a t i o n frequency of a nonlinear system.may be strongly amplitude dependent, and t h i s provides f u r t h e r motivation f o r the use of e l l i p t i c functions which have the same property. In 1955, Bogoliuboff and Mitropolsky [ 4 ] presented an extension of the K-B method, which w i l l be c a l l e d the B-M method, and i n p a r t i c u l a r presented a f i r m mathematical foundation f o r the approximation technique. S t a r t i n g from the basic equation 2 x + u> x = ^.f(x,x) , (0 < « 1 ) they assumed the following form f o r the s o l u t i o n x ( t ) : (1) 2 (2) x(t) = a cos ^ + - ^ i u v ; ( a , y ' ) + ^ . u ( a , y / ) + . . . . where 4 ,a = ^ A ( l ) ( a ) -f |U 2A ( 2 )(a) + .... y>=co +JLUBV ; ( a ) + ^ B v ; ( a ) + where the u^"^ are functions periodic i n 2n, and the A ^ and B^^ are functions to be determined. The major advantage of t h i s method i s that successive approximations of i n c r e a s i n g accuracy can be obtained r e c u r s i v e l y . The same basic approach may be adopted f o r non-autonomous systems of the form x + U) x = ji. f ( x , x , t ) where the u^"^ become functions of a, y and t, but the non-autonomous term can only be periodic,, i . e . perturbation terms of the form tx, tx cannot be handled. The B-M method i s e s s e n t i a l l y that of " v a r i a t i o n of parameters" [ 6 ] J but i t i s s t i l l only applicable to q u a s i - l i n e a r systems. Nonlinear, non-autonomous systems i n v o l v i n g no time delay may broadly be divided i n t o two categories: those which e x h i b i t resonance phenomena, and those which do not. Time delay phenomena, although of great i n t e r e s t , w i l l not be considered i n the present work. Resonance e f f e c t s i n quasi-l i n e a r systems lead to p a r t i c u l a r l y i n t e r e s t i n g and unusual phenomena, such as jump resonance [ 6 ] and sub-harmonic resonance [ 9 ], and i t i s possible to analyse these e f f e c t s by the B-M method. The research presented i n t h i s t h e s i s takes, as i t s s t a r t i n g point, perturbation s e r i e s s i m i l a r to those chosen by Bogoliuboff and Mitropolsky, i n an a n a l y s i s of the d i f f e r e n t i a l equation 2 2 ^ d x +ra (x + rx ) 4- /Og (x, dx, X ) = 0 (-1 < r < oo ) d t 2 ' d ^ using the Jacobian e l l i p t i c functions. This i s a non-autonomous equation, but, because of the a n a l y t i c a l d i f f i c u l t i e s associated with resonance 5 e f f e c t s i n grossl y nonlinear systems, the a n a l y s i s i s r e s t r i c t e d to a consideration of non-resonant systems. Equations of t h i s general form a r i s e i n the f i e l d s of, f o r example, astrophysics [ 12 ], c i r c u i t a n a l y s i s [ 35 ] - [ 45 ] and c o n t r o l systems [. 46 ] -[ 5 4 ] . 1.1 Thesis Outline Much of the' work presented i n chapter 2 appears i n a paper [ 2 ] where the first-approximation s o l u t i o n of the equation x + x + px + Bf ( x , i , t) •= 0 (p > 0) i s developed under the assumption that s o l u t i o n frequency remains approximately constant. Three examples, r e l a t i n g to mechanical systems and frequency-modulation c i r c u i t s , are chosen to demonstrate a p p l i c a t i o n of the approxi-mation method. I t i s shown that the approximate s o l u t i o n can be obtained i n a state-equation form where, f o r a p a r t i c u l a r d i f f e r e n t i a l equation, a l l the c o e f f i c i e n t s are known. This enables s o l u t i o n e r r o r to be predicted, and error r e s u l t s are shown f o r the three examples considered. Three s p e c i f i c d i f f e r e n t i a l equations, corresponding to the e a r l i e r examples, are solved using the first-approximation method. The graphical solutions of f i g u r e s 2.7 to 2.10 demonstrate the accuracy which can be maintained up to r e l a t i v e l y large values of the parameter (3, which determines the magnitude of the non-autonomous term. In chapter 3 the same basic equation i s considered, but the assumption that s o l u t i o n frequency remain approximately constant i s removed. The same three examples are considered, new error r e s u l t s are obtained and r e f i n e d solutions of the s p e c i f i c equations.are shown i n f i g u r e s 3*7 to 3 .9 . A p a r t i a l c a n c e l l a t i o n of c e r t a i n f i r s t order terms i s r e l a t e d to three d i f f e r e n t approaches f o r c a l c u l a t i n g the approximate s o l u t i o n , which enhance the accuracy of the r e s u l t i n g s o l u t i o n . Chapter 4 extends the r e f i n e d a n a l y s i s of chapter 3 to consider the equation 3 x + x - px- + B f ( x , i , t ) =0 (0 < p < 1 ) . Comparable examples to those of chapter 3 are chosen to d e r i v e e r r o r r e s u l t s and to f a c i l i t a t e comparison of the s o l u t i o n s shown i n f i g u r e s 4-7 to 4-9 w i t h the corresponding f i g u r e s of chapter 3-Two autonomous d i f f e r e n t i a l equations, a modified Van der P o l equation [ 30 ] and the damped D u f f i n g equation [ 34 ]» a r e considered i n chapter 5- Because n e i t h e r the K-B method or the B-M method co u l d be a p p l i e d to non-autonomous equations, a comparison of the r e s u l t s of the s o l u t i o n method w i t h those obtained from c l a s s i c a l methods was not p o s s i b l e • i n previous chapters. For the autonomous case, however, such a comparison i s p o s s i b l e and i s made i n s e c t i o n 5-5, the r e s u l t s being shown i n f i g u r e s 5.1 and 5-3. I t i s a l s o shown that use of the e l l i p t i c f u n c t i o n s o l u t i o n f o r the h e a v i l y damped D u f f i n g equation r e s u l t s , i n a f i r s t - o r d e r approximation of c o n s i d e r a b l e accuracy which i s , at the same time, i n a simple form. I n chapter 6 examples are chosen from the f i e l d s of mechanics, a s t r o p h y s i c s , c i r c u i t theory and c o n t r o l systems to demonstrate some areas of a p p l i c a t i o n of the general a n a l y s i s presented i n the e a r l i e r chapters. The t h e s i s i s concluded i n chapter 1, where some p o s s i b l e exten-s i o n s of the approximation technique are suggested. There are two s e c t i o n s i n an Appendix, the f i r s t c o n t a i n i n g some p e r t i n e n t r e l a t i o n s h i p s f o r e l l i p t i c f u n c t i o n s , and the second being a t a b u l a t i o n of two constants defined i n chapters 2 and 4-1.2 No t a t i o n I t i s u s u a l to express an e l l i p t i c f u n c t i o n as a f u n c t i o n of two 'variables, the f i r s t being the independent v a r i a b l e u (s a y ) , and the second being a'constant associated with the e l l i p t i c function,„called i t s modulus, k [ 5 ]- When the modulus i s zero the e l l i p t i c sine (Sn) and cosine (Cn) reduce to the c i r c u l a r functions s i n and cos r e s p e c t i v e l y . As k increases to 1 the e l l i p t i c sine and cosine depart from the c i r c u l a r functions and eventually become hyperbolic functions. The g e n e r a l i t y of e l l i p t i c functions can thus be seen to depend on t h i s modulus k. The e l l i p t i c f u n c t i o n s o l u t i o n of the equation x + x + px = 0 (p > 0) 2 1 /2 with x(o) = a and x(o) = 0 i s x(t) = a Cn [ ( 1+pa ) t,k] using the notation Cn[u,k]. However, the modulus k can also be expressed i n terms of ' p' and 'a' as / 2 ( 1 + P a 2 ) 2 Because the value of k i s expressible i n t h i s form, the notation Cn[u,k] w i l l subsequently be replaced by the simpler form Cnu. 2. THE ELLIPTIC FUNCTION APPROXIMATION, 2.1 Introduction Taking the approach of Bogoliuboff-Mitropolsky [ 4 ] as a basis f o r the present a n a l y s i s , approximate solutions to a c l a s s of non-autonomous, g r o s s l y nonlinear second order d i f f e r e n t i a l equations are developed, under the.assumption that s o l u t i o n frequency remains approximately constant and with the r e s t r i c t i o n that resonance e f f e c t s be n e g l i g i b l e or non-existent. The equation analysed i n t h i s chapter i s of the general form 2 2 3 d x + m (x + px ) + /Og(x, dx, f ) = 0 d t 2 / d r where p > 0, and i s a small constant c o e f f i c i e n t which may be e i t h e r p o s i t i v e or negative. This equation represents an o s c i l l a t o r y system with a "hardening" [ 6 ] c h a r a c t e r i s t i c , but i f the c h a r a c t e r i s t i c i s of a higher order than cubic i n the dependent v a r i a b l e i t may s t i l l be approximated by a cubic polynomial [ 14], [ 20]. E x i s t i n g techniques f o r the a n a l y s i s of non-autonomous d i f f e r e n t i a l equations include v a r i a t i o n of parameters [ 17 ] and the WKBJ method [ 6 ], hut are r e s t r i c t e d uO q u a s i - l i n e a r systems. The ana l y s i s presented i n t h i s chapter now makes possible the i n v e s t i g a t i o n of nonlinear non-autonomous, i n ad d i t i o n to g r o s s l y nonlinear autonomous, systems. 2.2 Development of the approximation Using well known techniques [ 3 ], [ 6 ], i t can be shown that any d i f f e r e n t i a l equation of the form it 2 3 x + m (x + px^) +yOg(x,x,t) = 0 where p > 0, yO i s a small constant parameter, and with a general i n i t i a l c o ndition x(0) = X Q and x(o) = 0, can be expressed i n the normalized form x + x H- px 5 + Bf(x,x,t) =0 (1 ) where x(o) = 1,0, x(o) = 0 , p > 0 and 6 i s a small constant c o e f f i c i e n t . The s o l u t i o n of equation ( l ) , when 8 = 0 , i s x(t) = Cn wt,where 1/2 2 Cn i s the e l l i p t i c cosine function [ 5 ], w = (Up)'' and 1c = _ p _ [ 6 ]. 2 . ( U p ) For p 7^,0, consider an approximate s o l u t i o n x ( t ) = a(t) Cn (cot + 9) , where a(t) i s the s o l u t i o n amplitude envelope, and 9 i s a phase-modifying term. For ease of notation i n the subsequent a n a l y s i s , x(t) w i l l be used instead of x ( t ) . In t h i s case 2 . 2 w = 1 + pa and k 2 = - P § _ 2 ( 1 + Pa 2) I f then y = wt + 0 = f(a,t) x = aCny/ = x(a, \j/) D i f f e r e n t i a t i n g equation ( 2 ) : x = a tix. + V 7 ^ JL ba ' and from equation ( 3 ) x = a ^ x + a. d da dt .ha. dt + yls^ x ' ay ( 2 ) (3) (4) The approximation of Bogoliuboff-Mitropolsky [ 4 ] i s based on the following polynomial representation of a and 9: ( 1 ) 2 , ( 2 ) ( i ) , = Ji-Av ; ( a , t ) + c^ A (a, t) + + jx. AK ; ( a , t ) + . 9 . 2 . ( 2 ) , p B V ' ; ( a , t ) +• ^ B ^ ; ( a , t ) + .... -r ^ i B ( i ) (a, t) + ( i ) -„A where |JL i s a small constant parameter, A and B (i=1 , 2 , 3 , • • •) are functions to be determined, and both polynomials are assumed to be conver-gent i n JJL and contain no small d i v i s o r s so that V i : 10 i+1 . (i+1) , i , ( i ) 1 + 1 B ( i + 1 ) « ^ B ( i ) D i f f e r e n t i a t i n g the expressions f o r a. and 9, and discarding terms of -0( |x J: A(1) 2 a t . ,(0 t 9 -^B-' + ^ A ^ B ^ + B ^ a t . where the subscripts denote p a r t i a l d i f f e r e n t i a t i o n with respect to amplitude (a) and time ( t ) . .,^er .If now the fo l l o w i n g assumption i s made: d (cot) = co dt i . e . the frequency of the s o l u t i o n x(t) i s assumed constant over the time i n t e r v a l of i n t e r e s t , then equation ( 4 ) may be written as: x.= a-Cny -'2a (w+O) Sny/ Dn.y/ - a Sny> Dny 9 - aR (co+9) (5) where R^ _S_(SnyDny< ) i . e . R = ( l - 2 k 2 ) Cny + 2k 2 Cr?f and hence R = 1 1 +pa Cn y + pa Cn y Equation (5) can now be expressed as a polynomial function i n JL/. plus one other term (which i s independent of jx ). As p r i m a r i l y only a f i r s t order approximation i s being considered, the polynomial function w i l l be concluded at terms of O(^). Note that, i n general, the expansion can be continued to any order of accuracy by i n c l u d i n g higher order terms i n |U. . The expression f o r x becomes: x = -aR (1+pa ) - j*. SnyDny/ a B ^ + 2A^ 1 ^ ( H p a 2 ) 1 / / 2 ] + p A ^ Cny - (Cr. y +-pa 2Cn 5y ) 2 a B ^ ( 1+pa 2) 1^ 2 ( 7 ) But (from equation (6)) 2 0 *7\ ' / t - ' R (1+pa") = Cr.y -i- pa d r y and from equation ( l ) : = -(aCny + pa^Cn^y/) - B f ( x , i , t ) x From equation ( 7 ) i t now follows d i r e c t l y that 6 f ( x , i , t ) = ^SnyDn^ " I * a B ^ + 2 A ^ ) ( 1 + P a 2 ) 1 / 2 ' (8) A ^ Cny - (Cny +• pa^n^y ) 2 a B ^ ( 1 + P a 2 ) l / 2 A convenient s i m p l i f i c a t i o n of equation (8") can be obtained by taking K J Cn^y K e j Cny/ d y -K (9) i . e . Cn y i s replaced by a term eCny having the same i n t e g r a l on the i n t e r v a l -K — \j/ — K , where K i s the complete e l l i p t i c i n t e g r a l of the f i r s t kind (see Appendix 1 ) . Now &_ (SnyDny/) = ( l - 2 k 2 ) Cny + 2 k 2 Cn 5y and hence 2 k 2 j Cn 5y dy = SnyDny - ( l - 2 k 2 ) j Cny d y By evaluating the i n t e g r a l of CnjC [ 5 ],.the expression f o r e may now be obtained as: 12 (2k -1 )20 4- sin20 4k 20 (10) —1 2 where 0 = s i n k. A t a b u l a t i o n of p, k , e and K i s given i n Appendix 2 , Table 1. I f the i n t e g r a l defined by equation (9) i s evaluated f o r cos y/ , instead of Cn y/, a value of e of 2 / 3 r e s u l t s . By applying l ' H o p i t a l ' s r u l e to equation (10), the same value i s obtained f o r k = 0, when Cn \j/ — c o s y / [5]. A s both jx and 8 are constant small parameters, no g e n e r a l i t y i s l o s t , i f , , i n equation (8), JJ. i s taken equal to 6 [ 4 ], [ 1 6 ] . S u b s t i t u t i n g f o r Cn Y » equation (8) then becomes: f ( x , x , t ) = Sny/Dn^ a B ^ + 2 A ^ 1 ) ( U p a 2 ) 1 / 2 Cn A ( l } - 2 a B ( l ) ( U s p a 2 ) _ (1+pa ) (11) Now from Milne-Thompson [ 16 ] the F o u r i e r s e r i e s expansion of S n y / D n y V may be derived as: Sn\^DnU/ = -d_ ( C n y / ) = -_d_ — . 00 " q ( e +l/2)" 271 Kk / , < (2B+1) 1+qv ' s=0 COS (2s + 1 ) 71 \u 2K 00 or S n v f D n y D y kK s=0 _ q 1+-q (34-1/2)-(2S+-1 ) (2s+1 ) s i n (2s4-1 ) • 71 u/ 2K where q = exp -71K' K , and K' i s defined i n Appendix 1. The expansion of SnyDny can thus be expressed as an i n f i n i t e sum of sine terms; and Cn^ as an i n f i n i t e sum.of cosine terms. This property suggests that C n y and S n y / D n y / may be treated i n the same way as are c o s y / and sin^/ i n the p r i n c i p l e of harmonic balance [ 6 ]. (1 ) ' I f equation (11) could be solved, then the unknown functions A and (1) B could be determined. Two terms prevent the d i r e c t s o l u t i o n of these equations, namely the p a r t i a l (1) (1) d e r i v a t i v e s with respect to time of A and B By d e f i n i t i o n : B ( 1 W 1 > ( a , t ) ' • .: and hence: a t B<'>-„<'>'. a t S u b s t i t u t i n g f o r a , to order u_ i<1>- M<'V1> +A<1> I a t i d ) I f the term i n (D B (D + B ( i ) = p. A ' B + B (D can be neglected i n each case, then A (1) (1) *0) '(1) and B\j/ can be approximated by A and B re s p e c t i v e l y . F i n a l l y equation (11) can be written as: . ' f (x,x , t ) = Sny Dny -Cny a B ^ + 2 A ^ ( U p a 2 ) 1 / 2 A ( 1 ) - 2 a B ( l ) ( l + s p a 2 ) ( U p a V ^ (12) For a p a r t i c u l a r f u nction f ( x , x , t ) , the unknown functions A (1) and B can now be evaluated. To c l a r i f y the a p p l i c a t i o n of equation (12), three examples w i l l now be considered. Complicated equations can be avoided 2 2 by taking a, (1+pa ), (1+epa ) etc. as constant [13 ] when evaluating the (1) (1) functions A and B . I f the amplitude i s slowly time-varying, then the error incurred by t h i s approximation w i l l be slight.. The examples each r e l a t e to systems i n v o l v i n g a time-varying component. The equation of example contains a l i n e a r l y - v a r y i n g phase term, and i s a nonlinear form of Hermite's 14 equation [10 ],. [ 18 ]. A l i n e a r l y time-varying damping term i s considered i n example 2, where the equation might describe o s c i l l a t o r y motion i n a f l u i d of time-varying v i s c o s i t y . Applications of the equations- of both examples 1 and 2 are treated i n greater d e t a i l i n chapter 5. The equation chosen f o r the t h i r d example describes the' steady-state operation of a nonlinear frequency-modulated negative-resistance o s c i l l a t o r , and, because of i t s p r a c t i c a l i n t e r e s t , a b r i e f d e r i v a t i o n i s given below. Cunningham [ 6 ] considers the nonlinear o s c i l l a t o r c i r c u i t shown i n f i g u r e 2.0. . . . R' N.R , F i g . 2.0 Negative - r e s i s t a n c e o s c i l l a t o r . where the negative resistance device has the current-voltage c h a r a c t e r i s t i c 3 i = -ae + be , a and b being p o s i t i v e constants. I f the r e s i s t a n c e R i s not n e g l i g i b l e , the d i f f e r e n t i a l equation d e s c r i b i n g t h i s c i r c u i t i s : e + e _ (1 - aR) + Rbe£ + e L C L C R - a + 3be L C C 0 Now suppose the capacitance C i s varied so that C = C Q (l+m cos wQt) [ 6 ] , where C Q i s the mean value and m i s a small constant. A f t e r some manipulation i t can be shown that the d i f f e r e n t i a l equation of the system takes the form *e +.e (1 - aR) + bRe£ - me L C , L C , 1 - aR + OJ L C 0 0 cos w^t + OJQ _R s i n ui^t L - ef(e,co Qt) + 0(mbR,m ) = 0 where f(e,wt) i s a function i n v o l v i n g constants, e and p e r i o d i c terms i n wQt. This i s a Van der Pol - type equation [ 6 ], and i f terms of order mBR, m may be neglected, i t s steady state s o l u t i o n may be obtained from the equation + e (1-aR) 4- bRe - me LC 0. LC 0 1-aR - co L L L C 0 0 cos co^t 4- U)Q R s i n co^t L = 0 A corresponding, and s l i g h t l y simpler, form was chosen f o r ana l y s i s i n example 3, i . e . 3 • X +• X 4- px - 6x cos co^t = 0 . The same equation i s also considered by Minorsky [17 ] and, i n connection wit a v i b r a t i n g s t r i n g problem, by McLachlan [15 ]• Example 1 Consider the equation X 4- x 4- px 4- 8tx = 0 where f ( x , x , t ) = tx = a t C n y . From equation (12), comparing c o e f f i c i e n t s of Cny and Sn^/Dn^: aB ( lU 2 A ( 1 ) ( H p a 2 ) l / 2 = 0 . A ( 1 ) - 2 a B ( l ) ( l 4 £ P a 2 ) = -at ( U p a 2 ) l / 2 D i f f e r e n t i a t i n g equation (14): a B ^ = ( A ^ 4 - a ) ( U P a 2 ) 1 / 2 2(l+epa ) id) S u b s t i t u t i n g f o r aB i n equation ( 1 3 ) , and rearranging: Xd) 4. 4 A d ) ( 1 4 £ p a 2 ) - a or, expressing t h i s as a Laplace transform: A ( 1 > ( s ) = a 4(l4-epa^ s " - 1 2 AfA 2 ^ s s -4- 4(1+epa ) (13) (14) assuming Ad)( 0-) = A ( 1 ) ( 0 ~ ) = 0. 16 (1 v A K ; ( t ) w i l l consist of a periodic o s c i l l a t i o n plus a p a r t i c u l a r i n t e g r a l term. By considering only the p a r t i c u l a r i n t e g r a l : A ( l > ( t ) 40+epa 2) Note that t h i s process i s equivalent to the "averaging p r i n c i p l e " employed by K r y l o f f and Bogoliuboff [13 ]• - -The expression f o r a. i s then: a. = -B a . 4(l+epa 2) and the s o l u t i o n f o r a(t) becomes: i(t) exp 4(Uep) (15) To determine the p a r t i a l d e r i v a t i v e of A ^ with respect to • (1 ) time i s substituted f o r A i n equation (14) ( t h i s removes one of the approxi-mations made e a r l i e r ) , g i v i n g B v ' = t(1+pa ) ' 2 ( U e p a 2 ) and f i n a l l y 6(t) = B t 2 ( U p a 2 ) l / 2 (16) Example 2 4 ( U e p a 2 ) x 4- x + px 4- Btx = 0 ? 1 / ? In t h i s case f ( x , x , t ) = tx = -at(l+pa ) ' Sn^Dny + O(fit). Note that terms of order Bt i n f( x , x , t ) are equivalent to terms of order B t i n Bf(x,x,t), and may consequently be neglected i n a f i r s t order a n a l y s i s . From equation (12): 2 / 1 ' ) ( l + P a 2 ) 1 / 2 = - a t ( U p a 2 ) 1 / 2 - 2 a B ( l ) ( H e p a 2 ) = 0 t< 2N1/2 (1+pa ) (17) (18) ,(0 D i f f e r e n t i a t i n g equation (18), and s u b s t i t u t i n g f o r aB i n equation (17) y i e l d s A ^ + 4 A ^ (1+epa 2) = -2at (1 + epa 2) Taking the p a r t i c u l a r i n t e g r a l of t h i s equation: and a(t) = exp |^-8t 2/4 A < " - -at , ( 'D , Example 5 S u b s t i t u t i o n f o r A ^ i n equation (18) gives f i n a l l y : 0 ( t ) = - B t ( U P a 2 ) l / 2 4(l+epa 2) x + x + px - Bx cos to^t = 0 ( 1 9 ) (20) where f ( x , i , t ) = -x cos to^t = -aCn^cos t o^t. To preserve the c o n d i t i o n that resonance e f f e c t s be n e g l i g i b l e , i t 2 2 i s necessary to require that co^ .csrl+pa . From equation (12): aB^'K 2 A ^ ( U p a 2 ) 1 / 2 = 0 A ^ - 2aB^1 ^ ( U e p a 2 ) 1 ^ 2 = a cos co.t (1+pa ) (21) (22) By following exactly the same procedure as before: a(t) = exp B(cos w^t - 1 4(l + ep) - to 18 or, i f the exponent i s small: a ( t ) - 1 + B(cos w Qt - 1) 4(1+ep) - w. 0 J (23) F i n a l l y 9 ( t ) = -26 ( l 4 - p a 2 ) 1 / / 2 s i n uiQt ' w 0 [ 4 ( l + e p a 2 ) - w Q 2 J (24) 2.3 E r r o r Analysis The approximation technique has been developed i n i t s e n t i r e t y , but with no dis c u s s i o n of the assumptions made during.the_ a n a l y s i s . As the usefulness of.any approximation depends on i t s range of v a l i d i t y , an i n v e s t i -gation of s o l u t i o n e r r o r to determine these regions of v a l i d i t y i s a necessary adjunct to the approximation method. From equations (3) and (4): x = aCn^ - aySnyDny •- (25) and . >^ 2 3 x = aCny - 2aySnyDny - aySnyDny - a y (Cn y +• pa C n y ) (26) 0+pa 2) Now x =- aCny i . e . , Cny — x a and from equation (25): aCny - x = ay/SnyDny aySnyDny = a x - x or S u b s t i t u t i n g f o r Cny and SnyDny i n equation (26): x = a x + a x - a. x a .y. + 2a y a •2 • / 3\ .y (x+px^) (27) ( U p a 2 ) A l l the c o e f f i c i e n t s i n equation (27) are known ( f o r a p a r t i c u l a r d i f f e r e n t i a l equation), and consequently i n t e g r a t i o n of the two.state equations 19 X 1 ~ X2 • • 2 / 3\ ' (1 + pa 2) y i e l d s the approximate s o l u t i o n x ^ ( t ) . The accuracy of the approximate so l u t i o n s f o r the three examples considered e a r l i e r was determined by comparison' with solutions obtained from a d i g i t a l simulation, using the fo l l o w i n g i n t e g r a l error d e f i n i t i o n : I = e Here x g x i s the s o l u t i o n obtained by numerical i n t e g r a t i o n of the equation •2 x + x 4-px 4- 8 f ( x , i , t ) = 0, which was assumed.for a l l p r a c t i c a l purposes to be the exact s o l u t i o n , and x i s the approximate s o l u t i o n obtained from app r r equation (27). The numerator of equation (28) then represents the area enclosed between the true and approximate s o l u t i o n s . The area under the exact s o l u t i o n was chosen as a normalizing f a c t o r , and I i s then a normalized erro r f u n c t i o n . E r r o r i n an approximation i s often the r e s u l t of phase inaccuracy, and t h i s e r r o r f u n c t i o n i s p a r t i c u l a r l y s e n s i t i v e to such a condition. The e r r o r i n t e g r a l I i s shown as a function of time (t) and 8. f o r the three examples i n f i g u r e s 2.1, 2.2 and 2.3, with the parameter p (which determines the degree of n o n - l i n e a r i t y ) held constant at 2.0. The t o t a l time of i n t e g r a t i o n i s 20 seconds, which corresponds to about four complete periods of the s o l u t i o n (depending on the p a r t i c u l a r f u n c t i o n f(x,£,t)). Figures 2.4, 2.5 and 2.6 show the v a r i a t i o n of the erro r i n t e g r a l with p, holding the parameter 8 constant. I t i s i n t e r e s t i n g to observe that the error f o r the second example i s almost independent of p (as shown i n f i g u r e 2.5). — a_ x.j + x^ — a a x yy 4- 2 a. f a J x - x dt ex app (28) I x dt ex L 0 Fig. 2.1 TIME (SECONDS) -Error integral as a function of time and 8 for f(x,x,t) = tx, with p = 2.0. TIME (SECONDS) Fig. 2.2 Error integral as a function of time and 8 for f(x,x,t) = tx, with p = 2.0. TIME (SECONDS) Fig. 2.3 Error integral as a function of time and B for f(x,x,t) '=- x cos t/m, with p = 2 . 0 Fig. 2.4 Error integral as a function of time and p for f( x , i , t ) = tx, with B = 0 .10 . 22 Fig. 2.6 Error integral as a function of time and p for f(x,x,t) = - x cos t/m, with 6 = 0.60. 23 The approximation derived f o r t h i s case i n d i c a t e s that Jhe amplitude envelope i s independent of p (equation (19))> and also shows the phase-variation term to be small; t h i s p r a c t i c a l r e s u l t i s consequently i n agreement with predicted s o l u t i o n behaviour. When determining the approximate s o l u t i o n to a p a r t i c u l a r d i f f e r e n t i a l equation (from equation (12)), an a u x i l i a r y second order l i n e a r d i f f e r e n t i a l equation must be solved. This equation i s always conservative, and i t s .unforced s o l u t i o n i s discarded. -This i s d i r e c t l y equivalent to the "averaging p r i n c i p l e " . I t might be a n t i c i p a t e d that t h i s approximation would lead to an o s c i l l a t o r y e r r o r i n the f i n a l s o l u t i o n , and indeed such behaviour i s r e a d i l y observed i n f i g u r e s 2 . 1 - 2 . 6 . 2.4 A p p l i c a t i o n of the method Now that the approximation technique has been f u l l y developed, and s o l u t i o n error determined, three equations corresponding to the e a r l i e r examples w i l l be solved to i l l u s t r a t e i t s a p p l i c a t i o n . . Taking s o l u t i o n time as i f t i s the independent v a r i a b l e of the normalized equations, then t = mt where m i s a constant. S u b s t i t u t i n g f o r t, the d i f f e r e n t i a l equations of examples 1, 2 and 3 become: d x + m \x + px ) + d^ 2 m 5 6 t x = 0 (Ex. 1) m 2Bt dx = 0 (Ex. 2) d t - m Bx cos mto^ TJ = 0 (Ex. 3) 3 2 For ease of notation, l e t - m B i n the f i r s t case, and | i = m p f o r the remaining examples. A f t e r s u b s t i t u t i o n f o r B and t i n equations (15] and (16), the approximate s o l u t i o n to the equation 24 d x + m (x + px ) + i*-%x = 0 d^ becomes: x exp 4m (1+e.p) Cn /. 2 n1/2^ (1+pa .) ' X m + X 4m(l+spa ) S i m i l a r l y , from equations (19) and (20), f o r the equation 2 2 3 d x + m (x 4- px ) 4- jxX dx =0 d t x( X ) = exp Cn dr d + p a 2 ) 1 / 2 r . m -4m(l+epa ) and, from equations (23) and ( 2 4 ) , the approximate s o l u t i o n to the equation 2 2 3 d x +• m (x +• px ) - ji x cos m 0JQ"£ = 0 where m =1.0 i s : x(T?) = 1 + jK.(cosV -1) m2|^4(l + e P) - w Q 2J Cn 0 + p a 2 ) l / 2 m -2 jusin m[4(l + £ P a 2 ) - u ) Q 2 ] The most d i r e c t method of obtaining a g r a p h i c a l s o l u t i o n i s f i r s t to plot the s o l u t i o n envelope, and then to determine the time i n s t a n t s at which the approximate s o l u t i o n i s e i t h e r at a l o c a l maximum, minimum, or zero. The f i r s t of these steps presents no problem; the second, however, u s u a l l y d i c t a t e s that a graphical approach be used to determine the phase. B a s i c a l l y , the fo l l o w i n g equation must be solved: y/( V ) = n K where y(x) i s the argument of the e l l i p t i c cosine function, n i s an integer (n=0,1,2 ) u/2 2^-1/2, K • J (1 - k 2 sin20)- l/2d0 , the complete e l l i p t i c i n t e g r a l of the f i r s t kind (see [ 5 ] and Appendix 1), 2 5 2 2 and k == pa ., 2 ( 1+ Pa 2) For the f i r s t example an allowance was made for the decrease i n K with amplitude by solving the following equation graphically: nK ( 5 . 0 ) = fit), where K ( 5-0 ) i s the value of K corresponding to the envelope amplitude at X. = 5 . 0 (the time for which the approximate solution was determined). If 'a' i s this amplitude: k 2 ( 5 . 0 ) pa 2 2(1+pa 2) and K ( 5 . 0 ) may be determined from Appendix 2 , Table 1 by finding the value of 2 K corresponding to this k . A solution of the equation d 2x + 5(x + 2 x 3 ) + 1 . 1 1 8 f x = 0 i s shown in figure 2 . 7 . -The phase correction term for the second example (where f(x,x,t) = tx) 2 1 / 2 i s small compared with m(l+pa ) X , and may be neglected.. The decrease i n K with amplitude must, however, be taken into consideration. The following algorithm was used to determine the approximate solutions of the equation d 2x + 5 ( x + 2 x 5 ) + utdz = 0 , d t 2 d r with ^ = 0 . 2 and 0 . 5 , shown in figures 2 . 8 and 2.9-1. Assume the average value of the amplitude to be 1 . 0 over the f i r s t quarter period, and let 'a' denote this average amplitude. 2 2 . Determine k from the relationship i 2 2 k = pa 2(1+pa 2) '•Mr - - - . 0-4-x(t) • exp- 0.024 t H 1- 1—I *—* 1 ¥• 1 —I h "—* • / \ 1.0 2.0 •) -0.4 -1.0 3.0' I •I I '• V}.0 \ • \ • \ so \ > \ . \ •' - V -Time (seconds) computer solution • • approximate solution . F i g . 2.7 Approximate and exact s o l u t i o n s 2 of d 2x + 5x + 10x 3 + 1 .118Tx = 0 dt 1.04- - — exp -0.05 f 0.4-<(t) n r - i ' — i - •H H 1.0 , . 2.0 N -0.41 -1.0-3.0 :, 4.0 \ s.o Tim e (seconds) F i g . 2.8 Approximate and exact s o l u t i o n s of d 2x + 5x + 10x 3 +- 0. 2 f c l x = 0 ~2 d f d t r computer solution. . approximate solution to 0-4-x(t) N -exp -0.125t2 v . •4 f— ;.o -+c 1 1 2.0 3.0 4.0. - 5.0 -0.4 -1.0 Time (seconds) computer solution . . . approximate solution. F i g . 2.9 Approximate and exact solutions of d 2x + 5x +1 Ox 3 +• 0.5tjLx=0 d*" dt *(t) -0.4 •1.0 • 1-0.1112 (1-cos t) 7 •I '•I 7 J J V. t — i — 2.0 i: i: 4.0'-t t \-V V ' v. fime (seconds] 6-0 V. 0.0 —M ;o-.o .v, computer solution . approximate solution F i g . 2.10 Approximate and exact solutions of d x + 5x + 10x^ - 5x cos V = 0 dV1 3. 4. Determine the value of K corresponding to k from Appendix 2, Table 1. Calculate the quarter period T from the r e l a t i o n s h i p 2V1/2 T = K to where to = m(l+pa ) Determine the average amplitude over the next quarter period (assuming i t to be of duration T seconds) from the s o l u t i o n envelope ( i . e . i n t h i s case exp - ^ 4 . J ), and return to step 2. For the f i n a l example: f\l) = : ( U P a 2 ) l / 2 m tr - 2 fx s i n i n X m^4(l + epa 2) to. = nK where K i s c a l c u l a t e d f o r the average value of the s o l u t i o n envelope (from Appendix 2, Table 1), i . e . an amplitude a = 1 m [i(Uep) - co 0 2 J The maxima, minima and zeroes of the approximate s o l u t i o n are then given by X - nK 2 jj-sint m (1 + pa 2) l / / 2 n^ H+O + spa^) - (oQ2 J where 'a' denotes the. average amplitude. This equation i s r e a d i l y solved g r a p h i c a l l y . An approximate s o l u t i o n of the equation 2 3 d x + 5(x + 2x ) - 5 x c o s t = 0 d ^ i s shown i n f i g u r e 2.10. 2.3 Discussion The major innovation of t h i s a n a l y s i s i s - a consequence of the (1 ) (1 ) def i n i n g equation f o r the unknown functions A and B (see equation (12)), which leads (for a s p e c i f i c example) to a l i n e a r , undamped second order 29 differential equation for A (t). The complementary solution of this equation i s discarded, leaving a particular integral term; this approach (as mentioned in section 2 .2) i s equivalent to the "averaging method" of Kryloff and Bogoliuboff. The K-B method and also that of Bogoliuboff-Mitropolsky, however, both reach this averaging procedure after a Fourier series expan-sion which assumes periodicity ..in the natural, i.e. unmodified, frequency of the unforced system differential equation. For systems where the solution frequency i s rapidly varying this approach might be expected to result in inaccuracy. Note that this problem does not arise i n the present analysis, as no such assumption about periodicity of the equations defining the amplitude and phase variation terms i s necessary to the derivation of the expressions A^ ^ and The anticipated oscillatory behaviour of solution error is shown i n figures 2.1 to 2 . 6 , but the graphical solutions (see figures 2.7 to 2.10) indicate that greater accuracy i s obtained in practice than would.be expected from the error analysis. This apparent increase in accuracy i s a consequence of the expedients employed,when deriving the graphical solutions, which could not be incorporated i n the d i g i t a l computer simulation. The main value of the integral error results i s the portrayal of qualitative behaviour, but i t i s not disadvantageous to have a pessimistic estimate of solution error. Actual error in the graphical solutions i s much less than that predicted, and hence the error results provide a conservative aid in determining regions of vali d i t y of an approximate solution. 2 .6 Conclusion A method for determining approximate solutions to a class of non-linear, non-autonomous differential equations characterized by 2 ^ x + m (x H- px ) + jxf(x,x,t) = 0 30 has been developed using the Jacobian e l l i p t i c functions; the approximation i s easy to apply, and, although only a\ f i r s t order approximation i s employed, the three examples considered demonstrate the accuracy which can be obtained. An expression f o r s o l u t i o n e r r o r makes i t possible to determine the accuracy of the approximate s o l u t i o n f o r any function f ( x , x , t ) once the amplitude envelope and phase r e l a t i o n s h i p s have been derived. This approximation method i s r e f i n e d .in the following chapter, to take i n t o consideration the v a r i a t i o n of s o l u t i o n frequency with amplitude. 31 .,3. REFINEMENT OF THE FIRST APPROXIMATION \ 3.1 Introduction In the previous chapter first-approximation solutions were developed under the assumption that the frequency of the non-linear oscillation remained constant (at least to f i r s t order). Although the accuracy of the solutions obtained demonstrated that this assumption was valid over a short time interval, more accurate solutions over longer time intervals might be anticipated i f allowance were, made for the variation of frequency with amplitude. In the present chapter, the analysis of chapter 2 i s extended to account for this variation of frequency. New error-integral results are shown, and refined solutions of the three examples considered i n chapter 2 demonstrate that an appreciable improvement in accuracy can be obtained. 3.2 The refined approximation The solution of the normalized equation •, ••- x + x -f px 3 + 8 f(x,i,t) = 0 , with x(o) = 1.0, x(o) = 0 and where p > 0, i s again taken as x(t) = a ( t ) Cn (cot + 9) = a Cn\j/ where . 2 to = 1 + pa y = cot + 9. However, in the present case, the assumption d [cot] = co i s no longer made. dt As before, writing '~ ' x = x(a, y) and differentiating with respect to time, the following equations are obtained: Also k = JJLA^1 ^ (a,t) + j/A^2\a,t) + 9 =|AB^(a,t) + jU 2B^ 2\a,t) 4-a ^ A ( ^ + 0(j,2) where the subscript denotes partial differentiation with respect to time. Now where and .• = CO + tco + 9 ( A J 2 ^ / 2 co = (1 + pa ) ' co ' = paa (1+pa ) Retaining only f i r s t order terms in j^, the expression for y becomes y/ = (1+pa ) + JU. B^-+ p a t A ^ " (29) Similarly y 2 = . U p a 2 +- JJL 2 B ( l ) ( U p a 2 ) l / 2 + 2 p a t A ( l ) ] (30) 2 2 2 2 where terms of 0 ( ^ L ) , 0(JU. t) and 0(|A- t ) are neglected. It should be notei1 that this assumption i s l i k e l y to cause a deterioration i n solution accuracy 2 2 for large t, i.e. when |^-t i s no longer negligible. For y/ : Y = 2co + t d _ ' d t paa /, 2x1/2 _(1+pa ) + 0 or >y = 9 + _2paa, ( H p a 2 ) 1 / 2 + Pt aa +• a ( 1 + P a 2 ) ' / 2 ( W p a 2 ) 3 / 2 _ .2 2 The term a. i s of order |U , and hence to f i r s t order in jU. d ) 2paA 0')' patA (1) ^ 2^/2 i. 2s 1/2 (1+pa ) (1+pa ) (31) 33 The r e f i n e d r e l a t i o n s h i p corresponding to equation ( 5 ) i s therefore 'x = a C n ^ - 2a^SnyDny - aSn^Trnfy - aR y (32) By appropriate s u b s t i t u t i o n i n t o equation ( 32 ) , the f i r s t - o r d e r approximation to x becomes: x =-(aCny + pa 3Cn 5y ) • jiSnyj/Dnxj/ A ^ Chy- - (Cn^ 4- pa2Cn3y/ ) 2aBv ' 2.(1) 2..(1) 2pa A v y4- pa tA v _ f« 2N1/2 . (1+pa ) . (1 )/ , 2x1/2 , 0 2 , . ( 1 ) v ' (1 +pa J ' 4- 2pa tA v ' pa ( H p a 2 ) 2 "5 ^ I f the term (Cn\j/ +- pa Cn ^ ) i s approximated by Cny/(l4epa ) (where the quantity e i s defined by equation ( 10 ) ) , then a r e f i n e d r e l a t i o n s h i p corres-ponding to equation (11) may f i n a l l y be w r i t t e n as: " (33) f (x,x,t) = SnyDny a B ^ 4 - 2 A ( 1 ) ( H P a 2 ) 1 / 2 + "n 2.(1-) 2.,(1) 2pa A v ' 4- pa t A ^ / ^ 2\l/2 ... (1 +pa ) ; Cn Y A ^ - ( U e p a 2 ) (Upa 2 ) . .L 2 a B ^ ( l + P a 2 ) l / 2 + 2 P a 2 t A ( 1 ) ' A comparison of t h i s equation and equation (11) reveals that the major modification i s to the term i n v o l v i n g SnyDny, but terms i n v o l v i n g t e x p l i c i t l y are now incorporated i n the c o e f f i c i e n t s of both CnyV and Sn^ Dny> To i n v e s t i g a t e the e f f e c t of t h i s refinement on equation s o l u t i o n s , examples 1, 2, and 3 of chapter 2 w i l l b r i e f l y be reconsidered. Example 1 Consider the equation 3 x 4- x 4- px 4- Btx = 0 with x (0) = 1.0 and i(o) = 0. For t h i s equation f ( x , i , t ) = tx = a t C n y . Applying equation (33 ) , and comparing c o e f f i c i e n t s of C n ^ and Sny/Dn^: 34 .B<;>+ 2 A < 1 > ( l t p a 2 ) , / 2 4 2pa 2A^ l )4 P a 2 t A ^ L d 4 p a 2 ) 1 / 2 ( 3 4 ) A ( J . } - (Uepa 2) ' (Hpa 2 ) 2 a B ^ ( l 4 P a 2 ) 1 / 2 4 - 2pa 2tA^^ = -at ( 3 5 ) Taking the partial derivative of equation ( 3 5 ) with respect to time: aB ( l t } = ( A { \ \ 4 a). ( 1 + P a 2 ) l / 2 -_ Pa 2(l4epa /, 2 x 1 / 2 (1+pa ) ( A ^ 1 ) + t A ^ ) (1 ) Substitution for &£> ^ ' i n equation ( 3 4 ) then yields A ( J ^ 4 - 4 A ( l ) ( l 4 e p a 2 ) 1 + pa 2(1+pa) -a . ( 3 6 ) At this juncture in chapter 2, the equation corresponding to equation ( 3 6 ) had been derived as A ( 1 )4- 4A ( 1 ) (Wepa 2) = - a , and this expression was then written as a Laplace transform. Inversion of the transform gave A ^ \ t ) as the superposition of a complementary solution and particular integral, of which only the particular integral was retained. The particular integral of equation ( 3 6 ) may be obtained by setting A^|| 0, and hence A 0 ) ( a , t ) = 4(Hepa 2) 1 + pa 2.(1+pa2) 2 2 Also, the term psi i s equal to k , where k i s the .modulus of the e l l i p t i c 2(1+ Pa 2) '••a--. (1) function. Integrating the expression for A , and assuming that the amplitude 'a' [ 13 ] may be treated as a constant without incurring gross 35 inaccuracy: i(t) = exp 4(l+£ Pa 2)(l+k 2) _ ,(1). S e t t i n g A .^ •,=• 0 i n equation (35) y i e l d s , a f t e r some manipulation: B ( l ) ( a , t ) = t ( U P a 2 ) l / 2 + 2. pa t 2(1+ Epa 2) 4 ( l + £ P a 2 ) ( H k 2 ) ( H p a 2 ) l / 2 and f i n a l l y : G(t) = p t 2 ( U p a 2 ) l / 2 + „ 2^2 Bpa t 4 ( H e p a 2 ) 8 ( l 4 - e p a 2 ) ( l + k 2 ) ( l + p a 2 ) l / 2 (37) (38) I t i s p a r t i c u l a r l y i n t e r e s t i n g to observe that the refinement 2 introduces terms i n v o l v i n g k , a parameter- which increases with the degree of n o n - l i n e a r i t y and which, f o r t h i s example, decreases the rate of decay of the s o l u t i o n envelope (by comparison with the. unrefined approximation). From fi g u r e s 2.7 to 2.10 i t might be a n t i c i p a t e d that a decrease i n the rate of decay of the amplitude envelope would improve s o l u t i o n accuracy s u b s t a n t i a l l y . This, i n se c t i o n 3-4, i s seen to be the case and the r e f i n e d approximation predicts s o l u t i o n behaviour accurately. Example 2. 3 x + x + px + Btx = 0 with x(0) = 1.0, x(0) = 0, and f o r which f(x,x,t) = t i = -at ( l + p a 2 ) 1 / / 2 Sn^Dny + 0(st) . 2 Note that terms of O(fit) i n f ( x , i , t ) are equivalent to terms of order 6 t i n Bf(x,x,t), and may consequently be neglected. From equation (33) : a B ^ + 2 A < 1 > ( U P a 2 ) 1 / 2 + 0 2.(1) 2,.(1) 2pa AK y+ pa_fcA_\j._ ( 1 + P a 2 ) l / 2 - a t ( l + p a 2 ) , / 2 (39) 36 (1+pa 2) _ W 1 ) ( l + P a 2 ) l / 2 + 2 p a 2 t A ^ ) ' (40) Taking the p a r t i a l d e r i v a t i v e with respect to time of equation (40) and s u b s t i -t u t i n g f o r aB^^ i n equation (39) gives: + 4A^ 1 ^ (l + E p a 2 ) ( l + k 2) = - 2 a t ( l + e p a 2 ) The p a r t i c u l a r i n t e g r a l of t h i s equation i s A<1> = -at , 2 ( l + k 2 ) and by i n t e g r a t i o n : t(t) exp 1.4(1+^) (41) (1) S u b s t i t u t i o n f o r A i n equation (40) y i e l d s : B ( l ) ( a , t ) = _ _ _ _ _ _ d + p a 2 ) l / 2 2 ( l + k 2 ) ( U p a 2 ) l / 2 4 ( H k 2 ) ( l + e p a 2 ) 2,2 pa t In t h i s example, i t i s no longer v a l i d to assume that the amplitude remains approximately constant, and the phase of the s o l u t i o n i s more e a s i l y deter-mined i n an approximate manner as described i n section .3.4. Example 3-x 4- x + px - Bx cos LO t == 0 with x(o) = 1 . 0 , x(o) = 0, and where f ( x , i , t ) = -x cos co^t = -aCnycos co^t. From equation (33): a B ^ ^ W 2 ) 1 7 2 - * - ' 2. (1) 2 . , (1 ) " 2pa A -+- pa t A ^ ' (A 2x l /2 (1+pa ) = 0 A ^ - (Uepa" (1+pa 2) 2 a B ( l ) ( l + P a 2 ) 1 / 2 4 - 2 P a 2 t A ( 1 a cos co^t (42) and hence: 37 A ( l ) ( a , t ) = -au) 0sinw 0t 4.(l+Epa 2)(l + k 2) - co Q 2 Taking- the p a r t i a l d e r i v a t i v e of with respect to time, and s u b s t i t u t i n g i n equation (42): B ( l ) ( a , t ) = L [ 4 ( l 4 E p a 2 ) ( l 4 - k 2 ) -VJL ( 1 + P a 2 ) l / 2 The expression f o r a(t) may f i n a l l y be obtained as pa tiu,,sinu)_t 2\1/2/.,.2\ , £ 0 0_ - 2(1+pa J ' (,1+k )cos coQt i(t ) = 1 -f 6(cos coQt -1 ) . 4 ( U e pa 2)(Wk 2) - LO q2 (43) The phase-modifying term 6(t) w i l l be discussed l a t e r ( i n section 3-4), but at t h i s juncture i t i s s u f f i c i e n t to observe that i n t e g r a t i o n of the equation (1 ) d e f i n i n g B (a,t) i s unnecessary. A comparison of these r e s u l t s with those obtained i n chapter 2 shows that the main consequence of the refinement i s the i n t r o d u c t i o n of f a c t o r s i n v o l v i n g k , the modulus of the e l l i p t i c function. As k increases with the parameter p (the c o e f f i c i e n t governing the degree of n o n - l i n e a r i t y ) , the e f f e c t of tha refinement should be most noticeable f o r high values of p. This hypothesis i s inv e s t i g a t e d f u r t h e r i n the f o l l o w i n g section. 3-3 Er r o r Analysis The equations d e f i n i n g the approximate s o l u t i o n i n algebraic terms, derived i n s e c t i o n 2.3, are s t i l l v a l i d , but the q u a n t i t i e s a, a, a, y/ and \j/ are now d i f f e r e n t . The major source of-'error i n any approximation i s u s u a l l y to be found i n the phase-modifying term, or, more generally, i n the argument of the periodic function. In the un-refined case considered i n chapter 2 t h i s argument (y/) was taken as y = cot + 9 and co -t 9 38 or, s u b s t i t u t i n g f o r to and 9: f = ( H p a 2 ) l / 2 + / * B ( l > In the present (refined) case, from equation (29): j/ = (1+pa 2) 1 / 2 + jl B ^ + . p a t A ^ ( H ' P T ) 1 7 " 2 -i . e . , the refinement has introduced a term U patA^ 1^ ( I 4 p a 2 772 i n t o the expression f o r \j/ . To demonstrate the e f f e c t of' t h i s modification, consider the values of y f o r example 1 (where f ( x , x , t ) — tx).. From chapter 2, the un-refined expression f o r y i s : • f = d+pa 2) l / 2+ B t ( l + p a 2 ) l / 2 • 2(l+epa 2) S u b s t i t u t i o n f o r A ^ ^ and B^^ into equation (29) gives the r e f i n e d expres-s i o n f o r \1/ as: f f - ( H p a 2 ) l / 2 + 8 " t ( U p a 2 ) l / 2 + 2(l+-epa2) 4 ( l + s p a 2 ) ( l + k 2 ) ( H p a 2 ) l / 2 p a 2 t 2, pa t 4 ( H s P a 2 ) ( l + k 2 ) ( l + P a 2 ) l / 2 which y i e l d s the i n t e r e s t i n g r e s u l t that, to a f i r s t - o r d e r approximation, the change i n \I/ due to _2_ (to) i s cancelled by the change i n y/ due to c> (6). ' .^a 5a • A s i m i l a r r e s u l t can r e a d i l y be obtained f o r the remaining two examples. Any improvement i n s o l u t i o n accuracy w i l l , as a consequence, only a r i s e from an.improved approximation to the s o l u t i o n amplitude and should be more apparent f o r high values of the parameter p. Families of error curves f o r the r e f i n e d case, corresponding to those obtained e a r l i e r , are shown i n f i g u r e s 3-1 to 3-6. Figures 3-1, 3-2 and 39 F i g . 3.2 Error integral as a function of time and 6 for f ( x , i , t ) = tic, with p = 2.0. 40 TIME (SECONDS) Fig. 3-4 Error integral as a function of time and p for f(x,x,t) = tx, with (3 = 0.10. 41 TIME (SECONDS) _ Fig. 3.5 Error integral as a function of time and p for f(x,x,t) = tx, with 6= 0.06. 42 3-3 refer to the variation of the error integral with 8 for the three examples, with p held constant at 2 . 0 ; figures 3-4, 3-5 and 3-6 show error integral variation with the parameter p, holding 6 constant. A comparison of figures 2.1 and 3-1 and 2.4 and 3.4 reveals that less error i s incurred by assuming to to remain constant for the f i r s t example, where f(x,x,t) = tx . The deterioration i n accuracy when using the refine-ment (which i s small, but not insignificant) i s the consequence of neglecting 2 2 terms of order u. t when deriving the approximation (see section 3-2). This particular case, i t should be noted, contains no energy-dissipative term (e.g.. x), so that any change in amplitude can result only from a change i n solution frequency. The time varying function tx alters the solution frequency or, more specifically, the total phase of the solution, with the result that substantial phase variation i s the dominant feature of the solution. This factor, coupled with the slow change in amplitude envelope with time, indicates that, of the three examples considered, the f i r s t w i l l be the most sensitive to phase error. The remaining two cases show that considerable improvement in solution accuracy i s obtained by using the refined approximation (see figures 2 . 5 , 2 .6 , 3-5 and 3-6) at high values of the parameter p, as had been anticipated. It should be emphasized that error predicted from the error-integral results tends to be pessimistic when compared with the f i n a l graphical solution,. as many expedients may then be employed which cannot be incorporated i n the d i g i t a l computer simulation. This becomes particularly apparent in the next section, where approximate solutions corresponding to the three examples of section 3-2 are obtained, and i t becomes evident that an improve-ment in solution accuracy results from use of the refined approximation for each example, including the f i r s t . 43 3.4 Applicaticn of the refinement For a non-autonomous function f(x,x,t) equation (29) w i l l always he of the form j/ = ( 1+pa 2)^ 2 + Bg(t), where g(t) i s a function containing t expl i c i t l y , i.e., the cancellation of terms of order 6 i s not complete. It w i l l be shown later (in chapter 5) that for a linear autonomous function f(x) this cancellation of f i r s t order terms i s complete. This i s an important result, but i t w i l l be treated i n detail i n section 5-2. In the non-autonomous case, however, the form of the expression for y/ can be used to indicate the method of approach when determining the graphical solution of a particular equation. By comparison with the phase-modifying term 9(t), the amplitude envelope a(t) i s , for a specific example, readily determined and hence the amplitude behaviour of the solution can be predicted. Also the function B^^ can readily be found once A^^ i s known, and the expression \j/ = (1+pa 2) 1/ 2 + B B ( 1 ) + _ P a t A ^ ' ( H p a 2 ) l / 2 (29) " 2 1 /2 can be evaluated tc datermine the function g(t) i n y/ = (1+pa ) + Bg(t), where g ( t ) = B ( 1 ) + ' patA ( l ) ( A 2X1/2 (1+pa ) The choice of approach should then depend on the nature of a(t) (the amplitude envelope), and the quantity g(t). a) If, over the time interval of interest, the amplitude envelope i s approxi-mately constant, then the phase of the approximate solution can be obtained (1) from 6(t). When integrating the expression for B to obtain 9(t) the amplitude may then be assumed constant, and the parameter K (the complete e l l i p t i c integral of the f i r s t kind) i n the expression f i t ) = nK , which i s us.ed to determine l o c a l maxima, minima and zeroes of the s o l u t i o n , takes a constant value corresponding to the i n i t i a l s o l u t i o n amplitude a(o). The evaluation of the quantity y ( t ) from 2 1 /2 ^ ( t ) = (1+pa ) t + 9(t) (where 9(t) i s given by, f o r example, equation (38)), should, however, take i n t o account v a r i a t i o n i n the amplitude. This p r o v i s i o n i s necessary to ensure that equation (29) i s not v i o l a t e d , i . e . the p a r t i a l c a n c e l l a t i o n of v a r i a t i o n s i n j/ caused by changes i n a(t) and e(t) must be taken i n t o consideration,- although i n a c i r c u i t o u s manner. Note that i t i s much l e s s tedious to integrate the expression (1) f o r B assuming the amplitude to be constant than i t i s to integrate equation (29) where t h i s assumption cannot be made. The s o l u t i o n of the normalized equation (e.g. equation 1) i s taken i n the form: x(t) = a(t) Cn [ ( H p a 2 ( t ) ) l / 2 t + e(t)] and the l o c a l maxima, minima and zeroes of the s o l u t i o n are determined from y(t) = (l + p a 2 ( t ) ) 1 / / 2 t + 0(t) = nK(0) . This approach i s adopted f o r the f i r s t example of t h i s s e c t i o n . I f the mean value of the amplitude envelope i s constant, and hence g(t) i s p e r i o d i c , then equation (29) may be integrated assuming the amplitude 'a' to take i t s mean value. The value of K i n y ( t ) = nK should corres-pond to t h i s mean-value of the. amplitude envelope. " '.'This approach makes i t unnecessary to evaluate ©(10'%, as y ( t ) may be obtained d i r e c t l y from equation (29). The approximate'solution i s then of the form x(t) = a(t) Cny . This method i s applicable to the f i n a l example considered. I f the amplitude v a r i e s s u b s t a n t i a l l y over the time i n t e r v a l of i n t e r e s t neither of these approaches can be adopted, and values f o r the l o c a l maxima, minima and zeroes of the approximation s o l u t i o n must be found by considering v a r i a t i o n i n the parameter K (see the algorithm of 4 5 • 2 1 / 2 section 2 . 4 ) . Note that the existence of g(t) i n \jJ = (1+pa ) -1- 6 g(t) for the non-autonomous case implies that the argument f of the e l l i p t i c function i s both time and amplitude dependent, which provides additional j u s t i f i c a t i o n for considering the variation of K with amplitude (again see the algorithm of section 2 . 4 ) . In.this case i t i s unnecessary to evaluate either ^ or 8(t) although, depending on the particular system, some inaccuracy may be incurred. The solution i s obtained i n the form x(t) = a(t) Cn [ ( U p a 2 ( t ) ) l / 2 ] . Exactly the same equations are those chosen i n section 2 . 4 w i l l now be used to demonstrate practical application of the refined approximation. For the f i r s t example, setting t = mX and jU=m 6 : 2 2 3 d x + m (x -+ px ) + uXx =0 2 At and from equations ( 3 6 ) and ( 3 7 ) : x( t ) = exp _4m2 (l+-ep)(l-t-k2) Cn (1+pa ) ' X m + 4m (1 +epa ^ P a 2 ^ 2 8m(H Epa 2)(l+k 2)(l + p a 2 ) l / 2 _ This solution i s shown i n figure 3 . 7 with m = 5 , p = 2 and p. = 1.118. The local maxima, minima and zeroes of the approximate solution were deter-mined by solving graphically the following equation: f ( X . ) = nK(0) (n=0,1 ,2 , . . . ) where K(o) i s the value of K (from Appendix 2 , Table 1) corresponding to the i n i t i a l amplitude. Note that this approach differs from that i n section 2 . 4 , because the change in frequency of the solution i s taken into account by the refined analysis; additional compensation for this frequency 46 t.o 0.4 aft) exp - 0.0183t i v ». V V V * X X i. H -1.0 /•0 2.0 G.O t 4-0 5.0 x V. Time (seconds) computer solution approximate solution — Pig. 3.7 Approximate and exact solutions of dfx + 5x + 10x 3 + 1 ,118t x = 0 d t 2 1.0 0-4 x(t) exp-0-0375 t v -+-1.0 20 -0-4 1.0] —• 3.0 * f f r \ f 4.0 SO Time (seconds) computer solution • • approximate solution Pig. 3-8 Approximate and exact solutions of d 2x +5x + 10x 3 + 0.2t & =0 dT 47 cos t ) 0.4 x(t) -0.4 / TP — i 6-0 K r t r to.o -/.oh-Time (seconds) computer solution approxim at e sotu tion F i g . 3.9 Approximate and exact s o l u t i o n s 2 3 of d x + 5x + 10x - 5x c o s t = 0 dtr2 48 change i s therefore unnecessary. For the second example, where ^ = m 6 and 2 2 3 d x '+ m (x +• px ) + uli_dx. = 0 d t 2 d r the approximate solution i s : x( X ) = exp .2 Cn .4(Hk 2). /, 2x1/2 ' m(1+pa ) ' V In this case phase cannot be determined e x p l i c i t l y in a simple form, and approach (c) i s necessary. In the analysis of the unrefined approximation the term 9(t) was found to be negligible, and the variation of phase with amplitude was deter-mined from the variation of K, the complete e l l i p t i c integral of the f i r s t kind. A suitable algorithm i s given in section 2.4. A solution of the equation , 2x + 5(x + 2r 5' d ^ dfx x 5) + 0.2 X dx = 0 It 2 d r obtained from this approach i s shown in figure 3-8. .'" For the f i n a l example 2 2 "5 d x -t- m (x + px ) - jU x cos mio^f = 0 2 . where jU = m 6 and mco^ = 1 .0, the solution envelope i s given by: a(x) = 1 + j>u(cos X - 1 ) where U = 1 ] m [4(l+epa )(l+k ) - ~ Q which can be obtained directly from equation (43). Note that U i s i t s e l f a function of 'a' (which is assumed to be constant), and that U defines the depth of modulation of the solution amplitude. The following short algorithm 49 may be used to determine U : 1. Set a = 1 2 2 2. Evaluate U, taking k = pa 2(H Pa^) -3. Calculate an improved value of the average amplitude 'a' from the relationship a = 1 - p-U. 4. Return to step 2, and repeat u n t i l the desired accuracy i s obtained. The algorithm converges rapidly, e.g. three iterations were sufficient to give three figure accuracy in the example subsequently considered. As the average solution amplitude may be assumed constant, the phase of the approximation solution ( can be determined from an integration of the expression for y/ .. From equation (29): B ^ ^ t A ^ (A 2N1/2 (1 + pa ) J f = ( H p a 2 ) l / 2 + B Substitution for A ^ and B^ ^ yields: • y, = ( 1 + p a 2 ) 1 / 2 _ 2B(1 + P a 2 ) l / 2 (l + k 2) cos cont [ 4 ( l 4 e p a 2 ) ( U k 2 ) - co Q 2] Note particularly that the term with a *t' multiplier has disappeared. The expression for y/( X) can now be obtained as: 2 ( U p a 2 ) l / 2 ( l 4 - k 2 ) s i n r y/(tr) = m(U Pa 2) 1' / 2t - ^ LU CO 0 where 'a' i s the average solution amplitude defined by a - 1 - JJ.U, and k corresponds to this value of 'a'. The local maxima, minima and zeroes of the solution are determined from a graphical solution of the equation nK = y/( X ) where n i s an integer (n=0,1,2, ...) and K is obtained from Appendix 2, Table 1 2 for the value of k corresponding to the average amplitude. A solution of the equation d 2x + 5(x + 2x 5) - 5 x cosX = 0 i s shown in figure 3.9-3.5 Discussion A comparison of the results of the unrefined and refined approxi-mations indicates that substantial improvement in both amplitude and phase accuracy can be obtained when the refinement i s employed. The f i r s t example (where f ( x , i , t ) = tx) i s the case most susceptible to phase error, and, as such, i s a severe test of an approximation method. It i s encouraging that the approximate solution for this example should demonstrate such a noticeable improvement in accuracy, but perhaps the most interesting aspect of the refine-ment i s the connection between amplitude and phase variation which was demon-strated algebraically in section 3.3 and discussed in section 3-4. 3.6 Conclusion The refinement of the f i r s t approximation of chapter 2 i s mathema-tically,no more d i f f i c u l t to apply, and yields approximate solutions which are significantly better than those derived by the un-refined method. The results of the error analysis show that the greatest improvement i s obtained at high values of the parameter p (which determines the degree of non-linearity There are three aspects of the refined approximation of particular significance. a) The use of e l l i p t i c functions i n the representation of.non-linear o s c i l -lations results in a zero-order approximation which i s versatile and inherently more accurate than a corresponding circular function approxi-mation. b) Consideration of solution-frequency variation with amplitude leads to the 2 introduction of terms Involving k , a parameter which varies with the degree of non-linearity and which can therefore accommodate, with accuracy, 51 wide variations in non-linearity. This feature is- unique to the present analysis, and i s an important contribution to the accuracy of the method. c) The equation defining ^ requires no assumption of solution frequency. A ^ i s obtained as the particular integral of a linear conservative dif f e r e n t i a l equation, and should therefore, on average, be an exact representation of the behaviour of the amplitude envelope. This is .in contrast to either the Kryloff-Bogoliuboff method [13 ] or the Bogoliuboff Mi.tropolski method [ 4 ], which both assume the. amplitude and frequency remain constant when deriving the expression for a(t). Although in the present method a constant amplitude may be taken when deriving.the phase-variation term, neither constant amplitude or frequency need be assumed when calculating the amplitude envelope. In the following chapter this refinement i s extended to consider equations of the form x+-x-px+-Bf(x,x,t) = 0, which relate to systems exhibiting saturation or limiting phenomena, e.g. electronic circuits involving saturating amplifiers, inductors and capacitors, negative-resistance devices, or oscillatory motion with a restoring force F 3 of the form F =.ax - bx . For the parameter p i n the range 0 < p < 1, and 6 =0, this equation has an exact solution in terms of the e l l i p t i c sine function. 52 4. SATURATING NONLINEARITIES 4.1 Introduction The a n a l y s i s , so f a r , has been concerned with the transient response of non-resonant nonlinear systems where the parameter p i n the equation ii + x + px +. Bf(x,x,t) — 0 i s p o s i t i v e , i . e . the n o n - l i n e a r i t y i s of the "hardening" type [ 6 ]. Many systems of p r a c t i c a l i n t e r e s t , however, contain non-linear elements of a "softening" type [ 6 ] (e.g. saturation e f f e c t s i n inductors and capacitors) [ 17 ], negative resistance devices i n e l e c t r o n i c c i r c u i t s [ 6 ] ) , and may be described by an equation of the form x +- x - px 5 +• Bf (x,x,t) = 0 (p > 0) This equation, f o r 6 = 0 , i s s t i l l s a t i s f i e d by an e l l i p t i c function, 2 but the form of the s o l u t i o n i s now dependent on the quantity pa , where 'a' i s the s o l u t i o n amplitude [ 21 ] . . I f pa i s l e s s than 1.0, then the s o l u t i o n i s o s c i l l a t o r y and i s given by x(t) -- a Sn(iot + K) when x(0) = 1.0, x(o) = 0; otherwise ( i . e . , i f pa > 1) the s o l u t i o n i s unbounded. The present a n a l y s i s w i l l be confined to a consideration of o s c i l l a t o r y motion, and, as before, resonance e f f e c t s are.required to be n e g l i g i b l e or non-existent. The parameter p (p > 0) i s retained, to f a c i l i t a t e compari-son with the r e s u l t s of chapter 3; when a more general d i f f e r e n t i a l equation i s considered i . e . , where the n o n l i n e a r i t y can be e i t h e r hardening or softening , the following form i s used: x + x + r x + Bf(x,x,t) = 0 where -1 < r < oo. 4.2 Development of the f i r s t approximation From Hsu [11 ], the exact s o l u t i o n of the non-linear d i f f e r e n t i a l •equation 53 •2 x + x - px = 0 (p > 0) with x(0) = 1.0, i ( 0 ) = 0 i s : x(t) = a Sn(wt + K) where co = 1 - pa_ 2 1/2 , K i s the complete e l l i p t i c integral of the f i r s t kind, and the modulus of the e l l i p t i c function, i s now given by ,2 2 k = psi 2 [ l - PJL 2 Accordingly, when p ^ 0 i n the equation x + x - px + 8f(x,x,t) = 0 let the solution for x(o) = 1, x(o)'= 0 be taken as x(t) = a(t) Sn('cot + 9 + K) = a Snj/ As before: and hence X = x(a,y/) x = a. 3x + \J/c> x ha ' hf x = a cbc + a d_ da dt ^x .3 a 4- y d 4- y 3x_ dt L dyj Sy Evaluating the derivatives i n this expression: 5x = Sn*J/ ha d_ dt d_ dt ax Sx" 2c„3 = yCny Dny •aCny Dny + a y J ^ (CnyDny) where _3_ (Cny Dny) = -(Sny - pa Sn^y ) ay and ^ x = a Cny Dny ay ' ' 1 ~PiL 2 By substitution, equation (45) may now be written as: . 2 2 *5 x = aSny 4 2ayCnyDny - a y (Sny -pa Sn y ) + ayCnyDny 1 - p a 2 2 (44) (45) (46) The polynomial representation of a and 6 i s unchanged from chapter 2, i.e. a = |*A ( 1 ) (a,t) + y 2 A ( 2 ) (a,t) ^ A ( i ) (a,t) 0 = j U B ( l ) (a,t) 4 - ^ 2 B ( 2 ) ( a , t ) + ^ B ( i ) ( a , t ) + ... and, r e t a i n i n g only terms of 0(jx): where the subscripts denote p a r t i a l d i f f e r e n t i a t i o n with respect to time. Following the approach of the r e f i n e d approximation (where the v a r i a t i o n of s o l u t i o n frequency with amplitude i s taken i n t o consideration): y/ = to + til) •+ 9 where co = 1 - pa_ 2 1/2 and To f i r s t co = -paa 1 - pa 2 . 1/2 order terms i n jJL, the expression f o r then becomes: r = 2 1 - pa 2 . 1/2 pat A (1) V 2 1-pa . 2 . 1/2 S i m i l a r l y : 1 - pa 2 . 2B (1) 1 - pa 2 . l/2 . A ( 1 ) ' - pat A x y (47)" (48) 2 2 2 2 where terms of 0 ( JA. ), 0(y t) and 0(yu t ) are neglected. For 55 or f e paa - pt, Vp2 2 1 - pa 2 J aa 1 - pa 2 TJ2 1 - pa 2 . 3/2 .2 Nov? the term a i s of order u. , and hence to f i r s t order i n f = r B ( l } - Pa A ^ 1 - pa_ 2 172 - pat A^ ^ 1 - pa 2 J 1/2 (49) • . 2 Substitution for y, y and y in equation (46) yields the following expression for x : x= -(aSny- -pa? Sn^y ) 4- JJL CnyDny 2A (1) 2 1 - pa_ 2 1/ 2 + a B ^ - (2pa 2 A< 1> t J>a 2t A ^ ) | 1 - pa_ ..2. 1/2 A^-'Sny/- (Sny- pa 2 Sn^y ) 1 - p a 2 2 . 2aB (1) , 2 1 - pa 2 1 7 2 - pa 2t A<1~>' (50) But, from equation . (44), x may also be expressed as x = -(aSny- pa?Sr?y + 8f(x,i,t)) and, setting 6 = jU, , i t now follows directly from equation (50) that f(x,x,t)= - Cn^Dny 2A (1) 1 - pa 2 1 / 2 + aB - (2pa 2 A^K pa 2t A ^ ) | 1 - pa_ 2 ,(0 2 e 3 1 - pa 2 ) W1> " 1 2 " 1 - pa 1 / 2 - p a 2 t A<1>" 2 (51) Following a similar approximation of chapter 2 let r 2K r o s r r J 2K 0 2.3, Sn 3ydy (52) Now d_ (CnyDny ) = -Sny (1 + k ) 4- 2k Sn^y 56 and hence 2k 2 J Sr?y dy = CnyDny -+ (1+k2) J S n y d y . By evaluating the integral of Sny/ [ 5 ] , the expression for o may now be obtained as S - 1 k 2 In (1 + k ) In [ 1 + k ] - k j~1 + k j L 2 L1 - k j (53) 2 A tabulation of p, k , o and K i s to be found i n Appendix 2, Table 2. Integration of equation (52) when k = 0, i.e. when Sny reduces to s i n y , gives the value of b as 2/3. An identical value i s obtained for equation (53) (when k = 0) after the application of L'Hopital's rule. Equation (51) may now be simplified to: f(x,x,t) = -CnyDny' V 1 ) 1 - pa 2 2 ' / 2 + aB - ( 2 p a 2 A < ' ^ pa 2t A<]>) 2ri--2T/2 pa 2 -Sn r A ( l } - ( i ^ i p a ! ) t r 2 . 1 - pa 2 2aB (1) (54) 1 / 2 2, .(1) ' - pa t A w From Milne-Thompson [ 16 ],fche x — i e r series expansion of Cny/Dn^/ may be derived as 00 CnyDny = d (Sny)= _d_ 2TI kK y Z , „ / L s = 0 Li - q . ( B +1/2) : ,(2s + 1) sin (2s -+ 1 ) 2K i . e . CnyDny = Tt"1 kK 00 s^O L1 - q (s + 1/2) (2s + 1) (2s +• 1) cos (2s + 1) _rc_y/ 2K where q = exp ^-liK' , and hence Sny and CnyDn^ may be expressed as infinite sums of sine and cosine terms respectively. As for the previous case (where r > 0 i n the equation 3 x +• x +• rx +• 8f(x,x,t) =0), this property indicates that Sny and CnyDny may be treated as analogous to s i n y and cosy in the principle of harmonic balance. The three examples of chapter 3 are again chosen to demonstrate the a p p l i c a t i o n of equation ( 5 4 ) . Example 1. x'+ x - px + 8tx = 0 with x(o) = 1.0, x(o) = 0 and where f ( x , i , t ) = tx = at Sny/. S u b s t i t u t i n g f o r f ( x , x , t ) i n equation ( 5 4 ) , and comparing c o e f f i c i e n t s of Sn^ and CnyDnij (55) a B ^ + 2A<1> . 2 1 - pa 2 1/2 _ ( 2 p a 2 A ^ + p a 2 t A ^ ) = Q. " 2 1 - pa 2 • 1/2 - (1 - o V ) 1 - pa 2 . 2aB (D , 2 1 - pa 2 . 1/2" 2. .(1)" ' - pa .t A v y = - at (56) Taking the p a r t i a l d e r i v a t i v e of equation (56) with respect to time: .(0 a B ( l ) = U ^ + a ) 2(1 -S pa 2) 1 - pa_ 2 ,0) 1 / 2 + p a f (A^K t A ^ ) 2 X 1 - pa_ 2 -1/2 and by s u b s t i t u t i o n f o r a B ^ ' i n equation (55) i t now follows that A ( J . | 4- 4 A ( 1 ) ( 1 - S p a 2 ) 1 - . pa 4 1 - pa 2 . - a (57) S e t t i n g A ^ = 0 (as i n s e c t i o n 3-2 f o r t h i s example), and noting that 2 2 pa =k_ , where k i s the modulus of the e l l i p t i c function Sn\ls, the • 21 2 ' 1 - pa. 2 p a r t i c u l a r i n t e g r a l of equation (57 ) may be obtained as A<1> = -a 4(1 - Spa" and then, by i n t e g r a t i o n : a ( t ) = exp 4(1 -S 1 - £ 2 - - P t 2 N pa 1 - £ 2 (1) (58) ,(D Taking A ^ = 0 i n equation ( 5 6 ) , the expression f o r B becomes 58 B < 1 > - . t r\ 2 i 1 - pa 1/2 2, - pa t . 2(1 - S pa 2) 2 . 8(1 - opa 2) r 1 - k 2 i 2 1 - pa_ 2 •1/2 and f i n a l l y : e(tj = B£ 4(1 - Spa2) .1 - pa 2 1/2 16(1 -6 2,2 a t ; 2 1 - pa S pa 2) n - k 2" 2 2 _ -1/2 (59) Example 2. x +- x - px +• B t i = 0 with x(0) = 1.0, i(o) = 0 and where f ( x , i , t ) = t i = at 1 _ E§_ 2 1/2 Cn^Dn^ + 0(Bt) As before (see section 3 .2) , terms of 0(pt) i n f ( x , i , t ) are equivalent to 2 • terms of order 6 t i n Bf(x,x,t), and may be neglected. From equation (54)., comparing c o e f f i c i e n t s of Sny and CnyDny : 1/2 1 - pa 2 2 J 1/2 2 .(1) 2, -*.(l)x - (2pa A v y+ pa t A y ) 1 -at 1 - Pa' 2 T/2 1r-pa2 (60) A ( l } - ( l ^ i p a ! ) ~ i 2 1 2 ~ 2 a B ^ 1 - p a 2 1 / 2 - P a 2 t A<1>" = 0 .. (61) 2... A f t e r taking the p a r t i a l d e r i v a t i v e of equation (61) with respect to .time, and substitution, f o r a B ^ i n equation (60), the fol l o w i n g equation i s obtained: A{1\ + 4 A ( 1 ) (1 -S pa 2) 1 -x. 2 = -2at (1 - & pa 2) . The p a r t i c u l a r i n t e g r a l of t h i s equation i s A<1>= -at and hence, by i n t e g r a t i o n : a(t) — exp 1 - kT 2 1 - k' 2 21 (62) 59 (1) S u b s t i t u t i o n f o r A ^ i n equation (61) then gives 21-1/2 1 - pa ' 2 T»0) 2 4 - 2 B = -Pa t 1 - £ 2 1 4(1 -S pa 2) h - k 2 l 2 1 - pa_ 2 1/2 The same d i f f i c u l t y i s found here as was encountered i n section 3.2, namely that a simple closed form expression f o r 9(t) cannot be found i f the amplitude v a r i a t i o n i s to be taken i n t o consideration.- However, as f o r the e a r l i e r case, i t i s s t i l l possible to determine s o l u t i o n phase by an approxi-mate method (see sections 3-4 and 4.4). Example 3. - 3 , x - f x - px - Bx cos to^t = 0 with x(0) =1.0, x(o) = 0 and where f( x , x , t ) = -x cos co^t = -aSn^cos co^t . Applying equation (54): 2 ' - t^pa A ' • -f 2" aB<}W1> 1 - pa_2 1/2 /_ 2 . (1 ) 2, . (1 ) N 0 (2 v ' -f pa t A v ') =0 11/2 t 2 A ( l } - (1_-J$paf) t 2 1 - pa 2 " 2aB<1> 2" 1 - pa 2 1 / 2 - P a 2 t A ^ ' = a cos co t . o (63) By the same procedure as before: A ( 1 ) _ - a co Q s i n c o Qt 4(1 -6 pa 2) .(0 1 - k: 2 - CO '0 and, a f t e r s u b s t i t u t i o n f o r A ^ i n equation (63): 4(l-«5pa2) [1-k 2] 21 2 1 - pa 2 _ _ 2 . _ 2 pa tco sinco t T/2 + 2 1 - pa 2 1/2 (64) 1 - £ 2 cos co^t F i n a l l y , i n t e g r a t i n g equation (64), and assuming 'a' to be constant: 60 i(t) = exp 6 (cos co^ t - 1 ) 4(1 -<5pa2) or, i f the exponent i s small: k_ 2 - to,. i(t)=-1 + r p ( c o s v -1) _4(1 - opa 2) |~1 - k*_ - co 0 2 (65) As discussed in section 3.2, i t i s unnecessary in practice to integrate the expression for B ^ and determine 6(t). Two limitations must be imposed in this particular case: 2 l™~ 2 2 a) To avoid resonance effects 4(1 - o* pa ) 1 - k_ »to L 2 J b) For solution s t a b i l i t y , the i n i t i a l value of the expression x - px - 8x cos to^t must be positive [ 11 ]. when x(o) = 1, this last condition may be written as 1 - p > 8. The equations developed above are applied to three specific examples i n section 4«4> following an analysis of error incurred i n applying the approxi-mation. One f i n a l detail, however, remains to be verified. In section 3.3 i t was shown that, to a first-order approximation in the case of r > 0 i n the equation x 4- x 4- r x 3 4- 6f(x,x,t) = 0, the change in j/ due to j^_(to) was 5a cancelled by that due to ^_(o). The first-order expression for u> i n the present case i s given by equation (47) as: f 1 - p a 2 1/2" + P B<1> - nat A<1> 1 - pa_ 2 Considering, for purposes of demonstration, example 3> and substituting for A<1> and B ( ' > : f 2 P§_ 2 1/2 [4(1 - p a 2 ) [ j co. pa t sin to^t I1 " ^ "T72 • on 2 l 1 / 2 4- 211 - pa | ' 1 - £ 2 , • pa t co- sin co_t cos to^t - 0 0 T72~ A similar cancellation to that observed in section 3 - 3 is thus also obtainable for the equation x + x - px + (3f(x,x,t) = 0. The above result may readily be deduced for the remaining two examples. In the following section the method introduced i n section 2 . 3 i s applied to determine error incurred when applying the approximation. 4 . 3 Error Analysis The analysis of the preceding section yields, for a specific example relationships defining the quantities, a, a, a, yv and + a^Cn^Lny/ ( 6 6 ) * 2 2 "5 > x - a Sny + 2ay»CnyDny/ - ay/ (Sny> - pa Sn y ) + a^CnyDny (67^ "l - p a f 2 Now x = a Sny i.e. Sny = x/a and from equation ( 6 6 ) : Substituting for Sn f x - a Sny = ayCnyDny Cny'DnU/ = k - ax . ' ' ay/ 2 • T a y> and hence, substituting for S n y and Cny/Dny i n equation ( 6 7 ) : x = a.x + a x - ax a _y£_+ 2 a a I (x - Px 3) 1 - pa_ 2 ( 6 8 ) A l l the coefficients of equation ( 6 8 ) are known (for a particular differential equation), and the approximate solution can now be obtained by integration of 62 the s t a t e equations X1 ~ X 2 x 2 = a x 1 -t a a_ x. _y) + 2 a c. I a _ y a 1 - pa 2 . These equations were i n t e g r a t e d (using.a d i g i t a l simulation), f o r the three examples of s e c t i o n 4.2, and.the i n t e g r a l e r r o r d e f i n i t i o n of equation (28) was again used to o b t a i n the r e s u l t s shown i n f i g u r e s 4.1 to 4.6. F i g u r e s 4.1, 4.2 and 4.3 r e f e r to the v a r i a t i o n of the e r r o r i n t e g r a l w i t h 6 f o r the three examples, h o l d i n g p constant at 0.5; f i g u r e s 4-4, 4-5 and 4.6 show v a r i a t i o n of the e r r o r i n t e g r a l w i t h the parameter p, w i t h 6 h e l d constant. An i n t e r e s t i n g f e a t u r e of f i g u r e s 4.4 to 4.6 i s that s o l u t i o n e r r o r i s shown to increase w i t h the parameter p ( c f . f i g u r e s 5-4 to 3-6 where the 2 2 converse i s seen to be t r u e ) ; note, however, that k = pa f o r t h i s case, 1 -pa 2 p tends to 1.0 as the modulus k tends to 1.0, and then the e l l i p t i c f u n c t i o n Sny tends to tanh y, which i s n o n - o s c i l l a t o r y . In s e c t i o n 3-4, where the equation under c o n s i d e r a t i o n was x + x + px +- B f ( x , x , t ) = 0 and p > 0, the r e s u l t s shown i n f i g u r e s 3-4 to 3-6 i n d i c a t e d that s o l u t i o n e r r o r decreased w i t h i n c r e a s i n g p, although i n that case the s o l u t i o n p e r i o d was a maximum f o r p = 0. The l i n k i n g f a c t o r of the two cases (which,, between them, cover the e n t i r e range of o s c i l l a t o r y s o l u t i o n s of the equation' x -+ x + px +Bf(x,x,t) = 0 where -1 < p 0 ) (70) This particular example i s also used to demonstrate application of the refined approximation of chapter 3 to systems where the i n i t i a l amplitude i s zero, but where x(o) i s non-zero, and where the amplitude envelope i s defined by an exponential function with a positive argument. The main d i f f i c u l t y encountered by previous workers when investigating solution behaviour of the unforced Duffing equation x + x + px 3 + Bx = 0 (71) was the prediction of solution phase (which, as might be anticipated from the discussion of section 4.6, i s far from being a linear function of time), particularly when calculated from a quasi-linear approach. Considerable ingenuity has been applied to the construction of suitable mathematical models 74 u s i n g c i r c u l a r functions,and Ludecke and Wagner [34] i n a recent paper adopted an approach which was c o n c e p t u a l l y s i m i l a r to that used i n chapter 3 ( i . e . , the v a r i a t i o n of s o l u t i o n frequency w i t h amplitude was taken i n t o c o n s i d e r a t i o n ) , except that c i r c u l a r f u n c t i o n s were s t i l l used i n the development of the approximate s o l u t i o n . In the f o l l o w i n g s e c t i o n i t i s shown that, f o r the damped D u f f i n g equation of equation (71), the phase of the f i r s t order approximate s o l u t i o n can be obtained e x p l i c i t l y i n terms of the amplitude envelope when e l l i p t i c f u n c t i o n s o l u t i o n s are considered. This r e s u l t makes p o s s i b l e a f i r s t order s o l u t i o n of equation.(71) i n a simple form, which i s accurate even f o r com-p a r a t i v e l y l a r g e values of the parameter 6. 5.2 L i n e a r d i s s i p a t i o n In s e c t i o n s 3-3 and 4.2 i t was shown that a p a r t i a l c a n c e l l a t i o n of f i r s t order terms i n the expressions f o r \j/ (equation (29) of s e c t i o n 3.2 and equation (47) of s e c t i o n 4 .2) was obtained f o r non-autonomous f u n c t i o n s 3 f ( x , x , t ) i n the equation x + x + r x + B f ( x , x , t ) = 0 where -1 < r < oo. Now consider an autonomous f u n c t i o n f ( x , x ) i n the equation -x + x - f r x 3 -+ 8f(x,x) = 0. (72) I f r > 0 ( i . e . the case considered i n chapter 3) the approximate s o l u t i o n i s of the form x ( t ) = a(t)Cny/, and f ( x , x ) may be expressed as 2 1 /2 f ( x , x ) = -g 1 a. (1+pa ) Sny/Dny/+- g ? a Cny, where f o r the moment no r e s t r i c t i o n i s placed on the f u n c t i o n s g^ and g^* except that they should not be f u n c t i o n s of time. Comparing c o e f f i c i e n t s of SnyDny and Cny when f ( x , x ) i s s u b s t i t u t e d i n equation (33) g i v e s (by the method of, f o r i n s t a n c e , example 1 of s e c t i o n 3-2): A ^ + 4 A ^ (l + e p a 2 ) ( U k 2 ) = -2g a(l+epa 2) 75 since -1 a t = 0, and the p a r t i c u l a r i n t e g r a l of t h i s equation i s : where . 2 2 k = pa 2 ( 1+ Pa 2) • I f f(x,x) i s a non-linear function then, although f(x,x,) can be expressed i n the l i n e a r i z e d form f(x,x) = - g l a(l + p a 2 ) 1 / / 2 S n ^ D n y / + g g aCn^ (see section 5-3), i t i s not s t r i c t l y true that g^ ^ g ^ ( t ) , as g^ w i l l be the approximation of a term i n v o l v i n g a mean value and pe r i o d i c component. I f , however, f(x,x) i s l i n e a r i n x and i , then g^ and g^ w i l l be constant and w i l l be zero. (1) (1) Evaluating B under the assumption A ^ — 0, and s u b s t i t u t i n g f o r A ^ ^ and B^ ^ i n equation (29) gives the f i r s t order approximation to \js as: f = ( U P a 2 ) 1 / 2 + P-g2 ^ - f 2 , ' 2 ( U e p a ) ' i . e . i f f(x,x) = g^x, where g^ i s a constant, then f i r s t order terms i n the expression f o r y/ cancel exactly. A corresponding a n a l y s i s f o r -1 < r < 0 (the case of chapter 4) where 1/2 f(x,x) = g 1 a 1 2 1 - pa 2 C n y / D n y / + g 2 a S n ^ / gives, under the same assumptions of linearity, 1 / 2 , Pg„ f , 2 1 - pa 2(1-Spa d ) 1 -pa 2 2 1 / 2 * 0 ( p 2 ) (from equation (47)). I t therefore follows that, f o r equation (72) and -1 < r < 00 , f i r s t order terms i n y/ cancel exactly; the s i g n i f i c a n c e of t h i s r e s u l t may be seen 76 from the fo l l o w i n g discussion. . Consider an approximate s o l u t i o n x(t)'= a(t) Cny ( i . e . the case of r > 0 ) . Then, i f y = ( l + p a 2 ) 1 / / 2 : d_ (Cny) = -Sny/Dny ]h£ dt & t and may be obtained by i n t e g r a t i n g yv = (1+pa ) t r e a t i n g 'a' as constant. I f the i n i t i a l value of the phase i s 9(o), then x(t) = a(t) Cn [ ( 1 + P a 2 ) t + G(o)] i . e . the phase-modifying term 9(t) i s time-invariant (to f i r s t order). This furt h e r implies that s o l u t i o n frequency i s only amplitude dependent, and the e l l i p t i c i n t e g r a l K, which defines the quarter-period T by T = K (1+pa ) may be regarded as i n v a r i a n t . This l a s t r e s u l t i s p a r t i c u l a r l y important when d e r i v i n g the approximate s o l u t i o n . A s i m i l a r argument may be applied to solutions of the form x ( t ) = a(t)Sny, f o r which the same conclusions are reached. . The f i r s t order s o l u t i o n of the equation 3 x+-x + rx + 8 x = 0 with -1 < r < oo , and an i n i t i a l amplitude defined by the i n i t i a l conditions x(o) and x(0) (see section 5.3), can then be obtained by c a l c u l a t i n g the amplitude envelope a ( t ) , and determining the l o c a l maxima, minima and zeroes of the s o l u t i o n from y (t) = nK(0) where K(o) i s the value of K corresponding to the i n i t i a l amplitude a(o) 2 • (derived from the value of k f o r a(o)), and y i s given by 77 f = ( l 4 P a 2 ( t ) ) 1 / / 2 t + e(o) f o r r > 0 and 1 - p a 2 ( t ) 1 / / 2 t + 9(0) f o r -1 < r < 0. 2 Although t h i s method r e l a t e s d i r e c t l y to the equation considered i n s e ction 5-4, note that, because of the nonlinear nature of the function (1 - ^ x ) i i n equation (70), i t cannot be applied to s o l u t i o n of the modified Van der Pol equation. 5.3 The Modified Van der Pol equation Consider the equation x + x + px 3 - 8(1 - # x 2 ) x =0 (p > 0) (70) where x(o) = 0. I f the approximate s o l u t i o n of t h i s equation i s taken as x(t) = a(t)Cny • where y = cot - K + 9 ( t ) , K i s the complete e l l i p t i c i n t e g r a l of the f i r s t 2 •] /2 kind, co = (1 + pa ) ' and 9(0) = 0, then: x = aCn y - aySnyDny/, 2 2 2 Now Sn(-K) = - 1, and since Dn u = 1 - k Sn u [ 5 ] i t therefore follows that - t „\ . /, . 2x1 /2 • ' , , 2 2 jjn^-K; = \\ - k ) where k = pa .. 2 ( H p a 2 ) I t i s now possible to define x(o) i n terms of the q u a n t i t i e s p, a and k 2 as: i ( 0 ) = a ( l + P a 2 ) l / 2 ( l - k 2 ) l / 2 (73) i . e . a choice of the i n i t i a l c ondition i ( 0 ) determines the i n i t i a l s o l u t i o n envelope anrplitude, the e l l i p t i c f u n ction s o l u t i o n of equation (70) then being defined by x(o) = 0 and x(o). Note that t h i s argument also applies to the a n a l y s i s of chapter 3, and may obviously be extended to that of chapter 4 where x ( t ) = a ( t ) Snf. These considerations of i n i t i a l conditions are necessary to obtain a f i n a l s o l u t i o n , but do not a f f e c t the d e r i v a t i o n of the equation de f i n i n g 78 the functions and B ^ ^ ( i . e . equation (33)) For equation (70): f ( x , x , t ) = -(1 - J-x 2)x and therefore f ( x , x , t ) = aSnyDny/ (1 - ^ a 2Cn 2 y/) ( 1 + p a 2 ) 1 7 2 + 0(p) (74) As before (see se c t i o n 3-2), terms of 0(B) i n f( x , x , t ) are equivalent to 2 terms of 0({3 ) i n 8f ( x , x , t ) , and may be neglected i n the f i r s t order a n a l y s i s . Equation (74) i s of a more complex form than has been considered i n 2 previous chapters, as i t contains a term Cn y Sny/Dny and may no longer be expressed i n the simple form f ( x , x , t ) = g(a,t)SnyDny> which had made possible a comparison of c o e f f i c i e n t s of Cn^ and Snyv Dny i n equation (33) (section 3.2, 2 Wow Cn y From Milne-Thompson [16]: 2 2 3 1 - S n y , i.e.. Cn yv Sny/Dny = Sny Dny - Sn f Dn\fJ co 00 3 3 Sn y Dny/ = 8 T T 3 3 / S s i n ( 2n+-l)7i\i/ v a ' n + 1 / 2 where S " q n 1 - q 2n+1 n=0 > d = _ q n ^ 7 1 + 2 7 1 d • cos ,.n TX uv 2K K n=1 and q = exp -1 + q 2n [~-7lK' L K _ I f k i s small only S Q and d^ w i l l be s i g n i f i c a n t . Using Hermite's expansion f o r q i n terms of k [ 7 ] and taking K = n 2 1 + k f + 3 k l + 4 8 (from Bowman [5 ] ) , the expression f o r Sn y Dny may be reduced to the form Sn y/DnU/ = ( l-A)sin TUJ/ ' 2K (75) where 1- 1 +• 0.9375k' L1 +- k 2 2. - 0.0625k 4 (76) i 4\ terms of 0(.k ) have been neglected, and only the fundamental component of the Fourier s e r i e s expansion has been retained. I f the following approximation i s made: S n 5 y D n y = (l-A)Sny^Dny/ (77) then f ( x , x , t ) can be expressed i n the form f ( x , x , t ) = f(a)SnyDny, and the method of sec t i o n 3-2 applied to determine ^ and From equations (74) and (77) : f ( x , x , t ) = TaSny/Dny ( 1 + P a 2 ) 1 / / 2 where T=1 - v A a , . ' (78) and from equation (33), comparing c o e f f i c i e n t s of Cny and SnyDny: n ,, 2x1/2 ,,(1) 0 , ( l ) / , 2N1/2 Ta(1+pa ) ' - aSr^ + 2k> y (1+pa ) 1 - 2pa 2A^ 1^ + p a 2 t A^ l ^ 1 + pa 2 1/2 (78) 0 = A ^ - ( U e p a 2 ) (1+pa 2) _ 2 a B ( l ) ( H P a 2 ) 1 / 2 + 2 P a 2 t A ( 1 ) (80) , Taking the p a r t i a l d e r i v a t i v e of equation (80), s u b s t i t u t i n g f o r (1 ) 2 2 aB i n equation (79) and s e t t i n g k = pa_ 2(1+pa) •+ 4A^ 1 ^ ( H e p a 2 ) ( l + k 2 ) = 2 T a ( n e p a 2 ) . The p a r t i c u l a r i n t e g r a l of t h i s equation i s : A<1>- P a 2(l+k 2) and therefore, s u b s t i t u t i n g - f o r P from equation (78.): a. - B_a + BxAa = 0 2(l + k 2) 2 ( H k 2 ) (81) Equation (81) i s i n the form of.a B e r n o u l l i equation [ 6 ], which may be integrated a f t e r the s u b s t i t u t i o n -2 z = a The f i n a l r e s u l t f o r a(t) i s : 80 i(t) g-A •+- exp L 1+k JL - *A 2 *0 -1/2 (82) where i s the i n i t i a l value of the amplitude envelope (from equation (71)). 0) As P i s , not an e x p l i c i t f u n c t i o n of time A ^ = 0 and hence, from equation (80): B<1L - P p a 2 t 2(l + k 2 ) ( U p a 2 ) l / 2 or, s u b s t i t u t i n g f o r T from equation (78): ,0) 2, - pa t + fl- pa tA 2 ( l + k 2 ) ( H p a 2 ) l / 2 2 ( H k 2 ) ( l + P a 2 ) l / 2 Integrating t h i s expression, and t r e a t i n g 'a' as a constant [13]: e(t) - p p a 2 t 2 4(l + k 2 ) ( U p a 2 ) l / 2 #Aa -1 (83) Equations (82) and (83) now define the f i r s t order approximate s o l u t i o n . To reduce these r e s u l t s to apply to Van der Pol's o r i g i n a l equation (see [30 ] and equation (69)), l e t ^ - 1 and p = 0. Then, i f x(o) = 0 and'x(o) = a, the approximate s o l u t i o n of equation (69) i s given by: :(t) = j_ +- exp(-Bt) 4 1 - 1 * 2 4 0 -1/2 s i n t (84) since 6(t) = 0 . The equation of the stable o s c i l l a t i o n i s obtained by l e t t i n g t tend to i n f i n i t y ( i . e . exp(-Bt) 0), and hence i n the l i m i t : x(t) = 2 s i n t. The r e s u l t 9(t) = 0 was also found by K r y l o f f and Bogoliuboff [ 1 3 ] , and produces a c i r c u l a r t r a j e c t o r y of the l i m i t cycle i n the phase plane (x-x) f o r equation (69). In p r a c t i c e the t r a j e c t o r y i s of a more complex nature [31], but i t should be noted that the present analysis (and the K-B analysis) produces averaged r e s u l t s which do not take into consideration the second-order e f f e c t 81 of amplitude v a r i a t i o n (e.g. the term 0(s) i n equation (74))- A n .analysis of Van der P o l ' s equation by the q u a s i - l i n e a r method of B o g o l i u b o f f - M i t r o p o l s k y can be found i n Minorsky [ 17 ]. A f t e r some manipulation i t can be shown that equation (84), i s e q u i v a l e n t to the f i r s t order approximate s o l u t i o n i n [ 1 7 ] . For the modified Van der P o l equation (equation (7 0 ) ) , the a p p r o x i -mate s o l u t i o n i s , from equations (82) and (83) :(t) #A + exp n 2 4 - 2 Bpa t -8±. U k 2 -1/2 Cn ( U p a 2 ) l / 2 t - K 4 ( H k 2 ) ( l + p a 2 ) l / 2 • Aa - 1 (85) where x(o) = 0, and i ( o ) i s de f i n e d by equation (73) w i t h a = a^, the i n i t i a l envelope amplitude. A s o l u t i o n of the equation x + x + 2 x 3 - 0.1 (1 - x 2 ) i = 0 , ( 8 6 ) u s i n g the f i r s t approach of s e c t i o n 3 - 4 , " i s shown i n f i g u r e 5.1. The i n i t i a l amplitude a n was chosen as 1.0, g i v i n g (from equation (73)) i ( 0 ) = 1.4.14. The parameters of equation (86) correspond to 6 = 0 . 1 , ^ = 1.0, p = 2.0, k 2 = 0.3333, K = 1 .734 and hence (from equation ( 7 6 ) ) : A = 0.283 . S u b s t i t u t i n g these values i n t o equation (82), the amplitude envelope becomes: t(t) = 0.283 + 0.717 exp (-0.075t -1/2 and the amplitude of the s t a b l e o s c i l l a t i o n a i s , t h e r e f o r e , s a s = ( y - A ) - l / 2 = ( 0 . 2 8 3 ) - l / 2 = 1.88 . 2 2 As a _ J , the term ( >fAa - 1) i n equation ( 8 3 ) vanishes when S ——— \j the s t a b l e o s c i l l a t i o n i s reached, and the s o l u t i o n of equation (86) then 82 0.4 ttt) • 1.4 ]o.?83 * 0.717 expt-0.07S f}J J H H 2.0 4.0 —I— 1.0 H f h --4-^ 1-(00 -a/ -i.o-4- - — ^ Time (seconds) computer solution approximate solution — — K-B solution a a a o F i g . 5.1 Approximate and exact s o l u t i o n s of . . x -f x + 2 x 3 - 0.1 (1 - x 2 ) x = 0 X(l) F i g . 5.2 Phase p o r t r a i t of the s o l u t i o n to x -f x + 2 x 3 - 0.1 (1 - x 2 ) i = 0 becomes- c ' . x(t) =• 1 .88' Cn[ (8.06) 1./ 2 t - I . 7 3 4 ] . (87) Differentiating equation (87) with respect to time, the maximum value of x(t) i s obtained from the equation x = a. Cny - ay/Sny/Dny . When the stable o s c i l l a t i o n i s reached, a. — 0 and hence . x = a y Sn \» Dn y max ' / ' max 2 2 For small values of k (i.e. k « l ) |Sny/Dny I w i l l be close to i t s maximum value when \j/ = (2n+1 ) K, and then a g j/ [ 1-k 2] (88) x max o 2 2 pa where k = -^T - s - (89) 2 [ l + P a 2 s _ Finally, from equations (87), (88) and (89): i 4 . 0 0 . max The phase portrait for equation (86), obtained from an analog computer simulation, i s shown i n figure 5 . 2 , where the limit cycle behaviour may clearly be seen. Actual values of x and x (from this simulation) max max were 1.92 and 4 . 1 5 respectively, which compare favourably with the predicted value of 1.88 and 4 . 0 0 . 5 . 4 The damped Duffing equation For the Duffing equation x 4- x + px 5 + 8x = 0 (p > 0 ) (71.) with x(o) = 1 . 0 and i ( o ) = 0 , the function 6f(x,x,t) i s , to a f i r s t approxi-mation Bf(x,x,t) =' Bx = -8a SnfDnf ( 1+pa 2) 1^ 2 . Following the method of section 3-2, and applying equation (33), 84 (1) the expression.defining A may be obtained as + 4 ( l + e p a 2 ) ( l + k 2 ) = -2a (Uepa 2) (83) where k 2 = pa 2(1+Pa ) Hence, taking the particular integral of equation (83): (1) - a 2(Uk 2) and t(t) = exp j:Bt_ 2(Hk 2) (84) From the analysis of section 5-2, the f i n a l solution for x(t) i s therefore :(t) = exp " -Pt Cn 1 + p exp ^ Bt_" 1/2 t -- 2 ( l + k 2X J+k 2. since 9(o), the i n i t i a l phase, is zero (for the choice of i n i t i a l conditions). The local maxima, minima and zeroes of the approximate solution can be obtained by a graphical solution of the equation 1/2 1 -t- p exp -P i . i + k s n K (0) (n = 0,1 ,2, ). A solution of the equation x + ' X + 2x 3 + x = 0 with x(o) = 1, i ( o ) = 0 i s shown i n figure 5-3, where although the parameter 8 i n equation (71) is equal to 1, solution phase accuracy, i s seen to be maintained. 85 1.0 | 0-4-x(t) V B . \ . exp -0.375t 0 * a _ • • _ i e <• a """" —. '. o{—• j o—+ o .2.0 p _ 1 --^..Kf ] £2- — — ~~8.0 10.0 a • % a Time {seconds) computer solution approximate solution — K-B solution o o a a F i g . 5-3 Approximate and exact s o l u t i o n s o f x + x + 2x^ + x = 0. 86 5.5 Comparison with the K-B Method It has not previously been possible to make a comparison between the e l l i p t i c f u n ction approximation and the K-B approximation, as the l a t t e r cannot be applied ,to non-autonomous systems of the type considered i n chapters 2-4 . For the two nonlinear autonomous systems of section 5 .3 and 5-4, however, such a comparison can be made. Consider f i r s t the modified Van. der Pol equation 3 2 x + x + px - B ( 1 - X ) X = 0 with x(o) = 0 , and assume (following the K-B method [ 1 3 ]) a s o l u t i o n of the form x(t) = a(t) s i n (cot + 0 ( t ) ) , where co = 1 . Now the function f(x,x) i n the equation 2 x + co x -+- e f ( x , i ) = 0 . . . . . . !1 assuming e =1 . S u b s t i t u t i n g f o r x and i : 3 2 i s given by f(x,x) = px - B(1 - x )x , f ( a s i n y , a y c o s y ) = pa 3 s i n y - s i n 3 y 4 L + 6 a 2 - 1 L4 c o s y - J L c o s 3 Y 4 ' _ where y = cot + 9 . The f i r s t order K-B approximation equations f o r a. and 0 may then be obtained as [ 13 ] : 2 a = pa 2 1 a. 4 9 = 3 pa 8 For the equation chosen i n sec t i o n 5 . 3 , where p = 2 , B = 0 . 1 , and assuming a(o) = 1 . 0 , the r e s u l t i n g K-B approximation s o l u t i o n i s : 87 :(t) = ^ 2 s i n ( l . 7 5 t ) . 0 . 2 5 + 0 . 7 5 exp ( - 0 . 1 t ) The f i r s t one and a h a l f periods of t h i s s o l u t i o n are pl o t t e d i n f i g u r e 5 .1» where the superior accuracy of the e l l i p t i c f u n ction approximation i s evident. . . • . Now consider the he a v i l y damped Duffing equation considered i n section 5 . 4 : 3 x + x + px + Bx = 0 , with x(o) = 1 and x(o) = 0 . Assuming an approximate s o l u t i o n of the form x(t) = a(t) cos (cot +• 0 ( t ) ) , where again co = 1 , the K-B approximations . f o r a and 9 are: a = -Ba 2 and 2 9 = 3pa . 8 For the equation 3 -x + x + 2 x + - x = 0 with x(o) = 1 . 0 and x(o) = 0 , the K-B approximate s o l u t i o n i s then x( t ) = exp ( - 0 . 5 t ) cos ( l . 7 5 t ) . The f i r s t period of t h i s s o l u t i o n i s shown i n f i g u r e 5-3, and again the e l l i p t i c f unction approximation i s seen to he considerably more accurate. In f a i r n e s s to the K-B method i t should be noted that both of the equations considered here are grossly nonlinear and, as such, are not s t r i c t l y amenable to an a l y s i s by a q u a s i - l i n e a r approach. This comparison does, however, serve to accentuate the misrepresentation of nonlinear system performance when l i n e a r i t y , or even q u a s i l i n e a r i t y , i s assumed. 5.6 Conclusion The a p p l i c a t i o n of the r e f i n e d approximation of chapters. 3 and 4 to non-linear autonomous equations of the normalized form 3 x + x + rx + Bf (x,x) = 0 (-1 < r < o o ) i s demonstrated i n t h i s chapter, by considering a modified Van der Pol equation [17 ] and a damped, forced Duffing equation [34 ]. I t i s shown, i n s ection 5-2, that an e x p l i c i t form of the s o l u t i o n to the equation 'x + x + rx +• Bx = 0 (-1 < r < oo ) can be obtained, making use of a complete c a n c e l l a t i o n of f i r s t - o r d e r terms i n the expression f o r ^ , where i s the argument of the e l l i p t i c f unction used i n d e r i v i n g the approximate s o l u t i o n . This i s a s i g n i f i c a n t r e s u l t , and represents an improvement both i n accuracy and s i m p l i c i t y over e x i s t i n g s o l u t i o n methods. I t i s also shown i n s e c t i o n 5.3, as the d e f i n i n g equations of the r e f i n e d approximation (equations (33) and ( 5 4 ) ) are independent of i n i t i a l conditions,, approximate solutions can r e a d i l y be obtained f o r x(0) = 0 when x(o) i s defined. A comparison of solutions obtained from the e l l i p t i c f u n ction approach and from the K-B method, demonstrates the superior accuracy of the e l l i p t i c f u n ction approximation. The a n a l y s i s of t h i s chapter completes the development and a p p l i c a t i o n of the f i r s t order e l l i p t i c f unction approximation to solutions of the normalized equation x + x + rx + Bf (x,x, t ) = 0 (-1 < r < oo ). In chapter 6 some physi c a l systems described by t h i s equation are considered to demonstrate possible areas of a p p l i c a t i o n of the approximation method. 89 6. SOME APPLICATIONS 6.1 Introduction The theory developed in the preceding chapters relates to second order d i f f e r e n t i a l equations exhibiting a cubic nonlinearity and small perturbing functions. It i s the purpose of this chapter to.show how such equations arise in practice by considering examples chosen from the fields of mechanics, astrophysics, ci r c u i t analysis and control systems. Detailed derivations have been avoided, the intention being to obtain differential equations corresponding to those chosen as examples in previous chapters. 6.2 Environmental studies of mechanical systems As a familiar example of an oscillatory mechanical system, consider a particle, of mass m, suspended on a non-linear spring of the hardening type L 6 J, which exerts a restoring force F = c (x + px ) (p > 0). Suppose that the system i s operating i n an environment subject to rapid tem-perature change, and that the effect of linearly decreasing temperature on system performance i s to be investigated. Contraction of the spring resulting from the temperature decrease [ 19] w i l l produce' an increase in spring stiffness with time. If x i s the displacement, of the particle from i t s equilibrium position, then this increase in stiffness can be represented by a term Btx i n the expression for the restoring force F of the spring: F = c(x + px + Btx) where c, p and B are constant positive parameters, and 8 is small. The equation of motion of the particle can then be written as x + _c (x + px + Btx) =0 m or, taking for simplicity c_ = 1 : m 9 0 x + x + px + ( 3 t x = 0 , which i s the equation of the f i r s t example of chapters 2 and 3. The f i r s t example of chapter 4 applies to a spring with a softening c h a r a c t e r i s t i c [.6 ] , hut subject to the same increase i n s t i f f n e s s as a r e s u l t of decreasing tempera-ture, and f o r which the r e s t o r i n g force F i s of the form F = c (x - px 3 + Btx) Then, i f o =1, the equation of motion of the p a r t i c l e i s m x + x - px 4 ptx = 0 . (0 <: p < 1 ) Now consider the same mass-nonlinear spring system which i s operating i n an environment of varying v i s c o s i t y , such as a l i q u i d which i s r a p i d l y c o o l ing [ 8 ]. I f i n i t i a l l y the v i s c o s i t y i s n e g l i g i b l e , and increases l i n e a r l y with time and decreasing temperature, then the damping force may be taken as A.tx (\> 0 ) . The equation of motion of the p a r t i c l e i s then . x -+- c_ (x -t- rx ) +• X tx = 0 . m m Taking c. =1 and X = B gives, f o r r > 0 , the equation considered i n the m m second example of chapters 2 and 3, and f o r -1 < r < 0 , the equation of the second example of chapter 4 . 6 .3 S t e l l a r pulsation Stars of the Cepheid type e x h i b i t v a r i a t i o n s i n luminosity which are of a periodic nature, and a possible explanation of t h i s phenomenon i s that the surface area of such stars i s varying. I f the s t a r i s assumed to be r a d i a l l y symmetric, qua n t i t a t i v e r e s u l t s f o r the periodic v a r i a t i o n i n radius can be obtained by observing frequency changes i n the emission spectrum of the Cepheid. As the v e l o c i t y of the l i g h t - r a d i a t i n g source i s r e a d i l y determined from such observations, i t i s usual to plot r a d i a l v e l o c i t y as a function of time. Such a c h a r a c t e r i s t i c i s c a l l e d a Cepheid v e l o c i t y curve. 91 The phenomenon of s t e l l a r p u l s a t i o n has at t r a c t e d the att e n t i o n of many .astrophysicists, and W.S. Krogdahl [12 ] proposed a p a r t i c u l a r l y i n t e r -e s t i n g mathematical model. Krogdahl based h i s analysis on the e x i s t i n g i d e a l i z e d models of Cepheids which, a f t e r c e r t a i n approximations, he derived i n the form 2 2 2 3 d_q = -q + 2. q + U q + + 2.(1 - q) d / 3 27 3 L d X . •7 •Here q i s given by q •= g(a , f ) where 'a' i s the mean or equilibrium radius 3 ~ of the star, g(a , t ) being a fun c t i o n which i s assumed small i n comparison 3 with a . From the assumptions that a) Cepheids are s t a t i c a l l y unstable stars b) Real thermodynamic systems are d i s s i p a t i v e and a somewhat h e u r i s t i c argument, he then proposed that t h i s equation be modified by the a d d i t i o n of a term 2 1 - q X* J dq , dtr where jX and X are small empirical constants. The r e s u l t i n g equation i s s e l f - e x c i t e d [ 6 ], and a f t e r a change of v a r i a b l e can be obtained i n the f i n a l form Q " +- Q - 2 X Q2 + H X2 Q3 - JU.(1 - Q2) Q' - 2X (1 - X Q ) ( Q')2 - 0 . 3 27 ' 3 When A. = 0 t h i s equation reduces to Van der Pol's equation (equation (69), and therefore should, f o r small X, e x h i b i t l i m i t - c y c l e behaviour i n the Q - Q' phase-plane. Krogdahl used a graphical technique to determine l i m i t - c y c l e t r a j e c t o r i e s which exhibited behaviour s i m i l a r to that of Cepheid v e l o c i t y curves. To quote H.T. Davis [ 7 ] (p. 371), "the general argument of the author seems to be confirmed". In section 5-3 i t was shown that the amplitude of the stable oscillation of the modified Van der Pol equation was a s ~ Now for Krogdahl's equation, and comparing the terms ^ x jx (from equation 70) and 6(1 - y x 2)i t can be seen that = 1 . - Q ) Q1, Taking, for purposes of approximation, A = 0 . 2 5 , the f i n a l amplitude of the solution Q w i l l be approximately 2. This implies 2 3 that the term J_4 X Q i s not negligible, and that Krogdahl's equation may 27 be written i n the form Q« + Q + J_4 \ 2 Q 5 - \ 27 2 Q +/* (1 - a n d control systems containing both linear and nonlinear com-94 ponents of the form shown in figure 6.2 are generically designated as "problems of Lur.'e" [ 51 ], provided ef(e) > 0 and f(o) = 0. An acceptable nonlinearity would be an odd function f(e) = e + re . x + ^ e f( ) f(e) G(s) y _ Fig. 6.2 The Lur'e type control system. G(s) represents a linear element in operator form, and f( ) represents a nonlinear element. - To investigate control systems which can give rise to nonlinear differential equations of the type analysed in chapters 2-5, consider a servo-mechanism preceded by a saturating amplifier. The servomechanism may be. represented by the transfer function G(s) = _ c 2 s -+- as +b where a, b and c are constant parameters, and the saturating amplifier represented by the function f(e) = e - pe (0 < p <1). Then, i f x(0~) = y(0") and x(t) =0 V t > 0, and since f(-e)= -f ( e ) , the di f f e r e n t i a l equation of the control system shown i n figure 6.2 w i l l be y +ay 4 y(b4c) - pcy = 0 . This i s of the form considered in section 5.4 for the damped Duffing equation, and hence certain nonlinear control systems are amenable to a transient analysis by the methods developed in chapters 2-5-6.6 Conclusion In this chapter some possible areas of application of the approxi-mation method have been indicated by considering four systems which result i n equations of the form x + x + rx + Bf(x,x,t) =' 0 (-1 < r < oo ). An exhaustive survey of practical applications was not intended, but the examples chosen serve to show some areas to which the general analysis pertains 7. CONCLUSION In this thesis a method has been presented for determining approxi-mate solutions to a class of grossly nonlinear, non-autonomous second order differential equations characterized by dfx + m (x +• rx^) + uf (x, dx, x) — 0 (-1 < r