UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Steady-state analysis of nonlinear networks. Agnew, David George 1971

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1971_A7 A45.pdf [ 2.89MB ]
Metadata
JSON: 831-1.0102065.json
JSON-LD: 831-1.0102065-ld.json
RDF/XML (Pretty): 831-1.0102065-rdf.xml
RDF/JSON: 831-1.0102065-rdf.json
Turtle: 831-1.0102065-turtle.txt
N-Triples: 831-1.0102065-rdf-ntriples.txt
Original Record: 831-1.0102065-source.json
Full Text
831-1.0102065-fulltext.txt
Citation
831-1.0102065.ris

Full Text

STEADY-STATE ANALYSIS OF NONLINEAR NETWORKS by DAVID GEORGE AGNEW B.Sc, University of Calgary, 1969 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE We accept this thesis as conforming to the required standard Research Supervisor Members of Committee Head of Department Members of the Department of E l e c t r i c a l Engineering THE UNIVERSITY OF BRITISH COLUMBIA March, 1971 In p r e s e n t i n g t h i s t h e s i s in p a r t i a l f u l f i l m e n t o f the r e q u i r e m e n t s f o r an advanced degree at the U n i v e r s i t y o f B r i t i s h Co lumb i a , I a g ree t h a t the 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 agree 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 purposes may be g r a n t e d by the Head o f my Department o r by h i s r e p r e s e n t a t i v e s . It i s u n d e r s t o o d tha 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 not be a l l o w e d w i thou t my w r i t t e n p e r m i s s i o n . The U n i v e r s i t y o f B r i t i s h Co lumbia Vancouver 8, Canada Depa rtment ABSTRACT In t h i s thesis the use of the d i g i t a l computer i n determining the steady-state response of nonlinear networks i s inve s t i g a t e d . Two classes of such networks, forced o s c i l l a t i n g networks and s e l f - o s c i l l a t i n g networks, are considered. For the f i r s t c l a s s , a method of obtaining the network equations i s developed, and an algorithm f o r c a l c u l a t i n g the numerical s o l u t i o n i s obtained. For the second c l a s s , the development i s r e s t r i c t e d to those networks which can be described by x + x = g(x, x) and the previous algorithm i s extended to solve t h i s problem. For both cl a s s e s , a s u b s t i t u t i o n i s developed which can expedite the analysis when the n o n l i n e a r i t i e s are exponential i n character. The eff e c t i v e n e s s of the methods i s demonstrated by examples. i i TABLE OF CONTENTS Page ABSTRACT i i LIST OF ILLUSTRATIONS iv LIST OF TABLES V ACKNOWLEDGEMENT v i I. INTRODUCTION 1 II. FORMULATION AND SOLUTION OF NONLINEAR NETWORK EQUATIONS 4 III. ON THE ANALYSIS OF OSCILLATORS 15 IV. EFFICIENT ANALYSIS FOR EXPONENTIAL NONLINEARITIES 24 V. PRACTICAL CONSIDERATIONS 28 VI. EXAMPLES " • 31 1. Nonlinear Inductor Example 32 2. Rayleigh Oscillator Example 37 3. A Transistor Amplifier Example 42 VII. CONCLUSION 50 APPENDIX A - Listing of Program for Example 1 52 APPENDIX B - Listing of Program for Example 2 55 APPENDIX C - Listing of Program for Example 3 58 REFERENCES 61 i i i LIST OF ILLUSTRATIONS Figure 2. 1 th The r Branch of a Network 2. 2 Graph and Corresponding A Matrix 2. 3 Branch Containing Nonlinear Element and Its Equivalent 5. 1 Example to Illustrate Network Behaviour with Respect to 5. 2 v as a Function of u for Circuit of Fig. 5.1. 5. 3 v as a Function of u for Modified Circuit of Fig. 5.1. 6. 1 Nonlinear Inductor Example 6. 2 Error Behaviour for Example One 6. 3 A Rayleigh Oscillator 6. 4 Error Behaviour of Rayleigh Oscillator 6. 5 Simple Transistor Amplifier 6. 6 Transistor Model 6. 7 Equivalent Circuit for the Transistor Amplifier 6. 8 Dependence of Gain on Amplitude 6. 9 Voltage Amplification as a Function of Input Voltage 6. 10 Harmonic Distortion as a Function of Input Voltage iv LIST OF TABLES Table Page 2.1 Nonlinear Dependent Sources 9 6.1 Results of Steady-State Analysis of R-L Network 35 6.2 Rayleigh O s c i l l a t o r Results 40 6.3 Behaviour of Transistor Amplifier 47 v ACKNOWLEDGEMENT I would l i k e to thank my supervisor, Dr. A. D. Moore, f o r h i s many contributions and assistance i n preparing the th e s i s . Also, Dr. A. C. Soudack's rapid reading of the thesis and h i s many suggestions are much appreciated. Thanks are also due to many others, but the assistance of Miss Linda Morris i n typing the t h e s i s , and Miss L. O'Leary and Mr. E. Tarnai i n proofreading i t , i s e s p e c i a l l y appreciated. F i n a n c i a l support from NRC i n the form of a 1967 Science Scholarship i s also g r a t e f u l l y acknowledged. v i I. INTRODUCTION In the past few years several general-purpose computer-aided algorithms f o r accurately p r e d i c t i n g the response of nonlinear networks have been developed (1, 2). However, i n s p i t e of the speed of modern computers, the time required f or such analyses i s f a r from n e g l i g i b l e , so more e f f i c i e n t algorithms continue to be developed. In t h i s t h e s i s , a new algorithm i s presented for obtaining the steady-state response of nonlinear networks, while avoiding some shortcomings of previous algorithms. To appreciate i t s advantages, i t i s necessary to consider some of the previous work on computer-aided network a n a l y s i s . There are now many general purpose programs a v a i l a b l e (3) f o r analyzing e l e c t r i c a l networks, mostly based on c l a s s i c a l analysis of l i n e a r systems. For l i n e a r networks, these programs give the d.c. response, the a.c. response, and the transient response. In contrast, c l a s s i c a l techniques f o r analysis of nonlinear systems have not been incorporated into general-purpose programs. Techniques such as the K r y l o f f - B o g o l i u b o f f method (4), and the perturbation method (5) are normally used only f o r f i r s t - or second-order almost-linear systems. Another common approach, the use of the p r i n c i p l e of harmonic balance, i s a general method, but as B a i l y (6) points out, "Though quite e f f e c t i v e i n p r i n c i p l e these methods s u f f e r , even i n r e l a t i v e l y simple problems, from gross algebraic complexities." The way these l i m i t a t i o n s have been overcome i n general purpose programs i s by use of a s t r i c t l y computer-oriented method: the numerical i n t e g r a t i o n of the nonlinear network d i f f e r e n t i a l equations. This approach to computer-aided nonlinear analysis has come to be almost u n i v e r s a l . B a i l y , r e f e r r i n g to algorithms for computer-aided analysis of nonlinar networks, says (7) 2 "The more general methods function v i a the numerical i n t e g r a t i o n of the network d i f f e r e n t i a l equations and the analysis r e s u l t s are usually presented as time pl o t s of the r e s u l t i n g network waveforms." As a r e s u l t , most general methods for computer-aided analysis of nonlinear networks are intended to provide the transient s o l u t i o n for the network. These approaches have some disadvantages. I f the steady-state response i s desired, i t i s frequently necessary to integrate f o r many periods of the e x c i t a t i o n before transients are n e g l i g i b l e . Furthermore, the equations of many networks contain widely separated time constants, and the shortest of these may l i m i t the s i z e of the time step used for numerical i n t e g r a t i o n , whereas only the longer ones may s i g n i f i c a n t l y influence the network behaviour. This drawback, while s i g n i f i c a n t i n analyzing l i n e a r networks, may cause more d i f f i c u l t y i n the nonlinear case, because terms which are e f f e c t i v e l y the time "constants" of the system may vary over a wide range. These l i m i t a t i o n s h a v e " r e s t r i c t e d the use of computer analysis of nonlinear networks. Three recent papers have reported on d i r e c t methods of obtaining the steady-state response of nonlinear networks, thereby avoiding the above problems. In 1968 B a i l y (7) developed a general approach based on the p r i n c i p l e of harmonic balance. He avoided the problems of algebraic complexity by using methods f o r multi-dimensional function minimization to adjust the parameters of the p e r i o d i c s o l u t i o n so as to minimize the res i d u a l s involved i n the harmonic balance method. In 1969 s i m i l a r r e s u l t s were published by N e i l l (8), and by Heywood and Moore (9), of i t e r a t i v e s o l u t i o n s to an i n t e g r a l formulation of the network equations. A l l these approaches have c e r t a i n features i n common. As opposed to i n t e g r a t i o n methods, which, once a step s i z e i s chosen, y i e l d a complete s o l u t i o n of f i x e d 3 accuracy a f t e r a given number of steps, the algorithms mentioned y i e l d a steady-state s o l u t i o n at each s t e p — t h e i t e r a t i o n s serve only to reduce the e r r o r . This has both advantages and disadvantages: while the s o l u t i o n process can be continued i n d e f i n i t e l y , and w i l l continue to improve the s o l u t i o n at each step making high accuracy p o s s i b l e , the improvement at each step may be so small that the time required to achieve reasonable accuracy may be excessive. C e r t a i n l y there i s s t i l l room for improvement. This thesis presents another general method for obtaining the > steady-state response of nonlinear networks. However, l i k e time-domain i n t e g r a t i o n , the method has associated with i t a step s i z e , which i s chosen i n advance, and the s o l u t i o n i s then obtained a f t e r a given number of steps. The s o l u t i o n has f i x e d error and there i s no p r o v i s i o n f or improving a s o l u t i o n already obtained. However, te s t problems t r i e d have shown that r e s u l t s obtained have reasonable accuracy, and require l e s s time than other methods to achieve the same accuracy. With that much j u s t i f i c a t i o n f o r the work, we s h a l l now b r i e f l y consider the content of the following chapters. In chapter I I , nodal v a r i a b l e s are used to develop a general method f o r w r i t i n g i m p l i c i t equations describing nonlinear networks. An algorithm f o r the numerical s o l u t i o n i s developed, i n v o l v i n g use of an extra parameter. This same technique i s then used i n chapter III as a basis f o r a general method for analyzing a class of o s c i l l a t o r s . A p a r t i c u l a r modi-f i c a t i o n to expedite analysis of exponential n o n l i n e a r i t i e s . i s the topi c of chapter IV. F i n a l l y , chapters V and VI contain, r e s p e c t i v e l y , some examples and the conclusion. 4 I I . FORMULATION AND SOLUTION OF NONLINEAR NETWORK EQUATIONS There are two parts to any general purpose algorithm f o r computer-aided network a n a l y s i s . The f i r s t involves a procedure for w r i t i n g equations describing the network. The second i s the numerical s o l u t i o n of the equations. Our algorithm for obtaining the steady-state response of non-l i n e a r networks w i l l be presented i n that order. Several such methods e x i s t , such as the nodal, mesh, and various hybrid d e s c r i p t i o n s , and t h i s chapter could have been developed i n terms of any of these; somewhat a r b i t r a r i l y , the nodal formulation was chosen as a ba s i s . Use i s made of Kron's standard branch convention, i n which the branch v a r i a b l e s are defined as i n F i g . 2.1. For the remainder of the chapter, a l l these v a r i a b l e s w i l l be frequency domain q u a n t i t i e s . Matrix methods of formulating the network problem are necessary. 6 The r th Branch of a Network F i g . 2.1 5 t i l The r branch contains an admittance Y , i n s e r i e s with a r r voltage source E and i n p a r a l l e l with a current source I , e i t h e r of r r . which may be zero. From Kirchhoff's laws, J = I + i (2.1) r r r V = E + e (2.2) r r r The branch current-voltage r e l a t i o n s h i p i s J = Y V (2.3) r r r r These r e l a t i o n s h i p s may be w r i t t e n f o r the whole network as J = I + 1 (2-4) V = E + e (2.5) and J = Y V (2.6) where Y i s a diagonal matrix with e n t r i e s Y » the self-admittances of the corresponding branches. L i n e a r l y dependent sources-may be accommodated by off-diagonal terms i n Y_. The other matrix needed i s the incidence matrix of the network, which w i l l be designated as A. The A matrix contains a column f o r every node but one (the reference node) and a row f o r every branch i n the net-work. Any given entry i s -1 i f the corresponding branch i s connected to and oriented toxjards the corresponding node; or +1 i f connected to but oriented away; or zero otherwise. F i g . 2.2 i s presented as an i l l u s t r a t i o n . 6 Reference IV Node I £/// 03 y Node 1 2 3 ; 0 ' 0 i 0 0 / -/ 0 i -i 0 0 ; -i A = 0 i -1 0 0 j 0 0 / 0 1 0 0 i 0 Graph and Corresponding A matrix F i g . 2.2 With the A matrix, two a d d i t i o n a l r e l a t i o n s h i p s may be wr i t t e n . The branch voltages (e) and node to reference voltages (e') are r e l a t e d by e = A e_' • (2.7) and, from Kirchhoff's current law, the branch currents are r e l a t e d by Afc 1 = 0 (2.8) Now s u b s t i t u t i n g (2.4) and (2.5) into (2.6), I+±= Y(E + e) (2.9) The s o l u t i o n i s as follows: I - Y E + i = Y e (2.10) i^(L ~ 1 _.) = Afc Y e^  , (2.11) At(I ~ 11) = At I A _' (2.12) e' = (Afc Y A ) " 1 ^ ( I - Y E) (2.13) Using (2.7) and (2.10), e_ and may be obtained. 7 A network containing nonlinear elements cannot be characterized i n the same way as a l i n e a r network. The A matrix i s w r i t t e n i n the same way as before, but the nonlinear elements cannot be described by s e l f -admittance or transadmittance. Consequently, when these elements occur the corresponding entries i n the Y_ matrix are not defined i n the usual fashion. Instead, the branches involved are put i n another form. the nonlinear element i s mathematically characterized i n such a way that the current through i t or voltage across i t can be ca l c u l a t e d from the s o l u t i o n . Then the element i s replaced by a current or voltage source (whichever i s known) whose value i s added to the e x c i t a t i o n already present. The equivalent sources representing nonlinear elements are indic a t e d by an a d d i t i o n a l subscript n, as i n F i g . 2.3. Actual Branch Equivalent branch - case 1. Equivalent branch -Assume the s o l u t i o n to the problem i s known. Assume also that case 2 Q 9 Q 6 6 6 Branch Containing Nonlinear Element and Its Equivalent F i g . 2.3 •8 The corresponding e n t r i e s i n Y can now be made. In case one, Y =0. In case two, the equivalent branch conductance i s i n f i n i t e , so a small impedance must be inse r t e d i n serie s with E and the corresponding admittance used. nr r o In doing so, one must be c a r e f u l not to perturb the network so far that the s o l u t i o n becomes u n r e a l i s t i c . th Consider the following example. Suppose the r branch consists of a nonlinear capacitor. At t h i s point we introduce the Fourier transform, and i t s inverse,3*^, where $t(f(t)) = F O ) = ^  r f ( t ) e j U t d t (2.14) J and ^ " 1 ( F ( c o ) ) = / 0°F(„)e j t b td„ (2.15) —OO Let the incremental capacitance, C , be a function of the instantaneous r ' nr' voltage,J5» ^{e^} , across the capacitor. The frequency domain current through the element i s therefore d t s T ^ e )] I =^{C [{* X ( e ) ] . T T — } , 0 nr nr r J dt J (2.16) Because time domain d i f f e r e n t i a t i o n i s avoided i n general (being prone to er r o r and r e q u i r i n g many numerical operations), frequency domain convolution i s s u b stituted. That i s , the transform current through the capacitor i s I n r = 5 i { C n r [ 3 r l ( e r ) ] } *{jtoer> . (2.17) where * denotes convolution. The equivalent current source, I , i s therefore nr I n r = ^ C n r [ ^ _ 1 ( e r ) ] } * {jcoer> (2.18) For purposes of numerical evaluation we require that i t be pos s i b l e to determine uniquely the value of C [??^(e ) ] . nr r 9 S i m i l a r l y , for an inductor whose incremental inductance i s a function of the current passing through i t , an equivalent voltage source i s s u b s t i t u t e d v ( 2 . 1 9 ) E n r = " ^ L n r [ > ( i r ) ] } * { j a 3V r - l , Again, i t must be p o s s i b l e to determine uniquely LiD* (i^_)}. This does not preclude h y s t e r e s i s c h a r a c t e r i s t i c s , but i f such a multiple-valued function i s used, there must be s u f f i c i e n t information a v a i l a b l e to choose the correct value at the time the function i s evaluated. Following the same argument, the current through a voltage and current dependent nonlinear conductance i s I = -ViC [t?i1('e ) , 0 ^ 1 ( i )]} * e nr nr r r r Nonlinear dependent sources are also allowed—they are summarized i n Table 2.1. (2.20) Table 2.1 Nonlinear Dependent Sources SOURCE Voltage-controlled voltage source Current-controlled voltage source Voltage-controlled current source Current-controlled current source V-I RELATIONSHIP ,e = D J { G [ ^ ( e )]} * e . r r q q (hybrid g) „-l, E = &{R [ 5 i x ( i )]} * i r m q q (mutual transresistance) \'^m[^\eqm * e q (mutual transconductance) I r =3f{Hh[3! ( i q ) ] } * ^ (hybrid h) 10 Eq. (2.13) can now be modified to account f o r the equivalent sources replacing the nonlinear elements. Since e^ , and i_ are functions of e 1 , E_, and I_, we may write e' = (At Y A ) ~ 1 A t { ( I - Y E) + L (e', E, I) - Y E (e', E, I)} (2.21) — — — — —n — — n — — — Now, (2.21) i s an i m p l i c i t equation which has no d i r e c t s o l u t i o n as i t s t a n d s — a common problem, unfortunately, with nonlinear equations. However, by m u l t i p l y i n g the nonlinear terms by a s u i t a b l e f a c t o r , a s o l u t i o n i s made po s s i b l e . We s h a l l adopt the notation f(y) f o r th i s f a c t o r . In the development which follows, i t w i l l be seen that y i s the independent parameter upon which the s o l u t i o n at each step depends. The fac t o r f(y) i s introduced i n t o (2.21) i n the following manner: e' = (Afc Y A ) " 1 A f c{(I - Y E) + f(y) [I (e', E, I) - Y E (e', E, I)]} (2.22) — — — — — — — 0 — — — n — — — The r a t i o n a l e f or introducing f(y) i s t h i s : by s e t t i n g f(y) = 0, the equations are l i n e a r i z e d and so are e a s i l y solved. Then by examining the dependence of the expression of Eq. (2.22) on y, a s o l u t i o n can be obtained for a network which i s perturbed s l i g h t l y from the l i n e a r i z e d network (the network represented by (2.22) with f(y) = 0 ) , i n the d i r e c t i o n of the o r i g i n a l network (the network represented by (2.22) with f(y) = 1). Then, with t h i s new s o l u t i o n as a s t a r t i n g point, the v a r i a t i o n with respect to y i s again determined, and the s o l u t i o n i s advanced s l i g h t l y further toward that of the o r i g i n a l network. This same step i s then performed repeatedly f o r a pre-determined number of times, to y i e l d the desired s o l u t i o n . After examining the choice of f ( y ) , we w i l l consider this process i n more d e t a i l . 11 To e x p l o i t this method i n obtaining a s o l u t i o n , some r e s t r i c t i o n s on f(u) are necessary:it must be zero for some value of u(say y^) so that the equations are l i n e a r i z e d and may be solved by conventional means; i t must be unity at some value of u(say y ) so that the equations represent the problem to be solved; and i t must be continuously d i f f e r e n t i a b l e on the range 0$f(y)£l, because the d e r i v a t i v e i n t h i s range w i l l be used i n obtaining a s o l u t i o n . In the examples which have been t r i e d , f(u) = V (2.23) was used i n almost every case and has given good r e s u l t s . However, i n chapter V an example i s given of a case for which i t would be more des i r a b l e to use a d i f f e r e n t f ( u ) ; consequently, the more general notation w i l l be retained f o r the r e s t of t h i s chapter. The f i r s t step i n the s o l u t i o n i s to set y = y„ (so f(y) . = 0) 0 'y=y o and solve the system of l i n e a r equations so obtained. Then, to advance the s o l u t i o n , i t i s necessary to d i f f e r e n t i a t e (2.22) with respect to y, which y i e l d Eq. (2.24). d l (e' , E, I) dE (e', E, I) —n — ' —' —_ _ —a. — — — dy dy (2.24) This i s an i m p l i c i t d i f f e r e n t i a l equation which cannot be solved for a r b i t r a r y nonlinear r e l a t i o n s h i p s but which can be integrated numerically. The following development uses Euler i n t e g r a t i o n ; provided accuracy require-ments are not too severe, i t i s the easiest to explain and understand. Eq. (2.24) i s put into the form of a d i f f e r e n c e equation by m u l t i p l y i n g by du and rewriting the d i f f e r e n t i a l equation as a d i f f e r e n c e 12 equation, using the d i f f e r e n c e operator A, Ae' = (A'Y A ^ V ^ f (vr) [Al^Ce' , E, I) - YAE^e' , E, I)•]' + [I (e', E, I) - Y E (e', E, I ) ] A f ( y ) } (2.25) —n — — — n — — — A l t e r n a t i v e l y , i f <^^ t J— i s e a s i l y evaluated, then Af(u) may be replaced , df(y) . b y ~ d 7 " A y ' At t h i s point a step s i z e must be chosen. We w i l l permit i t to have only c e r t a i n values, given by y l " y o Ay = , k integer (2.26) While t h i s may seem needlessly r e s t r i c t i v e , the f l e x i b i l i t y a v a i l a b l e i n the choice of f(y) obviates the need to consider other values. The number k i s the number of steps which w i l l be required to obtain a s o l u t i o n . I f k i s small, the s o l u t i o n w i l l be obtained r a p i d l y , but the e r r o r i n approximating the d e r i v a t i v e by f i r s t d i fferences may be excessive. Conversely, i f k i s too l a r g e , cumulative numerical errors may cause o v e r a l l accuracy to d e t e r i o r a t e . I t appears at t h i s juncture that optimum values of k must be chosen on the basis of experience with the method; the examples i n chapter VI should provide an idea of the range i n which to t r y . st Once Ay i s chosen, the method i s straightforward. At the j + 1 step, AI^(e' , E_, I) may be formed by subtracting I^(e_', E_, I) at the j ^  step st from I^CfL' > E_> I) a t t ^ e J + ^ step: Al j = I ( e , j , E, I.) - I ( e ^ - 1 , E, I) (2.27) —n n — — — n — — — 13 st In t h i s way a l l the differences are formed. We may now write the j + 1 step i n the algorithm. Ae' j = ( A S A ) " 1 A t { f ( y ( ) + j i n ) ^ ^ 3 , E, I) - I n ( e ' j ~ \ E, I) - Y(E ( e ' j , E, I) - E ( e , j _ 1 , E, I ) ) ] n — — — —n — — — + tl^ C e ' 3 , I , I) -YE^ C e ' J E, I ) ] [ f ( y Q + jAy) - f ( v Q ' + (j - l)Ay)]} (2.28) = e ' j + Ae' j (2.29) A l t e r n a t i v e l y , the term Af(u) = f ( y Q + jAy) - f ( y Q + (j -l)Ay) (2.30) may be replaced by Af(.y) = | i ( y Q + jAy)Ay (2.31) i f more convenient. The index j runs from 0 to k-1.. With j=0, one might i - l 1-1 wonder about evaluating I (e' J , E, I) and E (e , E, I ) ; however, no —Q — — — —n — — — problems a r i s e because f ( y Q + jAy) | = f ( y Q ) = 0 (2.32) 1-1 1-1 so the values of I (e , E, I) and E (e' J ,E,I) are immaterial In t h i s case —n — — — —n — — — This completes the development of the mathematical basis of the method, but i t s implementation as a computer method must s t i l l be considered. Returning to the example of the nonlinear capacitor, we have I n r = D i { C n r [ ^ 1 ( e r ) ) ] 3 r 1 ( j o ) e r ) } - (2.33) Now e' i s a s e r i e s of the form r N e' = a n6(w) + E {a 6 (co - nw.) + a '6 (u + ntO } (2.34) r U n n O n 0 n=l where - denotes complex conjugate. 14 The c o e f f i c i e n t s are stored i n the computer as a set of complex numbers. To perform the operation (jtoe^) each a^ i s m u l t i p l i e d by j u with the corresponding value of ui s u b s t i t u t e d . The operation D* (e^) i s then performed using the inverse f a s t Fourier transform on the set of numbers a^, y i e l d i n g a time-domain sequence. The evaluation of {C["3* ~^(ep 1 • 3» 1 ( J w e r ) } must then be performed for each point i n the time sequence, y i e l d i n g a new time sequence. The d i s c r e t e Fourier transform operation;.is then per-formed on the sequence using the f a s t Fourier transform. Problems also a r i s e i n evaluating Y_ and (AST A)-"'" because they are frequency dependent. These problems are avoided by computing the two matrices at each frequency of i n t e r e s t before commencing the algorithm. Then at each step, a f t e r evaluating I^Ce', E_, I) and E^(e', E_, I) , Equation (2.2 8) i s solved at the frequency of each harmonic of i n t e r e s t . The problem of truncation of the Fourier s e r i e s should also be mentioned. Obviously, use of an i n f i n i t e s e r i e s i s not f e a s i b l e ; but truncation may introduce e r r o r s . Using fewer harmonics saves computer time and storage, but may increase the e r r o r . This unfortunate l i m i t a t i o n , having to trade speed and accuracy f o r storage, cannot be avoided. Even worse, the analyst must choose the number of harmonics to be evaluated at the outset, when he usually knows l e a s t about the behaviour of the problem. Numerical experimentation may be necessary. This completes the discussion on analysis of nonlinear networks with p e r i o d i c e x c i t a t i o n . S t a r t i n g with nodal analysis f o r l i n e a r networks, the equations for nonlinear networks have been developed. Because they cannot be solved d i r e c t l y , a numerical method has been proposed. This technique y i e l d s the steady-state response of nonlinear networks e f f i c i e n t l y . 15 I I I . ON THE ANALYSIS OF OSCILLATORS Discussion of a general purpose method for "analyzing nonlinear networks would not be complete without considering s e l f - o s c i l l a t i n g networks. Unfortunately, the development of the previous chapter using nodal v a r i a b l e s i s , as w i l l be shown, unsuitable as a general purpose method. An a l t e r n a t i v e form i s considered which, though not completely general, i s u s e f u l for many problems. A numerical s o l u t i o n w i l l be presented, based on the a n a l y s i s from the previous chapter, but in c l u d i n g necessary extensions. Let us f i r s t put the problem i n perspective. The previously mentioned papers by B a i l y (6), N e i l l (7), and Heywood and Moore (8) do not develop the analysis of o s c i l l a t o r s , although B a i l y does state that with an extension of h i s method i t should be possible to do so. The reason these authors have neglected this aspect may be due i n part to the increased complexity involved. For example, one d i f f i c u l t y common to the above methods and the method of chapter I I , i s that an i n i t i a l s o l u t i o n i s necessary, but the amplitude and frequency are o r i g i n a l l y unknown. Using the method of chapter I I , t h i s means the Y_ matrix cannot be evaluated. Another d i f f i c u l t y i s that the frequency may vary i n the process of c a l c u l a t i n g the s o l u t i o n , leading to an a d d i t i o n a l parameter i n the s o l u t i o n . Any approach to the problem must take these factors i n t o account. Because of d i f f i c u l t i e s a r i s i n g from the increased complexity, i t i s necessary to consider a more r e s t r i c t e d approach f o r o s c i l l a t o r s than that considered i n chapter II for forced networks. The problem which w i l l be treated i s the analysis of those o s c i l l a t i n g networks described by the time-domain equation 16 x + x = g(x, k) (3.1) In this equation, the effects of a l l active elements, dissipative elements, and nonlinear elements are included i n g(x, x) . The coefficients of x and x are normalized for convenience. Unfortunately not a l l s e l f - o s c i l l a t i n g networks readily reduce to this form; for example, the R-C phase shift oscillator does not. However an L-C o s c i l l a t o r w i l l usually reduce to this form, which gives credence to our claim that this method is sufficiently general for many problems. Equation (3.1) is not new; i n fact, several classical techniques exist for analyzing i t , such as variation of parameters (10), perturbation methods (11), and Kryloff's and Bogoliuboff's averaging method (4). But none of these are general because they assume either that g(x, x) remains small compared to x or x, or else that g(x, x) must be large compared to either x or x at some point during the period, so that a different type of oscillation occurs (as in relaxation os c i l l a t i o n s ) . Alternatively, the principle of harmonic balance may be used, but the d i f f i c u l t i e s involved, especially the complexity of the nonlinear algebraic relationships which result, are even more formidable than with forced networks because frequency is also a variable in the oscillator problem. The generality of the computer-aided method which follows results because i t surmounts these d i f f i c u l t i e s . As the f i r s t step in the development of this approach, Eq. 3.1 is transformed to the frequency domain. (1 - co2)X(co) =y-{g[0T1(X), ^ ( j i o X ) ] } (3.2) where OO X(u) =.jC{x(t)} = E a 6(a> - nu^) (3.3) n=-°° 17 The function f(y) i s again introduced: (1 - u 2)X(u)) = f(y) O ^ g f O^CX), ^ _ 1 ( j o ) X ) ] } (3.4) The requirement on f(y) are s i m i l a r to those i n chapter I I : there must e x i s t y^ and y^ such that and f ( t t 0 ) = 0 (3.5) f ( y x ) = 1 (3.6) and f(y) must be r e a l and continuously d i f f e r e n t i a b l e i n the range 0*f(y)<l (3.7) Again, f( y ) = y (3.8) was used i n the examples t r i e d , but the notation f(y) i s retained for gen e r a l i t y . The next step, d i f f e r e n t i a t i n g Eq. (3.4) with respect to y, i s analogous to w r i t i n g Eq. (2.18) i n the l a s t chapter; the r e s u l t i s (1 - a)2) g =C?{g[Di" 1(X), )y- 1(ja)X)]} ^ H l + f ( y ) ^ {^tg(^" 1(X),Oi~ 1(ja )X))]} (3.9) The p a r a l l e l development i s continued by m u l t i p l y i n g (3.9) by dy and rewriting the d i f f e r e n t i a l equation as a f i r s t d i f f e r e n c e approximation using the d i f f e r e n c e operator A. The notation AX' = (1 - to2)AX (3.10) i s also introduced, giving 18 AX' =0Hg[3 i l(X),^ 1 ( j o J X ) ] ) A f ( y ) + f ( y ) A {3&[g(y 1 ( X ) , D^ _ 1(jcoX))]} (3.11) At t h i s point, the algorithm must deviate from that of chapter I I , even though Eq. (3.11) i s analogous to (2.24) i n the previous chapter. The d i f f i c u l t y a r ises because each step may involve a frequency change, which precludes a simple a p p l i c a t i o n of (3.11). However, i t turns out that a l l the a d d i t i o n a l information necessary to apply (3.11) can be obtained by ensuring that, at each step, (3.4) i s s a t i s f i e d at the funda-mental frequency f o r whatever value y has at that step. The following developments are based on th i s premise. The s o l u t i o n to (3.4) with f ( y ) = f ( y Q ) = 0 i s v n = a i 6 < u ~ 1) + ^ Ha + 1) (3.12) f(p).-0 1 1 where the value of a^ i s not uniquely determined. Further r e l a t i o n s h i p s are necessary f o r t h i s purpose. Because the time o r i g i n for the s o l u t i o n i s a r b i t r a r y , we may choose i t to obtain a given phase f o r the fundamental; further on, when y ^ y^, the phases of the other harmonics w i l l follow accordingly. For convenience, we a r b i t r a r i l y make r e a l , so x ( t > | f ( p ) = c j = 2 cos t (3.13) Having made a^ r e a l , i t s magnitude must now be determined. Assume the o s c i l l a t o r i s i n the steady-state. Let 00 3Kg(x, x)} = E b n 6(a) - no ) (3.14) n=~oo Then (3.4) becomes 19 2 2 E {(1 - n co 1)a 6 (co - nto, ) } = f(y ) E b 6 (co - nto..) (3.15) l n x n 1 n=-<» n=-°° Now, because a^ i s r e a l and f(y) i s required to be r e a l , then f o r (3.15) to hold, b^ must also be r e a l . Consequently, at each step i n c l u d i n g the f i r s t , the c r i t e r i o n for choosing a, (assuming a l l other c o e f f i c i e n t s a 1 n fo r n 4 1 corresponding to the current value of y are known) i s that Im {b^} must be zero. I f a s u i t a b l e value of a^ i s not known i n i t i a l l y , a numerical search must be made, using t h i s c r i t e r i o n . The choice of a^ to make Im{b^} = 0 does not guarantee that Re(b^} = 0; i n general i t does not, re q u i r i n g another r e l a t i o n s h i p i n order to s a t i s f y (3.15). The p h y s i c a l basis f o r the next step follows because, i f Retb^} 4- 0, there i s r e v e r s i b l e energy t r a n s f e r between the resonance producing elements and the elements described by g(x, x ) , the e f f e c t of which i s to modify the frequency. This i s e a s i l y accounted for mathematically, as the frequency of the fundamental component has not been s p e c i f i e d as yet except when f(y) = 0. The necessary r e l a t i o n s h i p i s obtained by considering the components i n (3.16) at the fundamental frequency. (1 - to 2)a 1[6(co - to1) + 6(to + u±)]= f(y)b 1 f 6(to - u^) + 6(to + o^) ] (3,16) I t follows that 2 b l 1 - = f(y) - i (3.17) 1 a± / b l cox = J 1 - f(y) (3.18) E s s e n t i a l l y , the reason f o r r e s t r i c t i n g the o s c i l l a t o r a n alysis to the form of Eq. (3.1) i s that the frequency r e l a t i o n s h i p , Eq. (3.18) i s a quadratic, '20 and therefore has a d i r e c t s o l u t i o n . A step s i z e , Ay, i s again necessary, and must be chosen on the same basi s as i n chapter I I , namely, y l " y0 Ay - _ } £ integer (3.19) where k i s the number of steps which w i l l be required to obtain a s o l u t i o n . The choice of k i s somewhat a r b i t r a r y , although the same considerations as i n chapter II apply: i f k i s small, the s o l u t i o n i s obtained r a p i d l y , but the er r o r i n approximating the d i f f e r e n t i a l s by f i r s t d i fferences may be excessive; or i f k i s too l a r g e , cumulative numerical errors may cause over-a l l accuracy to d e t e r i o r a t e . I t is. unfortunate that concise rules f o r choosing k have not been found—we can only r e f e r the reader to the examples i n chapter VI for an idea of the range i n which to try and the corresponding accuracy which may be expected. (An e r r o r estimate -is given f o r d i f f e r e n t values of k for the f i r s t two examples). st In terms of Ay, Eq. (3.11) becomes, at the j + 1 step, AX' =3i{g[Dr 1 ( x j),^" 1(ja)X j)]}{f [ y 0 + jAy] - f [ y Q + (j - l)Ay]} + f ( y Q + jAy) mg&hxh?'1^^))] - ^ [ g ^ " 1 ( X j " 1 ) - S ^ C j t o X ^ 1 ) ) ] } (3.20) With t h i s r e s u l t a v a i l a b l e , we may now consider the d e t a i l s of the algorithm. The changes i n the c o e f f i c i e n t s of the harmonics of the fundamental are determined using Eq. (3.20) to obtain AX*. Now AX' i s a Fourier s e r i e s of the form OO Ax' = Z c 6(u - ncO (3.21) n 1 21 Relationship (3.10) i s used to c a l c u l a t e the new values of the c o e f f i c i e n t s at each step. The r e s u l t i s c a = a J + -n n 1 2 1 - nco^ (3.22). n = 0,2,3,4,... Before t h i s r e l a t i o n s h i p can be applied, however, the fundamental frequency, s t co^, at the j + 1 step must be found. Fortunately, we are now i n a p o s i t i o n to evaluate (3.18), / b l co1 =y4 - f(y) -± (3.18) because the change i n the product f ( y ) b ^ at each step i s given by Re(c^} i n (3.21). Let v ='f(y)b ]_ (3.23) The value of v at the outset i s s t At the j + 1 step, v = f ( y Q ) b 1 = 0 (3.24) V J + 1 = v j + Re{cJ} (3.25) from which we may determine i o . ^ + ^ : a l where the l a t e s t values of y and a^ are used. With the new value of frequency known, new values for a l l harmonics but the fundamental can be calculated using (3.22). The one question s t i l l open i s that of determining the new value of a^, the magnitude of the fundamental component of o s c i l l a t i o n , at each 22 step. We again invoke the requirement that at each step Im(b^} = 0. At the beginning of each step t h i s requirement i s s a t i s f i e d from the previous step. Consequently, Im{c^} represents the deviation of Im{b-^ } from zero which would occur at that step i f a^ were not changed. However, the change i s e a s i l y c a l c u l a t e d , because new values of a l l other s o l u t i o n parameters may be calculated without the new value of a^, so numerical d i f f e r e n t i a t i o n may be used to evaluate for the new s o l u t i o n as obtained up to this point. Then a Newton-Raphson cor r e c t i o n i s used to obtain the new value of a^: •4-1 • Im{c1 }/ „ a ^ - a ^ - y ^ f - t l ^ ) } (3.27) This completes the development of the method. To l i n k the various steps of the algorithm together, a short summary follows: 1. I n i t i a l i z e the values of o>, and the c o e f f i c i e n t s a . 1 n Choose k and evaluate (3.20). 2. Evaluate A X ' from (3.21). 3. Calculate the new frequency from (3.29). 4. Evaluate new values of the c o e f f i c i e n t s a (except a.) n 1 using (3.23). 5. Using the new s o l u t i o n as determined so f a r , c a l c u l a t e 8 (inKc,)} 6. Use (3.28) to c a l c u l a t e a new value f o r a^. 7. If j^k - 1, Increase j by one and go back to (2); otherwise stop. This completes the discussion on s e l f - o s c i l l a t i n g networks, except for the example problem, an analysis of a Rayleigh o s c i l l a t o r , i n chapter VI. In t h i s chapter, the problems of analyzing o s c i l l a t o r s have been considered. The analysis has been c a r r i e d out for a r e s t r i c t e d c l a s s of o s c i l l a t o r s , those described by Eq. (3.1). A computer aided method for th i s class has been developed. An algorithm f o r the Fourier c o e f f i c i e n t s i n the steady-state s o l u t i o n i s the end r e s u l t . 24 IV. EFFICIENT ANALYSIS FOR EXPONENTIAL NONLINEARITIES Because of the profusion of b i p o l a r semiconductor components now used i n e l e c t r o n i c c i r c u i t s , i t would seem that e f f i c i e n t methods f o r analyzing such c i r c u i t s i n the steady-state should be p a r t i c u l a r l y d e s i r a b l e . Using the methods i n chapters II and I I I , a p a r t i c u l a r s u b s t i t u t i o n i s possible which w i l l , when using an exponential model for such devices, expedite the analysis of such c i r c u i t s . We s h a l l discuss the cases to which t h i s a p p l i e s , and develop the s u b s t i t u t i o n involved. The s i m p l i f i c a t i o n i s most e f f e c t i v e when using the Ebers-Moll model (12) f o r the p-n j u n c t i o n t r a n s i s t o r , because a l l the nonlinear terms i n t h i s model are exponential terms. When th i s model i s used i n (2.17) or (3.2) terms of the form ^(exp [O1 "*"(...)] } r e s u l t . The advantages occur because the exponential of a Fourier s e r i e s can be w r i t t e n as another Fourier s e r i e s i n terms of Bessel functions, so the repeated, transformations are avoided. The necessary function i s e a s i l y derived from the generating function for Bessel's c o e f f i c i e n t s , I 1 e 2 a ( u - ^ = E J n(z)u> n (4.1) n=-°° Sub s t i t u t i n g z = j x (4.2) and w = j _ 1 e j G (4.3) gives ,x,.-.l j9 . -j0, Ijb e J - je J ) -e 2 = E j n J n ( j x ) e : ] n 6 (4.4) n = - o o 25 j(eJ + e J ) * 2 - S I ( x ) e j n W (4.5) n n=-»o or xcos 9 „ T / , j n Q r\ e = £ I (x)e J (4.6) n n = - o o where ^ n ( x ) i s Bessel's function f o r imaginary argument, sometimes c a l l e d the modified Bessel function of the f i r s t kind. To consider the a p p l i c a t i o n to this a n a l y s i s , assume we have a frequency domain s i g n a l c o n t r o l l i n g a nonlinear element with exponential c h a r a c t e r i s t i c . For ease of explanation the analysis w i l l be c a r r i e d out r e t a i n i n g only the d.c. component, fundamental, and second harmonic at each step; therefore S = I a 6(co - noO (4.7) 1 0 n 1 n=-2 The desired s i g n a l , the response of the exponential-type nonlinear element to S^, i s S^, where S 2 =D*{exp[3i" 1(S 1)]} (4.8) The time domain s i g n a l S,>(t) i s S_(t) = exp(a_ + 2a cos to t + 2a 1 s i n to11 ^ U X c J. i s x where + 2 a 2 c c o s 2co1t + 2 a 2 s s i n 2to 1t) (4.9) a. = Re{a.. } = Re{a } l c i -1 Is = ~Im{a^} = Im{a }^ a 2 c = Re{a^} = Re-{a2> a 2 s = - Im{a2} = Im{a2> 26 L e t t i n g -1 Is ,. , n . <() = tan (4.11) 1 a. l c and _^ a 2 *2 = t a n " i f (4.12) 2c y i e l d s S (t) ='exp{a Q + 2 |a 1|cos(a) 0t - c^) + 2 |a 2 | cos(o) Qt -<J>2)} = exp{a Q} exp{2 (a^ | cos((jOQt - <j>^ )} exp{2 | a2|cos (2to Qt - <j>2)} (4.13) Sub s t i t u t i n g (4.6) into (4.13) and again r e t a i n i n g only terms up to the second harmonic, 2 S 2 ( t ) = exp(a Q){ E 1^(2 \ | ) (cosnc^ - j s i n n<j>1) exp (jnu) Qt) } . n=-2 .1 { E I n ( 2 | a 2 | ) (cos n<i>2 - jsinn<j)2)exp(j2nu) t) } (4.14) n=-l The -frequency domain s i g n a l S 2 may now be w r i t t e n d i r e c t l y as 2 S 2 = e x p ( a Q ) { n = E 2 I n ( 2 | a 1 | ) (cos nc^ - j s i n n ^ ) 6 (to - nu^) } * { £ I n ( 2 | a 1 | ) ( c o s n<j)2 - j s i n n<j>2) 6 (to - nto ) } (4.15) n=-l The convolution i s e a s i l y performed i n the frequency domain because only impulses are involved. In implementing t h i s r e s u l t i n a computer program, many fu n c t i o n a l evaluations can be avoided by using the f a c t that I_ n(x) = I n ( x ) ( f o r n integer) (4.16) Also complex conjugates occur and since only a sign change i s involved, the e n t i r e number need not be re-evaluated. 27 Eq. (4.15) i s an i n d i c a t i o n of the way an exponential n o n l i n e a r i t y can be treated i n a c i r c u i t analysis program without recourse to the fast Fourier transform or i t s inverse. As Gold and Rader (13) point out, the fast Fourier transform i s more e f f i c i e n t than summing products of s e r i e s only when the s e r i e s are longer than a c e r t a i n length. Thus i t may be expected that the method of this chapter should have i t s greatest advantage when only few harmonics are of i n t e r e s t . 28 V. PRACTICAL CONSIDERATIONS Before commencing the examples, i t seems worthwhile to consider some p r a c t i c a l aspects, the use of the f a s t Fourier transform, the choice of the function f ( u ) , and accuracy requirements. The f a s t Fourier transform has allowed considerable reduction i n the time required for numerical d i s c r e t e time-domain to Fourier domain transformations. However, i n a context i n which i t i s used to transform a s o l u t i o n of unknown character, a caution concerning i t s use may save some computing time. The f a s t Fourier transform i s most e f f i c i e n t when the number of time points i s an integer power of two. Since a large proportion of computer time taken i s often for performing the transformations, then i f the f i r s t fourteen harmonics are considered important, time may a c t u a l l y be saved by c a l c u l a t i n g using the f i r s t s i xteen. If a l l that i s known i s that the number of harmonics which are s i g n i f i c a n t i s between eight and sixteen, then i t i s even more l o g i c a l to c a l c u l a t e using sixteen. The second point concerns the choice of f ( u ) . The point can probably best be i l l u s t r a t e d by a l i n e a r example which can be solved e a s i l y . Consider F i g . 5.1. Example to I l l u s t r a t e Network Behaviour With Respect to u F i g . 5.1 By i n s p e c t i o n , v = (5 1 + lOOy P l o t t i n g v as a "function of y gives F i g . 5.2, V A •M v as a Function of y for C i r c u i t of F i g . 5.1 F i g . 5.2. The maximum value of |_dv| occurs at y = 0 and i s given by dy of y. I dv i ' dy 1 max = 100 Now suppose = lOOy . Once again v i s p l o t t e d as a f u n c t i o n 0.2 0.4 0.6 0.8 1.0 v as Function of y f o r Modified C i r c u i t of F i g . 5.1 F i g . 5.3 30 or y d v i dy max than an order of magnitude. This example has demonstrated that a choice of f ( y ) other than f(y) = y may reduce the maximum value of the de r i v a t i v e s with respect to y which are required i n the analysis method. By reducing without increasing e r r o r . The l a s t point concerns accuracy. For some purposes, p a r t i c u l a r l y network optimization, i t i s convenient to be able to trade speed f o r accuracy. The methods which have been presented allow t h i s . Assuming the step s i z e i s not so small that cumulative numerical errors are s i g n i f i c a n t , then using a l a r g e r step s i z e w i l l decrease the computer time required, but w i l l also increase the er r o r . This permits a reasonably well-defined choice to be made. aid the reader i n using the nonlinear steady-state analysis method i n an optimum manner. A f a c t o r a f f e c t i n g the choice of the number of harmonics to evaluate has been presented. Then, by example, i t has been demonstrated that the choice f(y) = y i s not n e c e s s a r i l y optimum. And f i n a l l y one advantage of the method, i t s a b i l i t y to trade accuracy f o r speed, has been explained. the maximum value of these d e r i v a t i v e s , a la r g e r step s i z e may be used The three points considered i n th i s chapter were discussed to VI. EXAMPLES The effectiveness of the steady-state nonlinear network analysis technique w i l l be demonstrated by three examples chosen to i l l u s t r a t e the most important points i n previous chapters. The f i r s t example i s a non-l i n e a r inductor i n s e r i e s with a r e s i s t o r and voltage source, the second i s a Rayleigh o s c i l l a t o r , and the t h i r d i s a b i p o l a r t r a n s i s t o r a m p l i f i e r . The examples were solved using Fortran IV programs on an IBM 360/67. For each example, the equations are set up i n the required form and numerical and graphical r e s u l t s are given. L i s t i n g s of the programs used are i n the appendices. 32 1. Nonlinear Inductor Example The c i r c u i t to be analyzed i s given by F i g . 6.1. B a i l y (14) has analyzed the same problem (also using an IBM 360/67) and has given s u f f i c i e n t data to enable us to compare performance of the two methods. v(t) 1(E) = 6 cos TTt 7\+ i = X + X" e = X By ins p e c t i o n , l e t t i n g gives where Nonlinear Inductor Example F i g . 6.1 f(y) = y e = E - y l I = 3={X + X } n X = C* 1{jcoe} (6.1) (6.2) (6.3) (6.4) st The i n i t i a l value of e i s E. At the j +' 1 step, the d i f f e r e n c e equation i s 33 e j + 1 = e J + j ( A y ) ( I n J - 1 ^ X) + I n J ( A y ) (6.5) where I i s given by (6.2) and Ay = 1/k (6.6) where k i s the number of steps to be used i n c a l c u l a t i n g a s o l u t i o n . The index j goes from 0 to k - 1. A program incorporating t h i s algorithm i s given i n Appendix I and r e s u l t s follow i n Table 6.1. As a measure of the e r r o r , the c r i t e r i o n chosen by B a i l y i s used: e = 1/T / T ( X + I + A 3 - 6 cos t ) 2 dt (6.7) 0 where T i s one period of the fundamental, and X i s the cal c u l a t e d s o l u t i o n . The best r e s u l t s were obtained using k = 110 and carrying 16 harmonics (the l a r g e s t number for which the program was intended). They -4 are given i n Table 6.1. The corresponding value of e was 1.019 x 10 For comparison, B a i l y ' s r e s u l t s are also given, the corresponding values i n each case being given immediately below i n brackets. The corresponding e i s 1.627 x 10~1J:. The r e s u l t s i n F i g . 6.2 give a reasonable comparison between the methods. The curves f o r B a i l y ' s method are the only two given i n h i s thesis so there was no choice i n s e l e c t i n g them f or comparison. Curves B and D show that B a i l y ' s ultimate s o l u t i o n has l e s s e r r o r , which i s expected since h i s method operates by reducing error u n t i l i t can be reduced no further. In contrast, the best s o l u t i o n obtained herein has considerably more er r o r . I t i s also apparent that the simplest approach to increased accuracy, that i s , taking more steps (k) with a smaller step-s i z e (Ay) i s not dependable: i f too many steps are taken, cumulative 34 numerical errors overwhelm the i n t e g r a t i o n e r r o r , and accuracy d e t e r i o r a t e s . It i s the behaviour with respect to computing time, however, which demonstrates the value of the method. As F i g . 6.2 shows, a good estimate of the s o l u t i o n can be obtained, i f accuracy requirements are not too severe, with a time saving of one to two orders of magnitude over that required by B a i l y ' s method to achieve the same accuracy. This example was chosen because i t permitted a good comparison between the method developed i n chapter II and B a i l y ' s method, a computer-aided harmonic balance method. The conclusions which may be drawn from t h i s example are that B a i l y ' s method can lead to more accurate r e s u l t s , but i f accuracy requirements are not too severe, the new algorithm saves computing time. Table 6.1 Results of Steady-State Analysis of R-L Network (Baily's Results are Given in Brackets) HARMONIC 1 3 5 7 9 11 13 15 Sine Coefficients 1.137 (1.139) 4.852x10"2 (4.834xl0 - 2) 4.216xl0~3 (4.154xl0~3) 4.180xl0"4 {4.007xl0~4) 4.331xl0~5 (4.007xl0~5) 4.612xl0~6 (4.050xl0~6) 4..816xl0"7 (4.082xl0~7) 5.312xl0"8 (4.049xl0~8) Cosine Coefficients 0.9123 (0.9089) 5.651xl0"2 (5.627xl0 _ 2) 5.775xl0"3 (5.715xl0~3) 6.382xl0~4 (6.315xlO~4J 7.292xl0~5 (7.240xl0~5) 8.497xl0"6 (8.474xl0 _ 6) 9.428xl0~7 (1.004xl0~6) 8.095xl0~8 (1.199xl0~7) Co 36 Error(e) 10"' --10' •2 10~3 + w~5 + io~7 + E r r o r Behaviour f o r Example One - Parameter Method Carrying 1 4 t b Harmonics B - Parameter Method Carrying 1 1 6 t b Harmonics s t „nd C D s t rd B a i l y ' s Method Carrying 1 and 3 Harmonic B a i l y ' s Method Carrying 1 & l l t b Harmonics st ,rd .th Numbers.beside c i r c l e s i n d i c a t e number of steps (k) used to obtain point. 0.1 0.3 1.0 3.0 Tims (sec.) 37 2. Rayleigh O s c i l l a t o r Example This example i s intended to i l l u s t r a t e the application of the theory derived i n chapter I I I for o s c i l l a t o r analysis. The c i r c u i t to be analyzed i s given i n Fig. 6.3. L=JH. :C=IF. ln=L01e 3 _e 3 A Rayleigh O s c i l l a t o r Fig. 6.3 The d i f f e r e n t i a l equation follows by inspection: dt R L n 3 (6.7) Eq. (6.7) i s rewritten i n terms of x, the inductor f l u x linkage. Cx + f + f =1 R L n (6.8) Substituting component values, x + 0.01 x + x = 1.01 £ -Upon reducing (6.9) to the form of (3.1), there results .3 (6.9) 38 L e t t i n g f(u) = y (6.11) the equation i s reduced to the necessary form. 3 £ + x = y {x - ~ } (6.12) This l a s t form i s known as Rayleigh's Equation, and i s the basis f o r the computer program i n appendix B. A s u i t a b l e quadratic er r o r c r i t e r i o n i s defined to be e where m = ^ n e = Z {[Re(x ) - R e ( x ' ) ] 2 + [Im(x ) - Im(x*)] 2]} (6.13) „ n m m m m=0 where h^ = no. of highest harmonic considered. The values x 1 were calculated using the Runge-Kutta fourth order method m to i ntegrate the d i f f e r e n t i a l equation, from which the fa s t Fourier trans-form was used to obtain the c o e f f i c i e n t s x'. The period was determined m from the time between zero crossings. The program was w r i t t e n i n double p r e c i s i o n , with a step s i z e that was progressively decreased. The magnitude of the even harmonics i s an i n d i c a t i o n of the accuracy obtained, as they should decay to zero. However, the transients were s t i l l s i g n i -f i c a n t , because with approximately u n i t amplitude of the fundamental -4 _5 the even harmonics were s t i l l i n the range of 10 . to 10 when the c a l c u l a t i o n was stopped a f t e r s i x t y seconds of computing time. The values x' are not intended to represent the ultimate s o l u t i o n to the problem m but only to v e r i f y the r e s u l t s of the method to wit h i n reasonable accuracy. Table 6.2 l i s t s the r e s u l t s obtained, along with the values x'. m As may be seen, a l l c o e f f i c i e n t s are i n error by l e s s than 0.02% of the •39 magnitude of the fundamental. However, these results do not demonstrate the virtue of the method as well as does Fig. 6.4, which gives error behaviour in terms of computer time used taking different values of k. This figure indicates that the solution which can be obtained in 0.1 second i s sufficiently accurate for engineering purposes in many cases. The results in Table 6.2 were obtained i n 2.4 seconds using k = 100. In contrast, for the f i r s t and fastest cycle of Newton-Raphson (fastest because a large step size was used for the f i r s t cycle) over two seconds were required, and after one cycle a l l coefficients were s t i l l in error by more than one percent. Continuation of the curve i n Fig. 6.4 for values of k>400 was not deemed feasible because of the limitations of the reference data. This completes the oscillator example. Besides demonstrating that the method developed in chapter III can solve this problem, i t has illustrated that a good solution i s possible for Rayleigh's equation with the parameter equal to unity, and that the solution can be obtained more rapidly than by the usual computer-aided technique of time-domain integration. Table 6.2 Rayleigh O s c i l l a t o r Results (Reference values are i n brackets) HARMONIC 1 3 5 7 Real Part 1.068 (1.068) 1.500 x 10" 2 (1.503 x 10~ 2) -3.901 x 10~ 3 • (-3.900 x 10~ 3) -6.053 x 10~ 4 (-7.000 x 10~ 4) Imaginary Part 0.0 (0.0) -3.935 x 10" 2 (-3.923 x 10~ 2) -3.241 x 10~ 3 (-3.268 x 10~ 3) 2.497 x 10" 4 (4.369 x 10~ 4) Fundamental frequency = 0.9430 rad./sec. (0.9430) i 1 1 1 1 Time 0.1 0.3 1.0 3.0 10.0 (sec.) 42 I I I . A T r a n s i s t o r Amplifier Example This c i r c u i t i s chosen to i l l u s t r a t e the use of the ideas i n chapter IV for the analysis of exponential n o n l i n e a r i t i e s . The c i r c u i t i s given i n F i g . 6.5. Vj cos t 0.25 <>vout (£=99 ) Simple Transistor Amplifier F i g . 6.5 It i s not d i f f i c u l t to conceive of a s i g n a l source and b i a s i n g r e s i s t o r combination which together would have a Thevenin equivalent s i m i l a r to that d r i v i n g the t r a n s i s t o r here. S i m i l a r stages without feedback are sometimes used i n simple t r a n s i s t o r a m p l i f i e r s , so the problem i s relevant to the r e a l world. The s i m p l i f i e d large s i g n a l model of the b i p o l a r t r a n s i s t o r of F i g . 6.6 i s used. 43 b o c ic=0.99,e = 99ib i n - 6 30vhp ie=iO e D e T r a n s i s t o r Model F i g . 6.6 It i s evident that a l l high frequency e f f e c t s are neglected, a l l ohmic resistances are neglected (they may be taken up i n the external r e s i s t a n c e s ) , £ i s assumed constant, and the diode c h a r a c t e r i s t i c i s approximated as an exponential with no c o r r e c t i o n terms. The parameters are chosen to be i n the range encountered i n actual devices. Using t h i s model, an equivalent c i r c u i t i s drawn, with v and i as defined i n F i g . 6.7. Woo A )0.25 + - Vj cos t 0.99 i S 10'-Se30v Equivalent C i r c u i t f o r the T r a n s i s t o r Amplifier F i g . 6.7 44 By i n s p e c t i o n , i = H f 6 e 3 0 v (6.14) v = 0.25 + v cost - 10 i (6.15) v Q = 10 - 5000(0.99) i (6.16) choosing f(y) = y (6.17) Eq. (6.14) i s w r i t t e n as a d i f f e r e n c e equation: • j I • / A \ T O " 6 3 C H R I I A ~ 6 3 0 v j _ 1 l = i J + 2 (Ay) 10 e - 10 e + 1 0 ~ 6 e 3 0 v J ( A y ) (6.18) A program to analyze this c i r c u i t i s given i n appendix I I I . I t c a r r i e s four frequency components at each step, and uses the methods of chapters II and IV. The information about the stage which i s of i n t e r e s t i s the voltage a m p l i f i c a t i o n and d i s t o r t i o n at various s i g n a l l e v e l s . Voltage a m p l i f i c a t i o n i s calculated as the output voltage at the fundamental frequency divided by the input s i g n a l voltage, and harmonic d i s t o r t i o n as the t o t a l harmonic power r e l a t i v e to fundamental output power. Table 6.3 gives output data for various values of the input voltage v^, and F i g . 6.9 and F i g . 6.10 display t h i s data g r a p h i c a l l y . I t i s i n t e r e s t i n g to note, as shown i n F i g . 6.9, that the gain increases as the input voltage i s increased. F i g . 6.8 demonstrates, i n an exaggerated fashion, why t h i s i s so. I t i s e a s i l y seen that although negative excursions of the out-put are l i m i t e d , small increases i n the p o s i t i v e excursion at the input may cause large increases i n the fundamental component at the output. F i g . 6.8 Dependence of Gain on Amplitude Note that F i g . 6.9 also indicates the point at which the t r a n s i s t o r would saturate. I f t h i s e f f e c t were also modelled, a sharp decrease i n gain could be expected beyond s a t u r a t i o n . F i n a l l y , i t i s worth noting that the time required to analyze the stage i s only about 0.02 second. While t h i s i s a very simple stage, such a short time does hold promise that this method could be used to r e p e t i t i v e l y analyze the c i r c u i t s u f f i c i e n t l y r a p i d l y to make i t s use i n a c i r c u i t optimization program not only f e a s i b l e but p r a c t i c a l . Table 6.3 Behaviour of Transistor A m p l i f i e r Input voltage Output voltage Harmonic D i s t o r t i o n (%) Voltage A m p l i f i -cation Time f o r preceding l i n e (se D-C F i r s t Harmonic Second Harmonic Third Harmonic 0.001 - 13.71 -0.1386 -0.0005578 -0.0006760 0.00400 ! -138.6. 0.027 0.003 13.71 -0.4158 -0.004996 -0.0003567 0.0145 -138.6 0.019 0.01 13.66 -1.387 -0.05533 -0.0003820 0.159 -138.7 0.019 0.03 13.22 -4.184 -0.4878 -0.008315 1.36 -139.5 0.019 0.04 12.85 -5.607 -0.8528 -0.01847 2.31 -140.2 0.018 0.05 12.38 -7.050 -1.305 -0.03284 3.43 -141.0 0.018 0.06 11.84 -8.514 -1.831 -0.04793 4.63 -141.9 0.018 0.07 11.25 -9.949 -2.380 -0.05722 5.72 -142.2 0.019 0.06+ Transistor saturated-model no longer v a l i d 48 Nega tive Voltage Amplification A I—€>Transistor saturating -model no longer valid 0.003 V: F i g . 6.9 Voltage A m p l i f i c a t i o n as a Function of Input Voltage 49 % Harmonic Distortion A 10 --F i g . 6.10 Harmonic D i s t o r t i o n as a Function of Input Voltage 50 VIII. CONCLUSION A general computer-aided method of e f f i c i e n t l y determining the steady-state response of nonlinear networks has been i n v e s t i g a t e d . The following contributions have been made. A method for w r i t i n g the nonlinear network equations for forced networks using nodal v a r i a b l e s has been discussed. I t i s f e l t that such a d e s c r i p t i o n i s often more convenient than the hybrid equations of state commonly used for nonlinear network a n a l y s i s . The r e s u l t i n g nodal equations are i m p l i c i t so have no d i r e c t s o l u t i o n . Instead, a computer-aided numerical technique, i n v o l v i n g an extra parameter, has been developed. Extensions have been given which enable analysis of a class of o s c i l l a t o r s i n the steady-state. I t has been demonstrated that nodal v a r i a b l e s are unsuitable, but a s a t i s f a c t o r y a l t e r n a t i v e has been developed. A numerical technique for t h i s case has also been developed, s i m i l a r to the one used for forced networks. The problem of analyzing c i r c u i t s containing exponential non-l i n e a r i t i e s has been considered. Bessel functions have been introduced as an a i d to s o l v i n g this problem. Use of these contributions i n analyzing networks has been demonstrated with examples. Further development of these methods i s c e r t a i n l y p o s s i b l e . On the basis of our experience with the algorithm, some d i r e c t i o n s f o r further work are i n d i c a t e d : 1. The use of f i r s t d i f f e r e n c e approximations has been found, i n most numerical methods, to be slow and inaccurate compared to higher order approximations. I t should be p o s s i b l e to develop the methods presented i n t h i s thesis to incorporate higher order approximations. 2. The method could be used i n conjunction with an optimizing algorithm. In t h i s way the value of the l i n e a r elements might be automatically c a l c u l a t e d so as to minimize the nonlinear e f f e c t s i n the network, or the power input might be minimized subject to gain or d i s t o r t i o n requirements. In conclusion, an a l t e r n a t i v e to e x i s t i n g computer-aided methods f o r determining the steady-state response for nonlinear networks has been developed, and has been shown to have some advantages over e x i s t i n g methods. APPENDIX A LISTING OF PROGRAM FOR EXAMPLE 1 ' C I H I S I S I-XAI-'IP Lb N U M o b R O N E — A S E R I E S C O N N E C T I •) N U F A V O L T A G t SUURCb,"53' C R E S I S T O R AND N U N L I N E A R I N D U C T O R C E I S T H E V O L T A G E A C R O S S T H E I N D U C T O R «~On D ' M 1>W c < Et) C ON P O N E NTS A R E C T H E R E A L P A R T S OF T H E C O E F F I C I E N T S S T A R T I N G W I T H T H E DC C O M P O N E N T , A N D C T H E E V E N N U M B E R E D C O M P O N E N T S A R E T H E I M A G I N A R Y P A R T S C I I S ' I N 1 F R O M E Q U A T I O N 6 . 1 > C I 1 I S I — A T T H E P R E V I O U S I T E R A T I 0 N — R E A L E ( 6 6 ) , 1 ( 6 6 ) , 1 1 ( 6 6 ) , O N E G A ( 6 6 ) , D I ( 6 6 ) , M O , L A M D A ( 6 6 ) I f Y ( 6 6 ) " ~ C " N ( ' l ) I S T W I C E T H E N U M B E R OF H A R M O N I CS C A R R I E D "AND M U S T ' B E A N I N T E G E R C POWER OF TWO I N T E G E R N( 1 ) ,N ( 1) ,,.|( 1 } = b i + C K I S T H E • N U M B E R OF S T E P S T O BE. T A K E N 8 0 RE AD ( 6 , 1 0 1) N , K I F ( K. E Q . 0) ST 0 P ~ " " C DMU I S T H E S T E P S I Z E D M U = 1 . / F L O A T ( K) r t r l FO KM AT ( k I 5 ) N 1=N( 1) N 2 = N 1+ 2 X'N'1 = F'L0 A T t N 1) " '"' " " ~ — P I = 3 . 1 4 1 1 3 9 2 7 C T H E N E X T 9 S T A T E M E N T S S E T U P I N I T I A L C O N D I T I O N S A N D G E N E R A T E T H E C F R E Q U E N C Y — A R R A Y O M E G A DO 1 0 J = l , 6 6 Y( J ) = 0 . •• • • L A N D AC J ) = 6. " • I 1 ( J ) = 0. E( J ) = 0. rn ON EG A ( J ) - F L O AT ( J/2H -P I Y ( 3 ) = 3 . b ( 3 ) = 3 . C I N I T I A L I Z E T H E C L O C K S = S C L O C K ( 0, ) C T H E F O L L O W I N G DO L O O P — I S T H E A C T O A L / A L G O R I T H M DO 2 0 J = 1 , K C I N T E G R A T E E ( D I V I D E BY J * O M E G A ) TO GET L A M D A DO 3 0 J 1=3 ,N 2 , 2 LAMD A ( J 1) = E ( J 1+ 1) / O ME G A( J 1) 3 0 L AM 0 A( J 1+ 1) = ™E ( J 1 ) / O ME G A ( J 1 ) L A n D A( 1 ) - 0. L A N D A( 2) = 0 . C T R A N S F O R M TO T I M E D O M A I N ~ " C"AL L " F 0 U R'2"( L AMD A, N , 1, + 1 T R 1) ' ~ DO 4 0 J 1= 1 , N 1 X = L A M D A ( J l ) *rti 1 ( J 1) - ••-'( X +X-i-X-:-X ) / X N 1 C T R A N S F O R M B A C K TO F R E Q U E N C Y D O M A I N C A L L F 0 0 R 2 ( I , N , 1 , - 1 , 0) "C" T H E F O L L O W I N G DO L O O P A D V A N C E S T H E S O L U T I ON " B Y ' l J N E S T E P DO 5 0 J 1 = 3 ,N 2 U I ( J 1 ) = I ( J 1 ) » I 1 ( J 1 ) . • 1 1 ( J 1 ) - I ( J I ) 5 0 E ( J l ) =E ( J l ) + M U * D I ( J l ) + 1 ( J l ) *DMU 2 0 . i-!U = MU+ DMU C : R E C O R D T H E T I M E V U S E D " " ' S= S C L O C K ( ' S ) C C R E A T E AND P R I N T T H E O U T P U T I N T H E S A M E FOR,-I AS E . B A I L Y U S E D 1 0 2 F D KH A T ( 4 OH I H A K M U N I C C O S I N E C O E F . S I N E C U f F . ) DO 6 0 J = 3 , N 2 , 2 L= J / 2 L A M D A( J ) = 2 . * E ( J+ 1 ) / O M E G A ( J ) L A m A ( J+ 1) = E ( J ) / I ' J r-l E G A ( J ) * 2 . W K i I t l b , 1 0 0 ) L , L M H J A l J ) , LA i ' lUA I J + 1 ) 6 0 LAiMO A( J+ 1) = « £ ( J ) /O i - iEG A ( J ) * 2 . 1 0 0 F O R M A T t I X , I 5 , i i X , E 1 2 . 6 , 3 X , E 1 2 . 6 ) WRITE (6, 103) S,K " " " ' " ~ ~ 1 0 3 F OR M AT ( 1 OT I M E I'J SE 0 = ' , F 1 0 . 6 , ' S E C . S T E P S T A K E N = ' , I 5 ) C C R E A T E AND P R I N T E . B/-\ I LY ' S E R R O R C R I T E R I O N LMI-ID At 1) = U. ; : LAMI) A( 2) = 0 . C A L L F O U R 2 ( L A N D A , K i 1 , + 1 , =» 1) C A L L " F 0 U R 2 ( Y , i - l , 1,+ 1 , - 1 ) ' " " " ~ ~ C A L L F O U R 2 ( £ , M , 1 , + 1 , - 1 ) E P S I LN = 0 . DU / (j j = i , o 4 X= L A MO A ( J ) / 2 . Z = E ( J ) + X + X * X * X " Y ( J ) "7 0 EPSI LfM = EPSI L I M + Z * Z / 6 4 . " W R I T E ( 6 , 1 0 4 J E P S I L N 1 0 4 F O R M A T ( 3 9 H 0 E . b A I L Y ' S E R R O R C R I T E R I O N E P S I L O N [ S , E 2 0 . 1 0 ) (30 " T u >3 u  E N D APPENDIX B LISTING OF PROGRAM FOR EXAMPLE C T H I S I S E X A M P L E 2=~AN A N A L Y S I S OF R A Y L E I G H ' S E Q U A T I O N D b C F I S D E L T A X P R I M E D IN THE N O T A T I O N OF C H A P T E R I I I C B l I S B FROM THE P R E V I O U S I T E R A T I O N R E A L M U , X ( 1.8) , B ( 1 8 ) , B 1 ( 1 8 ) , DB ( 18 ) , F ( 1 8 ) , E( 18 ) COMMON /W/WO I N T E G E R M( 1 )  COMMON / X / X , B , E C THE N E X T 16 L I N E S SET UP C O N S T A N T S AND I N I T I A L C O N D I T I O N S .. „_ N( 1) = 16 - . . -C K I S THE NUMBER OF S T E P S TO BE T A K E N K = 4 0 0 D M U = l . / F L O A T ( K) ; MU = 0 . C I N I T I A L L Y , THE A N G U L A R F R E Q U E N C Y OF THE F U N D A M E N T A L I S ONE WO= l . _ DO 10 J = l , 1 8 X( J ) = 0 . F ( J ) = 0 .  10 B 1 ( J ) = 0 . C THE F O L L O W I N G I N I T I A L C O N D I T I O N I S KNOWN FROM P E R T U R B A T I O N THEORY C BUT COULD HAVE B E E N D E T E R M I N E D IN T H E PROGRAM AS. WELL . .. _ X ( 3 ) = l . C I N I T I A L I Z E THE C L O C K S=SCLOCK( 0 . )  C THE A C T U A L A L G O R I T H M B E G I N S HERE DO 20 J = 1 , K C . M I S A S U B R O U T I N E WHICH T A K E S THE " S O L U T I O N , X , AND F R O M . I T C A L C U L A T E S C THE F O U R I E R C O E F F I C I E N T S ( B ) FROM THE N O N L I N E A R TERMS C A L L M DO 3 0 L = l , 1 8  DB( L) =B( D - B K L ) B i t L) =B( L) 3 0 F( L) =F( L ) +MU*DB( L) +B( L ) * D M U . .... _ C C A L C U L A T E THE NEW F R E Q U E N C Y OF THE F U N D A M E N T A L W O = S Q R T ( 1 . - F ( 3 ) / X { 3 ) ) C C A L C U L A T E NEW V A L U E S OF A L L H A R M O N I C S BUT THE F U N D A M E N T A L  DO 4 0 L = 5 , 1 8 F R = F L O A T ( [ L - D / 2 ) *W0 4 0 . X( L ) = F ( L ) / ( l . - F R * F R ) -C FROM HERE TO THE NEXT B L A N K COMMENT C A R D , A NEW V A L U E OF THE C F U N D A M E N T A L I S C A L C U L A T E D C A L L M . R = B ( 4 ) X ( 3 ) = X ( 3 J + 0 . 0 0 1 C A L L M ^ X( 3 ) = X( 3 ) - 0 . 0 0 1 - 0 . 0 0 1 * R / ( B ( 4 ) - R ) C 2 0 MU=MU+DMU  C RECORD THE T I M E T A K E N S = S C L O C K ( S ) C THE R E M A I N I N G S T E P S IN THE M A I N PROGRAM A R E FOR P R I N T I N G THE _ OUTPUT C AND C A L C U L A T I N G AND P R I N T I N G AN ERROR E S T I M A T E ES E S = 0 . . E ( 3 ) = 0 .  E{ 4 ) = X ( 3 ) * W 0 W R I T E ( 6 , 1 0 0 ) ...100 FORM AT ( ' 1 C A L C U L A T E D X ... TRUE +. X« / ' O F R E Q U E N C Y R E A L P A R T I M A G I N A R Y P A R T R E A L PA + RT I M A G I N A R Y P A R T ' ) DO 50 J = 3 t l 7 , 2 ;  F R = W O * F L O A T ( J / 2 ) ' " C READ IN THE ' T R U E ' V A L U E S R E A D ( 5 , 1 0 1 ) X R , X I 1 0 1 FO RM AT ( E 2 0 . 1 0 ) W R I T E ( 6 , 1 0 2 ) F R , X ( J ) , X ( J + l ) , X R , X I 1 0 2 FORMAT! 1 X , E 1 2 . 6 , 6 ( 6 X , E 1 2 . 6 ) )  ER=X( J ) - X R E I = X ( J + l ) = X I _„50 .. E S = E S + E R * E R + E I * E I . W R I T E ( 6 , 1 0 3 ) K , S , E S 103 FORMATt • OSTEPS - T A K E N = ' , I 5 , « T I M E T A K E N = « , F 1 0 . 6 , « S E C . , + E R R O R = 1 , E 2 0 . 1 0 )  STOP END . . S U B R O U T I N E M . . . . . . . .. C M I S A S U B R O U T I N E WHICH T A K E S THE S O L U T I O N , X , AND FROM I T C A L C U L A T E S C THE F O U R I E R C O E F F I C I E N T S ( B ) FROM THE N O N L I N E A R TERMS R E A L X( 1 8 ) , B ( 1 8 ) , E ( 1 8 )  I N T E G E R N ( 1 ) COMMON / X / X , B , E ..COMMON /W/WO .. N( 1) =16 C THE F O L L O W I N G DO LOOP D I F F E R E N T I A T E S X TO GET E DO 10 K = l , 1 8 , 2  F = W O * F L O A T ( K / 2 ) E( K ' ) = - X ( K + 1 ) * F 10 E ( K + 1 ) = X ( K ) * F : C GO TO THE T I M E DOMA IN C A L L F 0 U R 2 ( E , N , 1 , + 1 , » 1 ) DO 2 0 K = l , 1 6  C E V A L U A T E THE NONL INE A R I T Y — P U T THE R E S U L T IN B 2 0 B( K) = (E { K ) - E ( K ) * E ( K ) * E ( K ) / 3 . ) / 1 6 . ... . C . GO B A C K TO F R E Q U E N C Y DOMA IN . ' C A L L F 0 U R 2 ( B , N , 1 , - 1 , 0 ) R E T U R N END APPENDIX C LISTING OF PROGRAM FOR EXAMPLE 3 59 C THIS IS EXAMPLE NUMBER 3 ~» i\ TRANSISTOR AMPLIFIER C V AND I ARE AS DEFINED IN F I G . 6.6 C VOiIO ARE THE DC COMPONENTS C V I , II ARE THE FUNDAMENTAL COMPONENTS C E T C . IMPLICIT R E A L(I ,M)  C X IS AN ARRAY CONTAINING THE OUTPUT REAL X(4) WRITE(6,100) 100 FORM A T ( 5 X , 1 V I DC FUNDAMENTAL + SECOND HARMONIC THIRD HARMONIC HARMONIC DISTORTION') C THE FOLLOWING 13 STATEMENTS ARE INITIAL CONDITIONS  20 READ(5,101)VI 101 FORMAT!F20.10) I F( V I . LT . 0 . ) STOP . . . . . . MU=0. DMU=0.1 vo=o .25 . ; V 1 = V I V2 = 0. V 3 = 0 . _ 10=0. I 1=0. 12 = 0.  13=0. C INITIALIZE THE CLOCK S=SCLOCK( 0. ) . . ._; _ C THE ALGORITHM BEGINS NOW DO 10 J = l , 1 0 XV = V1*30. . ' " C THE A'S ARE EXP(30 .*V1) A0=BESSI0(XV) . _ . A 1= BESS I K XV) . „ _ ; . . „ A 2 = 2 . * ( A 0 - 2 . * A 1 / X V ) A3=2.*( ( 8 . / ( X V * X V ) + l . ) * A 1 » 4 . * A 0 / X V ) A 1 = 2 . * A 1  C THE B ' S ARE EXP(30 .*V2) B 0 = B E S S I 0 ( V 2 * 3 0 . ) B 2 = 2 i * B E S S I l ( 3 0 . * V 2 ) , C THE C'S ARE EXP(30.*V3) CO=BESSIO( V3*30. ) C3=2.*BESSI1 (30.-.--V3)  I D C =l . E » 6 * E X P ( 30. *V0) C THE NEXT 4 STATEMENTS ARE THE CONVOLOTION OF IDC, THE A ' S , B ' S , A N D C'S I N 0 = 1 D C * 1 A0 * B 0 * C 0 +A2 * B 2 * C 0 /2 .+A3 * B 0 * C 3 /2 . +A l* B 2 . * C 3 / 4 . ) IN 1= IDC*( A1^B0-:=C0+ A l * B 2 * C 0 / 2 . + A 3 * B 2 * C 0 / 2 . + A 2 * B 0 * C 3 / 2 . + A0* B 2 * C 3/2 .+ +A2 * B 2 * C 3 / 4 . ) IN2= IDC*(-A2*B0*C0+ A 0 * B 2 * C 0 + A l * B 0 * C 3 / 2 . + A 1»B2 : : ; C3/ 4. + A3*B2*C3/2 . ) I N 3 = ID C * ( A 3 * B 0 * C 0 + A1 * B 2 * C 0 / 2 . + A 0 * B 0 * C 3+A 2 * B 2 * C 3 / 2 . ) C CALCULATE DI 0 ' I 0 = I N 0 - I P 0 D I 1 = I N 1 - I P 1 D I 2 = I N 2 « I P 2 D I 3 = I N 3 ° I P 3  C STORE THE VALUES OF IN IN IP FOR THE NEXT STEP IP0=IN0 I P 1 = IN 1 ; „ „ . IP2=IN2 IP3=IN3 C CALCULATE NEW I 'S 60 1 0 = 1 0 + M U * D I 0 + I N 0 * D M U I 1= I 1+MU*D I 1+ I N 1 * D M U 1 2 = 1 2 + M U * D I 2 + I N 2 * D M U 13= I 3 + M U * D I 3 + I N 3 * D M U C C A L C U L A T E NEW V 1 S V 0 = 0 . 2 5 - 1 0 . * I 0 V1=V I - 1 0 . * I 1 V 2 = - 1 0 . * I 2 - V 3 = - 1 0 . * I 3 10 MU=MU+DMU C RECORD THE T I M E USED S = S C L O C K ( S)  X ( 1 ) = 2 0 . - 5 0 0 0 . * I 0 * 0 . 9 9 X ( 2 ) = - 5 0 0 0 . * I 1 * 0 . 9 9 X( 3) = - 5 0 0 0 . * I 2 * 0 . 9 9 - .. .. X ( 4 ) = - 5 0 0 0 . * I 3 * 0 . 9 9 H = ( X ( 3 ) * X ( 3 ) + X ( 4 ) * X ( 4 ) ) / ( X ( 2 ) * X ( 2 ) ) 1 0 2 FORM AT( 1 X , 6 ( E 1 5 . 6 , 5 X ) ) G = X ( 2 ) / V I W R I T E ( 6 , 1 0 3 ) S , G 103 F O R M A T l 1 T I M E U S E D = ' , F 1 0 . 6 , ' S EC GA IN = ' , F I 2 . 6 / ' ' ) GO TO 20 END 61 REFERENCES 1. Chua, L. 0., "Computer-Aided Analysis of Nonlinear Networks", i n Computer-Aided Design, F. F. Kuo and W. G. Magnuson, J r . , eds., Englewood C l i f f s , N. J . : P r e n t i c e - H a l l , 1969, p. 123. 2. Branin, F., "Computer Methods of Network An a l y s i s " , i b i d . , p. 92. 3. "Computer-Aided Design Goes Beyond F i r s t Generation Programs". Insert i n E l e c t r o n i c s , A p r i l 13, 1970. 4. K r y l o f f , N. and Bogoliuboff, N. Introduction to Nonlinear Mechanics. Princeton, N. J . : Princeton U n i v e r s i t y Press, 1947. 5. Cole, J . D., Perturbation Methods i n Applied Mathematics. Waltham, Mass.: B l a i s d e l l Publishing Co., 1968. 6. B a i l y , E. "Computer Analysis of Nonlinear Forced O s c i l l a t i n g Networks." Presented at 1970 I.E.E.E. Sixth Region Conference, S e a t t l e , Washington, May 26-28, 1970. 7. B a i l y , E. "Steady-State Harmonic Analysis f o r Nonlinear Networks". Ph.D. t h e s i s , Stanford U n i v e r s i t y , Stanford, C a l i f o r n i a , August, 1968, p. 1. 8. N e i l l , T. B. M., "Improved Method of Analyzing Nonlinear E l e c t r i c a l Networks". E l e c t r o n . L e t t . , 1969, 5, p. 13-15. 9. Heywood, D. R., and Moore, A. D., "Comment on an Improved Method of Analyzing Nonlinear E l e c t r i c a l Networks". E l e c t r o n . L e t t . , 1969, 5, p. 269-271. 10. Cunningham, W. J . , Introduction to Nonlinear A n a l y s i s . New York: McGraw-Hill, 1958, p. 63. 11. Bellman, R., Perturbation Techniques i n Mathematics, Physics, and Engineering. New York: Holt, Rinehart and Winston, 1964, p. 54. 62 12. Ebers, J . J. and M o l l , J . L., "Large Signal Behaviour of Junction T r a n s i s t o r s " . Proc. IRE 42, 1954, pp. 1761-1772. 13. Gold, B. and .Rader, C., D i g i t a l Processing of Signals. New York: McGraw-Hill, 1969, p. 204. 14. B a i l y , E., "Steady-State Harmonic Analysis for Nonlinear Networks". Ph. D. t h e s i s , Stanford U n i v e r s i t y , Stanford, C a l i f o r n i a , August, 1968, p. 102. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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-0102065/manifest

Comment

Related Items