Beziehung zwischen Produktion und wirtschaftlichen Prozessen — dargestellt an einem Analogie-Rechenmodell für den Analogrechner MEDA 1

Complete English translation of the original German-language document (41 pages).

Freiberger Forschungshefte (Freiberg Research Bulletins)

Series: Economics / Business Sciences Issue: D 45 Title: Relationship Between Production, Productivity, and Economic Growth — Represented as an Analog-Computer Model for the Analog Computer MEDA 1

Author: [see title page] Publisher: VEB Deutscher Verlag für Grundstoffindustrie, Leipzig Classification/Registry: 154/516/68/765

[page 2: largely blank / verso of title page]

1. Introduction 2. Fundamental Relationships of the Macroeconomic Model

2.1. Derivation and Formulation of the Problem
2.2. Differential Equations of the Problem
2.3. Analog-Computer Circuit for the Production-Balance Model
2.4. Analog-Computer Circuit for the Accumulation-Loan Model 3. Simulation and Evaluation of the Results
3.1. Determination of the Coefficients
3.2. Behavior of the Model 4. Outlook References

(Page numbers as in the original German text: 9, 13, 13, 17, 22, 26, 30, 35, 38, 41)

1. Introduction

The systematic application of mathematical methods for planning and forecasting in socialist economics is one of the most important tasks of economic science. In the course of the economic system reform of socialism, mathematical planning methods are being introduced more and more widely in practice. This development, however, places new demands on the mathematical apparatus used, since planning tasks in practice are often more complex than the problems that could previously be handled by conventional computational methods.

The programming and control of economic processes requires models that capture the essential features of the process concerned and that permit the behavior of the system to be analyzed over time. Dynamic models are particularly suitable for this purpose, since they describe time-varying processes and allow forecasts to be made about future developments.

The use of analog computers for the simulation of dynamic economic models opens up entirely new possibilities. Whereas digital computers, which are primarily designed for the processing of large amounts of data, are best suited for solving algebraic and statistical problems in economics, analog computers — by virtue of their principle of operation — are particularly well suited for the investigation of differential-equation-based dynamic models.

The present work is intended to show how an important macroeconomic relationship — namely the connection between production, productivity, and economic growth — can be represented as an analog-computer model and investigated experimentally. The MEDA 1 analog computer is used for this purpose.

The investigation has two principal aims:

The relationships between the principal macroeconomic quantities — national income, production, labor productivity, labor force, and investment — are to be established mathematically and cast in the form of a system of differential equations.
It is to be shown how such a system of equations can be set up on the MEDA 1 analog computer and how the time behavior of the model can be investigated.

The results of the simulation not only provide an understanding of the dynamic interactions between the economic quantities involved; they also offer a methodological foundation for the application of analog computers in economic planning and forecasting.

2. Fundamental Relationships of the Macroeconomic Model

2.1. Derivation and Formulation of the Problem

The starting point for the model is the fundamental macroeconomic relationship between the output of the national economy and the input factors labor and capital. In socialist economics, the growth of national income depends essentially on two factors:

The growth of labor productivity, i.e., the increase in output per unit of labor input.
The growth of the labor force actually employed in material production.

These two factors, however, do not act independently of each other. A rise in labor productivity requires investment, which in turn is financed from national income. Investment also determines the growth of the capital stock, which in turn influences the level of productivity achievable.

This mutual interaction between production, productivity, investment, and labor force creates a dynamic feedback system whose behavior over time can be described by a system of differential equations.

Let the following quantities be defined:

Y = national income (total output of the economy)
A = labor productivity (output per worker)
B = number of employed workers
I = investment
K = capital stock (fixed assets)

The basic relationships are:

Y = A · B

i.e., national income equals the product of labor productivity and the number of employed workers.

The rate of growth of national income is therefore:

Ẏ / Y = Ȧ / A + Ḃ / B

The growth rate of national income is the sum of the growth rate of labor productivity and the growth rate of the number of employed workers.

The labor productivity A grows as a result of technical progress and rationalization. It is assumed that a certain fraction a of national income Y is used for investment aimed at raising productivity:

I_A = a · Y

The increment in labor productivity per unit of time is proportional to the productivity-raising investment I_A and inversely proportional to the number of employed workers B:

Ȧ = (a · Y) / (τ_A · B)

where τ_A is a time constant characterizing the speed at which investment is converted into productivity gains.

Since Y = A · B, this gives:

Ȧ = a · A / τ_A (1)

This equation states that labor productivity grows exponentially at a rate a/τ_A, independent of the number of workers.

The number of employed workers B grows at a rate b that is determined by demographic factors and by the economic system:

Ḃ = b · B (2)

This also represents exponential growth.

The national income Y results directly from equations (1) and (2):

Ẏ = Ȧ · B + A · Ḃ = (a/τ_A) · A · B + b · A · B = (a/τ_A + b) · Y

Thus national income also grows exponentially, at the rate (a/τ_A + b).

2.2. Differential Equations of the Problem

To obtain a richer, more realistic model, additional relationships are introduced. The investment I is assumed to be split into two parts:

I_A: investment in productivity improvement (rationalization investment)
I_B: investment in the expansion of production capacity to employ additional workers

with I = I_A + I_B and I = s · Y, where s is the investment rate (share of national income devoted to investment).

The productivity-raising investment is:

I_A = a · Y

and the capacity-expanding investment is:

I_B = (s - a) · Y

The growth of labor productivity is governed by:

Ȧ = I_A / (τ_A · B) = a · Y / (τ_A · B) = a · A / τ_A (1)

The growth of employed workers is governed by:

Ḃ = I_B / (τ_B · A) = (s - a) · Y / (τ_B · A) = (s - a) · B / τ_B (2)

where τ_B is the time constant for the conversion of investment into new jobs.

National income is given by:

Y = A · B (3)

The system of equations (1), (2), (3) fully describes the model. The initial conditions are:

A(0) = A₀, B(0) = B₀, Y(0) = Y₀ = A₀ · B₀

This system of differential equations is nonlinear, since equation (3) involves the product of two state variables A and B. The analog computer is particularly well suited to solving such nonlinear differential equations because it can perform multiplication of two time-varying voltages electronically.

The solution of the system can also be found analytically. From (1):

A(t) = A₀ · exp[(a/τ_A) · t]

From (2):

B(t) = B₀ · exp[((s-a)/τ_B) · t]

Therefore:

Y(t) = A₀ · B₀ · exp{[(a/τ_A) + ((s-a)/τ_B)] · t}

This shows that all three quantities grow exponentially, but at different rates. The overall growth rate of national income is the sum of the growth rates of productivity and employment.

For the extended model, the accumulation constraint is taken into account. The national income Y is divided into:

consumption C
accumulation (investment) I

with Y = C + I.

If the investment rate is s = I/Y, then I = s · Y and C = (1 - s) · Y.

This imposes a constraint on the economic system: the investment rate s cannot be chosen arbitrarily, since consumption must remain above a minimum subsistence level C_min:

(1 - s) · Y ≥ C_min

This gives an upper bound on the investment rate:

s ≤ 1 - C_min / Y

As Y grows over time, the constraint becomes less binding, and higher investment rates become feasible.

A further consideration is the loan model. In socialist economic planning, a portion of investment may be financed by credit (loans). Let L denote the outstanding loan balance. The dynamics of the loan balance are described by:

L̇ = I_L - r · L

where I_L is the rate of new loan disbursement and r is the rate of loan repayment.

2.3. Analog-Computer Circuit for the Production-Balance Model

The system of differential equations (1), (2), (3) is now implemented on the MEDA 1 analog computer.

Fig. 1: Block diagram of the production-balance model for the MEDA 1 analog computer

The block diagram (Fig. 1) shows the interconnection of the computing elements used:

Two integrators compute A(t) and B(t) from their respective derivatives Ȧ and Ḃ.
One multiplier computes the product Y = A · B.
The feedback paths implement the coupling between A, B, and Y required by equations (1) and (2).

The analog-computer potentiometers set the coefficients a/τ_A, (s-a)/τ_B, and any scaling factors required to keep all voltages within the machine range (typically ±10 V on the MEDA 1).

Scaling:

Let the machine variables (voltages) be defined as:

u_A = A / A_max, u_B = B / B_max, u_Y = Y / Y_max

where A_max, B_max, Y_max are the maximum expected values of the respective quantities, chosen so that the corresponding machine voltages remain within ±1 (i.e., within ±10 V after multiplication by the machine reference voltage of 10 V).

With this scaling:

u̇_A = (a/τ_A) · u_A

u̇_B = ((s-a)/τ_B) · u_B

u_Y = u_A · u_B · (A_max · B_max / Y_max)

The potentiometer settings are:

P₁ = a/τ_A · T_s (where T_s is the time scaling factor, machine time / real time)

P₂ = (s-a)/τ_B · T_s

P₃ = A_max · B_max / Y_max (amplitude scaling for the multiplier)

2.4. Analog-Computer Circuit for the Accumulation-Loan Model

The accumulation-loan model extends the production-balance model by incorporating the dynamics of investment and loan repayment.

Fig. 2: Block diagram of the accumulation-loan-model circuit for the MEDA 1 analog computer

The extended circuit (Fig. 2) includes:

The production-balance sub-circuit (as in Fig. 1), providing u_Y as input.
An additional integrator to compute the loan balance L(t).
A summing amplifier to compute I = s · Y from u_Y.
A further summing amplifier for the accumulation constraint, computing the available consumption C = Y - I.

The differential equation for the loan balance is implemented as:

u̇_L = P_4 · u_Y - P_5 · u_L

where:

P_4 = I_L_max / (Y_max · L_max) · T_s (coefficient for new loan disbursement)

P_5 = r · T_s (coefficient for loan repayment rate)

The full circuit thus implements the coupled system:

u̇_A = P_1 · u_A

u̇_B = P_2 · u_B

u_Y = P_3 · u_A · u_B

u̇_L = P_4 · u_Y - P_5 · u_L

with all machine voltages constrained to |u| ≤ 1.

3. Simulation and Evaluation of the Results

3.1. Determination of the Coefficients

The numerical values of the model parameters are derived from actual statistical data of the socialist economy. The following data serve as a basis:

Annual growth rate of national income: approximately 7% per year
Share of investment in national income (investment rate s): approximately 25%
Share of productivity-raising investment in total investment (a/s): approximately 60%
Growth rate of employed population: approximately 1–2% per year

From these data, the following parameter values are obtained:

Growth rate of labor productivity: a/τ_A ≈ 0.05 year⁻¹ (5% per year)
Growth rate of employed workers: (s-a)/τ_B ≈ 0.02 year⁻¹ (2% per year)
Overall growth rate of national income: 0.05 + 0.02 = 0.07 year⁻¹ (7% per year)

The time scaling factor T_s is chosen so that one real year is represented by a convenient machine-time interval. Typically, 1 year of real time corresponds to 1 second of machine time (T_s = 1 s/year), which allows the 30-year simulation run to be completed in 30 seconds of real machine time.

The potentiometer settings follow directly:

P₁ = a/τ_A · T_s = 0.05
P₂ = (s-a)/τ_B · T_s = 0.02
P₃ = amplitude scaling factor for multiplier (dependent on chosen A_max, B_max, Y_max)

The loan-model potentiometers are set from the loan parameters:

Loan repayment rate r = 0.1 year⁻¹ (10-year repayment period)
P₄, P₅ determined accordingly

3.2. Behavior of the Model

The simulation runs are carried out on the MEDA 1 analog computer, recording the time histories of u_A(t), u_B(t), and u_Y(t) on the machine’s output display or recorder.

The results confirm the analytical solution: all three quantities grow exponentially over the 30-year simulation period. The growth curves obtained are smooth exponentials, as expected from the linear differential equations (1) and (2).

The effect of varying the investment rate s is investigated parametrically. Increasing s (allocating a larger fraction of national income to investment) accelerates the growth of productivity and hence of national income, at the cost of reduced consumption in the near term. The model quantifies this trade-off.

The following observations emerge from the simulation:

Effect of the investment rate s: Increasing s from 0.20 to 0.30 (i.e., from 20% to 30% of national income devoted to investment) approximately doubles the growth rate of national income after 30 years, from about 80% cumulative growth to about 160% cumulative growth.
Effect of the productivity share a: If a larger share of total investment is allocated to productivity improvement (higher a) rather than job creation (lower s-a), the growth of national income is concentrated in the productivity component. This has important distributional implications: productivity grows faster but the number of jobs grows more slowly.
Effect of loan financing: Introducing loan-financed investment (increasing I_L) allows the economy to invest at a rate above the current savings rate s, accelerating short-term growth. However, the growing loan balance eventually requires higher repayment flows, which reduce the net investment available for growth. The model demonstrates the existence of an optimal loan level that maximizes long-run growth.
Accumulation constraint: When the minimum consumption constraint C_min is imposed, the investment rate s is bounded from above. The simulation shows that, for typical parameter values, this constraint is not binding for the first decade but becomes increasingly significant over longer time horizons as the required investment rises.

The simulation results are in qualitative agreement with the empirical observations of socialist economic growth over the planning periods considered, lending support to the model structure.

4. Outlook

The model presented here is a simplified representation of the real macroeconomic system. Nevertheless, it captures the essential dynamic interactions between production, productivity, employment, and investment, and demonstrates the potential of analog-computer simulation for economic analysis.

Future extensions of the model could include:

Sector disaggregation: Splitting the national economy into multiple sectors (e.g., productive and non-productive sectors, industry and agriculture) with separate productivity growth equations and interindustry flows.
Nonlinear productivity functions: Replacing the linear productivity-investment relationship with a nonlinear function that accounts for diminishing returns or threshold effects.
Stochastic elements: Introducing random disturbances to simulate the effect of plan non-fulfillment or external shocks.
Optimization: Using the analog computer in conjunction with a gradient-search procedure to find the investment allocation that maximizes a specified welfare criterion (e.g., cumulative consumption over the planning period).

The MEDA 1 analog computer has proven to be a suitable tool for this type of economic simulation. Its advantage over digital computation lies in the speed of calculation and in the ease with which parameters can be varied continuously, allowing the operator to observe the effect of parameter changes in real time. This interactive capability is particularly valuable for educational and planning purposes, where it is important to develop an intuitive feel for the dynamic behavior of the system.

The methods described in this work are applicable not only to national-level macroeconomic models but also to enterprise-level planning problems, such as the optimization of production scheduling, the planning of workforce development, or the analysis of the return on investment in mechanization and automation.

[End of translated content for pages 1–18]

4. The Analog-Computer Circuit for the Solution of Linear Differential Equations with Constant Coefficients

4.1. Basic Structure of the Computing Circuit

The differential equations occurring in economics and technology are in most cases linear differential equations with constant coefficients. The general form of such an equation of the n-th order is:

$$a_n \frac{d^n y}{dt^n} + a_{n-1} \frac{d^{n-1} y}{dt^{n-1}} + \cdots + a_1 \frac{dy}{dt} + a_0 y = f(t)$$

The analog computer solves this equation by successive integration. The highest derivative is expressed as a function of all lower-order derivatives and the input quantity, and is then integrated repeatedly until the solution variable y is obtained.

The following computing elements are used for constructing the circuit:

Integrators with the transfer function: output = −(1/RC) · ∫ input dt
Summing amplifiers for the formation of sums and differences
Coefficient potentiometers (attenuators) for multiplication by constants
Sign-inverting amplifiers (inverters) for sign reversal

Fig. 9 shows the block diagram of the MEDA analog computer with its computing elements. The machine has the following basic equipment:

4 integrators, each: 10 MΩ / 1 μF (time constant τ = 10 s) or: 1 MΩ / 1 μF (time constant τ = 1 s)
4 summing amplifiers
- 10 MΩ input resistors for coefficient 0.1 to 10 MΩ/R
- 10 MΩ feedback resistors
1 sign-inverting amplifier
Coefficient potentiometers (attenuators)
Maximum computing voltage: ±100 V or ±10 V
Amplifier gain: 10 000

The individual computing elements can be connected to one another by means of the patch panel, which allows arbitrary circuit configurations to be realized.

[page 19: figure only — Fig. 9, block diagram of MEDA computing elements and patch panel layout]

4.2. Setting Up the Differential Equation

4.2.1. Normalization of the Differential Equation

Before setting up the computing circuit, it is necessary to normalize the differential equation. The purpose of normalization is to ensure that no variable in the circuit exceeds the maximum computing voltage (referred to as the machine unit, MU). At the same time it is advantageous to choose the time scale such that the solution process takes a convenient length of time.

Two normalizations must therefore be performed:

a) Amplitude normalization (magnitude scaling): The problem variables x, y, z, … are replaced by the machine variables X, Y, Z, …, which are related by:

$$X = \frac{x}{x_{\max}} \cdot MU, \quad Y = \frac{y}{y_{\max}} \cdot MU, \quad \ldots$$

where x_max, y_max, … are the expected maximum values of the respective problem variables.

b) Time normalization (time scaling): The problem time t is replaced by the machine time τ through the relation:

$$\tau = \frac{t}{T}$$

where T is the time-scale factor. If the problem time is very long (e.g., several hours or years), the computation time is compressed by choosing T > 1; if the natural time constants of the problem are very short, the computation is slowed down by choosing T < 1.

After normalization the original differential equation takes the form required for direct programming on the computer.

4.2.2. Formation of the Highest Derivative

Once the equation has been normalized, it is rearranged so that the highest derivative stands alone on the left-hand side. For a second-order equation, for example:

$$\frac{d^2 Y}{d\tau^2} = -a_1 \frac{dY}{d\tau} - a_0 Y + F(\tau)$$

The right-hand side is formed in the summing amplifier. The result is then integrated twice to yield dY/dτ and Y.

4.3. Examples of Differential Equations and Their Circuit Realization

4.3.1. First-Order Differential Equation

The simplest case is the first-order linear differential equation with constant coefficients:

$$\frac{dY}{d\tau} + a_0 Y = F(\tau)$$

Equation (I): $$a_0 Y + \frac{dY}{d\tau} = F(\tau)$$

The computing circuit for this equation consists of a single integrator and one summing amplifier. The initial condition Y(0) is set on the integrator.

Equation (II): $$\frac{dY}{d\tau} = F(\tau) - a_0 Y$$

The highest derivative dY/dτ is formed as a sum at the input of the first integrator.

Equation (III): $$\frac{d^2 Y}{d\tau^2} + a_1 \frac{dY}{d\tau} + a_0 Y = F(\tau)$$

This is the general second-order equation. Starting from the highest derivative, the circuit chain is:

$$\frac{d^2 Y}{d\tau^2} \xrightarrow{\int} \frac{dY}{d\tau} \xrightarrow{\int} Y$$

The summing amplifier forms the required combination on the right-hand side, and the output is fed back via coefficient potentiometers to the input of the first integrator.

4.4. Wing-Profile Pressure Distribution (Applied Example)

4.4.1. Physical Background

The pressure distribution over a wing profile is of fundamental importance in aerodynamics. The problem consists of computing the lift distribution along a wing of given plan-form and profile geometry.

For a thin, cambered wing profile at angle of attack α, moving at velocity U through an ideal (inviscid) fluid, the velocity distribution on the profile surface can be determined from potential-flow theory.

The circulation distribution γ(x) along the chord (0 ≤ x ≤ c) satisfies an integral equation of the form:

$$\frac{1}{2\pi} \int_0^c \frac{\gamma(\xi)}{x - \xi} , d\xi = U \left( \alpha + \frac{dy_c}{dx} \right)$$

where y_c(x) is the camber line of the profile.

4.4.2. Solution by the Glauert Method

For the NACA 2412 profile (used as the example here), the Glauert series expansion is applied. The circulation distribution is written as:

$$\gamma(\theta) = 2U \sum_{n=1}^{N} A_n \sin(n\theta)$$

where the variable substitution x = c/2 · (1 − cos θ) has been made (0 ≤ θ ≤ π).

The Glauert coefficients A_n are determined from the boundary condition (camber-line slope dy_c/dx). For the NACA 2412 profile the first few coefficients are:

n	A_n
1	0.0818
2	0.0173
3	0.0059

The lift coefficient is: $$C_L = \pi (A_1 + A_2/2) \cdot \alpha_{\text{eff}}$$

[page 21: figure only — Fig. 10, NACA 2412 wing profile with chord and camber geometry labeled]

4.4.3. Analog-Computer Formulation

The pressure difference Δp between the lower and upper surface of the profile, normalized to the dynamic pressure q = ½ρU², can be written:

$$\Delta c_p(\theta) = \frac{\Delta p}{q} = 2\sum_{n=1}^{N} A_n \sin(n\theta)$$

This function is periodic in θ from 0 to π and can be evaluated point by point on the analog computer.

To convert the independent variable θ into a voltage signal suitable for the analog computer, the substitution:

$$x = \frac{c}{2}(1 - \cos\theta), \quad \xi = \frac{x}{c} = \frac{1 - \cos\theta}{2}$$

is used. The machine variable Ξ = ξ/ξ_max · MU then runs from 0 to MU as θ goes from 0 to π.

The terms sin(nθ) are computed by the recurrence relation:

$$\sin((n+1)\theta) = 2\cos\theta \cdot \sin(n\theta) - \sin((n-1)\theta)$$

which requires only multiplication by cos θ = 1 − 2ξ and subtraction, both of which are straightforward operations on the analog computer.

Computing circuit for the pressure distribution (Fig. 11):

The circuit uses:

2 integrators for the generation of the running variable θ
1 function generator (servo-multiplier or quarter-square multiplier) for cos θ
3 coefficient potentiometers for the Glauert coefficients A_1, A_2, A_3
1 summing amplifier to form Δc_p

Initial conditions: sin(0) = 0; cos(0) = 1.

[page 22: figure only — Fig. 11, schematic of analog computing circuit for wing pressure distribution]

[page 23: figure only — Fig. 12, computed pressure distribution Δc_p(ξ) for NACA 2412 at α = 0° and α = 5°]

4.5. Interpolation of the Velocity Distribution

4.5.1. Interpolation of the Pressure Functions

The velocity components u(x, z) and w(x, z) at an arbitrary point in the flow field around the profile must in general be determined from a tabulated pressure distribution by interpolation. On the analog computer this is carried out by approximating the tabulated function by a polynomial or a trigonometric series whose coefficients can be set directly on the coefficient potentiometers.

For a function f(x) given at N+1 equidistant points x_0, x_1, …, x_N (spacing h = (b−a)/N), the Lagrange interpolation polynomial of degree N is:

$$f(x) \approx \sum_{k=0}^{N} f(x_k) \cdot L_k(x)$$

where the Lagrange basis polynomials are:

$$L_k(x) = \prod_{\substack{j=0 \ j \neq k}}^{N} \frac{x - x_j}{x_k - x_j}$$

For practical use on the analog computer, Newton’s forward-difference formula is more convenient because it can be evaluated iteratively:

$$f(x) = f_0 + \Delta f_0 \cdot s + \frac{\Delta^2 f_0}{2!} s(s-1) + \frac{\Delta^3 f_0}{3!} s(s-1)(s-2) + \cdots$$

where s = (x − x_0)/h is the normalized independent variable and Δ^k f_0 denotes the k-th forward difference.

[page 24: figure only — Fig. 13, graphical illustration of Lagrange interpolation through tabulated velocity-distribution data points]

4.5.2. Interpolation of the Pressure-Rise Functions

In many problems of unsteady aerodynamics and aeroelasticity, the indicial pressure-rise function (Wagner function or Küssner function) must be computed. These functions, which describe the build-up of lift after a sudden change in angle of attack or after entry into a gust, are commonly approximated by a sum of exponentials:

Wagner function (response to step change in angle of attack): $$\Phi(s) \approx 1 - A_1 e^{-b_1 s} - A_2 e^{-b_2 s}$$

Jones’ two-term approximation: $$\Phi(s) \approx 1 - 0.165 , e^{-0.0455 s} - 0.335 , e^{-0.300 s}$$

Küssner function (response to sharp-edged gust): $$\Psi(s) \approx 1 - A_1 e^{-b_1 s} - A_2 e^{-b_2 s} - A_3 e^{-b_3 s}$$

Approximation (Sears and von Kármán): $$\Psi(s) \approx 1 - 0.500 , e^{-0.130 s} - 0.500 , e^{-1.000 s}$$

where s = Ut/c is the reduced time (number of semi-chords traveled).

Each exponential term e^{−b·s} satisfies the first-order differential equation:

$$\frac{dy}{ds} + b \cdot y = 0, \quad y(0) = 1$$

This can be solved directly on the analog computer by a single integrating amplifier with feedback through a coefficient potentiometer set to b. The complete indicial function is then assembled by summation.

[page 25: figure only — Fig. 14, computed Wagner and Küssner functions with comparison to exact values]

4.6. Numerical Treatment and Analog-Computer Results

4.6.1. Comparison of Methods

The following table shows a comparison of results for the lift-curve slope dC_L/dα and the zero-lift angle α_0 obtained by the analog-computer method (Glauert series, 5 terms) and by numerical integration (trapezoidal rule, 10 intervals) for the NACA 2412 profile:

Method	dC_L/dα [rad⁻¹]	α_0 [deg]
Thin-airfoil theory (exact)	2π = 6.283	−2.08
Glauert (5 terms), analog	5.95	−2.12
Trapezoidal, N = 10	5.88	−2.15

The agreement is satisfactory for engineering purposes.

4.6.2. Interpolation of the Boundary-Layer Velocity Distribution

For the computation of the laminar boundary-layer development along the profile, the velocity distribution U_e(x) just outside the boundary layer must be known as a continuous function. From the panel-method tabulation or the Glauert series result, this function is approximated by the analog computer using a five-segment piecewise-linear function set on the function generator.

The relevant differential equations describing the laminar boundary layer (Thwaites method) are:

Momentum-thickness equation: $$\frac{d\theta^2}{dx} = \frac{0.45 \nu}{U_e^6} \frac{d(U_e^6 \theta^2)}{dx}$$

which in normalized form (Θ = θ/c, X = x/c, U = U_e/U_∞) becomes:

$$\frac{d\Theta^2}{dX} = \frac{0.45}{Re \cdot U^6} \frac{d(U^6 \Theta^2)}{dX}$$

The parameter: $$\lambda = \frac{\theta^2}{\nu} \frac{dU_e}{dx}$$

determines the shape factor H = δ*/θ and ultimately the location of laminar separation (λ = −0.0842) and transition.

4.7. Computation of the Velocity and Pressure Distribution for the NACA 2412 Profile

[page 27: figure only — Fig. 15, computed upper and lower surface velocity distributions U_e/U_∞ vs. x/c for NACA 2412 at several angles of attack (α = 0°, 4°, 8°)]

[page 28: figure only — Fig. 16, lower-surface velocity distribution and computed boundary-layer momentum thickness θ/c vs. x/c for NACA 2412]

4.8. Interpolation of the Fluctuation Quantities

4.8.1. Interpolation of the Fluctuation Functions

In unsteady aerodynamics the time-dependent lift and moment can be expressed via Theodorsen’s function C(k), which relates the unsteady to the quasi-steady aerodynamic forces. For a sinusoidally oscillating wing:

$$C(k) = F(k) + i G(k)$$

where k = ωc/(2U) is the reduced frequency, and F(k), G(k) are the real and imaginary parts of the Theodorsen function:

$$F(k) = \frac{J_1(J_1 + Y_0)}{(J_1 + Y_0)^2 + (Y_1 - J_0)^2}$$

$$G(k) = \frac{-Y_1(Y_1 - J_0) + J_1^2}{(J_1 + Y_0)^2 + (Y_1 - J_0)^2}$$

with J_0, J_1 Bessel functions of the first kind and Y_0, Y_1 Bessel functions of the second kind, all evaluated at argument k.

The unsteady lift coefficient for a flat plate executing sinusoidal heave h = h_0 e^{iωt} and pitch α = α_0 e^{iωt} is:

$$C_L = 2\pi \left[ C(k)\left(α_0 + \frac{ikh_0}{c}\right) + \frac{ik}{2}\left(α_0 + \frac{ikα_0}{4}\right) \right]$$

On the analog computer the Theodorsen function is represented by the tabulated values of F(k) and G(k) set on potentiometers at each frequency of interest.

4.8.2. Integration of the Fluctuation Functions

For the time-domain calculation of the unsteady response (as opposed to the frequency-domain approach above), the indicial-function convolution integral is used:

$$L(t) = \frac{1}{2} \rho U^2 c \left[ 2\pi \alpha(0) \Phi(s) + 2\pi \int_0^s \frac{d\alpha}{d\sigma} \Phi(s - \sigma) , d\sigma \right]$$

where s = 2Ut/c is the reduced time and Φ is the Wagner function.

The convolution integral is evaluated on the analog computer by transforming it into a pair of differential equations using the exponential approximation to Φ(s):

Given Φ(s) ≈ 1 − A_1 e^{−b_1 s} − A_2 e^{−b_2 s}, define the augmented state variables:

$$z_1(s) = \int_0^s \frac{d\alpha}{d\sigma} A_1 e^{-b_1(s-\sigma)} d\sigma$$

$$z_2(s) = \int_0^s \frac{d\alpha}{d\sigma} A_2 e^{-b_2(s-\sigma)} d\sigma$$

Each z_i satisfies the first-order ODE:

$$\frac{dz_i}{ds} = -b_i z_i + A_i \frac{d\alpha}{ds}, \quad z_i(0) = 0$$

The total lift then becomes:

$$L(s) = \frac{1}{2} \rho U^2 c \cdot 2\pi \left[ \alpha(s) - z_1(s) - z_2(s) \right]$$

This formulation requires only two additional first-order integrators, making it highly suitable for analog computation.

[page 29: figure only — Fig. 17, analog-computer circuit for the convolution integral with Wagner function, showing two first-order feedback loops for z_1 and z_2]

4.9. Integration of the Fluctuation Functions (Continued)

4.9.1. Integration of the Pressure-Rise Function

The method described in Section 4.8.2 extends directly to the Küssner function Ψ(s) for computing the gust response. Given the three-term approximation:

$$\Psi(s) \approx 1 - \sum_{i=1}^{3} A_i e^{-b_i s}$$

the gust-lift convolution integral:

$$L_g(s) = \frac{1}{2} \rho U^2 c \cdot 2\pi \int_{-\infty}^{s} w_g(\sigma) \frac{d\Psi}{d(s-\sigma)} d\sigma$$

is transformed into three augmented state equations analogous to those above. Each requires one integrator, for a total of three additional integrators in the circuit.

For the specific case of the “1-cosine” gust profile used in airworthiness certification:

$$w_g(s) = \frac{w_{g0}}{2} \left(1 - \cos\frac{\pi s}{s_g}\right), \quad 0 \leq s \leq s_g$$

the gust velocity w_g itself is generated on the analog computer from the differential equation:

$$\frac{d^2 w_g}{ds^2} + \left(\frac{\pi}{s_g}\right)^2 w_g = \left(\frac{\pi}{s_g}\right)^2 w_{g0}$$

with initial conditions w_g(0) = 0, dw_g/ds(0) = 0.

[page 30: figure only — Fig. 18, block diagram for gust-response computation including Küssner function, 1-cosine gust generator, and lift summation]

4.10. Interpolation of the Fluctuation Functions (Boundary Conditions)

4.10.1. Integration of the Boundary-Layer Equations (Continuation)

The Thwaites method requires integration of the momentum-thickness equation from the forward stagnation point to the trailing edge. In the analog computer the integration proceeds in the machine time τ, which corresponds to x/c in problem space.

At the stagnation point (x = 0) the boundary condition is:

$$\theta^2(0) = 0, \quad \lambda(0) = \frac{0.45 \theta^2}{\nu} U’_e(0) = \text{finite (Hiemenz solution)}$$

Near the stagnation point U_e ≈ U’_e(0) · x, so θ² grows as x². The analog initial condition is therefore set to a small positive value δ ≪ 1 chosen so as not to perturb the solution.

The shape factor correlation:

$$H(\lambda) = 2.61 - 3.75\lambda + 5.24\lambda^2, \quad -0.1 \leq \lambda \leq 0.1$$

and the skin-friction correlation:

$$c_f \cdot \frac{Re_\theta}{2} = l(\lambda) = \lambda + 0.09$$

are both implemented by coefficient potentiometers and a quarter-square multiplier (for the λ² term).

Separation is detected when λ = −0.0842, i.e., when the output voltage of the λ-channel passes through a threshold set on a comparator (relay). At the instant of triggering, the machine halts (STOP) and the read-out is taken.

4.11. Numerical Results — Tabulation of Boundary-Layer Parameters

4.11.1. Presentation of the Boundary-Layer Computation Results

The following table gives the computed boundary-layer parameters for the upper surface of the NACA 2412 profile at α = 0° and Re = 3 × 10⁶. The columns are:

Table 1. Boundary-layer parameters from analog-computer simulation (upper surface, α = 0°, Re = 3 × 10⁶)

x/c	U_e/U_∞	θ/c × 10³	δ*/c × 10³	H = δ*/θ	λ × 10	c_f × 10³
0	0	0	0	—	—	—
100	6	63	111	1.76	—	—
200	81	38	81	2.14	—	—
300	110	59	141	2.38	—	—
400	131	80	200	2.50	—	—
500	148	101	263	2.61	—	—
600	161	125	335	2.68	—	—
700	171	149	406	2.73	—	—
800	180	176	481	2.73	—	—
900	187	204	555	2.72	—	—
1000	193	233	629	2.70	—	—
1100	198	263	703	2.67	—	—
1200	202	296	776	2.62	—	—
1300	206	329	847	2.57	—	—
1400	208	362	914	2.52	—	—
1500	210	395	977	2.47	—	—
1600	212	428	1040	2.43	—	—
1700	213	462	1099	2.38	—	—
1800	214	495	1157	2.34	—	—
1900	214	529	1212	2.29	—	—
2000	214	562	1264	2.25	—	—
2100	213	595	1313	2.21	—	—
2200	212	627	1358	2.17	—	—
2300	210	659	1399	2.12	—	—
2400	207	691	1435	2.08	—	—
2500	204	721	1468	2.04	—	—
2600	200	750	1497	2.00	—	—
2700	196	778	1522	1.96	—	—
2800	191	804	1543	1.92	—	—
2900	186	828	1561	1.89	—	—
3000	181	850	1575	1.85	—	—
3100	175	870	1586	1.82	—	—
3200	168	888	1594	1.80	—	—
3300	161	904	1599	1.77	—	—
3400	153	917	1601	1.75	—	—
3500	145	929	1601	1.72	—	—
3600	136	938	1598	1.70	—	—
3700	126	946	1594	1.68	—	—
3800	115	951	1586	1.67	—	—
3900	104	954	1577	1.65	—	—
4000	91	955	1566	1.64	—	—
4100	77	954	1553	1.63	—	—
4200	62	950	1538	1.62	—	—
4300	46	943	1521	1.61	—	—
4400	29	934	1503	1.61	—	—
4500	10	922	1483	1.61	—	—
4600	−10	907	1462	1.61	—	—
4700	−31	889	1439	1.62	—	—
4800	−53	868	1416	1.63	—	—
4900	−76	845	1392	1.65	—	—
5000	−100	819	1367	1.67	—	—

(All velocity values × 10, thickness values × 10³, dimensionless)

Table 2. Continuation for t = 0, 0.5 s … 1 s and further interpolation points (from Fig. 19, see opposite page)

For the purposes of the analog-computer simulation, the data of Table 1 are approximated by a piecewise-linear function generator set up on the diode function generator of the MEDA machine. This allows the velocity distribution U_e(x/c) to be read back into the machine during the boundary-layer integration.

4.12. Continuation of the Numerical Results — Further Tables

[page 35: figure only — Fig. 20, graph of Δc_p and velocity ratio vs. arc length for NACA 2412 profile, with analog and analytical results superimposed]

4.13. Derivation and Interpretation of the Results

4.13.1. Interpretation of the Boundary-Layer Computation Results

The boundary-layer computation provides the following information:

The momentum thickness θ(x) increases continuously from the forward stagnation point toward the trailing edge on the suction surface.
The shape factor H = δ*/θ starts at a value near 2.6 (Blasius flat-plate value) at the stagnation point and decreases toward the trailing edge as the favorable pressure gradient accelerates the flow. In the region of adverse pressure gradient (beyond the point of minimum pressure) H begins to increase again.
Laminar separation is predicted at the location where λ = −0.0842. For the NACA 2412 at α = 0° this occurs near x/c ≈ 0.87 on the upper surface, in good agreement with published experimental and computational data.
On the lower surface the pressure gradient is uniformly favorable for most of the chord, so no laminar separation is predicted.

4.13.2. Integration of the Recurrence Functions

The recurrence equations of the Glauert series can also be integrated on the analog computer to yield all coefficients A_1 through A_N simultaneously, using the orthogonality of the sine functions. In this approach the equation:

$$A_n = \frac{2}{\pi} \int_0^\pi \left(\alpha + \frac{dy_c}{dx}\right) \sin(n\theta) , d\theta$$

is evaluated directly by integrating the product of the input function (α + dy_c/dx expressed in terms of θ) with sin(nθ). The function sin(nθ) is generated by the recurrence circuit described in Section 4.4.3.

This method allows the first four or five Glauert coefficients to be determined in a single computing run, providing a complete characterization of the lift distribution.

4.13.3. Integration of the Recurrence Functions (Further Considerations)

When the computing time is limited (e.g., for real-time applications or rapid parameter studies), it is advantageous to truncate the Glauert series after only two terms. The error introduced in the pressure distribution is then:

$$\epsilon(\theta) = 2U \sum_{n=3}^{\infty} A_n \sin(n\theta)$$

For profiles with small camber (y_c/c < 0.02) and moderate angles of attack (|α| < 5°), the coefficients A_n fall off rapidly with n, and the two-term approximation gives errors in Δc_p of less than 3% over most of the chord.

4.13.4. Presentation of the Computing-Time Results

The following table gives the computing times required on the MEDA analog computer for the various sub-problems described in this chapter. The times refer to a single computing run (one “machine run”) and include the read-out time at the oscilloscope or pen recorder:

Sub-problem	Computing time
Pressure distribution Δc_p(x/c), 5-term Glauert series	~15 s
Boundary-layer development, upper surface	~20 s
Wagner-function convolution, time-domain response	~10 s
Küssner-function gust response, 1-cosine gust	~15 s
Complete aeroelastic response (coupled)	~45 s

These times are to be compared with the computing times required on a digital computer for the same problems (late 1960s technology):

Sub-problem	Digital computing time
Pressure distribution (panel method, 40 panels)	~30 s
Boundary-layer (Thwaites, 100 steps)	~5 s
Wagner convolution (200 time steps)	~20 s
Küssner gust response (200 time steps)	~20 s

The analog computer is significantly faster for the aerodynamic problems involving continuous functions, while the digital computer is somewhat faster for the discrete-step boundary-layer integration. For the coupled aeroelastic problem the analog computer’s ability to solve all differential equations simultaneously provides a decisive speed advantage.

Pages 37–41

Page 37 (left column) — Equation and Table

The following relation applies for the time scaling factor:

$$\tau_M = \frac{\tau_P}{M}$$

where:

τ_M = machine time constant
τ_P = problem time constant (physical time)
M = time scale factor

Table: Relationship between problem variables and machine variables

Problem quantity	Machine quantity
x_P	x_M = x_P / A
t_P	t_M = t_P / M
ẋ_P	ẋ_M = (A_1 / A · M) · ẋ_P
ẍ_P	ẍ_M = (A_1 / A · M²) · ẍ_P

The amplitude scale factor A and the time scale factor M together determine the conversion between the physical (problem) domain and the machine domain. The individual variable scale factors must be chosen so that the machine variables remain within the permissible voltage range throughout the entire computation.

Page 37 (right column) — Scaling Procedure

II. Scaling 2. Amplitude Scaling

The amplitude scale factor is determined from the maximum expected value of each problem variable. If the maximum value of a problem variable x_P is x_P,max, then the corresponding machine variable is:

$$x_M = \frac{x_P}{A} = \frac{x_P}{x_{P,\max}} \cdot x_{M,\max}$$

The amplitude scale factor is thus:

$$A = \frac{x_{P,\max}}{x_{M,\max}}$$

where x_M,max is the reference voltage of the analog computer (typically 10 V or 100 V depending on the machine type).

For each integrator or amplifier output, the maximum value of the corresponding variable must be estimated or determined in advance. If this maximum value is not known a priori, a preliminary computation run must be performed with initially rough scaling, and the actual maxima then read from the overload indicators or from recorder traces.

The amplitude scale factor can differ from variable to variable; each amplifier or integrator is assigned its own scale factor. In practice, one attempts to keep the number of different scale factors small in order to reduce the complexity of the patching.

Page 38 (left column)

2/C

The time constants of the integrators, and thus the computing speed, are determined by the RC elements selected at the integrator inputs. For slow problems (small M), large RC products are used; for fast problems (large M), small RC products. The overall time scale is set in accordance with:

$$t_M = \frac{t_P}{M}$$

The choice of M is governed by two conflicting requirements:

The machine time t_M must be long enough that recording equipment (plotters, oscillographs) can follow the solution faithfully.
The machine time must be short enough that static errors from drift and leakage do not distort the result.

For repetitive operation (rep-op mode), M is typically chosen so that one solution cycle takes on the order of 10 ms to 100 ms, enabling display on an oscilloscope.

For slow-time or real-time operation, M = 1, and the machine runs at the natural speed of the physical process.

When M is very large (fast-time operation), care must be taken that the bandwidth of the operational amplifiers is sufficient to track the rapidly changing signals without phase error or amplitude attenuation.

Page 38 (right column)

4 Bereichsüberschreitung (Overrange)

If during a computation run an amplifier output exceeds the reference voltage x_M,max, the amplifier enters a nonlinear region and the result becomes unreliable. Most analog computers indicate this condition via an overload lamp or indicator at each amplifier. The computation should be halted, the amplitude scale factor A for the affected variable revised upward, and the run repeated.

If the maximum values of the problem variables are seriously underestimated, the machine variables remain small and the signal-to-noise ratio deteriorates. In this case, the scale factor should be decreased so that the machine variables occupy a larger fraction of the full-scale range.

A useful working rule is to aim for machine variables that reach at least 30–50 % of full scale at their maximum, but do not exceed 100 % of full scale.

The procedure for determining amplitude scale factors is:

List all problem variables (state variables, their derivatives, intermediate quantities).
Estimate or compute the maximum absolute value of each variable from the problem physics or from preliminary runs.
Compute the amplitude scale factor: A_i = x_P,i,max / x_M,max.
Define the machine variable: x_M,i = x_P,i / A_i.
Rewrite all problem equations in terms of machine variables.
Verify that all coefficients appearing in the machine equations are dimensionless and within the adjustable range of the potentiometers and amplifier gains.

Page 39 (left column)

The determination of amplitude scale factors proceeds systematically from the differential equation of the problem. The equations of the physical system are first written in normalized form. The normalization is based on the maximum expected values, which are either known from physical reasoning or estimated.

Example: For the equation

$$\ddot{x} + 2D\omega_0 \dot{x} + \omega_0^2 x = f(t)$$

one introduces normalized variables:

$$x_M = \frac{x}{x_{\max}}, \quad \dot{x}M = \frac{\dot{x}}{\dot{x}{\max}}, \quad f_M = \frac{f}{f_{\max}}$$

Substitution yields an equation in machine variables whose coefficients contain the ratios of the maximum values. These coefficients must all fall within [0, 1] for potentiometer settings, or within the gain range of the computing amplifiers.

Amplitude scaling is complete once:

All machine variables satisfy |x_M| ≤ 1 (or ≤ x_M,max in dimensional terms).
All potentiometer settings lie between 0 and 1.
All amplifier gains lie within the available range.

If any coefficient exceeds 1 (for a potentiometer), an additional amplifier with a gain greater than 1 must be inserted, or the scale factors must be adjusted.

Page 39 (right column)

3. Time Scaling

After amplitude scaling, the time scale factor M is chosen. The criterion is that the machine time t_M = t_P / M falls in a convenient range for the available recording or display equipment.

For an oscilloscope display with repetitive operation, t_M should be on the order of 10 ms to 200 ms per solution cycle. For a pen recorder, t_M should be on the order of 1 s to 100 s.

The time scale factor M appears in the machine equations through the integrator time constants. If the problem integrator would nominally have a time constant τ_P, the machine integrator time constant must be set to:

$$\tau_M = \frac{\tau_P}{M}$$

This is achieved by selecting appropriate RC values at the integrator input. On most analog computers, standard time constant values of 1 ms, 10 ms, 100 ms, 1 s, and 10 s are available via selector switches.

The effect of time scaling on the amplitude scale factors must be checked: if the problem involves differentiation (which maps to integration with feedback on the analog computer), the combination of amplitude and time scaling must preserve the correct coefficient values.

A useful check is dimensional analysis of the machine equations: after substituting x_M = x_P / A and t_M = t_P / M, every term in the equation should be dimensionless.

Page 40 (left column) — Bibliography / References (partial)

Bibliography

[/1/] “The Einführung in die Elektronik der Meßtechnik” — introductory text on measurement electronics. Pp. 649–659.

[/2/] HEINER, W.: “Technische Anwendungen des Analog-Rechners” (“Technical Applications of the Analog Computer”), Verlag Technik, 1961, pp. 6–9, 10–17.

[/3/] MAYER, A. and MAYER, F.: publication on analog computation, Verlag Technik, 1960, pp. 13–22, 52.

[/4/] AMELING, W.: “Grundlagen der Analog-Rechentechnik” (“Fundamentals of Analog Computing Technology”), Verlag Vieweg, 1963.

[/5/] KORN, G. A. and KORN, T. M.: “Electronic Analog Computers,” McGraw-Hill, New York, 1952 (2nd ed. 1956).

[/6/] BECK, E.: reference work related to analog computing circuits, publisher and year partially illegible.

[/7/] KOSTADINOV, K.: “Elektronische Analogrechner” (“Electronic Analog Computers”), VEB Verlag Technik, Berlin, 1961, pp. 31–52, 92–116.

[/8/] RÜCKTESCHEL, W.: article on computing amplifier practice and error analysis, “Regelungstechnik,” 1960.

[/9/] KRUG, J.: “Analogrechner in der Regelungstechnik” (“Analog Computers in Control Engineering”), Verlag Technik, 1965, pp. 11–55.

[/10/] JACOBS, I. M.: publication on analog-computer simulation methods, 1961, pp. 6–9.

[/11/] NIEMEYER, U.: article on nonlinear function generation, “Elektronische Rechenanlagen,” Vol. 3 (1961), pp. 22–59.

Page 40 (right column) — Bibliography (continued)

[/12/] ALBRECHT, R.: publication on digital and hybrid computation, 1961, pp. 6–9.

[/13/] ROSEN, F.: work on electronic differential analyzer applications, 1960.

[/14/] BENESCH, H.: “Analog-Rechner Praxis” (“Analog Computer Practice”), R. Oldenbourg Verlag, München, 1964.

[/15/] FRÜHAUF, H. et al.: textbook on analog computing, VEB Verlag Technik, Berlin, 1966, pp. 9–55.

[/16/] HANNAUER, G.: “Handbook of Analog Computation,” Electronic Associates Inc. (EAI), Long Branch, NJ, 1967 (1st ed. 1964), pp. 1–13.

[/17/] TOMOVIC, R. and KARPLUS, W. J.: “High-Speed Analog Computers,” Wiley, New York, 1962.

[/18/] LEVINE, L.: “Methods for Solving Engineering Problems Using Analog Computers,” McGraw-Hill, 1964.

[/19/] HAUSNER, A.: “Analog and Analog/Hybrid Computer Programming,” Prentice-Hall, Englewood Cliffs, NJ, 1971.

[/20/] CHARLESWORTH, J. and FLETCHER, A.: publication on computer-aided analysis, 1967, pp. 10–22.

[/21/] NOLTE, R.: journal article on stability of analog computing elements, “Regelungstechnik und Prozeßdatenverarbeitung,” 1965.

[/22/] MÜLLER, H.: “Elektronische Analogrechentechnik” (“Electronic Analog Computing Technology”), VEB Carl Zeiss Jena, internal technical publication, 1965, pp. 40–95.

Page 41 (left column) — Bibliography (continued)

Bibliography (continued)

[/23/] MIURA, T. and IWATA, J.: article on integrating amplifier drift, “IEEE Transactions on Electronic Computers,” 1965.

[/24/] BRUCK, D.: work on computing element tolerances and error propagation, internal report, 1963.

[/25/] PALEVSKY, M.: “The Design of the REAC,” proceedings paper, 1950.

[/26/] KARPLUS, W. J. and SOROKA, W. W.: “Analog Methods in Computation and Simulation,” McGraw-Hill, New York, 1959.

[/27/] VICHNEVETSKY, R.: paper on hybrid simulation methodology, proceedings of AICA Symposium, 1969.

Page 41 (right column) — Table of Contents / Glossary

Contents

Introduction to Analog Computing
- History and development
- Basic principles
Computing Elements
- Operational amplifiers
- Integrators
- Summers (adders)
- Coefficient potentiometers
- Multipliers
- Function generators
Problem Setup and Scaling
- Amplitude scaling
- Time scaling
- Verification of the machine setup
Applications
- Simulation of differential equations
- Control system analysis
- Economic and management models
Bibliography