LIBRARY OF THE 
 
 UNIVERSITY OF ILLINOIS 
 
 AT URBANA-CHAMPAIGN 
 
 510.84 
 
 nOel60-170 
 cop .3t 
 
The person charging this material is re- 
 sponsible for its return to the library from 
 which it was withdrawn on or before the 
 Latest Date stamped below. 
 
 Theft, mutilation, and underlining of books 
 are reasons for disciplinary action and may 
 result In dismissal from the University. 
 UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN 
 
 NOV 15 
 
 iff 1 
 
 m 
 
 1 1«» 
 
 L161 — O-I096 
 
Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/incrementalcharg166nusp 
 
r 
 
 gg DIGITAL COMPUTER LABORATORY 
 
 UNIVERSITY OF ILLINOIS 
 URBANA, ILLINOIS 
 
 .3 
 
 REPORT NO. 166 
 
 AN INCREMENTAL CHARGE METHOD FOR THE ANALYSIS OF NONLINEAR CIRCUITS 
 
 Stephen John Nuspl 
 
 August 6, I96I+ 
 
 'his work is being submitted in partial fulfillment of the requirements for 
 ■he Degree of Master of Science, August, 196k, and was supported in part by 
 he Office of Naval Research under contract Nonr-l834(l5) . ) 
 
DIGITAL COMPUTER LABORATORY 
 UNIVERSITY OF ILLINOIS 
 URBANA, ILLINOIS 
 
 REPORT NO. 166 
 
 AN INCREMENTAL CHARGE METHOD FOR THE ANALYSIS OF NONLINEAR CIRCUITS 
 
 by 
 
 Stephen John Nuspl 
 
 August 6, 1964 
 
 This work is being submitted in partial fulfillment of the requirements for 
 the Degree of Master of Science, August, 196k, and was supported in part by 
 the Office of Naval Research under contract Nonr-l834(l5) . ) 
 
ACMOWLEDGMEOTS 
 
 The author would like to express his sincerest thanks to his advisor,, 
 Professor W„ Jo Poppelbaum^ for his good counsel, support and encouragement „ 
 Thanks are also extended to the past and present Task 15 group at the Digital 
 Computer Laboratory with whom the ideas in this thesis were often discussed and 
 to Mrs. Phyllis Olson who prepared the final manuscript.. The excellent drawings 
 in this report are a result of the good work of the Digital Computer Laboratory 
 Drafting Department under Kenneth. Law, 
 
TABLE OF CONTENTS 
 
 Page 
 
 lo INTRODUCTION o ..„..,„.,....,.....„,,,, , 1 
 
 2, DEVELOPMENT OF THE INCREMENTAL CHARGE METHOD ...„,,..,, 3 
 
 2,1 Circuit Ordering and Topological Restrictions , o o « « „ » 3 
 
 2o2 The Charge Distribution Criterion , . . « = » , , . „ . » » 5 
 
 2„3 Node Voltage Calculation . o o » . . , . . » . . o , , „ „ 12 
 
 3. CALCULATION OF IMPEDANCE FACTORS . o . » . , » . » . , „ , « , „ I5 
 k. VOLTAGE -CONTROLLED SUBNODES .=.,=, o »,,,..„„,, „ 20 
 
 5. DISCUSSION OF CORRECTION FORMULA »,. o„o .„„.,„„., „ 26 
 
 5.1 Error Term for a Branch .0, o , o.» o .„„„„„„„ „ 26 
 
 502 Compensation of the Error « = o. = .<.„„.„.,„,„ , 28 
 
 503 A Predictor-Corrector Type Method o o = . . » « o o . . , . 30 
 
 ERRORS AND STABILITY 
 
 BIBLIOGRAPHY 
 
 33 
 
 60 1 Circuit Error and Response Error , » . » » , „ „ „ » „ „ o 33 
 
 602 Stability ..,„.,„„.,,„,,..„, ^S 
 
 RELAXATION OF TOPOLOGICAL RESTRICTIONS ...„.,.„,.„,„ 37 
 
 NONLINEAR AND TIME-VARYING PARAMETERS .,.,...,„.,,, kO 
 
 8.1 Sources o ,,» ..,..,„„.,..„.„.,,,„„ ^ I|.0 
 
 8.2 Nonlinear Elements ... o 00 ..,,,„ o, =,0 00 . 40 
 
 8.3 Semiconductor Diodes .,...<,, o <,,. o o , o., ,, 4l 
 80^ Transistors , ... o ... o,,,.„ ..^ ,,„„,„ „ ii5 
 8.5 Transmission Lines ...........„„„„,,„„ , 53 
 
 REALIZATION OF THE IQ METHOD . ...o ...,,,,,.„,„ , 58 
 
 9.1 The Program ..... .o,. ..... ... o., .o« =, 58 
 
 9o2 Sample Results .. .......«.„,.....,., „ 60 
 
 9.3 Efficiency .....,.......„,..,„„ „W' ° 6I 
 
 CONCLUSIONS .,..........„..„„„„..,„,„„ 63 
 
 65 
 
 APPENDIX .. o o .... „.,„„.„„,„,„.„.,„„„, , 66 
 
 ■IV- 
 
AW INCREMENTAL CHARGE METHOD FOR THE ANALYSIS OF NONLINEAR CIRCUITS 
 
 Stephen John Nuspl^ M.S. 
 
 Department of Electrical Engineering 
 
 University of Illinois, 196^^ 
 
 Conventional numerical methods for nonlinear transient circuit analysis 
 apply standard numerical procedures to the equations obtained from Kirchhoff's 
 lawSc In this thesis is developed a method of analysis in which only one of 
 Kirchhoff's laws is forced to hold and the other is used to make num.erlcal cor- 
 rections » The method developed here assumes that the current law holds and the 
 voltage law may he violated „ 
 
 To obtain sufficiently simple numerical formulas;, a circuit with a 
 tree-type topology is assumed. Relying on this the basic charge distribution 
 and node voltage calculation formulas are derived for the general branch which 
 may consist of a resistance, a capacitance, an inductance, a voltage source or 
 any series combination of these „ Since this method favors current controlled 
 nonlinear elements, an analogous but more restricted method for conveniently 
 including voltage -controlled devices is developed, A brief discussion of errors 
 and their compensation and of stability is included. 
 
 Because of the importance of diodes and transistors, the technique by 
 which they can be efficiently analyzed by the method is developed. The advan- 
 tages are that the method provides for the models the means for making numerical 
 corrections e A technique by which transmission lines can be handled is also 
 described, 
 
 A program based on this method is presented through a brief description. 
 To illustrate the use of the method four sample circuits and the waveforms they 
 produced are given. It is concluded that the method has some decided advantages 
 which should prove to make it useful, but not enough experience has been had 
 [With it to make any meaningful comparison with other efficient methods for non- 
 linear circuit analysis. 
 
1, INTRODUCTION 
 
 In the analysis of transients in nonlinear circuits one usually 
 attempts to obtain a first-order solution by the use of a linear transform 
 method which relies on a piecewise linear approximation of the nonlinear 
 elements. When this is no longer accurate enough or leads to cumbersome 
 algebraic manipulations, a more basic approach must be used. This normally 
 consists in the writing of the differential equations of the circuit obtained 
 from Kirchhoff 's laws and manipulating these until a set of first-order linear 
 equations in the form 
 
 X, 
 
 kjm, 1 L ^-'- 
 
 = a. 
 
 m, n 
 
 X. 
 
 1 
 
 n,l 
 
 is obtained. a is either a constant or a variable depending on any x. „ For 
 
 any but the simplest of circuits a numerical integration must be used. This 
 
 method has a number of disadvantages which become more and more severe as the 
 
 size of the circuit increases: finding the coefficients a, . is usually a time- 
 
 ki 
 
 consuming job; the evaluation of the coefficients at each step requires con- 
 siderable computer time since the coefficients in general contain some terms 
 related to the nonlinear elements; the addition of only one more element in the 
 circuit might require that new expressions for all coefficients be found and 
 programmed. These inconveniences severely restrict the practical use of this 
 technique for the solution of transients in larger nonlinear networks „ 
 
 Mechanizations of the above method which essentially separate the 
 ■ inear elements from the nonlinear are found in the literature. ^' ^^ The linear 
 'art of the circuit is usually described by a coefficient matrix which is auto- 
 latically assembled and inverted by program. The linear and nonlinear portions 
 •re then integrated separately and matched at convenient time intervals. This 
 lass of methods is excellent in that only a description of the circuit topology 
 
 ■1- 
 
-2- 
 
 and elements must be supplied. However, considerable programming effort is 
 required to make the methods functional and sometimes there are stability 
 
 problems . 
 
 The Incremental Charge (IQ) method described in this thesis in essence 
 goes back one step to more basic principles. In all of the above methods it 
 was assumed that Kirchhoff 's current and voltage laws hold and from these the 
 differential equations were obtained. In the Incremental Charge method we 
 impose one of Kirchhoff 's laws, allow the other to be temporarily relaxed and 
 use the deviations from this law as a means of making numerical corrections. 
 The process must always tend towards the reduction of the deviations. Programs 
 for the dc analysis of circuits ty the use of these principles have been success- 
 fully used'-^'^-' but the extension of these principles into the time domain does 
 not yet seem to have been attempted. 
 
2. DEVELOPMENT OF THE INCREMENTAL CHARGE METHOD 
 
 The basic philosophy behind the IQ, method is that one of Kirchhoff 's 
 laws will be allowed to be violated in order to permit numerical incremental 
 corrections o We can allow either the voltage around a circuit loop to be in 
 error by a small amount or we can relax the requirement that the algebraic sum 
 of the currents flowing into a node be zero„ It was found that allowing 
 Kirchhoff 's voltage law (KVL) to be violated and strictly enforcing Kirchhoff 's 
 current law (KCL) seemed to result in a more convenient formalization of a 
 numerical procedure. It is therefore this principle on which the IQ method is 
 based, 
 
 2=1 Circuit Ordering and Topological Restrictions 
 
 The first task is to devise a convenient scheme for presenting the 
 circuit topology to the computer and an order among the circuit elements which 
 the program follows during the execution of the method. To accomplish this each 
 aode of the circuit will be assigned a number which is also defined as the 
 "level" of that node. This notion of levels of nodes is introduced here only 
 for the reason that it conveys better the concept" of ordering among the nodes 
 than the word "number," The reference node is always assigned the number 1„ 
 Phese definitions are illustrated in Fig, 2,1, Between any two nodes there can 
 38 branches which may consist of a capacitance, a resistance, an inductance, a 
 Aoltage source or any combination of these. These branches also have numbers 
 issigned to them but since only the knowledge of which ones are under a partic- 
 |ilar node is required, the ordering can be arbitrary. 
 
 To define a sufficiently simple numerical procedure a number of 
 j-estrictions will be placed on the circuit topology. First it will be assumed 
 ■hat the topology is in the form of a tree as shown in Fig, 2,1, The forcing 
 
-k- 
 
 $ 
 
 B 
 
 FORCING FUNCTION 
 VOLTAGE SOURCE 
 
 BRANCH 2 
 
 >r-^ 
 
 NODE 5 
 
 B 
 10 
 
 B 
 
 I I 
 
 B 
 12 
 
 1— NODE L 
 
 EVEL 
 
 V 
 
 REFERENCE NODE 
 
 Figure 2a, Assumed Form of Topology for Development of Method. 
 
-5- 
 
 function to which the circuit must respond is in the branch at the apex of the 
 diagram and in the form of a voltage source. The second restriction is that 
 there be at most one branch between two nodes of which the reference node is 
 not one» These restrictions may seem drastic at present, but many nonlinear 
 circuit problems fall directly into this type of topology and many others which 
 do not conform to these assumptions can be handled but the error will be slightly 
 larger and other conveniences of the method are forfeited „ The second restric- 
 tion will eventually be eliminated entirely and the relaxation of the first will 
 be discussed later, 
 
 2.2 The Charge Distribution Criterion 
 
 There are essentially two steps in the numerical procedure. In the 
 first a current is injected into the highest node and passed to successively 
 lower level nodes via the branches under each„ When the current distribution 
 is completed, the currents in all branches are known, enabling us to calculate 
 the voltage across each element. The second step then consists in starting 
 from the reference node and determining what each successively higher level node 
 voltage will be. 
 
 To develop this method, let us consider the general node of Fig„ 2.2 
 which is assumed to have n branches immediately below it. In a time interval 
 At^ a current I or, equivalently, a charge 
 
 is injected into node k from higher level nodes. By a distribution process 
 which is to be developed branch j receives an amount of charge AQ.., The 
 assumption that KCL holds implies that 
 
-6- 
 
 AQ, 
 
 'i^T^" 
 
 ® 
 Iaq, 
 
 I"". 
 
 v^ *- O 
 
 ■ — »-0 
 
 11 I 
 
 ® 
 
 AQ 
 
 \ ^h 
 
 ■H. <.R. 
 
 X: 
 
 -6 
 
 NODE k 
 
 ,, TO LOWER LEVEL NODES 
 
 Figure 2„2o The General Node „ 
 
■7- 
 
 AQ^ = ZaQ (2.1) 
 
 The branches in Fig. 2.2 are not shown connected to node k to emphasize that 
 voltage V and all v.. 's may have different values if ICVL does not hold If we 
 
 K. J 1 
 
 required that it does hold, we have 
 
 + + + , , 
 
 V, . = V-, . = V^ . = „ „ „ = V (? p) 
 
 ki li 21 ni \^'^J 
 
 since the points (Y) to (n^ at the branch tops are connected to node k in the 
 actual circuits, If we perm.it this last set of equations to be violated, we are 
 allowing voltage differences 
 
 Av . . = V, . - V . . ( 2 o ^ ) 
 
 Ji ki ji v-o; 
 
 If this Av approaches zero, KVL will again hold. Hence an exact solution of 
 the circuit will be had if this condition holds over the entire network. There- 
 fore, the object of the numerical method will be to keep this voltage difference 
 as small as possible, or at least below a tolerable level. 
 
 The problem now is to distribute the charge AQl. among the branches in 
 such a manner that all branch voltages v. will experience an approximately equal 
 
 J 
 
 increase. To accomplish this we define a set of differences by 
 
 AQ.. = Q.. - Q., ^ . 
 
 A^Q . = AQ.. - AQ.,. . (2,^) 
 
 A^Q = A^Q.. - A^Q... ^. 
 Ji Ji j(i-l) 
 
 The voltages developed across the three types of elements are then given by 
 
-8- 
 
 Q.. 
 
 cji 
 
 
 V = li 
 
 2 
 
 ('A+ '> 
 
 >.2 
 
 (2.5) 
 
 Since we are allowing each branch to contain all three types of elements^ we 
 see ^hat only A^Q. . can he changed arbitrarily; the other charge differences 
 
 J 
 
 must be determined by equations 2.k. Also define 
 
 A + + + 
 
 Ji Ji j(i-l) 
 
 (2.6) 
 
 We are now in a position to define three series impedance factors between point 
 (T) of Figo 2,2 and the reference node: 
 
 Av ^..At. 
 cEjx 1 
 
 cj 
 
 AQ. 
 
 (2„7) 
 
 Ji 
 
 Av^„..At. 
 
 Rj A^Q. 
 
 (2.8) 
 
 Ji 
 
 Av^^. At. 
 LEji 1 
 
 ^'j a\. 
 
 (2.9) 
 
 Ji 
 
 where v , v and v^.^ are the voltages developed across these series impedance 
 cE RE -LiE 
 
 factors and which satisfy the relations 
 
 ^ji ^ ^cEji ^ ^REji + ^LEji 
 
 (2.10) 
 
 and 
 
 Av = Av . i- Av„^ . . + A\'-T ^, . . 
 "ji cEji RF.J1 LEji 
 
 (2.11) 
 
 We can also define a total branch impedance factor 
 
-9- 
 
 z . = z . + z^ . + z^ . 
 
 J cj Rj Lj 
 
 (2.12) 
 
 and two admittance factors 
 
 Y =-^ 
 
 (2a3) 
 
 n 
 
 Y, = Z Y. 
 
 (2aif) 
 
 Phese factors serve as a means of predicting branch voltage increases for a 
 5iven change in AQ. or in one of the higher differences » 
 
 Assume that the numerical procedure has been under way for a number of 
 :ycles and that after cycle (i-l) is completed there is a residual voltage error 
 ^^j(i-l) ^^"^^^^^ "°ds ^ ^^d branch j„ The first task is to raise the branch 
 roltage v to the node voltage v^^^,^) ^y injecting an appropriate amount of 
 charge which we shall define as A^Q.. into the brancho Assuming that the three 
 -mpedance factors are known for the branch, this condition requires A^Q.. be 
 such that 
 
 AQ 
 
 ^(i-D^^^Ki-D^^Si) , , ^^'S(i-i 
 
 At. ^cj ^ AT 
 
 ^Rj At, \.j 
 
 Av . / . . ^ 
 
 j(l-l) 
 
 his yields 
 
 (2.15) 
 
 ^ At. 
 
 J 
 
 Av - z i^^Jiiiill rz + z ) "^hilill 
 
 L j(i-l) cj \ At. )- ^^cj ^ ^Rj^ ^aE~~. 
 
 (2=16) 
 
 his is the basic "correction" formula in the method and will be discussed in 
 ore detail later. After determining 
 
-10- 
 
 for every branch and summing, we have distrihuted an amount of cb_arge 
 
 Wi - .| Vji 
 
 (2.18) 
 
 This leaves an additional charge to be distributed: 
 
 ^\i = ^\i - w^ ^'-^^ 
 
 Since it is assumed that at this point all branch voltages are at the level 
 
 V / X, this additional charge must be distributed so that all branch voltages 
 k(i-l)' 
 
 will increase the same amount. If branch j receives an additional A^Q^^ then 
 Since by definition 
 
 4Si = ''^(i-i) * 4si 
 
 ^«ji = ^'^Jd-D " 4Si 
 
 (2.21) 
 
 we get the result 
 
 4Si = 4Si = Vji '^■^^' 
 
-11- 
 
 The usefulness of this lies in the fact that the additional charge A Q will 
 
 A ji 
 
 cause an increase in branch voltage of 
 
 Vji ^cj At. ^ ^Rj At. ^ hj At. - At. ^^cj ^ ^Rj ^ ^Lj ^ = Y;^ 
 
 J 1 
 (2.23) 
 
 3r 
 
 A Q.. = Y.At.A.v^. (2.2^) 
 
 /hich gives the amount of charge which must be injected into branch j to raise 
 
 Lts voltage ^^v . We want the voltage of all n branches and that of the node 
 :o increase the same amount : 
 
 "^A^ji ^ ^A\i ^ constant for j = 1 to n (2„25) 
 
 fence 
 
 
 herefore 
 
 AQ.. Y. 
 A\i ""k 
 
 ^Si=^A\iXY:^ (2.28) 
 
-12- 
 
 This is the desired forimla for distributing the charge ^p9^^ among the branches, 
 The term 
 
 I 
 
 \ At.Av 
 
 1 A Kx 
 
 is in the form of an admittance and since the expression contains only terms 
 related to node k^ we could define it as the "admittance factor" for node k„ A 
 node impedance factor can therefore be defined by 
 
 2o3 Node Voltage Calculation 
 
 The above criterion is used to channel the charge injected into the 
 top node through all the branches and other nodes and into the reference node o 
 When this is completed, all the branch currents are known and therefore the 
 voltages across them can be calculated. The remaining step is to calculate the 
 
 node voltages v, . , 
 ^ ki 
 
 Through possible errors in the impedance factors all the branch 
 
 voltae-es v"^ will not have the same value „ Therefore some sort of mean value 
 ji 
 
 for V must be defined. There are a number of possible ways of determining 
 ki 
 
 this mean, the most primitive of which is the average value : 
 
 V . = i Z vt. (2o30) 
 
 ki n ji 
 
 A second formula may be obtained by imposing the conservation of energy conditioi 
 at node k, which implies that no energy is lost through the niMierical process o 
 This condition states 
 
-13- 
 
 ^i'^S.l = .^, ^Il-^Sl (2.31) 
 
 J=l 
 
 and hence we have 
 
 n AQ. . 
 
 ^-^5rji^ (-32) 
 
 This is effective when AQ and all AQ..'s have the same sign„ But, if 
 
 n 
 
 |ao . I « E Iaq.. I 
 
 which occurs if one of the branches under node k contains a dc source, numerical 
 instabilities may result. A modified version of the formula is 
 
 n 
 
 Z vT. Iaq. . I 
 
 j=i J" J" , , 
 
 \^-~ (2.33) 
 
 Z |aq.. I 
 
 When all ^Q . • ' s have the same sign energy is again conserved; otherwise the 
 node voltage will be the average between the branches delivering power and those 
 receiving it. 
 
 In a third method the impedance factors defined for the current 
 distribution are used: 
 
 n 
 
 Z Y.vT. 
 1=1 J J^ 
 
 k 
 
 ["his has the effect of latching the node voltage to the branch with the lowest 
 i-mpedance factor and therefore tends to stabilize the node voltage against 
 
-Ik- 
 
 numerical oscillations. Because of this stabilizing effect this last criterion 
 
 for the calc\ilation of the node voltages seems to he the best. However^ no 
 
 detailed analysis comparing the last two methods has yet been made. 
 
3. CALCULATION OF IMPEDANCE FACTORS 
 
 In the above discussion the knowledge of the values of Z , Z and Z 
 
 c"^ R L 
 
 for each branch was assumed. Since there are three parameters which must be 
 defined, we need three sets of equations for each branch. This makes impossible 
 the redefining of a set of impedances during each step of the analysis. Though 
 this would be desirable, it is not absolutely necessary. The initially defined 
 impedances will never change if all the elements are linear and will be 
 relatively constant over small ranges of voltages even if the elements are non- 
 linear. In the latter case the impedance factors must be redefined when the 
 actual branch voltage changes start to differ appreciably from the ones predicted 
 with this set of impedance factors. If we adhere to the restrictions in the 
 topology previously mentioned there are two methods which can be used to 
 calculate a new set of impedance factors. 
 
 Let us assume that we have the condition 
 
 or 
 
 Under this assumption it may be recalled that we will have for each branch 
 
 A^Q.. = A^Q.. = A 0.. 
 A Ji A ji A ji 
 
 Therefore 
 
 Av"!".At. 
 
 -7^-^^Ci*hi^W'^i (3.1) 
 
 ■15- 
 
-16- 
 
 After a complete numerical cycle all the terms on the left are known and 
 therefore the impedance factor Zj can be calculated. 
 
 To separate the three impedance factors two more equations are needed. 
 To obtain these we must first define an equivalent incremental capacitance, 
 resistance and inductance which are in series between point Qj and the reference 
 node : 
 
 AQ.. 
 
 ^^ (3 = 2) 
 
 ^ CEji 
 
 Av 
 
 R^. =-^ (3»3) 
 
 Ji 
 
 At. 
 1 
 
 1^=^^ (3A) 
 
 (At.)2 
 
 Combining these with the definition of Z ., Z^^. and Z we get 
 
 At. 
 
 ^RJ = % <3-'^ 
 
 Z, • = ^ '3.7) 
 
 For small changes in voltages, the equivalent elements C^ , R^ and L^ will be 
 constant , Hence 
 
 Z^. ^ At. (3-8) 
 
 Cj 1 
 
■17- 
 
 Z = constant (3»9) 
 
 ^Lj " i: (3.10) 
 
 From this we see that if we repeat the numerical cycle three times with three 
 
 different values of time increments k-|At., k^At. and k^At . we will be able to 
 
 1 1 2 1 3 1 
 
 get three values for each Z.. By using the variation of the impedance factors 
 
 J 
 
 rfith At. noted in equations 3-8 to 3«10^ we are able to write three equations 
 
 Z 
 
 k Z^. + Z^ — -^ = Z. m = 1, 2, 3 (3oll) 
 
 m Cj R. + k jm y > -J \-> J 
 
 //here Z. is the impedance factor found by using time interval k At.. When 
 jm mi 
 
 <:„ = 1.0, the solution of these equations yields 
 
 k k 
 Z. = - — -^ 
 
 L k^ - k3 
 
 "2 - \ VA' 
 . 1-k^- k3-l 
 
 (3^2) 
 
 ^c = -^^k^ (3a3) 
 
 Z^= \ -^L-^c (3.1V 
 
 The condition AQ . » A. Q, . can be realized by injecting into the 
 .op node a pulse which is much higher than the normal signal. Since this will 
 i'esult in large changes in voltages, all nonlinear elements must be temporarily 
 •eplaced by a linear equivalent whose parameter values are the same as the 
 ncremental parameters of the nonlinear element during the last normal cycle » 
 tecause of these much higher than normal responses, one of these impedance 
 actor calculation cycles cannot be used as a normal cycle in which time is 
 ncremented. 
 
-18- 
 
 The second method of determining the impedance factors also depends 
 on the assumed topology. If the circuit conforms to this topology, the various 
 impedance factors are nothing more than the equivalent impedances looking from 
 the point in question toward the reference node,, If we recall that a node 
 impedance factor Z was already defined, we can use this formula to calculate 
 three impedances for the th_ree time increments k-j^At^, -^t. and k At^, We can 
 then again use equations 3.12 to 3„1^ to calculate an equivalent capacitance 
 C J resistance R^ and inductance L^ for the node„ Since these are in series, 
 they can be added to the corresponding elements in the branch above this node. 
 In this way all the impedance factors may be calculated in a more direct manner 
 than that of the first method. Since this last method is also much easier to 
 program and much more conservative of computation time, in what follows It will 
 be assumed that this latter method is being used to calculate the impedance 
 
 factors. 
 
 At this point we might also introduce the conditions required to 
 change the time increment during the process. We must impose the condition 
 that during the time increment change the voltages across the three impedance 
 factor elements defined in equations 3 = 2 to 3,^4- remain constant; 
 
 v,^ . = C,^ . Q . = constant 
 
 AQ. 
 
 V _ . = R^ . —77- = constant 
 PZ j ^ j At 
 
 2 
 
 v^„. = Lv- . — ~^ = constant 
 ^J ^J (At)2 
 
 Since C , R_ and L^ are also constant during the increment size change, we 
 conclude that the terms 
 
■19- 
 
 AQ. A^Q. 
 
 Q , -^ and ^ 
 
 J ^^ (At)^ 
 
 must also remain constant o This gives us our desired information, 
 
k. VOLTAGE- CONTROLLED SUBNODES 
 
 In the previous formulation the current- controlled nonlinear element 
 has a decided advantage since all the currents are determined before the branch 
 and node voltages are csuLculatedo In order to allow voltage -controlled devices 
 to be handled just as efficiently it is desirable to formulate a similar process 
 which will be compatible with the above but in which the role of the current and 
 the voltage are interchanged. In this case we will allow KCL to be violated. 
 
 In the development of the relevant formulas we will assume that we are 
 working with only one branch extracted from the circuit. To be as general as 
 possible, we can postulate that this branch is made up of a series string of RLC 
 elements as shown in Fig, k,l. The current generator beside each parallel unit 
 represents the error current for that unit. Since the nodes between the parallel 
 units are associated only with branch j and do not have "levels" assigned to 
 them, we will call them 'voltage -controlled subnodes" to distinguish them from 
 the normal nodes. To derive the necessary relations we must turn to the dual 
 concept in which the flux is defined by 
 
 cp = r v(t)dt (i^.l) 
 
 
 
 Then for linear elements we have 
 
 cp = LI 
 
 ^ = RI (i^.2) 
 
 dt 
 
 2 
 dfcp ^ 1 -. 
 
 dt^"^ 
 
 -20- 
 
-21- 
 
 T- i-^iiA^o, 
 
 ^ ERROR CHARGE PATH 
 
 |AOl 
 
 ym ^m 
 
 m .i^QR 
 
 m 
 
 i— . 
 
 SUBNODE m 
 
 
 I 
 
 '** .-^N 
 
 7- ^-^ ^'Om 
 
 r 
 I 
 
 Figure 1^,1. A Branch Consisting of M Voltage-Controlled Subnodes, 
 
-22- 
 
 If we go to the difference notation we get for the mth parallel unit 
 
 L 
 . „ in 
 cp . = AQ . — — 
 mi mi At. 
 1 
 
 ACp . = AQ .R (4„3) 
 
 mi ml m 
 
 ,2. .„ ^*1 
 
 ^ V = '"Sai -C" 
 
 m 
 
 where 
 
 Acp . = cp . - cp . , 
 mi mi m(,i-l; 
 
 ^2, 
 
 A^cp . = Acp . ~ Acp .. ^x {h.h) 
 
 mi mi m(i-l; 
 
 and 
 
 ^22 
 A^cp . = A'^cp . - A cp 
 
 mi mi m(,i-l; 
 
 These equations suggest that we define a set of "admittance factors" in a 
 similar manner as the impedance factors were defined: 
 
 PL cp . 
 
 m rai 
 
 f = ^^ (4.6) 
 
 PR Acp . 
 m mi 
 
 AQ, 
 Y. 
 
 m A cp . 
 mi 
 
 'Cmi (i^„7) 
 
-23- 
 
 ^Pm ~ ^PLm ^ ^PRm + ^PCm (^"^^ 
 
 where AQ^, AQ^, AQ^ are the charges flowing through the parallel L^ R and C 
 elements respectively. Also define 
 
 A^Q . = AQ.. - AQ . (i+ Q) 
 
 where 
 
 AQ = AQ . + AO^ . + AQ^ . 
 mi Xjini ^mi Cmi 
 
 2 
 A Q corresponds to Av . . in the main method. 
 
 Let us suppose that during the time interval (i-l) an error of 
 
 2 
 ^ "^mfi-l) ^^sulted at subnode m. The condition for the compensation of this 
 
 error in this step requires 
 
 ^2, , , 
 
 •1 
 
 The amount of flux A^cp . required to do this is given by 
 
 3 1 
 
 Pm 
 
 _^ ^m(i-l) - ^PLm(%(i-l)) - (YpLm ^ ^PRm^^^fi 
 
 (i-l) 
 
 (^ai) 
 
 This correction formula is analogous to that of equation 21-6 and hence can be 
 treated with the same type of numerical procedure. Furthermore;, any future 
 comments on equation 2.16 can therefore also be applied with minor changes to 
 equation ^,11. Further derivations parallel to those of sections 2.2 and 2„3 
 could be continued here which would then make it possible to work directly with 
 such strings of parallel elements as shown in Fig. 4.1. However, in our method 
 we can be content with just being able to handle one set of parallel elements. 
 
-2k- 
 
 The subscript, m will be retained to distinguish charge in the voltage -controlled 
 subnode from the other charges. 
 
 Besides equation i4-.ll we also need a method of interconnecting the 
 voltage -controlled node with the rest of the nodes of the circuit. To obtain 
 this let us suppose that branch j which contains the voltage- controlled subnode 
 has a set of impedance factors Z^., Z^. and Z^^ defined for it and is regarded 
 by the charge distribution method of section 2.2 as a normal branch. After the 
 charge distribution cycle a charge AQ has been assigned to flow through 
 branch j . The increase in charge 
 
 ^ Si = ^Si - ^S(i-i) 
 
 can therefore be treated exactly the same as the error charge A Q^^_-l) and 
 hence the two can be summed to form a total effective error charge of 
 
 ^Sn(l 
 
 P 2 
 
 1-1) ^m(x-l) Ji 
 
 (i^.l2) 
 
 Hence the increase in <P is determined by 
 
 a3cp. 
 
 1 .2 
 
 mi Y. 
 
 Pm 
 
 [4«»(1-1) - ^Pl(^^(1-1)' - (^PL * ^PR)^\(i-l), 
 
 (l».13) 
 
 Since 
 
 Acp 
 
 V 
 
 mi 
 
 mi At, 
 
 we get the prediction formula 
 
 ^Sii(i-l) ^PLm ^^ 
 mi m(i-l) m(i-lj Y.^At 
 
 Pm i 
 
 ^Pm '"(l-^' 
 
 Y + Y 
 PLm PRm 
 
 ^Pm 
 
 Av 
 
 m(i-l) 
 
-25- 
 
 where 
 
 %i = \± - \il-l) {h.iy 
 
 With this formula the voltage across the parallel elements can be calculated 
 and consequently the currents through them can be obtained. The total charge 
 %i ^^"^ therefore be computed and the error charge for this cycle becomes 
 
 A Q^. = AQ - AQ (^^,6) 
 
 After these calculations have been performed for all voltage-controlled subnodes, 
 
 all branch voltages are known and therefore the node voltage calculation part 
 
 of the numerical cycle can be started. 
 
 From the above it should be apparent that branch j which contains the 
 voltage-controlled subnode can also have the full complement of series elements 
 
 in addition to this subnode. This results because we essentially have two 
 error-correcting mechanisms associated with branch j and therefore it can have 
 two sets of elements, one in series and another in parallel. 
 
 By the use of this voltage-controlled subnode technique, restriction 
 two of section 2.1 (only one branch can be connected between two nodes except 
 if one of the nodes is the reference) is eliminated since we now have a method 
 for dealing with parallel branches without violating the first restriction. Of 
 course this technique is also useful for voltage-controlled elements since the 
 voltage across the elements is first predicted and the currents are calculated 
 later. 
 
5, DISCUSSION OF CORRECTION FORMULA 
 
 In section 2.2 we derived equation 2,l6 which will compensate for a 
 residual error between the branch and node voltages from the previous step. 
 This formula was developed so that a circuit of the assiomed topology and with 
 correctly defined impedance factors will respond so as to have a zero resiilting 
 error between v . and vt. . If the impedance factors are not completely correct 
 due to the presence of nonlinearities, equation 2.l6 will not make the proper 
 correction. 
 
 5„1 Error Term for a Branch 
 
 To find the factors related to this correction formula, let us examine 
 the error that will be produced when the impedance factors are not correctly 
 defined. To be able to do this we make the following assumptions: 
 
 (1) The branch considered has impedance factors such that 
 V is relatively independent of v . 
 
 (2) Z', the impedance factor actually used, and Z, the 
 correct impedance factor, are related by 
 
 Z' = Z + AZ 
 
 where AZ is the error in Z' . Similar definitions hold 
 for 
 
 ^'C \' ^C ^R' ^R' ^R' ^' ^L ^^ ^ 
 
 2 
 (X) Q ,, AQ. -, and A Q. , are known. 
 
 ^~" 1-1-^ 1-1 1--L 
 
 iS) The error voltage after the last step in the numerical 
 process was 
 
 -26- 
 
-27- 
 
 (5) The node voltage v will be increased by an amount 
 
 Av, .. 
 ki 
 
 The IQ method will then produce the terms 
 
 and 
 
 4«i - f^ hl-l) - ^C ^ - (^H ^%^^] (5.1) 
 
 Av At 
 
 a3q' = a3q^ + a3q^ (5„3) 
 
 The correct third difference of charge should be 
 
 ^%=fK-X-c%^-(^C-K)%i| (5..) 
 
 4^ '- -^ (5.5) 
 
 The error is 
 
 a3q. = a3q. - a\' (5^6) 
 
 By the substitution of appropriate definitions this b 
 
 ecomes 
 
 3 AZ ^^n - ^r^ Z^ - Z^T o 
 
 ^\ = Z(^Taz) (^\-1 - ^\i)^t . -^xfr^ ^^i-l - Z(Z.AZ) ^^^i.-l 
 
 (5.7) 
 The new error voltage for step (l) is given by 
 
-28- 
 
 Av. = V, . - V. 
 X kx 1 
 
 By using previously derived relations we get 
 
 Av^ = Z^ (5.8) 
 
 Hence 
 
 ZAZ^ - Z AZ AQ. , Z^AZ - ZAZ^ A Q. -, ' 
 AZ /, , N C (J x-1 , _L L ^x-1 
 
 ^^i = zVaZ ^^^-1 ^ ^""ki^ ^ Z -kAZ aT "■ Z + AZ At ; 
 
 (5.9) ;| 
 
 The first term indicates that if AZ is small, the error due to i^^^^j^ + ^^ki^ 
 will be negligible. The second and third terms are directly related to the 
 current and the current increase respectively in the branch. It may be obserA-ed '1 
 that these can have much larger magnitudes than the first term. Nevertheless 
 they do exhibit the convenient property of being slowly varying functions of 
 time. This is a consequence of assuming that the impedance factors are rela- 
 tively constant even thovigh a nonlinear element may be present and that the 
 current will normally not vary greatly from one numerical step to the nexto 
 
 5.2 Compensation of the Error 
 
 The error caused by the last two terms in equation 5-9 is mostly 
 attributable to the fact that the impedance factors were developed so that the 
 circuit responds accurately to a pulse of width At which for linear elements 
 connected in the assumed topology may be of any height. If a step function is 
 applied, the circuit will respond accurately initially, but as the effects of 
 the transient decrease, new impedance factors corresponding to longer At's should 
 
-29- 
 
 be used and eventually when the circuit is again in a steady state the 
 impedance factors should merely he defined from the resistances of the circuit. 
 
 We note that as the effects of the transient are dying out the cur- 
 rents will become more and more constant from one step to another „ Therefore 
 the terms 
 
 ZAZ^ - Z^ AQ._^ Z^AZ - ZAZ^ A^Q. 
 Z + AZ aT + Z + AZ A^ 
 
 may be compensated by the introduction of a slowly time-varying factor into the 
 correction formula. Since the error impedance factors of the form AZ are not 
 known by the numerical process, this factor would have to be adjusted by 
 Observing the trend of the error voltage Av. and making changes to minimize Av 
 in the next few steps. The introduction of this adaptive variable eliminates 
 the need of a multiplicity of impedance factors which would be impractical to 
 use. 
 
 In the correction formula 2.l6 we replace the term Av, , bv 
 
 (i-1) -^ 
 
 ^""(i-l) + \ ''^^^^ \ is the new factor mentioned above. Whenever the circuit 
 is in a steady state, then at each branch 
 
 A^Q. » 
 
 1 
 
 Av. ~ 
 
 1 
 
 This requires that 
 
 ""a " ^C ^ (5.10) 
 
 When the transient analysis is started this condition must be forced and during 
 transients v^ should be adapted so that when a steady state is reached this 
 
■30- 
 
 condition will again hold. By the introduction of this correction factor the 
 response of the circuit to the transients is in no way impaired; the impedance 
 factors defined for At look after the fast transients and v^ adjusts the circuit 
 to the slower transients. 
 
 It should be emphasized at this point that the error voltage Av^ at 
 each branch is known after each numerical step and hence can be controlled „ If 
 an error at a branch is too large, one can always repeat the step using a 
 
 smaller time increment. Since Z varies directly with At and since the current 
 
 AQ. 
 
 — - remains constant during the step size change the error term 
 
 At 
 
 AQ. 
 ^C At 
 
 will also decrease with a decrease in the time interval. Hence the introduction 
 
 of V is not absolutely necessary, but its use will increase the speed of 
 a 
 
 computation considerably. 
 
 5,3 A Predictor-Corrector Type Method 
 
 To improve further the properties of the IQ method we can introduce 
 
 a predictor -corrector type procedure which is somewhat analogous to those used 
 
 [6] 
 
 in the integration of differential equations 
 eqixation 2„l6 as 
 
 We redefine the correction 
 
 4^--f 
 
 Av 
 
 m 
 Ei 
 
 'C At 
 
 ^\ ^ ^r) 
 
 At 
 
 m = 1, 2, 00,, M 
 
 (5.11) 
 
 where M is the number of times the process described in section 2 is repeated 
 for each increment in time. The procedure is to calculate a value for Av^^ 
 for each node, use this in the distribution of the currents and calculate all 
 
-31- 
 
 the node voltages. An error voltage Av. will result which is then used to define 
 
 2 
 a second value ^v^^. This cycle is then repeated M times^, the last "being the 
 
 one which will he accepted as yielding the charges and voltages for this time 
 increment . 
 
 The problem is to determine a set of expressions for the terms Av"^ . 
 
 Ei 
 
 This, in the general case, could become quite complicated but the results could 
 
 prove to be valuable here in the same way as those obtained by including more 
 
 terms in the predictor formulas in the integration of differential equations. 
 
 To simplify matters we will restrict the discussion to the case in which M = 2, 
 
 a value which would probably be close to the optimum when computation time 
 
 versus decrease in error is considered. The expression for Av„. then has the 
 
 El 
 
 form 
 
 Av^. = V + a Av. -, 
 El a 1 i-1 
 
 a^ can be either a constant for the branch or an adaptable factor which can be 
 slowly adjusted during the analysis so that 
 
 Av! ~ 
 1 
 
 2 
 The expression for Av„. becomes 
 
 El 
 
 Av2. = ^^ ^ a^Av._^ + a^Av^ 
 
 where a^ is again a constant or adaptable factor. 
 
 Since the three variable v , a and a all are directed towards 
 decreasing the error Av^, it is obvious that they are highly interrelated. 
 But, the role that they play is different: v is related to the algebraic sum 
 of the error over a number of time intervals and a and a are related to the 
 
-32- 
 
 numerical oscillations and damping. Since these parameters must be varied 
 slowly, it is apparent that the numerical process must be slow enough to allow 
 the adaptation of these parameters to the response. 
 
 This predictor-corrector type procedure was not studied further, but 
 an analysis of the factors introduced and their relation to stability would 
 probably prove useful. 
 
6c ERRORS AND STABILITY 
 
 6.1 Circuit Error and Response Error 
 
 Up to the present we have been talking about an error A^oltage Av. 
 at each branch or an error current A Q./At at a voltage-controlled subnode. 
 For purposes of discussion we will call these the "circuit" errors » There is 
 yet another type of error which we will name the "response" error,, The latter 
 refers to the difference between the theoretical value of a current or a voltage 
 at a given time and the value of the same current or voltage obtained by the 
 numerical process at the same time instance. 
 
 Previously it was shown that the circuit error can be controlled by 
 
 specifying a maximum allowable error between any branch and the node to which 
 
 it is connected. If we have K nodes connected in the worst case obtainable^ 
 
 namely all in series, and if at each node there is an error Av and since the 
 
 max 
 
 voltage at the branch connected to the highest level node must be within + Av , 
 
 — max 
 
 we can conclude that the circuit error of any node voltage is less than — Av 
 
 2 max 
 
 Normally this error would be much smaller and in the specific case the exact 
 
 upper bound can be determined from the topology of the circuit. We can also 
 
 conclude that whenever the circuit is in a steady state the response error will 
 
 be less than this maximum circuit error. 
 
 During transients the last statement does not hold, but we can define 
 
 a rough estimate of the upper bound for the response error during the transient 
 
 by considering the special case in which a linear capacitance is one of the 
 
 branches below the node whose response error we wish to find. Suppose that 
 
 the theoretical current flowing through the capacitance is given by l(t) and 
 
 AQ^ 
 that the current calculated by the IQ method is — r- where 
 
 iAt < t < (i + l)At 
 
 Also define 
 
 -33- 
 
-31^- 
 
 1 
 
 Then the error in the charge on the capacitance is 
 
 ^^ m=l m=l 
 
 If we assume Av. = (ioe,, v = v.), the response error becomes 
 
 1 i^l 1 
 
 1 . „ 1 ^ .2, 
 
 ^^E1 = C^\ = C ^AS 
 
 -, ^ m 
 
 m=l 
 
 But the maximum circuit error possible is Av and consequently 
 
 max 
 
 2 
 ^^i 
 C max 
 
 If there is no cancellation of the circuit error from one time interval to 
 another, we get the maximum response error 
 
 Av^. = iAv 
 Eimax max 
 
 In practical situations this would be much lower and fiirthermore a much closer 
 
 upper bound for a specific circuit can be calculated by observing the circuit 
 
 error at each time interval during the transient. From this it is also apparent 
 
 that the adaptive factors v and a^ should be adjusted so that small numerical 
 ^ a 1 
 
 oscillations will result. This will cause the circuit error to oscillate 
 between positive and negative values and therefore produce a much smaller 
 response error. It is also felt that a much lower upper bound for the response 
 error could be arrived at if the energy conservation criterion is used to cal- 
 culate the node voltage, but this was not examined. 
 
-35- 
 
 6,2 stability 
 
 The stability of the IQ method has not yet been examined in detail, 
 but from the discussion of errors we can make a few observations about 
 stability. First let us recall that in the normal predictor-corrector methods 
 of numerical analysis there is available an error term which is used to control 
 the increment. Consider the example 
 
 Predictor x . = x , ^^ + Atx , ^ 
 Pi c(i-l) c(i-l) 
 
 Corrector x . = x , , . + — (x + x , ) 
 
 ci c(i-l) 2 ^V ^ ''p(i-l)^ 
 
 where the term (x^. - x^.) is the error. It is to be noted that this is a 
 relative error since the corrected as well as the predicted value may be 
 erroneous. When both terms are in error but the difference between them is 
 small, the solution will possibly become numerically unstable. 
 
 On the other hand, in the IQ method the error which controls the step 
 size is the circuit error, an absolute error. We have already seen that we can 
 set an upper bound for this error. Since we can assume the circuit error is 
 initially within the allowed bounds, then for each succeeding step the error 
 will also be within these bounds. The last statement is justified if we assume 
 that there are no discontinuities in any v-i characteristics and if we recall 
 that a decrease in the time increment implies a decrease in the circuit error. 
 Therefore as long as we have continuous functions for the v-i characteristics 
 the solution will not be able to "blow up." 
 
 The problem of instabilities in the integration of the differential 
 equations of a circuit are usually related to the minimum time constants of 
 the circuits which will cause instabilities if the time increment is too 
 
-36- 
 
 largeo^"'"-' For the IQ method this does not appear to he a problem. Consider 
 for example the tranch in Fig. 6.1 which has a very low series inductance and 
 a high resistance: 
 
 L. 
 Lj At 
 
 Rj J 
 
 Z ~ 
 
 Figure 6.1. 
 
 If the node voltage varies slowly with time and a large At is heing used then 
 Z is insignificant compared to Z^^ and consequently Z^ is all but ignored by 
 the method. If a large transient (with a finite slope) appears at the node. 
 At will be decreased until the node voltage change is only incremental. Now Z^ 
 has a much higher value and might even dominate Z^. As the transient disappears. 
 At is again increased in size and Z^ again assumes an insignificant role. Thus 
 we see that the small inductance will only have an effect on the circuit during 
 fast transients. 
 
7. RELAXATION OF TOPOLOGICAL RESTRICTIONS 
 
 Suppose that it is desired to perform an analysis of a bridge type 
 circuit shown in Fig. 7»la. The elements marked L are assumed to be lower 
 impedances (e.g., 10 ohms) than those marked H (1000 ohms)o Two possible 
 redrawings of the circuit according to the method described appear in Figs. 7„lb 
 and 7»lc. Both forms of the circuit can be analyzed by the IQ method but for 
 fast transients considerably more difficulty will be experience with Fig„ 7olc 
 than with Fig. 7.1b. It may be recalled that the impedance factors were 
 derived from the "down-directed" equivalent impedances. The impedance factor 
 at the top node would therefore be approximately kO ohms for Fig. 7 .lb and 
 about 500 ohms for Fig. 7.1c, whereas the correct impedance is a little less 
 than kO ohms. If a voltage step is applied to the circuit of Fig. 7olc, the 
 injected current will be grossly underestimated and the process of establishing 
 equilibrium would take much longer than in the case of Fig. 7.1b where it would 
 be established in very few numerical cycles. 
 
 The conclusion from the above is that if a circuit does not have the 
 desired topology a wise rearrangement of the circuit will possibly make it 
 amenable to as fast a solution as that of a circuit which has the restricted 
 topology. There seem to be a few general rules which will usually lead to 
 convenient configurations if the original circuit does not have the desired 
 topology: 
 
 (1) If possible the highest level node influencing the 
 portion of the circuit with the undesirable topology 
 should have as low an impedance path to the reference 
 as obtainable . 
 
 (2) If possible the circuit should be arranged so that the 
 lowest impedance branches lead to the lowest possible 
 node levels . 
 
 -37- 
 
-38- 
 
 a) ORIGINAL 
 
 k3 
 
 L < R. 
 
 L<R. 
 
 R3<H 
 
 L<^K 
 
 Hi 
 
 L>R, 
 
 b) CORRECT REDRAWING 
 
 
 lSR, 
 
 L>. 
 
 R^^H 
 
 C) WRONG REDRAWING 
 
 Figure 7„lo Example of a Circuit Violating the Assumed Topology. 
 
-39- 
 
 By allowing these variations in the topology, we must also restrict the 
 forcing function voltage to incremental changes from one time interval to the 
 next since the impedance factors will no longer give the exact response to a 
 high pulse. In a specific circuit the amount of difficulty the method will 
 encounter can be predicted since it is directly related to how well the 
 impedance factors approximate the actual equivalent impedances. 
 
8o NONLINEAR AND TIME- VARYING PARAJvETERS 
 
 8.1 Soiirces 
 
 A voltage source is already Included in the IQ method and a constant 
 current source does not present any difficiilties since this requires only that 
 a charge I^At be subtracted from one node and added to another before the 
 start of the charge distribution cycle. With constant sources there are no 
 problems (except possibly during the initialization) but if they are time- 
 dependent or proportional to currents or voltages in another part of the circuit 
 the method might have trouble. If a voltage soaarce, for instance, is time- 
 dependent and during one numerical cycle experiences a large change, this 
 change will not be taken into accoiint by the first charge distribution cycle 
 since no impedance factors are changed. When the node voltages are then cal- 
 ciilated there will be a large error at the branch with the voltage source. 
 Since this error will then be diminished by the next numerical cycle through 
 the error-correction mechanism, we see that the response of the circuit to all 
 varying voltages will be delayed one cycle, A current source will also produce 
 the same result, but with the slight advemtage that the node to which it is 
 connected will respond immediately and the delay will appear in a higher 
 level node. This advantage can often be put to good use. The rule in general 
 seems to be that only incremental changes should be allowed from one time 
 interval to another for varying sources, 
 
 8.2 Nonlinear Elements 
 
 In general we can specify the characteristics of a nonlinear element 
 
 by 
 
 -^^ = f (v) or V = g \-^^ m = 0, 1 or 2 (8,l) 
 
 (Atr \(Atr/ 
 
 -ko- 
 
-41- 
 
 The second form is the one which falls naturally into the domain of the IQ 
 method but unfortunately it is not always the most convenient to use„ For 
 instance, in a tunnel diode the relation 
 
 is multivalued. In such cases the voltage-controlled node formulas developed 
 earlier can be used in conjunction with the first of the above form,ulas which 
 is in this case single valued. The important thing is that the prediction and 
 circuit-error correction mechanism is supplied by the method, therefore 
 avoiding the need of a repetitious process such as Newton's method to determine 
 the independent variable when the dependent variable is given. This can prove 
 to be an important time-saving factor during computation. 
 
 80 3 Semiconductor Diodes 
 
 Since the diode and transistor are the most important nonlinear 
 elements, they will be discussed here in more detail. The dc characteristics 
 of an ideal diode are given by 
 
 l = lQ[e^-l] A. = ^ (8.2) 
 
 q 
 
 The Incremental conductance is given by 
 
 V 
 
 (8,3) 
 
 ^D dv X 
 During the transients the incremental diode capacitance given by 
 
 ^ " dv 
 
 (8.i|) 
 
-42- 
 
 also becomes important « When in the nonconducting region the diode exhibits 
 a depletion layer capacitance of the form 
 
 C = -^ (8„5) 
 
 where v is the built-in potential, n is the grading constant which is usually 
 D 
 
 1 1 
 equal to — or - and C-. is a constant. 
 
 The diffusion type capacitance when in the conduction region is 
 
 C =^^==g ^ (8.6) 
 
 D dv ^D dl ^ ^ 
 
 The term — ^ is a time constant which can usually be assumed constant as a 
 dl 
 
 first approximation and equal to the diode constant T Hence 
 
 C ~ g^T^ (8.7) 
 
 Normally only one of the two capacitances will be important in a particular 
 application. Tn switching circuits C is most important and usually C^., is small 
 enough that it can be approximated by a constant. 
 
 If T cannot be considered a constant the relation Q = f (v) must be 
 obtained by niimerical integration during the transient: 
 
 Cj^(v)dv=j gj^(v)Tj^(v)dv (8,8) 
 
 
 
 If T can be assiimed constant we get the relation 
 
 v 
 Q = .,l/ (8.9) 
 
43- 
 
 When the depletion layer capacitance parameters v^ n and C can be assumed 
 constant, and if the diode is reverse "biased, we get the formula 
 
 Q = 
 
 n - n *- 
 
 
 (8a0) 
 
 From a programming point of view the expressions for 1, g , C and Q are in a 
 very convenient form since they are so highly interrelated » This gives the 
 very important result that no further simplifying assumptions which are usually 
 made in a normal analysis will be required. 
 
 To complete the model for the diode we can include the leakage 
 resistance R^ and if the reverse breakdown is important, a current of the form 
 
 ^B = 
 
 '0 
 
 1 - (^)- 
 
 (8.11) 
 
 where 
 
 1 
 
 1 - (^)- 
 ^B 
 
 is the multiplication factor and v_ is the reverse breakdown voltage o 
 
 Since all the above elements are in parallel and since the current 
 equation is in the form 
 
 
 we must use the voltage-controlled subnode methods The relevant parameters 
 are given by 
 
-kh- 
 
 ♦ * 9 
 
 > ^dS 
 
 ^D ^C, V 
 
 a) IN TERMS OF DIODE PARAMETERS 
 
 . AO: 
 
 6\ J 
 
 b) IN TERMS OF ADMITTANCE FACTORS 
 
 Figure 8.1<, Semiconductor Diode Model. 
 
c 
 
 (8.12) 
 
 y — ^ "V — 
 
 PC At P ~ Sj- + ^ + ^ 
 
 c 
 
 The prediction formula h.lk therefore reduces to 
 
 4^(i-l) ^PR 
 -i = -i-1 - ^-(i-1) - -Y^ - -f AV( ._^) (8„13) 
 
 The spreading resistance and the wiring inductance have not yet been included^ 
 hut these present no prohlem since we can still specify a complete set of 
 series elements for the branch containing the diode. 
 
 Q.h Transistors 
 
 Of the different large signal models available for the transistor j, it 
 was found that for the IQ method the Ebers and Moll model gives the best 
 balance between accuracy and efficiency in computation. For this reason the 
 technique available to handle this model will be developed in detail. 
 
 The equations of a PNP transistor for the dc case are 
 
 I ^ T ^ 
 
 ^ = -r^(e- -i)..^^^(e- .1) (8.ia) 
 
 I ^ T ^ 
 
 T - ry EO , kT ^ s CO , kT , . /o . n 
 
 ^c - ^n 1 -a^^ (^ - 1) - r^^^ ^^ - ^^ (Q"^5) 
 
 The notation and the sign convention is that used by Ebers and Moll^^ with the 
 exception that the symbol 9 is replaced by v. Therefore these will not be 
 
-k6- 
 
 discussed here in detail. For future convenience let us make the following 
 definitions °, 
 
 X = 
 
 kT 
 
 -1 
 
 ■^ 
 
 EO 
 
 1 - a 
 
 ^1 
 
 Therefore 
 
 E 
 
 ^=^o(- -^' 
 
 -I 
 
 CO 
 
 'CFO 1 - ajXt-j- 
 
 ^CF = ^CFO^^ - ^^ 
 
 ^ = ^r - "-^' 
 
 HF I CF 
 
 -■-C '^N-'-EF "^ """CF 
 
 (8a6) 
 
 A conductance for the emitter junction and one for the collector jiinction 
 are defined as follows °, 
 
 DTfl 
 
 S 
 
 E 1 5v, 
 
 "SiF "^ "S:: 
 
 FO 
 
 V =const 
 c 
 
 (8,17) 
 
 81, 
 
 gn = 
 
 C \5v. 
 
 ■'"CF '^ "'"CFO 
 
 Vw= const 
 ill 
 
-47- 
 
 When the junctions are reverse biased the leakage resistance R and R for 
 
 the emitter and collector junctions respectively become important and therefore 
 
 must be included in the equivalent circuit. 
 
 During transients either a diffusion capacitance or a depletion layer 
 capacitance come into play at each junction, depending on whether the junction 
 is forward or reverse biased. The two diffusion capacitances which are char- 
 acterized by two time constants T^ and T^ are given bv 
 
 E C "^ 
 
 C = T a; 
 DE E^E 
 
 ^CE ~ ^C^C 
 
 (8.18) 
 
 The time constant T^ and T^ can be assumed constant as a first approximation 
 and therefore similar charge versus voltage relations as equation 8„9 for 
 diodes can be obtained. The depletion layer capacitances are given by 
 
 o - °o= 
 
 TE / v„ \n^ 
 
 11-^ 
 
 V / 
 
 DE/ 
 
 (8,19) 
 
 c - 'oc 
 
 " ('-%r= 
 
 where C^^ and C^^ are constants of the transistor, v^^ and v are the built-in 
 potentials for the emitter and collector junctions and n^ and n are grading 
 constants. 
 
 Combining all of the aboA^e we have the model shown in Fig, 8,2a , 
 The reason for the unusual orientation will become apparent later. We note 
 two important properties of this equivalent circuit: the collector and the 
 emitter junctions have the same configurations and similar elements as the 
 
-kQ- 
 
 c 
 
 H 1 
 
 BO- 
 
 [ICF 
 
 Rr> ^cS 
 
 DC 
 
 'TC 
 
 'E S-> sA 
 
 !■. 
 
 'OE 
 
 'TE 
 
 E 6 
 
 ^^I^EF 
 
 (T I 
 
 N CF 
 
 a) IN TERMS OF TRANSISTOR PARAMETERS 
 
 E 6 
 
 b) IN TERMS OF ADMITTANCE FACTORS 
 
 Figure 8.2. IQ Model of a Transistor, 
 
-49- 
 
 equivalent circuit for the diode except for the current generator, and the 
 circuit is in the exact form required by the voltage-controlled subnode method 
 of section k. The first of the above indicates that we could easily use the 
 method developed for diodes in section 8.3„ But, the second property suggests 
 that a m.ore general technique designed specifically for transistors can be 
 formulated „ 
 
 Since we have no inductive element in the circuit, equation 4,11 
 reduces to 
 
 a\ 
 
 mi Y^ 
 Pm 
 
 Y 
 PRm „2. 
 
 .4Sr.(l-l)J - tf "^Vi-D <^-'°' 
 
 Pm 
 
 To put this in terms of voltages we invoke the relations 
 
 V 
 
 At ^^ = -A^ ^^ = ~^ (8-21) 
 
 Therefore 
 
 2 
 
 ^ Vi - -Y-aT- - ^7 ^%(i-l) (8»22) 
 
 Pm P ^ ^ 
 
 For the transistor we have two parallel units, unit 1 at the collector junction 
 and unit 2 at the emitter junctiono The admittance factors for these two units 
 are defined by 
 
 PRl ^C R 
 c 
 
 „ Sc ^DC /ON 
 
 ^PCl = -A^ -^ ^ (8-23) 
 
 Y = Y + Y 
 PI PRl PCI 
 
-50- 
 
 Y = a: + — 
 PR2 ^E Rg 
 
 PC2 At At 
 
 Y = Y + Y 
 P2 PR2 PC2 
 
 The voltages are given by 
 
 a2 
 V-, . = V . = V, / . , N + Av-, / . .X + A V 
 li ex l(i-l) l(i-l) li 
 
 (8.25) 
 
 ^2i ^ -^Ei ^ ^2(1-1) "^^2(i-l) ^^^2i 
 
 V . = -v^ . = V^ / . . \ + Av„ / . , ^ + A V^_. 
 
 The two current generators a 1 and « I can be absorbed into the 
 "error" current generators A^Q^/At and A^Q^/At^ It may be recalled that A Q^ 
 was the error charge as defined in equation 4„9 and AQ was the sum of the 
 charges flowing through the parallel R, L and C elements. For the transistor 
 the relations become 
 
 a\ = AQ^^ - AQ^ 
 
 A^Q^ = AQ^^ + AQ^ - AQ^ 
 
 where AQ and AO^^ are defined in Fig, 8,2b and AQ^ and AQ^ are given by 
 
 v^ Av^ Av^ 
 
 
 V Av Av 
 
-51- 
 
 There is still one remaining inconA/enience which concerns timing., 
 If the normal IQ procedure is followed, the current will first be distributed 
 and therefore at the transistor we get the terms 
 
 '^CTl = '^CT(l-l) + ^^%Tx 
 
 V-^(l-l) ^^\l 
 
 The error term which is to be used in equation 8„22 then becomes 
 
 ^^li = ^'^Ki-D - ^'^cTi 
 
 4^i=^%(i-l) --^'Si-^^^Tl 
 
 8.28) 
 
 If the new voltages v^^ and v^^ are calculated from this, there is a possibility 
 of a larger error in AQ^ and AQ^ since the correction formula does not take 
 into account the fact that OC^I^^t and aj^I^^^At will have different values in 
 time interval i from those in time interval (i-l)„ Consequently the two cur- 
 rent generators will always be lagging behind the voltages by one numerical 
 
 cycle during transients „ To compensate for this, values for v, and v can 
 
 li 2i 
 
 be predicted from errors in Interval (i-l) and these may then be used to predict 
 A^alues for 0!^I^.^t and Qj^IgpAt for the time interval i„ The changes from 
 interval (i-l) to i can then be included in the error formulas 8o28o It is 
 also apparent that a predictor-corrector type method such as described in 
 section 5,3 could be used here to further reduce this numerical lag effects 
 
 To see how the transistor equivalent is incorporated into the rest of 
 the circuit we must recall the assuro,ption that only one branch of the main cir- 
 cuit is involved when dealing with a set of voltage-controlled subnodes. Hence 
 the entire equivalent for the transistor except the base is included in one 
 
-52- 
 
 branch. The collector bulk resistance R and any possible lead inductance 
 
 cc 
 
 can, as with the diode, be specified in the series elements associated with 
 this branch. A second branch is required for the base and therefore it would 
 contain the base resistance R,-, o Since the voltage-controlled subnodes such 
 as that in the center of the transistor model are normally not available to 
 the main IQ method, a special set of nodes, one per transistor, must be 
 defined so that the channeling of the charge which flows through the base 
 branch to the internal transistor subnode can be accomplished. 
 
 The orientation of the transistor equivalent shown in Fig„ 8.2 is 
 mainly for purposes of being able to determine the approximately correct 
 impedance factors and to help prevent the current generator lags mentioned 
 above. The parallel admittances of the transistor equivalent can be converted 
 to series impedance factors by the method described in section 3« These are 
 illustrated in Fig. 8.3- Normally Z would be fairly small and capacitive 
 
 11 
 
 \i 
 
 Node below 
 emitter 
 
 Figure 8.3 
 
-53- 
 
 because of a large C^^, and Z^^ would be large except when the transistor is 
 in saturation „ Thus if Z^ is small the values determined by the impedance 
 factor calculation method of section 3 will be 
 
 ^2'\*h^ "CT " ^cT (8.29) 
 
 S'^l'\' \o (8.30) 
 
 Because of the normally high Aralue of Z^^ as compared with R , the defined 
 impedance factors are a good approximation of the actual AQ versus Av behavior 
 at the transistor terminals. Therefore this orientation of the transistor 
 equivalent will normally be used„ As an example of how a transistor circuit 
 would be redrawn before analysis by the IQ method, an asymmetric flipflop 
 circuit is shown in Fig.. 8„4, 
 
 8„5 Transmission Lines 
 
 When working with high-speed circuitry one invariably meets a situation 
 in which a signal delay is involved. Since an equivalent transmission line will 
 usually represent the phenomenon adequately, a technique by which the IQ method 
 can handle transmission lines will be briefly discussed. 
 
 Consider the lossless transmission line of Figo 8.5a. It was found 
 that the m.ost convenient equivalent circuit is that shown in Fig, 8o5b. Z is 
 the characteristic im_pedance of the line and the forward and reverse currents 
 at the sending and receiving end are defined as follows: 
 
-54- 
 
 4+25 
 
 I — V^ 
 
 <.I.5K <,l.8 
 
 .+ R- 
 
 •ih-0--v^ 
 
 •Ih- 
 
 i* >2 5 
 
 ^- 
 
 ] 
 
 2K 
 -50 
 
 4.9< 
 
 <' 3 3K 
 . + 2 2> 
 
 2i •'^-/ 
 
 <> <► <► 
 
 ■vHi- 
 
 >860 "-3 
 
 4 <t 
 
 t 
 
 i 
 
 * 
 
 -5 
 
 2K< 
 
 T T 1 
 
 50 
 
 4,9K 
 
 •® 
 
 .860 
 
 ♦-* f- C9 
 
 ■CD 
 
 • •■ 
 
 ^ 
 
 2 5K 
 
 .8K< 
 
 ® 
 
 i^bb 
 
 J — 
 
 3.3K< <Jr 
 
 J. 5 I 25 I 31 I 25 I 2.2 | 25 [ ^ J. 
 
 i i i • i 4 i • 
 
 ® 
 
 Figure 8.4. Asymmetric Flipflop Example, 
 
-55- 
 
 U I. 
 
 I 
 
 ■o 
 
 >z< 
 
 CHARACTERISTIC IMP = Zq 
 VELOCITY OF PROPAGATION =V, 
 
 >Zl 
 
 T =L 
 Q) ACTUAL TRANSMISSION LINE 
 
 b) EQUIVALENT CIRCUIT 
 
 [I NODE 10 
 
 I 
 
 NODES 
 
 l'-©h a,p' 
 
 ? ? 
 
 I I 
 
 I I 
 
 I I 
 
 [* NODE 9 
 
 J 
 
 ws 
 
 ? ? 
 
 I I 
 
 I I 
 
 I I 
 
 NODE 5 
 
 ? 
 
 ? 
 
 NODE 4 
 
 I i 
 
 I I 
 
 i I 
 
 C) A TRANSMISSION LINE IN A CIRCUIT 
 
 Figure 8.5« Transmission Line and Equivalent 
 
■56- 
 
 3 
 1^=1^ at the sending end 
 
 ■p 
 T = T at the receiving end 
 F F 
 
 ■□ 
 
 I = I at the receiving end 
 
 g 
 I - I at the sending end 
 R R 
 
 (8.31) 
 
 We therefore have the following relations 
 
 13=1^. I^ and 1, = 1^^ . 1^ (8.32) 
 
 I^(t) = l|(t -T) 
 
 4(t) = l^(t 
 
 (8.33) 
 
 By imposing boundary conditions so that the reflection coefficients 
 are correctly defined we get 
 
 S ^S 
 
 R ' ^OS 
 
 i; = i: . I 
 
 (8.3^) 
 ^ = I? - I 
 
 R "T OR 
 
 q "P 
 
 From the previous set of definitions we see that I and I^ are the delayed 
 forward and reverse currents respectively and therefore are known quantities 
 at any time instance. If the equivalent impedances of the line at the sending 
 and receiving ends are treated as normal branches by the IQ method, I^g and I^^^ 
 
 will be known after the current distribution cycle. Consequently the currents 
 
 S R 
 I^ and 1^ which are to be entered into the two delay tables are easily calcu- 
 
 F R 
 
 R S 
 lated. These emerge after a time T as the values for I^ and I^ respectively. 
 
-57- 
 
 The above can also be adapted to lossy transmission lines^ the 
 difference being that now the data in the delay tables must be processed in 
 some manner Instead of being merely stored. Since this is out of the domain of 
 this report, the subject will not be discussed further. 
 
 Although the above relations for the transmission line are very 
 simple;, two complications enter during implement ati on „ Since the delay tables 
 must be of fixed size there is a problem of matching the fixed At of the line 
 to the variable time increment which is possible in the IQ method. This problem 
 is fairly easily solved by integrating the charge entering the line when the 
 line At is larger than the IQ time increment and by "filling" the spaces in 
 the table which are not used when the line At is smaller than the IQ time 
 increment. The second problem is one of initializing the delay tables. Since 
 points 1 and 3, and 2 and k in Fig. 8.5 are parts of the same node at dc, they 
 should also be treated by the IQ method as the same nodes but this presents 
 some problems which have not yet been completely resolved and which are directly 
 related to how the main IQ method is programmed. However, this problem also 
 does not seem to be insurmountable. 
 
9. REALIZATION OF THE IQ METHOD 
 
 9.1 The Program 
 
 To verify some of the ideas contained in the above, a program based 
 on the IQ method was written. Although the time element prevented the inclusion 
 of all of the techniques developed,, the basic ideas of sections 2 and 3 were 
 all included. In addition a method for handling nonlinear two-terminal devices 
 such as tunnel diodes and the method for analyzing transmission lines were 
 incorporated into the program. Since some major changes in the organization of 
 the program would be required to include efficiently the voltage-controlled 
 subnode method^, the addition of this and hence the model for the diode and 
 transistor has not yet been completed. The program was written in FORTMN II 
 for the IBM 709^ system at the Digital Computer Laboratory, University of 
 Illinois , 
 
 The input part of the program is written in such a manner that the 
 sequence in which the nodes are read automatically defines their level. For 
 each node the numbers of all the branches below the node are also read from the 
 tape. The branch data read in by the program consists of the node number below 
 the branch, the four constants of the branch and if it has a nonlinear element, 
 the type and number of the nonlinearity. The rest of the data is read by the 
 program in a standard manner. 
 
 The main part of the program consists of routines for inputting and 
 outputting data, incrementing circuit variables and for the actual calculation 
 of the response, A simplified flow chart of the latter is shown in Fig, A-1 
 of the appendix. After initializing of the table, the dc initial conditions 
 are evaluated. Since this uses the same routines as the transient analysis 
 with minor changes, the program is put in "Mode l" during the dc calculations. 
 While in this mode, the subroutines will replace all capacitances by open 
 
 -58- 
 
"SQ- 
 
 >9 
 
 circuits and all inducatances by shorts, which leads to a relaxed condition 
 when the dc analysis is completed,. Other initial conditions could be easily 
 specified by making minor changes in the program, but these have not yet been 
 programmed „ Initially all sources have zero values; then the program "turns 
 them on" by successively increasing the values of the sources until they have 
 their desired values, meanwhile continually solving the circuit to prevent large 
 discontinuities. This process does not seem to take excessive time; yet it 
 allows bistable circuits to be set as desired. This part of the program in 
 itself provides a means for the dc analysis of nonlinear circuits. 
 
 When the dc initialization is completed, Mode 2 is set and the 
 transient analysis is started. The control of the printing, time incrementing 
 and impedance calculation is shown in the right half of Fig. A-l, QV is the 
 name of the subroutine which distributes the charge and calculates the node 
 voltages, IMP calculates the impedance factors and ALTERT changes the time 
 increment and checks if the present set of values is to be printed. The sim- 
 plified flow charts of IMP and ALTERT are found in Figs, A-k and A-5„ The terms 
 found in Fig, A-5 are defined by the following^ 
 
 DT --Present time increment 
 
 TMAX --Time at which analysis is stopped 
 
 DTM --Minimum time increment allowed 
 
 DVM --Maximum node voltage increase during cycle 
 
 DW --Maximum node voltage increase allowed 
 
 LVL --If DVM is less than this the time increment can be doubled 
 
 It is to be emphasized again that the flow charts do not represent the subrou- 
 tines exactly; their purpose is merely to present some of the major computational 
 steps in the subroutines. 
 
-60- 
 
 Figures A-2 and A- 3 contain the charge distribution and node voltage 
 calculation loops which will require the most computer time during the analysis. 
 The flow charts in this case describe the steps discussed in section 2 fairly 
 completely, but again minor points are not included. As soon as QV is entered, 
 it calls a subroutine ADJERR (Fig, A-6) which changes the adaptive factor v 
 discussed in section 5„2, At present a slightly larger than linear extrapola- 
 tion formula is used to minimize the error in future numerical cycles. This 
 is done only when there is a definite increasing trend in the error over the 
 last three cycles; otherwise an average is used and if the error is oscillating 
 between positive and negative values no correction is made at all. 
 
 At present the decrease of the time increment is based on how large 
 the node voltage changes are rather than how large the circuit error is. This 
 is a result of earlier indecision about which node voltage calculation criterion 
 to use. Experience with the program has indicated the time interval change 
 should rely more on the circuit error, 
 
 9 , 2 Sample Resul ts 
 
 To demonstrate some of the properties of the IQ method, four circuits 
 and their response will be presented here. The first, that of Fig, A-7, illus- 
 trates the effect of v on the circuit error. The maximum allowed error was 
 
 a 
 
 set arbitrarily high so that the method would not decrease the size of the 
 increments. When v = 0, we see that the error in the initial two cycles after 
 the application of the 10~volt step is very small but increases slowly^ as 
 discussed in section 5,2, The vast improvement by correctly changing v^ is 
 noted, confirming the discussion of that section. 
 
 In the second example a tunnel diode monostable (Fig, A-8) is analyzed, 
 The tunnel diode curve used is composed of three cubics for the tunneling, 
 valley and diffusion regions of the diode which were found by a least-square 
 
■61- 
 
 method of curve fittings The same data was processed by the IQ method and by 
 the second-order Runge-Kutta method, both of which gave the same response, 
 agreeing within three millivolts at any given time. The circuit error in the 
 IQ method was always less than one millivolt during the fast transients and 
 even less at other times o 
 
 In the third example a tunnel diode pulse amplifier and waveform 
 reshaper is considered^ The tunnel diode on the input side had a peak current 
 of two milliamperes and that on the output a f ive-milliampere peak. The two- 
 milliampere tunnel diode triggers the five-m.illiamper e diode which then supplies 
 the output power „ The effect of the coupling resistance R is observed in the 
 two response diagrams „ Since the average circuit error during the response was 
 less than O.5 millivolt, the waveforms can be expected to be fairly accurate „ 
 
 As a final example, the reflections from an unbiased tunnel diode 
 terminating a transmission line (Fig, A-IO) are analyzed by the IQ method. The 
 input pulse and the reflection appear in the v versus t plot. The initial 
 negatively-reflected voltage occurs because of the low tunnel diode resistance 
 in the tunneling region. After the peak is passed the resistance becomes nega- 
 tive, then very high, resulting in a positive voltage reflection coefficient. 
 The final spike results from the capacitance discharging into the line. Since 
 the sending end is correctly matched, it produces no reflections, 
 
 9-3 Efficiency 
 
 The program that is presently operative has been written in a manner 
 so as to allow easy changes; efficiency was not a primary objective and con- 
 sequently not much information can be gained by comparing the speed of computa- 
 tion of the present program with that of other methods, A more significant 
 approach, an examination of the relative time spent in different parts of the 
 program, indicates that the calculation of the branch currents and voltages in 
 
-62- 
 
 QV by far uses the most time. This then leads to the conclusion that the time 
 required to analyze a circuit will "be proportional to the number of branches 
 it haso To get an order of magnitude estimate of the speed of the IQ method 
 on a computer we need therefore only consider the operations required in the 
 branch loops of the program. 
 
 For ILLIAC 11^ for instance, it was found that about 20 additions^ 
 30 multiplications and 50 access times will be required per branch calculation. 
 Based on times of 3-5 [J^sec, 6 (asec^, and 1,8 lasec respectively for these three^ 
 this comes to a total of 350 |-isec per branch. Since branches with nonlinear 
 elements will require more time^ we can roughJ-y say that the average time per 
 branch calculation is 500 lasec. Therefore a 20-branch circuit requiring 1000 
 steps for the solution should take about 10 seconds. 
 
 Two possible methods of increasing the speed of the IQ method should 
 be mentioned. Firsts in most circuits there are many branches which are con- 
 nected to ground;, the reference node for the IQ method. Since these often 
 contain only single elements such as resistances or capacitance S;, more simplified 
 versions of the correction formula^ equation 2, l6^ could be used for these special 
 cases. Since sometimes as many as half of the branches will be connected to the 
 reference node^ this represents a significant increase in speed. 
 
 Second, for large circuits it is possible that one end is undergoing 
 transients but the other is still in a steady state. Therefore the two parts 
 of the circuit can be analyzed with different time increments. Although this 
 will lead to problems in matching the time increments the result should be 
 especially rewarding during the analysis of large transistor switching circuits. 
 
10 „ CONCLUSIONS 
 
 At present not enough experience has been had with the IQ method to 
 draw very definite conclusions about it. But;, to summarize, we can restate 
 some of its advantages and disadvantages „ Among the advantages are the 
 following: 
 
 1) An absolute error is available which the IQ method uses 
 to make corrections „ 
 
 2) An error correction mechanism is provided by the method 
 for nonlinear elements, which eliminates the need for 
 repetitive relaxation formulas o 
 
 3) The method is particularly suited for analyzing diode 
 
 and transistor circuits since the provided error correction 
 mechanism eliminates the need to integrate the transistor 
 and diode equations separately and match the solutions 
 with the rest of the circuit periodlcallyo A second 
 merit Is that no further assumptions, aside from those 
 used to get the Ebers and Moll equations, need to be made 
 because of the convenient manner in which the relevant 
 transistor quantities can be calculated. 
 h) The efficiency seems reasonable and the time required to 
 analyze a circuit is linearly dependent on the number of 
 branches in the circuit „ 
 
 5) Programming is not too complicated since no matrices must 
 be prepared and the initial conditions and transients can 
 be calculated by the same routine. 
 
 6) The IQ method is very much directed toward giving the 
 user insight into his circuit. The method deals with 
 
 -63- 
 
-6k- 
 
 quantities he can visualize and allows him to control 
 more directly how the program handles his circuit. 
 
 The chief disadvantage is that the method will not handle all circuits 
 with equal ease since it was developed for a circuit with a restricted topology. 
 This will require some work from the user before he has his circuit in a form 
 in which the method will handle it most efficiently. How serious these limita- 
 tions can be is not yet know. 
 
 It is felt that at present the IQ method is still in too much of an 
 untested state to allow a detailed comparison with other numerical methods for 
 circuit analysis. But, considering the above advantages of the IQ, method^ it 
 is also felt that it will certainly prove to be a useful tool and therefore 
 deserves further analysis and refinement. 
 
 As a final note, it should be pointed out that the IQ method developed 
 in this report is a consequence of the chosen topology and set of differences. 
 It is possible that other incremental charge methods with better characteristics 
 can be formulated by applying the same basic principles and deriving parallel 
 formulas , 
 
BIBLIOGRAPHY 
 
 [1] Branin, F. H„, Jr., "D-C and Transient Analysis of Networks Using a 
 
 Digital Computer/' IRE Int. Conv. Rec , Vol„ 10, pt. 2, pp, 236-256, I962, 
 
 [2] Dewitt, D. and A. L„ Rossoff, Transistor Electronics „ McGraw-Hill Book 
 COo, Inco, New York, New York, Chaps „ 2, 3, 1957„ 
 
 [3] Domenico, Ro J„, "Simulation of Transistor Switching Circuits on the 
 IBM 70^/' Trans, PGEC , EC-3, pp„ 2^2-2^7, Dec, 1957 = 
 
 [k] Ebers, J„ J„ and J. L„ Moll, "Large Signal Behavior of Junction 
 Transistors," ProCo IRE , vol, k2, pp. I76I-I772, Dec, 195^, 
 
 [5] Fosdick, L„ D„, "A Program for Calculating Node Voltages and Branch 
 
 Currents in Circuits with Nonlinear Elements," Univ, of 111,, Dig, Comp, 
 Lab, File No, 287, 1959, 
 
 [6] Hamming,, R, W,, Numerical Methods for Scientists and Engineers , McGraw-Hill 
 Book Co,, Inc, New York, New York, Chap, I5, I962, 
 
 [7] Leichner, G, H,, "General Tolerance Analysis Program (1206)," Univ, of 
 111,, Dig, Comp, Lab, File No, 255, 1958, 
 
 [8] Poppelbaum, W, J, and N, E, Wiseman, "Circuit Design for the New Illinois 
 Computer," Univ, of 111,, Dig, Comp, Lab, Report No, 90, I959, 
 
 -65- 
 
APPENDIX 
 
 -GG- 
 
■67- 
 
 INITIALIZE 
 
 TABLES 
 SET MODE 1 
 
 SET DC 
 TURN ON 
 LOOP ( i ) 
 
 INCR. DC 
 SOURCES 
 STOP AT 
 FULL VALUE 
 
 SET MODE 2 
 PREPARE FOR 
 TIME INCR. 
 ETC. 
 
 SET TIME 
 LOOP { i) 
 
 DECISION COMPLEX 
 REPEAT CONTINUE PRINT 
 
 IMP 
 
 I 
 
 SET 
 
 CONTINUE 
 
 INDICATOR 
 
 \ 
 
 END 
 
 Figure A-lo Flow Chart for Main Routine, 
 
ADJERR 
 
 SET NODE 
 
 LOOP 
 
 (k) 
 
 SET BRANCH 
 
 LOOP 
 
 (J) 
 
 Qji = Qj(i-i)+^Qj ♦ 
 
 Q. = Cv- 
 
 J J J 
 
 CALC. Vcj.Vrj.Vlj. 
 AQj TO LOWER NODE 
 
 CALC. A^Qj (2.16) 
 
 CALC. A_Q. (2.17) 
 K J 
 
 AQ^ = AQ^-A^Q. 
 
 YES 
 
 YES 
 
 YES 
 
 SET BRANCH 
 
 LOOP 
 
 (J) 
 
 CALCULATE 
 A^3 Q. ( 2.28) 
 A 2 Q. (2.20) 
 AQj (2.20) 
 
 Figure A-2 , Flow Chart for the Charge Distribution Cycle. 
 
■69- 
 
 '^VK = Vki-Vk(i-i) 
 
 Vk=0 
 SET BRANCH 
 LOOP 
 
 YES 
 
 
 Vk=Vk+T- 
 
 YES 
 
 V|_ = 
 ^C = Vk(i-l)-Vj"{i-i) 
 
 ACCEPT 
 CALCULATIONS 
 i.e. UPDATE 
 TABLES 
 
 Figure A-3. Flow Chart for Node Voltage Calculation. 
 
■70- 
 
 
 
 YES 
 
 NO 
 END 
 
 SET NODE 
 LOOP (k) 
 
 i 
 
 Ykm=0 
 
 SET BRANCH 
 
 LOOP (j) 
 
 ADD LOWER 
 NODE C^, R^, Lg. 
 TO Cj, Rj.Lj 
 
 YES 
 
 CALCULATE 
 
 Y = Y + — 
 km km 2 
 
 jm 
 
 CALCULATE 
 CEk'f^Ek.LEk(3l2) 
 
 Zjm=Rj+RE 
 
 Figure A-4. Flow Chart for Impedance Factor Calculations 
 
-71- 
 
 SET 
 "CONTINUE" 
 INDICATOR 
 
 DT=2DT 
 
 SET "END" 
 INDICATOR 
 
 DT= DT/2 
 
 SET "REPEAT" 
 INDICATOR 
 
 Figure A-5. Flow Chart for At Alteration Subroutine 
 
-72- 
 
 SET 
 BRANCH 
 LOOP(J) 
 
 YES 
 
 INITIALIZE 
 
 - 7 AQ 
 c A t 
 
 LINEAR^\^YES 
 
 NA(J) = NA(J )-!- I 
 
 YES 
 
 
 
 v«c= v^-2-5vi+ Vi_2 
 
 NO^^LL SAME ^^^YES 
 SIGN ? 
 
 Vj + Vi_i+Vi_2 
 
 NA(J) = I 
 
 Figure A-6o Flow Chart for Adjusting v to the Error „ 
 
-73- 
 
 > 
 < 
 
 o 
 
 I 
 
 3 
 
 o 
 
 c! 
 o 
 o 
 
 o-\ 
 
 
 o— ^^ 
 
 lO 
 
 II 
 
 > 
 
 cr 
 
 UJ 
 
 (M <-• 
 1 1 
 
 -'T 
 
 — rj- — 
 
 > 
 
 > 
 tr 
 
 
 LU 
 
 
 CC 
 
 
 UJ 
 
 Q.X 
 
 LlI 
 
 1- 
 
 
 (£ 
 
 > 
 
 < 
 
 o 
 
 1- 
 
 •-* 
 
 V) 
 
 
 1 
 
 UJ 
 
 
 
 
 > ■" 
 
 
 
 
 UJ 
 
 ^iM 
 
 
 > 
 
 
 , 
 
 
 < 
 
 Q _ 
 
 
 
 in 
 
 
 
 
 CJ 
 
 < 
 
 ■"" 
 
 ) 
 
 1 
 
 3 
 
 _l — 
 
 i 
 
 
 > 
 
 < — 
 
 > 
 UJ 
 
 
 
 
 II 
 
 
 
 1- 
 
 > 
 
 
 - 
 
 
 
 
 
 
 
 
 
 
 
 - 
 
 
 
 
 1 
 
 ^B 
 
 
 
 1 ' 
 
 1 ' 
 
 MM' 
 
 ro 
 
 Q 
 
 Z 
 
 o 
 (J 
 
 UJ 
 
 CO 
 
 > 
 
 < 
 
 "3- 
 O 
 O 
 
 cvj 
 O 
 O 
 
 C\J 
 O 
 O 
 
 O 
 O 
 
 < O 
 
 o 
 
 o 
 
 u 
 o 
 u 
 u 
 
 -p 
 
 •H 
 
 ;^ 
 o 
 
 •H 
 O 
 
 O 
 
 S 
 
 o 
 
 +^ 
 o 
 
 (D 
 
 05 
 
 0) 
 w 
 a 
 o 
 ft 
 
 CQ 
 0) 
 
 « 
 o 
 
 K 
 
 t 
 < 
 
 •H 
 
 en 
 
 CO 
 
 OT 
 
 o 
 > 
 
 o 
 
 > 
 
 o 
 
 > 
 
-7h- 
 
 + 150 MV. 
 
 a) CIRCUIT 
 
 - 500 
 
 -400 
 
 - 300 
 
 - 200 
 
 - 100 
 
 NANOSECONDS 
 
 b) INPUT AND RESPONSE 
 
 MA. 
 
 Figure A-8o Tunnel Diode Monostable Response by the IQ, Method, 
 and by the Second-Order Runge-Kutta Method, 
 
■75- 
 
 3 
 O 
 
 •—I 
 
 o 
 
 > 
 
 
 
 ^ 
 
 
 
 CL) 
 
 
 
 P4 
 
 
 
 cd 
 
 •'</) 
 
 
 CO 
 
 Q 
 
 
 <u 
 
 Z 
 
 
 K 
 
 o 
 
 
 
 o 
 
 
 0) 
 
 Ul 
 
 
 > 
 
 </J 
 
 
 cd 
 
 o 
 
 
 :5 
 
 z 
 
 
 
 < 
 
 
 '(d 
 
 z 
 
 4 
 
 
 •H 
 
 •H 
 rH 
 
 
 UJ 
 
 H 
 
 
 to 
 
 z 
 
 £ 
 
 
 o 
 
 (U 
 
 
 a. 
 
 TJ 
 
 
 </) 
 
 o 
 
 
 UJ 
 
 •H 
 O 
 
 
 a: 
 
 (L) 
 
 
 
 
 
 JD 
 
 EH 
 
 c6 
 
 O 
 0) 
 
 
 
 o 
 
 ^ 
 
 
 CO 
 
 o 
 
 
 (U 
 
 o 
 
 
 K 
 
 UJ 
 
 
 
 CO 
 
 
 
 o 
 
 z 
 
 
 jl 
 
 < 
 
 z 
 
 
 0) 
 
 •H 
 
■76- 
 
 I CURRENT 
 I PULSE 
 
 U 
 
 50n<. 
 
 Zq = 50 /V 
 
 T 
 
 DELAY = 1 NSEC. 
 
 a) CIRCUIT 
 
 5^/xf 
 
 LOAD LINE 
 AT T=2 30 
 
 T 1 1 ■ I I 
 
 100 200 300 400 500 M V. 
 
 b) TUNNEL DIODE CHARACTERISTICS 
 
 200 - 
 
 100 - 
 
 300 
 
 200- 
 
 iOO- 
 
 T 1 r— I 1 r 
 
 I 2 
 
 C) SENDING AND RECEIVING END VOLTAGE WAVEFORMS 
 
 Figure A-IG, Reflections from an Unbiased Tunnel Diode „