PART I SCALING OF LINEAR DIFFERENTIAL EQUATIONS FOR ELECTRONIC DIFFERENTIAL ANALYZERS PART II SCALING OF NONLINEAR DIFFERENTIAL EQUATIONS FOR THE ELECTRONIC DIFFERENTIAL ANALYZER Robert M. Howe The University of Michigan

TABLE OF CONTENTS PART I Page 1. Introduction..............1............. 1 2. Example of a Second-Order System......................... 2.1 Equation Describing the Second-Order System..... 2 2.2 Computer Circuit Without Regard for Scaling..... 2 2.3 Computer Circuit With Proper Scaling........... 6 2.4 Change of Time Scale........................... 9 2.5 Use of Dimensionless Time..................... 10 35.. Example of a Two-Degree-of Freedom System............... 12 3.1 Equations of Motion...............15 3.2 Preliminary Examination of the Problem.......... 15 3.3 Circuit for Solving the Ride System............. 17 3.4 Computer Circuit Including Nonlinearities....... 17 4. Scaling of Eigenvalue Problems......................... 19 4.1 Equation for Heat Flow................. 20 4.2 The Eigenvalue Problem.......................... 4.3 Computer Solution of the Eigenvalue Problem..... 24 5. Scaling of Difference Equations...................... 26 5.1 Basic Difference Equations for Heat Flow........ 27 5.2 Electronic Differential Analyzer Circuit........ 29 6. Summary.29 6. Summary.............~..................................... 29 APPENDIX A.1 Summary Amplifiers.............................. 52 A.2 Integrating Amplifiers.......................... 55 PART II 1. Introduction............................................... 2. Example of a Control System, Including Nonlinearities..... 57 2.1 Circuit Using Direct Representation of the Problem Variables............................... 59

TABLE OF CONTENTS (CONT'D) Page 35. The van der Pol Equation................................ 43 4. Solution of the Satellite and Ballistic-Missile Trajectory Problem................................... 47 5. Solution of Perturbation Equations for the Trajectory Problem..................................9................ 49 5.1 Determinination of Initial Conditions for the Example Problem................................. 55 6. Summary of Nonlinear Scaling Techniques.................. 57 APPENDIX A DERIVATION OF THE EQUATIONS OF MOTION FOR A POINT MASS MOVING IN A CENTRAL-FORCE GRAVITATIONAL FIELD.... 60 APPENDIX B TRIGONOMETRIC RESOLUTION IN ANALOG COMPUTERS BY MEANS OF MULTIPLIER ELEMENTS......................... 63

PART I SCALING OF LINEAR DIFFERENTIAL EQUATIONS FOR ELECTRONIC DIFFERENTIAL ANALYZERS 1. Introduction One of the most important problems in successful utilization of electronic differential analyzers involves the scaling of the computer, i.e., the proper choice of dependent and independent variable changes from original problem to computer problem. This scaling procedure is necessary for both analog and digital computers, but it is particularly important that it be done reasonably well in an analog computer because of the limited range of variables over which accurate computing is possible. These notes contain a discussion of methods for scaling mathematical problems for the electronic differential analyzer (analog computer). This is accomplished by considering several typical example problems among which are the common mass-spring-damper system, a two-degree-of-freedom system representing an automobile suspension system with both a linear and a nonlinear shock absorber, and finally the problem of heat flow in a slab, solved by two different methods. The scaling problem generally involves the consideration of two basic categories: 1. Scaling of the dependent variables (e.g., the problem output variables). 2. Scaling of t.he independent variables (e.g., the time scale in a dynamic problem). This second category must be considered in both linear and nonlinear problems; the approach for dynamic problems depends on whether or not it is necessary to solve the problem in real time. The first category must also be considered in both linear and nonlinear problems, but the degree of consideration is much different for linear as opposed to nonlinear problems. In linear problems it is the scaling of dependent variables relative to each other that is important while in nonlinear problems the absolute scaling of each variable is important. 2. Example of a Second-Order System As a first example, consider the mass-spring-damper system shown in Figure 2.1. Let the mass be constrained to move in the vertical direction. Denote the mass by m, the spring constant by k, and the viscous damping constant by c, and let the vertical displacement of the mass from its equilibrium position be denoted by y. The external force is f(t), where t is the time variable. -1

-2m, f(t) Figure 2.1. Mass-Spring-Damper System 2.1 Equation Describing the Second-Order System Summing vertical forces acting on the mass gives y+ cy + ky f= (t) (2.1.1) To completely specify the problem it is necessary to denote the initial conditions. These are the values of the displacement y and the velocity y at the beginning of the problem. Let y(o) - YO, y(o) - yo (2.1.2) We will now consider the possible electronic differential analyzer circuits for solving this problem. 2.2 Computer Circuit Without Regard for Scaling It will be assumed that the reader is familiar with the basic element of the electronic differential analyzer, namely, the operational amplifier. A brief description of the basic characteristics of operational amplifiers is given in the Appendix. To obtain the analyzer circuit for solving (2.1.1) we first solve for the highest-derivative term.

-3Thus from (2.1.1) my = -c -ky + f(t) (2.2.1) Starting with a voltage assumed equal to my, we synthesize the righthand side of Equation (2.2.1) as shown in Figure 2.2.1. The voltage my is first multiplied by 1/m by sending it through an attenuating potentiometer set at 1/m. The resulting signal y is integrated in amplifier 1 to obtain -y, which is in turn integrated with amplifier 2 to obtain y. These voltages, multiplied by the appropriate constants through potentiometers 2 and 3, are summed with f(t) at the output of amplifier 4. This voltage represents the right side of Equation (2.2.1) which, according to the equation, must equal the left side my. This is assured by completing the connection shown by the dashed line in Figure 2.2.1, and the analyzer circuit is now complete. Following the release of initial-condition voltages on integrators 1 and 2, the output voltages -y and y will represent the dynamic solution to the problem, which will, of course, depend upon the forcingfunction voltage f(t).. |II I I I/m. I- I k I -cy-ky+f(t)-k l ~ ~ tyf(t) f(t) o — / V' All Resistor Values are Megohms All Capacitor Values are Microfarads Figure 2.2.1. Synthesis of the Differential Equation by Means of Operational Amplifiers.

-4One advantage of the circuit in Figure 2.2.1 is that there is a one-toone correspondence between potentiometer settings on the computer and parameters in the physical problem; thus the setting of potentiometer 1 represents 1/m, while potentiometers 2 and 3 represent c and k respectively. But for almost any combination of m, c, and k normally encountered, the circuit of Figure 2.2.1 is almost certain to cause scaling difficulty. Consider the following example: m = 600 slugs, c = 2500 lb.sec/ft., k = 50,000 lb./ft. According to the circuit of Figure 2.2.1, the settings on potentiometers 1, 2 and 3 would be 1/600, 2500, and 50,000. Since one cannot obtain gains of greater than unity through a potentiometer, amplifiers 3 and 4 would need to have larger gains than the unity gains shown. For example, we might obtain a gain of 10,000 with amplifier 3 by means of a 0.01 megohm input resistor and a 100 megohm feedback resistor. Potentiometer 2 would then be set at c/10,000 or 0.25. Similarly, amplifier 4 might be given a gain of 100,000 by using a 0.01 megohm input resistor and a 1000 megohm feedback resistor. Potentiometer 3 would then be set at k/100,000 or 0.5. The resulting circuit is shown in Figure 2.2.2. I 1 1000 1/60 |. 0 oo2 600"'y ^ ""y Y -2500y -50000y+f(t 1000 100 2f(t)t Figure 2.2.2. Poorly Scaled Circuit for an Underdamped Second-Order System.

-5This circuit is an example of very bad scaling for the following reasons: 1. Summing amplifiers should never be used with gains very much larger than unity except in some unique special circuits which require it (there is no such unique requirement here). This is because amplifier noise and zero drift at the output is essentially magnified by the closed-loop gain of the amplifier, which here is 10,000 for amplifier 3 and 100,000 for amplifier 4. Also at these large gains the wide-open gain of the amplifiers (i.e., the gain before feedback) is not sufficient to insure an accurate closed loop gain. For example, if amplifier 4 has an open-loop gain of 500,000 (this is the gain from summing junction to output), then the gain with a 1000- meg feedback and a 0.01 meg input resistor will be 100,000/(1 + 100,000/500,000) or 83,300 instead of 100,000. Mathematical derivation of this relationship is given in the Appendix. Although there is no example here, it can also be stated that summing amplifiers should never be used at gains very much smaller than unity if they are part of a main computing loop. This is because no matter how small the amplifier gain, there is a certain minimum noise and drift level at the amplifier output (see the Appendix). Thus decreasing an amplifier gain to well below unity at one point in the circuit will not decrease the noise or drift appreciably at that amplifier output, but the corresponding increase in amplifier gain to many times above unity at some other point in the circuit will increase the noise or drift at that amplifier output by the factor of the amplifier gain. Thus it is best to keep all summing amplifier gains near unity when they are in a main computing loop. 2. When two or more amplifiers are used as integrators, they should never be used at time constants which are extremely different from each other if they are part of a main computing loop. The time constant of an integrator with an input resistor R and a feedback capacitor C is equal to RC seconds. Thus an amplifier with a 0.1 microfarad feedback capacitor and a 0.1 megohm input resistor should be considered as having a time constant of 0.01 seconds. In the circuit of Figure 2.2.2 integrator 1 has a time constant of 1/600 second if potentiometer 1 is considered part of the integrator, while integrator 2 has a time constant of 1 second. This is bad because in general this will mean very unequal voltage levels at the outputs of the integrators. 3. Potentiometers cannot be set accurately at values which are very small compared with unity. Thus it would be difficult to set potentiometer 1 at 1/600 to closer than 10 or 20% accuracy.

-64. Amplifier output voltages are limited to maximum positive and negative values. For most computers this limit is the order of + 100 volts. Thus the output of amplifier 4 should not exceed + 100 volts which means that the.output y of amplifier 2 must not exceed + 1/500 volt. Hence the maximum amplitude of solution we could hope to obtain with the circuit of Figure 2.2.2 is + 1/500 volt for y. This would not only be inaccurate due to noise and drift effects, but it would also be difficult to record. 2.3 Computer Circuit With Proper Scaling Our statements of the previous section said that all summers in a main computing loop should have a gain somewhere near unity, while all integrators should have about the same time constant. In the circuit of Figure 2.2.2, it should then be apparent that the output of amplifier 4 must be the order of [-cy -ky + f(t)] x 10-5 in order to have a gain of unity and a potentiometer setting of 0.5 for potentiometer 3. Thus the input to potentiometer 1 becomes my x 10-5. With potentiometer 1 set at 500/m or 5/6 the input to integrator 1 becomes y/200. This is shown in Figure 2.3.1. 500/m o0. |I s 0.05 rlS m10 K 1 c/5000. 10'5 f(t) ~ t)1 Figure 2.3.1. Correctly Scaled Circuit for an Underdamped Second Order System.

-7We should now choose the time scales of the integrators 1 and 2 to be roughly equal and. such that we pick up a gain of 200 in going from y/200 to y. In Figure 2.53.1 this has been done by using a time constant of 0.1 seconds for integrator 1 and 0.05 seconds for integrator 2. The time constants shown could equally well have been obtained by using 0.1 mfd feedback capacitors instead of 1 mfd, in which case the input resistors to integrators 1 and 2 would have been 1 and 0.5 megohms respectively. The output of integrator 1 is now -y/20, and potentiometer 2 is set at c/5000 or 0.5 in order to produce an output of [cy - f(t)] xlO 5 at amplifier 3. Note that the input to amplifier 3 is now f(t) x 10-5. The circuit of Figure 2.35.1 has achieved the objective of near unity gain for the summers and time constants about equally distributed among the integrators. All potentiometer settings are reasonably close to unity. Since the problem considered here is linear, the choice of relationship between problem units and computer units for the dependent variable, y, is rather arbitrary. If the maximum expected mass displacement y is one foot, then we might let one foot = 100 volts, or full scale output at amplifier 2. A 100 volt output of amplifier 2 would then represent a displacement of one foot; a 100 volt output of amplifier 1 would represent a velocity y of -20 ft./sec. A 100 volt output of amplifiers 3 or 4 would represent a force of 105 lbs. A step function input force f(t) of 104 lbs. would consist of a step voltage of 10 volts applied to the input of amplifier 3. Several times we have referred to the "main computing loop." By this we mean the chain of amplifiers which has the primary influence on the solution to the problem. In more complicated problems the main loop or loops are often difficult to recognize, particularly if the general nature of the solution is not known. In the simple mass springdamper problem of Figure 2.3.1 the main loop consists of amplifiers 1, 2, and 4 for an underdamped system (~ = 0.228 in the example considered, where - is the fraction of critical damping). If, on the other hand, c = 50,000 lb. sec./ft. instead of 2500 lb. sec./ft., as given earlier, the system is overdamped. ( = 4.57) and amplifiers 1, 3 and 4 constitute the main computing loop. We might then denote the output of amplifier 1 as -y, in which case, potentiometer 2 would be set at c x 10-5 as 0.5; this would maintain the output of amplifier 3 at [cy - f(t)] x 10-5 as before, and the output of amplifier 4 at [-cy -ky + f(t)] x 10-5 as before. (See Figure 2.3.2).'The input to potentiometer 1 is my x 10-5. If we set potentiometer 1 at 500/m as before, the time constant of integrator 1 must be 1/200 in order to provide an output of y. This can be achieved with a 0.05 megohm input resistor and a 0.1 mfd feedback capacitor as shown in the figure.

-80.1 1 I 10 y o / M I >\ [- cy-ky+f(t)]ly. 3 I 10(5 f(t) VA Figure 2.3.2. Correctly Scaled Circuit for an Overdamped Second-Order System. In the steady state following a constant input force f(t), the velocity y will be negligible, as will the acceleration y. Hence the input voltages to amplifier 4, one coming from y, one coming from f(t), must be equal and opposite. Thus we expect y values the same order as 10-5 f(t), so that we shall let the output of amplifier 2 be y as before, rather than some constant factor times y. Reference to Figure 2.53.2 shows that this is a case where the integrator time constants are nowhere nearly equal for a properly scaled circuit. The important point is that amplifier 2 is not part of the main computing loop, at least as far as the initial transient operation is concerned. The main loop consists of amplifiers 1, 3 and 4, as explained earlier. Had we used the circuit shown in Figure 2.3.1 for the case considered here namely, C = 50,000 lb. sec./ft., amplifier 3 would have needed a gain of at least 10, and we would have been amplifying a small signal -y/20 from amplifier 2 by a large gain in amplifier 3 in order to produce a nominal output cy x 10-5. The superiority of the circuit in Figure 2.3.2 is evident. Here all amplifier voltage levels will be the same order of magnitude.

-92.4 Change of Time Scale We saw in Figure 2.3.1 for the case where m = 600, c = 2500, and k = 50,000 that integrators 1 and 2 have time constants of 1/10 and 1/20 second respectively. The resulting dynamic solutions will exhibit damped oscillatory transients of about 1-1/2 cps because on, the undamped natural frequency, equals fk/m = 9.12 radians/second. Although such solutions will not tax the bandwidth capabilities of the operational amplifiers, they may be too fast for some recorders, e.g., X - Y type recorders. If the problem had one or more nonlinearities for which servo-driven multiplying potentiometers or function generators were used, oscillations at 1.5 cps might also be too fast. Then, too, some dynamic systems may have oscillatory transients which are so high in frequency that they are too fast for the operational amplifiers themselves. Or they may be so low in frequency that low-frequency limitations such as amplifier drift or the patience of the computer operator become a factor to consider. In any case it is frequently desirable to either speed up or slow down the computer time scale relative to the original problem time scale. About the only situation where this cannot be done is when a "real time" simulation is necessary because of the tie in of physical hardware with the computer. Probably the simplest argument to explain the method of time scale change for the circuit in Figure 2.5.1 begins with the following realization. In discussing the circuit shown it was assumed that 1 second of computer time equaled one second of problem time. Thus an integrator with a 0.1 megohm input resistor and a 1 microfarad feedback capacitor has a time constant of 0.1 second in problem time. But if we assume that 10 seconds of computer time equals one second of problem time, then an integrator with a 1 megohm input resistor and a 1 microfarad feedback capacitor still has a time constant of 0.1 second in problem time, i.e., the computer will run ten times as slow as the actual problem. Thus by increasing the input resistors of integrators 1 and 2 in Figure 2.53.1 from 0.1 and 0.05 to 1 and 0.5 megohms respectively, we can slow the problem down by a factor of 10. This might then allow us to use a slower responding recorder. On the recorder trace 10 seconds along the abscissa would equal one second for the actual problem. Thus the computer can be slowed down or speeded up relative to "real time" simply by changing the input resistor or feedback capacitor of all integrators by the same factor. For the case discussed in the previous paragraph we slowed the computer down by a factor of 10. It is equivalent to introducing a new time variable T given by T = 10 t (2.4.1)

-10from which d = i dcT = 10 L (2.4.2) dt dr dt dT and =100 (2.4.3) dt2 dT2 Written in terms of the new time variable T, Equation (2.1.1) becomes 100 + + lOc ky f(t) (2.4.4) dT2 dT or -2 100 m d = - lc d - ky + f(t) (2.4.5) dT- dr The computer circuit for solving Equation (2.4.5) is shown in Figure 2.4.1. It is precisely the circuit suggested in the previous paragraph in order to slow the problem down by a factor of 10, but here it was arrived at formally by rewriting the equation with a change of time variable, whereas before we simply took the circuit of Figure 2.53.1 and increased all input resistors to integrators by a factor of 10. Either method is suitable, and it becomes more a matter of personal preference as to which way one prefers to look at the problem of time scale change. 2.5 Use of Dimensionless Time Another approach to the time-scale problem, and probably the most preferable approach, in many cases, is to introduce a dimensionless time in the differential equation to be solved. Consider again Equation (2.1.1). By dividing through by k we can write m'y +r + y = 1 f(t) (2.5.1) k k k This equation can be rewritten in terms of the undamped natural frequency On and the damping ratio ~. Thus 1 — y + y+y =-1 f(t) (2.5.2) c n C>n k

-11C C I R R.. d2y | - t - [ — 2 - d d rt k- f(z-) o- - Va Figure 2.5.1 Circuit for Underdamped Second Order System Using Dimensionless Time. - 2; F C R dr" ~ f()0-^Ar Figure 2.5.2 A Three-Amplifier Circuit for an Underdamped Second-Order System Using Dimensionless Time.

-12where con =^ 3(2.5.3) m; = =c (2.5.4) 2 -fmk Let us now introduce a dimensionless time variable T given by T = ant (2.5.5) from which d = d2 2 d(2.5.6) dt dT dt2 dT 2 Equation (2.5.2) then becomes d + 2 + y= 1 f(t) (2.5.7) dT2 dT k The analyzer circuit for solving this equation is shown in Figure 2.5.1. Note that the time constant for each integrator is RC seconds. Our choice of RC will dicate the number of seconds of computer time equivalent to one unit of dimensionless time T. If R = 1 megohm and C = 1 microfarad, then 1 second of computer time equals one unit of dimensionless time T. The advantage of introducing the dimensionless time T = Ont for this problem is obvious; the scaling problem for reasonable ~ more or less solved itself, and the change of computer time scale is particularly easy to visualize. Note also that the forcing function is actually a displacement, ~ f(T) so that all voltage inputs and outputs represent displacements. Hence a single choice of relationship between volts and feet is adequate. A three amplifier circuit equivalent to the circuit of Figure 2.5.1 is shown in Figure 2.5.2. 3. Example of a Two-Degree-of Freedom System We have already examined in Section 2 the many scaling problems which arise in determining the electronic differential analyzer circuit for solving a mass-spring-damper system. A variety of solutions to these scaling problems was discussed. In this section we turn to a

-13two-degree-of-freedom system. Here we have all the scaling problems of a one-degree of freedom system with a few more thrown in. M A Y. X t x(t) Figure 5.1. Simplified Automobile Suspension System. Consider the system shown in Figure 53.1, which represents a simplified automobile suspension system. The upper mass M exhibits a vertical displacement y1, with respect to a fixed equilibrium reference, and represents the sprung iass. It is coupled to the lower, unsprung wheel of mass m with a spring constant K and viscous damping constant C (shock absorber). Displacement of the lower mass is denoted by y2, and k represents the spring constant of the tire. Finally, x(t) is the vertical displacement of the support (road). The following constants represent approximately one-half the front end of a 1956 Chevrolet. M = 22 slugs, C = 100 lbs.sec./ft., K = 1600 lbs./ft. The wheel mass and tire spring constant are given by m = 5 slugs, k = 12,000 lbs./ft.

-14 —'~~o~ ~Go ~~~0~~~~~~~0 ~~~1 L 7 a>~1, lo _ i2 T / - ~~~~~~~\ <~~- p S~~~~~ e c

-153.1 Equations of Motion The equations of motion describing the system are obtained by summing the forces acting on the upper mass and on the lower mass. The following two equations result: Myl + C (Y - Y2) + K (y1 - y2) = 0 (5.1.1) My2 - C ( - y2) K (Yl - Y2) + k (Y2 - x) = 0 (3.1.2) 3.2 Preliminary Examination of the Problem Before attempting to formulate the computer circuit, including scaling, it is well to make some preliminary calculations concerning the expected behavior of the system. This will often aid in such decisions as time scale, relative magnitude of displacements, etc. First let us calculate the natural frequency of each of the mass systems assuming that the other is held fixed. For the upper mass M 1 o = EI-K = 6 - 8.52 rad/sec. 1J M 22 For the lower mass m -2~n = - =" = 67.5 rad/sec. c- m 5 On the basis of these frequencies, we would guess that the time constants of the two integrators which will be part of the upper-mass system should be the order of 0.1 seconds (this will give a natural frequency of about 10 rad/sec), while the time constants of the two integrators for the lower system should average about 0.02 seconds (this would give a 50 rad/sec. frequency in conjunction with a unity gain summer). These integrator time-constants are, of course, based on the assumption that we wish to solve the problem in real time. Let us proceed on this basis; later we can easily slow the problem down by changing all four integrator time constants by the same factor. We might also note that although the displacements yl, Y2 and x will be the same order of magnitude, the relative displacement Y2 - x will be much smaller (this represents the compression of the tire). This means that even though Y1, Y2 and x are scaled to cover the full-scale voltage range, Y2 - x can be followed by a gain factor of 4 or 5 without danger of amplifier saturation.

-16^N 0 _ oil K 10C 0 0 <~ E-i a) S0, ~ 3 X _ \ IN >l ~ ^0 <( 0 I --------— 0A^-A-^ A^ ------- ^" 0 0' I 03 d*~MC C M * 0 0, 0 C5 0~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~.o0' d <lo n' ^^ ^Z,<,- <3 o~o~d~d~o^

-173.3 Circuit for Solving the Ride System The circuit for solving Equations (3.1.1) and (.1.2) on a one-to-one or real time scale is shown in Figure 3.3.1. It was arrived at by keeping in mind the preliminary considerations of the previous section, and by applying the basic method of approach presented in Section 2.3. Note that all potentiometers in Figure 35.53.1 have settings between about 1/3 and unity. The time constant of the different integrators is the same order of magnitude as the predictions of the previous section. Even though the output -Y2 of amplifier 4 is scaled so that 100 volts (full scale) equals the full excursion of Y2 (say 0.5 ft.), and amplifier 8 has a gain of 4, its output 4(y2 - x) is not likely to saturate. Note also that we have solved for +yl and for -Y2The opposite signs allow us to take the necessary differences yl - Y2 and Yl - Y2 without additional inverters. By taking these differences we require only a single potentiometer to set K, and a single potentiometer for C. If we had collected terms in Y2 in. Equation (3.1.2) and computed (K + k)y2, this would not have been possible. As the size of the problems considered becomes larger, it is always more difficult to perform a near-optimum job of scaling before the problem is set up on the computer. What one tries to do is to make reasonable estimates of the scaling based on prior knowledge of the problem. It will often be necessary to readjust the scaling once the problem has been set up and run. Frequently one finds that one or more amplifiers saturate, or that the output level of one or more amplifiers is too low. When this is the case, the modification of gains necessary to correct the situation is apparent. We have already pointed out that if the time scale of the circuit in Figure 3.3.1 is too fast, the problem can easily be slowed down by changing all of the integrator time constants by the same factor. For example, by increasing all feedback capacitors from 0.1 mfd to 1 mfd, 10 seconds of computer time will equal 1 second of problem time. 3.4 Computer Circuit Including Nonlinearities Up to now -we have discussed only linear problems. Actually, one of the chief advantages of the electronic differential analyzer (or any high-speed computer) lies in its ability to solve nonlinear differential equations. For example, in the automobile ride problem we might wish to represent the viscous damping constant C as a nonlinear function of the velocity (yl - Y2). This can be accomplished on the computer by the use of a diode function generator (other methods

-18are also available but will not be considered here). The diode function generator can be set up to provide a series of segmented straight line approximations to the function under consideration, where the break point between segments can be set anywhere over the full voltage range of ~ 100 volts, and where the slope of each segment can be set so that the output voltage of the function generator falls anywhere within ~ 100 volts. Because of limited function-generator accuracy it is best to use as much as possible of the + 100 volt range for both the abscissa and ordinate. cyl- ) (LBS) +100 -/ I' /' +1 Y1Y2 (ft/sec), /...__ ~ _ _ - - ~J - #100 Figure 3.3.2. Nonlinear Shock Absorber Characteristic To deal with a specific example, let us assume that the shook absorber saturates at 100 lbs. when the relative velocity is 1 ft./sec. In Figure 3.3.2 is shown the resulting plot of viscous damping force C (Yl - Y2) versus Y1 - Y2. If we let 0.5 ft. = 100 volts for our basic scaling in the problem, this allows a maximum excursion for yl and Y2 of + 0.5 ft. or 1 ft. total. Since the output of amplifier 7 in Figure 35.5.1 is -c(yl - y2)/4000, it will represent viscous - damping force in lbs. x 4000. Assuming that 0.5 lbs. = 100 volts, a viscous force of 2000 lbs. will produce 100 volts at the output of amplifier 7. Let us now assume that the diode function generator is inserted following amplifier 7. From Figure 3.35.2 it is evident that we want the

-19function generator to pass the viscous damping force output without attenuation until it reaches + 100 lbs., at which point the output must saturate. From the previous paragraph we see that this will correspond to a + 5 volt output, since 2000 lbs. = 100 volts from amplifier 7. Hence the diode function generator must be set up as shown in Figure 3.353. VOLTS OUTPUT +5.I ~/"/~ ~ -^ ~~~~DIODE -5 ____,,U_^ ___ C f <\ GENERATOR / +5 VOLTS 4000 YY! / |' INPUT TO AMPLIFIERS I AND 3 - -- -5 <- I Figure 3.3.3. Diode Function Generator Simulation of Nonlinear Shock Absorber. Actually, it would be better to utilize more of the full-scale output of the function generator than + 5 volts. For example, if the function generator output were fed into amplifiers 1 and 3 through 1.0 meg resistors instead of 0.25 meg resistors, the function would need to be magnified by a factor of 4, and the characteristic of Figure 3.3.4 would apply. This would use more of the full-scale output range of the function generator. 4. Scaling of Eigenvalue Problems Since the electronic differential analyzer can integrate only with respect to time, it can solve only ordinary differential equations. In order to solve partial differential equations, we must first convert them to one or more ordinary differential equations. This can be done either by separation of variables, in which case an eigenvalue problem results, or by finite-difference techniques. In this section we will consider the eigenvalue problem.

-20VOLTS OUTPUT +20 __ + 5 VOLTS INPUT -20 Figure 35.53.4. Improved. Diode Function Generator Simulation of Nonlinear Shock Absorber. 4.1 Equation for Heat Flow As a simple example of a well-known partial differential equation, let us consider the heat flow problem shown in Figure 4.1.1. Here we wish to find the temperature u in an infinite conducting slab as a function of distance x through the slab and time t. The left boundary at x = 0 is held at constant temperature uo, while the right boundary is an insulator. If C is the heat capacity per unit volume of the slab, and K is the thermal conductivity, then the equation of heat flow balance is C. - = K x (4.1.1) at x with boundary conditions u(O,7) = 0, (L,) = 0 (4.1.2) In solving partial differential equations it is almost always wise to introduce dimensionless time and space coordinates. First consider a

-21FiNgur 4.T / e TEMPERATURE aSb N...,. / INSULAT x-O x=L Figure 4.1.1 Heat Flow Through a Slab.

-22dimensionless space coordinate x given by x = - (4.1.5) L Then a = adx _ = a (4.1.4) bx bx d5 L ax In general the heat capacity C and thermal conductivity K may be functions of x. Let C(x) = Co Oc (x) (4.1.5) and K(x) = Ko X (x) (4.1.6) Where C0and K~represent the maximum values of C and K, and where 0c(x) and OK(x) are dimensionless functions which represent the variation of heat capacity and conductivity with x, and which have a maximum value of unity. In terms of (4.1.4), (4.1.5), and (4.1.6) Equation (4.1.1) becomes O OC(x) au = 1 [k(x) u] (4.1.7) Ko t L 6 x Next we introduce a dimensionless time variable t given by Kot t= (4.1.8) In terms of t Equation (4.1.7) becomes 0C(x) U = [k(x) ] (4.1.9) with boundary conditions u(O,t) = 0,.(l,t) = 0 (4.1.10) ox Let the initial temperature distribution be denoted by u(x,o) = U(x) (4.1.11) Equation (4.1.9) with boundary conditions of (4.1.10) and the initial condition of (4.1.11) defines our heat-flow problem.

-234.2 The Eigenvalue Problem To solve Equation (4.1.9) by separation of variables we assume that u(x,t) = X(x) T(t) (4.2.1) Substituting (4.2.1) into (4.1.9) we have X(x) oc(x) aT(t) = T(t) d [0k(x) X(x)] dt dx dx or 1 dT(t) = 1 d [k(x) aX(x)] T(t) dt 0,(x) X(x) dx dx (4.2.2) Since the left side of (4.2.2) is a function only of time t while the right side is a function only of distance x, to be equal for all t and x they must both equal a constant, say -p. Thus 1 dT = _- or dT + T = 0 (4.2.3) T dt dt In the same way [ik(x),] + P c(x) X = 0 (4.2.4) dx dx The solution of Equation (4.2.3) is simply the exponential decay T = Ae-t (4.2.5) where A is an arbitrary constant. The boundary conditions (4.1.10) of the original problem will be satisfied only if X(o) = (1) (l) = (4.2.6) Equation (4.2.4) can be solved with the boundary conditions of (4.2.6) only for discrete values of the constant A. These discrete values, (n' are called eignevalues, and the corresponding solutions Xn(x), are called eigenfunctions or normal modes. Once the Pn and Xn are computed, the complete solution can be written as the infinite series of functions. 00 u(x,t) = Z AnXn (x)e-nt (4.2.7) n=l where the constants An are determined by the initial temperature distribution U(x).

-244.3 Computer Solution of the Eigenvalue Problem We will now consider the computer circuit for solving Equation (4.2.4). To simplify the problem, assume that the thermal conductivity is constant throughout the slab, so that XK(x) = 1. Then Equation (4.2.4) becomes d2X + c(x) X = 0 (4.5.1) dx2 with boundary conditions X(o) = dX (1) = 0 (4.5.2) dx Since the electronic differential analyzer can integrate only with respect to time, we must let x correspond to time on the computer. Thus (4.3.1) is analogous to a mass-spring problem with a time-variable spring constant. The boundary conditions become initial and final conditions; we must find a computer solution which, starting with zero initial X, have zero X1 one unit of computer time later. We will, of course, find such computer solutions only for discrete values Pn of the parameter P. The lowest such value of P we call Pl, the next lowest 23, etc. In determining the computer circuit we must try to make the time constants of the two integrators about equal and the gain of any summers about unity. With this thought in mind the circuit of Figure 4.3.1 has been drawn. In this circuit RC seconds equals one unit of computer time and hence unity in x. Thus if we wish the problem to run for 5 seconds, corresponding to x running from 0 to 1, we might let R = 5 megohms, C = 1 mfd. Since potentiometer 1 cannot exceed unity in its setting but should not be too low in value, the additional factor o is incorporated with P, so that the potentiometer is actually set at c23. In this way we can always choose a such that potentiometer 1 is reasonable in value. The method of solution is as follows: We select a convenient time length for the problem} say 5 seconds. This sets RC for the integrators at 5 seconds. We then guess at a value for P and hence a reasonable setting for a. Starting with zero initial X and arbitrary dX/dx (initially the output of amplifier 1 should be close to 100 volts to use up full scale) we run the problem and by means of a recorder determine whether the output of amplifier 1 (i.e., dX/dx) vanishes 5 seconds later. If it doesn't (and it never will the first time), we try a new value of f and repeat the run. After a moderate number of trials it is possible

-25C c I L/5 a 1/5 a/ Em aX /dX v V a -a a2, stc dx 0c -100 VOLTS R/5 FUNCTION GENERATOR (x) Figure 4.5.1. Circuit for Determining Eigenvalues. to converge on the correct p necessary for an exact solution which satisfies the end condition. During the convergence process it may be necessary to change the setting of a in order to keep potentiometer 1 at a reasonable setting. A sketch of some typical solutions along with the ( and ce settings is shown in Figure 4.3.2. For each solution the initial condition on amplifier 1 should be adjusted as necessary to make X, the output of amplifier 2, reach a maximum of near 100 volts. This will allow accurate use of the multiplier. Note also that the input x to the function generator is obtained from an integrator with the same time constant as integrators 1 and 2 and with a -100 volt input. This generates an output starting at zero and sweeping in RC seconds to 100 volts. This represents the input necessary for the function generator to use up its full scale. Since.c(x) was designed to have a maximum value of unity, this will correspond to a maximum output voltage of 100 volts from the function generator.

-26/3- 0.12 a= 2 0? /30 o.9l a = X2! /a = 2.20 a - 0.5 X337? Figure 4.3.2. Typical Computer Determinations of Eigenvalues. When using nonlinear equipment it is particularly important to scale inputs and outputs so that they use up nearly the full available voltage range, in most cases + 100 volts. This is because the nonlinear equipment, either function generators or multipliers, is the weakest link with regard to accuracy, particularly zero drift. This is not quite such a problem when servo-type nonlinear equipment is used, however. 5. Scaling of Difference Equations It was pointed out at the beginning of Section 4 that partial differential equations must be converted to ordinary differential equations before they can be solved with the electronic differential analyzer. In Section 4 this was accomplished by separation of variables, resulting in an eigenvalue problem. In this section it will be accomplished by replacing derivatives with respect to the distance variable by finite differences.

-275.1 Basic Difference Equations for Heat Flow As in the previous section, the problem of heat flow in a slab will be used to illustrate the solution of a partial differential equation by difference techniques, and, in particular, the scaling of such a problem. Consider the problem illustrated in Figure 4.1.1. After introducing a dimensionless distance variable x = x/L and a dimensionless time variable t = (Ko/CoL2)~E we have for the equation describing the temperature U in the slab c(x) = A [() (4.19) with boundary conditions u,(o,t) = o, xU (l,t) = 0 (4.1.10) and initial condition u(x,o) = U(x) (4.1.11) Instead of considering the temperature u at all values of x, let us consider u only at certain discrete stations in x,. as shown in Figure 5..1.1. xFo Dse Ax S x=l Figure 5.1.1. Slab Divided Into. Discrete Stations.

-28Let the separation Ax between stations be a constant, and let uo denote the value of u at x = 0, u1 denote the value of u at x = Ax, u2 denote the value of u at x= 2Ax, un denote u at x = nAx. Then clearly a good approximation to the heat flux X au/8x at the n-l/2 station is Kxx in K/ 2 1 (un - Un-1) (5.1.1) Here Kn /2 is the value of the dimensionless heat-conduction parameter X at x = (n-l/2)Ax. In the same way we can approximate 6x 6x at the nth station as -[K -] = OK{[K-] - [OK-] } ax ax n Ax x n + 1/2 x n - 1/2 (5.1.2) Thus the equation for heat flow balance at the nth station becomes ocn dt = 1 [)Kn + 1/2 (un+l - un) - ~Kn - 1/2 (un - un-l)] (5.1.3) Note that since un is the temperature at a fixed value of x (e.g., x = nAx), 6un/6t becomes the total derivative dun/dt with respect to time. Equation (5.1.53) is iterated for each of the stations across the slab. If there are N stations in all, then a set of N simultaneous first-order equations results. The boundary condition at the left edge which says that u(o,t) = 0 is replaced by uo = 0 (5.1.4) and the boundary condition at the right edge which says that a (lt) = 0 is replaced by UN = uN + 1 (5.1.5) The initial condition u(x,o) = U(x) becomes an initial condition on each of the station temperatures un.

-295.2 Electronic Differential Analyzer Circuit From the discussion of the boundary conditions in the previous section it should be clear that for N active temperature stations, there are actually N + 1/2 distance increments Ax making up the slab thickness of unity, i.e., the right boundary occurs at the N + 1/2 station, and Ax = 1/(N + 1/2). For reasonable values of N, say N = 10, the coefficient 1/(Ax)2 in Equation (5.1.3) will be very large. This will make scaling of the computer circuit more difficult, so the (Ax)2 term is often incorporated into the time variable. Thus let a new time-variable T be given by T = (N + 1/2)2 t (5.1.6) Equation (5.1.3) then becomes Cn dT = + 1/2 (un+l - un) - Kn- 1/2 (un Un-1) (5.1.7) where the fc or K + 1/2 are the order of unity. Thus the scaling of Equation (5.1.7) is particularly simple. The electronic differential analyzer circuit is shown in Figure 5.2.1. Note that RC seconds equals unity in T. Often RC = 0.1 second is chosen for the integrator time constants when solving heat-flow problems. This is because for large N the dimensionless time variable T must take on large values to represent reasonable excursions in the original dimensionless time variable t [see Equation (5.1.6)]. Hence it is convenient to speed up the computer by at least a factor of 10 in relation to the variable T. To summarize, in scaling a partial differential equation for the electronic differential analyzer when using difference techniques, the original Equation (4.1.1) was simplified by introducing dimensionless distance x and dimensionless time t. After writing the difference Equation (5.1.3) we eliminated the Ax from the equation by introducing a new dimensionless time variable T. The resulting Equation (5.1.7) was then in such a form that the computer circuit of Figure 5.2.1 satisfied the basic aims of scaling, i.e., summers with unity or near unity gain, integrators with equal time scales, and potentiometer settings near unity. 6. Summary The basic points to remember in scaling mathematical problems for the analog computer are as follows:

-30(=0here) c r lI R<O I I N ~~~~~~~~~ — Fige 5.2.1 CrcfrOoa aFo R I I!. *k, -I 2' %1 C-_d U' Fiur 5.. i

-311. Set the gains of all important summing amplifiers close to unity. 2. Set the RC time constants of all important integrators to approximately the same value. 3. Arrange the circuit so that all important potentiometers are set between 0.100 and 0.999. 4. Arrange the time scale in the machine so that the bandwidth capabilities of the operational amplifiers, nonlinear computing elements and recording equipment are not exceeded. Calculation of system natural frequencies will help in choosing a time scale. Use dimensionless time wherever possible. Although all these rules must be broken from time to time, it will generally be true that when a mathematical problem has been properly scaled for the analog computer, the output voltages of all the important amplifiers in the circuit should approach full scale at some time during the solution.

APPENDIX BASIC CHARACTERISTICS OF OPERATIONAL AMPLIFIERS A.1 Summing Amplifiers The basic building block of the electronic differential analyzer is the operational amplifier. This consists of a high-gain d-c amplifier with a feedback impedance and one or more input impedances. For summation the impedances are all resistors, as shown in Figure A.l.l. Ra ia ea(; Rf if INPUT VOLTAGES ebOUTPUT e \ J? HIGH GAIN DC AMPLIFIER, GAIN = -, Figure A.i.1. Operational Amplifier as a Summer Here ea. eb, and ec are the input voltages, and eo is the output voltage. If we neglect any current into the amplifier itself (this neglects the grid current in the first stage of vacuum-tube amplification and is normally justifiable), then the sum of the input-currents equals the feedback current or ia + ib + ic = if (A.I.l) -52

-33But from Ohm's law, i ea - et 1a Ra where e' is the input voltage to the amplifier proper. Using similar expressions for ib, ic, and if, we have ea - e+ eb - e e + e'T e' - e(A.1.2) Ra Rb Rc Rf But eo = -te'l, where 4 is the gain of the amplifier (the minus sign takes care of the phase reversal). Thus et = -eO/, and Equation (A.1.2) can be solved for eo, obtaining 1+(1 Rf Rf Rf E (1 + + —a +- eb + e) Ra Rb Rc Normally the amplifier gain ji is quite large (between 30,000 and. 100,000,000) and hence > ~> 1 + Rf/Ra + Rf/Rb + Rf/Rc for a properly scaled summer. Thus the denominator in (A.1.3) is essentially unity, and the summer output is seen to be proportional to the sum of the input voltages, the constant of proportionality in each case being equal to the ratio of feedback to input resistance. Next let us consider the effect of noise in the amplifier. If the noise is very low in frequency,we refer to it as drift. In any case it is convenient to designate the noise as the equivalent voltage en which would need to exist at the amplifier input in order to produce the same effect on the output as does the noise. If this is the case, thenthe equation relating e to eo is eo = -L(e' + en) (A.1.4) We can eliminate e' between Equations (A.1.2) and (A.1.4), obtaining [R ea + Rf eb + Rfec] - [1 + Rf + Rf + Rfen Ra Rb Rc Ra Rb Rc (A1. 5) Here we have assumed p. ~ 1 + Rf/Ra + Rf/Rb + Rf/Rc

-34The second term is the output due to the noise en referred to input. Clearly the higher the gain of the summer, the larger the noise component in the output. For example, for a summer with three inputs representing a gain of 10 each, the output noise = 31 en. For en = 10-3 volts offset due to drift, the output drift will be 31 x 10-3 or 0.03 volt. For en = 5 x 10-3 volts rms of high-frequency noise, the rms noise in the output will be 0.15 volts rms. On the other hand, no matter how low the summer gain is made, the output noise will always be at least en. This may be seen by setting Ra = Rb = Rc = ~o in (A.1.5). Finally, let us consider the dynamic error of a summer resulting from the limited bandwidth of the dc amplifier. This is best illustrated by calculating the phase shift of the summer for sinusoidal inputs of frequency c. One might realize, of course, that the d-c amplifier gain p. is. actually a function of frequency, and for reasons of stability will fall off approximately as 1/w for large frequencies co. In general let (Xo) = 1(coX) eeJl(W) where.Ll is the magnitude of the amplifier transfer function and 01 is the phase shift of the amplifier transfer function. Then from (A.1.5) it is apparent that the summer output voltage due to the input voltage ea = EeiJt is given by - Rf Eet = Ra -----— (A.1.6) 1+.. (1 +-!+ L+ 1 ( Ra Rb Rc Often the amplifier phase shift 01 is the order of -90~ or -t/2 radians. If the gain R_.ff Rf Rf ~1 >> (1 + + + R) Ra b Rc which is in general still a good.assumption, then the summer output eo for the sinusoidal input ea can be written as a where the summer phase shift E is given by where the summer phase shift is given b+ + (A..9) p1 Ra Rb Rc

-35As in the case of the noise or drift effects, the summer phase shift becomes proportional to one plus the summation of the summer gains. Thus the larger the gains, the more the phase shift. However, even for zero gain the phase shift is 1/1. Thus the maximum bandwidth is obtained by distributing gains evenly through summers. A.2 Integrating Amplifiers When the feedback resistor Rf in Figure A.l.1 is replaced by capacitor C, the operational amplifier becomes an integrator. Since the capacitor has an impedance l/Cp, where p is the operator d/dt, Equation (A.1.3) becomes ea + R1 eb + _ e_ eO = — (RaC e-+-RbCp e — + RcCp (A.2.1) 1i[l+ 1 + 1 1 + [1 + RaCp RbCp RcCp] if A>> [1 + 1 + 1- + 1 ], RaCp. RbCp RcCP then 1 1 1 o = ff [lRC ea + RbC eb + R-C ec] dt (A.2.2) where 1/p has been replaced by the integral. Thus the output is proportional to the integral of the sum of the input voltages. Next let us consider the effect of noise or drift voltage en referred to the integrator input. Replacing Rf by 1/Cp in Equation (A.1.5) we have f [RC. ea + R eb + R e]'dt - [en + / ( + 1 + ) endt (A.2.5) RaC RbC RcC) enCt] Here the second term in the brackets is the output resulting from the noise en. It not only appears directly but is also integrated, as if en applied to each of the input resistors. This second error can become

-36quite serious for open ended integrations when en represents a drift voltage. For example, if en = 10-3 volts and U1 + 5 + = 3 = 30 sec-1, RaC RbC RcC 0.1 after 20 seconds the integrator output would be 20 (30) 10-3 = 0.6 volt. Up to now we have neglected the grid current ig flowing into the amplifier. This may sometimes not be possible for integrators. One can show that grid current produces an integrator output voltage given by f(ig/C)dt which, of course, is added to any other integrator outputs. Note that the effect of grid current is independent of the input resistors. For example, if ig = 10-10 amps and C = 0.1 microfarads, the output of an open-ended integrator will increase at 10-3 volts per second due to grid current. Finally, consider the bandwidth problem with integrators. We have seen that the ideal integration given in Equation (A.2.2) is achieved only if the amplifier gain 1 1 1 > [1 + RaCp + RbCp + RCp At very high frequencies Ct will fall off until this approximation no longer holds. In general the additional integrator phase shift at a high frequency X is approximately -1/p1, where il is the amplifier gain at that frequency. Generally this is not a severe limitation compared with summer phase shifts. At very low frequencies (1 + 1 + 1) becomes Ra Rb kRc Jc which can become quite large, so that again the approximation A~ [1+ no+ gr) i] is no longer valid. Also, leakage resistance of the feedback capacitor is a limiting factor at low frequencies. Thus integrators cannot be operated at too low a frequency. However, since the leakage resistance of the capacitors is generally 1012 ohms or higher, and since p may be 107 or greater for drift stabilized amplifiers, this low-frequency limitation is seldom met in practice.

PART II SCALING OF NONLINEAR DIFFERENTIAL EQUATIONS FOR THE ELECTRONIC DIFFERENTIAL ANALYZER 1. Introduction In Part I we treated a number of linear examples, including a mass-spring damper system, a two-degree-of-freedom automobile ride problem, and the equation of heat flow solved first by separation of variables and secondly by the difference method. In one of these examples, the two-degree-of-freedom ride problem, the inclusion of a simple nonlinearity in the viscous damper was considered. At the end of Part I the following basic points to remember in scaling problems for the analog computer were made: 1. Set the gains of all important summing amplifiers close to unity. 2. Set the RC time constants of all important integrators to approximately the same value. 3. Arrange the circuit so that all important potentiometers are set between 0.100 and 0.999. 4. Arrange the time scale in the machine so that the bandwidth capabilities of the operational amplifiers, nonlinear computing elements, and recording equipment are not exceeded. Calculation of system natural frequencies will help in choosing a time scale. Use dimensionless time whenever possible. It is the purpose of Part II of these notes to consider the problem of scaling of nonlinear differential equations by considering a number of examples. 2. Example of a Control System, Including Nonlinearities Consider the feedback control system shown in Figure 2.1. Here 9 represents the output of the system and can be considered an angle of shaft rotation. The input is 9i, the desired shaft angle. The purpose of the controller is to make the output angle -0 equal the input angle Gi. To accomplish this, the output is fed back to the input and subtracted from it, producing the error signal e = 9i- 9o. To this signal is added a quantity proportional to the output rate, i.e., -CdGo, in order to provide damping. The resulting signal e - C Pd is -37

-38TACHOMETER d 1-, Cd$ CCP _ COULOUMB FRICTION D C d - TORQUE' C INERTIA PLUS ~ ~~ d C d~ AMPLIFIER VISCOUS DAMP D - Go Figure 2.1 Schematic of Feedback Control System. fed into an amplifier which produces an output torque To = p(e - Cdo), which is applied to the output shaft. In addition there is a torque + C due to the coulomb (dry) friction.* If we assume that the system has inertia I and a viscous-damping constant f referred to the output shaft, then the transfer function relating the total torque Ta applied to the output shaft and G9o the output angle, is simple l/(Ip2 + fp) as shown in the figure, where p is the operator d/dt. Next consider the torque amplifier as a linear device which saturates at the output torque level + Tm. Figure 2.2 shows the input-output characteristic. In order to have a specific problem, let us assume the following values for the parameters in the problem: I = 0.05 slng ft2 Tm = 5 ft# f = 0.25 ft# sec C =0.2 ft# P= 50 ft#/rad Cd = 0.03 sec The equations describing the system can be written from Figure 2.1. Actually, in this case it is just as easy to proceed directly from the figure itself in determining the computer circuit. Without regard for scaling and as a point of departure, consider the differential analyzer circuit shown in Figure 2.3. Here the diode function generator (DFG) following amplifier 2 is used to represent *Actually, the coulomb friction torque is equal and opposite to To for |To| < C, and is equal to TC for G0 positive or negative, respectively, when \To1 > C.

-39OUTPUT +Tm ---- I/' -SLOPE -,u __________"'~!._____/________________INPUT 7 /L -Cd o I / I / Figure 2.2 Torque Amplifier Characteristic the amplifier torque saturation characteristic of Figure 2.2 Pot 1, set at 1/i, is connected in the feedback loop of amplifier 2 to provide an amplifier gain of p for the case where p > 1. For p < 1 the pot would be inserted between amplifier 2 and DFG no. 1. Amplifier 6 with no feedback resistor drives a second DFG the output of which represents the coulomb frinction torque + C.* 2.1 Circuit Using Direct Representation of the Problem Variables Let us redo the circuit in Figure 2.3 with proper scaling. One of the most popular methods of scaling nonlinear problems is to let a convenient computer unit, say one volt, equal unity in each of the problem variables. Scaling is then accomplished by appropriately relabeling each amplifier output. For example, in Figure 2.3 the input 91 and output 9o represent angles, so that we might let 1 radian equal one volt. Then if the full-scale range of input and output. * Actually, a pair of biased diodes can be connected across amplifier 6 so that the amplifier has no feedback (i.e., wide-open gain) until. the output voltage reaches the diode bias voltage, at which point one of the two diodes conducts and the amplifier does not overload.

-400 ~_ _. I_. CO m I -0 _J I I I i i <F U.)

-41angles is + 2 radians, this will represent + 100 volts providing we label the input as 50 91 and the output as 50 GO (see Figure 2.4). In the same way the output of amplifier 1 is 50e. Thus if amplifier 1 has a 20 volt output, 50e = 20 volts = 20 radians, and e = 0.4 radian. We have already seen that summer gains should be kept near unity. But amplifier 2 in Figure 2.3 has a gain of p where p. = 50. We can reduce this to a gain of 2 by relabeling the amplifier output so that it is 1/25th what it otherwise would have been. In Figure 2.4 we see that the output of amplifier 2 becomes (1/25)[-p.(50e - 50Cd6o)] = 2p(e - Cd9o)- Again 1 volt equals one unit of the output variable, i.e., 1 ft lb of torque. Thus an output of 5 volts from amplifier 2 would mean that - 2.(e - Cd9o) = 5 volts = 5 ft lbs, and the unsaturated torque output p.(e - CdGo) = -2.5 ft lbs. If we let the diode function generator have a gain of unity in the unsaturated region, then it should saturate at an input voltage level corresponding to p.(e - Cdo) += Tm = + 5 ft lbs, or 2.(e - CdAo) = + 10 volts. Thus the DFG has unity slope for inputs smaller than 10 volts in magnitude, and saturates at + 10 volts for larger inputs. The output of the DFG represents -2To. We could have let amplifier 2 have a gain of unity by reducing its gain by a factor of p. rather than p./2. But by keeping pot 1 in the circuit we can make small adjustments in p-without repatching the circuit. Since the input to amplifier 3 from DFG no. 1 is -2To, so also is the input from DFG no. 2 given by + 2C, representing the coulomb friction torque. Since C = 0.2 ft lbs,- + 2C = + 0.4 ft lbs = + 0.4 volts. Because the DFG is not very accurate at such low voltage levels, in Figure 2.4 we have let the output be + 4 volts and then added a pot set at 0.1 to attenuate the output. The output of amplifier 3 in Figure 2.4 is 2Ta, where again 1 volt = lft lb. Note that we need a total gain of(l/2I)(50) = 500 through integrators 4 and 5 in order that the output of 5 will be 50 9o. Following the rule of distributing gains (i.e., time constants) approximately equally through integrators we have let integrator 4 have a gain of 25 and 5 have a gain of 20, as shown in Figure 2.4. The setting on pot 2 is then 1/801 = 0.250 while the setting on pot 3 is f/8I = 0.625. The input to pot 4, which sets the derivative feedback constant CD (CD =.0.04 sec) is now -2.5 9G, while the output of pot 4 should be -50 Cd^o. Hence the setting of pot 4 is 20 CD = -9..

-42-I? 0 LO t ~~~~~~~~~~~~~~~~~?I) 1 A~i or 0) o _< 0 0 l\ __..H.~I- - r r^-rO)^ ^ t~ i ~S / \^ ui F ~Z'^ ( L^ )' I WLr<^ -Jg | 1 ~~ d i+ H o$ M a - -<c' - * k'S~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~Lp -^ ^~~~~~~~~~~~~~~~~~~~~~~~I U. I I U~~~~~~~~~~~~~~~~~~~~~~-1

The circuit of Figure 2.4 is now adequately scaled. All amplifier gains are reasonable and amplifier outputs are labeled so that one volt equals unity in the particular output variable. If the time scale is not correct, e.g., if the solution is too fast for the type of recorder being used, the RC time constant of integrators 4 and 5 can be changed. Thus the solution can be run at 1/10 real time by replacing the 0.1 mfd feedback capacitors with 1.0 mfd capacitors. When scaling nonlinear problems using the technique of this section, or, for that matter, using any technique it is important to know the expected full-scale range of the variables. If this is not known in every case, then a reasonable guess must be made, subject to correction after the problem has been set up and several trial solutions run off. In general the scaling has been successful if, during the course of a solution, the output voltage of each amplifier in important computing loops reaches some reasonable fraction of full scale output, normally + 100 volts. This is particularly important at the input to nonlinear elements, since their accuracy is usually the limiting factor in the overall computer accuracy. 3. The van der Pol Equation An example of a famous nonlinear differential equation which represents an interesting scaling problem is the van der Pol equation, which represents the behavior of certain electronic oscillator circuits. After choice of a dimensionless time variable, the van der Pol equation results:* x'- (1 - x2) + x = 0 (5.1) where p is a constant. This can be thought of as a mass-spring-damper equation with unit mass and spring constants and a damping constant which is negative for lxi < 1 and positive for x|J > 1. It turns out that regardless of the initial conditions the solution x reaches a stable limit cycle of amplitude approximately 2. By a stable limit cycle in this case we mean a periodic solution which is always approached regardless of the initial conditions. This is evident in the curves of Figure 3.1, which show differential analyzer solutions for O = 0, 1, and 5. This problem is an example of one where prior knowledge of the general nature of the solution is extremely helpful in scaling the problem. If such knowledge is not available, then one must make some reasonable guesses. Since we know in this case that the solution x * J. J. Stoker, Nonlinear Vibrations, Interscience Publishers, Inc., New York, 1950; pp. 247-252.

-442.0 x =0 5 0 5 0 -1 A =1 Figure -.1 Solutions to v-an der Pol's Equaton

-45will not normally exceed + 2 in magnitude, we can make a choice in voltage scale for x. If we were to follow the procedure used in the previous section, we would let unity in x be one volt. Then 50 x should be the quantity we compute, since x = + 2 would correspond to + 100 volts. However, in this section let us use a slightly different scaling technique. Let x be the quantity we compute rather than 50 x. To have x = + 2 correspond to + 100 volts, then, we should let unity in x be 50 volts. Reference to Equation (3.1) shows that all terms are linear except the term p. x2x. This can be obtained by multiplying x by x and the resulting product by x again. Almost all analog-computer multipliers, whether electronic or servo-type, produce 100 volts output when the two inputs are each 100 volts. Thus a multiplier with inputs x = + 100 volts and x = + 100 volts produces an output of 100 volts. But if unity in x = 50 volts (similarly in x), the product x x should in this case be 2 x 2 = 4, which corresponds to 200 volts if 50 volts equals unity. Thus the actual analog multiplier output of 100 volts corresponds to x x/2.. In general, if unity in each of the problem variables equals V volts, then the output Z of an analog multiplier with inputs X and Y is Z = XY~ —-0. This assumes the computer reference voltage is + 100 volts. This is illustrated in Figure 3.2, which is the circuit for solving van der Pol's equation. Multiplier no. 1 has inputs of -x/5 and x, and since unity = 50 volts, the output is -xx/10. Similarly, the output of multiplier no. 2 is -x x/20. The scaling of this equation is made particularly difficult for large values of p. because the velocity x and acceleration x' go through much larger peak values than for a simple-harmonic oscillator (p = 0). This is evident in Figure 3.1. This is the reason x/10 and -x/5 are computed in Figure 3.2 rather than x and -x, respectively, which would have been more reasonable were it not for the p.(l - x2)x term. Actually, when p. = 5 the circuit of Figure 3.2 is fairly close to being scaled in an optimum fashion, i.e., the outputs of all amplifiers exhibit near full-scale voltage excursion during the solution. For higher values of p. the gains of integrators 2 and 3 must be increased even further, with an accompanying reduction in gain of amplifier 1 for the x input. Had we not known the nature of the solution to van der Pol's equation we would probably have scaled the circuit with unity gain in amplifier 1 for the x input and equal gains (time constants) in amplifiers 2 and 3. For large p. we would have observed overloads in amplifiers 1, 2, and 5, and would have had to rescale the problem in a

L i z N I ~ 0 H o 0 Cd 0 rd i~ LO 0) I x~~~~~~~ -i Ac

-47manner similar to Figure 3.2. Note that pot 4 is placed in the feedback loop of amplifier 5 and is set at 1/2 p. This gives amplifier 5 a gain of 2p.. An alternative approach would have been to place pot 4 between the output of amplifier 5 and the 0.1 meg input resistor to amplifier 1. However this would limit p to a maximum value of 0.5. By giving amplifier 5 a gain of 10, the maximum value of p. could be raised to 5 (pot 4 would be set at p./5); however, by placing pot 4 in the feedback loop of amplifier 5, we can set p. at any value above 0.5, being limited only by the possible voltage saturation at the amplifier output. For the p = 0 case, the output of amplifier 5 into amplifier 1 is disconnected. If the multipliers used for this problem are of the servo type, then the high velocity and acceleration in x for large p. (see Figure 3.1) may exceed the servo capability. Under these conditions the computer time scale can be slowed by decreasing pots 1 and 2 by the same factor. In fact, if these pots are available ganged together on a common shaft (a servo multiplier with multiple ganged pots works well for this purpose), then the time-scale pots always move together and we can slow computer time down or speed it up right in the middle of a solution. If an XY recorder is used with an integrator output on the x axis to give a time sweep, then a third ganged time-scale pot at the integrator input slows down or speeds up the recorder sweep simultaneously with problem time-scale changes. Thus we can slow the problem down in the areas where x changes rapidly and speed it up in other areas so as not to exceed the dynamic capabilities of multipliers 1 and 2, if they are servo driven. 4. Solution of the Satellite and Ballistic-Missile Trajectory Problem An interesting problem which illustrates several techniques in scaling is the problem of the trajectory of a satellite or ballistic missile. We will consider the vehicle a point mass m and will assume that the earth provides a pure central-force gravitational field. As long as we neglect any external forces other than gravity (i.e., neglect any power plant or aerodynamic forces), then the fact that the angular momentum of the body remains constant requires the motion to take place in a plane which contains the center of the earth. Hence it actually becomes a two-dimensional problem, and we choose to use polar coordinates r, 9 to describe the location of the point mass m in the plane of motion. Thus let r be the distance of the body from the center of the earth and O be the polar angle, measured from an arbitrary reference line, as shown in Figure 4.1. Note that the plane of the motion (the xy plane in Figure 4.1) is non-rotating with respect to inertial space, i.e., the earth's surface rotates with respect to this plane.

-48In Appendix A the following y/i i familiar equations are derived: VELOCITY 2 VECTOR i - r I' ) -. rEAT = -cr2 (4.1)./^ ^and r G + 2r (4.2) (. Here go is the acceleration due to gravity at the surface of the earth, and Ro is the radius of the earth (we assume a purely spherical earth for Y > EARTH simplicity). Equation (4.2) actually represents the condition of constant angular momentum, since d/dt (r2 G) = r G9 + 2r r~. Thus Equation (4.2) can be Figure 4.1 integrated to give m r2G = p = constant angular momentum (4.3) From Equation (4.5) -9 = A^...... (4.4) Mrr Equations (4.1) and (4.4) must be solved to determine the trajectory of the point mass m. They are clearly nonlinear. Often the scaling of nonlinear problems is simplified by introducing a dimensionless dependent variable. Here we let the ratio of r/Ro (actual radial distance to earthts radius) be denoted by p. Then Equation (4.1) becomes p -2 2 (4.5) We have seen previously that the scaling is also frequently simplified if we define a dimensionless independent variable, here time. The choice is evident if we write Equation (4.5) as - 4 p 1 (de)2 go dt go dt p B0^

-49where C = Nv2 pO v^/vE. = p-20 Vo/V~ Again all quantities are dimensionless including VgO/vE, the ratio of initial tangential velocity to the escape velocity. Now that all variables, both independent and dependent, are expressed as dimensionless quantities, the scaling is relatively simple. The differential analyzer circuit is shown in Figure 4.2, where we have let 50 volts on the computer be unity. Since the solution is in terms of the polar coordinates p, 9, we need to convert to rectalinear coordinates x and y to record the trajectory on an XY plotter. This can be accomplished using the conventional servo-driven sine-cosine pot shown at the bottom of the figure, where from Figure 4.1 it is evident that X/Ro = p cos 9, y/Ro = p sin 9. In Section 3 we saw that if unity equals V volts for the V input variables X and Y to a multiplier, then the output Z = XY 100. In Figure 4.2 servo-driven pots are used in the feedback circuit of amplifiers 4, 5, and 6. Since the shaft angle is in each case proportional to p, the gain of each amplifier is proportional to 1/pj in effect we have divided the amplifier input voltage by p. In general, if unity equals V volts and the amplifier has an input X which is divided by Y with a servo driven pot in the feedback loop, then the amplifier output voltage Z = X 100. Y V Pot 9 in the feedback loop of amplifier 7 is used to convert 9 in radians to 9 in degrees, where 50 volts = 100~ for the particular servo-driven resolver pots used in this computer. A separate pot 8 is used to set the parameter c independently of pot 9. If a resolver servo is hot available, a circuit employing only multipliers and integrators can be used to obtain p cos 9 and p sin 9 from d9/dr and p (see Appendix B). 5. Solution of Perturbation Equations for the Trajectory Problem The electronic differential analyzer is limited in accuracy and probably can solve the trajectory problem of the previous section to 0.1 - 1% accuracy. But this accuracy is not adequate for detailed trajectory studies, e.g., determination of miss distances for ICBM missiles resulting from small errors in velocity at cutoff. For this reason it is convenient to introduce a perturbation vector A which represents the deviation from some norm trajectory, perhaps the desired trajectory.

-50or 2 2 - = d- (4.6)d 2 p 92 t) la -t p The new dimensionless time is clearly given by T = t (4.7) from which Equation (4.6) becomes d2 de2 1 X K)dT p (4.8) In terms of p and T Equation (4.4) becomes di3 = (4 9) dT m R 2p22 g Next consider the initial velocity vow which makes an angle ac with the line parallel to the earth's surface and has tangential and radial components Vg and vr, respectively, as shown in the figure. The angular momentum p is given by p = mrovoCOs C = m PoRovgo (4.10) - \ \,/V where ro is the initial radius / \.\ / (po = ro/Ro). It can be shown V90 x 0\ ~ so \.> fb Vr that the escape velocity vE at 0\\ / the surface of the earth is arm simply r^or vE = gRW *. (4.11) In terms of vg /vE Equation (4.10) becomes p = m poRo Hf2go Ro E 0 (4.12) which, when substitued into Equation (4.9) yields dG = -2 o _v = c (4.153) dr P- VE p

-512/p f BINS dLP dP dr dr IN.S2 I I I 0.25 1 p tim e s one seon o. _' AceWrT ~ 1/4 (.573) 2 d6e /2 s I pCOs6 = X/Ro NOTE: The pots labeled s set the S4^ __-__ * computer time scale; s seconds of,,-e y Bcos problem time equals one second of computer time. P go,/ Bsin p sing = y/Ro Figure 4.2 Circuit for Trajectory Equations

-52Thus let -4. -4 -4 R = r + A (5.1) where r is an exact and known solution to the trajectory equation for specified initial conditions. Thus r satisfies by definition the equation of motion, which can be written as 2-4.2-, d2r _ gon0 rr (5.2) dt2 r3 </ -\ *<' where we have written er/r2 as r/r3. The equation for R, the position vector from the earth's center to the actual location of the point mass m -is d~2-4~~~ 2-*-~ -*~-d R go~o R + f dt2 R3 where we have included an additional perturbing force f which could include aerodynamic forces, vernier thrust force, forces due to the fact that the earth is not perfectly spherical, etc. Let us calculate first the term R/R3 in terms of r and A. Assume that the perturbation vector A has components Ar along er and A along eQ, i.e., = 4 er + eQ (5.4)

-53Since |A |R| R - r + r A and R = (re~ + Ar er + L Q) (r + A)-5 = (rer + Ar er + eQ) (r5 3r- 4 Ar +...) =r + 1 [-2 Ar e + ] (5.5) r3 r3 Substituting Equations (5.1) and (5.5) into (5.3), we have d2r+ 2 goR2r _goo (-2 r + e) +f dt2 dt2 r3 r 3 (5.6) But from Equation (5.2) d2r goRo r and Equation (5.6) becomes.2 _ g- (-2 er + e) + f (5.7) dt2 r3 This is the perturbation equation, where all of the terms are small compared, with the original terms in Equation (5.3). Next let us calculate d.2/dt2 in terms of Ar and Ag. As noted in Appendix A, the coordinate system with er and eg is rotating with angular velocity Q eZ. Thus dA = A er + e + Q ez x (Ar er + A ) e) dt = (4- A O) er + (g + Ar ) eg d = (Ar A G - O ) er + (' + 4r +r Q) eg + ez x [(. - ) e + ( + ) e] " = (4r - 2 4A -Z5 _' - _ r 02) er

-54If we let the perturbing force f = fr er + fg er, then Equation (5.7) becomes from (5.8) and Ar = 2 A- G + AQ - + tr 2 + 2 0goRo + f (5) adA = - 2 Ar 9 - r G + Ag 2 + f (5.10) r3 These, are the equations to be solved with the electronic differential analyzer. Their solution will give the deviations Ar and. A from the norm trajectory. The equations are linear, since the variables r, 0, and 0 are prescribed functions of time, namely, the solution of Equation (5.2). The equations have time-varying coefficients for this reason. To solve Equations (5.9) and (5.10) on the electronic differential analyzer the circuit of Figure 4.2 will be used to generate the time varying coefficients. Thus the solution of both problems will proceed simultaneously, the solution of Equation (5.2) acting as an input in the solution of Equations (5.9) and (5.10). Equations (5.9) and (5.10) must be rewritten in terms of the dimensionless radius p = r/Ro and time T = go t. We will also assume that fr = fg = 0 in our illustrative problem. For economy of notation we wil write dA/d.T as Ar, dAr/dr as, etc. i.e., from now on all derivatives are with respect to T. Then Ar = 2 + AG @Q + Ar 2 +2A (5.11) 7P *.~~ p Also note from Equation (4.2) that 0 is given by 0e e - a c The electronic differential analyzer circuit is shown in Figure 5.1. Note that unity is equal to 50 volts as before. However, since the problem is now linear as far as Ar and A9 is concerned, the scaling of these variables can be changed at will without changing the circuit.

-555.1 Determination of Initial Conditions for the Example Problem As an example problem we will consider a 6200 statute mile minimum-energy ballistic trajectory (one fourth the circumference of the earth). If Q = 0 at launch, then 0 = 90~ at impact. From other calculations it can be shown that vgo = 0.595 VE at launch and Vr = 0.245 vE at launch, where VE = escape velocity = 2Rogo = 56,700 ft/sec. Now c =-\T2 p0 vQGvE = 0.842 po If we assume for simplicity that the velocity cutoff occurs at the earth's surface, then po = 1 and c = 0.842. - Vr ~2 Vr o0.347 Rogo VE Po = 1, o =~ 0. This provides all of the parameters and initial conditions for the circuit of Figure 4.2 which is used to compute the reference trajectory. As stated earlier the outputs from this circuit (e.g., p, i, p) are used as inputs for the circuit of Figure 5.1 which solves the perturbation equations. An initial value of 2r corresponds to an error in vr at launch, while Ag corresponds to an error in va. Computer solutions for initial Ag and A4 values of 1 volt are shown in Figure 5.2. The resulting Ar and Ag values at impact are tabulated below. Ar = +5.25 volts Initial AG of + 1 volt —- At impact, A = -1.67 volts Ar = +2.71 volts Initial Ar of +, 1 volt - At impact, = -5.15 volts Let us arbitrarily set 1 volt = 1000 ft. Since 43 =^Ji ^ =804 A, dT go dt dt the initial A of 1 volt corresponds to 1000 or 1.24 ft/sec. Similarly the initial A. corresponds to 1.24 ft/sec. Next consider the actual miss at impact corresponding to the Ar and Ag values found from the computer. A sketch of the impact

-56-, _ 2P ), 7 f A I 0.5. 2 A-S7A.t \8 /2 S6A 6 -----— A^A~r-l- -. 1A'C^ "" IAl _o-,. 8 ( O7B (S7A).-Atri 2 I.-ws/- IN S6 I S6C.. s I II 303 L" t-^/S/^-1 *^,r r i Ar/2 I FAigue 5.1 CircAt for Pera Ibation Equations

-57Ae= I volt -2Ar 8 so I volt -4A8 / Ar = I volt 4Ar | Ar= I volt -4Ap -r Figure 5.2 Computer Solutions for Perturbation Equations

-58geometry is shown in the figure, and it is apparent that miss = - + tan of N^~~~~~ " <'Ge __ _ _ _Nl1, N ei -i- - Ar/tana fFor the example, here tan-1 -' = 22.5~ voo and the following miss distances are computed: NA = 1.24 ft/sec, miss = 52 50 - 1670 = 6160 ft tan 22.50 A.r = 1.24 ft/sec, miss = 2710 3 5150 = 3380 ft tan 22.5~ These compare with 6950 ft and 2860 ft respectively when calculated theoretically. Had we chosen initial conditions on AQ and Ar of, say, 50 volts instead of 1 volt, then an order of magnitude increase in accuracy might have been expected. 6. Summary of Nonlinear Scaling Techniques We have seen that in general the scaling of nonlinear problems is more difficult than linear problems, although all of the rules for linear scaling should, whenever possible, be followed. Thus the four rules given in Section 1 still apply. In addition, the following rules are helpful and have been illustrated in the preceding sections. 1. In general it is best to redefine a dimensionless dependent variable or variables, e.g., the dimensionless radius p = r/Ro in Section 4. Under these conditions some convenient voltage level, say 20 volts or 100 volts can be set equal to unity in the dimensionless dependent variable.

-59Another technique is to represent the variables directly in dimensional form and to let a convenient voltage, say one volt, equal unity in each variable. Amplifier outputs are then relabeled to insure proper full-scale voltage excursions, e.g., the example in Section 2. 2. Amplifier gains should be selected so that voltage inputs to nonlinear elements such as function generators and multipliers exhibit full-scale excursion during the problem. 5. If V volts equals unity for the inputs X and Y to a multiplier, then the output Z = XY 1. For a divider the X 100 output Z = Y = 1 assuming that Z = 100 volts when X = Y = 100 volts. 4. When using a diode function generator to generate f(x) it is frequently convenient to let f(x) = fo (x), where fo is a constant equal to the maximum absolute value of f(x) and,6(x) is a dimensionless function with a maximum absolute value of unity. When f(x) is a pure linear function, then 6(x) is simply a linear function of x with a value of unity at the maximum value of x. If f(x) is a particularly simple function, such as the torque-saturating function of Figure 2.2, this procedure is not worthwhile. 5. Often we are concerned with small deviations from a particular solution to a linear or nonlinear differential equation as a result of small changes in initial conditions or forcing functions. In this case the perturbation equations can be solved with the electronic differential analyzer, as in Section 5. 6. When nonlinear elements with limited bandwidth (such as servos) are used, it may be useful to use a "time throttle" consisting of ganged pots ahead of each integrator (see the end of Section 5). This allows the solution to be speeded up when the variables are changing slowly or to be slowed down when the variables are changing rapidly. 7. In general a nonlinear equation cannot be scaled perfectly the first time. It is usually necessary to rescale the circuit in a number of places after running the first solutions.

APPENDIX A DERIVATION OF THE EQUATIONS OF MOTION FOR A POINT MASS MOVING IN A CENTRAL-FORCE GRAVITATIONAL FIELD Consider the equation of motion of a point mass m in the earth's central gravitational force field. The motion will take place in a plane, and we will use polar coordinates r and Q to describe the location of the point mass. First we calculate the expression for the components of acceleration along the direction of increasing r and Q. Let er be the unit vector aligned in the direction of increasing r, and eg be the unit vector aligned in the direction of increasing 0. Then the vector R which describes the position of the point mass with respect to the origin (center of the earth) is given by R = r er. (1) Next consider the rate of change R of the position vector R with respect to the er, eg coordinate system. Since der/dt = 0 in this coordinate system, R is simply = r er. (2) But the inertial forces on m will depend on the acceleration with respect to inertial space. Hence, we must calculate dR/dt and then d2t/dt2, where the derivatives are with respect to inertial space. If we denote eZ as the unit vector perpendicular to er and eg in the figure, and hence, directed out of the plane of the figure, then the er, e coordinate system -60

is rotating with angular velocity QeZ. Thus* dR ~ -4 -4.- t =- R + Oez x R = rer + re (3) = B + ve (3) Similarly d2f dk.. d-4 - = + QeZ x dit dt dt - (e-i + roQeQ + rQeQ) + (rQeQ - rQ2r) or 2 = (r - r42)e + (rg + 2r@)e9 (4) dt2 Here 2rQ is the Coriolis acceleration.** -4 The gravitational force F acting on the point mass will be Fr = 2 7- er (5) r where 7 is the gravitational constant and M is the mass of earth. When r = Ro, where Ro is the radius of the earth, F = -mgoj where go is the gravity acceleration, not including the centrifugal component due to the earth's rotation (the usual g includes this). Thus, Equation (5) can be rewritten as -mgoRo2, (6) r r2 * This equation, Equation (3), states that the velocity with respect to inertial space measured in a moving system of coordinates is equal to the time rate of change of R relative to the moving reference frame plus the cross product of the angular velocity vector 6eZ and the position vector R. This second term 6eZ x R represents the velocity of the coordinate system with respect to inertial,space due to the rotation of the reference frame. ** Equation (4) can also be derived from Lagrange's equation for the point mass m (e.g., see Synge, J. L. and Griffith, B. A. Principles of Mechanics. New York: McGraw-Hill, 1949, 458-462.

-62d2i -, The equation of motion is simply -m - + Fr = 0, or from Equations (4) and (6). dt r rQ2 = - goR () r2 and rQ + 2rQ' = 0 (8) which are the equations which must be solved to determine the trajectory of a ballistic missile or satellite in free flight (no external forces such as thrust or drag)..

APPENDIX B TRIGONOMETRIC RESOLUTION IN ANALOG COMPUTERS BY MEANS OF MULTIPLIER ELEMENTS Many of the problems solved by analog computers require the generation of trigonometric functions. In present general-purpose electronic differential analyzers this is usually accomplished by servo-driven sine-cosine potentiometers or by diode function generators. The purpose of this appendix is to describe an alternative method of generating sine and cosine functions which involves only multipliers as nonlinear elements. Furthermore, this method allows continuous angular travel through an indefinite number of revolutions, whether servo or all-electronic multipliers are used. The new resolving technique is based on the solution of the equation dG2 dQ2 + X =00 ( with initial conditions x(O) = 0, ud (0) = 1. (2) The solution to (1) with the conditions in (2) is simply x = sin @ (3) y = c= os 0 (4) which can then be used to provide coordinate resolutions involving the angle Q. In (l) the independent variable is 0, whereas the electronic differential analyzer (electronic analog computer) can integrate only with respect to time t. But by multiplying the integrand by dQ/dt before integrating, the integration with respect to ~ is effectively accomplished. Thus f z (d-) dt = f zdQ. (5) The electronic differential analyzer circuit for accomplishing the solution to (1) is shown in Figure 1. Note that the input t.o each integrator is multiplied by dO/dt (i.e., 0) before being integrated with respect to time. The multipliers shown schematically in-the circuit could be either servo or all-electronic.

One unit in t - RC seconds Solution. One difficulty with the circuit as it is shown in Figure 1 is the tendency of the solution drift so that cos2 g + sin2 0 / 1. For example, if Q = X = constant, the theoretical solution would be x = sin wt, y = cos ct. Due to power loss in the feedback capacitors the solution may damp somewhat. The net result is that after several cycles x2 + y2 is no longer unity, but perhaps a percent too low. Other computer errors such as dynamic multiplier errors, integrator drifts, etc., can cause similar errors. This error in amplitude can be calculated using multipliers and a summer to obtain 1 - (x2 + y2). Thus define the error by E = 1 - (x2 + y2). (6) The magnitude of x and y can be adjusted automatically upward or downward as needed by providing each integrator with negative or positive damping proportional to. This is shown in Figure 2 where is shown in Figure 2 where is computed by (6) above, multiplied by x and y, and fed back into integrators 2 and 1, respectively. When servo-driven potentiometers are used for multiplication, this resolving circuit requires two servos, each with at least 3 multiplying potentiometers in addition to the reference potentiometer. The shaft angle of each servo is proportional to cos 0 and sin 0, respectively. This means that additional ganged potentiometers on each shaft can be used to multiply cos 0 or sin 0 by other voltages. The output voltages of integrators 1 and 2 are also proportional to cos Q and sin 0, respectively.

-65Figure 2. Complete Computer Circuit fd8or Trigonometric slaves for each master. Additionalt d to multiply cos sinR by other voltages. + 4 and - Q voltages as inputs, providing output voltages and shaft EC-6, No. 2, (June, 1957), 86-92. *.'Hw n[.GlbeWrt,"rgnmercRsluto in,>-X -Analog Corn — pute r b Man o M u ipieEemnts"R ras lctCm r;-3-~~~~~~~~~~~~~~~~~~~~~~.^ EC-6,No.2,(une 95 86-9. -