Die Lösung und Darstellung von partiellen Differentialgleichungen auf dem elektronischen Analogrechner

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

Solution and Representation of Partial Differential Equations on the Electronic Analog Computer

Electronic Data Processing — Technical Reports on Programmable Machines and Their Applications, Issue 5/1962

Editors: Dr. H. K. Schuff (chief editor), Dr. Hans Christien (Hamburg), Dr. Günter Doetsch (Freiburg i. Br.), Prof. Dr. W. Haack (Berlin), Dr. W. Händler (Erlangen-Bochum), N. D. Hill (Hayes, England), Prof. Dr. C. Knott (Loughborough, England), Dr. G. Leitler (Neubiberg/Obb.), Prof. Dr. W. v. Murjahn (dir., Amsterdam)

Working group (Fachausschuss): Prof. Dr. W. Händler (chair), Dipl.-Ing. K. W. Braunbeck, O. Haubner, K. Manne, K. Meyberg, Dr. K. Schwarzelbach, Dipl.-Ing. A. Wohlmuth

The Solution of Partial Differential Equations and Their Representation Possibilities on the Electronic Analog Computer

M. St. Albring, Aachen

Abstract: A brief general survey of the state of the electronic analog computer art is given, and the programming procedures and methods of solution for partial differential equations on the electronic analog computer are treated, because the properties and limitations of the analog computer make it necessary, when dealing with such problems, to apply special problem-adapted numerical methods. The solution of some example problems illustrates the discussion.

1. The Programming of the Electronic Analog Computer

Whereas previously the solution of technical problems was accomplished mainly by graphical methods or mathematical analysis, numerical evaluation, using mechanical computing aids, has increasingly appeared alongside these. At the same time, much more complex problems and more exact results can be achieved, and the computing times required can be substantially shortened. The development of the electronic analog computer has brought about a further essential advance in this area. The analog computer makes it possible in many cases to solve problems arising from the theory and simulation of technical and physical processes both rapidly and elegantly. Beyond this, several variants of a problem can be examined in rapid succession and, if needed, the results can also be processed directly.

The programming of an analog computer consists of developing an analog circuit from the system of equations describing the problem. This circuit comprises interconnected building blocks, which are arranged in such a way that they perform the required mathematical operations. The individual building blocks are described in the following section.

2. The Programming of Differential Equations on the Electronic Analog Computer

2.1 Programming Procedure

The basic programming procedure consists of the following steps:

Formulation of the mathematical relationships describing the physical problem (differential equations, boundary conditions, initial conditions) in normalized form.
Setting up of the analog-computer circuit by translation of the mathematical relationships into the circuit diagram. For the conversion from mathematical operations into circuit elements, the symbols in Table 1 are used.
Verification of the circuit diagram for correct implementation of the mathematical equations (in particular, attention to sign).
Scaling for amplitude and frequency.
Setting the potentiometer coefficients and the initial conditions.
Running the computation.

3. The Building Blocks of the Analog Computer

Table 1 (page 3) gives an overview of the elementary building blocks. The columns show: designation, symbol, circuit diagram, and mathematical description.

Designation	Symbol	Circuit	Mathematical Description
Summing amplifier / Summing inverter	[symbol]	[circuit]	y = −(x₁ + x₂ + … + xₙ); in general: y = −∑ xᵢ, with 0 ≤ k ≤ 1
Potentiometer	[symbol]	[circuit]	y = k · x, with 0 ≤ k ≤ 1
Coefficient setter	[symbol]	[circuit]	y = −k · x, with 0 ≤ k ≤ 1
Inverter (sign changer)	[symbol]	[circuit]	y = −x
Summing integrator	[symbol]	[circuit]	y = −∫(x₁ + x₂ + … + xₙ) dt + y₀; or y = −∫∑xᵢ dt + y₀, with 0 ≤ k ≤ 1
Tracking/storage integrator	[symbol]	[circuit]	[as summing integrator, with additional hold mode]
Multiplier	[symbol]	[circuit]	z = W · [f₁(x₁, y₁) + f₂(x₂, y₂)]; also: z = W · [f(x, y)]; electronically or mechanically
Function generator	[symbol]	[circuit]	a.o. relay networks, diode circuits, electronic construction; solution via Diode switch or similar; y = f(x)
Switching element / Comparator	[symbol]	[circuit]	y = f(x)

Fig. 1 — Symbols for the building blocks (Table 1)

[page 3: figure only — table of building-block symbols with circuit diagrams]

4. Programming the Variables

4.1 The Method of Treating the Variables

After this overview of the building blocks, the method of treating the variables (integration variables) in a circuit must be explained. On the analog computer, the solution of differential equations is achieved by integration, which is distinguished from differentiation because it does not amplify noise.

For an n-th order differential equation, one needs n integrators. As a starting point for building up the circuit, it is assumed that the highest derivative d^n y/dt^n is available. From this, through successive integration, all lower derivatives and finally y itself are obtained. The circuit is closed by feeding the computed lower-order terms back according to the differential equation.

Example: For the second-order differential equation

d²y/dt² + a₁ dy/dt + a₀ y = f(t)

one assumes that d²y/dt² is available. Then:

dy/dt = −∫(d²y/dt²) dt + (dy/dt)₀

y = −∫(dy/dt) dt + y₀

Feeding back according to the equation gives:

d²y/dt² = −a₁ dy/dt − a₀ y + f(t)

The complete circuit uses two integrators, two coefficient setters or potentiometers, and one summing amplifier.

4.2 The Method of Treating the Eigenvalues

As a supplement to the method of treating the variables (described above), the method of treating the eigenvalues should be mentioned. In this method, the constants of the differential equation are not set directly, but rather the characteristic roots (eigenvalues) of the equation are set. This is particularly advantageous when it is desired to change the type of oscillation (aperiodic, critically damped, oscillatory) in a systematic manner.

For a second-order system with the characteristic equation

λ² + 2δω₀ λ + ω₀² = 0

the parameters δ (damping ratio) and ω₀ (natural frequency) characterize the behavior. The eigenvalues are:

λ₁,₂ = −δω₀ ± ω₀√(δ² − 1)

The circuit parameters can be set directly in terms of δ and ω₀ instead of in terms of the equation coefficients. This simplifies parameter studies.

5. Scaling

5.1 Amplitude Scaling

On the electronic analog computer, the voltage at any point in the circuit may not exceed a certain maximum value (the machine unit, typically ±100 V or ±10 V), since otherwise overloading (clipping of the signal) occurs and the computation results become falsified. On the other hand, the signals should be as large as possible to achieve good signal-to-noise ratio. A good rule of thumb is that all quantities in the problem should be scaled to lie within the range of ±1 machine unit, so that maximum excursions approximately reach but do not exceed the full machine range.

For the variable xᵢ of the problem, the machine variable Xᵢ is defined by:

Xᵢ = xᵢ / xᵢ,max

where xᵢ,max is the anticipated maximum value of xᵢ. This scales all variables to the range [−1, +1].

5.2 Time Scaling (Frequency Scaling)

Besides amplitude scaling, time scaling is also necessary in many cases. If the time constant of the physical process is very small (fast process) or very large (slow process), a time-axis transformation must be performed:

τ = t / tₘ

where tₘ is a suitable time-scaling factor (machine time unit).

For a fast process one stretches the time axis (tₘ < 1), and for a slow process one compresses it (tₘ > 1). This allows the computation to run within a convenient time window on the machine.

Table 1 (referred to above) — detailed table with scaling formulas for typical circuit elements and transformations for rescaling by factors of 10:

n	Differential expression	Scaled to a smaller time constant	Scaled to a larger time constant
1	dy/dt	−(a₀/β) · Y	− (a₀ · β) · Y
2	d²y/dt²	+(a₀/β²) · Y	+(a₀ · β²) · Y
…	…	…	…

(Full table on page 8.)

6. Partial Differential Equations

The handling of partial differential equations (PDEs) on the analog computer requires special treatment because, unlike ordinary differential equations (ODEs), they involve derivatives with respect to more than one independent variable.

6.1 Formulation

In general, the problem at hand is a PDE of the form

F(x, y, t, u, ∂u/∂x, ∂u/∂y, ∂u/∂t, ∂²u/∂x², …) = 0

with associated boundary conditions and initial conditions.

The following discussion concentrates on the two-dimensional heat equation and similar parabolic PDEs, as well as on elliptic boundary-value problems, since these occur most frequently in technical applications. Hyperbolic equations (wave equation) are also treated.

6.2 Replacement of the Partial Differential Equation by Ordinary Differential Equations

Finite-difference method (spatial discretization)

The most commonly applied method for solving PDEs on the analog computer is the method of lines (also called the spatial finite-difference method). In this method, the spatial derivatives are replaced by finite-difference approximations, reducing the PDE to a system of ODEs in the remaining independent variable (usually time), which can then be solved by the analog computer in the usual way.

For the heat conduction equation

∂u/∂t = a² · ∂²u/∂x²

the spatial derivative is approximated by:

∂²u/∂x² ≈ (uᵢ₊₁ − 2uᵢ + uᵢ₋₁) / (Δx)²

This leads to the system of ODEs:

duᵢ/dt = (a²/(Δx)²) · (uᵢ₊₁ − 2uᵢ + uᵢ₋₁), i = 1, 2, …, N−1

with boundary values u₀ and uₙ specified.

This system of N−1 ODEs is then programmed on the analog computer in the usual manner. Each internal grid point requires one integrator and the corresponding summing and coefficient-setting elements.

For a two-dimensional PDE

∂u/∂t = a² (∂²u/∂x² + ∂²u/∂y²)

the spatial discretization in both x and y gives a system of ODEs:

duᵢⱼ/dt = (a²/(Δx)²)(uᵢ₊₁,ⱼ − 2uᵢⱼ + uᵢ₋₁,ⱼ) + (a²/(Δy)²)(uᵢ,ⱼ₊₁ − 2uᵢⱼ + uᵢ,ⱼ₋₁)

The number of ODEs required is (N−1)×(M−1) for an N×M grid.

6.3 Accuracy of the Finite-Difference Method

The error introduced by replacing spatial derivatives with finite differences is of order (Δx)² for central differences. Specifically, the truncation error in the spatial discretization is:

ε ≈ (Δx)²/12 · ∂⁴u/∂x⁴

Thus, reducing the grid spacing Δx by a factor of 2 reduces the discretization error by a factor of 4.

However, reducing Δx also increases the number of grid points (and thus integrators needed), and it stiffens the ODE system (the largest time constant is proportional to (Δx)²), which may require a smaller time step and slower computation on the analog machine. There is therefore an optimum balance between accuracy and computational complexity.

7. Representation of Partial Differential Equations

[page 11: figure only — circuit diagrams and block diagrams for the heat equation with spatial discretization]

The circuit for the one-dimensional heat equation with N−1 internal grid points is set up as follows:

One integrator per internal grid point i.
Each integrator receives inputs from the potentiometer-weighted values of uᵢ₋₁, uᵢ, uᵢ₊₁.
The coefficient (a²/(Δx)²) is set on the potentiometers.
Boundary conditions are set as fixed voltages (for Dirichlet conditions) or through additional circuitry (for Neumann conditions).

The solution proceeds in real time (or scaled time), and the output at each integrator represents the temperature (or other field variable) at the corresponding grid point as a function of time.

8. Estimation of the Error

Before treating further solution methods, the sources and estimation of error in the analog-computer solution should be examined. Errors arise from the following causes:

a) Static computing errors — arising from amplifier offset voltages, finite open-loop gain, and resistor/capacitor tolerances. b) Dynamic computing errors — arising from the finite bandwidth of operational amplifiers and phase errors at high frequencies. c) Truncation (discretization) errors — arising from the finite-difference approximation of spatial derivatives. d) Quantization and setting errors — in the potentiometer settings and initial conditions.

For the evaluation of the discretization error, the following consideration applies:

For the difference quotient

(uᵢ₊₁ − 2uᵢ + uᵢ₋₁) / (Δx)²

the Taylor series expansion gives:

uᵢ±₁ = uᵢ ± Δx · u’ + (Δx)²/2 · u” ± (Δx)³/6 · u''' + (Δx)⁴/24 · u'''' + …

Adding and rearranging:

(uᵢ₊₁ − 2uᵢ + uᵢ₋₁)/(Δx)² = u” + (Δx)²/12 · u'''' + …

The leading error term is thus (Δx)²/12 · u'''', which is second order in Δx.

This error can in principle be reduced by using higher-order difference schemes. For example, a fourth-order accurate approximation is:

∂²u/∂x² ≈ (−uᵢ₊₂ + 16uᵢ₊₁ − 30uᵢ + 16uᵢ₋₁ − uᵢ₋₂) / (12(Δx)²)

However, implementing this on the analog computer requires more complex wiring and additional potentiometers.

9. The Approximation of the Error

To estimate the influence of the truncation error on the computation result, the following approach is used. The exact solution u(x,t) and the computed (discrete) solution uᵢ(t) differ by an amount that depends on the grid spacing Δx and on the higher derivatives of the solution. If the solution is smooth, the error is small; if u has large fourth derivatives, a finer grid is needed.

For a specific example — the heat conduction equation with a sinusoidal initial condition u(x,0) = sin(πx/L) — the exact solution is known:

u(x,t) = sin(πx/L) · exp(−a²π²t/L²)

The discrete solution with N grid points gives an approximate eigenvalue:

λ_h = −(4a²/(Δx)²) · sin²(πΔx/(2L))

which for small Δx/L approaches the exact eigenvalue −a²π²/L². The relative error is:

(λ − λ_h)/λ ≈ (πΔx/L)²/6

This shows that with N = 10 grid points (Δx = L/10), the relative eigenvalue error is approximately 1.6%, while with N = 5 the error is about 6.6%.

9.1 Representation of Linear Functions

Every first-order ordinary differential equation can in principle be programmed on the analog computer. It should also be noted that a first-order ODE can arise in the context of PDE solution, representing the evolution of the field variable at one grid point.

[page 14: figure — block diagram for a first-order system with labeled inputs and outputs]

For each internal grid point, the following circuit is required:

One integrator
One summing amplifier (for combining contributions from neighboring grid points)
Two potentiometers (for setting the coefficient a²/(Δx)²)

[page 14: additional figure — block diagram showing three coupled integrators for a three-point discretization]

The three integrators are coupled through summing amplifiers such that each integrator output feeds as input into its neighbors, multiplied by the appropriate coefficient.

From the approximation of the functions one gets, for example for y = A · x:

y / y_max = (A · x_max / y_max) · x / x_max

i.e., the potentiometer must be set to:

k = A · x_max / y_max

9.2 Representation of Staircase Functions

Each staircase function can be approximated arbitrarily closely by a combination of shifted step functions. On the analog computer, step functions are generated using switching elements (comparators, relays). For each step in the staircase, one switching element and one coefficient potentiometer are required.

[page 14: figure — circuit for a staircase function with three steps]

9.3 Linear Interpolation between Boundary Values

For Dirichlet boundary conditions in which the boundary value varies linearly with time — e.g., the temperature at the boundary increases linearly — the boundary-value generator need only be a linear ramp, which is obtained from a constant-input integrator.

For time-varying boundary conditions of more general shape, a function generator is required to produce the specified boundary-value time history.

[page 14: figure — schematic showing time-varying boundary value with function generator]

10. Error Estimation for the Particular Differential Equation

The behavior of the solution — and hence also the error — depends essentially on the character of the PDE (parabolic, elliptic, or hyperbolic). For parabolic equations (heat equation), the solution is smoothed by the diffusion process, so that discretization errors tend to decrease over time. The step-by-step checking of whether truncation errors remain small at all grid points during the computation is strongly recommended.

For the estimation of the maximum allowable grid spacing Δx, the following practical rule is available: If the solution u(x,t) is known to have significant spatial frequency content up to a wavenumber k_max, then the grid spacing should satisfy:

Δx ≤ π / (5 · k_max)

to keep the discretization error below approximately 2%.

10.1 Generation of Lösung Functions (y = A · x^m)

For the solution of nonlinear PDEs, or for the calculation of nonlinear source terms, function generators on the analog computer can produce the approximation y = A · x^m for m ≠ 1. These are typically realized using a combination of diodes and resistors (diode function generator, or DFG), which piecewise-linearizes the nonlinear curve.

[page 13: figure — diode function generator circuit with breakpoints]

The accuracy of the DFG depends on the number of diode segments used; with n segments, the approximation error is of order 1/n² in the worst case.

11. Parabolic Differential Equations — Heat Conduction

The heat-conduction equation is the prototype parabolic PDE:

∂u/∂t = a² · ∂²u/∂x²

with initial condition u(x,0) = f(x) and boundary conditions (Dirichlet or Neumann).

Circuit setup:

The spatial interval 0 ≤ x ≤ L is divided into N equal subintervals of width Δx = L/N. The internal grid points are at x₁, x₂, …, x_{N−1}. Each grid point xᵢ corresponds to one integrator on the analog computer.

The ODE for grid point i is:

duᵢ/dt = (a²/(Δx)²) · (uᵢ₋₁ − 2uᵢ + uᵢ₊₁)

Amplitude scaling:

The variable uᵢ is scaled to Uᵢ = uᵢ / u_max so that Uᵢ ∈ [−1, +1]. The coefficient becomes:

a²/(Δx)² · u_max / U_max = a²/(Δx)²

(unchanged if all variables are scaled uniformly).

Time scaling:

The natural time constant of the heat equation is τ = L²/(a²π²). If this is, say, 0.1 s, a time scale factor of β = 100 stretches the computation so that one machine second corresponds to 0.01 s of physical time.

Initial conditions:

At t = 0, the initial voltage on each integrator capacitor is set to Uᵢ(0) = f(xᵢ)/u_max using the initial-condition setting facility.

Boundary conditions:

Dirichlet: The boundary voltages U₀(t) and U_N(t) are fed as fixed (or time-varying) voltage sources into the neighboring integrators’ summing inputs.
Neumann: A zero-flux boundary condition ∂u/∂x = 0 at the boundary is enforced by setting u₋₁ = u₁ (image condition), which requires an additional summing connection.

[page 11: figure — complete circuit for the 1D heat equation with 4 internal grid points]

11.1 Two-Dimensional Heat Equation

For the two-dimensional case the computation requirement grows rapidly. With a 5×5 grid (25 internal points), 25 integrators are needed, which is at the limit of a medium-sized analog computer.

The coupling equation for interior point (i,j) is:

duᵢⱼ/dt = (a²/(Δx)²)(uᵢ₊₁,ⱼ + uᵢ₋₁,ⱼ − 2uᵢⱼ) + (a²/(Δy)²)(uᵢ,ⱼ₊₁ + uᵢ,ⱼ₋₁ − 2uᵢⱼ)

For Δx = Δy this simplifies to:

duᵢⱼ/dt = (a²/(Δx)²)(uᵢ₊₁,ⱼ + uᵢ₋₁,ⱼ + uᵢ,ⱼ₊₁ + uᵢ,ⱼ₋₁ − 4uᵢⱼ)

[page 16: figures — schematic and block diagram for the 2D heat equation on an 8×8 grid system using the analog computer with plotters]

11.2 Hyperbolic Differential Equations — Wave Equation

The one-dimensional wave equation

∂²u/∂t² = c² · ∂²u/∂x²

is handled similarly, but now requires two integrations in time instead of one for each grid point. The semi-discretized system is:

d²uᵢ/dt² = (c²/(Δx)²)(uᵢ₊₁ − 2uᵢ + uᵢ₋₁)

This requires two integrators per grid point (for uᵢ and duᵢ/dt).

Stability: The explicit time-stepping implied by the analog computer integration is stable for the wave equation provided (for the continuous-time ODE system) all eigenvalues are purely imaginary, which is the case. However, errors in the potentiometer settings and amplifier offsets can cause slow drift.

12. Control and Display of Functions

12.1 Representation of Functions y₁ = y₁(x), y₂ = y₂(x)

In many cases it is not enough to solve the PDE numerically; the solution must also be displayed and recorded. The following output devices are available on the analog computer:

Oscilloscope (X-Y or time-sweep)
X-Y plotter (curve plotter)
Y-t recorder (strip-chart recorder)
Printer (point-by-point)

For displaying the spatial distribution of the field variable at a fixed instant of time, the X-Y plotter is used: the x-coordinate is swept and the output voltages of the individual integrators are connected to the Y-input one after another (manually or via a multiplexer).

For displaying the time evolution at a fixed spatial point, the Y-t recorder is used, with the integrator output connected directly.

[page 16: figures — photographs of recorded outputs; strip charts showing temperature distributions as a family of curves for different times]

[page 17: figures — three-dimensional surface plots of the solution u(x,t) showing evolution of a heat distribution over a two-dimensional spatial domain, obtained by composite recording on the X-Y plotter]

12.2 Representation of Functions y₁(x), y₂(x), y₃(x) for Parameter Studies

For a parameter study (varying a coefficient such as a² or boundary conditions), the computation is run repeatedly, once for each parameter value. The family of solution curves is recorded on the same plotter sheet by allowing successive runs to overprint on the same chart. This produces a family of curves parameterized by the chosen variable.

[page 17: figure — family of curves showing parametric variation, with curves labeled by the parameter value]

[page 18: figure — three-dimensional mesh plot showing the solution u(x,y) of the 2D problem, obtained by raster scanning]

13. Summary

After this survey of the analog computer and its possibilities for the solution of partial differential equations, the following further methods should be mentioned:

The methods described here — in particular the finite-difference method (method of lines) applied to the spatial coordinates — are the most widely used on analog computers. The methods were illustrated through the heat-conduction equation as a typical example.

The second Heft (volume) will be concerned with the solution of elliptic differential equations by various relaxation and iteration methods suited to the analog computer.

The approximation method for treating PDEs on analog computers can be used generally: for linear and nonlinear PDEs alike, provided the nonlinear terms can be handled by the available function generators and multipliers.

The principal advantages of the analog computer for PDE solutions are:

Rapid setup and computation of families of solutions (parameter studies)
Direct display of the solution in real time or scaled time
Intuitive relationship between the physical process and the computing circuit

The principal limitations are:

The number of grid points (and hence accuracy) is limited by the number of computing elements available
Amplitude and time scaling must be carefully set to avoid overload and to keep the computation within the machine’s voltage and time ranges
The inherent computing errors (amplifier drift, component tolerances) impose a lower bound on the achievable accuracy

Despite these limitations, the electronic analog computer has proven to be a valuable tool for the rapid exploration of the solution behavior of partial differential equations, particularly in the early stages of a design study when many parameter variations need to be examined quickly.

[page 18: figure — additional three-dimensional surface-plot outputs demonstrating the visualization capabilities of the analog computer for 2D field problems, showing u(x,y,t) for various t values]

An iterative error-elimination method is applied by passing the result of each computational step back through the system. The greater the number of iterations, the more accurately the individual differential equations can be solved. A single integration step is sufficient for simpler differential equations; for more demanding problems, several passes may be required.

Through iterative computation, it is possible to eliminate errors that arise during the integration of differential equations. An integration step can be executed multiple times, and on each pass the accumulated error is reduced further.

Various iterative approaches to solving differential equations are available on the analog computer; the most commonly used are the iterative methods associated with the functional relationships of the solution variables.

Figures 8, 19

[Page 19: The page consists primarily of a triangular/pyramidal diagram (Fig. 8, 19) illustrating the iterative error-elimination scheme, accompanied by a references list at the bottom of the page.]

References

[1] C. L. Johnson: Analog Computer Techniques. McGraw-Hill Book Co., New York 1963.

[2] M. Korn and T. M. Korn: Electronic Analog Computers. McGraw-Hill Book Co., New York 1956.

[3] G. W. Soroka and R. C. Wood: Principles of Analog Computation. McGraw-Hill Book Co., New York 1954.

[4] A. S. Jackson: Analog Computation. McGraw-Hill Book Co., New York 1960.

[5] H. Herschel: Analogie-Rechenmaschinen (Analogy Computing Machines), in “Taschenbuch der Nachrichtenverarbeitung” (Pocket Reference for Signal Processing). Springer-Verlag 1962.

[6] G. Doetsch: Introduction to the Theory and Application of the Laplace Transformation (Einführung in Theorie und Anwendung der Laplace-Transformation). Birkhäuser Verlag, Basel 1958.

[7] The Solution of Partial Differential Equations by Differential Analyzers (anonymous). McGraw-Hill Book Co., New York 1953.

[8] K. Brunner-Schwer and H. Yanagawa: A Study of Simultaneous Partial Differential Equations Using the Analog Computer. Internal report, 1965.

[9] A. Ralston, H. S. Wilf, and G. Levy: Mathematical Methods for Digital Computers (Mathematische Methoden für Digitalrechner). R. Oldenbourg Verlag, Munich 1965.

(Submissions received on 31. 5. 1958)