APPROXIMATIONS TO THE FREE RESPONSE OP A DAMPED NON-LINEAR SYSTEM by PAUL TSANG-LEUNG CHAN BrA-..Sc;i, University of B r i t i s h Columbia, 1962 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE In the Department of E l e c t r i c a l Engineering We accept this thesis as conforming to the standards required from candidates for the degree of Master of Applied Science Members of the Department of E l e c t r i c a l Engineering THE UNIVERSITY OF BRITISH COLUMBIA March 1965 I n p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f t h e r e q u i r e m e n t s f o r an a d v a n c e d d e g r e e a t t h e U n i v e r s i t y o f • B r i t i s h C o l u m b i a , I a g r e e t h a t t h e L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and s t u d y , I f u r t h e r a g r e e t h a t p e r -m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y p u r p o s e s may be g r a n t e d by t h e Head o f my D e p a r t m e n t o r by h i s r e p r e s e n t a t i v e s . I t i s u n d e r s t o o d t h a t . c o p y i n g o r p u b l i -c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l n o t be a l l o w e d w i t h o u t my w r i t t e n p e r m i s s i o n . D e p a r t m e n t o f E l e c t r i c a l E n g i n e e r i n g The U n i v e r s i t y o f B r i t i s h C o l u m b i a , V a n c o u v e r 8, C a n a d a D a t e M a r c h 3 0 . 1 9 6 5 A B S T R A C T I n t h e s t u d y o f many e n g i n e e r i n g s y s t e m s i n v o l v i n g n o n -l i n e a r e l e m e n t s s u c h a s a s a t u r a t i n g i n d u c t o r i n a n e l e c t r i c a l c i r c u i t o r a h a r d s p r i n g i n a m e c h a n i c a l s y s t e m , we f a c e t h e p r o b l e m o f s o l v i n g t h e e q u a t i o n „ . 3 x + 2 e x + x + a x = 0 w h i c h d o e s n o t h a v e a n e x a c t a n a l y t i c a l s o l u t i o n , . B e c a u s e a c o n s i s t e n t f r a m e w o r k i s d e s i r a b l e i n t h e c o u r s e o f t h e s t u d y , we c a n a s s u m e t h a t t h e i n i t i a l c o n d i t i o n s a r e x ( 0 ) = 1 a n d x ( 0 ) = 0 w i t h o u t l o s s o f g e n e r a l i t y . T h i s e q u a t i o n i s s t u d i e d i n d e t a i l b y u s i n g n u m e r i c a l s o l u t i o n s o b t a i n e d f r o m a d i g i t a l c o m p u t e r . When e a n d \i a r e s m a l l , c l a s s i c a l m e t h o d s s u c h a s t h e m e t h o d o f v a r i a t i o n o f p a r a m e t e r s a n d a v e r a g i n g m e t h o d s b a s e d o n r e s i d u a l s p r o v i d e a n a l y t i c a l a p p r o x i m a t i o n s t o t h e e q u a t i o n a n d e n a b l e t h e e n g i n e e r t o g a i n u s e f u l i n s i g h t i n t o t h e s y s t e m . H o w e v e r , w h e n e a n d a a r e n o t s m a l l , t h e s e c l a s s i c a l m e t h o d s f a i l t o y i e l d a c c e p t a b l e r e s u l t s b e c a u s e t h e y a r e a l l b a s e d o n t h e a s s u m p t i o n t h a t t h e e q u a t i o n i s q u a s i - l i n e a r 0 T h e r e f o r e , t w o new a n a l y t i c a l m e t h o d s , n a m e l y : t h e p a r a b o l i c p h a s e a p p r o x i m a t i o n a n d t h e c o r r e c t i o n t e r m a p p r o x i m a t i o n , a r e d e v e l o p e d a c c o r d i n g t o w h e t h e r e < 1 o r e ^ l , a n d a r e p r o v e n t o b e a p p l i c a b l e f o r v a l u e s o f e a n d u f a r b e y o n d t h e l i m i t o f c l a s s i c a l m e t h o d s . ACKNOWLEDGEMENT Acknowledgement i s due to a l l who have helped during the course of this work. In pa r t i c u l a r , the author wishes to thank Dr. A.C. Soudack, supervisor of the project, for his invaluable guidance and much-needed encouragement, Dr 0 F„ Noakes, Head of the Department of E l e c t r i c a l Engineering, University of B r i t i s h Columbia, and the staff of the Computing Centre of the University of B r i t i s h Columbia, without whose help this work would have been impossibles The author also wishes to express his gratitude to the National Research Council of Canada for a Bursary and a Student-ship awarded him i n 1962 and 1963 respectively. TABLE OF CONTENTS Page L i s t of I l l u s t r a t i o n s « • • • • » > « • o o o * o <> <> o <> o a « • • • « « o o e • * v AcklioWl edgemeil"fc . . . . • . • • • o . o o e o s o e o o o o e o o o o . o o o o o o . o . YX1 1 s IH"trOdllC"biOIl . > . > . » » , » , » , « » » . . . o o o o o o o o o o o o o e o o o o o o o . * . 1 1»1 Mathematical Models . o o . o o . s o o o a o o o o o o . o a . o . 1 1.2 Ana l y t i c a l Approximations o o o « » . . . 1 1.3 Derivation of the System Equation . . . . 3 2. Study of the System Equation ..«.<>.....<>..•...»«• 9 2.1 Normalization •<> • • . . o o o « . « « . « . « < > < > « . . « . « o . . . • 9 2.2 Phase-Plane Analysis . . a . . . . . . . . . . . . . . . . . . . . 11 2.3 Investigations i n the Time-Domain .......»•• 13 2.4- Conclusions «>«>•*>.«>•« »«««»•»«>•••«> «. .•.•..»» •.« 20 3. Approximate Solutions to the System Equation .... 21 3 » 1 Mo t X V at X Q.Xl • » b * > * A « o e o e « o o o o o o o o o c * o o o o o o o o o t > » 2 1 3 » 2 C c X S G I ™ £ 1 a « * o o « o o o o o * o o o o o * o o o o o o * o o o » * 2 3 3 « 2 . 1 Choice of Approximant » 0 * . > o o o o o o & o * * * 2 3 3 . 2 o 2 The Angle. C r i t e r i o n and the Determine atiOn 0f t / - . . » • *•-. . o o o o o o o o o e o o o o o o o . 2*4 m 3.2.3 Determination of A(t) and ••••• 36 3.2.4 Determination of P and 0 q ... 44 3.2.5 Exampl© A . O . O O O O O O O O O O O O ' O . O O O O O O O O O O . . 45 3.2.6 Refinements i n the Approximation .... 49 3.2.7 Summary and Examples of the Refined Parabolic Phase Approximation . . . < , . 5 4 3.2.8 Errors and Limitations ......<>.<>...•• 60 3.3 Case II £ 1 * . 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 . 0 65 3.3.1 Choice of Approximant . . . . o . 6 5 3.3.2 Determination of n . . . . . 0 . 0 . 0 0 . 0 . . O . . 68 i i i Page 3.3.3 Determination of g ..««».....oa......*' 73 3.3.4 Summary and Example of the Correction Term Approximation .................. 80 3.3.5 Errors and Limitations .............. 83 3.4 Summary «ft* + .....»«o . .oo9««o».o.»..fte.« 87 4. Co nclUSion *^.»..ooooe»oooo* « o 9 o o.oo«..o»«* 92 Appendix A On Computation •.......»••........•...<>•• 94 Appendix B The Kry l e f f and Bogoliuboff Approxi-ITl3 /~bX0Xl » « a « o « * o o o o o o « o o « o « o o « « o e o o « o e o « o « 95 Appendix C A Measure of Closeness between Linear and Non-linear Isoclines • i > e . l > o . i e o . . . o a . . 98 References «»oooo»o*«*oo40o.o.oooooo.oao.ooo.oo.o».«« 102 iv LIST OF ILLUSTRATIONS F i g u r e • Page 1»1 N o n - l i n e a r c h a r a c t e r i s t i c s ao©(&o»soooooooo«» 4 l e 2 Non=l inear RLC c i r 103 R e s t r a i n e d s l i d i n g mass a o o o o o o e o o o B o o o o a n a 6 1«4 T o r s i o n a l pendulum »ooooooo*e.»o©«>oo 2©X Phase-^plane diagrams a o o « o © e * 8 o o » & o o o « o o 8 o « » 14 2«2 Dependence of overshoot on s and gx « 0 oooooa 15 2«3 S o l u t i o n curves to the equa t i on x + 0 a4x * x + J A X 3 = 0* x(0) = lg x(0) — 0 o o o o o o o o o o B o u o o o e a o o o o e o e o a 17 2©4 S o l u t i o n curves to the equa t i on x + 2 e4x + x • f*x3 = O s x(0) = 1; x(0) — 0 oosooooeooooe oeoooeoooeoee 19 3»3L D i f f e r e n c e between the l i n e a r and non -l i n e a r i s o c l i n e s of the same slope s m »««»« 26 3»2 E f f e c t of the n o n ^ l i n e a r term o»eooo«e<><>o<><» 27 3»3 E f f e c t of the n o n - l i n e a r term of ampl i tude 32 3o4 Dependence of t Q on A ( t m ) * i» e«o«ooooe<o«oo»» 32 3<>5 Dependence of t on JX , 34 3Q6 Dependence of t on e j>oooooe«oo»o»oooooo«.o 35 o 3o7 Comparison between approx imat ing methods s » 46 3«8 Comparison of the t rue f requency and phase and t h e i r approx imat ions ©ooo©/aoo«»oooo«>oo» 51 3»9 Approx imat ion by the r e f i n e d p a r a b o l i c phase method f o ? e = 0»4p |A=3 g o i n i » « « « > i t « i > i i > i t i < o 58 3M10 Approx imat ion by the r e f i n e d p a r a b o l i c phase method f o r equa t i on - + x + 5 x + 1 0 x 3 = 0 . 61 3 a l l - C o r r e c t i o n term z ( t ) f o r c = l » 2 and u,=2 a*.** 67 3©12 Contours of cons tant n from exper imenta l r 6 SUl t S a o « O 3 C i a o G £ i o a t > 0 f l b i i i U 9 « a v « * ' i ' « » a c o « u a « & 70 3«*13 Approx imat ion to the contours of constant 71 Page 3 < > 1 4 Determination of c as function ofna9e»»0o 7 2 3 « 1 5 Comparison between z(t) and z(t) for e = l » 2 and U/ — 2 o o o o o o o o o o o o o o o o o o o o o e o o o e o o o o o o o o * 7 4 3 * 1 6 Contours of constant g •« o • • • • • • • « « < > « o o o o • •> 7 6 3 « 1 7 g as a function of a for constant e <>« <>•»<>• 7 7 3 « 1 8 l ° § i o ^ v s " * " ° ^ 1 0 ' x " ^ o r c o n s " * ' a 1 1 ^ e 7 8 3 . 1 9 Approximation !to F i g . 3 , 1 8 • . • • • • • . . . . . . < > • • 7 9 3 o 2 0 Approximation/by the correction term method for e=1.4, u= 3 . . • • • • • • • • » * o o . o « « » » » 8 4 3 » 2 1 Magnitude of the deviation for e = l , 4 , u= 8 • 8 6 C*l A measure of closeness between l i n e a r and non-linear isoc l i n e s of the same slope m «• 9 9 v i APPROXIMATIONS TO THE FREE RESPONSE OF A DAMPED NON-LINEAR SYSTEM lo INTRODUCTION 1.1 Mathematical Models In the study or analysis of physical systems, i t i s common practice to represent them by mathematical models. In order to make the mathematical models more tractable, certain simplifying assumptions are usually made„ For example, the d i f f e r e n t i a l equations evolved are usually l i n e a r i z e d , so that they may be solved by well established techniques used for linear d i f f e r e n t i a l equations 0 Most physical systems, however, behave i n a manner which i s far from l i n e a r , for example, a triode amplifier with large signal inputs or a mass, restrained by a non-linear spring, o s c i l l a t i n g with large amplitudes. Therefore, non-linear analysis i s required in order to y i e l d results closer to r e a l i t y <> 1o2 A n a l y t i c a l Approximations Exact solutions to non-linear d i f f e r e n t i a l equations, are usually d i f f i c u l t y i f not impossible, to f i n d i n closed form 0 Techniques for solving the equations vary according to the types of-equations involved, and are very l imitedo Y i t h the aid of d i g i t a l computers, numerical solutions, to almost any degree of a c c u r a c y t o any non-linear d i f f e r e n -t i a l equations are available» Using analogue computers, solu-tions to ordinary d i f f e r e n t i a l equations can be obtained. The solutions obtained from these computers, however, do not furnish a l l the information concerning the physical system of interest, i . e . they usually reveal only the behaviour of the system under certain p a r t i c u l a r conditions. Exhaustive tests are needed i f some insight into the system i s required, and an engineer cannot necessarily predict from them how the system w i l l behave i f some parameters i n the system are changed. The problem of cost and a c c e s s i b i l i t y i s another disadvantage i n using computers as a means to solve a non—linear d i f f e r e n t i a l equation. For these reasons, an a l y t i c a l approximations to the solutions of ordinary non-linear d i f f e r e n t i a l equations are developed. These approxi-mations are obtained i n algebraic or transcendental form without the necessity of introducing numerical values for parameters or i n i t i a l conditions during the process. Though some degree of accuracy i s s a c r i f i c e d , an over-all insight into the system i s often obtained at a low cost. For instance, the dependence of the solution on a certain parameter may be e x p l i c i t , thus y i e l d i n g useful information for system design. A few well established approximate ana l y t i c a l m e t h o d s a r e (a) Perturbation method, (b) Variation of parameters, (c) Averaging methods based on residuals, and (d) Principle of harmonic balance. Though these methods are developed to cover a very large class of non-linear d i f f e r e n t i a l equations, they have a common weakness in that they are incapable of dealing with equations exhibiting gross n o n - l i n e a r i t i e s . This l i m i t a t i o n i s due to the general approach to solving the equation, namely, making the grossly non-linear equation only s l i g h t l y non-linear, or quasi-linear, i n an attempt to get more insight into the behaviour of the system, using linear theory. In order to break through this l i m i t a t i o n a bolder approach i s i n order, i . e . a direct attack on the non-linear equation i n question. To this end, a study, for the purpose of obtaining approximate solutions to a certain type of grossly'non-linear equation which arises from many engineering systems, was undertaken* 1»3 Derivation of the System Equation A large class of physical systems contain a non-linear element whose characteristic i s represented by an odd cubic polynomial with positive c o e f f i c i e n t s , for example, a hard spring characterized by F(x) = a^x + a^x where F(x) ^ restoring force i n spring, x = displacement, and a^ and a-j are positive c o e f f i c i e n t s . This odd cubic polynomial i s often an approximation to a grosser non-linearity such as an odd polynomial of higher order, i 0e„, k = n F(x) = ^ a kx k, n > 3. k odd The general shapes of th i s odd polynomial and the "odd cubic" characteristic are shown i n F i g . 1.1. Now consider the following systems containing "odd cubic elements: 4 non-linear inductor F i g . 1.2 Non-linear RLC c i r c u i t (a) RLC c i r c u i t . The p a r a l l e l RLC c i r c u i t i n Pig. 1.2 consists of a line a r r e s i s t o r R, a lin e a r capacitor C, and an inductor which i s non-linear because of saturation. Neglecting hysteresis, the inductor current can be approximated by 3 = a-^ A + a-^ A where A = flux linkage, and a-^ and a^ are positive c o e f f i c i e n t s . Applying Kirchhoff's current law, we obtain Cv + Jj + & X + a 3X 3 = 0 * where v i s the voltage across the p a r a l l e l elements. Since, by Faraday's law, v = X, this equation can be re-written as CA + | A + a ] A + a 3A 3 = 0 ( l . l ) which i s , then, the equation of the RLC c i r c u i t . (b) Hard spring with pure viscous damping. F i g . 1.3(a) shows a simple mechanical system involving a mass s l i d i n g on a surface with pure viscous f r i c t i o n and res-trained by a hard spring whose characteristic i s given by F = b-j^ x + b 3 x 3 where F = restoring force i n spring, and b-^ and b^ are positive c o e f f i c i e n t s . From the free-body diagram shown i n F i g . 1.3(b), • A dv V = dt 6 \ non-linear spring (a) Bx (a) Pig. 1.3 Restrained s l i d i n g mass D I non—linear torsion bar I c-Pig. 1.4 Torsional pendulum 7 we obtain F + Bx = -Mx where M = mass. and B = c o e f f i c i e n t of viscous f r i c t i o n , and hence Mx + Bx + b ^ + b 3 x 3 = 0. (1.2) (c) Torsional pendulum with pure viscous damping, A simple torsional pendulum i s i l l u s t r a t e d i n F i g . 1.4. It consists of a disc with moment of i n e r t i a J, a support with viscous f r i c t i o n a l c o e f f i c i e n t D, and a non-linear torsion bar whose characteristic i s given by T = c 1Q + c 3 © 3 where T = restoring torque i n the torsion bar, 0 = angular deflection, and c^ and c 3 are positive c o e f f i c i e n t s . Consideration of 7 torques gives T + DO = ^J0 and the system equation becomes JO + DO + CjO + c 3 © 3 = 0. (1.3) The systems described above are, i n f a c t , analogous to each other, because th e i r equations a l l have the following form ay + Py + yy + S y 3 = o (1.4) where a corresponds to C, M, or J, P corresponds to jj^ , B, or D, "y corresponds to a-^ , b^, or c-^ , 3 corresponds to , b-j, or c^, and y corresponds to 0, x, or Oo Thus, a s o l u t i o n to equation (1.4) w i l l provide a s o l u t i o n to a l l the above systems. I f both the terms Sy and § y ^ are r e l a -t i v e l y s m a l l , a c l a s s i c a l method based on v a r i a t i o n of para-meters, such as the K r y l o f f - B o g o l i u b o f f method, gives s a t i s -(2) (3) f a c t o r y a n a l y t i c a l approximate s o l u t i o n s . However, i f e i t h e r py or §y^ i s l a r g e , t h i s method f a i l s to y i e l d good r e s u l t s . A d e t a i l e d study of the equation where |3y and Sy"^ are not n e g l i g i b l e and a d i r e c t approach to f i n d i n g a n a l y t i c a l approximate s o l u t i o n s was t h e r e f o r e attempted. 2. STUDY OF THE SYSTEM EQUATION 9 2.1 Normalization Equation (1.4) ostensibly contains four arbitrary c o e f f i c i e n t s , which make i t d i f f i c u l t to study. However, two of these c o e f f i c i e n t s can be made i m p l i c i t i f the following normalization i s performed. Dividing through by a, equation (1.4) can be rewritten as y + ay + by + cy 3 = 0 (2.1) where a = 6/oc, b = 'Y/a., and c = c5/cx.. Letting T = Jb t, obtain y _ dt _ dy_ dT ~ dT ' dt = fb y' we where y' -dT' 2 , d y and y = — i r dt* dT v d t ; = Vb ^ ( VF y ) = b y " , ' d 2 where y" = — ^ . Substituting y' and y" into equation (2.1), ar 2 we have by" + a Vb y' + by + cy =0, or y " + ^ y ' + y + f y = o . a „, , , c „3 10 In order to f a c i l i t a t e subsequent work, this equation i s re-written i n the form y » + 2ey» + y + ay 3 = 0 (2.2) where e = a _ , 2 fb and u = ^ « Moreover, since a consistent framework i s desirable, t h i s equation i s assumed, to have the following i n i t i a l conditions g y(0) 1 P y'(0) = 0. To show that this assumption does not affect the generality of the approach, l e t the i n i t i a l conditions be y(0) = Q, y»(0) = 0, Replacing y by x = ^— y i e l d s *^o x(0) := 1, x' (0) = 0, and equation (2«2) becomes Q Qx» + 2eQ0x» + Q Qx + U.Qq3 x 3 = 0 y' (0) i s assumed to be zero i n this work for physica.l reasons. For examplep i n the study of the s l i d i n g mass referred to i n 1.3, one usually displaces i t from i t s neutral position by the i n i t i a l amount of Q O ? then releases i t . Very seldom does one incorporate an i n i t i a l v e l o c i t y because i t i s d i f f i c u l t to obtain accurately. Moreover, i n the case where the system i s o s c i l l a t o r y , one can a r b i t r a r i l y f i x t ~ 0 at the peak of the o s c i l l a t i o n j , i . e . y f(0) — 0. 11 2 3 or x" + 2ex' + x + uQQ x = 0, x(0) = 1, x'(0) = 0 (2.3) A comparison between equations (2.2) and (2,3) reveals that the substitution of x = ¥j— leads to an equation with a d i f — *^o ferent c o e f f i c i e n t i n the non-linear term i f Q ^ 1; but since o the method to be used for solving the equation i s not altered by the values of this c o e f f i c i e n t , no generality i s l o s t . 2.2 Phase-plane Analysis Consider the equation 3 x + 2ck + x + ux =0, x(0) = 1, x(0) = 0. (2.4) If u = 0, i t degenerates to the linear equation x + 2ex + x = 0, x(0) = 1, x(0) = 0. (2.5) which w i l l be referred to as the complementary linear equation of (2.4), Exact solutions of this equatioti depend on the values * of e , namely, i f (a) e< 1 (underdamped) f then X = 1 e ~ e t cos ( Jl - e 2 t + 0 )F 2 0 Jl - e' (2.6) Because the system to be considered contains only passive elements, i . e . no energy sources, the value of ,e w i l l either be zero or p o s i t i v e . where 0 = tan -1 t e 2 (b) e = 1 ( c r i t i c a l l y damped), then x = (1 + t) e t , (2.7) (c) e >1 (over-damped), then 1 (! + e ) e ( - E + Ve 2- l ) t x te" - 1 2 ' j ! * + | (1 + ~ = ) e ^ 6 ^ " ^ ' * . (2.8) If [i ^ 0, solutions are also dependent on the damping factor e For example, i f (a)' e = 0 (conservative system), then equation (2.4) becomes x + x + ux = 0, x(0) = 1, x(0) = 0. The solution i s the Jacobian e l l i p t i c cosine x = Cn (k,<ot)(4) where k^ = ^ " • 2(l+u) 2 and (0 = 1 + [j,« (b) e >0 (non-conservative system), then no exact solutions to equation (2.4) are available i n closed form. Though the solutions are similar to those of equation (2.5), they exhibit a greater number of o s c i l l a t i o n s i n the same length of time. In order to c l a r i f y the picture, a phase-plane analysis i s most 13 (5) h e l p f u l . Using the method of isocl i n e s , phase—plane diagrams of equations (2.4) and (2.5), as shown i n Pig. 2.1, are obtained. In F i g . 2.1(a) where e < l , the two tr a j e c t o r i e s suggest damped o s c i -l l a t i o n s . In the case where u>0, the period of o s c i l l a t i o n i s shorter because x decreases with a greater slope. In F i g . 2.1(b) where e=l, the trajectory for u=0 represents a solution without "overshoot", i . e . x never going negative, while tr a j e c t o r i e s for u>0 show one or more overshoots, the number of which increases as u increases. Similar results are observed i n F i g . 2.1(c) where e > l . If under—damping and over-damping are defined respectively by the presence and absence of overshoots, one sees, therefore, that consideration of e alone i s not s u f f i c i e n t to predict whether the system i s under-damped or over-damped, as i n the l i n e a r case, for the value of a i s also an important factor. An extensive i n v e s t i -gation of t h i s aspect was undertaken, using the d i g i t a l computer to provide numerical solutions. The re s u l t , as i l l u s t r a t e d i n F i g . 2.2, i s a curve i n the e-u plane showing regions where the system has overshoot and where i t has not. Contrary to li n e a r theory, over-shoots may be observed for e > l , i f a is high enough. This result i s not surprising as equivalent l i n e a r i z a t i o n also predicts possible overshoots. 2.3 Investigations i n the Time-domain The independent variable, time, i s i m p l i c i t i n phase-plane diagrams, and solutions to equation (2.4) are, therefore, not readily available as functions of time. Using the method of equi-(6) valent l i n e a r i z a t i o n and numerical solutions obtained from the d i g i t a l computer, the effect of u on the solution i s as followss _ , . _ _ See Appendix A for computational d e t a i l s . (a) 0 < e < l (b) e = 1 (c) e > l u=--u3>u2 Pig. 2.1 Phase-plane diagrams 15 F i g . 2.2 Dependence of overshoot on e and u 16 (a) 0 < e < 1 Pig. 2.3 shows a t y p i c a l example i n this case. A numerical solution to the equation x + 0.4x + x + 2x 3 = 0, x(0) = 1, x(0) = 0 (2.9) i s displayed together with the solution to i t s complementary linear equation, i . e . x + 0.4x + x = 0, x(0) = 1, x(0) = 0 (2.10) Any difference between these two curves r e f l e c t s the effect of the non-linear term. The solution to equation (2.10) i s x =' — ± — e ~ ° * 2 t cos (7o796t - 0 ) yo~96 where 0 = tan * • It represents a damped sinusoid having 0 V0.96 an envelope — g-0.2t a n < ^ a phase of JO .96 t - 0 . Although JO.96 0 the solution to equation (2.9) resembles a damped sinusoid, i t s amplitude decays at a slower rate than — 1 e -^' 2^ and i t s 70.96 phase increases i n a non—linear manner, i . e . the phase i s retarded as time increases. Equivalent l i n e a r i z a t i o n of equation (2,4), based on (6) v a r i a t i o n of parameters , yields , .2 2 x + 2ex + (1 + x = 0 (2.11) -A where A = — e = amplitude. Here, the value of A at t = 0 must not be taken as x(0), for i t Envelope for u,=2 Envelope for \i=0 A 1 -0.2t A = — = e 1 -0.2t cos(,/0.96t -0.2) F i g . 2,3 Solution curves to the equation x + 0.4x + x + [xx3 = 0, x(0) = 1, x(0) = 0 18 represents the envelope of the solution and i s always greater than x(0) , though i n most cases the difference i s small. Using linear theory, the effect of u, on the frequency i s evident. This equation represents an o s c i l l a t i o n with both amplitude and frequency varying with respect to time. The amplitude A decays —e t exponentially according to e , and the frequency s i m i l a r l y decreases as the amplitude decays with increasing time. As time progresses, A ultimately becomes small enough so that 3uA 2 — ^ — i s negligible compared to unity. Then, the frequency I 2 becomes e f f e c t i v e l y Jl - e , which i s the frequency of o s c i l l -ation of the complementary linear equation of equation (2.4). Therefore, equation (2.4) degenerates to i t s complementary linear equation as time increases. (b) e > 1 A t y p i c a l example i n this case is i l l u s t r a t e d i n F i g . 2.4 by x + 2.4x + x + 7x 3 = 0, x(0) = 1, x(0) = 0. Although this equation has an over-damped complementary linear equation, an overshoot i s observed i n the solution. The presence of the non—linear term i s responsible for this over-shoot as already shown i n the phase-plane analysis. Consistent results are also predicted from equivalent l i n e a r i z a t i o n . Consider now the equivalent linear equation (2.11) where e ^ 1 . If u, has a value such that 1 + > e , o s c i l l a t o r y solutions may be obtained. It must be noted, however, F i g . 2.4 Solution curves to the equation x + 2.4x + x + u.x3 = 0, x(0) = 1, x(0) = 0 20 that this inequality only predicts q u a l i t a t i v e l y how overshoots may possibly occur, and i s not necessarily capable of y i e l d i n g accurate r e s u l t s , for small values of e and a have been assumed in the equivalent l i n e a r i z a t i o n . If accurate results are required, then the curve i n Pig. 2.2 may be used. 2.4 Conclusions In conclusion, the general behaviour of the system has been studied, using the d i g i t a l computer and the method of phase plane analysis* Solutions to the non-linear equation are compared with solutions to i t s complementary linear equation. In damped o s c i l l a t o r y systems, the amplitude of the solution decreases more slowly i n the non-linear case, and i t s frequency, being i n i t i a l l y greater, approaches that of the complementary linear equation as time progresses. In the case where the complementary linear equation i s c r i t i c a l l y or over-damped, the presence of the non-linear term may lead to over-shoots . Because the system i s damped, x w i l l eventually disappear, 3 and the non-linear term ux w i l l become i n s i g n i f i c a n t compared to x, when x becomes small. Hence, an approach to finding the a n a l y t i c a l solution i s suggested by neglecting the non-linear term at a point where x has become s u f f i c i e n t l y small. 3. APPROXIMATE SOLUTIONS TO THE SYSTEM EQUATION 21 3 .1 Motivation Consider equation (2.4) x + 2ex + x + ux 3 =0, x(0) = 1, x(0) = 0. As mentioned previously, the solution to this equation i s not available i n closed form. Therefore, attempts were made to approximate the solution. In the case where e < 1, the method of Kryloff and *(3) Bogoliuboff was used. From Appendix B, the solution i s given by x(t) = e " e t cos (1 + 2|)t, (3.1) where e, a ^ 1. By comparison with x(t) = , 1 2 e " e t cos ( Jl - e 2 t + 0 ) (3.2) VI - e 0 which i s the solution to the complementary linear equation (2.4), one observes that solution (3.1) cannot be extended to higher values of e and u t because (a) contrary to the resu l t i n the l a s t chapter, the frequency in (3.1) remains constant as time progresses, and (b) this frequency approaches unity i f a becomes zero, which does not agree with (3.2) i f e i s not ne g l i g i b l e . For large values of e and a, therefore, a new method must be _ This method w i l l hereafter be referred to as the K-B method. See Appendix B. 22 developed. In the case where e ^ 1 , the K-B method i s no longer applicable because the solution i s not o s c i l l a t o r y i n nature. (7) Here, the Ritz method , which i s an averaging method based on residuals i s used i n conjunction with i n i t i a l condition matching. Since the solution of interest i s either monotonically decreasing or exhibits one overshoot, the approximate solution i s assumed to be of the following form: A at , -n bt x(t) = A e + B e where A and B are constants and both a and b are dif f e r e n t negative numbers. Substituting this solution into the o r i g i n a l d i f f e r e n t i a l equation, we have the residual given by 0~(t) = x -f 2ek + x + u,x3 = (a 2+ 2ea+l)Ae a t + (b 2+ 2eb+l)Be b t+ u A 3 e 3 a t _L T J 3 3bt _ A 2 B (2a+b)t A T a2 (a+2b)t + \iB e + 3uA Be v ' + 3uAB e v ' . The Ritz c r i t e r i a are 'CD C T ( t ) e a t dt = 0, and / ( X ( t ) e b t dt = 0. 0 Now, from the i n i t i a l conditions x(0) = 1, x(0) = 0, we also have A + B = 1 a A + b B = 0 Hence we obtain four equations in four unknowns. After integrating and eliminating A and B from these equations, we have 23 3 1. i^(a 2+ 2ea+l) + -%r(h2+ 2eb+l) - — ^ 2 a v T - a + b ^ T — / 2 ,4a(a-b) 3 2 2 ixa + 3ub a 3aa b (3b+a)(a-b) 2 (3a+b)(a-b) 2 2(a+b)(a-b) 2 3 2o =|(b2+ 2eb+l) + - T u(a 2+ 2ea+l) - — ^ = 0 2 b r a + b „ 4 b ( a - b ) 2 3 2 2 ub + 3aa b _ 3aab (3a+b)(a-b) 2 (3b+a)(a-b) 2 2(a+b)(a-b) 2 0 . Without the aid of a d i g i t a l computer, solving the above two equations simultaneously i s very laborious. Therefore, from an engineer's point of view, this method i s highly impractical. As a r e s u l t , two new approximating methods were developed, depending on whether e i s less than unity or greater than unity. The rest of this work w i l l be devoted to the development of these methods. 3.2, Case I - c < 1 3.2.1 Choice of Approximant In this case, where e < 1, both phase-plane analysis and computer solutions from Sections (2.2) and (2.3) have indicated o s c i l l a t o r y solutions resembling damped sinusoids. The phase has also been shown to increase non-linearly with time. The approximant, therefore, w i l l assume the form x(t) = A(t) cos H ( t ) , where A(t) = amplitude, and .Q(t) = phase. Consider equation (2*4) with e < 1. Because the system i s damped, x w i l l eventually vanish and the non-linear term 3 u,x w i l l become negligible compared to x, when x i s suf-f i c i e n t l y small. Let this point of n e g l i g i b i l i t y occur at t = t m < Beyond this point, then, equation (2.4) es s e n t i a l l y degenerates to i t s complementary linear equation (2.5) and w i l l be treated as such. Therefore, the approximant w i l l have the following form: for 0 < t ^ t m x(t) = A(t) cos f l U ) (3.3) for O t x(t) = P e ~ e + cos ( V l - e 2 t + 0 ) ^ m v ' v o (3.4) where P and 0 are constants, o IJere, i t must be noted that P ^ . s=- and 0 ^ tan ^ e e as they are i n equation (2.6), because i n i t i a l conditions for (3.3) must be adjusted to match (3.2) at t m . And so, the problem i s now to f i n d the functions A(t) and .Q(t), and the constants P, 0 , and t . ' o' m 3.2.2 The Angle C r i t e r i o n and the Determination of t ;_m In order to determine t , some information about the m' 3 point at which [xx may be neglected i s necessary. Consider, again, the phase-plane diagram. Usually the f i r s t step i n the construction of the phase-plane diagram i s the construction of i s o c l i n e s , i . e . curves of constant slope i n the x-x. plane. Therefore, i f two systems have almost i d e n t i c a l sets of 25 i s o c l i n e s , they must have almost the same phase-trajectories. Furthermore, i f two systems have almost i d e n t i c a l phase-t r a j e c t o r i e s , i t i s reasonable to assume that their solutions as functions of time are almost i d e n t i c a l . Hence a measure of "closeness" between the i s o c l i n e s of two systems may be regarded as a measure of how close their solutions are to each other. Now, consider the two systems represented by the equations (2.4) and (2.5)« Typical i s o c l i n e s of the same slope m for these systems are shown in F i g . 3.1. In order to have a measure of "closeness" between these i s o c l i n e s , a c i r c l e of radius R i s constructed, intersecting the linear i s o c l i n e at point P. Through P, a v e r t i c a l straight line i s drawn, i n t e r -secting the non-linear i s o c l i n e at point Q. Then, the angle So between the lines OP and OQ can be regarded as a measure of "closeness" between the two i s o c l i n e s . If a slope d i f f e r e n t from m i s chosen, and the same construction performed, the resulting angle § 0 may be d i f f e r e n t . However, the maximum value ( So) of these angles, as m varies, i s a function of ^ ITlcLX aR . In p a r t i c u l a r , ( So) decreases as uR decreases, but m ax since the amplitude A of the solution follows R quite closely, (So) decreases as uA decreases. Therefore, the value of 2 uA can also be regarded as a measure of "closeness" between the two i s o c l i n e s , or a measure of the effect of the non-linear term. For example, i f uA = 0.2, ( So) i s approximately 3 , or 0.05 radian. Therefore, as the amplitude A decays from i t s 2 i n i t i a l value and reaches a value such that uA =0.2, the two sets of is o c l i n e s may be considered coincident for a l l p r a c t i c a l See Appendix C. 2 6 x=y F i g . 3.1 Difference between the l i n e a r and non-linear is o c l i n e s of the same slope, m 27 purposes . T h i s p o i n t then determines a p o i n t at which the term 3 |xx may be n e g l e c t e d . Hence t may be chosen such tha t I t must be remembered t h a t t h i s r e l a t i o n i s an a r b i t r a r y c r i t e r i o n based on the c o n s i d e r a t i o n of the angle ( § © ) , and ° max i t serves on l y the purpose of o b t a i n i n g a p o i n t where the non-l i n e a r equa t i on can be r e p l a c e d by the l i n e a r one. In choos ing 2 such a p o i n t of t r a n s i s t i o n , how the va lue of uA v a r i e s w i t h time must, a l s o be c o n s i d e r e d . Maximum - t i m e , t P i g . 3.2 E f f e c t of the N o n - l i n e a r Term F i g . 3.2 shows the genera l shape of the va lue of uA as t ime v a r i e s . S ince the ampl i tude A decays almost e x p o n e n t i a l l y , 2 as mentioned b e f o r e , uA drops v e r y q u i c k l y a t the beg inn ing and approaches zero a s y m p t o t i c a l l y as t ime i n c r e a s e s . T h i s 28 curve also represents the effect of the non^linear term. In 2 order to show that consideration of the value of u,A alone may not necessarily lead to a wise choice of the point of t r a n s i t i o n from the non-linear equation to a l i n e a r one, the following two choices of such a point are comparedo F i r s t , 2 . l e t the t r a n s i t i o n occur at t = t , , The value of ixA i s m ml r small as shown in F i g . 3*2, indicating that the equation i s e s s e n t i a l l y linear for t = t = t and therefore the " l i n e a r " m ml part of the approximant, i . e . equation (3.4): for t ^ t m x(t) = P e ~ e t cos ( J l ~ e 2 t + 0Q) 9 i s very close to the exact solution. As a second choice l e t b = t 0, m m2 the t r a n s i t i o n occur much e a r l i e r , at t_ _ o 0 The value of 2 uA i s now larger, and the " l i n e a r " part of the approximant i s therefore not as good as the f i r s t choice„ This does not necessarily mean, however, that the second choice i s poorer, for i t may y i e l d a better "non-linear" part of the approximant, i.e , equation (3.3): for 0 < t ^ t m x(t) = A(t) Cos O ( t ) . In f a c t , a better "non-linear" part i s usually expected because i t s range of approximation i s now greatly reduced. As a r e s u l t , one must consider a compromise i n choosing the point t , so that both the " l i n e a r " and "non-linear" parts of the approximant are reasonably accurate. To this end, one must * i 2 choose t as small as possible and at the same time, avoid an m unduly large value of ufA(t )1 . Using Galerkin's method, which (8) i s an averaging method based on residuals ', i t seems that the 29 value of t may be optimized i n the sense that the integral <T 2(t) dt i s a minimum, where C7"(t) = x + 2ex + x + ax. However, this i s impractical because (a) the functions A(t) and .Q(t) are not known, and (b) even i f they are known, solving the set of equations r o-2(t) dt = o cr 2(t) dt = o CO c r 2(t) dt = o w i l l be a formidable task due to the presence of the non-linear term. Therefore, experimental results are used to obtain an empirical c r i t e r i o n for choosing an acceptable t . For example, numerical solutions of the equation x + 0.4k + x + ux 3 = 0, x(0) = 1, x(0) = 0 with various values of a have indicated that the phase-retard-2 ation becomes negligible when uA = 0.2. This means that the 2 equation behaves p r a c t i c a l l y l i k e a linear one when aA = 0.2, and therefore t can be chosen such that m u[A(t )] = 0.2. (3.5) From this empirical c r i t e r i o n , as u becomes larger, A ( t m ) becomes smaller, y i e l d i n g a greater t m« This i s reasonable 3 because with a larger a, the non-linear term ux must take a longer time to become negligible compared to x. It must also be noted here that the r e l a t i o n (3.5) i s obtained empirically with e = 0.2. If equations with a larger e are considered, however, a di f f e r e n t empirical c r i t e r i o n may be obtained. In fac t , study of the numerical solution of the equation x + 0.8 k + x + ax 3 = 0, x(0) = 1, x(0) = 0 has revealed that the phase-retardation becomes negligible at tm S^en by From F i g . 3.2, this indicates that the effect of the non-linear term i s greater at t m , but i t should be noted that a larger e results i n a faster decay of the amplitude A. Because the 2 ef f e c t of the non—linear term varies as A , this means that 3 ux becomes negligible compared to x i n a much shorter time i n t e r v a l . Therefore i t i s conceivable to relax the c r i t e r i o n . Many examples with various values of e have been solved numerically and the result has suggested that this c r i t e r i o n can be assumed to depend on e i n the following l i n e a r manner. a [ A ( t m ) ] 2 = | . (3.6) Note that this r e l a t i o n i s only an empirical c r i t e r i o n to be used as a "rule of thumb" for choosing t wisely, and i s not necessarily the best c r i t e r i o n , i f one exists at a l l . Now, the next step i s to evaluate t , using this c r i t e r i o n . 31 Although this c r i t e r i o n gives a value of A ( t m ) when e and u are specified, i t does not provide the value of t d i r e c t l y . It i s necessary to f i n d a relationship between A(t ) and t • Again, consider equation (2.4) x + 2ex + x + ax 3 = 0, x(0) = 1» x(0) = 0. It has been shown e a r l i e r i n section 2.3(a) that the envelope A(t) of the solution decreases at a slower rate with a ^ 0 than with u = 0, as i l l u s t r a t e d i n F i g . 3.3 by curves I and I I . A horizontal l i n e i s drawn through A(t) = A ( t m ) , i n t e r -secting curves I and II at F and G respectively. Therefore, the abscissa for G i s t , and i f the abscissa for F i s denoted by t , we have J mo' , -e t . / , \ 1 mo which gives Letting t be the i n t e r v a l between F and G, then t = t + t m mo o "7" log [/l - e 2 A(t m)] + t o . (3.7) Hence the problem becomes finding t i n terms of A ( t m ) , u and e, which are a l l the known quantities. Here, i t may seem that the introduction of t does not help i n solving the problem at a l l , because the o r i g i n a l problem was to find t also i n terms of these three known quantities. But this i s not the case, because by finding t , one i s looking only for that part m I. uj^ O I I . u=0 mo Pig. 3.3 Envelope of the solution to equation (2.4) i x + 2ex + x + fxx3 = 0 Fig« 3.4 Dependence of t on A ( t m ) of t due to the non-linear term, i . e . only the effect of the m non-linear term on the envelope. ' The approach i n finding t i s to investigate how t depends on the three quantities A ( t m ) , u, and e respectively. Here, experimental results are again used. F i r s t , the dependence of t on A ( t m ) i s i l l u s t r a t e d by F i g . 3.4, i n which the numerical solutions to the equations x + 0.4 x + x = 0, x(0) = 1, x(0) = 0 and x + 0*4 x + x + 3x 3 = 0, x(0) = 1, x(0) = 0 are used. As A(t ) decreases from i t s i n i t i a l value, t m 7 o increases f a i r l y l i n e a r l y and reaches a maximum at A ( t m ) = 0. Hence we have t QcX 1 - A ( t m ) (3.8) Secondly, i n order to reveal the dependence of t on u,, the numerical solutions to the equation x + 0.4 x + x + u.x3 = 0, x(0) = 1, x(0) = 0 are shown i n F i g * 3.5(a) as u. varies from 0 to 3 i n increments of 1.0. A fixed value of A(t ) i s chosen, and values of t , m corresponding to particular values of u, are found. If t i s now plotted against u,, a straight l i n e i s obtained, as i n F i g . 3.5(b)* Thus, t varies approximately l i n e a r l y as u,, or o t 0CC [x (3.9) The i n i t i a l values of the envelopes are assumed to be unity here for the purpose of finding an approximation. Their true values, however, are greater than unity. 34 A(t) (b) F i g . 3.5 Dependence of t on |x 35 Now, the dependence of t on e can be seen from the numerical solutions to the equation x + 2ex + x + 3x 3 = 0, x(0) = 1, x(0) = 0 as e varies frpm 0 to 0.9. Corresponding values of t Q , as e increases, are obtained with A ( t m ) fixed, for example, at 0,4, Then t i s plotted against e as i n F i g . 3*6. The curve obtained resembles a hyperbola, suggesting that t varies inversely as e , or t cc -o e (3.10) Therefore, from equations (3.8), (3.9) and (3.10), we have a[ l - A(t m)] t = k o where k = constant,. From numerous examples with various values o 2.5 2.0--1.5-1.0--.5-0 a=3.0 A(t m)=0.4 .1 .2 .3 .4 .5 .6 .7 .8 Fig* 3*6 Dependence of t on E 36 of e, a and A ( t m ) , k i s empirically found to be about l / l O . Ve then arrive at a [ l - A(t )1 t = LL -J2-=L , ( 3 . 1 1 ) 0 10 e which enables us to calculate ± n when E , a, and A ( t m ) are known. Hence, from equation (3.7), we have * m - - ^ o 8 e [ y T T ^ A ( t + ( , . i a ) Note that t i s always positive because both e and A(t ) are less than unity« In short, t can be calculated from equations (3«6) and (3.1l) when e and a are specified. Equation (3.6) i s an arbitrary c r i t e r i o n based on the consideration of the angle (So) between the linear and non-linear is o c l i n e s i n the max phase-plane, and equation (3.1l) i s obtained empirically using numerical solutions obtained on the d i g i t a l computer* It must also be remembered that the t thus calculated i s not necessarily m an optimal choice of the point of t r a n s i t i o n from the non-linear equation to the linear one, but rather, i s a judicious choice of such a point for the purpose of approximating the exact solution. 3.2.3 Determination of A(t) and Q(t) In the determination of the amplitude function A ( t ) , consider f i r s t the amplitude of the solution to the comple-mentary linear equation, i . e . x + 2 E X + x = 0, x(0) = 1, x(0) = 0, which i s the case for u, solution i s given by = 0, The amplitude A(t) of the 37 A(t) = , 1 e ~ e t / l - e In the case where u, ^ 0, graphs such as F i g , 3.5(a) have shown that the amplitude A(t) also resembles an exponential (13) but decreases at a slower rate. As suggested by Tuttle, ' A(t) can assume the following form: A(t) = A Q e " ^ , where A and p are constants, o Since A(0) i s assumed to be unity as mentioned previously, A =1, and therefore o A(t) = e ~ p t . (3.13) But at the point where the non-linear term becomes negl i g i b l e , -pt t = t . Hence m A(t ) = e m m 'or p - - ( l / t m ) l o g e A ( t m ) (3.14) Note that p i s always positive for A ( t m ) less than 1* "..Having calculated t f f l and A ( t m ) from equations (3 .6) and (3.12), p i s now easily obtained, and equation (3.13) becomes A(t) = exp log A(t ) &e m t m (3.15) 38 The next step i s to determine the phase function Q ( t ) . Previous study of the equation with u. ^ 0 has revealed that the phase increases i n a non-linear manner,, or the frequency of o s c i l l a t i o n Varies with time. As a f i r s t approximation, we consider that the frequency varies l i n e a r l y with time and therefore | r [ n ( t ) ] = 2<o2t + 0 ) 1 (3.16) where to^ a n < i are constant. Integrating once j we have O(t) = <o 2t 2 + a>1t + « o (3.17) where to = constant, o To f i n d « q . »^, and to^i consider equation (3.3) which now become s x(t) = e ~ p t cos (« o + w-j^ t + " 2 t 2 ) Since x = 1 at t = 0, we have 1 = cos to o Therefore to = 0. (3.18) o (9) In. order to f i n d <o^ and to^, the method used by Soudackv ' i s * In the linear case, the frequency of o s c i l l a t i o n i s constant, and equals the f i r s t derivative of the phase w.r.t. time. Therefore, as a generalization to the non-linear case, the f i r s t time derivative of the phase i s referred to as the frequency of o s c i l l a t i o n . 39 adopted. They are obtained by considering that both the phase and the frequency i n the approximants (3.3) and (3.4) should be matched at the point of t r a n s i t i o n . This means that at t = t , m* H ( t ) = Jl - e 2 t + 0 , m * m o and 0 ( t ) = Jl - e 2 . From equations (3.16), and (3.17), we have <o + <ont + <0 ot 2'= Jl - e 2 t + 0 , (3.17a) o 1 m 2 m ^ m 'o 7 \ > / and tt + 2»2 tm = «/l - e 2 . (3.16a) From the l a s t equation, we then obtain Jl - e - » 1 w0 = A — . (3.19) 2 * m Since ttg i s now expressed e x p l i c i t l y i n terms of tt^, a l l that remains to do i s to determine tt^ independently. But before we do so, l e t us see whether the parameters p and (6^ determined are consistent with the case where e —»-0. The l i m i t i n g values for both p and (a^ a r e expected to be zero, because i f e = 0, we have (a) A(t) = 1 and (b) frequency = constant, i . e . no phase retardation. Let us f i r s t consider equations (3.6) and (3.12) as e —»-0. We have lim A(t ) = lim /rr- = 0, m' V 2u s c—0 e —0 and lim t = lim ' 0 E-*-0 1-, f f. 2 xi . - A ( tm }] - l o g |V 1 - e A ( t ) J + — lOe 40 = CO These results could also be derived from the following argument: As e becomes smaller, the envelope w i l l decrease at a slower rate and i t w i l l take longer time to reach the point of t r a n s i t i o n to the li n e a r equation, i . e . t w i l l become larger. Eventually, as z gets very close to zero, t w i l l approach i n f i n i t y and A ( t m ) w i l l approach zero, for the envelope i s a lways decreasing so long as E / 0, Now, from equation (3.14) lim p = lim [- i - l o g e A ( t )] e —»- 0 e — 0 m lim l o g e A ( t ) _ e —- 0 lim t e —— 0 - lim log A(t ) lim e = e — 0 e m e — 0 e —»- u = lim e e^-0 = 0 . F i n a l l y , from equation (3.19) V 1 - e -co, lim = l i m ~2~t— e—*-0 e-*-0, t --co m ' m = 0 , provided ^ CO , However, we are guaranteed that ^ oo , for i f i t i s , we have an i n f i n i t e frequency, which i s not l i k e l y to occur i n the physical systems with which we are dealing. Thus, the l i m i t i n g values for both to^ and p are consistent• 4 1 Soudack then proposed a method of finding by con-(9) sidering these l i m i t s as e—*-0, - As e—*-0, the d i f f e r e n t i a l s equa t i on becomes x + x + j ix 3 = 0f x ( 0 ) = 1 , x ( 0 ) = 0 * ( 3 . 2 0 ) Now, s i n ce t m ~* - a? f P - * * 0 , and °* "the approximant ( 3 . 3 ) degenerates to the form x ( t ) = COS fiO^t and the f requency of o s c i l l a t i o n i s <o^. The exact s o l u t i o n to ( 3 * 2 0 ) i s , however, g i ven by the J a c o b i a n e l l i p t i c cos ine x ( t ) = Cn (k, ttt) ^ where k 2 « • • > ^ 1 11 t 2(1 + jO and Prom books o f t a b l e s , f o r example, Jahnke and Emde ,^^^ the qua r te r p e r i o d K(k) of the o s c i l l a t i o n can be f ound . Now can be chosen such tha t the degenerate s i n u s o i d a l case has t h i s same q u a r t e r p e r i o d » S ince K(k) = « £ =yrrj z and w l 4 = 2 f o r the co s ine t o be zero at ^ , we r e q u i r e *1 = 2 K & ) ^ " + 7 < (3-21) 42 To complete the proof of consistency of the solution i n the degenerate cases, we need to consider the l i m i t i n g cases of the functions A(t) and O(t) as \i ——Oo In order to do t h i s , the angle c r i t e r i o n must be re-examined* Prom equation (3.6), A(t ) = / - £ - , m' J 2|x which i s not v a l i d for very small u,, because then the value of A(t f f l) i s very much higher than 1. However, since A ( t m ) i s the amplitude at the point of t r a n s i t i o n to a linear solution, and this point occurs at t = 0 i f u,—*-0, we have lim A(t ) = A(0) = 1. — 0 m Now, from equation (3.1l), since ii [i - A(t )] lim t = lim m | i — 0 0 u. — 0 , A ^ m ) ^ 1 1 0 e = 0, the amplitude curve coincides with that of the linear solu-t i o n As a re s u l t * the amplitude function A(t) of the non-linear solution degenerates properly to the amplitude of the linear solution as u,—*-0. Considering the phase, since t f f l = 0 as [i 0, 0("t) =J1 - E 2 t + 0 q for a l l time. There-fore <o.t2 + tt.t + « =Jl - e 2 t + 0 (3.22) 2 1 o * o for a l l time. Since t°, t"'", and t 2 are l i n e a r l y independent, the only consistent solution of equation (3.22) for a l l time i s 43 • 2 = 0 2 e to = 0 = 0 o o Hence as a ——0 the parameters co^ and ft)^ degenerate properly to the linear solutions However, the proper value of 0 Q i s —1 E -tan , 2 i a n& n o " t 0° This discrepency arises from the J1 - e assumption that the i n i t i a l value of the amplitude i s unity, which leads to <0 q = 0. Since the solution i s only approxi-mate and the c r i t i c a l parameters are and (Og, this discrepency w i l l be toleratedo The error thus introduced w i l l be small because for a> as high as 0„2 radians, costt o 5 1 o = 0»98o This completes the proof of the consistency of the non-linear solution with the known solutions i n the degenerate cases where E —»- 0, or u—» - 0 . In conclusion* the functions A(t) and O(t) are obtained i n the forms A(t) = e - p t, and A ( t ) = <o + «Lt + < 0 o t 2 . o 1 2 The parameters p» ( 0 q > <O^, and can eas i l y be calculated from equations (3*14)$ ( 3 » 1 8 ) , (3ol9) and (3 D 2l) for specified values of E and u* The value of p i s obtained by making the amplitude equal to A ( t m ) at t = t m , and the value of <0 q i s obtained from the i n i t i a l condition x(0) = 1» Values of co^ See section 3 . 2 „ 6 for i n i t i a l value correction of the ampli-, tude*. 44 and U)^ a r e found by l e t t i n g e— » - 0 , and by matching the phase to the f i r s t derivative, which can be regarded as the frequency* Consistency with the known solutions of the degenerate cases where e — 0 or jj, —»-0 has been shown from the li m i t s of A(t) andO(t), except for a small error i n <Dq, which arises from the assumption that the i n i t i a l value of the ampli-tude A(t) i s unity. 3.2,4 Determination of P and 0 q After A(t) and O(t) are determined, i t i s a r e l a t i v e l y simple matter to f i n d P and 0 . In fac t , 0 can be calculated c o o di r e c t l y from equation (3.17a), i . e . 0 o = Vm + V i / " V l ~ e 2 K + " o = (^ ~Jl - e 2 ) t + <*0t 2 + (a m 2 m o But from equation (3«16a) , e - 2 <0 ot . 2. m Therefore 0 = (/1 - e 2 - 2<o_t ~J\ - e 2 ) t + » 0t 2+ •„ o 2 m ^ m 2 m u = - < 02 tm 2 + t to ( 3 ' 2 3 ) P i s found by matching the amplitude of (3.3) and (3.4) at t = t . Thus, -et A ( t m ) = P: e m . et Hence P = A(t ) e m . (3.24) m 45 This completes the determination of the parameters of the approximant to the solution of the non-linear equation (2.4) where e <C 1• An example using this approximating scheme i s worked out i n the next section. 3.2.5 Example In order to i l l u s t r a t e the approximating scheme just developed and to see how good i t i s . an example i s worked out. Because t h i s approximating scheme has a t o t a l phase which i s a quadratic i n t* i t w i l l hereafter be referred to as the parabolic phase approximation, a notation f i r s t used by ( l l ) Soudack, ' Consider,then, the equation x + 0.8 x + x + 3x 3 = 0, x(0) = 1, x(0) = 0 which has e = 0.4 and a = 3.0. Both the magnitudes of e and u are inadmissible i n the K-B method as we s h a l l see when we compare the solution with the true numerical solution. Using the parabolic phase approximation, however, a much better solution i s obtained. These three solutions, i . e . , the true numerical solution, the K-B approximation, and the solution obtained from the parabolic phase approximation, are shown i n Figo 3.7* F i r s t , from Appendix B, the K-B method yiel d s the following approximation: x(t) = e ~ e t cos (1 + 2|)t (B.7) or x(t) = e ~ ^ 0 ^ cos (2.125 t) . Using the parabolic phase approximation, on the other F i g . 3.7 Comparison between approximating methods 47 hand, we have f o r 0 < t < t m x ( t ) = e ~ p t cos (» 1t + ^ 2 ^ 2 ) f o r t > t m x ( t ) = p cos ( Jl - e 2 t + 0 q) The unknown parameters are then obtained from the equations developed i n the l a s t four s e c t i o n s as f o l l o w s : (3.6) 0.4 2T3T = 0.258 m e l o g e Jl - e2 A ( t m ) all - A ( t m ) ] lOe (3.12) = - ^ l o g e [ y o - 8 4 - (0.258)] + 310^7421 = 4,17 P = " T~ l0Sa A ( t m ) m (3.14) = - 47YT l o g e ( 0 t 2 5 8 ) = 0.325 (o = 0 o (3.18) -n; <°i = 2 K T k T y i + * (3.22) 0 3 (lO) where k = 2 ( i + n ) = 2(1+3) = C ^ 3 7 5 ' a n d f r o m J a h n k e a n d E m d e > 48 K(k) = 1.761* Therefore ten <6. 3,142 2(1.761) V 4 = 1-783 2 t. m ./0.84 - 1,783 2(4.17) .0,104 P - A(t ) m m = 0.258 e 0 ° 4 < 4 * 1 7 > (3.19) (3.24.) = 1 = 37 0. = =W-t + 0) 2 m o = 0,104(4.17) ' = 1 a 81 a (3,23) F i n a l l y , the complete approx imat ion becomes f o r 0 < t < 4 . 1 7 x ( t ) = e " 0 e 3 2 3 t c o s ( l . 7 8 3 t - 0 „ 1 0 4 t 2 ) f o r t > 4 » 1 7 x ( t ) = 1.37 e " 0 , 4 t cos (0 .916t + 1.81) Nov, a comparison between the two app rox ima t i on s v as shown i n F i g « 3.7? i n d i c a t e s tha t the K=B method i s exceed ing l y s imple to c a r r y out* but the r e s u l t i s poo r . The f requency of o s c i l l a t i o n i s too l a r g e , and there i s no phase r e t a r d a t i o n . A l s o , the ampl i tude of o s c i l l a t i o n decays too r a p i d l y . On the o ther hand, the p a r a b o l i c phase approx imat ion r e q u i r e s a few more s imple computat iona l steps but the e x t r a e f f o r t i s w e l l rewarded* The f requency of o s c i l l a t i o n i s now sma l le r than 4 9 that obtained by the K-B method and there i s phase retardation i n the f i r s t part of the solution, where the effect of the non-linear term cannot be neglected. Also, i n t h i s f i r s t partj the amplitude decays more slowly than e as already predicted from previous studies of the equation. Moreover, examples with s t i l l higher e and/or u w i l l show that the parabolic phase approximation i s far superior to the K—B-method in dealing with these types of non-linear equations. 3.2.6 Refinements i n the Approximation Although examples such as the one considered i n the l a s t section have indicated that the parabolic phase approxi-mation yields far better results than the K-B method, close examination of these examples suggests that some refinements i n the method would make the results even better. The f i r s t refinement involves no extra labour and i s essentially a modification based on the consideration of the phase term. The second one i s a correction of the i n i t i a l value of the amplitude. Prom various examples, the amplitude function A(t) obtained i s i n f a i r l y good agreement with the numerical solution. It i s the phase function Q(t) that contributes mainly to the discrepency in the solution. Examination of many examples reveals that the phase O(t) always leads that of the true solution, indicating that either i s too large or the phase retardation too small. Before we attempt to improve the phase function, l e t us review how i t i s obtained and compare i t to i t s true value. F i r s t , the angle c r i t e r i o n i s used to determine A ( t m ) and consequently t m , at 'which the phase retardation i s considered n e g l i g i b l e . Theny the f i r s t time derivative of the phase, i . e . the frequency, i s assumed to decrease l i n e a r l y as time increases for t <Ct m For t = t m , the frequency i s assumed to be Jl — e" which i s the frequency of o s c i l l a t i o n of the solution to the complementary lin e a r equation* Thus, as shown in F i g . 3.8(a), the graph representing the frequency begins at the point M(0, <6^), drops l i n e a r l y to the point N ( t m j J1 - e 2 ) and becomes level thereafter. The area under this graph then represents the approximate t o t a l phase, as shown i n F i g , 3.8(b)* But since the approximate phase i s always leading the true phase as already pointed out* the true phase may be represented by the dotted curve i n Fig* 3.8(b). It i s always below the approxi-mate phase and approaches an asymptote with a slope of Jl - e 2 . The true frequency may, therefore, be represented by the dotted curve i n F i g . 3.8(a). It becomes clear now how the discrepency arises. Apparently, the straight lines MN and NH do not approximate the true frequency too well for t < t , and consequently do not give a p a r t i c u l a r l y good approxi-mation to the tota.l phase. In an attempt to improve the approximation, l e t us consider the point L(t /29 Jl -* e 2 ) as shown in F i g . 3*8(a). The lines ML and LH would give a better approximation to the true t o t a l phase because the area KLN under the dotted curve would compensate for the area JMK ¥ : = ~ A small error i s present here, because the i n i t i a l value of the phase i s s l i g h t l y different from zero. However, for the purpose of finding an approximation, we have assumed, that A(0) = 1, which leads to to = 0, or 0(0) = 0. 2 51 (a) Frequency (b) Phase F i g , 3.8 Comparison of the true frequency and phase with their approximations 52 over i t . The point L i s so chosen partly because the approxi-mating scheme already developed can be adopted with p r a c t i c a l l y no change. We need only match the li n e a r and the non-linear parts of the approximant at t /2 instead of at t • Following the above' argument* i t might be noted that this new matching point could have been di f f e r e n t from t /2, such as ir t , r m 3 m' 3 2 7 t , or 7 t » etc. and that matching at t 12 does not 4 n r 5 m7 ' 5 nr necessarily give the best r e s u l t . However, we must not forget that the objective of this section i s to modify the parabolic phase method i n order to give improved results i n general, and not optimum results i n particular cases. Since various numerical examples have shown better results by matching at t m / 2 , we now replace t i n the equations previously developed by t /2 and obtain the following equations: Jl - e 2 - «, From equation (3*19), « 2 = 1 * (3.19a) m From equation (3.23). 0Q = f t + «>o . (3.23a) t | t From equation (3.24), P = A (-f) e* m 2 m = e Therefore, the complete approximation becomes for 0 < t < t /2 x(t) = e ~ p t cos(w + «. t + < 0 o t 2 ) . ^ nr ' o 1 2 for ±^tm/2 x(t) = P e ~ e t cos (Jl - e 2 t + 0 ) . (3.24a) 53 The other refinement i s to improve the amplitude function A ( t ) . So f a r , the i n i t i a l value of the amplitude has been assumed to be unity, which i s not quite correct. Prom various examples, i t has been observed that the true i n i t i a l amplitude i s greater than unity and the difference between the true i n i t i a l amplitude and unity decreases as u, increases. For most cases where e i s not too large or [x i s greater than 3, this difference i s n e g l i g i b l e . However, i f e becomes close to unity or \i becomes smaller than 3, this difference w i l l be appreciable and a correction added onto the assumed i n i t i a l amplitude of the approximant w i l l d e f i n i t e l y improve the r e s u l t . Since t h i s difference i s greatest for \i = 0, and becomes negligible for [x = 3, we may assume as a f i r s t approximation* that i t drops l i n e a r l y as [x increases from 0 to 3. Knowing that the true i n i t i a l amplitude i s —====. i f / l - c u, = 0, we therefore obtain the following relations 1 I n i t i a l Amplitude Correction = — ^ e - l Thus, for ,u,^3, A(0) = 1 + - l (3.25) Denoting this value of the i n i t i a l amplitude by A Q, we have A(t) = A Q e " p t , which from equation (3.14) leads to 1 A ( tm } P = - — l o g e -j . (3.14a) m o Prom the i n i t i a l condition that x(0) = 1, ¥e also have 54 1 = A cos a> , o o # or w = c o s ~ 1 ( l / A ), » <0 (3.18a) o ' o o ^ F i n a l l y the non-linear part of the approximant becomes for 0 < t < t /2 x(t) = A e cos(w + «, t + at,) m o o JL w As a wholej these refinements i n the approximation w i l l y i e l d better results and can be made with almost no extra e f f o r t , for the general development of the procedure i s not changed. Two examples w i l l be worked out i n the next, section to i l l u s t r a t e t h i s refined parabolic phase approximation. 3.2.7 Summary and Examples of the Refined Parabolic Phase Approximation In order to f a c i l i t a t e the application of the refined parabolic phase approximations, a summary of the computational procedure i s given as follows: (1) Normalize the equation into the form . 3 x + 2ex + x + ux = 0. (2) Compute A ( t m ) from the angle c r i t e r i o n , i . e . equation ( 3 . 6 ) 2u (3) Compute t m from equation (3.12) Jl - e 2 A(t ) t = - — log .„ _ m e toeL m'j - A(t m)] 106 ft Positive values of w are not used, as w i l l be explained i n section 3.2.8* (4) If a <31 compute A from equation (3.25) 55 A = 1 + 2=lt o 3 and assume A Q =1 i f a ^ 3 . (5) Compute p from equation (3.14a) p = log m e A. (6) Compute » q from equation (3<>18a) -1 ,1 s (0 = cos {-.— J o VA o w <0. 0 ^ (7) Compute k from k 2 = 2(i+a) 9 a n d "k n e n obtain the quarter period K(k) from tables of e l l i p t i c functions. (8) Compute co^ from equation (3.22) <°1 = 2KTkT" J 1 + * • (9) Compute from equation (3.19a) J l - e 2 -<°2 " m (10) Compute P from equation (3.24a) t | t i , , / m\ Z m P = A(—2) e ( l l ) Compute 0 q from equation (3.23a) 0 = 7 t 2 + <o 'o 4 m o The complete approximation f i n a l l y becomes for 0 < t ^ t /2 x(t) = A e ~ p t o x(t) = P e •et for Since one can quickly arrive at answers within slid e - r u l e accuracy, the method seems to be very suitable f o r p r e l i m i -nary engineering analyses, and as an added benefit i t w i l l give the engineer some useful insight into the behaviour of the system. The same equation considered i n section 3.2.5 i s again used so that the advantage of the refined method can be i l l u s t r a t e d * The equation has e = 0.4 and \i = 3. Following the steps just outlined, we obtain Example x + 0.8 x + x + 3x 3 = 0, x(0) = 1, x(0) = 0, t m - fa l o g e [J0^4 (0.258)] + ZlS^ml = 4.17 A 1 o P o l o g e ( 0 „ 2 5 8 ) = Oo325 K(k) = 1.761 3*142 Jl + 3 = lo783 1 ~ 2(1.761) ^ =/b-84- i1.783 = _ 0 o 2 0 8 P = e ^ 0 e 4 ~ °-325)(2.085) = x 1 ? 0 = 0±|08 ( 4 e l 7 ) 2 = O o 9 0 5 Hence the complete approximation i s as follows! for 0<t<2.09 x(t) = e-0.325t c o s ( l o 7 8 3 t _ o.208t 2) for t>2«09 x(t) = 1.17 e ~ 0 o 4 t cos(0.9l6t + 0.905) This approximate solution i s plotted i n Pig. 3.9 together with the numerical solution, the K-B approximation, and the unrefined parabolic phase approximation, which are ob-tained i n section 3.2.5. It i s observed that the refined parabolic phase approximation i s the closest to the numerical solution. Example As another example of the refined parabolic phase approximation, consider the equation x + k + 5x + 10 x 3 = 0, x(0) = 1, x(0) = 0. Since normalization i s required i n this example, l e t 7" = ,/~5t , to obtain k = 2.236 x» x = 5 x 1' cLx cL^ x where x' = -r— and x " = — ~ . dT d T 2 Substituting x' and x°' into the o r i g i n a l equation and x + 0.8x + x + 3x 3 = 0, x(0) = 1, x l(0) = 0 KB approximation Parabolic phase approximation (unrefined) True numerical solution Parabolic phase approximation (refined) F i g . 9 Approximation by the refined parabolic phase approximation for the equation x + 0.8x + x + 3x 3 ='0, x(0) = 1, x(0) =0 00 59 d i v i d i n g through by 5, we have x 1' + 0.4472 x 1 + x + 2x 3 = 0, •which has e = 0.2236 and u = 2. Carrying out the computational steps, we obtain *m - " 0^6 l o g j y i - 0.05 (0.236)] + = 7.26 A = 1 + ~"\-±— - l l = 1.008 0 3 lj^95 J * = - 7 T 2 6 l 0 g e [ M l ] - 0 - 2 0 0 tt0 = c o s~ 1 i r f e s = - °-13 K(k) = 1.734 3.142 " i = 2(1.734) V T T T = . U 5 6 8 . B ybT9T-21.568 a _ 0 > Q 8 1 7 p = e(0.2236 - 0.2)(3.63) = l e Q 9 0 = 0.0817 ( 7 ( 2 6 ) 2 - 0.13 = 0.945 0 4 Hence the complete approximation i s for 0 < t ^ 3 o 6 3 x ( T ) = 1.008 e 0 o 2 T c o s ( l . 5 6 8 T - 0.0817T 2 - 0 . 1 3 ) for t > 3 . 6 3 x(r) = 1.09 e " o 2 2 3 6 T c o s (0.975T + 0.945) But since 7~ = 2.236t, the f i n a l approximation becomes for 0 < t ^ l . 6 2 x(t) = 1.008 e " 0 , 4 4 7 t c o s (3. 508t - 0.409t 2 - 0 . 1 3 ) for O l , 6 2 x(t) = 1.09 e~°" 5 t c o s (2..l8t + 0.945) Now, l e t us compare this so lut ion with the K-B approximation. Prom Appendix B, the K-B approximation i s given by or x (t) = e" e t cos (1 + 2%) t ( t ) = e " 0 o 5 t cos ( 3 . 9 1 3 t ) . (B.7) These two solutions are displayed i n F i g . 3.10 together with the numerical so lut ion . Results are, again, i n favour of the ref ined parabolic phase approximation. 3.2.8 Errors and Limitat ions An important uncertainty inherent i n a n a l y t i c a l (12) approximating methods i s the error i n the so lut ion obtained. Because of the presence of the non-l inear term, i t i s usual ly not a simple matter to make an error ana lys i s . After the approximate so lut ion x i s obtained by using the parabol ic (8) phase approximation, a common c r i t e r i o n for the system error i s given by the in tegra l F i g . 3.10 h - 1 62 •t. co m J = o - 2 ( t ) d t + CT 2(t)dt 0 m where cf(t; = x + 2ex + x + ux . As already discussed i n section 3.2.2, the evaluation of thi s integral i s a formidable job, and leads to no immediate insight into the accuracy of the solution. It must be noted that t h i s type of error analysis does not require the knowledge of the true solution. I f , on the other hand, we knew the true solution,, i t would be very simple to measure the error. The absolute deviation, which i s defined as the magnitude of the difference | between the approximate and true solutions, could be computed and plotted against the independent variable. Ve could then have at our disposal a number of quantities as measures of error, for examples (l) the maximum deviation or (2) the area under the deviation curve* In this work, the numerical solution obtained from the d i g i t a l computer i s considered to be the * true solution and i s compared with the refined parabolic phase approximation* The maximum deviation i s chosen as the measure of error and i s evaluated for a large number of examples with a wide range of s and u. If the maximum deviation i s found to be s u f f i c i e n t l y small, we can conclude, because of continuity, that the approximating method i s satisfactory for every £ and u within the range. The following table shows a few numerical results as compared to those obtained from the K—B method. The machine solution i s accurate to three decimal places. See Appendix A. 63 Maximum Deviation e K-B method Refined Parabolic Phase Approximation 0.2236 2 9 3 3 9 3 3 0.64 1.04 0,66 0.44 0.80 0.48 0.50 0.13 0.23 0.10 0.07 0.16 0.03 0.10 0.3 0.4 0.6 0.6 0.8 0.9 As suggested from the above table, the l i m i t s of e and \i are not very d e f i n i t e , because they depend on each other as well as on the specified accuracy. For example, higher values of u. may be accepted i f the value of e i s higher, while smaller values of e make the allowable value of u. lower. This i s reasonable because as the damping becomes l i g h t e r , the non-linear effect takes a longer time to become ne g l i g i b l e , and the parabolic phase i s not suf-f i c i e n t to ensure.a good phase f i t . For a maximum devi-ation of about 0.1, various experimental results have shown that \i may be as high as 5 i f e < 0.5 and as high as 10 i f e >0.5. In most cases, better accuracy can be expected i f u- i s not so large. Thus, the above error consideration gives us an idea of the upper l i m i t of [x. The lower l i m i t of J A , however, i s determined by one of the approximating steps. Early i n the approximating procedure, the amplitude at the point of t r a n s i t i o n from the non-linear to the linear solution i s determined by equation (3.6), i . e . 64 Since the amplitude cannot be gr e a t e r than 1, the r e l a t i v e value of e and a should be such that or 'a ^ 2 . Therefore, i n the case where a < ^ , the p a r a b o l i c phase approximation becomes i n a p p l i c a b l e and we can conside r that the equation i s e s s e n t i a l l y l i n e a r because the angle c r i t e r i o n has i n d i c a t e d that the n o n - l i n e a r e f f e c t i s n e g l i g i b l e from the very s t a r t . Thus, we conclude that an acceptable lower £ l i m i t of u i s 2 = Another e r r o r that i s also i n h e r e n t i n the r e f i n e d p a r a b o l i c phase approximation i s that the i n i t i a l c o n d i -» , \ **** t i o n x(0) = 0 i s u s u a l l y not met. I f x i s d i f f e r e n t i a t e d we have x ( t ) = ~pA e- p tcos(<o + » 1 t + » 0 t 2 ) - 2tt 0t)A e" p tsin(«> 0 0 1 2 1 2 o 0 2. + <o1t + <o 2 t ) Therefore, x(0) = -pA cos w - <o,A sin w ' r o o 1 0 o - _p _ to, A s i n 0) . 1 o o Assuming p^?jto^A osin <OQ we know t h a t the maximum e r r o r i n the i n i t i a l slope w i l l never exceed —p i f we do not allow p o s i t i v e values of a>o from equation (3.18a). This i s the reason why i O Q i s e i t h e r negative or zero. Thus, we. see that an e r r o r i n the i n i t i a l slope i s inherent i n the method, but t h i s e r r o r can be ignored because we are p r i m a r i l y i n t e r e s t e d i n the solution of x(t) for a wide time range. This completes the development of the approximation for the case where e < 1. 3.3 Case II - e >1 3.3.1 Choice of Approximant In this case, where e ^ l , the solution may or may not have overshoots depending on the r e l a t i v e values of e and u., as already discussed i n section 2.3. Since numerical solutions obtained from the d i g i t a l computer for s = 1.1 have shown that \x has to be higher than 25 be-fore a second overshoot appears, and a s t i l l higher (j, w i l l be required for a second overshoot i f e i s larger than l o l , we need only consider the case with at most one over-shoot i f we l i m i t our interest to p, ^ 1 0 . Because the solution i s not o s c i l l a t o r y i n general, we can no longer assume an approximant of the form x(t) = A(t) cos O ( t ) , and consequently, both the parabolic phase approximation and the c l a s s i c a l K-B method are not applicable. There-fore, a new method of approximating the solution must be developed. To this end, consider the "complementary" li n e a r equation x + 2ex + x = 0, x(0) = l j x(0) = 0 where e ^ 1 . As indicated i n chapter 2, i f e =1, the solution of this equation i s given by 66 x(t) = (1 + t)e'" t , (2.7) and i f e > l , the solution becomes x(t) = \ (1 + e ) e ( ~ E + V / ^ r i ) t + 1 ( 1 ( _ e _ y e 2 _ i ) t E ~ 1 (2.8) Let us denote, hereafter, the solution to the complementary linear equation by x ( t ) , which i s given by equations (2.7) c and (2.8), and examine the effect of introducing a non-3 lin e a r term [i.x; to the equation. Consider then the equation x + 2ek + x + px 3 = 0, x(0) = 1, x(0) = 0 whose solution i s x(t) = x (t) - z ( t ) , (3.26) where z(t) i s a correction term to account for the effect of the non-linearity. Thus, the problem now i s to approxi= mate z ( t ) . In order to determine the form of z(t) for the approximation, solutions x(t) to various examples have been obtained numerically from the d i g i t a l computer and z(t) i s then calculated from equation (3.26)* i . e . z(t) = x (t) - x ( t ) . F i g . 3.11 shows a t y p i c a l example of z ( t ) . Close exami-nation, of the general shape of z(t) has revealed i t s characteristics from which possible approximants are sug-gested as follows: (1) Before z(t) reaches i t s maximum* i . e . for t < t , i t may be approximated by the function t n f where n >1. (2) For t > t p 9 since z(t) decreases as t increases, the possible approximants are t m or e , where m < 0. (3) z(t) i s always positive and vanishes at t = 0 and t = 00 . Z(t) 0.4 — 0 1.0 t 2.0 3.0 4.0 5.0 F i g . 3.11 Correction term z(t) for e = 1.2 and a = 2 Let us therefore examine the function F(t) = t n e _ t , n >1, 68 to see i f i t has a l l the above features; (1) F(t) has only one maximum, i . e . at t = -n. For small t, F ( t ) ^ t n . (2) For large t, e ^ i s the dominating factor. (3) F(t) i s always positive and vanishes at t = 0 and t =00. Since this function has a l l the features z.(t),hasj we w i l l assume that z(t) takes the form z(t) = g t n e~* where g i s a constant. Our problem now i s to determine the constant parameters g and n i n terms of e and p. Discrepencies may arise from the assumption that m = -1, but from the accuracy of the results which w i l l be developed l a t e r , they may be either too small to be of significance or else may have been taken up by the other factor t n . As a result, the form of the approximant of z(t) becomes z(t) = g t n e""5 (3.27) and the rest of this work w i l l be devoted to the determina-^ t i o n of g and n as functions of e and J A . 3.3.2. Determination of n Empirical results are used to determine both the parameters n and g. In order t o determine n, consider the equation (3.27). The maximum value of z ( t ) i s given by dt - u ? , n - 1 - t , n - t „ o r g n t e - g ' t e = 0 . S i n c e g , t , a n d e ^ a r e n o t z e r o f o r f i n i t e t , we o b t a i n t = n , ( 3 . i . e . , n i s n u m e r i c a l l y e q u a l t o t h e t i m e a t w h i c h z ( t ) i s a m a x i m u m . T h e r e f o r e , i f we i n s i s t t h a t z ( t ) h a s a maximum a t t h e same t i m e as z ( t ) . we c a n f i n d n f r o m e x p e r i m e n t a l r e s u l t s b y s i m p l y n o t i n g t h e t i m e a t w h i c h t h e maximum o f z ( t ) o c c u r s , as s h o w n i n F i g . 3 . 1 1 . F r o m e x a m p l e s w i t h v a r i o u s v a l u e s o f e a n d a , t h e f o l l o w i n g t a b l e i s o b t a i n e d ; N u m e r i c a l V a l u e s o f n 0 . 5 1 . 0 1 . 5 2 . 0 2 . 5 3 . 0 1 . 0 1 . 6 2 1 . 5 5 1 . 4 9 1 . 4 4 1 . 4 0 1 . 3 7 1 . 2 1 . 7 1 1 . 5 8 1 . 5 0 1 . 4 7 1 . 4 0 1 . 3 8 1 . 4 1 . 7 2 1 . 6 2 1 . 5 5 1 . 4 9 1 . 4 2 1 . 4 0 1 . 6 1 . 8 7 1 . 7 3 1 . 6 0 1 . 5 4 1 . 4 3 1 . 4 0 1 . 8 2 . 0 2 1 . 7 8 1 . 6 8 1 . 5 9 1 . 5 3 1 . 4 2 2 . 0 2 . 1 3 1 . 9 0 1 . 7 6 1 . 6 9 1 . 6 0 1 . 4 9 2 . 2 2 . 3 9 2 . 1 0 1 . 8 7 1 . 7 4 1 . 6 6 1 . 5 5 2 . 4 2 . 5 0 2 . 2 4 1 . 9 7 1 . 8 7 1 . 8 1 1 . 7 1 A n e r r o r o f t h e o r d e r o f 0 . 0 5 may be e x p e c t e d i n some o f t h e s e f i g u r e s b e c a u s e some c u r v e s o f z ( t ) h a v e a r a t h e r f l a t p e a k a n d i t i s d i f f i c u l t t o l o c a t e t h e maximum a c c u -r a t e l y . A t a n y r a t e , t h e s e f i g u r e s g i v e a g o o d p i c t u r e o f how n v a r i e s w i t h b o t h e a n d a . I f t h e c o n t o u r s o f c o n -s t a n t n a r e p l o t t e d i n t h e e - u p l a n e , F i g . 3 . 1 2 i s o b t a i n e d . T a k i n g t h e p o s s i b l e e r r o r o f n i n t o c o n s i d e r a -t i o n , we may now a p p r o x i m a t e t h e s e c o n t o u r s b y a s e t o f p a r a l l e l s t r a i g h t l i n e s a s s h o w n i n F i g . 3 . 1 3 . The s l o p e o f t h e s e s t r a i g h t l i n e s i s f o u n d t o be 0 . 4 7 2 . T h e r e f o r e , F i g . 3„12 Contours of constant n from experimental results F i g . 3,13 Approximation to the contours of constant n they can be represented a n a l y t i c a l l y by the simple r e l a t i o n e = 0.472u + c (3.29,a) where c i s the intercept and depends on n, Now, the i n t e r -cept c i s plotted against the corresponding n as shown in Pig. 3.14. 2.3 2.0 --1.5 -1.0 --0.5 --0 o Experimental Intercepts c = 2.3(1 - e - 2 - 0 ( n - U 4 ) ) n 1.4 1.6 1.8 2.0 2.2 2.4 F i g . 3.14 Determination of c as function of n From this diagram, we observe that the curve exhibits sat-uration at about c = 2.3 and intercepts the abscissa at n = 1.4. Considering also the shape of the curve, we are led to believe that c can be expressed i n terms of n as 73 follows: c = 2.3 (1 - e ^ ( n ~ U 4 0 ) ) where q i s a constant. By t r i a l and error, a q of -2.0 has been found to give good results as i l l u s t r a t e d i n Pig. 3.14 where the function c = 2.3 (1 - e- 2'°( n " i ' 4 0 * ) i s also plotted. Now substituting c into equation (3.29a), we have e = 0.472a + 2.3 ( l - e " 2 - ° ( n - U 4 ) ) which yie l d s n = 1.4 + (0.5) log — 2-2^ (3.29) e 2.3 + 0.472u-e Hence n can be computed when e and a are spec i f i e d . 3o3.3 Determination of g After n has been obtained, g can be determined by making zf(t) equal z(t) at the maximum, i . e . at t = n, from equation (3.27). The maximum i s chosen as the matching point because we have determined the parameter n, such that the peak of the approximant z(t) occurs at the same time as that of z(t)« Since the maximum occurs at t = n, we have g = z(n) n- n e 1 1 (3,30) Using the same example as i n F i g . 3.11, where e = 1.2 and a = 2, we have n = 1.4 + 0.5 log — — — 6 2.3 + 0.944 - 1.2 74 = 1.46 and g = 0.325 ( l . 4 6 " l o 4 6 ) e 1 * 4 6 = 0.805 • Hence z(t) = 0.805 t 1 " 4 6 e~ + . This approximation i s now shown i n Pig. 3.15 together with the true z(t) obtained numerically. The result i s encouraging because very good agreement between z(t) and z(t) i s observed, indicating that the form assumed for z(t) i s a good one. 75 Examples with various values of e and (j, are then invest igated i n a s imi lar manner, and the corresponding values of g are shown i n the fol lowing table: Numerical Values of g 0.5. 1.0 1.5 2.0 2.5 3.0 1.0 .275 .500 .680 .840 1.000 1.14 1.2 .260 .455 .651 .805 .946 1.07 1.4 .238 . .420 .584 .731 .862 .978 1.6 .209 .385 .541 .681 .802 .912 1.8 .184 .351 .493 .629 .742 .848 2.0 .159 .314 .454 .578 .691 .793 2.2 .131 .274 .408 .533 .640 .743 Our next task is to determine ,'g as a function of E and n from this tab le . The f i r s t attempt was to use the contours of constant g i n the e-p, plane as we d id before i n the determination of n . The contours of constant g were then plot ted as shown i n F i g . 3.16. From this diagram, we observed that the contours could not be represented by a set of simple functions such as p a r a l l e l s tra ight l ines because there i s a def in i te trend showing that the slope of each l i n e i s d i f ferent from the r e s t . In another attempt, g i s now plot ted against u for constant values of e, as i n F i g . 3.17. From the shape of the curves obtained, g i s seen to have the form b g = a u. for constant E , where a and b are constants. Therefore, for constant E , the p lot of log^Q g against log^Q^ i s a set of s tra ight l ines as shown i n F i g . 3.18. Hence we can write l o g 1 0 g = b l o g l O [ A + l o g l O a (3.31) 76 e , s t»o 1-^ 2-0 2>5 3 0 u Pig. 3.16 Contours for constant g g(e,u) 77 1.2 1.0 --0.8 0.6 0.4 0.2 --0 0.5 1.0 1.5 2.0 2.5 3.0 F i g . 3.17 g as a function of u for constant e Since these straight lines are almost p a r a l l e l and equally spaced, we may approximate this set of lines by F i g . 3.19 in which the straight lines are p a r a l l e l and equally spaced. Therefore, the constant b i s the common slope of a l l these straight l i n e s and i s found to be 0.794. Since the v e r t i c a l distance between the. straight lines for e e 1.0 and e = 2.0 i s 0.15, and the intercept for 1.0 i s -0.325, equation (3.3l) becomes l o g 1 0 g = 0.794 l o g 1 0 a - [o.325 + 0.15(e-l)] . (3.32) Pig. 3 . 1 8 log, ng vs log, 0u. for constant e F i g . 3.19 Approximation to F i g . 3.18 80 Thus, g can be computed when e and a are given. Let us now show that the result i s consistent i n the l i m i t i n g case where a —»-0. As a—*-0, the non-linear equation approaches i t s complementary lin e a r equation whose solution i s x ( t ) . Prom equation (3.32), we have c lim ( l o g 1 0 g ) = -CO u — 0 . or lim g = 0 . a—-0 Therefore lim z(t) = 0, a — 0 which yie l d s lim x(t) = x (t) * a-^0 c Hence, our approximate solution degenerates to the correct solution for the complementary equation* An example w i l l be worked out i n the next section to i l l u s t r a t e method* 3*3.4 Summary and Example of the Correction Term Approximation This approximating method, then, i s es s e n t i a l l y the determination of the solution x f i(t) to the complementary linear equation and the approximation of the correction term z(t) 3 due to the presence of the non-linear term ax . The pro-cedure i s summarized as follows; (l ) Determine the solution to the complementary lin e a r equation by equations (2.7) and (2.8), i . e . i f e = .1, x (t) = (1 +' \)<T% , 81 and i f e >1, x c ( t ) =|d + ^ r T ) e ( - ^ i ^ r T ) * + |(1 e (2) Approximate z(t) by ( _ e _ y e 2 _ i ) t z(t) = g t n e * , where 2.3 n = 1.40 + 0.5 l o g e 2.3 + 0.472a - e l o g 1 Q g = 0 . 7 9 4 l o g 1 0 u - [ o . 3 2 5 + 0 . 1 5 ( e - l)] The complete approximate solution i s then given by x(t) = x c ( t ) - g t n e~\ . An example i s now worked out to i l l u s t r a t e the method. Example Consider the equation + 2.8 x + x + 3x 3 = 0, x(0) = 1, x(0) = 0, x i n which e = 1.4 and a = 8. Following the steps just outlined, we obtained x (t) = I(! + 1*4 ^(-1.4+yo.96)t+ 1 ( 1_ 1.4 ) e(-1.4-yo.96)t c 2 70.96 2 ybT96 = 1.21 e " ° - 4 2 t - 0.21e- 2* 3 8 t 82 n = 1.40 + 0.5 log — 6 2.3+0,472(3)- 1.4 = 1.4 l o g 1 0 g = 0.794 l o g 1 0 3 - [0.325 + 0.15(1*4 - l ) ] = -0.006 g = 0.987 . F i n a l l y , the complete approximation i s x(t) = 1.21 e - ° * 4 2 + - 0.21 e - 2 - 3 8 t - 0.987 t ^ V * . This approximate solution i s compared -with the numerical solution i n F i g . 3.20 and i s seen to be quite satisfactory, the maximum deviation being 0.04. In the same Figure, the approximation by using the Ritz and i n i t i a l condition matching method i s also shown. Referring to Section 3.1, thi s solution i s obtained by solving numerically the equations 3 3 - M*2+ 2ea+l) + ~(b2+ 2eb+l) ^ ~ + ^ 5 a a + b 4a(a-b) 2 (3b+a)(a-b) 2 2 2 3 u-b a _ 3nab _ Q (3a+b)(a-b) 2 2(a+b)(a-b) 2 3 3 £(b^+ 2eb+l) + -rr(a + 2ea+l) - 0 -r 0 b a + b 4b(a-b) 2 (3a+b)(a-b) 2 2 2 + 3|xa b _ 3u-ba _ Q (3b+a)(a-b) 2 2(a+b)(a-b) 2 where e =1.4 and \i = 3.0. By extensive t r i a l s , a d i g i t a l computer produced the following roots: 83 a - -1.11 b = -4.44 Hence the solution i s ~ m . - l . l l t . , n -4.44t x^tj = A e + B e . Applying i n i t i a l conditions of x(0) =1, and x(0) - 0, we have A = -=£ = 1.333 a-b B = -—- = -0.333 a-b Therefore, the Ritz method and the i n i t i a l condition matching gives ~tj.\ -i o i l - l . l l t p. 0 0-, -4.44t x(t) = 1.333 e - 0.333 e As i l l u s t r a t e d i n Pig. 3.20, this solution has a maximum deviation from the numerical solution of 0.05, which i s not as good as that obtained by the correction term, method just developed. If we consider the p r a c t i c a b i l i t y and the e f f o r t required i n applying the Ritz and i n i t i a l condition matching method, i t i s evident that the correction term method i s much more tractable. 3.3.5 Errors and Limitations Following the same reasons as i n the case where e < 1, we again use the deviation between the true and the approximate solutions as a measure of accuracy to j u s t i f y the v a l i d i t y of the correction term method. The following F i g . 3 .20 Approximation by the Correction Term Method for e = 1.4. u = co table i s the result obtained from a large number of examples: e Maximum Deviation 1.0 1.0 0.02 1.0 3.0 0,07 1.0 7.0 0.24 1.4 1.0 0.02 1.4 3.0 0.04 1.4 8*0 0.21 1.8 1.0 0.04 1.8 2.0 0.06 1.8 3.0 0.08 2.2 1.0 0.05 2.2 2.0 0.09 2.2 3.0 0.11 Prom this table, ve see that the maximum deviation i s rather high for a. = 8. But i t must not be forgotten that the deviation i s generally much smaller than i t s maxi-mum. An example w i l l help c l a r i f y this point. The magni-tude of the deviation between the approximate and true solutions of the equation x + 2.8 x + x + 8 x 3 =.0, x(0) = 1, x(0) = 0 i s shown i n F i g . 3.21. Thus, we see that the deviation i s f a i r l y large for small t, then drops off quickly and remains well under 0.08. As a r e s u l t , we may allow a to be as high as 8 for e = 1.4. It is also suggested from the above table that the allowable a i s lowered i f e becomes larger. This i s reasonable because as e increases, one of the exponents (— — /z2— l~)t involved i n x c ( t ) , i . e . e « , becomes negligible compared to the other and the choice of e~^ i n the correction term w i l l not be a very good one. At any rate, for e as high EL S 2 • 2 j £L a of 5 may s t i l l be allowed. Thus, we obtain an idea as to the upper l i m i t of a from the above error con-86 Deviation sideration. The lover l i m i t of |x i s zero of course, because.-" we have shown i n Section 3.3.4 (page 80) that the approxi-mation x(t) degenerates to the solution x c ( t ) of the complementary equation. Regarding the value of e, however, there i s an inherent l i m i t i n the method. Consider equation (3.29): n = 1.40 + 0.5 log — , 6 2.3+0.472u-e F i r s t , n must not be i n f i n i t e . Therefore, we have the r e s t r i c t i o n that 2.3 + 0.472^ 1 - e / 0 . Secondly, since the logarithm of a negative number i s not allowable, the r e s t r i c t i o n then becomes 2.3 + 0.472u. - e > 0 . Knowing that the lowest value of jx i s zero, we see that for \i — 0, the upper l i m i t of e i s 2.3. However, since we are dealing with non-l inear equations, |x i s always greater than zero and the upper l i m i t of e i s , therefore, usually higher than 2.3 o F i n a l l y , i t may be worth mentioning that i n th i s approximation, the i n i t i a l conditions x(0) = 1 and a x(0) = 0 are met because z(0) = z(0) = 0 which lead to x(0) = x (0) = 1 c ' and x(0) = x (0) = 0 . This concludes the correct ion term approximating method which has been developed for e ^ 1 . 3 04 Summary In th i s chapter, two methods have been developed to. approximate d i r e c t l y the so lut ion to the equation x + 2ex + x + u-x3 = 0, x(0) ~ 1, x(0) = 0 where both e and jx are not small numbers. F i r s t , i n the case where e < 1, the parabolic phase approximation was developed and re f ined . A heur i s t i c argu-ment was given for the use of the form for 0 < t ^ t x(t) A(t) eos Q( t ) for t ^ t x(t) = P e~'e t eo s ( / l - A + 0 ) ' m w o as the approximant. The value of t where the non-l inear effect- becomes neg l ig ib le was determined after A(t ) had been obtained using an a r b i t r a r y c r i t e r i o n based on the consider-at ion of the angle between phase plane i soc l ines for the l i n e a r and non-l inear equations. The amplitude A(t) was assumed to have the form A(t) ~ A Q e-P* , and a parabol ic phase of the form o f l ( t ) = toQ + w^t + tt>2t . The value of co was determined from the i n i t i a l con-o d i t i o n x(0) = 1, then and were found by l e t t i n g r.—*- 0 and by matching the phase to the f i r s t d e r i v a t i v e . The value of A was f i r s t assumed toi be unity and p was o obtained by making the amplitude equal to A ( t m ) at t = t^. The method of obtaining these parameters was l a t e r ref ined by using t = "^ m / 2 a s ^ n e point of t r a n s i t i o n i n s -tead of t = t . and by correct ing the i n i t i a l amplitude. The approximation then became for 0 < t ^ t /2 x(t) = A Q e~ p t c o s ( « + » x t + » 2 t 2 ) for t ^ t /2 x(t ) = P e" e t cos{Jl - s 2 t+0 ) . ^ m * o Hence, the parameters P and 0 Q were determined by matching the two parts of the approximate so lut ion at t = t m/2. Con-sistency with the known solutions of the degenerate cases where e — » - 0 and a — 0 was also shown from the l i m i t of t r m and the l i m i t s of a l l the parameters i n the approximant. Examples using th is ref ined method of parabol ic phase approximation were worked out, and compared with the true numerical so lut ion obtained from the d i g i t a l computer, as 8 9 shown i n F i g . 3 . 9 and F i g . 3 .10. The approximations were very close to the true solutions but could be improved by further reducing the value of This suggests a project for future research, since this work has i l l u s t r a t e d the v a l i d i t y of the approach. F i n a l l y , the K-B approximations were also plotted i n F i g . 3 . 9 and F i g . 3 . 1 0 for comparison. It is obvious by inspection that the K-B approximations were not as good as the parabolic phase approximations for this type of equation because the phase retardation appeared i n the parabolic phase approximations and did not i n the K-B approximations. In the case where e ^ .1, both the parabolic phase approximation and the K-B approximation f a i l to y i e l d ac-ceptable results because the solution i s no longer o s c i l -l a t o r y . Therefore, an e n t i r e l y d i f f e r e n t method was deve-loped. The solution x (t) to the complementary l i n e a r equation was f i r s t computed, and a correction term z(t) was then defined by z(t) = x c ( t ) - x(t) where x(t) was the solution to the non-linear equation. Thus, the problem of approximating x(t) was reduced to approximating z ( t ) . From various numerical examples, i t was suggested that z(t) could be approximated by ,n - t z(t) = g t e where g and n are constants depending on e and \i0 Using F i g . 3 . 1 2 and F i g . 3 . 1 3 , which were contour diagrams of n i n the £ p - ( i plane, n was empirically determined to be given by 90 n = 1.40 + 0.5 log 2 - ^ 2 — . 6 2.3+0.472uu - E Plots of log-^Q g against l o g ^ u for constant values of E were used, and g was then found to be given by the empirical formula l o g l 0 g = 0,794 l o g 1 0 a - [o.325+0.15(e~l)] Hence we were able to compute z'(t) when e and a were speci-f i e d , and the approximation to the solution x(t) was x(t) = x c("k) - 8 e~*^ • Consistency of this approximant with the known solution of the degenerate case where a »—0 was also shown. An example was worked out to i l l u s t r a t e this cor-rection term method, and the result was compared with the numerical solution i n F i g . 3.20. V i t h the aid of a d i g i t a l computer, the Ritz method i n conjunction with i n i t i a l con-d i t i o n matching was used to obtain another approximate solution which was also shown i n F i g . 3.20. When accuracy, e f f o r t , and p r a c t i c a b i l i t y were considered, the correction term method was much preferred. F i n a l l y , we see that by using these two approximating methods, values of a up to 10 and s and high as 2 may be accepted. Therefore, unlike a l l the c l a s s i c a l methods, they are good for f a i r l y gross n o n - l i n e a r i t i e s . The essential difference between these methods and c l a s s i c a l methods i s t h e i r d i r e c t approach i n attacking systems which are not qua.si—linear. In ^ conclusion, both the parabolic phase approxi-mation and the correction term approximation have strong 91 potential i n approximating the solutions to non-linear equations with f a i r l y large non-linearities whose characteristics can be represented by a d d cubics such as the flux-current r e l a t i o n of a saturating indicator, or the force—displacement r e l a t i o n of a hard spring. \ 4. CONCLUSION 92 As stated i n the Introduction, the purpose of this work has been to f i n d a dir e c t method of approximating the solution of the non-linear d i f f e r e n t i a l equations of the type x + 2ex + F(x) =0 where F(x) i s , or may be approximated by, an odd cubic with positive c o e f f i c i e n t s * Without loss of generality, we have studied i n d e t a i l the equation x + 2ex + x + fxx3 = 0, x(0) = 1, x(0) = 0, and then two methods have been developed to approximate the solution, according to whether e < 1, or e ^ 1. In the cake where e < 1, the parabolic phase approxi-mation was developed* The approximant was f i r s t derived to be of the following form for 0 < t ^ t x(t) = e ~ p t cos (o^t + « 2 t 2 ) for t > t x(t) = P e " e t cos (Vl - e 2 t + 0 q) where a l l the parameters were determined i n terms of e and u« Then, a refinement of the method changed the approximant to for 0 < t ^ t /2 x(t) = A e ~ p t cos (tt) + t + tt$9t2) in o o _L for t > t /2 x(t) = Pe" e t cos (Jl - e 2 t + 0 q ) , which yielded better r e s u l t s . In the case where e ^ 1, the solution was approximated by subtracting a correction term z(t) from the solution x c ( t ) of the complementary linear equation. The correction term z(t) v 93 was of the form z(t) = g t n e _ t , where g and n were computed from formulae involving e and a only. Therefore, x(t) = x (t) - g t n e~* . c Since the values of e and a are not limited to small values, we have found a dir e c t method of approximating the solution without resorting to quasi-linearization of the equation. The l i m i t of e i s s l i g h t l y above 2 and the l i m i t of a i s close to 10* These values are far too large to be handled by any c l a s s i c a l method* Although this method cannot handle e and a beyond their l i m i t s , i t has i l l u s t r a t e d the v a l i d i t y of the approach, and. further research along this l i n e i s encouraging. For example, similar methods may be developed for more general types of grossly non-linear equations. F i n a l l y , the goal of finding d i r e c t l y an approximate solution to the type of grossly non-linear equation has been achieved, and valuable insight into the free response of many engineering systems with odd cubic non-linear c h a r a c t e r i s t i c s , such as the hard spring and the saturating inductor, can readily be obtained. APPENDIX A ON COMPUTATION 94 The 4-th order Runge-Kutta method was used to ob-tai n the numerical solutions used throughout this work. In our case, the formulae for the computation of x and x are as follows: 9 2 x ,, = x + hx + l/6.' h (k, + k«+ k.,) n-f-1 n n 1 2 3 * n + l = *n + 1 / 6 h < k l + 2 k 2 + 2 k 3 + k4> * 3 where k^ = -(2ex n + x n + ax n ) k 2 = ^-[2e(xn+ 1/2 h kj_) + (x n+ | h + a(x n + \ k 3 = - [ 2 e ( x n + | hk 2) + ( x n + | hx n + I h 2 k x) + a ( x n + \ h x n + I h 2 k l ) 3 ] k 4 = - [ 2 e ( x n + hk 3) + ( x n + h x n + \ h 2k 2) + ^ ( x n + h * n + 2 h 2 k 2 ) 3 ] and h — t , ~, — t « n+1 n The University of B r i t i s h Columbia's IBM 1620 computer was used to carry out the computations to eight s i g n i f i c a n t figures* The program was written i n Fortran I I . Since the error of this method i s of the order of 5 h j, and h = 0.2 was used, the error expected was of the order of (0.2) , or 0.0003, which i s negligible for a l l p r a c t i c a l purposes« APPENDIX B THE KRILOFF AND BOGOLIUBOFF APPROXIMATION 95 The Kryloff-Bogoliuboff, or K-B method i s con-cerned with the transient solution of equations of the type x + co2x + k g(x,x) = 0 where g(x,x) i s arbitrary and K i s smallo An approximate (3) solution i s developed by Kryloff and Bogoliuboff and is e s s e n t i a l l y a va r i a t i o n of parameters technique• The solution i s assumed to have the form x(t) = A(t) cos 0(t) (B.l) where 0(t) = cot + 0(t) . Di f f e r e n t i a t i n g once, we have x(t) = A(t) cosQ(t) - A(t)[co+0(t) ] sin©(t). Since we have introduced one more variable, we can impose a constraint such as A(t) cosO(t) - A(t) 0(t) sinO(t) = 0 (B.2) so that x(t) = -A(t) co sinO(t), (B*3) and therefore x(t) = -A(t) co sinO(t) - coA(t) O(t) cos©(t). Substituting x ( t ) , x(t) and x(t) into the o r i g i n a l equation, we obtain -A(t) co sin©(t) - coA(t) 0(t) cosO(t) + K g(x,x) = 0 (B.4) Solving equations (B.2) and (B„4) simultaneously, we have A(t) = K sinO(t) (B.5) 96 and A(t) 0(t) = K cosO(t) . (B.6) The approximation i s made by averaging (B„5) and (B.6) over one period of o s c i l l a t i o n , assuming that A(t) i s constant over this period and can be taken out from under the integral sign. This means that i f A(t) i s slowly varying, the approximation i s a good one. In our case, g(x,x) = 2ex + (xx3 3 3 = -2eA(t)<osinO(t) + uA cos O(t), tt = I , K = 1. Substituting i n (B»6), and averaging over a period of 2TZ, we obtain A(t) = -cA(t), and 0(t) = ^ f 2 - ^ - . But A(0) = 1, from our framework of i n i t i a l conditions, and since A(t) i s assumed to be constant over the period of o s c i l l a t i o n , A(t) =* A(0) = 1. Therefore, we have A(t) = e ~ e \ and 0(t) = ^ . F i n a l l y , (B.l) becomes x(t) = e ~ e t cos (1 + 2g.) t, (B.7) 97 which i s the approximant used i n this work as a comparison to the methods developed i n Chapter 3. Note that the K-B method f a i l s to y i e l d acceptable results i n the case where e and a are not small, because i n such cases, A(t) does vary considerably over one period of o s c i l l a t i o n and the assumption required i n the averaging pro-cedure i s not a reasonable one, However, i f we do not make this assumption removing A(t) from under the integral signs, the integrations become very d i f f i c u l t , i f not impossible, to handle. Thus, we may not expect good results from equation (B.7) when e and a are r e l a t i v e l y large, as already i l l u s t r a t e d by various examples. i 98 APPENDIX C A MEASURE OF CLOSENESS BETWEEN LINEAR AND NON-LINEAR ISOCLINES The method of isocli n e s i s often used to construct phase-plane diagrams of 2-nd order d i f f e r e n t i a l equations. Con-sider the equation x + 2ex + x + ax = 0 ( C l ) and i t s complementary linear equation x + 2ex + x = 0 . (C.2) The i s o c l i n e s for (C«l) and (C.2) are respectively given by a n d ' y = = - f • ( 0 - 4 ) where y = x dx m = TJ^- = Slope of trajectory, and a = m + 2e. If these two sets of iso c l i n e s are close to each other i n the phase-plane, then equation (C.2) i s a good approximation to equation ( C . l ) . This means that the non-linear effect i n equation (Col) may be neglectedo In order to obtain a measure of closeness between the two sets of i s o c l i n e s , F i g . 3ol has been constructed, part of which i s reproduced i n F i g . C.I for convenient reference. In F i g . C . l , the c i r c u l a r arc of radius R intersects the lin e a r and non-linear i s o c l i n e s , for the same slope m, at points P and S respectively. The angle subtended by 99 y=x F i g . C.l A Measure of Closeness between Linear and Non—linear Isoclines of the Same Slope m the arc PS i s therefore a measure of the closeness between the i s o c l i n e s . However, because the co-ordinates of the point S are d i f f i c u l t to obtain, the point Q i s chosen instead, by dropping a v e r t i c a l l i n e from P to meet the non-linear i s o c l i n e . Then, the angle S© can be regarded as a measure of the closeness between the i s o c l i n e s . From equations (C.3) and (C.4), the co-ordinates of the point P are x = R cos 0 P v _ -R cosQ _ a and those of the point Q are x = R cos© Q 100 --R cos© /, , D2 2 A \ yg = — ( l + [XR COS ©) Therefore, we have tan© = — = — , x a P 2 2 and tan(© + So) = i f l = ^ " ^ R c o s Q i D2 2 i ( i + — a - ) . 1 + a Hence O© i s given by tan 8 © = tan [(© + 8 © ) - 8©] = . (c.5) (1+a 2) 2 + uR 2a 2 This shows that the angle 8© i s a function of R and a . As a varies, the angle 8© varies, and i t s maximum value, for a con-stant R, i s given by X ~ PfcanS©! = 0 3 a J which can be reduced to a 4 - (uU 2 + , 2)a 2 - 3 = 0, Therefore, we have 2 a = ~(2 + uR 2) + y | (2+uR 2) 2 + 3 . (Co6) This i s , then, the value of a 2 that w i l l give the maximum S© for a given R. From equations (C.5) and (C.6), therefore, the maximum value of §Q, i . e . (§©) m a x» i s dependent on uR2 only, o and uR becomes a measure of the closeness between the linear and non-linear i s o c l i n e s . For example, i f uR =0.2, (OO) Illct.A-i s approximately 3 , or 0.05 radian. Thus, for uR = 0.2, the two is o c l i n e s are almost coincident, suggesting that the effect of the non-linear term may be neglected at this point. F i n a l l y , since for equation (C.2), x(t) = 1 e " e t cos(Jl - e2 t - 0 ) JTT72 where 0 = tan" 1 g , we have A(t) = 1 e " e t x + x = . 11 + e sin (2 v l - e t Jl - eZ Thus, R o s c i l l a t e s about A(t) with a smaller and smaller ampli-tude as t increases. In other words, A(t) i s a very good 2 approximation to R, and therefore, uA (t) i s also a measure of the closeness between the two sets of i s o c l i n e s . This has enabled us to establish the angle c r i t e r i o n i n the develop-ment of the parabolic phase approximation. 102 REFERENCES 1. Cunningham, V.J*, "Introduction to Non-linear Analysis", McGraw-Hill Book Co., New York, 1958, pp. 121-170, 2. McLachlan, N.Y., "Ordinary Non-linear D i f f e r e n t i a l Equations in Engineering and Physical Science", Oxford, at the Clarendon Press, Second Edition, 1958, pp. 91-92. 3. K r y l o f f , N. and Bogoliuboff, N., "Introduction to Non-linear Mechanics", Princeton University Press, Princeton, 1947. 4 . Cunningham, W.J., op. c i t . , pp. 7 5 - 7 8 . 5 . Cunningham, V.J., op. c i t . , pp. 3 2 - 3 6 , 6 . Cunningham, V.J., op. c i t . , p. 1 6 8 . 7 . Cunningham, V. J. , op. c i t o , p. 1 5 7 o 8 . Cunningham, V.J., op. c i t . , pp. 1 5 1 - 1 5 4 . 9. Soudack, A.C., "Jacobian E l l i p t i c and Other Functions as Approximate Solutions to a Class of Grossly Non-linear D i f f e r e n t i a l Equations", Technical Report No* 2054-1, A p r i l 24, 1961, Stanford Electronics Laboratory, pp. 73-76. 10. Jahnke, E. and Emde, F., "Tables of Functions", Dover Publications, Fourth Edition, New York, .1945, p. 95. 11. Soudack, A.C., op. c i t . , pp. 75-76. 12. Cunningham, V.J*, op. c i t . , p. 123, 13. Tuttle, D.F., Unpublished notes on Nonlinear Analysis, Stanford University, 1960. 14. Antosiewicz, H.A* and Gautschi, V., "Numerical Methods i n Ordinary D i f f e r e n t i a l Equations", A Survey of Numerical Analysis, edited by J. Todd, McGraw H i l l Book Co. Inc., New York, 1962, p. 321.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Approximations to the free response of a damped non-linear...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Approximations to the free response of a damped non-linear system Chan, Paul Tsang-Leung 1965
pdf
Page Metadata
Item Metadata
Title | Approximations to the free response of a damped non-linear system |
Creator |
Chan, Paul Tsang-Leung |
Publisher | University of British Columbia |
Date Issued | 1965 |
Description | In the study of many engineering systems involving nonlinear elements such as a saturating inductor in an electrical circuit or a hard spring in a mechanical system, we face the problem of solving the equation ẍ + 2εẋ + x + μx³ = 0 which does not have an exact analytical solution,. Because a consistent framework is desirable in the course of the study, we can assume that the initial conditions are x(0) = 1 and ẋ(0) = 0 without loss of generality. This equation is studied in detail by using numerical solutions obtained from a digital computer. When ε and μ are small, classical methods such as the method of variation of parameters and averaging methods based on residuals provide analytical approximations to the equation and enable the engineer to gain useful insight into the system. However, when ε and μ are not small, these classical methods fail to yield acceptable results because they are all based on the assumption that the equation is quasi-linear. Therefore, two new analytical methods, namely: the parabolic phase approximation and the correction term approximation, are developed according to whether ε < 1 or ε ≥1, and are proven to be applicable for values of ε and μ far beyond the limit of classical methods. |
Subject |
Differential equations, Nonlinear |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-09-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0302279 |
URI | http://hdl.handle.net/2429/37525 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1965_A7 C5.pdf [ 4.33MB ]
- Metadata
- JSON: 831-1.0302279.json
- JSON-LD: 831-1.0302279-ld.json
- RDF/XML (Pretty): 831-1.0302279-rdf.xml
- RDF/JSON: 831-1.0302279-rdf.json
- Turtle: 831-1.0302279-turtle.txt
- N-Triples: 831-1.0302279-rdf-ntriples.txt
- Original Record: 831-1.0302279-source.json
- Full Text
- 831-1.0302279-fulltext.txt
- Citation
- 831-1.0302279.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0302279/manifest