IBM_Scientific_Computation_Forum_1948 IBM Scientific Computation Forum 1948
IBM_Scientific_Computation_Forum_1948 IBM_Scientific_Computation_Forum_1948
User Manual: IBM_Scientific_Computation_Forum_1948
Open the PDF directly: View PDF 
.
Page Count: 127
| Download | |
| Open PDF In Browser | View PDF |
PROCEEDINGS
Scientific
Com putation
Forum
1948
H. R. J.
GROSCH, Editor
WATSON SCIENTIFIC COMPUTING LABORATORY
INTERNATIONAL
BUSINESS
NEW
YORK
+
MACHINES
NEW
YORK
CORPORATION
Copyright 1950
International Business Machines Corporation
590 Madison Avenue, New York 22, N. Y.
Form 22-5904-0
P
R
I
N
TED
I
N
THE
U
NIT
E
D
S
TAT
E
S
o
F
AMERICA
FOREWORD
A
SCIENTIFIC COMPUTATION FORUM, sponsored by the
International Business Machines Corporation, was
held in the IBM Department of Education, Endicott,
New York, from August 23 to August 26, 1948. The
Forum concluded with two sessions held in New York
City on August 27.
Earlier meetings in this senes, which began in
1940, were devoted largely to statistical procedures. In
the 1948 Forum, for the first time, an attempt was made
to cover many of the fields in which large-scale computing methods have proved important. The exchange
of ideas between workers in fields as diverse as aerodynamics and physical chemistry proved fruitful from
the very beginning, yet specialists in the same field also
found time for intensive discussions.
It is hoped that the contributions printed here
will prove of value not only to the participants but to
other members of the growing group engaged in technical calculations on punched card equipment.
CONTENTS
Evaluation of Higher Order Differences on the Type 602
Calculating Punch
-FRANK M. VERZUH
9
Differencing on tbe Type 405 Accounting Machine
-GERTRUDE BLANCH
14
The Use of Optimum Interval Mathematical Tables
-H. R. J. GROSCH
23
Punched Card Techniques for the Solution of Simultaneous
Equations and Other Matrix Operations
-WILLIAM
28
D. BELL
Two Numerical Methods of Integration Using Predetermined Factors
Integration of Second Order Linear Differential
Equations on the Type 602 Calculating Punch
Integration of the Differential Equation
Type 601 Multiplying Punch
~~
-LELAND W. SPRINKLE
32
- N . ARNE LINDBERGER
34
=P·F(r) Using the
r
-PAUL HERGET .
Some Elementary Machine Problems in the Sampling
Work of the Census
- A . ROSS ECKLER
IBM Applications in Industrial Statistics
-CUTHBERT
Some Engineering Applications of IBM Equipment at the
General Electric Company
-FRANK
Planning Engineering Calculations for IBM Equipment
c.
.
.
.
.
.
42
45
HURD
49
52
J. MAGINNISS
-BEN FERBER
.
39
.
.
.
A Survey of the IBM Project at Beech Aircraft Corporation
Aerodynamic Lattice Calculations Using Punched Cards
-JOHN KINTAS
54
-HANS KRAFT.
60
. Dynamics of Elliptical Galaxies
67
-JACK BELZER, GEORGE GAMOW, GEOFFREY KELLER
Application of Punched Cards in Physical Chemistry
-GILBERT
70
\v. KING.
Application of Punched Card Methods to the Computation of Thermodynamic Properties
of Gases from Spectra -LYDIA G. SAVEDOFF, JACK BELZER, H. L. JOHNSTON
74
Calculation of the Equilibrium Composition of Systems of
Many Constituents -STUART R. BRINKLEY, JR., ROBERT
77
Punched Card Calculating and Printing Methods in the
Nautical Almanac Office
W. SMITH, JR.
-FREDERICK H. HOLLANDER
Programming and Using the Type 603-405 Combination Machine in the
Solution of Differential Equations
-GEORGE S.
Applications of Punched Card Equipment at the
Naval Proving Ground
..
FENN
.
.
.
.
95
99
-CLINTON C. BRAMBLE
Use of the IBM Relay Calculators for Technical Calculations at
Aberdeen Proving Ground
-JOSEPH
83
.
101
Simultaneous Linear Equations
-FRANCIS J. MURRAY
105
Computation of Shock Wave Refraction on the Selective
Sequence Electronic Calculator
-HARRY POLACHEK
107
Computation of Statistical Fields for Atoms and Ions
- L . H. THOMAS.
H. LEVIN.
.
•
123
PARTICIPANTS
ABRAMOWITZ, MILTON,
Mathematician
Computation Laboratory, National Bureau of Standards
New York, New York
ANDREWS, THOMAS B., JR.,
Aeronautical Research Scientist
National Advisory Committee for Aeronautics,
Langley Aeronautical Laboratory, Langley Air Force Base, Virginia
ANDREWS, T. G.,
Scientific Advisor
General Staff, Department of the Army, Washington, D. C.
ARNOLD, KENNETH
J.,
Assistant Professor of Mathematics
University of Wisconsin, Madison, Wisconsin
Project Engineer
ARNOLDI, W ALTER,
Vibration Section, Hamilton Standard D.ivision,
United Aircraft Corporation, East Hartford, Connecticut
BAILEY, EDWARD W., Statistician
Carbide & Carbon Chemicals Corporation, Plant Y 12,
DUNWELL, S. W.
Future Demands Department
International Business Machines Corporation, New York, New York
ECKERT, W. J.,
Director
Department of Pure Science,
International Business Machines Corporation, New York, New York
Deputy Director
ECKLER, A. Ross,
U. S. Bureaq of the Census, Washington, D. C.
FENN, GEORGE S.,
Research Engineer
Northrop Aircraft Incorporated, Hawthorne, California
FERBER, BEN,
Research Engineer
Consolidated Vultee Aircraft Corporation,
San Diego, California
FRANTZ, WALTER,
Engineer
Boeing Airplane Company, Seattle, Washington
Oak Ridge, Tennessee
GEAN, JAMES,
BELL, WILLIAM
D., Vice-President
Telecomputing Corporation, Burbank, California
BELZER, JACK,
Supervisor
GROSCH, H. R. J.,
Computing Section, Cryogenic Laboratory,
Ohio State University, Columbus, Ohio
BISCH, PAUL E.,
Engineer
Chief of Structures
Beech Aircraft Corporation, Wichita, Kansas
BLANCH, GERTRUDE,
Head
Department of Research in Physical Chemistry,
Mellon Institute, Pittsburgh, Pennsylvania
BRAMBLE, C. C.,
Director
HARDGRAVE, E.
J., JR.,
Research Group Engineer
Analysis Section, Ordnance Aerophysics Laboratory,
Daingerfield, Texas
BRINKLEY, STUART R., JR.,
HARMAN, H. H.
Personnel Research Section, Office of the Adjutant General,
Department of the Army, Washington, D. C.
HF;RGET, PAUL,
Director
University of Cincinnati Observatory, Cincinnati, Ohio
HOLLANDER, FREDERICK H.
Computation and Ballistics,
U. S. Naval Proving Ground, Dahlgren, Virginia
Physical Chemist
Explosives Branch, U. S. Bureau of Mines,
Pittsburgh, Pennsylvania
CALDWELL, S. H.,
Staff Member
Los Alamos Scientific Laboratory,
Los Alamos, New Mexico
Mathematician
Institute for Numerical Analysis,
National Bureau of Standards, Los Angeles, California
BOWMAN,JOHN R.,
Senior Staff Member
Watson Scientific Computing Laboratory,
International Business Machines Corporation, New York, New Yark
HAMMER, PRESTON C.,
Special Projects, North American Aviation, Inc.,
Los Angeles, California
BLAKELY, W. R.,
Chief Analytical Engineer
Aircraft Division, Parsons Corporation,
Traverse City, Michigan
.
Professor of Electrical Engineering
Massachusetts Institute of Technology,
Cambridge, Massachusetts
CHEYDLEUR, B. F.
Numerical Analysis Section, U. S. Naval Ordnance Laboratory,
Washington, D. C.
COSTELLO, GEORGE R.
Flight Propulsion Research Laboratory,
National Advisory Council for Aeronautics, Cleveland, Ohio
Nautical Almanac Office, U. S. Naval Observatory,
Washington, D. C.
HOPKINS, RALPH,
Special Representative
International Business Machines Corporation,
Washington, D. C.
HURD, CUTHBERT C.,
Technical Research Head
Mathematical Statistics and Automatic Computing
Carbide & Carbon Chemicals Corporation,
Oak Ridge,. Tennessee
JAKKULA, A. A.,
Executive Director
Texas A. and M. Research Foundation,
College Station, Texas
JOHNSTON, JOHN,
Staff Member
Los Alamos Scientific Laboratory,
Los Alamos, New Mexico
KING, GILBERT W.,
Physical Chemist
Arthur D. Little, Inc.,
Cambridge, Massachusetts
Department of Pure Science,
International Business Machines Corporation, New York, New York
Structures Research Engineer
KINTAS, JOHN,
Beech Aircraft Corporation, Wichita, Kansas
Assistant Head
KORN, GRANINO,
Research Section, Curtiss-Wright Corporation,
Columbus, Ohio
KRAFT, HANS,
Aerodynamicist
Turbine Engineering Divisions, General Electric Company,
Schenectady, New York
KROLL, C. W.
KUSNER, JOSEPH, LT. COL.
Military Requirements Division, Munitions Board,
Office of Secretary of Defense, Washington, D. C.
Chief
LEVIN, JOSEPH H.,
Machines Branch, Ballistics Research Laboratories,
Aberdeen Proving Ground, Maryland
IBM Supervisor
LIGON, LEONARD,
Analytical Group, Aircraft Division, Parsons Corporation,
Traverse City, Michigan
LINDBERGER, N. ARNE,
Fellow
Royal Institute of Technology,
Stockholm, Sweden
Senior Engineer
Research and Development Laboratory,
International Business Machines Corporation, Poughkeepsie, New York
MACK, CHARLES
E., JR., Chief of Research
Grumman Aircraft Engineering Corporation,
Bethpage, New York
MAGINNISS, FRANK
J.
MARSH, H. W.,
Staff Member
U. S. Navy Underwater Sound Laboratory,
New London, Connecticut
C.
F.,
Staff Research Engineer
Scintilla Division, Bendix Aviation Corporation,
Sidney, New York
McNAMARA, WALTER
J.,
Institutional Representative
International Business Machines Corporation,
New York, New York
MORRISSON, REEVES,
Head
Analysis Section, Research Department,
United Aircraft Corporation, East Hartford, Connecticut
MURRAY, FRANCIS
J.,
Professor of Mathematics
Columbia University, New York, New York
MYERS, FRANKLIN G.,
Design Specialist
Glenn L. Martin Company, Baltimore, Maryland
POLACHEK, HARRY,
Mathematician
U. S. Naval Ordnance Laboratory,
White Oak, Maryland
Mathematician
Machine Development Laboratory, Applied Mathematics Division,
National Bureau of Standards, Washington, D. C.
D., Mathematician
RINALDI, LEONARD
Wind Tunnel Department, Cornell Aeronautical Laboratory, Inc.,
Buffalo, New York
SAVEDOFF, LYDIA G.,
Associate Supervisor
Cryogenic Laboratory, Ohio State University,
Columbus, Ohio
E.
Special Representative
C.,
International Business Machines Corporation,
New York, New York
Mathematician
SCHULTZ, OSCAR T.,
Sperry Gyroscope Company,
Great Neck, New York
Senior Staff Member
SEEBER, R. R., JR.,
Department of Pure Science,
International Business Machines Corporation, New York, New York
SELVIDGE, HARNER,
Director
Special Products Development, Bendix Aviation Corporation,
Detroit, Michigan
SIMMONS, LANSING G., Chief
U. S. Coast & Geodetic Survey,
Mathematician
Washington, D. C.
SPRINKLE, LELAND W.
Acoustics Division, Naval Ordnance Laboratory,
White Oak, Maryland
STANLEY,
J.
P.
Computation Centre, McLennan Laboratory,
University of Toronto, Toronto, Ontario, Canada
C.
STEVENSON,
Analytical Division, General Electric Company,
Schenectady, New York
MAYLOTT,
RHODES, IDA,
SCHROEDEL,
Computing Division, Army Map Service, Department of the Army,
Washington, D. C.
LUHN, H. P.,
REID, WILLIAM H.
H.,
Chief
Strength and Weight Section, Douglas Aircraft Company,
El Segundo, California
THOMAS,
L.
H.,
Senior Staff Member
Watson Scientific Computing Laboratory,
International Business Machines Corporation, New York, New York
THOMPSON, OSCAR N.,
Structures Engineer
Product Engineering Department,
Consolidated Vultee Aircraft Corporation, Fort Worth, Texas
E., Mathematician
TILLITT, HARLEY
U. S. Naval Ordnance Testing Station,
China Lake, California
TUCKER,LEDYARD,
Head
Department of Statistical Analysis,
Educational Testing Service, Princeton, New Jersey
TUKEY, JOHN W.,
Associate Professor of Mathematics
Princeton University, Princeton, New Jersey
VERZUH, FRANK
M.
Electrical Engineering Department,
Massachusetts Institute of Technology, Cambridge, Massachusetts
WILLIAMS,
J.
V., JR., Technical Engineer
International Business Machines Corporation,
Poughkeepsie, New York
Evaluation of Higher Order Differences on the
Type 602 Calculating Punch
FRANK
M.
VERZUI-I
Massachusetts Institute of Technology
SO M E of the uses of the higher differences of a given
function are:
1. Numerical integration using finite differences,
2. Numerical differentiation using finite differences,
3. Subtabulation or interpolation to a fixed interval,
4. Location of errors in a given set of data,
5. Smoothing the irregularities in experimentally obtained data.
ence by an amount E times the binomial coefficients of
(a - b )n. An error of + 1 affects five fourth differences
by amounts of + 1, -4, +6, -4, + 1. An error of +.1
affects three second differences by + 1, - 2, + 1. This
characteristic variation of magnitude and sign serves to
locate errors immediately.
This method of evaluating higher differences on the
The tabulation of a function and its higher differences
shown in Figure 1 illustrates the derivation of the equations used'for the computation of the fourth differences.
Mathematical Basis of This Differenc,ing 1\11ethod
p02 was originally utilized to detect the presence of errors
in a computed table of the function x tanh x.
In this circuit, one card containing the value of the function is used for each given value of the argument. As
each function card passes through the 602, it has' the
higher differences as ~ell as the recomputed function
punched upon it. Since the values of the function and its
higher differences at preceding ordinates are required in
each computation, the cards must be arranged in their
proper sequence. At the beginning of each set of cards,
the higher order differences are not directly available. We
merely assume that they are equal to zero and therefore
the first few values ot these differences will be in error.
The recomputed function provides a complete check
upon the entire operation of the Type 602 Calculating
Punch. The recomputed function will be identical with the
given function only if all the cross footing operations required to solve the set of algebraic equations are correct.
This agreement is definite proof that all the calculations
are correct.
The comparison circuits of the accounting machine are
used to indicate any discrepancy between the given function and the recomputed function. An asterisk is printed
adjacent to the error each time a discrepancy occurs. Any
machine errors in the differencing process are immediately
detected by visually scanning the printed record for the
presence of asterisks.
The presence of errors in the given function are also
quite evident. An error of magnitude E in the function
will affect (n
1) consecutive values of the nth differ-
FUllction
First
Differences
Third
Second
Fottrth
!-2.0
6,'-1.5
!-1.0
6,i_ 0. 5
fo.o
tl O. 5
f1.0
fl1.5
f2.0
l1 i1.0
6,i~.o
l1~.o
fl~.o
l1ii-0.5
l1ii0.5
tlv0.0
i
6, y'o
flii
1.5
6,i2.5
f3.0
FIGURE
1.
TABULATION OF A FUNCTION AND
ITS I-IIGHER-ORDER DIFFERENCES
The first difference may be defined, for instance, by
6,i2..., = f3.0 - f2.0 ;
(la)
the second difference by
(lb)
.6i~.O = 6,i2.5 - fl1.5
the third difference by
(Ic)
.6'1~5 = .6to - .6i~.o
and the fourth difference by
(Id)
Higher order differences may be expressed by equations
of several different forms which are some combination of
+
9
10
SCI E N TI FIe
Figure 1 and equations (1) and (2) have illustrated the
origin and derivation of the formulas which are necesspry
in the computation of the fourth order differences. However, the principle of this differencing technique may be
obtained by a consideration of the simpler equations used
for the computation of second order differences which are
indicated below. The succeeding description will be devoted to a discussion of these second order difference
equations:
the formulas given in equations (1). The derivation of an
equation for the fourth difference may be given in terms
of the function at two given points and the backward
diagonal differences from the first point:
13.0 - (/2.0 + 6 i1 . 5)
i
6~.5 = 6 ro + 6\.5
13.0 = 6~.5 + 12.0
6i~.0 =
In addition to the basic equation given for the fourth
difference, the corresponding equations for the differences
of lower order may be given:
6 ii.o + 61~~5
6i~.0 = 6i~~5 + 6 i1.o
6 i2 . 5 = 6 iJ.o + 6\.5
13.0 = 6 i 2 . 5 + 12.0 .
L~H~5 =
(2b)
(2e)
Discussion
~ultiplier
Prob. No.
F VA /
11.4 T/()N
TYPE
?tJ2
Multiplicand
1 - 2
en.. J
I
I
I
I
~
+--- ~
/lDo
~
RO-RC
!::/.'Z
----
$U81:'
CrL.Z
/J'2S
-=--
PROG.
CTL.3
TRAIIISFER
TO.
STOR.llG£
tf'O-RC
-
Summary
Resul t Storage
13 - 16 14 - 15
~
RC
RC
I
~
_I
I
I/I)D
I
f1"Z
./
I
I
~r--t-RO
~
;
I
I
--+
/100
.1' Z.5
-;:-----
---1--- ~ r-...+
I
I
I
£,3
- --r--
I
I
t-RO-
2.
THE
I
I
SUST.
FIGURE
lAI/TH
I
/1'
PR06.
4../II)/4A
I
/
:I
13
CYCL.E
Ptf'OG.
{)RrJFJ? D/EF£RFNr;.-.<,
RHC
LHC
ADO-SlleT.
/rEADING
(3c)
r; PL rIlL ArnR
[f2.+~/.5] /1' 1.5
RC
RC
MJlST£/?
MULTIPL.Y
2@
9 - 12 10 - 11 5 - 6 - 7
3 - 4
X 60
Date
Ed. No.
()F
(3b)
The time-sequence chart shown in. Figure 2 illustrates
the flow of information in.the Type 602 during a computation of the second order differences. Assuming that the
value of the function and its first difference at the
preceding argument are available, i.e., if 12.0, 6\.5 and
- ({2.0 +6\.5) are in the. indicated counters, the second
difference is computed as follows:
(2d)
~~ZtLt:f.
(3a)
Operation of the Differenci,l'l.g Circuit
(2c)
These equations are so arranged as to utilize the results
of the preceding equation in the computation of the succeeding equation.
Customer
COM PUT A T ION
I
FLOW CHART OF
RO
---602
liDO
£.3
---
RO-RC
=-=------------ :::::----
COMPUTATION
f.t...
FORUM
11
PROCEEDINGS
CALCULATING PUNCH-TYPE 602-CONTROL PANEL
234567.
--------5--C~O~N~~~~~~
23
24
TORS
25
26
27
C
28
29
30
31
32
33
5------CARD
0000000000
1
25
0000000000
-.. . . . . . . . -·R
40
3.
42
43
~
35
65
75
44
20
40
o
0
0
0
o
0
0
0
45
0000000,00
FIGURE
41
1 _
£
OA
6O~
ON
10 G
o I
WIRING DIAGRAM FOR EVALUATION OF SECOND-ORDER DIFFERENCES
1. The value of 13,0 is added to - (/2,0 + 6\,5) in
counters (9, 12) to provide 6~.0 in accordance with
equation (3a).
2. This value of 6~.0 is stored in the MP counters
(1, 2) and simultaneously added to 6\.5 in counters
(10, 11) to provide the new first difference 6 12 . 5 111
accordance with equation (3b).
3. The value of 6 12 . 5 is subtracted into counters (9, 12)
where it will be available for the computation of the
next card. In addition, 6 12 . 5 is added to 12.0 in
counters (13-16) to provide the recomputed value
of 13.0 as indicated in equation (3c).
4. The value of 13.0 is subtracted from - 6\!.5 in
counters (9, 12) to provide - (/3.0
6 i2 . 5 ) . This
+.
quantity is required for the computation of the next
card; 13.0 is stored in summary counters (14, 15).
5. During the transfer-to-storage cycle the values of
6 i to, 6 i 2 .'J and the recomputed /3. 0 are punched on
the card containing the original function 13.0'
Since this computation is a sequential operation, i.e., the
terminal values of the computation of the first card become
the initial values for the computation of the next, counters
(9, 12), (10, 11), and (13 -16) are not reset during the
card feed cycle. A master X60 card is used to clear the
machine after the differencing of an entire set has been
completed. The actual 602 wiring diagram used in the
computation of the second order differences is shown in
Figure 3.
11
SCIENTIFIC
SERIAL
1
689
t:!.'
1Ja'
f.
#"
~/l
COMPUTATION
i'
~v
6245
689
6245
689
6245
689
6245
689
6245
689
6245
0143
-683
6102
.... 373
2347
-062
8592
695
6388
683
6114
56
8461
701
6543
89
-683
6025
707
6799
0269
-000
0358
713
6887
538
719
7076
2
695
6388
6
3
701
6543
6
0155
12
0256
101
-4
707
709
6
5
713
6887
6
0088
168
0189
101
-000
269
6
719
7076
6
7
125
7278
6
0202
13
-000
0088
0357
725
7278
0213
11
-000
0002
86
731
7491
2
737
7715
1
1
743
7951
8
731
7491
6
9
737
7715
6
0224
11
743
7951
6
0236
12
749
8199
6
1260
0
1024
12
755
8459
5
9248
0
2012
13
761
8736
6
0277
0
1029
14
767
9013
6
0277'
15
773
9308
6
0295
18
16
779
9614
6
0306
11
17
785
9932
6
0318
12
18
792
0262
6
0330
12
0341
11
10
___ 11
19
798
0603
FIGURE
4.
6
-000
1012
1011
749
9211
-000
3036
-000.4048
755
8459
6077
761
8736
-000
1029
4070
767
9013
-000
0007
3041
-000
18
1047
773
9308
-000
00.25
779
9614
8
785
9932
-000
0001
792
0262
-000
0001
798
0603
1
-000
0001
ResuLTs OBTAINED DURING EVALUATION OF FOURTH-ORDER DIFFERENCES
The accounting machine is used to provide the printed
record shown in Figure 4. This record illustrates the information present on each punched card after it leaves the
602. This record contains a serial number, the given
function, the higher order differences in increasing order,
and the recomputed function. The two errors present in
the functional data are quite evident because of the characteristic appearance of the higher differences.
The flow chart diagram shown in Figure 5 illustrates
the flow of information within the 602 during the evaluation of the fourth order differences. The solution of the
set of simultaneous equations (2) is effected by theindicated crossfooting operations.
Lim#ations and Advantages of This Differencing Method
The Type 602 Calculating Punch has been wired to
obtain fourth order differences as well as second order
differences. Whenever the given function is quite irregular, the higher differences are of appreciable magnitude.
In such cases, the available twenty-four positions in the
result storage counter are not adequate to punch all of the
computed differences. Occasionally, the punching of the
first and second differences is omitted, and the available
result storage positions are used to punch the large-size
higher order differences.
The second order difference board is used for many
applications. Whenever the second differences are too
large for proper interpretation, the cards are reinserted
and the operation is repeated. In such a case, the second
order differences are used as an input function, and the
fourth order differences are then obtained.
This method of differencing is quite rapid, as differences
may be obtained at a rate of 1500 an hour. This number
varies slightly with the number of columns punched.
In addition to the error detection properties of this
method, it may be easily used for subtabulation, i.e., interpolation to fixed intervals. The 602 may be wired to
perform integration using the given value of the function
and its backward differences.
Since this method does not require any blank cards for
operation, it is capable of a higher operating speed than
the methods which do use blank cards.
DISCUSSION
[This paper and the following one by Dr. Gertrude Blanch were
discussed as a unit.]
lL.E t£Z.UI::l..
Customer
Prob. No.
Discussion
EVALIJATI{)N
TYPE
Multiplier
Multiplicand
1 - 2
3 - 4
X-GO
4
OF
r
~n2
.4 L
LRC
fJ."1
0.5
T#
r: {II
9 - 12 10 - 11
[fZ. +11'1.5+/J.'+/1'
ID 0.5
MASTER
Ed. No.
()IFF£R£NC:~.")
ORDER
ATJ Nc,
Date
Summary
5 - 6 - 7
13 - 16 14 - 15
I
I
1.0
I
/
I
7i-t~
PIINr:1-I
RRC
iii
lNlTI-I
4/J.5/4-R
Result Storage
~'1.5
£2
RC
RC
I
RC
RC
I
RC
I
i
I
I
/lDD -JUST
RERD
DETAIL
.... [$
fJ.'V,..
I
.,
I
I
.~
PROG.
CTL..1
P~OG-.
CTL.2
I
I
~ --..+
/lOD
----
I
.6'V
1.0
L1'V
1.0
RO-RC
I
ADD
I
~
~
~
I
~+
RO
11'"0.5
C-rL.3
-=---===
I
ADD
SUST.
I
Ll.'"0.5
~
LJ.~O
PROG.
I
I
~
r--
I
L
-r--t----r--.-+
I
I
SUBT.
RO
/1"2.0
I
I
I
IIDD
£J.II
2.0
~
~z....5
...............
I
FROG.
OUST
CTL.4
.1'2.5
I
I
I
I
~+
/IDD
RO
Ak.5
-;:-----'
.!;,3
I
I
PI('OG.
en. 5
I
SUBT.
-
£3
I
I
--=
I
~
-
RO
--..
~
~ ~::: ~- t-=--r---I
TRIINSFER
To
5TO~IIG£
I
RO-RC
RO
I
I
RO
flv
6/1'
I
I
I
I
I
FIGURE
S.
EVALUATION OF FOURTH-ORDER DIFFERENCES WITH THE
13
602
--
f3
Differencing on the Type 405
Accounting Machine
GER TR UDE
BLANCH
Institute for Numerical Analysis, National Bureau of Standards
D IFF ERE N C E S have to be taken frequently, not
only for the purposes of integrating and differentiating
but for the purposes of checking data in the process of
computation. I f functions are given at uniformly spaced
intervals, the process of differencing at strategic stages of
the computing process offers a very satisfactory check on
the operations.
N ow, since differencing has to be done so frequently, it
is important to be able to do it on many machines.
Sometimes the 602 is tied up on other work and you
want to do it on the accounting machine. The 405 was
used for differencing long before any other IBM machine.
Everybody knows how to take a first difference, summary
punch it, and then take the second difference. In this
manner successive differences of any order can be built
up. But you can also take a sixth difference or fifth difference or fourth difference without taking intermediate
differences. This is valuable if you want the differences
mainly for checking data. For when you take a high
enough difference-the sixth difference for instance, if it
is small-you can generally detect the errors by the difference pattern and actually takeout the card where the
error occurs.
In the summer of 1947, Dr. E. C. Yowell spent his time
in our New York laboratory. Dr. Abramowitz had picked
up Comrie's paper on getting higher order differences on
the National bookkeeping machine, and he asked Dr.
Yowell to do something similar on the accounting machine.
Dr. Yowell wired a control panel which has been used
very successfully, and I will try to give you the wiring of
that panel.
(Since the rest of the talk depended heavily on blackboard diagrams and slides, r have taken the liberty to substitute in its place Dr. Yowell's own lucid description of
the wiring. It goes into greater detail than my own talk
did, and will be much easier to follow for anyone who
wants to reproduce the wiring.)
Sixth Differences on the 405
EVERETT C. YOWELL
Institute for Numerical Analysis,
National Bureau of Standards
THE MET HOD used. in computing the sixth difference is that given by Comrie. 1 The first six functions are
used to compute the first through the fifth difference.
There is then available in the machine
The sum of these six quantities is an approximation to f7 ,
and the difference between f7 and the sum of these terms
will be 6.64 • This can be verified by writing
/7 = /6 +
1
6. 1312
Hence
/7
=
/6 + 6i1/2 +
6.~ ,
Hence
This process can be continued until the fifth difference is
written as the sum of the previous fifth difference and a
sixth difference. Then transposing the equation gives
6 64 = f7 - (f6
+ 6.~112 + £S5 + !i~/2 + 6.~ + 6.57/2)
.
The machine process is as follows: we assign one
counter to the function and one to each difference-seven
counters in all. Let us suppose the machine is set up with
6.~ in Counter 1, 6.~/2 in Counter 2, 6.44 in Counter 3,
6.39/2 in Counter 4, 6.~ in Counter 5, 6.\1/2 in Counter 6
and f7 in Counter 7. In the next eight card cycles we will
compute 6.65 ,
14 ,
FORUM
15
PROCEEDINGS
During the first card cycle, 6.64 is rolled out of Counter 1
into Counter 2. This addition of 6.~ t.o 6.~/2 gives 6.~/2 in
Counter 2. At the end of this cycle, a total is taken, printing
6.64 out of Counter 1 and resetting this counter. During the
second card cycle, 6,59 / 2 is rolled out of Counter 2 and
entered positively int.o Counter 3 and negatively into
Counter 1. This addition of 6.~/2 to 6.44 gives 6.45 in
Counter 3, while - 6.5912 stands in Counter 1. During the
third card cycle, 6.~ is rolled out of Counter 3 and entered
positively into Counter 4 and negatively into Counter 1.
This addition of 6.45 to 6.~12 gives t111/2 in Counter 4,
while - (6.5912 + 6.\) stands in Counter 1. During the
fourth card cycle, t111/2 is r.olled out of Counter 4 and
entered positively into Counter 5 and negatively into
Counter 1. This addition .of t11112 to 6.; gives 6.26 in
Counter 5, while - (6. 5912 + 6.45 + 6.311 12) stands in
Counter 1. During the fifth card cycle, 6.26 is rolled out
of Counter 5 and entered positively into Counter 6 and negatively int.o Counter 1. This addition of 6.26 to 6.\1/2 gives
6.113 / 2 in Counter 6, while - (6.~12 + 6.45 + 6.31i/2 + 6.26 )
stands in Counter 1.
During the sixth card cycle, only the negative transfer
into Counter 1 is needed. In the previous cases, we have
had to build up our difference order of n from a higher
order difference and the previ.ous order of difference n.
But the function has been read from the card, so that the
function does not have to be built up from the previolls
function and first difference. Hence, during the sixth card
cycle, 6.113 /2 is rolled out of Counter 6 and entered negatively into Counter 1, thus giving - (6.59 / 2 + 6.~ + 6.31112
+ 6.2 + 6.\312) in Counter 1. During the seventh card
cycle, i1 is rolled out of Counter 7 and entered negatively
into Counter 1. This gives - (6.59 / 2 + 6.~ + 6.311 / 2 + 6.26
+ 6.i312 + i1) in Counter 1. During this cycle, i1 is reset
as it rDlls, leaving the counter .open for receiving the next
function. This process will be explained later.
During the eighth cycle, is is read positively intD Counter
1 and Counter 7. This leaves is in Counter 7 for use in computing the next difference and completes 6.65 in C.ounter 1,
fDr 6.65 = is - (6.59 / 2 + 6.~ + 6.\112 + 6.26 + 6.\312 + i1)·
During the eight-cycle operation, each counter has advanced from one difference to the succeeding difference .of
the same order. Hence we are ready to repeat the cycle
again and compute 6.66. Since .one function card is read
every eight cycles, seven blank cards must f.ollow every
function card.
The first thing to be done is to set up an eight-cycle
control panel. This is done by lacing together seven X
distributors (Figure 1).
($
As the first card is fed into the machine, the VCI hubs
emit impulses at every digit position. The X impulse is
selected by the digit selector and passed on t.o the common
SI!:LECTIVE LIST CTRL.
~
~
X
0
NX C
010
FIRST CARD'CONTROL
o-.!!...o
o-!!!!-o
---43
o
0
0
~
0
SELI!:CTORS
NX
C
X
0 10
0
NX C
01 0
lC
X
0'1 o---Q
0
0 2 0
0
020
0
020
020--0
o-!Lo
0
030
0
030
0
030
030---0
0
040
HMR.LCK.
o!Q1O
0
040
0
040
000
oJLo
o
050
0
050
o
0
--43
0
45
..
:: 6
~:
'-I
:: S :,
l-CTRL
o
0
up- Sp
"
0 7 ~:
0
~:
o
P
:: 8 ::
( : : : 8 ::
0
0
0
0
"
7 c;
()
~:.
1~11~I~su''''CT.oT
CO.P.. lltING "!L.'f'-~UG I'ttOM
.
:: 8 0--0
__~1
u"." MUSMU - 1 1 - - - - - - - ,
00000000000000
+
8.:;
lONE. tfl.rCTIOM (0." • I')
000
,
00000000000000000000
00000000000000000000
T""'I......tc::--_~_O_CT._~T_.L~_T.L~
o
0
o
0
0
43
0
r
0
0
0
c•• o
0
or-o--o-~~9--0
EEl
COONT
FlGURlt
1
hub .of Distributor 7. Since 7 is not energized, the impulse
comes out of the NX hub and into the common hub .of 6.
Since this is not energized, the impulse comes .out of the
NX hub of 6 and enters the common hub of 5. Since 5, 4,
3, 2 are all unenergized, the impulse finally goes from the
NX hub of 2 to the common hub of 1. This distributor is
unenergized, so the impulse comes out .of the NX hub and
picks up Distributor 1. Once picked up, the distributor
holds for one cycle.
As the second card is fed, the VCI hubs emit impulses
fDr all digits, and the digit selector passes the X impulse
along to the distributor chain. It passes along the lacing,
as the first impulse did, until it reaches Distributor 1.
Since this is still energized, the impulse is shunted t.o the
X hub .of the distributor, and from here to the pickup .of
Distributor 2. Once picked up, this als.o holds during the
next card cycle. As the third card is fed, the VCI hubs
emit impulses for all digits, and the digit selector iSDlates
the X impulse and passes it along tD the distributor chain.
It passes along the lacing until it reaches Distributor 2.
Since this is energized, the impulse is directed to the X
hub .of DistributDr 2, from whence it picks up DistributDr
3. Notice that the impulse never reaches Distribut.or 1,
16
SCIENTIFIC
COMPUTATION
AUXILIARY WIRING DIAGRAM
000
0
25
0
0
FOR X DISTRIBUTION AND COUNTER CONTROL
0
65
0000000
x-
COL. _ _
0000000
OSC-I
UCI
0-0
0-0
o-ll,o
X- COL. _ _
X-COL. _ _
0000000
x-
-----25---.
0000000
6B O
0
0
I0
0
0
0 60 0
0
0
0
0
0
0
COL. _ _
X- COL. _ _
X- COL. _ _
x-
------25-----
rrI I I rr
COL. _ _
)(- COL.-.oj
-----65---
! III 1 I 1
80
----co c,---
+
0
9
o
SUP
0
~
.SA~
0--0-=-0--0
0-=-0-0--0
-------+------.
I6B
CDUNTU £ N T R V - + - - - - - - - - - - - i
1r I I I I leAI I I I I I I ISBI 1 I I
oC~C,
0
9 SUP
o
0
~-'-I ~
9 SUP
0
o
g
s~
o
FIGURE
since its path is broken by Distributor 2. As the fourth
card is fed, the VCI hubs and the digit selector pass an
X impulse along to the distributor chain. The impulse
passes along the lacing until it reaches Distributor 3. Since
this is energized, the impulse is diverted to the X hub of
3, from whence it picks up Distributor 4.
In a similar manner, the fifth card picks up Distributor
5, the sixth card picks up Distributor 6, and the seventh
card picks up Distributor 7. As the eighth card is fed, the
VCI hubs emit impulses for all digits, and the digit selector passes the X impulse on t.o the common hub of
Selector 7. Since this selector is energized,! the impulse is
shunted to the X hub of this distributor. As this hub is not
wired, the impulse has no effect on the machine, and all
selectors are unenergized as the ninth card starts to feed.
This is the same condition that we had when the first card
started to feed. Thus we have wired a sequence of events
which repeats itself every eight cycles.
Seven counters are necessary for a sixth difference computation. This permits handling ten-digit numbers. We
shall assign Counter 4A-6A to the function, 4B-6B to the
sixth difference; 4C-6C to the fifth difference; 2A-8A to
the fourth difference; 2C-8C to the third difference; 2B-8B
to the second difference; and 2D-8D t.o the first difference:
All counters are balance coupled, that is, t~e CI from the
highest position is wired back into the C hub of the units
0
o
0
000
o
Q
000
NX
0
0
0
0
c
2
position and the "hot 9" is jackplugged to the SVP hub.
An analysis of the desired counter additions and subtractions shows that the following counters must be activated on the indicated cycles.
Roll
Negatively
Into the
6,6 Counter
Cycle
Roll
Positively
Onto
1
6,6
6,5
2
6,5
6,4
6,5
4C-6C Subt.
2A-SA Add (2)
4B:-6B Subt.
3
6,4
6,3
6,4
2A-8A Subt.
2C-SCAdd (3)
4B-6B Subt.
4
6,3
6,?
6,3
2C-SC Subt.
2B-SB Add (4)
4B-6B Subt.
5
6,2
6,1
6,2
2B-SB Subt.
2D-SD Add (5)
4B-6B Subt.
6,1
2D-SD Subt. (6)
Counters
4B-6B Subt. (1)
4C-6C Add
6
H.oll
7
Roll the function 4A-6A Subt. (7)
4B~6B Subt.
S
Read the new function into the
4A-6A Add
(8)
function counter and the 6, 6 counter 4B-6B Add
FORUM
PROCEEDINGS
One cycle of Plug to C impulses is sufficient to control
all but the 6 6 counter. Since any difference counter adds
as the next higher difference counter subtracts, only the
subtraction needs to be wired to the Plug to C. The addition can be taken care of by wiring 2A-SA add to 4C-6C
subtract; 2C-SC add to 2A-SA subtract; 2B-SB add to
2C-SC subtract; 2D-SD add to 2B-SB subtract. The subtraction impulses are generated as shown in Figure 2.
It is at this stage of the wiring that the differencing
cycles must be correlated with the card feed cycles. Since
we wish to read from the first card, this step must take
place as the first card is under the lower brushes. This is
the last step in the differencing cycle. As the first card
feeds, an X impulse picks up Distributor 1. As the first
card passes the upper brushes, the X impulse picks up
Distributor 2. Hence this distributor is energized as the
first card passes the lower brushes, and the Plug to C
impulse causing Counters 6A and 6B to read must come
from the X hub of Distributor 2. This fixes the end of the
differencing cycle, and the remaining Plug to C impulses
are wired in sequence from this point.
Selector F is used for algebraic sign control. Only the
reading of the function depends on the punched sign. Once
it is entered into the counters as a complement, if negative,
or a true figure, if positive, it will always be transferred
as a complement or a true figure. Another position of the
selector will be used to indicate the sign of the function in
the listing process.
The entry of digits into all but the function and sixth
difference counters is made by using the card cycle total
transfer device. We have already wired the counters to
subtract on the cycle when they transfer, and to add on the
cycle when they receive. The wiring of the counter exit
and counter entry circuits is as follows:
from
to
4B-6B
4C-6C
Counter total exits.
Counter entries.
from
to
4C-6C
2A-SA
Counter total exits.
Counter entries.
from
to
2A-SA
2C-Se
Counter total exits.
Counter entries.
from
to
2C-SC
2B-SB
Counter total exits.
Counter entries.
from
to
2B-8B
2D-SD
Counter total exits.
Counter entries.
The function counter receives impulses only from the
brushes:
from
to
4A-6A
Lower brushes.
Counter entries.
The sixth difference counter must receive impulses from
all difference counters, the function counter, and the
17
brushes. This demands the use of several selectors. Only
five positions will be shown in each selector, but each is to
handle the full ten digit field.
Selector A is picked up on the card reading cycle by an
X control shot from the pickup hub of Distributor 2 to
the X pickup hub of the selector. Selector B is picked up
by an X control shot from the pickup hub of Distributor 1
to the X pickup hub of the selector. Selector C is picked
up twice; once by an X control shot from the pickup hub
of Distributor 4 to the X pickup hub of the selector, and
the second time by an X control shot from the pickup hub
of Distributor 6 to the D pickup hub of the selector. Selector D is picked up twice; once by an X control shot
from the pickup hub of Distributor 5 to the X pickup hub
of the selector and once by an X control shot from the
pickup hub of Distributor 7 to the D pickup hub of the
selector. These impulses are taken from the pickup hubs
of the distributors. It must be remembered that the pickup
hubs of the selectors and distributors are double hubs so
that an impulse wired into one hub is automatically emitted
from the other. Hence the selectors and the corresponding
distributors will pick up at the same time.
Wiring from the counter list exits is done to avoid the
use of split wires from the counter total exit hubs. Whenever a counter is impulsed to add or subtract, the counter
list exits and the counter entry hubs are internally connected. Thus, whenever Counter SA is impulsed to add or
subtract, any number reaching the counter entry hubs will
also reach the X hubs of Selector C by ~ay of the SA
counter list exits.
The read-in to Counter 6B will proceed as follows: as
Distributor 2 picks up, the first card feeds under the lower
brushes. Selector A also picks up, so that the entries to
6B are connected by way of the double entry hubs of 6A
to the brushes. This reads f into the counter. As Distributor 2 picks up, Distributor S picks up also (Figure 2)
and the new f is added or subtracted into 6B according to
the sign of the function. On the next cycle, Distributor 3
is picked up, but no selector is energized. Hence the 6B
entry hubs are connected to the 8B total exit hubs. As only
the 6C counter is impulsed on this cycle, no impulses are
emitted from 8B into 6B. Also 6B subtracts, sO' that the
difference in 6B is transmitted out through the total exit
hubs intO' the entry hubs of 6C which adds the sixth difference onto the previous fifth difference. On the next
card cycle, Distributor 4 is picked up. Hence 6C subtracts
and SA adds, and 6B subtracts through the NX hub of
Distributor 8. Selector C is also picked up. TherefO're, any
information reaching the SA entry hubs is emitted from
the SA list exit hubs and transferred to the 6B entry hubs.
The SA entry hubs are connected to the 6C total exit hubs.
As 6C subtracts, the fifth difference is emitted from the
total exit hubs and added into the previous fourth differ-
18
SCIENTIFIC
.------5-- ...
-------20-
a~~~~~~----------40·
o
0
0
0000000000000000000
o
0
0
0 6B O
0
0
1
M.... &.,CAL "'1" •• "IT INfRy
0
00
0 2A
o
0
0
028 0
I
0
0
0
I
0
0
4A o
-,,-r
00000
9 SUP 9 SUP
9
o
0
0
0
0
0
5 - - - "UIOl.'eAL
0
0
0
0
0
0
" " . . . . LIlT UTAY
0
0
0
0
0
--------20-
000
-----.-::~----40-
CDUNT£II .isT Elf" - - - - - - - - - - -
0
I
0
0
0
T
0
I
0
0
0 6A
o
0
0
I
0
0
~I~
00
SUP
9
0
0
c
c
0
SUP
0
o·
9
0
~~~
0-=-0
0-=-0
COMPUTATION
8A
+
28
14A
I
0
+
0-=-0---0-0
cr--o-=-o--<>
I2A II I I I I
SLf>
88~
I I I lSAI I I I I I I 1 l I I I
8B
c
"-'1'-'1'
~~
~
+
0-=-0
~
0
00000
9 SUP 9 SUP
o 000
9
0
0-=-0
0 69 0
0
60
0
SEL.4
o
0
o
0
0
o
0
000
0
0
Be
Sltt
80~
+
0-=-0-0-0
0
0
0
0
0
0
0
0
\i:,s
x
X
0
0
9
0
SEL. 6
SEL.5 SEL. 7
fff:.os EL E CTOR'C ."..ZF"'--O'->:..:()....",T.,..,OT".--.,.~
E LEe TO '0 ZF
-/1::;,
o
SUP
0
0
0
0
0
0
0
NX
NX
x
o
0
o
a
0
0
0
0
0
0
NX
C
0
10
0
0
o
000
------------20-
0000000
00000
000
0
0
40-
00000
o
0000000
0
0
0
000
0
0
0
0
0
0
0
0
20
0
000
00000000000000000
FIGUR~
ence in SA. It also passes through SA and is subtracted in
6B. On the next card cycle, Distributor 5 picks up. Hence
SA subtracts and SC adds and 6B again subtracts through
the NX hub of Distributor S. Counter SA transfers the
f.ourth diff~rence into the SC entry hubs where it is added
to the previous third difference.
The fourth difference also is emitted from the list exit
hubs and reaches the X hubs of Selector D. Since this distributor is energized, these digits reach the entry hubs of
6B and subtract on top of the fifth difference. On the next
cycle, Distributor 6 and Selector C are energized, SC subtracts, and SD adds; 6B still subtracts through the NX
hub of Distributor 8. Since se subtracts, it transmits the
third difference to SB, where it adds t.o the previous second
difference. The third difference is emitted from the SB list
exit hubs. From there, it passes through the X hubs of
Selector C into Counter 6B. Note that SA is not active in
0
40
0
3
this cycle. Hence the SA list exit hubs are dead, and no
confusion can result from double wiring the SB and SA
list exit hubs to the X hubs of Selector C. On the next
card cycle, Distributor 7 and Selector D are energized.
Counter SD adds, while SB and 6B subtract. Counter SB
emits the second difference, which enters SD and adds onto
the previous first difference. It is also emitted from the SD
list exit hubs and passes through Selector D into 6B.
Again the double wiring of the SC and 8D list exit hubs
t.o the X hubs of Selector D causes no confusion as 8e
and 8D are not active at the same time.
On the next card cycle, as no distributor or selector picks
up, SD subtracts and 6B subtracts (Figure 2). Thus SD emits
the first difference from its total exit hubs, and this passes
through the NX hubs of Selector D to Counter 6B. On the
next cycle, Distributor 1 and Selector A are energized and
FORUM
19
PROCEEDINGS
SELECTORS
X NX C
NX
0
C
0
0
0
0
0
0
0
0
0
X- OOL. _ _
0
0
0
0
0
0
0
0
0
0
X
)(-OOL. _ _
0
X- OOL. _ _
0
X-OOL. _ _
0
0
X- COL. _ _
0
0
0
0
0
0
X- OOL. _ _
0
0
0
0
0
0
X- OOL. _ _
0
0
0
()
0
0
X- OOL. _ _
o
0
0
0
0
0
PLUG TO"C'"
au
0---0--0--0---0
~
+
~
+
~
~
~
+
~
"1-
~
"1-
~ ~
+
+
~
0-<>-=-0--<>
0-0-=-0--<>
~
~
~
~
+
~
+
~
"1-
0--0-=-0--0
0---0-=-0----0
0-0-=-0----0
+
~
+
~
FIGURE
6A and 6B subtract. The function is emitted from the
total exit hubs of 6A and passes through the X hubs of
Selector A into Counter 6B. The next cycle feeds another
function card past the brushes and starts the operation all
over again.
All counters are to accumulate the numbers read into
them except the function and the sixth difference counters.
The sixth difference counter should reset after the difference has been rolled out. This is controlled as shown in
Figure 4.
The function counter is to be reset after it has transferred its value to 6B. This is best done by rolling the
function out of 6A back into 6A (Figure 5).
4
The pickup hub of Selector E is wired to the pickup hub
of Distributor 2. This is to prevent a reset when a negative function is read into 6A. The resetting principle is as
follows: 6A is subtracting in order to transfer the information. Hence all its counter wheels are turning. As each
wheel passes 9, it emits an impulse which enters the
counter through the list exit hub. When an impulse enters
a counter position while that counter is subtracting, it
stops the wheel. Consequently, the counter wheel is stopped
as soon as it emits an impulse-when it stands at 9. A
counter with 9 in every position is a counter containing
zero, if a "nines complement" system is used. Thus the
function counter is reset as it transfers~
f'·i'''r~I 11' aI 9I loI i zI
~~1~-~1~-9S~I~-95:1~
0000
00
00
~
0-=-0
0--<>--=-0---0
0---0-=-0--0
0--<>--=-0---0
1 111 111 luI III
24
28
0
9 SUP
00
III I 16A! 1111 r
9 SUP
00
~~~~
0..=-0
0-=..0
0---0-=-0---0
:
o
oNXo
Q
0
UUl'-------
FIGURE
5
0!2!0
x
000
oNXo
C
~-~l~-~r'-~I~~""oo--~ ~9 SUP 9 _ SUP
0000
x
000000
-------
14s1
ffrofIT1
~ELECTOR-FrJ!20
o!2Io o1Lo
x
~ s~p
~o¥o~~
0-=-0
1---::--~7o~LECTOR-E~
'1"-
o
~
C
0
0
Q
0
20
SCIENTIFIC
To print the sixth difference and list the argument and
function, the following wiring is used (Figure 6) :
COMPtJTATION
the comparing relays, and connecting the unequal impulse
to the minor control hubs.
When there is one intermediate and one final total
each cycle, the paper is spaced three times. In order to
bring all the listing and total printing for a single step
onto a single line, the upspace suppress is wired from first
card control intermediate. This wiring kills spacing on the
intermediate total cycle and allows ,only the single upspace
on the minor listing cycle.
In order to compute an nth difference, an (n + 2)-cycle
panel is needed and (n
1) blank cards must be inserted
between function cards. Higher differences than the sixth
can be computed on this same scheme if counter capacity
and distributor capacity are available. No extra selector
capacity is needed, as Selectors C and D can be multi wired
on the same scheme as indicated here.
If the function and argument are read into counters
instead of being listed, they can be summary punched
together with the highest order difference. An intermediate summary punch control will be sufficient. Since the
argument and function can be total printed at the same
time as the difference, the listing cycle is not necessary
and the minor control break can be eliminated. The intermediate control can then be shifted to minor control, saving one cycle for each difference.
The machine is cleared on the last card by wiring all
differencing counters to major total. Since a last card
causes all these total breaks, this will reset all counters
and clear the machine.
+
FIGURF;
6
Counters 4B-6B are wired for balance conversion (Figure 7):
o:r'
-;-~l"-~o~ o_---.;~_4_+_-_-~-,' ~!o-o
0---+--0
~
- ._
0 ._0---0-=-0--0
0--0-=-0---4
_'
"'"'''a'.'I1' _ _ _
_ _ __
0
1.1 I I I I
1
-
o ..o·
-~
• SUP
o
~
0
9
0
supi
0
:
0
'rl"-
o
-:'
..-
o 0
-~
0
9 SUP
_
,
0
FIGURF;
0
7
RF;I<'F;RF;N CF;
The argument and function are listed (Figure 8) :
1. L. J. CoMRU:, "On the Construction of Tables by Interpolation,"
Mon,thly Notices Royal Astroft. Soc., 88 (1928), pp. 447-59 and
518-22.
_ _ _ CONTROL 8RUSHES - - - - - - - - ,
o
0
0
0
o
0
700
0
_ _ .......*.. "'.'t" .. ''' ..'~.a.
0
0
0
:>
0
~
0
7~
0
0
4~1
o.
80
'
u.-----4O'''''l
ZOtl, .....
000000000000000
o
0
0
x
0
0
oNXo
0
C
OQOOOOO
FIGURE;
8
This will list an N after each negative function. If a
numerical type bar is used for the sign of the function,
read one Subtract Units Position into Selector F instead
of a "hot 5."
To get the listed numbers printed, a list cycle must be
introduced before the function card reads, This is done by
wiring the pickup of Distributor 2 to a control position of
DISCUSSION
Mr. Ferber: Do you insert blank cards?
Dr. Blanch: Yes, seven blank cards behind every function card, otherwise you couldn't do the operation. The
cards are put in by the collator.
};Ir. Bisch: Are they absolutely blank cards?
Dr. Blanch: Absolutely blank, and they come out blank.
You can use them over and over again. You are really
programming the accounting machine.
Dr. Abramo'lvitz: I would like to mention that the sixth
difference control panel described by Dr. Blanch may be
used for the computation of differences of lower order if
so desired. To illustrate this property it is instructive to
consider the following operations which take place to produce the successive differences. Let us consider a table of
values tv ,,~, ... in which we have interspersed seven blank
cards between successive function cards. I f the cards have
been fed into the machine and the 11 card is at the lower
brushes, we then compute as shown on page 21.
6"
6'
6'
63
6'
6'
Coullter
Coullter
Counter
Counter
Counter
Counter
I, card
+1,
Blank card 1
Clears
°+1,
°
°
°
°
Blank card 2
-I,
Blank card 3
-21,
Blank card 4
-31,
Blank card 5
-41,
Blank card 6
-51,
Blank card 7
I. card
Blank card 1
Blank card 2
- 1.+ 51,
-212+ 91,
Blank card 4
-31.+121,
Blank card 5
-41.+141,
Blank card 6
-51.+151,
Blank card 7
-61.+ 151,
I, card
Clears
Blank card 2
- 13+ 5/.- 101,
Blank card 3
-213+ 912-16f,
Blank card 4
-3f.+ 12f2- 19f,
Blank card 5
-4f.+ 14f2- 2Of,
Blank card 6
- 5f3+ 15f.- 2Of,
Blankt:ard 7
Blank card 2
- f.+ 513- 1Of.+1OI,
Blank card 3
-2f.+ 9f3-16f2+14I,
Blank card 4
--3f.+12f3-19f2+15f,
Blank card 5
-4f.+ 14f,-20f.+ 151,
Blank card 6
- Sf. + 15I,-20f,+ lSI,
Blank carel 1
Clears
- f5+ 5f.- 1Of.+lOf.-5f,
-21,+ 91.-1613+141.-61,
Blank card 4
-3f5+121.-1913+15f.-6f,
Blank card 5
-4f,+ 14f. -20f3+15f,-6f,
Blank card 6
-51,+151. -20f3+ 1512-6f,
Blank card 1
12- 21,
1.-1,=6',*
Clears
+1.
13- 5/2+ 10/,
13- 41.+61,
I, ·-3f.+3I,
f.- 2f.+f,=62,
il,+6', = 6'.
Clears
+1.
f.-5f.+1012- lOf,
f.-4f,+6f.-4f,
1.-3f,+3f.-f, =6',
IS, +
/'i" = 6'.
6'2+62.=t};.
Clears
+1.
f5- 5f.+1013- 1Of2+ 5f,
f,-41.+6f.- 4f2+I, = 6\
6"+6',=6'.
IS, + 6', =
6~,
6\+6'3=6~
-6f,+15f. -2°f, + 151.-6f,
Clears
- f.+ 5f.-lOf.+ 1Of3- 5f.+f,
Blank card 3
-21.+ 915-16f.+14f3-6f.+f,
Blank card 4
-3f.+ 12f,-19f.+15f3- 6f.+f,
Blank card 5
-41.+ 1415-20f.+ 15f3-6f.+f,
Blank card 6
-51.+ 15f,-20f.+ 15f3-6f.+f,
I, card
1•.,--31,
Clears
f. -6f,+ 15f. -20f3+ 15f.-6f,
Blank card 2
Blank card 7
1.- 41,
-6f. + 15f3 -2°f. + 15f,
Blank card 2
Blank card 7
+1.
I.-51,
f,-6f.+ 15f3-20f2+ lSI,
Blank card 3
f. card
Clears
-6f.+ 15f2- 2Of,
Clears
Blank card 7
+£,
f. -6f3+ 1512- 2Of,
Blank card 1
f, card
+£,
13- 61.+151,
Blank card 1
f. card
+1,
-61,
Blank card 3
+1,
+1,
1.-61,
Clears
Function
Counter
+f.
f. -5f5+ lOf. -lOf3+ 5f2-f, = IS,
6"+6",=6'.
1S.+6;=IS,
6 23 + 6 33 = 6'.
6'.+ IS. = 6'.
-6f.+ lSI. - 2°f, + 15f3 -6f2+f,
Clears
f,-6f.+ 15f, -20f,+15f3-6f,+f, = 6·,
+f,
*The subscript convention used here differs from that of Yowell in Dr. Blanch's paper.-Ed,
21
22
To explain the foregoing techniques, let us confine
our attention to the sequence of operations occurring
with the fs card. When the fs card is at the lower brushes,
it is added into the 6. s counter and function counter. On
blank card 1 the 6. s counter adds to the 6. 5 counter, the
amount in the 6. s pririts, and the counter is cleared. On
blank card 2 the amount in the 6. 5 adds to the 6. 4 counter
and subtracts from the 6. s counter. A similar process takes
place on blank cards 2 to 5. On blank card 6 the 6. 1 counter
subtracts from the 6. 6 ~ounter. On blank card 7 the function subtracts from the 6. 6 counter. The function counter
is cleared on this card by having the amount stored in the
counter subtract from itself. This method of clearance
eliminates the necessity of having to stop to take a total.
From this point on the pattern described above continues
to produce the successive values of the sixth difference.
I have only indicated the changes which take place in the
various counters. It is clear that if the add impUlse to the
6,1 counter is eliminated we will transmit a zero balance
on blank card 7 and the order of the differences in the
remaining counters will be reduced by one. Similarly, if
the add impulse to the 6. 3 counter were removed, the
process described would generate second differences only.
I f only fourth differences are desired a minor modification
of the above process (using six blank cards) will print all
the columns of differences 'as true figures. In the sixth
difference control panel only the quantities in the 6. 6
counter print so that the differences in the other counters
are carried as complements. If one wishes to list all columns of differences it is necessary to i~troduce additional
counters from which the amounts may be printed as true
figures with appropriate sign indication.
In both cases just described it is possible to difference
ten-digit function values taking account of algebraic signs.
If only fourth differences are desired, the capacity of the
counters may be extended to sixteen digits. It is also possible to take second differences of three functions (using
four blank cards), or third differences of two functions
simultaneously (using five blank cards).
It is clear that the Type 602 Calculating Punch is a
superior machine for the calculation of tabular differences
since it is faster and the results are punched on cards.
Although the results can be punched on cards when differencing on the accounting machine, this procedure involves summary punching. However, when only a printed
record of the differences are wanted for the purpose of
checking tables, the. accounting machine method is desirable. Comparison of the speed of the accounting machine
method with that necessary for differencing oh a National,
Burroughs or Sundstrand machine shows that it is approximately four times as fast except for time consumed
in punching the cards. However, once a card file has
been prepared, further computations can be made with
SCIENTIFIC
COMPUTATION
it, and this usually compensates for the time spent in key
punching. The time necessary for interspersing blank
cards is never appreciable, and this operation can usually
be done at the same time that the differences are being run.
I would like to mention that at the present time we can
compute tables on IBM equipment and type them on the
card-controlled typewriter. From our experience, this typewriter made one error in typing 35,000 ten-digit numbers.
This is work of a high order of accuracy, but for a table
maker there is no compromise with perfection. The resulting manuscript must be subjected to various tests. It would
be highly desirable if there were a means of preparing a
card file from the typed manuscript. This would obviate the
necessity of checking by hand, proofreading or repunching
a new set of cards to be compared with the original files.
111r. Hollander: I would like to suggest something about
the notation on the diagrams. There ought to be more detail
about where the information is coming from. When a
counter is indicated by a column, it requires little more
writing to indicate in that column that the counter is being
impulsed plus or minus. The flow of information through
the machine can then be more easily determined, because
one knows something is happening in that column.
Dr. Herget: I would like to point out that the chain of
X distributors Dr. Abramowitz showed might better be
activated by a punch on one of the function cards, because, if at any time there is a machine failure, the 405
will go on with an out-of-phase cycle of eight. If started
by the card, each cycle would be independent.
lV/r. J-lollander: The blank cards could all carry control
punches.
Mr. Bell: Another advantage of that is that the board
can be used as a general purpose difference computer with
the order of the highest difference determined by the control cards used.
- Dr. Eckert: There is one comment I would like to make
about Dr. Abramowitz' point that we should have automatic means of reading tables so as not to have to key
punch them. If you consider the value of a table and the
amount of work you put into computing, it doesn't seem
excessive to have such a machine. But to key punch the
figures from a paper you publish for posterity is a comparatively cheap method of proofreading. Even though
you check the plates perfectly the printing may have been
bad, so it is very useful to have at least one copy of the
tables looked at by individuals as it comes from the press.
Dr. Abramowitz: The difficulty in the printing doesn't
detract from the value of a machine for reading back and
checking. You want to be able to read back accounting
machine records, too.
Dr. Grosch: I think many people here already know it is
possible to take second differences on a two-brush accounting machine by ~ip-flops, without interleaving blank cards.
The Use of Optimum Interval Mathematical Tables
H.
R.
J.
GROSCH
Watson Scientific Computing Laboratory
ures," but a detailed analysis of the permissible error as a
function of the argument. Thus a typical specification
might be "no error greater than 2.6 X 10-8 for x < 1, nor
greater than 2.6% X 10-8 for % > 1; average error of
random interpolated values to be zero." The average number of values to be taken from the table at each use should
be known, as should the probable number of times the
table will be used (on the work to which its construction
cost will be charged).
In addition, the speed and operating cost of the machines
involved both in using and in constructing the table must
be given, and their storage and sequence capacities. Finally, some sort of estimate may be made of the value of
the table designer's own time; this is the factor which
makes it unwise to spend a month planning how to save a
total of three or four hours of 604 time!
The most satisfactory type of table, given requirements
permitting its use, is the critical table. No interpolation is
required; final answers are obtained by sorting or collating
and gang punching. In printed critical tables, one line is
required for each possible value of the function (more if
the function is not monotonic in the range tabulated).
Values of the argument are constructed so as to correspond
exactly to values of the function midway between those
actually printed. In the optimum interval methods I have
developed, a general form of critical table arises from
setting the degree of polynomial approximation equal to
zero~general in the sense that the maximum tabular error
need no longer be exactly one-half in the last place.
Critical tables, however, are indicated only if the number of values of the function required from each use of
the table is large compared to the size of the table. Thus a
critical table of the sine function from 0° to 80°, with
maximum allowable error 1.5 X 10-'*, will consist of 3283
cards. It obviously would be a good choice if more than a
thousand values were required each time the table was
used, and it obviously would be a poor choice if less than
a hundred values were required.
A table requiring linear interpolation usually is considered next. It is far easier to construct than higher order
tables, and much smaller than a critical table. But machine
THE ART of constructing printed tables of mathematical functions is not by any means static. Indeed, of the
half dozen great table makers, two have flourished in our
time: the late Jean Peters and L. J. Comrie. The requirements of a good printed table are not often explicitly formulated, but most of the Forum members have worked
with "good" and "bad" specimens. Not only must a computer constructing a printed table worry about interpolation methods, tabular intervals, and the detection and
elimination of errors, but he must consider very carefully
the typography, page format, paper quality, and binding.
Discussion of these latter items is not often found in the
literature; aside from notes in MTAC reviews and the
introductory material in the new Chambers tables, l the
only discussion I have referred to recently is in the Napier
Memoria1. 2
Naturally the aim of the table maker is to facilitate the
use of his product. In specifying a figure of merit for usefulness, however, one runs some risk of controversy in
assessing a "design" in terms of speed of use, reduction
of ocular fatigue, and protection against misreading; the
relative weighting of these estimates is even more uncertain. The human element is the vital one in hand computing, and it is not su,rprising that attempts to predict what
computers wi1llike have led to rather varied results!
The situation is far different when we turn to automatic
digital computing equipment on the level, say, of the Type
602. For a given system of input, storage, and output the
variables of typography, format, and binding disappear;
estimates of the usefulness of a particular table can be
made with almost the precision of cost accounting; and
machine characteristics become more important than
human foibles. It is not surprising that machine methods
of calculatimi have led to a very different formulation of
the problem of table design, nor is it surprising that the
formulation can be much more exact than opinions about
printed tables.
A problem involving table lookup should be specified as
completely as possible. As always, one needs to know the
function and range of arguments required. An exact statement of error requirements is needed: not just "five fig-
23
24
SCIENTIFIC
characteristics may be an important factor; if the multiplications involved in the interpolatory process are performed on the 604, for example, quadratic or cubic
formulas take no longer to evaluate than linear ones, and
it is possible to save table bulk subject to limitations imposed by the storage capacity of the 604 and by the
extra cost of computing more complicated tables. On the
602, extra multiplications slow down the operational
speed, and a different balance must be struck. On the
601, extra card passes are required; for quadratic interpolation Herget's card reversal technique eliminates the
extra control panel, but not the second pass.
There is no point in belaboring this subject of economic
decisions too far. I will close this section of my paper by
giving one very detailed example of the process, and then
pass on to more technical matters. Suppose the estimation
procedures to be explained later have been applied to a
certaIn problem, and that the following approximate table
sizes resulted:
18,000 cards
Critical
1,400 cards
Linear
340 cards
Quadratic
110 cards
Cubic
Quartic
60 cards
Quintic
30 cards
(12)
(17)
(23)
(29)
(35)
(42)
The numbers in parentheses indicate the number of digits
punched on each card of the table, including the argument.
The arguments on the detail cards are six-digit numbers.
Further suppose that the equipment available includes
sorter, collator, reproducer; and 602A and 604 punches.
The relative cost of using these machines is taken as
2, 3, 4, 8, and 16 including operator and overhead (these
figur.es will of course differ from installation to installation, and also will be changed for different machine models
and extra accessories). Finally, suppose that an average of
250 values is needed for each use of the table.
In the critical case, we shall do best to use both sorter
and collator. First we sort the detail cards on six columns,
then we collate these on all six columns with the 18,000card table, selecting out unnecessary table cards. The
merged deck, not much under 500 cards, is passed through
the reproducer and gang punched. A final collator run removes the. detail cards and reassembles the big table. The
costs, in arbitrary units, are respectively 0.3, 3.9, 0.4, and
4.0; the total, 8.6. Because of excessive card handling
problems we will adopt 9.0 as our reference figure.
In linear interpolation, collating will still pay in spite of
the drastic reduction in table size. We sort the 250 detail
cards on {our columns, collate on four columns with the
1,400-card table, selecting unnecessary table cards, run the
merged deck of say 400 cards through a gang punching
operation, separate and restore the table as before, and
COMPUTATION
finally make a multiplication of simple A X B + C form
on each- detail card. The costs are 0.2, 0.4, 0.4, 0.4, and 1.4
(Type 602A) or 0.8 (Type 604). Choosing the 604, the
total cost is 2.2 units, a considerable saving over the previous critical table.
In quadratic interpolation, it is no longer economical t.o
select out unwanted table cards, since it is cheaper to pull
apart the table and detail cards on a single pass through
the sorter than to reassemble the table deck by collating,
while the extra cost of running a few unnecessary table
cards through the reproducer is negligible. Hence, the
operations are sorting on four columns, collating without
selection (340
250 cards), gang punching the whole
590-card merged deck, sorting once for separation, and
repeated multiplication of the form (A X B + C) X
A
D on either 602A or 604. The costs are 0.2, 0.2, 0.5,
0.1, and 2.2 (Type 602A) or 0.8 (Type 604). The total is
1.8 if the 604 is used.
In the cubic case the collator is no longer used, and this
will be true for the still smaller tables in higher-order interpolation. The costs are 0.3 for sorting, 0.4 for gang
punching, 0.1 for separating, and 3.0 (Type 602A) or 0.8
(Type 604) for interpolation. The total cost is therefore
1.6 units.
The quartic case costs 1.4 units, using the Type 604,
since the initial sorting cost drops to 0.2 again and the
gang punching drops to 0.3. This total cost of 1.4 arbitrary
units will not be substantially reduced by going to higher
.order interpolation, and in fact the Type 604 has to stop
here, as its capacity for storing the multiple coefficients of
the quintic and higher tables, reading the detail card and
punching the answer is exceeded.
N.o consideration was given to passing the merged deck
through the 604, omitting the gang punching operation. Passing a table card through the 604 costs four
times as much as passing either a table or a detail card
through the repr.oducer; therefore gang punching should
be omitted only if the size of the table is less than onethird the size of the detail deck. This just begins to be the
case for the quartic, and the costs for that case figure
0.2
1.0 + 0.1 = 1.3, undoubtedly the very best that can
be done.
If the number of ~imes the job is to be done is great
enough to warrant constructing the very complicated quartic table, we may claim that this latter case is the most
economical. If figures amortizing the cost of constructing
the various tables are added to the above, a final choice
can be made.
So much for the economics of special tables; now I
want to tell you about the methods we use to design and
construct optimum interval tables of various order;;.
The idea of expanding the interval is not new. 3 ,4 The
exigencies of hand computation prevented adoption in the
+
+
+
FORUM
25
PROCEEDINGS
past, but the advent of punched card equipment in technical computation immediately brought the matter to the
fore. The linear univariate case was discussed in some
detail, 5 but even at this level it is possible to increase the
permitted interval by 40 per cent. G The material for the
general case is new.
Let us consider the Besselian interpolation formula
f=
i£ +
.
n . LYi+1/2
n(n-1) ,,00
2
. L::/l+1/2
+
+ ...
involving odd and mean even differences. The error of
neglecting the third difference is
n (n - 1/2) (n - 1)
6
.
"iii
Ui+1h •
This error is zero at n = 0 and n = 1, and has two extrema of equal size and opposite sign at n = 1/2 ± vT7f2~
The extreme error is y3 iii/216, or about ~i/125. If this
kind of quadratic interpolation is adopted, the rule for
interval would be L1 ii = 125 £, where £ stands for the maximum allowable error due to neglect of third and higher
differences.
If we define the error of approximation of the jth line
of a table as
L
p
f(x) -
Aij Xi,
_t'j <
X
< .t'j+1
-
<,
Li
0
aij n~ ,
0< n
<1,
(2)
i=O
with
n=
.1:"j
.1:" Xj+1 -
Xj
This is ordinarily a polynomial of degree
out £ we write it as
p + 1; taking
(3)
Of course £ may vary from line to line of the table. In the
example of p = 2 Besselian interpolation, Ep(n) was
n/12 - n2/4 + n 3 /6.
Without further ado, I will simply state that the maximum possible interval is obtained when the error polynomial Ep (n) is chosen from
Eo(n)
=
1 - 2n
E1 (n) = 1 - 8n + 8n2
=
1 - 18n
Ea (n) = 1 - 32n
E2(n)
or in general
+ 48n
+ 160n
2
2
-
32n 3
256n3
+ 128n
4
L~i
P+1 (
1)"
1.
(p+i)! • ni
(2i-1)!·(P-i+I)!
(4a)
( -1 )P+1 cos [2 (p + 1) cos-1 yin]
(4b)
2::!H
The identification with Chebyshev polynomials is due to
Dr. C. C. Bramble.
The error of approximation is therefore + £ at n = 0,
passes through p extrema (alternately - £ and + £), and
is ( -1 )P+1 .£ at n = 1. This is the error distribution which
maximizes the interval w = Xj+l - Xj' Note that the tabular
values ~Aij:vJ are no longer' exact except for rounding
error. Instead they are wrong by the maximum allowable
amount + £.
This material is sufficient to handle certain simple cases.
For example, consider the sine function near the origin,
and quadratic interpolation. We can replace the sine by the
first terms of its series, and write,
(x - x 3 /6) - (Ao
+ A1_'l: + A2x2)
£ •
Then using
where the A's are the coefficients of the approximating
polynomial of degree p, we can write this in n-measure as
+ 1)
'i=o
(1)
,
i=O
f ( n)
1 + (p
.t'
=
l1w
(1 - 18n
=
+ 48n
2
-
32n3 ).
and equating like powers of n we find
-1/6n3 w 3 = -32n 3 £
or
the remaining equations can be solved for the A's and the
first line of our table is complete. Next we could write the
Taylor series expansion about -1:"1 and repeat the above to
get the next interval -1~2 - x] and the next set of A's; this
would go along finely until the term +.1,-5/120 began to
bother us.
H one tries this process with a more slowly convergent
series for f (x), even the first line of the table will give
trouble. The square root function is a case in point; I tried
very early in my experiments to construct a table of
~ to be used in going from sines to cosines, but
found thatfC.'\:) = 1 - .1,'/2 - _1:":l/8 - .1,'3/16 - ... converged too slowly for even the first line of a ten-decimal
p = 3 table.
The clue to a more general approach comes from picturing the error polynomial after it is slightly distorted by the
presence of terms of degree greater than p+ 1. The shape
is different; the extrema are no longer exactly -lor + 1 ;
the positions of the extrema have shifted, but not much!
Therefore, if we insist that the error be + £ at n = 0,
(-ly£ at 11. = nk, and (-1)p+1 £ at 1l = 1, where the nk are
the p roots of
dEp(n) = 0 ,
dn
26
SCI E N T IF I C
errors very slightly larger than € may occui' in the neighborhood of the nk; this excess error may be taken care of
by making the interval very slightly smaller than theory
might indicate.
In general, then, after an interval has been adopted, we
compute the function for the p+2 arguments corresponding to 11, = 0, the nk's, and n = 1. This is done with the
usual table-making precaution: two or three extra decimal
places are carried. \Ve then write p+2 simultaneous equations in the unknowns ao, a l , • • • , ap, and €:
ao
a o + n 1 a1 + n 2 G2 + ... +
ao +
J:.G p
a1 + n~G2 + ... + n~ap
1'i 2
+
€
= f(O)
-
€
= f(n 1 )
+
€
= f(n 2)
+ npa1 + 1t~a2 + ... + n~ap + (-l)p € = f(np)
ao +
a1 + a2 + ... + Gp + (-l)P+l f = f(l)
Go
These equations are solved for the a's and for €. If € is less
than but nearly equal to its desired value, the interval has
been chosen correctly.
Due to the nature of the polynomials Ep (n), we are able
to write a general rule for the positions of the extrema:
1lli,
= sin 2 2(;~1) ; k::::; 1,2,3, ... p.
[
, +1/2
+1/2
+1/2
-1/2
]
[
I
+1/2
+3/4
~~/4
o
-1/4
+1
-1/2
+1/4
+5/6
+1/3
-1/3
+1/6
-10/3
+2
+10/3
-2
+8/3
-8/3
-8/3
+8/3
+1/6
-1/3
+1/3
-1/6
+7/8
-7
+1/4
-1/4
+ 4y2
+14
...c..12y2
+4
-4
-8
+ 8y2
0
+1/8
-1/4
+1/4
-J
-1/8
+1/4
- 4y2
+12y2
+3
-10
-
+8
8y2
-1/4
+1/8
A very important feature of the matrix presentation is
that the bottom row of N-t, which is the one used in calculating € from F, has the general form
(5)
Thus there are no extrema for p = 0; for p = 1 there is a
minimum at 11,1 = 1/2; for p = 2 there is a minimum at
n 1 =1/4 and a maximum at n 2 = 3/4; for p = 3 there are
minima at n 1 = 1/2 - Y 1/8 and 113 = 1/2 + Y 1/8, and
a maximum at 112 = 1/2.
It is evident that if we write the above system in matrix
notation
N·A=F,
N is a square matrix of order p+2 whose elements depend
only on the choice of p. The inverse of N may therefore
be computed once and for all, for the various values of p,
and the unknown column matrix A formed from
A = N-1. F .
TABLE
CO M PUT A T I ON
(6)
The bottom element of A is €; from the top down, the
others are ao, av a2 , • • • , ape Having obtained the a's, the
linear transformation to the desired A's is obvious but
laborious, since each line of the table requires a different
transformation. Values of N-1 for p = 0, 1, 2, and 3 are
given in Table 1.
[ +2(P~1)
1
+ p+1
( -l)p
p+1
...
( -1)P+l]
2(p+ 1)
(7)
.
This furnishes the definitive estimation process required
in choosing an interval, and is a powerful tool in polynomial approximation theory.
While these matrix methods serve well in the final
stages of table design and in the actual construction,
rougher procedures are valuable in the preliminaries. If
we consider the term i == P + 1 in (4a), put n = 1, and
equate it to the (p+ l)th term of the Taylor expansion for
f(x+w) we get
w P+1
f
(P+l)
/(p+ I)!
= (-1) PH • 2:W+
1
• €
•
Absorbing the sign' into €, and remembering the relation
between differences and derivatives, we can write
FORUM
27
PROCEEDINGS
The more usual cases are as follows:
p=o
P= 1
p=2
p=3
L,1 =
L,ii =
2£
16 £
L,iil = 192 £
L,iv = 3072 £
Suppose we wish to make a linear table of In x, 1.0 <
x < 5.0, with an overall accuracy of 1.8 X 10-6 throughout, and no requirement on the mean error. First we have
to consider the form of the table. I have found from experience that in this sort of problem there is a real gain in
making the Ao's seven-decimal numbers, but not more.
The A/s will also be seven decimals. The rounding error
of the worst Ao will not exceed 5 X 10-8 , the rounding
error of A1x (assuming % to be exact) will not exceed
5x X 10-8 • The final rounding of an interpolated answer f
may introduce errors as large as 5 X 10-7 if the answer is
given to six decimals only. Therefore £, the maximum allowable error due to the degree of approximation, must be
£ = (125 - 5x) X 10-8
•
From (8) we have
00
2
fii = (20 - 0.8x) X 10-6
;
if it is desirable to squeeze the final table down to the very
smallest size, this equation may be used to calculate the
intervals. Usually, however, we would take the worst case
(x = 5.0) and simply use
2
fii = 16 X 10-6
00 2
= 0.004(fii)-1/2
00
or
as the interval rule. For the natural logarithm, fi = -x ~
and hence the very simple result, giving £ a negative sign, is
w
= 0.004%
.
As a test we take the last interval, with argument 4.980
and interval 0.020. The column matrix F is
0.80 X 10-6 , so the above value of £ is exactly as it should
be. The actual tabular values turn out to be
Ao = +0.6074339
Al
and application of (7) gives £ = -1.005 X 10-6 • At x = 5.0
the maximum rounding error from all sources will be
+0.2004010 ,
and the errors at n = 0, n 1 = 1/2, and 1Z, = 1 are respectively -0.99,
1.02, and -0.99 in units of the sixth
decimal. These errors are small because there happens to
be no error in Al (at least to eight decimals), the rounding
error of Ao was only 1.5 X 10-8 , and the final results were
not rounded at all. This favorable combination of circumstances will not hold throughout the table.
As for the size of the whole table, detailed examination
using values of 00 from 0.004 to 0.019 will show that 430
cards are required. For rough estimates, I use
+
fd~/w = 250 f d:/x
J1.O
J1.0
= 250 In 5.0 = 403
This, of course, assumes a smooth change of interval, and
is usually ten or fifteen per cent low.
It is not practical to describe here the extension of these
methods to other problems. I have worked out approaches
to such variations as cases with error terms of order p+2
instead of p+ 1, cases where the mean error must be zero
and p is odd (for p even the above methods suffice), and
even the bivariate problem. In general, it is fair to say that
optimum interval methods in the latter case are warranted
only if the table is to be referred to many hundreds of
times, for the work of design and construction is enormous
even for the linear ca~e.
REF'ERENCES
1. 1. J. COMRIE, SLr-Figure Mathematical Tables (Chambers,
194~),
2.
Vol. II.
J. R. MILNE, "The Arrangement of Mathematical Tables,"
j\fapier Tercentenary M etl10rial VOht1ne
3.
4.
1.60542989 ]
1.60743591
[ 1.60943791
=
5.
6.
(Longmans Green,
1915), pp. 293-316.
]. E. A. STEGALL, "On a Possible Economy of Entries in a Table
of Logarithms and Other Functions," lac. cit., pp. 319-28.
P. HERGET, "A Table of Sines and Cosines to Eight Decimal
Places," Astron. I., 42 (1932), pp. 123-25.
P. H~:RGET and G. M. CLEMENCE, "Optimum-Interval Punched
Card Tables," MT A C, I (1944), pp. 173-76.
J. c. P. MU,LER, "[Note on] Optimum-Interval Punched Card
Tables [of Herget and Clemence]," MTAC, I (1944), p. 334.
DISCUSSION
[Discussion of this paper was omitted because of time limitations.]
Punched Card Techniques for the Solution of Simultaneous
Equations and Other Matrix Operations
WILLIAM
D.
BELL
Telecomputing Corporation
ELECTRICAL punched card equipment has been used
for matrix calculations of various sorts for some time.
There have been wide discrepancies in the operational
methods, efficiency and general utility of the procedures
being used. From the widespread interest in this subject,
it is evident that there is a genuine need for a good basic
approach to the problem, particularly in terms of actual
machine operations.
The method explained below has been successfully used
for a wide variety of problems over an extended period
of time. Among the advantages of the system are the following:
Given a set of simultaneous linear equations
3X
+ 12Y -
+
7Z
8Z
=
=
6
2X -
8Y
10
6X -
2Y -- 3Z = -7
write the matrix
3
12
-7
6
2
-8
-2
8
-3*
10
6
-7 .
The first reduced matrix, by "starring" on the A33 term,
becomes
1. Both card handling and machine operations have
been reduced to a minimum. This results in a definite
time advantage, as well as simplicity from the operator's standpoint.
3 _ (-7)(~
(-3)
-12 _ (-7)(-2)
(-3)
(-7)(-3)
-7---(-3)
2-~-~
_ 8- (8)(~
(-3)
8- (8)(-3)
(-3)
(-3)
2. Among the basic matrix operations to which the
same procedure is directly applicable, are the following:
a. Solution of systems of simultaneous equations.
b. Computation of determinants.
c. Calculation of inverse matrices.
(-7)(-7)
6-<=3)
10 _ (8)(-7)
(-3)
,
which reduces to
-11
18
6
50/3
-40/3
2
0
0
-3
67/3
-26/3
7 .
Repetitions of this pattern, working each time with the
transformed matrix from the previous operation and
starring on one of the remaining main diagonal terms, will
produce a diagonal matrix from which the unknowns or
the determinant can be computed directly. The basic operation is thus of the form
3. The same method, procedure, and control panels can
be used for either real or complex numbers.1
The equipment described is that made by the International Business :Machines Corporation.
A .. _
'tJ
Mathematical Basis of M etlzod
A'im A'III} -
Arnrn
-
'A ..
'tJ
where' Aij is the A'ij term in the transformed matrix and
Ainrn is the starring term in a matrix of order n.
There is a systematic method of operating on a matrix
in such a way as to make all of the elements of a given
column, except one, equal to zero. 'I'he method consists of
subtracting some multiple of a given row from each of the
other rows in such a way as to make one column zero in
'all except one element. The process has been well described
in the literature 2 , 3 and is often called the "starring process." A simple numerical example follows.
IBA1 M etlzod
The major problem in matrix work is that a given card,
representing a certain element of the matrix, must be
handled in a different manner at different steps in the
28
FORUM
29
PROCEEDINGS
procedure. As an example, what is an answer from a previous transformation might become a multiplier in the
next reduction. An obvious solution is to use a master deck
which knows the factors to be used, the operations to be
performed, and the identity of the resulting term. The
actual numerical values for a particular problem are then
transferred to the master cards by a relating and gang
punching process. A method of this type has been published. 4 The disadvantages of such a process are the necessity of constructing the master files, and the number of
unnecessary operations involved.
The key to the procedure explained here lies in a very
simple way of causing a machine to differentiate automatically between two identical cards. The usual method
of identifying a type of IBM card is with an X punch in
all cards of that type. This is a high punch placed in some
column of the card and falling at the upper edge. It can
be used for controlling the various functions of the different machines through which the card passes. If a card
goes through a machine face up rather than face down,
all punching in the card is transferred to a mirror image
position with right and left transposed. It becomes, to all
intents and purposes, a different kind of card. The fact
that this requires the machine to read each amount from
different columns and from right to left instead of left to
right imposes no restrictions. The card form used is shown
in Figure 1. The card is in its normal reading position;
when the card is in its reversed or mirror image position,
the machines sense X58.
I
I
GANGED
~
OOOOOOO~OOOOOOOOOOOOOOO
123456789WngU~mffinm~w~nn
The control panels required are as follows:
1. A gang punching control panel which transfers the
factors of the starring row into all related cards. The
starring row becomes gang punch master cards
which are used in the "turned" position. This function is performed automatically by a Reproducing
Punch, Type 513 or Type 519, at a constant speed
of one hundred cards per minute.
2. A Type 602 Calculating Punch control panel which
performs the operation
'Ai' - A .. _ Aim Am;
J IJ
Amm
The factors Aim and Amm are in "turned" cards
which are recognized by the machine by a common
X58. The ratio of these terms is computed by a
division operation and the result stored within
counters in the machine for use on all terms in a
particular row of the matrix. This eliminates any
redundancy in the division operation, which is the
slowest operation performed, requiring approximately four seconds per card.
On an X23 card (normal position) the machine
multiplies the ratio by Amj and subtracts it from Aij
to form the new term' A ij , which is punched into
each detail card. The order number of the new
matrix term is simultaneously calculated and re-
REPRO. FROM PREVIOUS ORDER
513
AMOUNT
REVERSED
Control Panels
""~
~~
v:>
• ""v:>
"" ...... "" ""
...... CD ......
......
z
......
~
CD
!~~ ~~ ~~~ ~~
0000000000
0
ct:
::It
AMOUNT
-'
:; :!:
602
"" ""
......
~§
AMOUNT AFTER
REDUCTION
c:>
""""
...... ::>
~~
o0 00 000 000 o0 00 00001000 0 0 0 000 000 010000 0 0 00 00000000
4 25 26 27 2829 30 31 32 33 3435 3637 383940 414243 4445 4647 ~8 49 50 51~2 5354555657 585960 61 62636416566.6768 69 70 7172 737475 76 77 78 79 80
11111111111111111111111 1111111111 11 11 1 1 1 1 1 1 1 1 1 1 11,1:,1,,11 111 1 1 1 1: 1 1 1 1 1 1 11 1 1 1 1 1 1 1 1
I
I
22222222222222222222222 2222222222 22 2 2 222 222 2 2 22 2 2 2 2:2 2 2 2 2 2 222 2 2 2 212 2 2 2 2 2 2 2 22222222
I
I
33333333333333333333333 3333333333 3 3 3 3 333 333 33 3 3 3333:333 333 333 3 3 3 3:3 3 3 3 3 3 3 3 33333333
I
44444444444444444444444 4444444444 44 44 444 444 44 44 4 4 4 4:4 4 4 4 4 4 444 4 4 4 4:4 4 4 4 4 4 44 44444444
I
I
55555555555555555555555 5555555555 5 5 55 555 555 55 5 5 5555:555 5 5 5 555 5 5 5 515 5 5 5 55 55 55555555
I
I
66666666666666666666666 6666666666 66 66 666 666 66 66 6666:666 666 666 66 6 6:6 6 6 6 6 6 66 66666666
I
I
77777777777777777777777 7777777777 7 7 7 7 777 7 7 7 77 7 7 777 717 7 7 777 777 7 7 7 7:7 7 7 7 7 7 7 7 77777777
I
1
I
I
I
8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8888888888 88 88 888 888 88 88 8 8 8 8 8 8 8 8 8 8 888 8888188 8 8 8 8 88 88888888
9999999999999999999 9 9 9 9 9 9 9 9 9 9 999 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9999:9999 9 9 9 9 9 9 9 9 9:9999999 999999 9 9 9
1 2 I ~ ~ 5 6 7 8 9 10 n 12131415 16 n 18 1920212223242526272829303132333435363138 3940 41424344454647484950 51152535455565158 59 6~ 61626364165666768697011 7z 1314757677 78 79 80
FIGURE
1
-
30
SCIENTIFIC
corded. The speed here is about twenty-five cards
per minute. All of these mathematical manipulations
are completely automatic, with the 602 sensing the
type of card, determining what is to be read, and
performing the proper operations.
3. A reproducing control panel which at a hundred
cards per minute will make a new deck with the
newly computed transformed elements returned to .
the original locations in the new cards in readiness
for the next starring operation.
Procedure
The detailed steps followed by the operators are exceedingly simple. A complete operating procedure is given
below. This covers the calculation of the diagonal matrix.
This method, as explained above, will give the solution of
simultaneous equations, inverse matrices and determinants.
1. Sort the cards to column 44-45.
2. When column 44-45 equals 59-60, reverse the selected cards, place in front and sort:
46-47 minor
41-43 major
3. Gang punch using control panel (1).
4. Select the X58 (turned) cards and hold aside.
5. Sort remaining cards to 46-47.
6. When 46-47 equals 59-60, reverse the selected cards,
place in front and sort:
44-45 minor
41-43 maj.or
7. Calculate 'A ij , the elements of the transformed
matrix, using board (2).
8. Select X58 cards and hold aside.
9. Reproduce the balance of the cards on control panel
(3).
This completes one starring operation. Using the new
cards from (9), begin again at step (1) and repeat.
Usually it will be adequate to select the "starring" terms
by· going up the main diagonal.
RE;F-E;RENCES
1. W. D. BELL, "A Simplified Punched Card Approach to the
Solution of the Flutter Determinant," f. Aeron. Sci., 15 (1948),
pp~ 121-22.
2. R. A. FRAZER, W. J. DUNCAN, and A. R. COLLAR, Elementary
Matrkes (Macmillan, 1947), chap. IV.
3. H. MARGENAU and G. M. MURPHY, The Mathematics of Physics
and Chemistry (Van Nostrand, 1943), pp. 482-86.
4. P. D. JENNINGS and G. E. QUINAN, "The Use of Business Machines in Determining the Distribution of Load and Reactive
Components in Power Line Networks," AlEE Preprint 46-195
(1946).
,""
COMPUTATION
DISCUSSION
Mr. HaY11'tan: When you use the cards face up, don't
you have trouble with curvature?
Mr. Bell: No trouble at all.
Dr. Caldwell: What kind of climate do you have?
,lllr. Bell: Los Angeles.
Dr. Herget: We could do it at Washington.
Dr. Grosch: We did it in New York.
Mr. Bell: This approach of reversing the cards will
work in many problems. When you reverse the cards in
the matrix work, the identification field is reversed also,
so it is necessary to have another field where the identification is punched as a mirror image. Then when you turn
them over and sort you do not have to invert the sorting
procedure or change the sorter brush setting.
MI'. Ferber: Are your elements real in this case?
Mr. Bell: Complex.
.Alr. Ferber: In finding characteristic roots, you make a
guess from previous experience?
Mr. B ell~' Right.
Mrs. Rhodes: You have been very lucky with fiftieth
order linear equations. What did you do about loss of
accuracy, go on faith?
Mr. Bell: Yes, I went on faith.
Mrs. Rhodes: Did you position for division? If so, how?
Mr. Bell: You are talking about a very real problem.
Occasionally we get a matrix that doesn't want t.o come
out. I have yet to find a way by which you can easily and
quickly predict that you are about to get into that situation. And we have been working with relatively small
problems.
In working large problems the gang punch master cards
can be deleted at each stage. Then a straightforward back
solution is used after the first unknown is obtained. This
does not work well if you have just a few equations. It is
more work to do the special back solution than to carry
the extra terms along as you go.
Dr. Grosch: There are still only a few installations
which have had experience with problems of very large
order. We did a problem with forty-one unknowns, carrying the full eight figures permitted by the 601. Our
procedure for size was simply that we never starred a
number which was less than one-tenth as large as the
largest element in its row. We looked at a printed record
produced by the 405 while applying the usual check
sums. In mass production one could suppress listing for
all but the starred row to save paper.
Following that rule, we had to rearrange only once in
forty reductions. We lost only two significant figures. How
many figures did you lose in your fiftieth order problem?
Mr. Bell : We had four or five decimal accuracy.
Dr. Caldwell: That is about our experience in large
cases.
FORUM
31
PROCEEDINGS
Mr. Ferber: If you want to solve a small set of equations with seven or eight unknowns, the machines aren't
adapted to working with such a small number of cards
and one comes to the point where it is quicker to do them
some other way. Then again, since the work goes up very
rapidly with increasing order, you soon come up against a
blank wall in that the work is prohibitive. How do you
handle the problem?
Mr. Bell: We do enough matrix work to keep the basic
602 panels wired up ; so we need not allow time for that.
Our operators are familiar with the job, and they handle
it efficiently even with only seven or eight unknowns.
Everything is kept ready for them. To give you some idea
of time, it requires between one and two hours to evaluate
an eighth order complex determinant. I f there are enough
problems, you can split them up so that some are mUltiplying while others are being sorted and gang punched, with
no idle machines.
We have tried to handle special cases by special procedures, but have lost on it. I feel it is better not to use
special procedures. Instead we stick to methods that all of
our people understand; then we are more certain of coming out with a good answer.
Mr; Harman: I was wondering if you could get around
the division step in your formula by the expansion of a
second-order minor, which involves the difference between
two multiplications.
Mr. Bell: We have tried doing that, and here is the.problem we got into: the numbers change tremendously in
size, and it is necessary to stop for inspection. That takes
longer than the regular method.
Dr. Tukey: The thing is to get away from division by
small numbers.
Mr. Bell: Yes. 1\1:ost of our work is brought in to us.
Sometimes we know what is behind the problems, and
sometimes we do not. Most of the work is in engineering
fields where extreme accuracy is not required. But we do
have to calculate with expanded accuracy occasionally.
Dr. Tttkey,' As I understand the situation, the pessimists
thought-and I was one of them-that you lost, roughly
speaking, a constant number of significant figures for each
new equation. Now it has been proved! that it doesn't go
like that; the loss goes like the logarithm of the order.
If you are thinking of big equations that is a tremendous
improvement. I think that size ought not to be taken as
grounds for pessimism.
REFERENCE
1. J. VON NEUMANN and H. H. GoLDSTINE, "Numerical Inverting
of Matrices of High Order," Bull. Amer. Math. Soc., 53 (1947),
pp. 1021-99.
Two Numerical Methods of Integration Using
Predetermined Factors
LELAND
W.
SPRINKLE
Bethesda, Maryland
V AR IOU S systems have been proposed for numerically
integrating expressions for which no formal method has
been found. Although many of these systems are highly
accurate, the adaptability of some of them to mechanized
methods of computing has not been fully explored.
The central-difference formula subsequently referred to
offers an unusually accurate method of determining the
definite integral but at the same time presents practical
difficulties in the calculating and handling of differences
through the seventh order. The following discussion, with
corresponding examples, will show two ways of overcoming the practical difficulties involved.
A simple, accurate, flexible method of determining the
value of the definite integral may be described as follows:
1. Determine the ordinates between and including the
specified limits of integration, maintaining a constant
interval.
2. Multiply the first nine ordinates respectively by the
following nine factors:
Example:
x
.40.0.
.425
.450.
.475
.50.0 .
.525
.550
.575
.60.0.
.625
.650.
.675
.700
.725
.750
.775
.800
.825
.850.
.875
.900
Find the value of
Ordinate
.0.0.163840.0.0.0.0.
.0.0.250.450.850.2
.0.0.3736694531
.0.05455760.125
.0.0.781250.0000
.01099297205
.01522435234
.0.2078140.531
.02799360.0.0.0
.... .03725290.298
.0.490.2227891
~
~
.06384492921
V,
"- .0.8235430.0.00.
.1052848896
.1334838867
.1679236701
.2097152000
.2601223326
.3205770883
.3926959038
.47829690.0.0
""'
~
..,
~
~
"-
S~'
X'dX
Factor
.33380.470.13
1.3289260.91
.687220.8443
1.261894566
.8333333333
1.071438767
.9794458223
1.00.4407242
.9995286325
""'
to:>
""'
..,
~
~
"-
.9995286325
1;00.4407242
.9794458223
1.071438767
.8333333333
1.261894566
.687220.8443
1.328926091
,33380.470.13
Pro.dltct
.0.00.54690.56
.0.0.332830.67
.0.0.25679344
.00.68845941
.0065104167
.0117782964
.0.149114283
.0.208729940
.0.27980.40.47
.0372529030.
.0490222789
.0638449292
.0823154809
.10.57489056
.1307402352
.1799199300
.1747626667
.3282469580
.2203072573
.5218638324
.1596577538
2.149064412 (Step 4)
0.3338047013
1.328926091
0.6872208443
1.261894566
0.8333333333
1.071438767
0.9794458223
1.004407242
0.9995286325
Step 5. The sum in Step 4 mUltiplied by the constant
interval (.025) of the example equals .05372661030 which
is the value of the integral
1:
X'dX, with an error of
only .00000012905.
In the event that four extra ordinates can be accurately
determined at each end of the original group of ordinates,
the following even more accurate method is applicable:
3. Multiply the last nine ordinates respectively by the
same factors in reverse order.
4. Add all the eighteen products formed in the above
steps to the unmultiplied remaining ordinates.
5. Multiply this sum by the constant interval; the resulting product is the definite integral between the originally specified limits of integration.
1. Determine the ordinates between and including the
specified limits of integration, maintaining a constant interval; then determine four more ordinates at each end of
the group, still maintaining the constant interval.
2. Multiply the first nine ordinates of the entire group
respectively by the following nine factors:
32
FORUM
33
PROCEEDINGS
Step 5. The result of Step 4 multiplied by the constant
interval (.025) of the example equals .05372648125 and
0.0004713679453
- 0.004407242063
0.02055417763
-0.07143876764
0.5
1.071438767
0.9794458223
1.604407242
0.9995286325
is the correct value of the original integral,
3. Multiply the last nine ordinates of the entire group
respectively by the same factors in reverse order.
4. Combine all eighteen products with the unmultiplied
remaining ordinates.
5. Multiply the aboye sum by the constant interval ;,the
resulting product is the definite integral between the originally specified limits .of integration.
Example: Find the value of
.300
.325
1.0002187000000 ~
~
1
.350
.375
. 400
.425
.450
.475
.500
.525
.550
.575
.600
.625
.650
.675
.700
.725
.750
.775
.800
.825
.850
.875
.900
r5~
.950
.975
1.000
::;~
~
V)
'-
.0003829865538
.0006433929688
'.001042842865
.001638400000
.002504508502
.003736694531
.005455760125
.007812500000
.01099297205
.01522435234
.02078140531
.02799360000
.03725290298
.04902227891
.06384492921
.08235430000
.1052848896
.1334838867
.1679236701
.2097152000
.2601223326
.3205770883
.3926959038
.4782969000
1.5794181954
.6983372961
.8375915935
1.000000000
X' dX
P,actor
Ordinate
X
S:'
!
2.149059250 (Step 4)
X'dX,
to the number of places shown.
In either of the foregoing numerical integration systems
the number of ordinates chosen may be either even .or odd.
Also, advantage is gained in an extended calculation since
the majority of the ordinates are simply added without
change. Further, there is no limit to the number of ordinates that can be used although only nine at each end of
the group are multiplied by factQrs.
In view of these advantages, the establishment of two
permanent factor files of only eighteen cards each is all
that is necessary to make both systems readily adaptable
to punched card equipment.
The two preceding systems of integration are based on
a central difference formula carried through the fi fth
differences by Scarborough. 1 • 2 The author is indebted to
Dr. J. W. Wrinch of the David W. Taylor Model Basin
for deriving an additional term containing the seventh
differences. 3
Product
.0000001031
.0004713679453
- .004407243063 -.0000016879
.0000132244
.02055417763
-.000-0744994
.07143876764
~ . 0008192000
.5
~
.0026834275
'" 1.071438767
~
'. 0036598898
.9794458223
. 0054798050
1.004407242
.0078088174
.9995286325
.0109929721
.0152243523
.0207814053
.0279936000
.0372529030
.0490222789
.0638449292
.0823543000'
.1052848896
.1334838867
.1679236701
.2096163471
.Y995286325
.2612687547
1.004407242
.3139878899
.9794458223
'%""'
.4207496150
1.071438767
~
2391484500
.5
'"
~
-.0413929218
.07143876764
'- .0143537488
.02055417763
- .004407242063 -.0036914689
.0004713679
.0004713679453
S:
REFERENCES
1.
J. B. SCARBOROUGH, Numerical Mathematical Analysis (Johns
Hopkins Press, 1930), pp. 128-41.
2. E. T. WHITTAKeR and G. ROBINSON, The Calculus of Observalions (1st ed.; Blackie, 1924), pp. 146-47.
3. L. M. MILNE-THOMSON, The Calculus of Finite Differences
(Macmillan, 1933), pp. 184-87.
DISCUSSION
Dr. Blanch: These factors are related to those 111 our
volume of Lagrangian coefficients. 1
Dr. Grosch: It is possible to generate integrati.on coefficients of this sort based on other than polynomial approximation. Unfortunately in most cases the weighting factors
become functions of the limits .of integration, but if the
latter are fixed for a large group of problems one may find
it economical to derive factors based on an approximating
function that suits the physical situation .
REFERENCE
1. MATHEMATICAL TABLES PROJECT, Tables of Lagrangian Interpotation Coefficients (Columbia University Press, 1944), pp.
390-92.
Integration of Second Order Linear Differential Equations
on the Type 602 Calculating Punch
N.
ARNE
LINDBERGER
Royal Institute of Technology, Stockholm
T R I SPA PER is concerned with the numerical
solution by means of punched cards of the equation
In order to eliminate the fourth derivatives, a modified
function y can be defined
= g(%)S + p(%) with initial value conditions. In the
dd2~
second section, the special case of p = 0 is treated. By
y=S-l2"Sii
h2
%~
(4)
Insertion of (2) and (3) gives the value of y at the
p.oint %n+1
approximation to a difference equation, a step-by-step
procedure is derived which has been set up for automatic
integration on the Type 602. The machine operation is
described in the third section. In the fourth section, the
method is extended t.o the general second-order linear
differential equation which can be reduced to the case
above, where P =F O. This can be set up on the Type 602
with a slight change of the control panel described in the
third section.
and by reversing the sign .of h
Appro%i11wtion by Finite Difference$
(6)
Consider the second-order linear differential equation
The central difference of y around the point
(1)
~2Yl~ =
where 9 is a given function over the integrati.on interval.
Suppose that in the neighborhood of an arbitrary point
Xn of the integration interval, a solution S and its second
derivative can be developed into a Taylor series. Let the
constant step in % be h. The value of the function S at the
point X n+1 = Xn + h, is then
}'1I+1 -
2Yn
+ Yn-1
Xn
is now
.
(7)
•••
(8)
Insertion of (5) and (6) yields
~.,
O~Yn
= h')~ Sii
n -
6
h svi+
240
n
Terms of the sixth and higher orders are going to be
neglected and thus form the truncation error. The approximation gives
(9)
Insertion of (1) in (9) and (4) yields
~2Yn =
gn h 2 Sn
(10)
and
(11).
S
Here S,! = (dd )
x
Xn
and so forth.
After rewriting equati.on (11)
··
S~~1
2
where
3
.. 2T
h Siv
h sv
= S"~ + hSJl+
n + TI n + ...
+ P-n) Yn
(12)
gn h2 / 12
1 - gnh2/12
(13)
Sn = (1
The same development of the second derivative gives
(3)
p,n=
34
FORUM
35
PROCEEDINGS
Now the approximation of the difference equation
(10) will be
(14)
with
(15)
Yn = 12p.n
The derivation given here* is essentially the method
used in a paper by Feinstein and Schwarzschild. 1 The
main machine used in their work was a special multiplying
punch.
In order to apply equation (15) for computation on the
Type 602, it will be written in a different manner. With
the advancing difference notation
~'Yn
= Yn+l
(16)
- Yn ,
equation (14) is
~'Yn = .6.}'n-1
+ YnYn
.
The initial values are YI and D.Yo' In (17), putting n
~YI
=
.6.Yo
+ YIYt
.
(17)
= 1,
(18)
From (17) it is evident that every difference .6.Yn is computed from the previous one by adding YnYn. This procedure will give
~Yn
=
n
2:
~'Yo +
YiJ'i •
(19)
i=1
After having used (19), the next value of }' is computed
from (17)
(20)
The Jvl achine Set-up
In order to describe how the equations (12), (19) ,and
(20) are set up for automatic integration on the Type 602,
a flow chart will be used. In this, the first eight columns
represent the eight counters MC, MP, RHC, LHC, Summary Counters 13, 14, 15 and 16. The chart is divided into
ten rows, indicating machine cycles, which are shown in
the ninth column.
Suppose that the integration has proceeded up to the nth
integration cycle (not to be confused with machine cycles).
This means that n -1 cards have been run through the machine, each one of them carrying the functions Y and p..
Every card carries the integration one step forward and
is punched with the computed values of y and S.
The machine has computed Yn during the (n - 1) th integration cycle and keeps the absolute value of this number
stored in the MP counter during the reading of the nth
card. In order to perform its function correctly, MP must
contain a true number. The sign of Yn has to be stored elsewhere, for which purpose SC 13 is reserved. I f the con*From lectures by Dr. L. H. Thomas.
tents of this counter are zero or a true number, it is interpreted as a positive sign in Yn' If the contents are a complement number, it means that Yn is less than zero.
Moreover, the machine stores the progressive differences of Y in SC 14 and 16 which are coupled together.
In the beginning of the nth integration cycle, the stored
difference is
fI, 1
.6.YIH = .6.Yo
+
2:
YSi
(21)
i=1
as indicated in the flow chart (Figure 1, page 36).
The essence of the procedure is the following. During
the nth card operation, the machine is going to compute
the product YII)'n and add it to .6.Yn H thereby yielding .6.YM
according to equation (19); .6.Yn is added to Yn, giving
l'
• Yn is multiplied with 1 + rlin which gives Sn by (20)
.)'1£+1,
and (12). Finally Sn and Y,HI are punched out on the nth
card.
The details of the operation appear in the flow chart.
During the card reading cycle, the precomputed numbers
"'In and 1 -I- p'n are fed from the card into MC and SC 15.
The numbers in the MC and MP counters are then multiplied. The rounded product comes out in LHC. This
product is transferred to SC 14 and 16 on the first program cycle. As it appears positive in LHC, the sign has to
be taken care of by adding the product intO' SC 14 and 16
if it is positive, and subtracting it, if negative. This is
accomplished in the following manner: an eventual sign
in "'In appears as an X punch in the card, wired to pick up
an entry control selector. A negative sign in Yn appears as
previously mentioned in SC 13 which therefore is balancetested during the first program cycle. If Yn is negative, the
test impUlse will get through the balance control. From
there it is wired to pick up a selector. By coupling the
above-mentioned selectors in series, a plus or a minus
shot will be available from them, due to the sign of the
product. The impulses are then fed into the appropriate
control hubs of SC 14 and 16 which will take care of entering the product into the counters with the correct sign.
The details of the wiring are shown in the control panel
diagram.
The quantity 1 + P.n is transferred into lVIC on the second program cycle and multiplied with J'n during the third
program cycle. The product Sn will appear in LHC where
it is stored until the transfer-to-storage cycle. The sign is
taken care of by balance-testing SC 13, as the sign of Sn
is always the same as that of Yn. That sign will now be
stored in LHC; SC 13 can be used for other purposes
and is reset on the fourth program cycle.
On the fifth program cycle, MC reads out Yn to SC 13
and 15, coupled together. lVIC is reset on the same cycle.
As MC always stores true numbers, the sign is taken care
of by balance-testing I~HC, picking up a selector and feed-
36
MC
tn
SCIENTIFIC
MP
i
RHC
Contains sign
of Yn
IYnl
{ Yl}
sc 14
sc 13
LHC
sc 15
sc 16
~Yo
1 + fn
+ t)'iYi
Multiply
RO to SC 14 &- 16,
RC SignCtrl. by Bal.
Test of SC 13 and
Entry Ctrl.Sel. 2
Balance Test
Yn
Yn
RO to MC
RO to sc 13
& 15, RC
RC
(1 + j1h) Yn
= Sn
Balance Test
3rd Program
Cycle
RC
RC
4th Program
Cycle
Yn
Yn
5th Program
Cycle
RO to sc 13
RO to MP
Balance Test
Iyn+d
R~o
Storage,
1st Program
Cycle
2nd Program
Cycle
Balance Test
Yn + b,.Yn
RO Sn
Read card
[A Yo}
1 + j-lh
I
SEQUENCE
1"'1
ot1n
RC
I
+ r'tiYi
1=1
[A Yo]
I
RC
~Y
COMPUTATION
RO Yn + 1 to
Storage
Yn
FIGURE
ing an appropriate impulse to the plus or minus hubs of
SC 13 and 15. As distinguished from MC and MP, the
summary c.ounters and LHC are used to store negative
numbers as complements.
On the-sixth program cycle, 6Yn is read out from SC 14
and 16 and accumulated in SC 13 and 15. This means that
the operation (20) is carried out.
During the seventh program cycle, Y1£+1 is read into MP.
For reasons previously mentioned, it has to be subtracted
into MP if negative or added if positive. This is again
taken care of by balance-testing SC 13.
The program cycles are now finished and the machine
transfers to storage. Sn is read out and punched from
LHC,3'n+1 from SC 13 and 15. By using the balance punch
feature, negative signs will appear as X punches.
LHC and SC 15 are reset, as distinguished from SC 13
which has to store the sign of Yn+l for the next card. If the
6th Program
Cycle
RO to MP
7th Program
Cycle
to
RO Yn + 1
Storage RC
Transfer to
Storage Cycle
1
Dn+1
RO to sc 15
+ b,.Yn
+ j!h+1
Read Next
Card
1
sign is negative, there will be a nine in the left-most position of SC 13 that will be detected during the balance tests.
The allowed number .of digits in the function Y must be
limited to nine as there are ten positions in SC 13 and 15
together and the left position has to be left free.
The initial values Yl and 63'0 are fed into MP and SC
14 and 16 from the first card. In the flow chart they appear within braces. As all counters are reset before the
starting of a new integration, the zeros in SC 13 will indicate the positive sign .of the initial value Y1The control panel has been wired as shown in the wiring diagram (Figure 2) and the procedure has been found
to work satisfactorily_ The coefficients of integration,
yand 1 + p, are conveniently computed with punched card
machines_ In the mentioned set-up p was calculated in .one
run of the 602 ; y and 1
p came out together in one run
on the 601.
+
....-_ _..,....c<=::;;....._ _
Sig~fJn
2/ 3
__
4
o
CALCULA TlNG PUNCH - TYPE 602 - CONTROL PANEL
5
6
7
8
9 10 )J J2 J3 J4 J5
5-CONTROL--10-READING--15
0
0
0
0
0
J6
..
J7
J8
J9
20 2J 22
20-rSELEC
23 24
TORS
25
i
26
27
28
29
30 3J
11
32 33
CARD
34 35 36
READING
37
38
39 40 4J 42 43
15-------
n'
35
oC6
0
o
oTo
6
7
8
ONO
0
ONO
0
10
0
9
VJ
'-.J
o
~
o
FIGURE 2
0
0
0
0
0
0
0
000
0
0
000
~
0
0
0
0
0
~
38
SCIENTIFIC
The General Case
Any linear second-order differential equation can be
reduced to the form 2
After insertion of z, the equations (23), (24) and (20)
will be
(26)
Sn = (1 + /J-n) Zn
(22)
where the first derivative is lacking. This equation, inserted
in (9) and (4), will give recurrence formulas corresponding to (12), (19) and (20) :
L.(
COMPUTATION
~zn
=
~zo +
L.(
n
"Ii Zi
+ Pi)
(27)
i=l
(28)
Here
(23)
(29)
(24)
Equations (26) and (28) are basically the same as (12)
and (20) of Section 2. Equation (27) is taken care of in
the same manner as (24), by feeding the additional term
P n from the nth card into the counters that accumulate the
difference in z. As P can be precomputed in one run on
the 602, the extra run of equation (23) is saved.
Thus it has been shown that the general second-order
linear differential equation can be automatically integrated
on the 602. Reduction to the form (22) has not been
treated here. It will probably imply the same amount of
work as the rest of the procedure.
fI,
~Yn == ~yo +
"tiYi
+ q'i)
i=l
(20)
The addition of qn = h 2 Pn (1 + /J-n) in each step is what
makes (24) differ from (19). It can be done in the setup
of Section 3 by feeding qn from the nth card into SC 14
and 16 which accumulate the difference ~Yn-l'
As for (23), there is no counter storage left for the
additional term qn/12. This can be overcome, however, by
letting the machine compute (1 + /J-n) Yn in the integration
run and then add the above-mentioned term in an extra
run. The calculation of S can be split in this manner because it is not part of the progressive computation of Y;
qn and qn/12 can be computed in one run on the 602.
Another way of relieving the situation is the following~
Define an auxiliary function
Z
= Y + h P/12 .
2
(25)
REl"ERENCES
1. L. F~INS'fl':IN and M. SCHW~RZSCE}LD, "Automatic Integration
of Lmear Second-Order DIfferential Equations by Means of
Punched Card Machines," Rev. Sci. htst., 12 (l941), pp. 405-8.
2. Ibid., p. 405.
DISCUSSION
[This paper and the following one by Dr. Paul Herget were discussed as a unit.]
Integration of the Differential Equation ~:f = P'P(r)
Using the Type 601 Multiplying Punch
PAUL
HERGET
Cincinnati Observatory
ferenced, it becomes possible to compute 6. iip and build
up the 6. ip and P columns by addition. It is necessary to
proceed step by step and by successive approximations.
In the solution of the original, more general, equation
it was possible to employ a Type 601 MUltiplying Punch
equipped with sign control and a net balance summary
counter. The board wiring may be illustrated schematically:
d 2P
THE EQU ATION dr~ = P·F (where F is a predetermined function of r), arises in the computation of
the wave equations for atoms in various stages of ionizaE),
tion. In that case it is necessary to replace F by (F
where E is such a constant as will cause P to vanish when
r approaches infinity. In practice the correct value, of E
must be determined by trials, and hence it is necessary to
run through this kind of a solution many times.
From the calculus of finite differences we have the following relationships:
+
d 2P
.
( 6.r)2 :1r~
=
(6.r)2 P (F
6. iip = f
+ 6. ii fJ12
Reading hrushes
+ E) = .f
- 6. iv f1240
+
Multiplier
Multiplicand
LH Counter
(5) = 0 Pickup
Sum. Counter
The ultimate objective of our computations is to obtain a
table of numerical values of P (r) which satisfy these two
conditions and which may be illustrated by the following
arrangement of the intermediate results:
r
P
0.0
1.000 000
0.1
0.927 288
f::.ip
f::.iip
f
-0.000 012
0.000 000
+0.000 915
0.000 927
-0.072 712
0.2
0.855 491
0.3
0.785 393
0.4
0.717 640
0.5
0.625 748
-0.071 797
-0.070 098
+0.001 699
0.G01 711
+0.002 345
0.002 356
-0.067 753
+0.002 861
0.002 871
-0.064 892
+0.003 255
0.003 255
f::.if
Punch
The first position in the punched field receives an X punch.
The group multiplier switch is wired oFF and ON at the
same time. The OFF switch permits the multiplier to reset
on every card. The column in which the X is punched is
wired to read as if it were the group mUltiplier master
card indication. This has the effect that when any card is
punched and thei1 fed through the machine a second time
it will be skipped out as a master card, without punching.
The field A is always cross footed into the LHC and it
transfers to the SC with sign control. The only multiplier
is (6.r) 2F, which is prepunched into a set of salmon cards.
The multiplicand is wired reversed from the positions
where the punched field is "reflected" through the center
of the card. All cards have their index numbers prepunched. The P's are manila cards, the 6. ip'S are green
cards, and the 6. iif112's are blue cards. These must be
obtained from a previous approximation and have the
above mentioned X punch.
f::.iif f::. iiif
-146
+927
+784
+645
+515
+393
-143
-139
-130
-122
-110
+ 3
+4
+ 9
+8
+12
This illustration represents the solution of the simplified
d 2P
equation (1;:2 = Pr, where P(O) = 1, and P (00) = 0; the
problem is to find dP / dr at r = 0 such that the condition
P ( 00) = 0 will be fulfilled.
The only numbers which can be entered directly into
the table are in the f column, when they are computed
according to the first equation above. vVhen these are di f-
39
40
SCI E N T I FIe
One cycle of operations' consists of the following (the
index i refers to the ith hori:;onial line in the numerical
table):
Card Color
Salmon (i)
Blue (i)
Green ( i - 0)
Green (i+ 0)
Manila (i)
Manila (i+ 1)
Salmon (i 1)
+
1\,[ achille Operation
Multiply F by P (i.)
Add 6,iif/12
Add 6,lP (i - 0)
Blank
Add P (i)
Blank
Blank and reversed
Result -in Sum. Counter
Punch f
NOll punch 6,liP (i)
Non punch 6,lP (i
0)
Punch b,iP (i +0)
Non punch P (i
1)
Punch P (i
1)
Punch P (i
1)
+
+
+
+
The operator has the cards of different colors piled
separateIybefore him, each pile in order of the index i.
On one cycle he performs the following sequence of
operations:
punched Xon the salmon card when it is fed in the reversed position. This automatically clears the counter at
the end of each cycle of operations.
I f we undertake to apply these 'principles to the solution
of the first order differential equation, (6r )~p
works out as follows:
r
P(i) == 'f(i) - D,if(i)/12
3. The top card from the blue pile is picked up and
held in one hand.
4. The first green card is allowed to fall into the
stacker.
5. The second green card is placed behind the blue
card being held in one hand.
6. The top card from the green pile is picked up and
placed behind the other cards being held in one
hand.
7. The first manila card is allowed to fall into the
stacker.
8. The second manila card is placed behind the other
cards being held in one hand.
9. The top card from the manila pile is picked up
and placed behind the other cards being held in
one hand.
10. The top card from the salmon pile is picked up
and placed in the 1'eversed, position behind the
other cards being held in one hand.
11. The last salmon card is placed (in the direct position) ahead of all the cards being held in one hand,
and
12. This deck is now placed in the feed hopper to begin
the next cycle.
The operator may be illiterate, so long as he is not color
blind! The work proceeds at the rate of thirty seconds per
step in the table, which is nearly the speed at which the
cards can pass through the machine.
The SC does not reset except under control of the class
selector. The selector transfer is obtained only from a pre-
+ 11 6 ii1f(i)/720
= f,
it
- ...
where it will be noted that all the quantities on the right
side are "on the line" in odd difference columns, so that
they are actually the means of the quantities on the half
lines above and below. Then
6 i P(i + 1/2)
1. The salmon card is allowed to fall into the stacker.
2. The blue card is allowed to fall into the stacker.
COM PUT A T 1·0 N
+
= 4f(i)
4f(i + 1) -
-
2~ 6
H
2~ 6 uf(i +1)
f(i)
+ 1;;o6 ivf(i) .. .
iV
+ 1!!0 6 f(i +1) .. .
Now, f may consist of the algebraic sum of any number
of cards, and if the higher order difference cards are
already available from a previous approximation,it is only
necessary to include one control card in each control group
and to use two counters in order to build up the table of
numerical values of P. All the cards representing quantities on the ith line are entered into both counters. The
control change causes an intermediate (progressive) total
This gives the value of the integral, P ( i), on the ith line.
As the next control group starts through the machine, the
first card is the control card, and this rolls the second
counter into the first counter, then causes a minor total
which clears the second counter. This enters all of
the quantities from the ith line which are needed for
6 iP(i + 1/2) into the first counter. Then the remaining
cards in that. control group enter the quantities from the
(i + 1) th line into both counters, as before, and the resultant progressive total in the first counter is P (i + 1).
DISCUSSION
Dr. Grosch: I would like to make the general remark
that Dr. Herget has a big point in using the human element
in his cycle. We can all make use of the tricks that Mr.
Bell and Dr. Herget have described. If you reduce the
number of 601 or 602 control panels that are kept wired
up by thirty or forty percent, you have effected a sub~
stantial saving. We did an optical calculation two years
ago which required twenty-eight· boards. I don't think
we had heard of the reversed card at that time. That
one idea would have released eight 601 control panels for
other work.
Dr. Stanley: I have one question. Notice that in Mr.
Lindberger's equation (22)
FO RUM
PRO C E E DIN G S
g(x) is a function of one variable only. The equation is
linear and of the second order. Now a more general equation of the second order in the normal forn) may be
written
2
d S = g(xJ S) S
dx 2
+ p(x)
This equation is non-linear. You might conceivably treat
it in the same manner as the speaker has suggested, except
that you would obviously run into the difficulty of .computing the quanties Yn. We might construct some sUltable
program beforehand and use it to estimate the Yn. I wonder
if either of you gentlemen have tried such a method. .
Dr. Thomas: This is exactly the thing that Hartree dId
in his so-called "Self-consistent Field" computations.
There are two ways you can do it. vVith the notation:
where
41
He used V and any numerical approximation to t/ln to get
V nJ then solved the differential equation to get t/lnJ these to
get new V and VnJ and so on until you come out with
what you put in. A somewhat different trick was one we
tried a few years ago. Instead of assuming Vn we assumed
V and put o/n to get Vn continuously as t/ln was being computed. I don't think you gain anything by that, except that
every answer you get is a solution of the differential
equation.
Dr. Caldwell: It might be possible to do that kind of
thing with the 602 provided the functions were not too
complicated.
Dr. Thomas: It was the double integral that I had in
mind. You can do two of them simultaneously on the
405 as well as the constant. You go all through an
integration to get preliminary values for t/ln. These must
be normalized. These integrals must be obtained to get an
"energy" to put in on the right-hand side before repeating
the integration.
Some Elementary Machine Problems in the
Sampling WJrk of the Census
A.
ROSS
ECKLER
Bureau of the Census
The second way in which we use sampling is to carty
out independent field surveys based upon a sample of the
population from which we can estimate the total population of the country and the population in various economic
and social groups.
The third way in which we use sampling is in connection with measuring or controlling the quality of statistical
operations. I will refer to each of these uses very briefly
in some of the applications I will mention.
First of all, I should like to refer to' an application of the
machines which is a very happy one. This use is in connection with drawing samples of blocks for certain types of
surveys. We want to determine certain blocks in which we
are going to collect information. We have put in punched
cards certain facts relating to each block in all of aUf cities.
That information, among other subjects, includes the
number of dwelling units, the number of stores, and the
number of various types of institutions. As we take our
population sample, we want· to select certain blocks in
which we will do our sampling.
We have determined that under many circumstances an
efficient procedure of drawing the sample of blocks is to
draw it on the basis of probability proportionate to size,
i.e., the number of dwelling units, or number of stores.
VvTe have been able to develop a procedure for selecting
blocks by the use of the Type 405 whereby we run through
the cards, accumulate, and select every nth dwelling unit.
The machine can be wired so that if in a very large block
there are two or more units which are to be included in
the sample, this fact will be indicated by the machine. If
any of you are interested in that, we would be glad to have
you write to inquire about the method.
Another area in which we have made use of sampling is
in connection with the processing of data. Weare particularly interested in the development of better equipment to
handle sample materials because it will give us a possibility
of increasing the use of sampling, thereby taking greater
advantage of the benefits it offers. We are anxious to extend the use of this tool as far as possible, and in certain
PER HAP S I have a definite advantage over most of
you in being able to recognize the value of this kind of
meeting. I do not come here as a mathematiCian nor as an
expert in machine accounting; so perhaps I am peculiarly
able to see the advantages of bringing together these two
types of people. In my opinion the International Business
Machines Corporation is to be commended for its vision
in making possible meetings of this kind. The advantages
for both groups are very great, and I have been much impressed with the gains from this sort of meeting even
though much of the material is highly technical.
Most of you are familiar with the long-run interest of
the Census in large scale accounting equipment. Weare
very proud of the fact that in the early years men like
Hollerith and Powers were employees of the Bureau of
the Census, and we have for many years used equipment
specially developed for our needs. We have used that as
well as very large quantities of the different types of IBM
equipment.
I shall speak primarily of our work in the field of sampling, which involves certain applications of equipment
somewhat different from what we get in our complete
tabulations, and which illustrates some areas in which the
present equipment fails to meet the requirements that we
would like to see met.
It is unnecessary to inform this group about the advantages of sampling. Most of you are familiar with the
theoretical work to a far greater degree than lam. You
doubtless know that through the application of sampling
we have been able to save very large sums of money in our
tabulating work. Moreover, we have been able to speed up
results so that we have been able to carry out many types
of detailed tabulations which would be far too expensive
to carry out on the basis of complete coverage.
There are several directions in which we apply sampling.
One· is the use of sampling to serve as a supplement to a
complete census, asking certain questions on a sample basis
only. In this way, we have been able to increase very
greatly the number of subjects covered.
42
FORUM
PROCEEDINGS
areas the mechanical equipment is a limiting factor. We
could go further with it if we had equipment which fitted
the needs more precisely than the present equipment does.
Just as we depend upon equipment to expand the use of
sampling, we also use sampling to improve the use of
equipment. We are carrying out a great many of our
processing operations on the basis of sample verification.
This takes place in a great many fields; one example is
our foreign trade statistics, which involve tabulations of
information on imports and exports by country of origin,
country of destination, etc. vVe have developed a system
of sample verification, which usually provides for a sample
of one card in fifty. We continue with that sample as long
as the operator is making fewer than sixteen errors per
400 cards sample verified. When she exceeds that rate of
err:or, we shift over to 100 per cent verification for a short
time. Then, when the evidence is available to show that the
person has come back down to a lower error rate, we shift
back to a five per cent sample and after a period of that,
if the rate still continues low, we go back to the two per
cent sample.
Through the use of that type of sample operation, we
have been able to have the verification of the work of
thirty operators handled by three. This achieves a very
considerable saving in the verificat~on operation, and still
provides control of the work so that we can be sure that
our error rate is under two per cent.
N ow I should like to mention the major applications of
sampling which take place in a number of fields in the
Bureau. In the field of current population surveys we are
carrying out samples on a monthly basis. Weare doing
somewhat the same kind of work in the field of business
statistics for retail and service firms, and generally similar
work in government statistics, where we collect employment data for state and local government units.
In the first field I mentioned, our current population
surveys, we interview a sample of about 25,000 households
once a month. vVe get information from them on the
number of people who are employed, the number unemployed, the hours worked, the occupation, industry, and
so forth. The households are selected by the use of area
sampling, a method probably familiar to most of the people
in this room. It is based upon units which are selected
from sixty-eight different sample areas scattered around
the country, scientifically determined so as to give a good
cross section of the country as a whole. We insist that all
of our samples have measurable accuracy; in other words,
that the design be such that we can determine the degree
of error in the results.
In this current survey of 25,000 households we estimate
that the figure on the total labor force will be within one
per cent of that which would be obtained from a complete
census nineteen times out of twenty. vVe achieve that very
43
high degree of accuracy partly by virtue of the fact we
have control totals for various groups to which we can
adjust the sample results. Obviously, in a sample survey
of this sort giving monthly information, speed is of greatest importance. These data are highly perishable and it is
important we make them available as rapidly as possible
because they are widely used. The information we get for
these 25,000 households is punched in about 65,000 cards
for individuals and those cards are weighted according to
the sampling ratios that were used. As each card gets a
weight which depends upon the age, sex, and residence
group of the population from which it is drawn, a considerable number of weights must be applied. There is
nothing about that job which can not be handled by
standard equipment. The difficulty is that it takes quite
a while to carry it out.
It is necessary first of all to sort the cards into these
sub-groups, to determine the total number in each, and
then to determine the weight which has to be applied to
each type of card. Even after considerable experience, we
found ourselves unable to do the whole job in less than
about thirty days, which meant that the data collected in
one month would appear fairly late in the next month.
After some experimental work we shifted over to a
Census-built machine, the unit counter, which has the
advantage of somewhat greater speed of operation; for
this particular kind of job, it can count in a considerable
number of categories-sixty-at the rate of 400 per
minute.
This procedure has certain disadvantages, however; in
order to use this machine we are forced to use a less precise weighting system than we could use on the 405. We
accomplish the weighting by classifying our cards in 144
different groups and then rejecting some cards by random
methods and duplicating some others by random methods
so as to get our results weighted according to the predetermined weights.
With that slight sacrifice of precision we were able to
speed up results a great deal and I believe, at the present
time, we use about fifteen days. I think it is pretty clear
we don't yet have ideal equipment for meeting this particular problem.
One other area of work I will mention briefly in closing.
In connection with our sample work, we attempt to establish a very careful measure of the degree of accuracy of
the results so that we have to compute large numbers of
variances and, as you well know, that involves calculating
very large numbers of sums of squares and sums of products. In fact, for the measurement of accuracy of just one
item, it is necessary to get the squares of more than three
/thousand numbers, to weight them by certain factors, and
then to combine them. It is possible to get these sums of
squares and Stul1S of prodl1cts through a rather compli-
44
cated series of operations, but the time required for that
is considerable. It is not a .very efficient procedure. Some
consideration has been given to the extent to whiCh some
of the new high-speed multipliers will meet the problem,
but study of the situation so far indicates that we are still
considerabl y handicapped in the direction of getting these
measures of the accuracy of the results of sample data.
There is a need for further development which will increase the use of sampling by making it possible to measure the accuracy of the results more speedily.
I have brought several types of exhibits in which some
of you may be interested, a number of pamphlets and
bulletins which show cases in which we have made direct
use of machine tabulation sheets for publication: some of
our housing reports, some of our foreign trade reports,
SCIENTIFIC
COMPUTATION
and a job we did for the Air Forces, on all of which we
printed the 405 sheets directly. If there are questions you
want to ask about them, feel free to write in. There are
also several copies of the booklet on work of the Census
Bureau-Fad-Finder for the Nation. If any of you want
a copy of that, I would be very glad to furnish it.
DISCUSSION
Mr. Tillitt: At one time the Bureau of the Census
turned out a little sheet called "Tab Tips." Is that distributed any more?
Mr. Eckler: I think I have the first forty. If anything of
that sort is being distributed now, it has not come to my
attention.
IBM Applications in Industrial Statistics
CUT H B E R T
C.
HURD
Carbide and Carbon Chemicals Corporation*
a linear relation between x and y. That is, we assume that
there exists a relation of the form
THE N U:NI B E R of IBM applications in Oak Ridge is
great, partly because of the variety of activities in which
Carbide and Carbon Chemicals Corporation engages there,
and partly because of the fundamental. import~n?e ~f
methods of probability and of mathemabcal statistics 111
the atomic energy field. Thus, on the one hand, Carbide
and Carbon, an Atomic Energy Commission contractor,
operates the Oak Ridge N ati.onal Laboratory, the Magnetic Separation Plant, and the Gaseous Diffusion Plant,
the first of these being devoted to fundamental research,
the second and third to development and to production.
On the other hand, atomic energy problems in whatever
state of development are problems which require the most
careful kind of statistical analysis both in the design of
experiments, in the interpretation of experimental results,
and in maintaining statistical quality control in various
production and measurement programs. \iVhen the size of
a problem, as measured by the number of measurements
to be made, has been great, we have f.ound IBM methods
advantageous and in some cases almost indispensable.
Originally I had thought to give a survey of the type of
statistical problems encountered in Oak Ridge and to describe the machine methods which are used in their solution. However, in view of the previous papers .of this
Forum and particularly in view of informal discussions
with other participants I have now decided to describe
only two types of problems and on one of these I should
like to ask your advice. These problems are: one, that of
curve fitting; two, calculating approximate solutions to
partial differential equations.
y=aX+f3 .
I f the assumption is now made that measurements o? x
can be made with perfect precision, that is, that x 1S a
fixed variate, we set x at the values xi' ~t·:!, . • . , Xn, say,
and measure the c.orresponding values of y: .Yl1 )'2' ... , Yn.
We realize that the set of n value pairs represents only a
sample of all possible value pairs and consequently we are
confronted with this problem in statistical inference:
On the basis of the 11, sample values (Xl' 3'1), (.r 2, Y2)'
... , (.rnJ )'n), required to estimate the parameters a and f3
of the equation above. N.ow, mathematical statistics is a
pure science in the sense that conclusions follow inescapably from assumptions; and, in certain fields, investigators
believe that a reasonable assumption for the present problem is that of normality of distributi.on with fixed variance
in y corresponding to each fixed value of x. More specifically it is assumed that if x is fixed at Xl' say, and the
corresponding y measured repeatedly then these y values
will be normally distributed about a mean value given by
a + f3 Xv and a variance of (J'2, say. Similarly repeated
measurements of y corresponding t.o fixing X at X 2 would
result in a normal distribution about a mean value of
a + f3 ~t'2 with the same variance, (J'2.
Under these conditions it is not difficult to show that the
optimum estimates ~, f3" of a and f3 are given as solutions
.of the simultaneous equations
THE PROBLEM OF CURVE FITTING
The Straight Line (One Fixed Variate)
The pr.oblem of curve fitting is as old as experimental
science and it is familiar to each of you in one connection
or another. As a beginning, suppose that we are able to
measure each of two variates x and y and that we postulate
I will not define optimum, but will only say that in the
case considered this procedure leads identically to the
least squares solution and the maximum likelihood solu-
*Since appointed Director of Applied Science, International Business Machines Corporation.
tion. As such the estimates ~ and f3" have the properties of
45
46
SCIENTIF1C
consistency, efficiency, and sufficiency. More.over they are
unbiased. 1
N ow it is clear that the requisite sums and sums of
products n, ~x, ~X2, ~Y, ~xY in the above equations can
be computed on a desk calculator. Also, it is well
known that the accounting machine provides, with the
use of progressive digiting, 2 an extremely fast meth.od
for computing these sums. A third method would be
to use the 602, punching individual products as the calculation proceeds. I will not attempt to describe completely
the conditions under which we ch.oose one of the three
methods described above. 'these general observations can
be made, however: because of the increased opportunity
for checking we generally prefer to use a punched card
method even when the amount of data is small; second,
since the calculation in the linear case is frequ'ently only
preliminary to later work, we are inclined to use either the
602 method or a combined 602 and 405 method. In this
way we calculate individual products and save them in the
card, and toward this end we have permanently wired 602
control panels for some of the curve fitting problems
which we encounter frequently.
Several Fixed Variates
I will now discuss a more general case of curve fitting
and in so doing I will indicate that in an important class of
statistical problems it is necessary not only to compute the
value of the inverse of a matrix but to compute explicitly
the elements .of the inverse matrix. We shall see, then, that
the Type 602 methods of equation solving which Mr. Bell
described are important in statistics and we wi11note also
that it is convenient to augment the matrix of coefficients
with the unit matrix in order that the elements .of the inverse matrix may be obtained explicitly.
One example of curve fitting in several dimensions arises
in the plastic industry. We denote by y the molecular
weight of a plastic. We denote the operating conditi.ons
of a production process as follows: temperature, Xl;
pressure, X 2 ; amount of agitation, X 3 ; time, X 4 ; amount of
stirring, Xo; amount of monomer,x s . We assume that
Xl> X 2 , • • • , Xs can be measured with perfect precisi.on
and that they can be fixed at prescribed values. We
then perform an experiment in which n sample values
(Yj; Xlj, X 2 j) ••• ,xsj), j = 1, 2, ... , n > 6, are obtained.
If a linear relati.on of the form
Y =
ao
+ al + ... + as X
·1'1
6
is postulated, then the problem becomes that of estimating
the parameters ao, aI' .•• , as'
Another example arises in the field of industrial medical
statistics in which it is desired to predict the size of the
dispensary staff and the amount of dispensary equipment
needed in a new plant on the basis of knowledge of characteristics of employees; age, sex, race, occupati.on, educa-
COMPUTATION
tion, salary, kind of chemicals worked with, etc. The
medical aspects of this problem are under the direction of
Dr. J. S. Felton, the mathematics under Dr. A. S. Householder, both of Oak RidgeN ational Laborat.ory, and the
calculations include the inversion of a matrix of order
fifty-five, these calculations now being performed on IBM
equipment under J. P. Kelly and B. Carter at the K-25
Plant.
In general, we suppose that we have fixed variates
.1'1' X 2 , • • • , Xk, in which, for convenience of notation,
we define Xl to be identically equal to one. I will not formulate the problem in detail but rather will refer you to an
excellent discussion of the problem by vVilks. 3 We make
certain assumptions about normality of distribution of repeated measurements of y corresponding to fixed values of
the %'s, and we assume that there is a relation of the form
y =
a1 Xl
+ a 2 X + ... + ak Xl.;
2
•
On the basis of a sample of size n
the maximum likelihood estimates and likewise the least
.
" "a 2 , ••• , ak
" 0f t he parameters are then
squares estImates
al>
given as s.olutions of the simultaneous equations
L.
fI,
+
a2
L.
fI,
Xlj X 2j
+ ... +ak
}=1
=
}=l
L.
fI,
j=l
L.
fI,
Xkj X
Ij
+
a2
akL. %2j Xlcj
j=l
L.
L.
fI,
j=l
al
Xlj Yj
j=l
fI,
+ ... +
L.
fI,
'xl} Xkj
=
%2j Yj
j=l
fI,
Xkj X 2j
+ ... +
ak
X1Zj
j=l
Computationally, the first problem is that of computing
the k(k + 3)/2 sums of products which are exhibited
above. Second, we wish to solve the equations themselves.
We thus arrive at estimates for the parameters which
are optimum in a certain sense but these estimates are
single estimates or point estimates. I f pressed for a single
estimate we would give ~l as an estimate of al> etc. But
statisticians prefer to give also a range of values for aI'
a range of values for a 2 , etc., and to associate this range
of values with a probability statement. For example, ~e
might choose a confidence level of 95 per cent and then
arrive at statements such as:
FORUM
47
PROCEEDINGS
accounting machine equipment applies also to the present
case.
vVe assert with 95 per cent confidence that
(a~,
aD
a2 lies in the interval (ag,
a~)
a1 lies in the interval
ak
ApPROXIMATING SOLUTIONS
To PARTIAL DIFFERENTIAL EQUATIONS
lies in the interval (aIL al) .
I do not intend to describe either the theory of confidence intervals as developed by Neyman, Pearson and
others or the theory of fiducial probability as developed
by Fisher. Rather, I will say only that, in order to arrive
at statements of this kind, it is essential to c.ompute explicitly the elements of the inverse of the matrix of the
coefficients of the equations exhibited above. In concluding
this section, let me say that all of us have rule of thumb
methods for deciding when to use hand methods and when
t.o use IBM methods for solving linear equations. These
rules must be modified in curve fitting problems because
we must not only solve the equations but also compute the
coefficients of the equations and we must compute explicitly the elements of the inverse matrix. Consequently,
IBM methods are efficient for a lower value of k (number
of equations) in the statistical problem than in the straightforward linear equation solving problem.
More Ge'neral Forms
It should be mentioned that the mathematical models
assumed in the previous sections can be thought of as
including many, but not all, of the situations which one
encounters in science and industry. For example, we note
that parabolic forms can be obtained from
by setting
Xl
=
1, .t'2 = x 2,
X3
=
xi, ... ,Xk
=
X£,;-1 •
However, there is another type of curve fitting problem
whose solution requires a more complicated procedure
than that given above but for which accounting machine
equipment can still be used advantageously. I refer to a
problem such as that of fitting an equation of the form
y
fi1%
= a1
e
+a
fi2%
2
e
fik%
+ ... + ak e
In this case the equations to be solved in order to minimize
the sum of the squares of the deviations of the observed
from the calculated are not linear. Hence the method
of successive approximation described, for example, by
Deming4 must be applied. But by this method we arrive at
a set of linear equations whose solution leads to the numbers by which to adjust our first appr.oximations. Consequently, what has been said above concerning the use of
I shall only describe a problem on which we are working
for Dr. A. S. Householder and Dr. B. Spinrad at the Oak
Ridge National Laborat.ory. I should then like to ask for
comments arising from your experience with problems of
this sort.
vVe have the Laplacian partial differential equation in
two dimensions with boundary conditions given on two
squares, one square being l.ocated within the other. Because of certain symmetry conditions we need to consider
the problem over only a triangular region. We then construct a lattice and replace the differential equation by a
difference equation which relates the value of the solution
at one point to the value of the solution at each of four
neighboring points. By punching cards on which are designated the coordinates of the lattice points and the boundary values 'we then proceed in a combination of the
operati.ons of collating, reproducing, and calculating. The
process is that of Gauss-Seidel and we find that, for one
thousand points, we can perform an iteration in about two
hours and a half.
Now the questi.on which I have is whether any of you
have a criterion as to when sufficient convergence has been
obtained. Do you continue to iterate until there is no
change in the final digit carried in your machine? Do you
prescribe in advance an upper bound to the sum of squares
of deviations from one iteration to the next? What is a
good criterion?
In concluding this section I should say that many problems in industry can be attacked either by curve fitting
methods or by analysis of variance methods. F. C. Uffelman and E. W. Bailey of the Y-12 Plant in Oak Ridge
have devel.oped efficient machine methods for applying
analysis of variance techniques. I am sure that l\fr. Bailey,
who is at the Forum, would be glad to discuss that particular problem with you.
REFERENCES
1. H. CRAMER, Mathematical Methods of Statistics (Princeton,
1946), chap. 32.
2. W.]. ECKERT, Puttched Card Methods of Scientific C011'lputatio11t
(Thomas J. Watson Astronomical Computing Bureau, Columbia University, 1940), p. 31.
3. S. S. WILKS, Mathem,atical Statistics (Princeton, 1944), chap~
VIII.
4. W. E. DEMING, Statistical Adjustment of Data (Wiley, 1943),
chap. IX.
DISCUSSION
Mr. Stanley: In regard to that last problem and the
criterion required, I notice that the problem seems to be
one which could be very easily adapted to the relaxation
48
method of Southwell. I might mention a paper which gives
a very good criterion in that case. 1 Offhand, I can't state
whether the criterion can be easily adapted to the GaussSeidel method, but in view of the close similarity between
the two methods it is quite likely that it can be.
In reference to Dr. Hurd's first problem: in that fitting
of exponentials, I suppose by using least squares, the idea
was to take a number of ordinates and, with the assumption that there were perhaps three terms in its series, to
fit a least squares system to them. You arrive at a cubic
equation and solve it for the three unknowns, whence the
required values. The alphas are relatively easily derived.
However, I wonder how you attack that problem when
the number of dimensions is very great-say twenty or
thirty. Then least square fitting becomes very laborious.
Dr. Hurd: Your first comment, as I understand it, is
about whether the Gauss-Seidel method is a relaxation
method. They are usually talked of as being one and the
same.
Dr. Thotlws: I should like to say that while many people
use the names indifferently, I think it is well to make the
distinction that if you iterate regularly you are using the
Gauss-Seidel method; if you always improve the worst
point you are using the relaxation method.
Dr. Hurd: In reply to Mr. Stanley's second question,
we wouldn't have ten ordinates. We would probably have
a thousand or perhaps two thousand. Actually the least
squares method is the way we would approach the problem.
Dr. T'ltkey: With regard to this problem of fitting the
exponentials, there is an old paper in the Vienna Academy
Proceedings,2 when people were first starting to untangle
radioactive chains, which I think probably isn't as well
known as it should be. It amounts to this, that if you take
your % values, as if you had taken a square network and
projected it down in a slanting manner, then you can solve
these in a particularly simple fashion.
Now where you have a thousand values, you ought to
be able to sort out one hundred which approximately meet
these conditions, and get a rather simple and accurate first
set of values.
Afr. Stanley: As for this relaxation method, I wonder if
you have had any experience in trying to adapt it to
punched card equipment.
As Dr. Thomas has pointed out, the relaxation method
is not well suited to this equipment, as it goes from one
point of the net to another in what might be considered
a haphazard manner. Possibly punched card machines can
be used in a satisfactory manner in conjunction with a
system of block or group relaxations.
Dr. Thomas: In regard to using relaxation methods, it
is more advantageous, at least for convergence, to go over
the network by hand, make the adjustments directly, and
then compute the residuals on the machine. That is the
SCIENTIFIC
COMPUTATION
main advantage of the relaxation method; adjusting the
worst points gives more rapid convergence.
Mrs. Rhodes: Did you try the harmonization scheme?
I have never tried it on the IBM equipment, but I have an
idea with the Type 604 we might save a great deal of time
using this scheme. Moskovitz 3 and Frocht 4 gave this
scheme for just such areas as Dr. Hurd described; it
yields a first approximation in no time at all. Emmons 5
gives a problem very similar to the one Dr. Hurd showed,
as an illustration of the superiority of the relaxation
method over the Liebmann transit method of solution. As
I recall, the times for the respective solutions were given
by him as 20 and 11 hours.
Dr. Grosch: To do single point relaxations, one can use
both feeds of a collator to search for the largest residual
in a deck of cards. That means it is possible to survey
eight residuals a second: faster than a human operator by
almost a whole order of magnitude. A twenty by twenty
mesh could be reviewed in two minutes. The catch is that
the machine does not skip right to the troublesome areas
of the mesh in the way the human eye does.
.lllr. Stanley: It really seems that, because of the time
consumed, a system of single point relaxation is out of the
question. Group relaxation might do for a part of the
problem.
Dr. Grosch: Unfortunately, group relaxations don't
work out too well in terms of standard equipment. I am
afraid it requires a sort of human thinking that the 604
and the collator won't do. The SSEC, perhaps, might.
Dr. Thomas: I would like to remark that just as in the
problem of ordinary differential equations we talked about
earlier, you can get a simple formula. The error term is in
the sixth order and actually computing is no more complicated. The same thing is also true of partial differential
equations in two or three or more dimensions; in two, if
you take an improvement formula that involves only the
four values at the corners of the square in addition, you can
get a formula which is accurate to error terms in the sixth
differential coefficients; it will enable you to get accurate
results with a much coarser mesh and will save an enormous amount of time where you have a large number of
points.
Mr. Bell: Getting the inverse of a sixty by sixty matrix
is a big problem any way you look at it.
REFERENCES
1. G. TEMPLE, "General Theory of Relaxation Methods Applied to
Linear Systems," Proc. RoyaISoc., A169 (1938), pp. 476-500.
2. F. AIGNER and L. Fr,AMM, "Analyse von Abklingungenskurven"
Sitzber. K. Akad. Wiss. Wien, 121 IIa (1912), pp. 2033-44.
'
3. D. MOSKOVITZ, "The Numerical Solution of Laplace's and Poisson's Equation," Q. Appl; Math., 2 (1944), pp. 148-63.
4. M. FROCHT, "Numerical Solution of Laplace's Equation in Composite Rectangular Areas," J. Appl~ Ph~,'s., 9 (1946), pp. 730-42.
5. H. W. EMMONS, "The Numerical Solution of Heat Conduction
Problems," ASME Trans.} 65 (1943), pp.607-15.
Some Engineering Applications of IBM Equipment
at the General Electric Company
FRANK
J.
MAGINNISS
General Electric Company
calculations could be made in a satisfactory length of
time by using IBM machines. Accordingly, a set of these
machines was installed, and for the past two or three years
every new turbine set which has been built at the General
Electric Company has been previously analyzed for the
location of the critical speeds of the shaft. Moreover, all
sets which have been built in the past and are now in
service are being analyzed. In all cases in which test results are available, they agree very closely with the results
of calculation.
Briefly the method of calculation is this: the shaft (and
by shaft is meant the shaft and the units on it) is assumed
to be made up of a definite number (about thirty) of concentrated masses. Equilibrium equations are written for
each of the masses for the forces and moments, on it due
to the reaction of adjacent masses and the centrifugal
force due to rotation about the unstressed axis of the set.
A speed of rotation is assumed, and starting at one end
the deflection, slope, shear, and moment are calculated for
succeeding sections until the far end of the shaft is
reached. Unless the speed assumed is exactly a critical
speed, it will not be possible to have both the shear and
moment zero at this point. I f the shear is made zero, there
will be a residual amount which may be positive or negative. When a number of such calculations as described
above has been made for a number of speeds of rotation
(at about 200 rpm intervals for a 3600 rpm machine) the
curve of residual moment vs. shaft speed will indicate by
its zeros the desired critical speeds.
W HEN I was asked to give this talk, it was suggested
that I describe one of the problems for which we had employed our IBM machines. Since we have two sets of these
machines being used for engineering calculations on which
some variety of problems has been worked out, I thought
it might be of interest to give a more comprehensive picture of what use we are making of these machines, without,
however, going into much detail for anyone application.
Our two groups of machines are operated by the Turbine Engineering Division and by the Central Station Engineering Divisions. The former set is used for the calculation of shaft critical speeds and for the solution of a
problem which will be described by Mr. Kraft in his talk.
The other set has been used for a greater variety of problems. I shall try to describe these problems (with the exception of Mr. Kraft's) to indicate why we felt IBM
machines would be useful to us at the General Electric
Company.
Critical Speeds of Turbine Shafts
One of the important problems facing the designer of
large turbine sets is that of determining the critical speeds
of the set. Critical speeds are those rotative speeds which
coincide with the natural frequencies of transverse oscillation of the shaft. A better set from a vibration standpoint
is obtained if the critical speeds are not too close to the
running speed.
In the past, these criticals were approximated by considering that each unit between bearings (the high pressure turbine, the low pressure turbine, and the generator)
was (l separate entity, and the lowest frequency critical
speed of each of these single spans was calculated under
this assumption. This method was used only because no
other method was available which would give results in a
practical length of time.
Recently, our Turbine Engineering Division decided
that a more accurate determination of the critical speeds
of the multiple-span shaft involving a large number of
Correlation of Data on Electrical Steel
Obviously a factor of ,great importance to the manufacturer of electrical equipment is the magnetic characteristics of the steel he receives from various fabricators.
Any significant departure from a standard may have serious consequences in its effect upon the performance of a
rotating machine or a transformer. It is therefore a question of some interest to the people in our laboratories
49
50
what factors are responsible for such deviations from
standard. Are they caused by small variations in the
amounts· of carbon, manganese, cobalt, silicon or other
elements which enter into the composition of the steel?
Or are they caused by slight differences in the heat treatment the steel may be given by different manufacturers or
by the same manufacturer at differe,nt times?
An attempt was made to answer the first of these questions by determining the correlation between the percentage of each of the elements present in the composition of
the steel with certain magnetic properties and by making a
multiple correlation between overall composition and the
magnetic properties. This particular study showed only a
small correlation between chemical composition and magnetic properties.
This problem afforded an excellent opportunity to use
the technique of progressive digiting on the accounting
machine.
Multiple Cond1utor Circuits
A recent trend in the field of transmission of electric
power has been in the direction of higher voltages. A
reason for this is the increased distance between power
source and load center. Tests up to 500 kv are being conducted to determine the feasibility of higher voltage transmission.
When high voltages are to be used, there appear to be-t
certain advantages in using more than one conductor per
phase. For example, the line inductive reactance will be
lowered as multiple conductors per phase are used and the
capacity susceptance increased. These result in an increased load-carrying capacity of the line for the same
voltage. It is of interest to determine the spacing and arrangement of conductors in one phase and between phases
which will prevent excessive corona loss. In order to determine these it is necessary to know the maximum value of
the voltage gradient which will occur at any of the conductors. This in turn depends on the instantaneous value
of the charge on the conductors.
The problem is to set up and solve the equations relating
the charges appearing on each of several conductors with
the voltages on the conductors. In the general case there
will be as many such algebraic equations as there are conductors. Moreover, in order to determine the maximum
gradient it will be necessary to know the charge distribution for several values of the voltages. For a two conductor
per phase, three phase line with two ground wires, we are
faced with eight algebraic equations. Each variation of
conductor size or spacing which is to be investigated will
require a new set of equations with different coefficients.
The coefficients are Maxwell's potential coefficients and
the equations may be readily set up. It then becomes a
SCIENTIFIC
COMPUTATION
question of matrix inversion· for the reason mentioned
above. An interesting and significant result of the study
was the fact that the maximum voltage gradient can actually be decreased for a two conductor per pha~e line over
that for a single conduetor of the same total area per
phase.
Incompressible Fluid Flow
For many field problems it is possible to obtain a very
close approximation to the stream lines and equipotential
lines by using some sort of network calculator such as the
AC Network Analyzer or the DC Calculating Board. As
an example, Concordia and Carter! determined the fluidflow pattern in a centrifugal impeller by this method.
However, most network devices are inherently limited in
range and accuracy (as are all simulator type calculators) .
I f, therefore, it is required to determine flow. lines to a
high degree of accuracy in a region of great curvature, a
digital method of solution is necessary.
The reported results of the problem of the centrifugal
impeller study just mentioned were further refined, using
IBM machines to improve their accuracy. A small region
around the trailing edge of the blade was subdivided into
a fine mesh of points, values of the stream function of the
boundaries of this region were assumed to be fixed at the
values obtained from the DC Board study, and values of
the stream function at interior points of the region were
improved by a relaxation method.
Incidentally, we believe that the only way to study a fine
mesh relaxation problem is to start out by setting it up
first on the AC or DC Network Analyzers using as many
units as possible to give reasonably close starting values.
If such a device is unavailable, it would be necessary to
start with a very coarse network and gradually increase
the number of points until the required fineness of the
mesh is obtained.
In addition to the work done on the centrifugal impeller,
a study is at present being carried out to determine the
flow lines around a turbine bucket. In principle the two
problems are identical.
Some other problems we have been or are at present
studying by means of our IBM machines are:
1. A study of the distribution of flux density in the
interior of a steel lamination under conditions of iron
saturation. Such a study could be valuable in predicting
transformer losses to be expected from a new steel. This
problem involved solving one of Maxwell's equations in
the steel subject to the condition that the total flux in the
steel be sinusoidal in time. This is a rather tricky condition
and difficult to satisfy since the best we could do was to
assume the flux density at a given point, for example at
the center of the lamination, as a function of time and by
FORUM
51
PROCEEDINGS
means of a difference equation determine the flux density
at successive points from the center out to the edge. If the
total flux was not sinusoidal in time, a different central
flux density was assumed and the procedure repeated.
Actually, of course, we did a number of cases at the same
time. This pr.oblem also led us to develop a method for
harmonic analysis for odd harmonics only.
2. A study to determine crystal structure is being carried on at the present time. The method we are using has
been described by Shaffer, Schomaker and Pauling. 2
3. Finally, we endeavored to solve a small but complicated circuit problem involving complex voltages, currents and impedances t.o which several different frequencies and values of load were applied. Our conclusion
from this study was that the ratio of the number of steps
involved in the solution to the number of cases considered
was so large as to make this particular problem unsuited
t.o IBM machine solution.
I think I may state our conclusions as to the use of these
machines as follows. Our set of machines on which the
critical speeds are being calculated is being used very efficiently. Although there are not many cases (that is, values
of frequency) for each turbine set studied, the successive
steps are identical, control panels are permanently set up,
the operators are th.oroughly familiar with the procedure,
and the work grinds out day after day.
The other set 0 f machines on which the remaining
problems I have mentioned have been worked out cannot,
obviously, be used so efficiently, since we must be prepared
to solve many different kinds .of problems. Of the types of
work we have carried out we believe the most fruitful use
of the machines compared to other means of performing
the same calculation has been in the field of solving algebraic equations or matrix inversion.
REFERENCES
1. C. CONCORDIA and G. K. CAR'l'ER, "DC Network Analyzer Determination of Fluid Flow Pattern in a Centrifugal Impeller"
ASME Trans.,69 (1947), pp. A113-18.
'
2. P. A. SHAF'F'ER, JR., V. SCHOMAKER, and L. PAULING, "The Use
of Punched Cards in Molecular Structure Determinations:
1. Crystal Structure Calculations," f. Chem. Phys., 14 (1946),
pp.648-58.
DISCUSSION
Mr. Stevenson: On the turbine blade problem, we have
worked a somewhat similar problem using basically the
method of Theodorsen and Garrick 1 for arbitrary airfoils.
By means of a prepared table, we are able to get pressure
distribution on an arbitrary airfoil in approximately fifteen minutes' work with a 405 only. That meth.od can be
further expanded to give you either the velocity potential
or any streamline you wanted to specify. I would be glad
to let you know about the method as I think it would reduce your machine time.
Dr. Caldwell: I w.ould like to mention two points. Mr.
Maginnis briefly touched on a subject which I believe has
considerable possibility in many of these jobs, and that is
the possibility of using a DC board to get a start on a
numerical solution. Likewise, one can use electronic analog
machines of lew precision but high speed to get a rough
approximation, and then jump off from there. I think that
is a very powerful type of combination. I f any of you are
doing this, and making it work, I would be interested in
knowing it.
Dr. Fe1'l1't: In some work which I may describe later at
these sessions, precisely that was done. We have converted
a gun director into an analog machine, and sets of suitable
approximate solutions for certain iterative digital processes are actually worked out in large numbers. It is, as
you say, a very powerful method and saves a tremendous
amount of time.
Dr. Caldzoell: An.other thing I would like to see somebody do some day is to work some fundamental improvement on the method of crystal structure analysis that has
been described. It isn't too happy, the way you have first
to draw out a lot of cards-picking them by hand is about
the only practical scheme-and run aro.und the place
juggling little groups of cards. Unless you have a very
conscientious clerical job done, you are going to get something mixed up. At least I haven't seen it work yet without
depending a great deal on the human element to keep
things straight.
Dr. Th01nas: I would like to say we have been developing much smaller files of cards for these various harmonic
analysis calculations, and we are proposing to do it entirely
by sorting, collating, gang punching and tabulating. It can
be done that way, and in many cases is quite as fast as
using many cards. If you have a really fast multiplier, such
as the 604, you could probably do it quicker with that, but
with only the 601 or even the 602 available, it is quicker
to do. the mUltiplying by progressive digiting.
REFERENCE
1. T. THEODORS:EN and 1. E. GARRICK, "General Potential Theory
of Airfoils of Arbitrary Shape," N ACA Tech. Report 452
(1940),
Planning Engineering Calculations
for IBM Equipment
BEN
FERBER
Consolidated V ultee Aircraft Corporation
Many of our jobs use relatively small quantities of
cards but many steps. For these jobs it is often a time
saver to wire a few simple control panels rather than a
single one that is very complicated. Of course, we can
afford to spend more time on the wiring of a complicated
board if it is to be a repetitive job, or if it is very large.
We have very good cooperation with our accounting
machine installation. Having our installations together
does have certain advantages. We can have our multiplying done on one machine and checked on another.
Although the actual procedure used to solve a given
problem would to a large extent be a function of the type
of machine that is available, the following approach to a
common structural problem may prove interesting to some
of you. Let
(MyIilJ - MilJK) x (JlfilJ Iy - My K) Y
IilJIy - K2
IilJIy - K2
where
fb = stress at any point whose coordinates are
x and )'J measured from· any pair of rectangular axes passing through the centroid
of the cross section
CON SOL I D ATE D V U L TEE started back in 1942
to use standard IBM equipment to help solve engineering
problems. The first use was on electrical load analysis. The
solution of electrical load problems continued until about
1944, when the use was extended to structural problems.
Because of the saving in time and money, the equipment
has been in use continuously since that time, and has been
applied to a wide variety of engineering problems.
We have found that it is a good idea to get an estimate
from the person whO' wants a job done as to how long it
would take to do it manually. Then we put it on IBM
machines if we are sure that we can at least hetter his
cost estimate; Another important point to consider is the
advantage of having the engineer cooperate and assist in
getting the job on the equipment. We must replan his. job
to suit the equipment. In fact, with his assistance we might
discover new techniques and provide him with a better
service than he requested or anticipated.
For a particular problem the first two columns of the
card are used to identify the problem. The card layout
form with the card code is filed for possible later use on
a new problem of the same type. The card code was very
helpful to us in looking back to find the cards that represented a particular job. It could also be used to sort cards
into sequence, or for control on various machines both as
to its digit and as to its zone punch.
After the war we had a large supply of control panels,
so we used a large number of permanent control pan~ls.
When the time came that we ran short of permanent
panels, we have gone to semi-permanent panels. For a
602 where we have a great deal of basic wiring that does
not change, we use a few manual wires. Perhaps some of
you have considered such a possibility for a 405 differencing control panel, using manual wires to make minor
changes.
Another use of manual wires combined with permanent
wires might be in wiring a panel for the first time. Testing
the panel' after a few steps and substituting permanent
wires for the manual wires enables us to see what we are
doing.
M ilJ and My = bending moment about the x-axis and
y-axis, respectively,
Ia: and Iy = moments of inertia of the cross-section
about the %and yaxes
K
= product
about the
of inertia of the cross-section
% and y axes.
If
1
C1
= IilJ1 y _
C2
Cg
C4
Dm
Dy
= C1 K
K'J
= C1IilJ
= C1I y
= xC2
= Y Cli
-
-
then
fb = MmDm
Y C4
Cg
%
+ J.l1yDy
•
What we did was simply to put % and y inside the
parentheses and take out the bending moments. Most of
the time we wish to hold the section constant and solve
52
FORUM
PROCEEDINGS
stress for various combinations of bending moments. If
we wish to continue and solve for shear flow, the basic
formula for shear flow q due to vertical and side shear 'in
an open section is
where
VilJ and V y
= shear
forces perpendicular to the x and y
axes
QIl) and Q'Y = static moments about the x and y axes
QIl) =
~AiYi
i
Q'Y = ~A'iXi
i
Ai
= area of the cross-section of the ith item
53
DISCUSSION
Dr. Caldwell: I especially appreciate Mr. Ferber's idea
of getting the customer to estimate how long it would take
to do the job. It is very common to find a customer arriv- .
ing at one of these installations that are supposed to do
mass production computations with his eyes bigger than
his stomach. He says, "Well, we can get a lot of work
done here for practically nothing," so he specifies a lot of
work. Later he says, "vVhy do you say this is going to
cost thousands of dollars? I could do it for a hundred
dollars. " You investigate a little further and find he is
going to do one per cent of the work for one hundred
dollars 1
Mr. Ferber: I might add that when a problem is worked
on standard IBM equipment with no extra devices other
than those normally available, one operator does about
twelve machine-hours of work per day.
Mr. Bell: The indirect saving resulting from the elimination of rework in production can be a very important
result of using IBM equipment for engineering design.
This follows directly from the ability to investigate many
design conditions quickly and completely, so that all design
parameters are well in hand before construction begins.
A Survey of the IBM Project at
Beech Aircraft (7orporation
JOHN
KINTAS
Beech Aircraft Corporation
THE PRESENT PAPER is intended to describe
some important problems being solved by the IBlYI installation in the engineering department of the Beech
Aircraft Corporation. Emphasis is placed .on the various
types of jobs processed by the IBM group for our engineering and sales departments. Consistent with the stated
objective of the Forum, particular emphasis is placed on
problems which arise frequently in structural engineering.
It is hoped that the ideas contained herein will help stimulate discussi.on and thereby foster a mutual exchange of
ideas.
particular emphasis will be placed on technical prablems
of aircraft engineering.
"VEIGHT CONTROL PROBLEMS
Of impartance is the current usage.of machine methods
in weight contral calculations. In much of this type of
work the machines are required to shoulder intermittent,
heavy work-loads, and have demonstrated appreciable savings in time and labor. In a general sense, our weight control graup determines the c.onditions of weight and balance
of the variaus Beech airplanes. In particular, this group
investigates the weights and center of gravity locations of
campasite airplanes far various loading canditions. Such
informati.on usually is obtained by cansidering the weights
and centroid lacatians of all fixed and movable items of
mass in the airplane. Punched card methods are readily
adapted to handle the large volumes of weight and balance calculations associated with this type .of wark.
INTRODUC1.'ION
At the present time the computing installation in Beech's
Engineering Department includes one each of the following International Business :Machines:
Type 601 Multiplier
Type 513 Reproducing Punch
Type 080 Sorter
Type 405 Accounting Machine
Type 552Interpreter
Type 031 Alphabetic Key Punch
SALES RESEARCH PROBLEMS
The need for such an installation was envisioned during
the war as a practical solution to some problems arising
from critical manpower shortages. It was considered that
the IBM group would function in the following ways:
Our computing installation was used recently in a sales
research pr.ogram ta determine potential markets far Beech
praducts. This pragram was essentially a statistical survey
.of numeraus parameters relating to sales patentials in
various geographical areas. The labar entailed in the
mathematical manipulatian of these parameters was considerable. However, many basic operations were repeated,
and therefore adapted to punched card meth.ods of solution. The machines relieved the sales research group of
many haurs of "donkey wark."
1. Accammodate certain types .of periodic and intermittent wark-loads of our engineering department.
2. Alleviate shortages and losses toO the armed forces of
skilled technical personnel.
3. Help relieve engineering personnel of routine calculations, thereby providing time far more important
duties.
STRUCTURAL PROBLEMS
To date our IBM group has handled a cansiderable
amount of w.ork relating to airplane weight cantrol, structural analyses, sales research, time-labor studies, field
service engineering, library records, and other engineering
problems. Some of these projects are discussed herein.
Consistent with the stated area of discussion, however,
Weare set up to handle the fallawing structural engineering problems by machine methads:
Three-dimensional flutter analyses .of aircraft
structures.
Harmonic analysis.
54
FORUM
55
PROCEEDINGS
SQlution of linear simultaneous equations.
Solution of complex matrix equations.
Computation of section properties of aircraft
structures.
For each of th~se problems we have provided our IBM
group with a master deck .of precoded cards, wiring diagrams, and a set of written instructions. The master decks
are considered as permanent equipment, that is, they are
not processed in the solution of a particular problem. For
specific problems, the master decks are reproduced to
obtain working decks of cards which are prQcessed in
accordance with the written instructions. The instruction
sheets we use at Beech avoid, wherever possible, reference
to technical significances of the steps being performed.
This divorcing of engineering aspects of a problem from
the required machine operatiQns enables the IB:M operator
to concentrate mainly on manipulation of the cards.
Flutter
The punched card method of flutter analysis we use is
based on the theory given in ACTR No. 4798. 1 Standard
procedures have been set up for the following basic flutter
modes:
1. Fixed surface bending vs. fixed surface torsiQn vs.
rotation of control surface (with or without geared
tab).
.
2. Fixed surface bending vs. fixed surface torsion vs.
airplane roll or vertical translation.
3. Fixed surface bending vs. fixed surface torsion vs.
rotation .of control surface (with or without geared
tab) vs. airplane roll or vertical translation.
These flutter modes are solved entirely by the machines,
except that a few manual operations are required during
the final stages of solution.
Usually, flutter analyses are conducted to determine the
critical flutter speeds and frequencies and the associated
values of damping coefficients. In SQme cases, additional
information such as mode shapes and amplitude ratios of
the component degrees of freedom also may be required.
The problem of determining these items may be resolved
into the two rather distinct phases of formulating and
solving the stability determinant. We formulate the determinant by straightforward operations on matrices, then
solve the determinant by trial-and-error iteration. Some
discussion on these important steps is considered desirable.
For any mode of flutter, the stability determinant is
composed of complex elements. The numerical value of
each element may be determined by evaluating and summing up certain aerodynamical and mechanical integrals.
Careful study has shown that:
1. Each of these integrals may be expreSsed as a summation of products of finite quantities. This would
be equivalent to considering that the wing is divided
into a finite number .of chordwise strips.
2. Multiplication and summation of these products may
be accomplished readily using methods of matrix
algebra.
3. The llumerical value of each element of any stability
determinant may be expressed as the product of four
matrices.
On the basis of these observations, we form stability determinants entirely by machine·methQds.
As previously mentioned, we solve the flutter determinant by trial-and-error iteration. It has been determined
that this process will converge at a practicaJ rate, since the
preponderant elements usually lie on the principal diagonal.
In most cases the iteration stabilizes satisfactorily in two
or three trials. However, in some cases fQur or more trials
may be required. In general the iteration process is carried
out in the following way:
1. First a trial value of w is substituted in all but one of
the elements along tli.e principal diagonal.
2. The determinant is reduced to the third order, if
necessary, then expanded to obtain a linear equatiQn
in one unknown.
3. The linear equation is solved to obtain the second
trial value 0 f w.
4. Steps 1, 2, and 3 are repeated until the trial value
of w agrees with the solution.
5. The entire process must be repeated a number of
times equal to the number of degrees of freedom
considered.
Practically all of the labor of solving the stability determinant is done by the machines. However, some manual
operations are required to estimate initial trial values of
and calculate flutter speeds, frequencies and damping coefficients at the final stage of solution.
Estimates of the time required for the solution of a
three-degrees-of-freedom flutter mode by manual methods
as compared with our punched card method, based on the
assumption that five values of (v/bw) are investigated f.or
each mode, are as follows:
(J)
IBM
lIfanual
Time
Time
Operation
35 hours
Computation of
elements of the
stability determinant
30 hours
26 hours
Solution of the
stability determinant
44 hours
61 hours/mode
(1 operator)
Tot~l
74 hours/mode
(1 person)
time
56
SCIENTIFIC
Wave/orm. Analysis
Complex periodic waveforms frequently occur on vibration records of structural investigations such as flight testing, fatigue testing and vibration tests of power plant
installations. It is impossible to analyze many of these
waveforms by ordinary inspection methods. However, any
complex periodic wave may be represented by the superposition of a number of simple .sine and cosine waves. It
might also be mentioned that aperiodic curves, continuous
in a finite interval, also can be represented by assuming
that the given curve represents a single cycle of variation.
The problem of waveform analysis is primarily that of
determining the amplitudes and frequencies of the sine
and cosine components present in the synthesized wave.
At any point along the reference axis, the ordinate of a
composite waveform is equal to the sum of the ordinates
of the component harmonics.
At Beech, we have expanded the Fourier series to obtain
the general equations corresponding to five, eleven, twentythree and forty-seven harmonics. Particular solutions may
be obtained by substituting into these expansions the nu- '
merical values of the ordinates of a given curve.
We have transferred the trigono.metric. constants of
these expansions into. a master deck of 4704 coded cards.
Solution of a, specific problem may be obtained by punching into a working deck (reproduction of the master deck)
the measured values of the ordinates. The cards are then
pro.cessed in accordance with standard instructions. Solution is accomplished almost entirely by th~ machines; some
divisions and extractions of square roots must be performed manually.
A comparison of time required by manual and punched
card techniques for waveform analysis is of interest.
Numberof
Ordinates
to Curve
Corresponding
Number of
Harmonics
r~
5
11
23
47
24
48
96
Est. llfallual
Time_
Hour{
2.0
7.5
28.0
110.0
Est. Machine
TimeHours,
1.0
2.5
7.5
26.0
COMPUTATION
of coefficients is operated on by rows and by columns until
all terms below the principal diagonal are zeros and each
term along the diagonal is unity. During operations on
rows, the· column matrix on the right side of the equation
is also modified. It is clear that the value of the nth unknown is immediately given on completion of these operations. The value of the nth unknown then is employed in
a back-tracking process to determine the (11-1) th unknown, and so on.
It has been determined that time saved by machine
methods over manual methods increases appreciably with
increasing numbers of equations and unknowns. This may
be seen from the following comparisOli:
Number 0/
Equations
and Unknowns
10
30,
50
ApproJ;imate
IBM
by Decimals
TI~me
for Soluti011i-7ManHouys
IBM
by Powers
Manual
Approximately Same
50
39
78
148
8.5
74
418
Matrz'J,' Manipulation
Weare equipped to handle certain types of matrix
equations by the Kimballmethod. 2 It is particularly useful
,in performing the basic operations of matrix algebra on
matrices with complex elements. The method is also useful
in manipulating matrices with variable elements when the
solutions are approximately known or when the preponderant elements lie along the principal diagonal. This latter
feature is a reasonable guarantee that trial-and-error iteration will converge at a practical rate.
Other Projects
At the present time we are investigating the feasibility
of solving the following structural problems by punched
card methods:
1. Spanwise airload distribution for monoplane wings.
2. Natural torsional frequencies of crank..,systems using
Holzer's technique.
3. Natural uncoupled frequencies of beams using Stodola's iteration procedure.
4. Analysis of shear lag in aircraft structures.
Linear Simultaneous Equations
Our IBM group now is set up to solve systems of fifty
linear equations in fifty unknowns. While this number
accommodates our preseilt needs, it can readily be expanded to any practicallil11it.
In general, we utilize a modified Gauss method. Here,
the linear equations are converted into matrix fo.rm. The
matrix equation, on the left side, contains a square matrix
of constant coefficients postmultiplied by a column matrix
of the unknowns. The right-hand side of the equation has
only one column matrix of constants. The square matrix
TIME ASPECTS
The time estimates previously given for several structural problems were based on actual performances. They
were determined by solving given problems manually and
by machine. Now the approximate average rates of, our
machines are as follows:
Multiplier-15 cards per minute
Reproducer-lOO cards per minute
Sorter-450 cards per minute
Accounting Machine-' 80 cards per minute (detail print)
150 cards per minute (group print)
Interpreter-60 cards per minute
FORUM
PROCEEDINGS
Processing time for most problems can definitely be improved through usage of faster calculating punches. We
probably will consider faster machines when the need for
more rapid processing becomes manifest.
RECOMMENDA'tIONS
Usually, punched card procedures can be adapted to a
given engineering problem in a number of ways. From
time and labor standpoints some procedures will be more
efficient than others. Among other factors, a determination
of the optimum procedure depends on a knowledge of the
full capabilities of the machines available for our use. For
this knowledge we rely to an appreciable extent on the
Wichita staff of International Business Machines Corporation. We have always found them highly cooperative.
However, in some cases they were unable to provide us
with enough information on specialized capabilities of the
machines, particularly our Type 601 Multiplier.
We recommend that local IBM offices be provided with
up-to-date information on the full computing capabilities
of the machines in their region. It may be possible to establish, on a current basis, the flow of such information from
the various IBM research laboratories and computing
centers to branch offices. This would help people like ourselves to realize the maximum potential utility of IBM
installations and avoid needless duplication of effort.
REFERENCES
1. B. SMILG and L. S. WASSERMAN, "Application of Three-Dimensional Flutter Theory to Aircraft Structures,"Materiel Division, Air Corps, ACTR 4798 (1942).
2. E. KIMBALL, JR., "A Fundamental Punched Card Method for
Technical Computations," Bureau of the Census, U. S. Department of Commerce (undated).
DISCUSSION
Dr. Eckert: I might point out that the use of IBM machines for technical computation grew very slowly for a
number of years. It was very easy for two or three people
to keep in touch with each other. Within the last two or
three years, there has been such a sudden cloudburst that
to have people everywhere supplied at the right time with
the right information is a little difficult. Steps are' being
taken; IBM has now, among other things, special representatives in the Sales Department who understand what
you are trying to accomplish. They know what is available
in IBM, and they are at your call.
With respect to the local manager, he also has a very
tough assignment when you ask him for methods ,of doing
things he has never heard of. That gap has to be bridged,
and I am sure it will be in a very short time.
There is still another way which is open at the moment,
and that is to call on us at the Watson Laboratory or write
us a letter. Weare always glad to hear from you. I think
about half of the people in this room have already done
that.
57
Mr. Kintas: You probably know Mr. K~mball of the
IBM office in Dallas. Would it be possible for him, or
others like him, to circulate periodically as good-will ambassadors? A person like that could not only convey information fro111 your laboratory to our installations, but
could pick up ideas from us to be passed on. Such representatives would be very valuable.
Dr. Eckert: That is the intention. The difficulty is to get
the right man, get him trained in a very difficult field, and
get him to you.
Dr. Kortz: Why not send Dr. Grosch?
J.\:lr. Schroedel: 1\ir. Kimball, Dr. Grosch, and I participated in a rather successful experiment along those
lines recently. Lectures were given at the Cornell Aeronautical Laboratory in Buffalo, and people from other
computing groups in the vicinity also attended. We learned
a good deal from the meeting. But there are fifty or sixty
installations in various parts of the country, and it is not
easy for all of us to visit every locality and work with you
long enough to really make a contribution. We hope to
have more technically trained representatives in the Sales
Department as time goes on.
Mr. Bisch: Several points, which could stand some comment and emphasis, were picked out of the very interesting talks of our aircraft representatives, Messrs. Ferber,
Bell and Kintas. Instead of taking one point at a time, it
seems more constructive and brief to present all my comments as a whole.
Many engineering organizations are looking today for
powerful means of calculation. What the Engineering Department of North' American did four years ago was to
look to IBM for the main answer to this need, after a
general survey of the field. Fast progress was desired,
and to this end production methods were used in the creation and the development of our Engineering Section of
accounting machines.
Such methods call mainly for extensive specialization
and perfect coordination. To be concrete, we select the
structures section, although the following would be true
for the aerodynamics section as well.
An engineer with a long acquaintance with problems of
structures, company policy, and organization methods was
asked to select the problems and among their various solutions those which result in speed and efficient use of the
IBM machines. This engineer familiarized himself with
the various functions of the IBM equipment and the pattern of calculation most suitable for it, but he made no
attempt to learn how to operate it. As a result, he was able
to increase the contribution of the machines and he
promptly reached the conclusion that for a given standard
problem, the machines can do everything from the punching of the initial data to the final report printing. It is
58
important to remark that he was the only engineer in direct
contact with the accounting machine section.
On the other hand, an IBM operator with several years
of experience and a mathematical knowledge equivalent to
a master's degree was assigned to program for the machines the problems offered by the structural engineer, to
propose profitable changes in the mathematical processes,
to suggest further use of the equipment, and to select the
type and the number of IBM machines. In no case did
he Concern himself with the engineering aspects of the
problems.
When this work was sufficiently under way to call for
more than two operators, an experienced operator was
selected as supervisor of this group. This group is an independent section of the main accounting section with its
own machines, and it derives many obvious advantages
such as readily available service for the machines and the
incidental facilities of a large. installation, by not being
separate from the main IBM body. The two men at its
head, the mathematician and the supervisor, form a perfect team for the dual purposes of large volume and continuous improvement.
The Engineering Section is therefore made of two parts:
the engineering part and the accounting machine part,
through which a perfect coordination is possible by a oneman contact. It is now opportune to detail further the
duties and achievements of each component.
The engineering part, which is also assigned to seek
new technical and experimental solutions of aircraft engineering problems, to write reports on the algebraic and
tabular solutions of problems and to conduct experimental
work, keeps an accounting of the approximate number of
arithmetical operations required by every job sent to the
IBM section. This number is easily arrived at for standard
problems, using simple algebraic formulas. On the other
hand, an extensive survey has shown that an engineer can
perform an average of one hundred arithmetical operations for each remunerated hour. As shown by two years
of coverage, an IBl\1 operator performs on the average of
one thousand such operations per hour; therefore his
speed is ten times greater.
Finally, last year's operations show an average of seven
thousand engineer-hours per month performed by an IBM
Section of four operators, which would cost a minimum
of twenty-one thousand dollars per month, if performed
by engineers, against an over-all cost of six thousand dollars by the IBM group, thus effecting a saving of fifteen
thousand dollars per month. This approximate saving of
seventy-five per cent which was evidenced at the very beginning was the real selling point to our management.
However, the engineering management has become conscious of other less tangible although equally important
advantages. They are:
SCIENTIFIC
COMPUTATION
Unprecedented dependability from the point of view
of time and accuracy.
An average speed ratio of ten resulting also in cutting
down waste in fabrication, as peinted out by
Mr. Bell.
Availability of solutions of problems, previously
prohibitive on account of cost and time.
I should like to conclude by giving a list of several questions frequently asked, and brief answers.
What is the accuracy generally 'used? Eight, sometimes
ten, significant figures.
T¥hy would slide rule accttracy not be a factor in th4
selection! Because in extensive calculations,even when
performed in tabular form and on desk machines, there is
a need for mathematical and engineering checks, which
would be meaningless with slide rule accuracy. In some
problems, such as the solution of systems of -linear equations, more accuracy sometimes means the only chance of
getting a correct answer.
~Vhy
not ask every engineer the esti'mate of engineering
ti'N'te, as practiced by Mr. Ferber? Because it is easier to
gather in one room all the supervisors and experienced
stress men and have them propose once and for. all a sober
rounded-off average.
Why not teach engineers how to operate the machines!
Although they are all welcome to visit the installations
and to see the machines performing, they only are interested in the relief in their task provided by use of the
equipment. Maximum efficiency is obtained by showing
them concrete types of calculations suitable or unsuitable
to these machines; even this is not necessary for the standard problems. Finally, all the past experience being available to the contact engineer, it can be poured into every
new problem.
What is the criterion for a solution by the l11a.chines!
The IBM section can solve a mathematical problem in any
way selected for efficiency, provided the final answers
would coincide with the answer found by engineering
methods.
Is the most efficient method reached the first time? Far
from it. First of all, any substantial increase of efficiency
is worthwhile at any time because the gain will be repeated
many times and will outweigh the cost of the changes.
Moreover, the survey of the first numerical application of
a method always suggests other improvements.
Wouldn't the IBjI! section soon be confttsed by a maze
of various methods? On the contrary; thanks to the specialization of the engineering group, it has been possible to
restrict them to a few general solutions from which many
problems can be derived as simple cases without change of
FORUM
PROCEEDINGS
procedure. For instance, the bending and shear stress distribution of shell structures can be obtained by the same
process, regardless of the type of cross-section and the
part of the airplane. Several obvious advantages are attached to such generalization.
What types of solutions are preferred, the type which
is short but inefficient on the machines, or the long but
efficient type? The first type is generally adapted to desk
machines. Our experience shows that the last type is preferable, because it taxes the operator less, and usually the
answer is reached sooner.
What aeronautical problems were found to be solved
efficiently!
1. Conventional and elastic determinations of stresses
and rigidity of shells.
2. Weight balance, static and dynamic.
3. Airplane performances.
4. Determination of design loads, magnitude and distributions.
S. General problem of arch analysis, stresses and deflections.
6. Redundant structure analysis.
What mathematical problems?
1. Solution, by iteration methods, of characteristic
equations such as in vibration problems.
2. Direct solution of systems of linear equations up to
seventieth order.
3. Solution, by Galerkin's method, of systems of linear
differential equations with variable coefficients.
What is our trend with respect to wiring control panels?
It may be that the nature of our work calls for flexibility.
In any event, more efficiency is attained by breaking the
process of solution into simple steps which are made with
standard control panels, than with special panels aimed at
shortening such processes.
59
Are we trying to use a minimum of caras! In general,
one card is used for one machine operation. The cost and
inconvenience of more cards is outweighed by the advantages of simpler procedures.
vVhat is the fut'ltre of such machines in aeronautical
engineering? As long as this industry secures experimental
contracts an increasing need for their use will exist.
Has use of IBM machines contributed to the progress
of aeronautical engineering? Very much. Experimentally,
thousands of stresses are daily recorded directly in cards
from a specially wired stress rosette recorder. Technically,
modern airplane structures can now be correctly analyzed
at reasonable cost, whereas before the advent of the machines, this would have been impossible. M-oreovcr, the
mentioned recorded stresses can be rapidly transformed,
through matrix calculations, into information valuable for
future designs.
Could the progress in 'ltse of IBM equipment for aeronautical engineering be readily used for other engineering
specialties ? For problems with identical mathematical
equations the answer is obviously yes. As to engineering
problems such as all general structure calculations (like
shell and arch analysis), their IBM solution can readily be
used in civil engineering, when presented under the general
form which is desirable in aircraft.
Finally, I might remark that none of our existing or
planned reports mentions the use or the wiring of the control panels.
Dr. Kortz.: Do you have those figures available? I have
to work a little on administration, too, in that respect.
Mr. Bisch: Every month we add to our little book several pages about new jobs; I can show you that. We make
entries for every job we send to the IBM machines which
involves at least a hundred man hours of engineering
time. Many jobs are in the thousands of hours. We keep
a close record all the time, and every two or three months
I send in a progress report.
Aerodynamic Lattice Calculations
Using Punched Cards
HANS
KRAFT
General Electric Company
W HAT I am going to present to you is by no means a
finished product. It is riot ev·en something I am exceedingly
pr.oud of. I would like to show it to you in some detail,
and hope that I will receive from you some criticisms and
help on how we could have this computation done in a
much shorter time than it takes us now. It is possible that
we are on the wrong track entirely.
The problem of the turbine engineer with all its complications is essentially as shown in Figure 1.
modern turbine-and it performs very well already-will
be made when, and only when, we are able to follow by
calculation the flow through this nozzle and bucket combination with the buckets moving at high speed.
This means that we have to compute a flow through a
row of nozzle profiles. In aerodynamic language such a
row of equally spaced profiles is called a lattice. We will,
in addition, need to kn.ow the flow through the bucket lattice. Furthermore, there is an interaction between these
nozzles and buckets. This interaction appears as a time
variation. As the buckets move past the nozzles, different
configurations of the available flow space result.
\Ve cannot rely much on the well-known approximati.on
of the flow by that of an incompressible fluid. Our velocities are such that we always have to consider the fluid
as compressible. Thus, we must first of all learn to compute compressible flow through a stationary, two-dimensionallattice ; later on we must study interference between
the two lattices as one passes by the other. Theoretically,
we think We know more or less how to handle the problem,
although as far as actual computation is c.oncerned, we
still have a long distance to go.
I should like to discuss some of the initial work which
we have done to describe a simple flow through a row of
buckets. It has been performed for flow of an incompressible fluid, but was done in a manner identical to that to
be followed to give us the compressible counterpart of
this incompressible calculation. The compressible computation still awaits the completion of a set .of input functions before it can be performed.
To solve the incompressible problem Laplace's equation
must be solved. We do not attempt, however, to solve a
boundary value problem. \Ve try to learn to build up profiles from given functions and accept the resulting shape
if it seems to be one which we actually do want, i.e., a
shape which will perform well.
We use the representation of a flux function t/t given as
a function over a field with the coordinates..X' and y. In the
compressible case we will not have this simple Laplace's
Momentary
Streamlines
Relative
Bucket Velocity
Nozzle Velocity
FIGURE
1
Most of you know that a turbine is essentially a windmill. A very fast flow .of steam or gas issues from stationary passages which are called nozzles. The moving blades
we of the General Electric Company call buckets. The
steam flow streaming from the nozzles passes between
them and is deflected. This process generates mechanical
power which is removed by the rotating shaft.
We have had a long history of experimentation. We
have experimented very intensively since 1920. We would
like to do some theoretical computations in addition. We
feel that we are somewhat against a blank wall with only
experimental approach. It is my own honest, private opinion that further improvement in the performance of the
60
FORUM
61
PROCEEDINGS
equation to deal with. Instead we must solve a rather forbidding non-linear equation if the problem is written in
terms of the physical coordinates % and y. To evade this
situation the equation is written in other variables, those
of the "hodograph" plane.
Figure 2 shows a streamline in the %-y plane. At any
point along any streamline there exists a velocity vector
given by magnitude and direction. The sequence of these
velocities along a streamline can be represented ina coordinate system such that the ends of vectors existing
along a streamline are connected. The result is a map 0 f
the true streamline. The map of all streamlines so pictured
is called the hodograph representation. It describes the
flow as well as does the original picture.
In the case of incompressible flow Laplace's equation
describes also the field in this hodograph plane. In the case
of compressible flow the differential equation applying to
the hodograph plane is linear. The important consequence
is that in this plane, solutions can be superimposed.
In our actual computations we are not using the true
hodograph plane. \iVe use the logarithm of the vector and
t
v
VORTE)C.
1
SOU~CE
4
o
A- ()
Physical Plane
PATTER.N
WITH SINGULARITIES
o
FIGURE
Hodograph
FIGURt
2
3
thus have a Cartesian system. Figure 3 depicts these representations for the case of relative flow through a row of
turbine buckets. In the physical ~t"-y field the flow comes
from infinity to the left and disappears to infinity at the
right. In the hodograph map all flow must come from the
upper singularity within the closed curve and disappears
into the lower singular point. A number of singularities
correctly placed outside the closed curve is needed to generate this pattern. These singularities are not shown here.
The streamlines within the closed curve describe all flow
62
SCIENTIFIC
through the bucket lattice. The two saddle points· at the
left infinity of the logarithmic hodograph represent the
entrance and exit of the buckets.
The object of the computation is to generate closed figures of this topology. They are generated by a vortex
source representing upstream infinity and a vortex sink
representing downstream infinity, combined with logarithmic singularities located outside the closed curve. Their
location and strength distribution will ultimately determine the contour of the profiles.
The mathematical background of the method is shown
for the more complicated case of compressible flow. The
basic ingredients for the final differential equation are
rather well known: continuity, irrotationality, gas law (in
this case for isentropic flow). Flux function", and potential cp appear as mathematical tools.
Equations (6) show the important conversion from hodograph variables to the physical flow picture. The linear
hodograph differential equation is Equation (8).
The linear fields so represented can be considered as
arrived at by the intelligent addition of logarithmic singularities. Additional non-logarithmic singularities may also
be added. In general, the desired results can be obtained
from logarithmic singularities alone. So much for compressible flow.
Here we are considering the far simpler case of incompressible flow. It serves as a pilot process for the really
desired more complicated compressible case. The difference appears in the computation of the logarithmic singularities which serve as input functions. This computation
which is desperately difficult in the compressible case is
much simpler for incompressible flow.
The only input function needed is given by the equation
Phys£cal Lam!s
+ a(pv~ = 0
Continuity
aept{,)
ax
Irrotationality
au _ av = 0
ay
Gas Law (isentropic) Py
P
ay
F
(1)
(2)
ax
= constant
(3)
Mathel1wtical Tools
=
It
Potential
Stream Function -
acp
ax
,v =
Pl V =
a",
ay
~"',
(JX
(4)
a",
Pltt = ay
(5)
H odograph Coordinates
= dcp cos 0 -
dx
~u
d", sin 0
Pivu
dcp .
d",
dy = - sm 0 + - cos 0
W
(6)
PlVU
dcp = ds
~u
Relations Betvueen cp and '"
avu
Pl
_W2)~
a",
a
ao
ocp
1
0'"
ocp = _ ~(1
ao = -
2
w
(7)
Pl ~u OVU
Differential Equation in H odograph Variables
1 -.,M2 02~
P-I
dO-
+ ~~(W at/!) = 0
PI dVU
PI dW
(8)
COMPUTATION
= cp + i", = In
(1n vu + iO) ,
(9)
where
(10)
tan 0 = viu .
( 11)
F is the logarithmic singularity called a source. Multiplied
by i it is called a vortex. This is well known to everybody
reasonably familiar with conformal mapping. The source
in the hodograph plane represents the axial flow component emanating from infinity of the .r-y plane. The vortex
furnishes the tangential component of the flow jlt infinity.
Added together in proper proportion, i.e., with one multiplied by the correct factor to depict the flow vector at
infinity, they furnish in the hodograph plane a map of the
physical infinity conditions. What has been said here for
the upstream infinity repeats itself at downstream infinity.
Here the source is negative, in other words, a sink.
We must realize that what is desired in the end is a
picture of the flow in physical space. The conversion from
the hodograph map to physical flow is linear, as is clear
from Equation (6). Hence x and y values can also be
superimposed. If they are known over the whole field of
the logarithmic singularity they can be added in exactly
the same manner as can the cp jlnd t/! values themselves.
In other words, if a field is composed of CP1 and o/H
CP2 and "'2; its physical coordinates are added from the
X 1 ,X2 and YVY2 of the individual singularities.
We calculated this x,y field for one singularity, but we
did not obtain very good accuracy (Figure 4). The equation for this field is shown on the figure. One-half of this
field has been computed to great detail and accuracy by the
National Bureau of Standards Computation Laboratory in
N ew York. The field as shown is for a vortex. To find
x and Y for a source, the relations are
(12)
Ie
)(
~
...
~
G)
:0
'iii
-4
«I)
...
G)
~
''''''
E
0
~ ..,I~
'--.:... i
(,)
....=
0
.. at
G)
...0
ii:
(.)
s::
>-
...
't:J
0
»
~
«I)
G)
0
.!:
cD
.J
~
~
Ie
II
-<
CJ
\
\
...E
\
I\V'"
\
\
N
"0
Ie
CJ
\
~\
\
~
_\
CJ
"
01
0
..1'<";......
't:J
"'i\
....\. \
-.!.!.\\
~,
*' ',\
'"..........
0
X :::E:
\
\
'\
\
\
,\
'\
I
,
I
- - $t:~ii,( - - ....
I
\
\
\
/
/
I
\
,"
\
/
I
I
\
\
\
"\
"'\ ,
/
I
\
\
/
,
\
\
\
/
I
\
\\
/
I
/
I
/
~\ ~I
I
'.. \ :.\
/
.,:;\~
\ \
I
,
\
\
\
,
I
,
o
I
>
63
64
SCIENTIFIC
COMPUTATION
Hodograph of Simple Turn
with Constant Velocity
on Suction Side
Simple Turn with Constant Velocity
on Suction Side
FIGURE
We now have four sets of numbers cp, 1/1, x, y known at
every point of the In w-() plane. We cover this plane by
a close square mesh of numbers, always four at every intersection. This number system we can reproduce as often
as needed. These systems then can be moved bodily with
respect to each .other. The numbers cp, 1/1, ~t", Y then can be
added at every location and a new combined solution
emerges. Obviously, it is best to move always in straight
multiples of the mesh interval.
Before addition each system can be multiplied by a
common multiplier. In view of what has been said before
such multiplication is necessary in the case .of the singularities representing the flow of infinity. Another multiplication is needed for the .;r-y system every time the
singularity is moved. A motion in the In w direction represents a relative shortening of physical dimensions, and
one in the () direction means a rotation .of the physical
system.
The equations for thiscorrectian are:
.;r = e-1n w (.t" cos 0 - y sin 0)
(13)
+ x sin 0)
(14)
y
= e-1n w
(y cos 7i
5
where .r,Yare the physical coordinates of the singularity
displaced by In 'W,. if in the hodograph plane. These multiplications and additiDns must be performed for each displaced singularity before it is added to the .others.
Procedure of Computation
One master stack of cards holds the four-value .table.
It can be reproduced onto as many stacks as singularities
are needed. This reproduction can already be guided in
such channels that the character of the singularity as
source, sink, positive .or negative vortex is taken care .of.
Nothing but a control of the signs is needed for this.
N ext the singularities must be moved to their predetermined place. This means a displacement of the origin of
each stack, in other words a change of identification. A
constant is added to either independent variable. This new
identification .orients the singularities with respect to each
other. Then we must multiply for strength and CDrrect the
.;r's and y's as shown above.
The field of interest for a particular computatian will
be smaller than that of the master table. The overflow
cards are now remaved. This involves .one sDrting.
FORUM
65
PROCEEDINGS
Here are some of the more simple examples which we
have done. We are collecting experience about the most
promising combination of singularities. Figure 5 shows· a
very simple flow turn which has a constant velocity on the
inside of the turn. By the addition of two singularities, one
source and one sink, we arrive at turning streamlines, one
of which has constant velocity throughout. The other
streamlines decelerate first and then accelerate. None of
these streamlines generates a closed profile.
The simplest closed curve we could generate is shown
in Figure 6. A vortex source flows into a vortex sink and
both are 0pP9sed by equal and opposite singularities. This
results in a closed curve in both hodograph (circle) and
physical plane. Here we have a saddle point which is not
situated at minus infinity. As a result, the profile does not
have a finite entrance angle at its nose. The velocity there
is not zero. There is, by the way, a very simple condition
which assures that if you have a closed curve in the hodograph map you will also get a closed physical flow curve.
In this manner we finish with as many stacks as there
are singularities. They are properly coded with respect to
each other. Next they are fed through the collator. Here
all cards belonging to the same coordinates are stacked together. This new total stack enters the accounting machine
where the values cp, tf' x, yare added together for each
coordinate. A new card is summary punched for each addition. The new resulting stack is the solution. Other singularities can be added to it if a modification of it is desired.
The solution is not yet in the fon11 in which we need it.
We must find the points for a number of (equally spaced)
constant values of tf (stream lines). The .t' and y values
appearing along these lines furnish coordinates of the
physical stream lines. One of these is the profile. What is
needed is a fast and simple inverse interpolation for tf and
a direct interpolation for .t' and y. In the absence of a fast
method, we use the old and time-honored method of cross
plotting. We hope to be able some day to do this part by
machine. As long as machines become faster there is hope.
60· - 60· Bucket
Turn
FIGURE
6
66
All that is necessary is that the strength of the infinity
singularities matches the components of the velocity at
which they are situated, and that the residue in the closed
hodograph figure be zero.
I may add that the fate of the compressible counterpart
of this computation is now entirely in the hands of IBM.
IBM is willing to perform the very complicated calculation
of the compressible singularities on the Selective Sequence
Electronic Calculator. After these tables are available we
can calculate subsonic compressible flow. Supersonic flow
is of small interest to us. \Ve leave that to the artillery !
Nothing will help us much, except being able to do two
things. One is to compute with a great deal of accuracy.
If we could sacrifice accuracy, we could proceed along
cheaper ways by experiment alone. The other is to be able
to mass. product theoretical results, because what we primarily are after are not solutions for production, i.e., for
the immediate turbine that goes into the shop. VYe want
series of solutions for experimental purposes. vVe want to
get experimental parameters which are related to blade
shapes, to the interaction between blade shapes, and to a
number of additional variables which now rather obscure
a clear conception of the working process of a turbine. If
this computation can help here we shall be able to produce
a still better performing machine than we have now.
SCIENTIFIC
COMPUTATION
DISCUSSION
Dr. Fenn: How justified is your irrotational flow with
the blades moving past each other?
Mr. Kraft: The flow can be considered completely irrotational as long as it is considered as a two-dimensional
problem. Even in three dimensions, when the turbine is
designed for constant circulation, you still would have an
irrotational problem. It becomes rotational only when you
take the boundary layer into account. This must come
necessarily after we can calculate irrotationally. We have
to follow the procedure which the great'masters of aerodynamics have laid down. \Ve do not think we know anything better.
Mr. Stevenson: Have you ever tried the classical aerodynamic scheme using conformal transformation?
Mr. Kraft: As I emphasized at the beginning, if it were
only the incompressible solution we were after, we would
not do this computation as I described it. We do it in this
manner only because we can replace it by the compressible
calculation as soon as the additional logarithmic singularity
tables for compressibility are available. We are well aware
that incompressible lattice computations can be performed
by simpler procedures.
Dynamics of Elliptical Galaxies
JACK
BELZER
Ohio State University
GEORGE
GAMOW
George Washington University
GEOFFREY
KELLER
Perkins Observatory
I N HIS C LAS SIC A L investigations on the dynamics
of elliptical galaxies, Sir James Jean has said, "We have
seen that the density of matter in the central lenticular
masses of nebulae is of the order of 10-21 g/cm 3 • The free
path in a gas of this density is 1014 cm., whereas the
diameter of the central mass of Andromeda nebula is
about 1.6 X 1021 cm . . . . It follows that the various
nebular configurations may legitimately be interpreted as
configurations of masses of rotating gas."1 He further
calculates the free path of a star, allowing for its gravitational interactions with other stars of the cloud, and finds
that for the same mean density it is equal to 10:.l9 cm.
(about 50,000,000 times the diameter of the nebula). He
concludes that the concept of gas pressure cannot be legitimately used in connection with nebular dynamics when
the nebula is supposed to be a cloud of stars, and that
clouds of stars should not assume the special shapes of
observed nebulae.
Today we know that the elliptical galaxies and the central bodies of the spirals are made up entirely of stars.
This fact, first suggested by the stellar type of their spectra, was established beyond any doubt by the recent work
of Baade 2 who was able to resolve the celestial objects into
a multitude of individual stars.
To -reconciliate Jeans' theoretical consideration with the
observed fact, George Camow suggests that the regular
shapes of elliptical galaxies were established in some past
epoch when they were entirely gaseous, and are now retained as "dead skeletons" after all the original gaseous
material was condensed into stars. This sttggestion agrees
with the recent theory of the evolution of the expanding
universe proposed by Dr. Camow,3 according to which the
masses and sizes 0 f galaxies can be predicted on the as-
sumption that they originated as the result of "gravitational instability" of the expanding primordial gas. These
gaseous clouds would tend to assume spherical or ellipsoidal forms depending upon the amount of angular momentum, and their internal density distribution was presumably that of the rotating isothermal gas-spheres. The
linear velocities of gas-masses at each point were proportional to the distance from the rotation axis. As the
star forming condensations took place under the forces of
gravitation and radiation pressure, the net outward force
due to the difference in gas (and radiation) pressure between the surfaces toward and away from the center of
the galaxy was reduced as a result of the reduction in the
volume and surface area. Consequently, the stars were
accelerated toward the center under the influence of unbalanced gravitational force, assuming that at no time during the history of the galaxy do the stars exchange any
appreciable amount of energy with their immediate neighbors as a result of close encounters.
The only overall forces acting on a star will be: 1. The
resistance of the gas to its motion, which will be a radial
force as long as the motion is radial and which will diminish in importance as the gas is consumed in the star formation process. 2. The radial force exerted by the smoothed
gravitational potential of the remainder of the stars. Then,
the stars will tend to oscillate radially through the center
of the galaxy. The amplitude of oscillation or the points
of maximum elongation of the newly acquired elongated
elliptical orbits of the stars should correspond to the distances at which they were originally formed. As more and
more stars were formed, the major axes of original elliptical orbits were gradually changing due to the gravitational action of other stars which have originated outside
67
68
SCIENTIFIC
of them, but are now penetrating during a certain fraction
of their motion into the interior regions. As a result of
these interactions, the stars acquired an altered distribution. It is expected, however, that the stellar orbits were
not shuffled; i.e" the stars which were formeda-t larger
distances from the center and therefore with larger angular moments' will also have their points of maximum elongation farther out from the center. In observing mean
radial velocity of stars at different distances from the
original rotation axis, with respect to an observer, we will
essentially obtain the tangential velocities of those stars
which pass through their point of maximum elongation at
that distance. Since the stars whiCh move beyond that
point will be much smaller in number and will also present
us only with the projection of their actual velocity, it can
be expected that the observed rotational velocities will increase with the distance from the center (linearly in the
first approximation), giving the impression that the entire
galaxy is rotating approximately as a solid body. This
"solid body" rotation is actually observed in elliptical
galaxies and in the central bodies of spirals, and was considered an unexplainable phenomena in view of the large
free passes of individual stars.
In the present work we. plan to analyze the observed
density distribution in the elliptical galaxies which was
very carefully measured by Dr. Hubble.'
Then the total mass between rand r
_
p ( r )dr Therefore
_~
r
= the present mass in the shell between radii
+
rand l'
dr
f3(a)da = the initial mass between the shells of radii
a and a + da
per) = the present potential at r.
A stellar mass starting at a and dropping inward t.o r will
lose the potential energy
1'1'llP(a) -.; P(r)]
= 11tV2/2
,8(a) da dr
T(a)v
R
dA
1
-=lM·da
a
dQ
1
= ..~:f·-dr
r
-
= modulus = .434
M
then
pea)
=
2y'2. M T(a)
7r
a
~[S~~(r)
. dA
A.
-
pea) dPd(g) dQJ. (1)
For integration, the Gregory formula is suggested, since
the errors of integration are readily estimated:
(1
.
1)
IJ.a+nw
W
af(x)dx = 2 fo + fl + f2 + ... + fn-l + 2,fn .
1
1
,
- 12
(I::::.l!n_l
- 1::::. 'f
0) -
1 (A
Aiif0)
24
Wiifn-2 - .W
19 (A
iiifn-3
- 720
W
Aiiif',0 ) _
-
W
•••
For obtaining derivatives near the end of the tabular sequence we use the following:
'lV/'(a) = I::::.if(a) -
~ I::::.iif(a) + j I::::.iiif(a)
-*
I::::. iv f (a)
+ ...
Away from the ends of the tabular sequence the following
central difference formula converges more rapidly:
wf'(a) = ~ (I::::.ia 112+ l::::. ia1/ Z)
1 ( WA
- 12
~O (I::::. va_ 1:l +
1
The time which the star will spend in the shell of thickness dr will be dr Iv (on the inward trip), and the fraction of all its time which it spends in this shell will be
dr/uT(a), where T(a) is the time required for the mass
to drop from a to the center. The contribution to the mass
observed now between rand r + dr arising from the mass
which started to drop from between a and a + do. is
r p(a)da
dr
T(a)v .
Jr
pea) da
pCr) = y'2Jr T(a) [pea) -- P(r)]I/Z'
and hence acquire the velocity
v = y'2[P(a) - per)] .
+ dr is
R
Computation of the initial density distribution in a spiral
galaxy before formation of stars:
Let p(r) dr
COM PUT A T I ON
I::::. va 1/2) - 2!0 (I::::. vila_liz
By definition
T(a) '"
f~
=
T
+
+
Ai")
W Ila1/2
I::::. vii a1/2 ) +
+
...
~2f [Pea) ! rp (r)]lI'
dr
Q
".
llla 1/~
1
r
r
= log1or ;dQ = dQ/dr = M = .434
a _ _1
__
( ) - My'2
r-a [Pea) rdQ
- Per) ]1/2
JA
(2)
FORUM
69
PROCEEDINGS
Near r = a, P (r) becomes nearly equal to P (a) and the
integrand diverges. To overcome this a series approximation in the neighborhood of r = a is suggested.
Let h = (a - r). Hence the first term is
~/-2-
-
a'JMaX yh
(3)
For small a the first term approaches zero, and in the
second we may set 3Ma = %-a3 p (0), and pea) = p(O).
Let h = a; then to the first approximation
of the problem is possible. It is the hope of the authors
that in the near future they will be in a position to publish
results of these investigations which will conclusively decide whether this theory is valid.
REFERENCES
1. SIR J. H. JEANS, Astronomy and Cosmogony (Cambridge,
1928), p. 335.
(4)
2. W. BAADE, "The Resolution of ... the Central Region of the
Andromeda Nebula," Astrophys. J., 100 (1944), p. 141.
3. G. GAMOW, "The Origin of Elements and the Separation of
Galaxies," Phys. Rev., 74 (1948), pp. 505-6.
4. E. HUBBLE, "Distribution of Luminosity in Elliptical Nebulae,"
Astrophys. J., 71 (1930), pp. 231-76.
The third term of the series should be calculated if greater
accuracy is desired.
The equations (1), (2), (3) and (4) with which we are
here concerned, lend themselves readily and very conveniently to solutions on standard IBM equipment. It is
gratifying that a computing laboratory of this nature is
available for these investigations. Although to date no conclusive results are available, a more thorough investigation
DISCUSSION
Mr. Hollander: Do we have any proof that stars are
moving radially, either in or out?
Mr. B elzel': It would take very many years to notice any
motion in the stars. I think that according to Dr. Gamow's
theory of the evolution of the universe, in a period of
about 200,000,000 years the stars might have made seven
or eight oscillations.
1~ I
1
TCa) = 4'J 6rrG p(O) .
Application of Punched Cards
in Physical Chemistry
GILBERT
W.
KING
Arthur D. Little, Incorporated
THE F 0 L LO WIN G concep~s of the application of
punched card machines are the conseq.ucmce of the particular problems in which the author is interested, but probably the examples are typical of present-day theoretical and
physical chemistry.
The Hollerith machine was invented for scientific work,
but reached the present stage of refinement because of its
application to commerce, and in the field of chemistry it
will return to science because of its bookkeeping facilities.
The chemical literature has reached such an enormous
volume that the tremendous burden of bookkeeping is a
serious detriment to efficient research. The application of
punched card techniques to indexing is the center of interest in chemistry today, and outweighs the other possibly
more important applications. One is tre recording, storage
and handling of experimental data. Next, there is as a rule
an opportunity for machine methods in the calibrating and
correcting of these data. After this there is the analysiS of
the data, various correlations to be searched for and other
statistics to be extracted. Finally, there are purely theoretical calculations.
parameters of the instrument the deflection can be converted to an absorption coefficient and the distance along
the· paper converted to wavenumbers. The absorption coefficients at a given wavenumber are invariant to the particularexperimental conditions; it is these quantities, not
the original record, which have universal interest. This
conversion and calibration can only be done after a first
step of visual "reading off" the charts, point by point. This
is a process of conversion from a continuous physical
variable to an abstract digital number. The conversion and
calibration now proceeds by digital calculations. Clearly, it
would be most desirable to record the original physical
variable (often a voltage )as a digital· number. At first
sight, it might seem that such a process is less accurate.
However, this is what the eye ultimately does in "reading
off" the chart. Furthermore, every good experiment should
be done with equipment that has a definite "noise" level
(in a general sense) just visible. Even measuring a distance should be done with a scale, or a micrometer eyepiece, or a cathetometer whose accuracy is just sufficient
to meet the requirements, or if a limiting factor, at its
extreme of sensitivity, so that each reading has a few
percent of "noise." The magnitude of the noise can be
taken as the unit in the digital scale; and a digital reading
is as accurate as the experiment, even though the continuous trace may look more precise.
If th~ readings are to be made digitally, the binary system is the· pest. First, a punched card works on a binary
system. There are two and only two "digits"-a hole or
no hole. Secondly, the binary system is most economical in
number of required digits per number (above eight). Finally, as a practical matter only one punching time is required to (multiply) punch a binary number <1024 in a
single column of a standard card, in contrast with the four
punching times for a four-digit decimal number.
The problem of recording data has been solved by a
Digital Reader, which reads a voltage, converts it to a
digital binary number and punches this in a card. An accuracy of one percent has been achieved with a very simple
circuit. The resulting seven digit « 128) binary number
Recording Experimental Data Digitally
Today, the recording of data is extremely widespread.
Almost always continuous variables are used to measure
quantities in an experiment or industrial process; the most
primitive is the visual reading of a meter scale, and recording of the nearest number on the scale. Recording
potentiometers are commonplace. Photographs of oscilloscope patterns are used in transient phenomena. Such records are rarely useful in themselves. Even if the primary
data on the record is sufficient by itself (for example, the
temperature recorded as a function of time), the study of
more than a few dozen charts by eye is almost impossible.
Usually the primary recorded data have to be converted to
some absolute quantities. An example is the. recording of
infrared spectra. The actual conventional record consists
of the distance a pointer is deflected for various angular
displacements of the paper roll. From a knowledge of the
70
FO RUM
PRO C E E DIN G S
is punched in one column of a card at the fastest punching
speed of IBM equipment, namely twenty-five a second.
Processing of Data
With the original measurements directly punched 011
cards as primary data with the original "noise" retained,
the conventional processing of the data can proceed with
machine methods. The example of infrared spectra brings
out the power of punched card methods in extracting more
information from experimental data than can be done by
continuous record. Hundreds of new spectra are being
taken every day. For the most part they are compared
with spectra taken in the same laboratory, and in almost
all cases, with spectra already extant in the literature. Several thousand infrared spectra are now generally available,
but their comparison is a difficult matter. In addition to
the various scales used in reducing drawings for publication, there are several conventions of plotting wavelength
or wavenumber, increasing or decreasing to the right. Finally, even with two records on the same scale, there are
some rather subtle comparisons to make. In general, the
problem is not whether two spectra match at every point,
but whether they have certain peaks in common, and to
what extent. This is similar to the problem of weather prediction carried out by punched card methods, based on the
matching, within tolerances, of today's weather map with
a file of the last forty years' maps. This is an ideal situation for a collator. Methods of calibrating the data, and
extraction of statistical information, are relatively conventional and are to be found in the literature.
There is a converse of the above problem, namely the
plotting of results of calculations done on cards. There is
no question that a plot is more readily understandable
than a table of values. It would be very desirable to have
an instrument which would take the impUlses from an
accounting machine and convert them to a continuous
quantity, presumably a voltage, which would drive a recorder. Lacking such equipment, we have derived a method
of using the type bars of a 405 to plot points. 1
Theoretical Calculations
At first sight, one might be' sure there are many opportunities for punched card methods in theoretical chemistry.
There are, however, several.arguments against their use.
The cost, if not the very presence, of a battery of IBM
machines in a chemical laboratory will have an undue influence on the type of work done there. There will be
emphasis on the ponderous data collecting and large scale,
poring over of material notable for its quantity rather than
quality, away from the elegant simple experiment. On a
more theoretical level, the interests and, ultimately, training of the research worker will be away from analysis and
71
be likely to get in a rut of conventional punched card approaches. There is no time to read Whittaker and Watson
when machines stand idle!
Another difficulty in the introduction of standard machine methods is the fact that IBM equipment is a parallel
type of computer. The scientific approach inherited from
pre-machine days is sequential. For example, many problems are developed with the aim of finding the solution as
the root of a polynomial. Now there are no analytical ways
of finding roots of polynomials. But IBM machines are not
suitable for the numerical solution of a single polynomial,
because they are efficient in parallel and not in sequential
calculations. If, however, the problem, perhaps by generalization, requires the solution of a few thousand polynomials, the process could be run in parallel even though
it required very many sequential steps. The' best known
numerical method of solving polynomials involves continued fractions, which are known to converge, a.nd hence
behave better than many approaches in which the root has
to be raised to a high power. The recent development of
the 602 and 604 which can divide has opened up the field
of extraction of roots to machine methods.
The repetitive feature of punched card machines is a
possible advantage to physical chemists, especially in the
construction and tabulation of functions. One example
would be the tabulation of the free energies of all substances for which information is available, from empirical
constants. This could be done to great advantage at every
degree from -273°C. to 5000°C. and thus save a great
deal of interpolation. Another application is the direct
calculation of the thermodynamic functions from fundamental frequencies from spectroscopic analysis, without
the forcing of such data into simple empirical formulas in
order to sum or integrate.
The purely theoretical problems of chemistry lie ina
field in which the role of machines, although certain, is
quite obscure. Chemists deal with molecules, and except
for hydrogen, a rigorous quantum mechanical approach is
beyond the power even of the most advanced machines
today. One difficulty with the quantum mechanics of a
molecule is that it is'many-dimensional, and integration
even in three dimensions in general involves astronomical
numbers of unit operations. It may be possible to approach
these problems with some sort of statistics, as Dr. Thomas
did so well for atoms. Another difficulty in theoretical
chemistry, as contrasted with fundamental physics perhaps, is that the solutions of problems do not involve
analytical functions. This appears even in the simplest
problems, such as the interpretation of infrared spectra.
The spectrum of water, for example, consists of several
thousand lines which are differences of a fewer number of
energy levels. The latter, however, are not spaced according to any elementary function.
72
They are roots of polynomials, and are non-elementary
functions of the moments of inertia, so can only be expressed as tables. Thus, advanced spectrum analysis can
only be approached by machine methods.
Still a further difficulty arises in the real problems of
physical chemistry from the nature of experimental data
themselves. The unravelling of the rotational structure of
an infrared spectrum has been described above as being
very complex. One might infer that, with machine methods
and enough time, a unique interpretation could be made.
This is true when the spectrum is completely resolved, as
in the photographic infrared. In the main infrared region,
lines are not completely separated. The spectrum consists
of a continuous record consisting of a number of peaks,
each of which contains a number of more or less resolved
"lines" which are not mathematical lines, but Gaussian
curves. Thus, the conventional analysis by finding a set of
constant differences between the various lines fails, because the overlapping of lines in peaks leads to uncertainties in their position of the order of magnitude of
significant differences in the differences.
There is an analogous situation in many x-ray and electron-diffraction spectra. The classical approach of deducing a unique structure from its observational data cannot
be carried through. Modern spectrum analysis almost
always is achieved by the stochastic process. A structure
of the molecule is assumed. Certain theoretical expressions
are used to calculate the appearance of the spectrum in
which such a structure would result. The calculated spectrum is compared with the observed. Successive trial structures are assumed until a satisfactory agreement between
calculated and observed spectra is achieved. There is no
proof the final fit is unique, but as a rule the observed
spectra have sufficient complexity that the probability of
any other structure giving rise to the same or better fit is
remote. The application of the stochastic method to the
analysis of rotational structure of infrared spectra has
been described in detail elsewhere. 2
It is appropriate to review briefly the wiring of the various IBM machines in this work.* Many of the calculations
are straightforword. The only unusual feature in the
collator is the" effective use of splitting up the comparing
magnets into nine fields. This allowed comparing, and
hence merging on three fields of quantum numbers simultaneously, two of which increased and one decrea~ed.
The whole process consisted of forty steps, of which
seventeen were done on the accounting machine. It might
be interesting to point out that all calculations were done
on the same control panel. The most noteworthy single
calculation done on this control panel was the calculation
*Complete wiring diagrams will be presented in a report to the
Office of Naval Research who supported some of the spectrum
analysis with punched card equipment.
SCIENTIFIC
COMPUTATION
of the expected spectrum with a finite slit. A deck of cards,
one for each small absorption line in the spectral region
( 300 cm-l wide) is made as the final step in the purely
theoretical treatment. Each line is represented by a card
giving its position in wavenumbers and its intensity. These
cards are summary punched over a small interval (0.5 cm-l )
of the spectrum, so that lines less than 0.5 wavenumbers
apart are combined in intervals. This deck of summary
cards is made complete with blanks, to give a final deck, a
card appearing for every 0.5 cm-l interval in the whole
region. This deck then represents the appearance of the
spectrum if the resolving power were 0.5 cm-l. In actuality
it is higher. A region of wavenumbers is seen by the slit
of the spectrometer at anyone time. It is desirable then
to compute the actual transmission at each wavenumber v
2:
+2
T ( v) =
P(}" •
1 (v
+ (}") ,
(}"=-2
where P-2' P-l1 Po, Pv P2' are weight factors describing how
much on either side of v the slit sees. This summation is to
be carried out for each l', proceeding by intervals of 0.5 cm-l
along the region of 300cm-l . It was found that values of
P = 1,2,4,2,1 or 1,2,1 satisfactorily expressed the experimental slit shape function. Summary cards T(v) were obtained .in one pass through the accounting machine, using
progressive totals in five counters, where counter entries
were fed the number 1 (v) by a permuting switching arrangement. A minor total cycle occurred every card,. One
counter has tn be cleared every card. This was done connecting the counter exit to entry through the permuting
class selectors.
It should be mentioned that some calculations were
made possible on the single accounting machine control
panel by the use of external manually controlled gang
switches which connected different fields of the cards via
brushes to different counters, which were too few in number. The well-known "Octal" tube socket was used for
additional hubs, the standard IBM wire fitting the holes
in these sockets somewhat better than the hubs in the control panel!
Another large field in physical chemistry where large
scale computing machinery will playa great role is statistical mechanics. The properties of a single molecule have
to be determined from experimental measurements by a
stochastic process as outlined above. But the macroscopic
behavior of materials is an average of certain quantities
of an extremely large number of molecules, each one moving, vibrating, rotating with different velocities, or in a
different configuration. A very simple example is the
theory of rubber-like elasticity as based on a simple theory
of the statistical mechanics of high polymers. Rubber
molecules are long chains of ten thousand or more seg-
FORUM
73
PROCEEDINGS
ments. A molecule can therefore exist in very many different configurations, and as a consequence has high
entropy. On stretching, the molecules become less random,
and the chains tend to become parallel. The decrease in
entropy is responsible for the strong force in a piece of
stretched rubber which attempts to contract it to its original state of maximum randomness. Clearly the force, and
hence the modulus of elasticity of a plastic made of long
chain molecules, could be calculated if we could enumerate
the possible configurations of such an assembly of chains.
This can be done in a simple way and leads to results in
rough agreement with experiments, at least for highly
elastic substances. The simple theory has a number of defects, principally the fact the molecules have volume and
can only occupy the same volume of space once. The problem of this effect of "excluded volume" on the number of
configurations is a topological one of great difficulty. We
have set out to solve this topological problem by a straightforward enumeration of the configurations of chains. In
essence, we have investigated the famous "random walk"
problem, in a tetrahedral lattice accounting for the effect
of excluded volume.
This, then, is an example of the use of punched cards in
sampling a Gibbsian ensemble, in which each system is
described appropriately on a set of cards. The required
statistical averages can be very readily made by arithmetical means by conventional processing of these samples
on cards.
RBI<'BRENCES
1. G. W. KING, "A Method of Plotting on Standard IBM Equipment," MT AC, III (1949), pp. 352-55.
2. G. W. KING, P. C. CROSS, and G. B. THOMAS, "The Asymmetric
Rotor. III. Punched Card Methods of Constructing Band Spectra," J. Chem. Physics, 14 (1946), pp. 35-42.
DISCUSSION
Mr. Bell: What mechanical or electrical system is between your prism and the punch?
Dr. King: A whole bunch of relays; it took me a day to
wire them up. The method is not very complicated.
Dr. Caldwell: This problem of converting a voltage to
a number on punched cards is going to occur more frequently and in more difficult form than is described here.
One case involves the fact that the voltage is coming from
a device over which you have no control. I think the general approach is that you can apply the pulse code modulation system. Within the last year and a half, the Bell
Technical Journal has contained several papers on that
subject, the general principle being to compare voltage
patterns against a standard voltage. That means that you
can easily produce these binary digits; the real problem is
that they are coming too fast for any punch to record
directly.
Dr. King: vVe don't have a storage problem because we
get all the digits out more or less simultaneously. We take
a new reading as soon as the old one is punched.
Dr. Grosch: lVIay I make a remark about this question
of finding the roots of a polynomial? There is a class of
problems which can be handled in parallel fashion-polynomials of degree, say, five to twelve. If you have several
to do at once, so much the better. The process is to evaluate the nth degree polynomial for n
1 equal spaced values of the variable x, without rounding. This is a severe
limitation, since it involves large multiplications toward
the end; the answers will usually be about n + 6 significant
figures for n<20. These exact answers are then differenced n times; the single nth difference should equal n !wn
times the original coefficient of x n , where w is the interval;
this is used as a check. Further values of the polynomials
can then be built up from the constant nth difference on
the 405 three or four orders at a time, with summary
punched intermediate results. The tabulation of the polynomial can be extended to cover all real roots, or restricted
to the neighborhood of a single root. Inverse interpolation
for the exact root is the final step; usually this is hand
work.
I f the above procedure seems too complicated, another
trick worth trying on the 602 or 604 is to prepare a deck
of n + 1 cards carrying the coefficients of the polynomial,
highest degree coefficient first. A card carrying a guess at
the root is placed in front of the deck, and one pass through
the calculating punch evaluates the polynomial. A new
guess is made, key punched on a fresh card, and the
process repeated. The 602 or 604 is here used as a specialized desk calculator.
+
Application of Punched Card Methods to the Computation
of Thermodynamic Properties of Gases from Spectra*
LYDIA
G.
JACK-
SAVEDOFF
HERRICK
BELZER
L. JOHNSTON
Ohio State University
THE P R INC I P L E S of the statistical calculation of
thermal functions of gases from spectroscopic data have
been used for many years with a large degree of success.
The basis for much of this work is the so-called "summation method" which was first used by Hicks and Mitchell. 1
From the distribution of the molecules of a gas among
the various energy levels of the molecule, the energy and
other thermal properties of the gas can be found through
the evaluation of the partition functions. In the summation
method, this summing takes place over the actual energy
levels of the gas. When there are large numbers of energy
levels, as in all but the simplest molecules, the direct application of this method becomes impractical and labor saving devices must be used. In the past, these have taken the
form of substitution of integrals. for the series involved as,
for example, in the extended Mulholland treatment. 2 The
present paper introduces a method whereby the original
summation method is adapted to use of punched cards.
The energy levels of the molecule may be obtained
through a study of the band spectrum of the gas. In the
equations for the thermodynamic properties for the gas
from these energy levels, the contributions due to translation may be separated from the effects due to rotation
and vibration, thus simplifying the expressions to be evaluated. Johnston and Chapman3 have introduced the following notation whereby the contributions of the rotational
and vibrational energy of the molecule may be expressed
simply in terms of three basic quantities:
!A = !Pi e-f.i/kT
For convenience, these three quantities may be redefined
in terms of base ten exponentials:
!A* = !Pi 10-f.i/zkT
!S*:;:: !Pi f.i/ zkT . 10-f.i/zkT
= !Pi f.i e-f.t/ kT
(2b)
i
!c* = !Pi f.l / zkT . 10-f.-ijzkT
(2c)
i
In the above equations f.i is the energy of the molecule in
the ith quantum state, Pi the statistical weight of this state,
k the Boltzmann constant, T the absolute temperature of
the gas, and z the numerical constant In 10~
The thermodynamic functions are given by the following equations.:
EO -
E~
= 3/2 RT + zRT !B*/!A*
C$ = 5/2 R
(3a)
+ z2R [!C*/!A*
- (!B*/!A*)2]
(3b)
SO = 3/2zR 10gA-f + 5/2 zR log T - 2.3140
+ zR [log!A* - !B*/!A*] (3c)
FO - Eg = 5/2 R - 3/2 zR log M - 5/2 zR log T
+ 2.3140 - zR log !A*. (3d)
The final term in each of the equations (3) gives the con;.
tribution of the rotational and vibrational states, while the
preceding terms are due to the translation of the molecules.
New quantities appearing in these equations are the gas
constant R, and the molecular weight of the gas M.
This laboratory, with the assistance of Dr. Thomas of
the Watson Scientific Computing Laboratory, has set up
a punched card method for obtaining the quantities !A*,
!B*, and !C*, once a punched card table of the energy
(la)
i
!B
(2a)
i
(lb)
i
(lc)
*This work was supported in part by the Office of Naval Research
under a contract with the Ohio State University Research Foundation.
74
FOR UM
75
PRO C E E DIN G S
levels is available. These are readily computed in many
cases. Since the sums are to be computed for a number of
different temperatures, it becomes c.onvenient to transform
the exponentials of equations (2) into double exponentials. The sums A*, B*, and C* now have the form
~A* = ~Pi,
10-10$£
(4a)
i
~B*
= ~Pi 10»i 10-10'l'i
(4b)
i
~C*
=
~Pi 10::: x i
10-10:"i
i
(4c)
'i -
where Xi = log f.i/zk - log T =
l.og T. The advantage
of this will be apparent when the calculations of these
sums at different temperature is discussed. For the present,
the temperature will be set at 10 K., so that Xi = Ct. A
master deck of the functions 10-10 ; lOX 10-1°:'10:::x 10-10Ji was
prepared for the argument x = - 8.00 (.02) + 1.00. This
set includes all tht; significant values of the three functions
(which vary between zero and one). Before thi~ table can
be utilized, it is first necessary to obtain the logarithms
of the quantities f.il zk and reduce these to the arguments
of the master table. To achieve the former of these objects,
a six place optimum interval logarithm table is used,
whereby Xi is obtained to six decimal places and one whole
number. Four point Lagrangian interpolation coefficients 4
are used to distribute the statistical weights P'i of the
energy levels f.i to the appropriate arguments of the master
deck.
The. punched card Lagrangian interpolation table built
up for this purpose consists of 5000 cards at intervals of
0.0001 in the argument, P = 6x/h, for the range 0.00010.5000. Since the table reflects itself for the range 0.50001.0000, the table is punched in a symmetrical manner
around the centeri.of the card, with the right half reading
in reverse from right to left.. Tumbling the cards puts
them in position for use in the upper half of the interval
range. A secondary argument, q = 2p, is also punched on
the table so that the table may be collated directly with
the last four digits of Xi. To accomplish this most readily,
the detail cards. are first divided int.o odd and ev.en groups
on the second decimal digit of .~'i to determine which half
/of the table is to be used.
The use of the .interpolation coefficients is shown most
readily by the following illustration. Suppose that for an
energy level f.i of statistical weight Pi, the log f.i/ zk is equal
to 0.498306. The distribution .of the statistical weight
among the four adjacent arguments is given in Table 1.
TABLE
I
Xn
= 0.46
= 0.48
P~-l = PiA-l
p~ = PiAo
X n +1
= 0.50
X n+2
= 0.52
P~+1
P~+2
X n -1
= P·iA
= piA:
l
The four coefficients, A_i' Ao, Ai' and A2 are those corresponding to the argument q = 1.8306. A prepared table of
~.
~.....
and X 1H2 with the argument Xn is also repro""n-H
·~1h """+1
duced .onto these same detail cards so that the proper
argument may be matched with the new weights P~t. In
general, there will be a considerable amount of overlapping, and so the P~t's are next summed for each value
of Xn such that a total weight, W,t = ~p~, is obtained for
each different value of x". In practice, it has been found
most convenient to make up three additional sets of cards
so that the arguments ,t'n and the products P:t = p.tA always
appear in the same fields and the new weights are obtained
easily on the accounting machine.
The advantages of the procedure described above are
twofold. In the first place, there is a great reduction in the
total number of terms over which the sums are to be
taken. Even in a relatively simple gas, the total number of
thermodynamically important energy levels lies in the
thousands, whereas the master deck of functions for these
new arguments contains at most 400 cards. Secondly, the
functions 10-10~ lOX 10-10: 102 .v 10-1O'v forthe arguments X at
intervals of 0.02 are available in a permanent file and need
only be reproduced into the detail cards to .be multiplied
and summed.
Up to this point, the temperature has been neglected.
However, by shifting the table cards one value (0.02) of
x, the detail cards may be prepared for a new temperature.
The entire temperature range may thus be covered in this
way in logarithmic increments of 0.02 in T. This is the
major advantage of the double exponential form adopted
for the problem. The range of log T from 0.00 to 3.84 is
sufficient to c.over all temperatures up to 6OOO 0 K. The
manipulation of the sets of cards to take account .of the
temperature shifts is more clearly shown by the following
numerical example, taken from an actual calculation. The
symbol 'n has been introduced for the Xn used earlier to
av.oid confusi.on. The smallest value of
is 0.24, and the
first set of cards prepared is that for the highest value of
log T which is equal to 3.84 if the computations are to be
carried up to 6000 0 K. Successive cards in this set will
have arguments and weights as given in Table II.
'n
TABLE
II
W1
'1 = 0.24
w2
'2
= 0.26
.1:'2
'3
= 0.24 = 0.26 -
= 0.28
Xa
= 0.28 - 3.84 = -3.56 = 4.44, etc.
Wa
Xl
=
3.84 =
3.84
= 4.40
-3.58 = 4.42
-3.60
The last card in the set should have the argument X = 1.00.
At this point, the functions have dropped to zero and any
left .over cards may be discarded.
The second set of cards will be for the next lower temperature, log T = 3.82. The master set is shifted up one
76
SCIENTIFIC
hydrogen by the authors. The major error in the above
procedure is due to the use of the four point interpolation
coefficients. The maximum error in each of the sums is of
the order of 2X 10-1 "$.Pi. In cases where the "$.Pi and hence
this absolute error might seem to become unduly large, the
values of the sums themselves are large ·and the relative
error remains within reasonable limits. It might also be
pointed out that this error is approximately half the maximum error attained by the hand computing procedure used
previously for the direct summation method. The agreement of the results as given in Table IV is quite satisfactory.
The authors wish to express their thanks to Dr. W. ].
Eckert, Dr. L. H. Thomas, and other members of the staff
of the Watson Scientific Computing Laboratory for their
invaluable assistance in getting this program started.
value, thus wand the X 2 are matched. The last value of Wn
will not match any table card and may be discarded. This
is shown in Table III.
TABLE
W1
'1
'W2
'2
Wg
'3
= 0.24
= 0.26
= 0.28
"t'2
Xg
"t'4
= 0.24 = 0.26 = 0.28 -
COMPUTATION
III
3.82
= - 3.58 = 4.42
3.82 = - 3.56 = 4.44
3.82 = -3.54 = 4.46, etc.
The process is repeated until the entire temperature range
has been covered. It has been found possible to condense
the number of individual temperature sets to one third by
including all the necessary data for three temperatures on
one set. The sums A*, B*, and C* may now be obtained
according to equations (4). Since the individual products
were not required, this was done by progressive digiting.
There may be as many as thirteen digits in W n , so that this
is quite a prodigious task. It is particularly unpleasant to
handle since the individual groups are so small.
Once the sums have been obtained for the arbitrary
temperatures, they in turn are interpolated using the four
point coefficients, yielding the results for the desired temperatures. From these, the thermodynamic properties
may be readily computed by using the relations given in
equations (3).
A comparison of the hand and machine computed
results is made in Table IV. The data is quoted from a
forthcoming paper on the thermodynamic properties of
REFERENCES
1. H.C. HICKS and A. C. G. MITCHELL. "The Specific Heat and
Entropy of Hydrogen Chloride Derived from Infra-Red Band
Spectra," 1. Am. Chem. Soc., 48 (1926), pp. 1520-27.
2. H. L. JOHNSTON and C. O. DAVIS, "Heat Capacity Curves of
the Simpler Gases. IV. Extension of the 'Free Energy' Formula
of Giauque and Overstreet ...•" J. Am. Chem. Soc., 56 (1934),
pp.271-76.
3. H. L. JOHNSTON and A. T. CHAPMAN, "Heat Capacity Curves
of the Simpler Gases. 1. Heat Capacity, Entropy, and Free
Energy of Gaseous Nitric Oxide from Near Zero Absolute to
5000 o K.," 1. Am. Chem. Soc., 55 (1933), pp. 153-72.
4. MATHEMATICAL TABLES PROJECT, op. cit., pp. 104-203.
DISCUSSION
[Discussion of this paper was omitted because of time limitations.]
TABLE
IV
COMPARISON OF THERMODYNAMIC FUNCTIONS
COMPUTED BY HAND AND MACHINE METHODS
Temperature
1500 0 K
Function
'2A*
};B*
'2C*
EU-E~
Cpo
SJ
-(F o_ Eg)IT
5000 0 K
'2A*
};B*
};C*
EO-Eg
COl'
SO
-(Fo-E~)/T
Hand
Computed
Machine
Computed
PerCent
Deviation
18.76256
8.86564
9.07780
7711.5
7.7103
42.702
35.574
18.76255
8.86569
9.07783
7711.51
7.71025
42.7018
35.5743
.00005
.0006
.0008
95.79684
69.60681
90.85673
31516.1
9.39459
53.0818
44.7921
95.80611
69.62325
90.87847
31518.3
9.39447
53.0823
44.7923
.0096
.024
.024
.0070
.0013
.0013
.0004
Calculation of the Equilibrium Composition
of Systems of Many Constituents
ROBERT W. SMITH, JR.
STUART R. BRINKLEY, JR.
U. S. Bureau of Mines
THE COMPOSITION of a system at chemical
equilibrium is easily calculated when there is only a single
reaction to be considered. In this case, the concentrations
of each constituent can be related to a single variable, "the
degree of reaction," and the solution of the mass-action
equation is straightforward. Difficulties are encountered
if this method is extended to a consideration of two simultaneous equilibria, and when the number of such simultaneous equilibria becomes large, the ordinary methods
become very laborious.
There is need for a systematic procedure designed to
provide a method for writing down the necessary relations
in the form most appropriate for numerical computation.
When the number of constituents is large, the relations
must usually be solved by an iterative procedure. In the
course of an extended program of such calculations, it is
usually necessary to formulate a number of computational
procedures, in order to assure sufficiently rapid convergence. If the calculations are to be carried out by punched
card methods, it is desirable that the smallest possible number of arithmetical operations of different kinds be involved in order to minimize the number of different control panels required. In a recent publication,! a systematic
procedure for calculating the equilibrium composition of
a system of many constituents was presented. This method
presents a simple rule for formulating the work program
of such calculations, with the result that very little time is
required for setting up a particular problem. The systematic nature of the computational procedure makes
the method well-adapted to punched card methods. The
method has been routinely employed in this laboratory in
a long series of such calculations. The method is as easily
applied to a mixture with a very large number of constituents as to a mixture with a small number of constituents, although the time required to obtain the solution
would be greater for the former case than for the latter.
In the publication cited, the method was developed for
systems of a very general nature. In the present communication, we restrict application of the method to the calcu-
lation of the equilibrium composition of mixtures consisting of a single homogeneous gas phase, and we assume
that the gas phase is adequately described by the ideal gas
equation of state. By taking advantage of these restrictions, we are able to formulate a computational method,
applicable to these particular cases, which is substantially
simpler and more systematic than the more general
method. A large number of systems of industrial and
academic importance is included in this category.
THE COMPONENTS
In a system containing many constituents, it is possible
to select certain constituents which are sufficient to describe the composition completely. By this it is meant that
if the system is conceived to consist of the selected constituents only, its gross composition (in terms of the
amounts of each chemical element present) is completely
defined. The constituents thus sufficient to describe the
composition are called the components of the system. An
analytical criterion has been published 2 for the choice of
the components. In terms of this criterion, the conditions
, for equilibrium may be written in a form which has a high
degree of symmetry and is particularly well adapted for
formulating a computational method for the calculation of
the equilibrium composition.
The number of constituents of any system depends upon
the accuracy with which it is desired to describe its composition. The constituents to be considered must be chosen
a priori, and this choice usually will imply the neglect of
certain equilibria that may be expected to exert a negligible effect on the composition of the system at equilibrium.
Consider a closed system containing s different substances, which are assumed to be in chemical equilibrium.
The molecular formula of the £th substance may be represented by
(1)
i
= 1,2, ... s, where X(k) is the symbol of the kth element,
ail,;
77
is the subscript (which may be zero) to this symbol in
78
SCIENTIFIC
the formula of the ith substance, and m is the total number
of elements represented in the system. For every i, the
array of sUbscripts ail., k = 1, 2, ... , ln, may be said to
define a vector
(2)
which may be called the formula vector of substance i. If
the rank of the matrix of the vector elements ail. is c, it
follows from a well-known theorem of algebra that there
are c linearly independent vectors, and if c0 ~.r:- }>0
,,0
CI
0
I
0
13m
,
I
o 0.0
0
1
2
3
4
3 00.0
3 00.3
3 00.5
3 00.8
301.0
300.5
3 00.7
3 01.0
301.2
301.5
251.8
252.0
252.3
252.5
252.8
1
2
3
4
5
6
7
8
9
301.3
3 01.5
301.8
3 02.0
3 02.3
301.7
302.0
3 02.2
3 02.5
302.7
253.0
253.2
253.5
253.7
253.9
10
11
12
13
14
3 02.5
3 02.8
303.0
303.3
3 03.5
303.0
3 03.3
3 03.5
3 03.8
3 04.0
15 ' 3 03.8
16 3 04.0
17 3 04.3
18 3 04.5
19 3 04.8
20
21
22
23
24
3
3
3
3
3
•.r:-
f....
,,0
;]>0
CI
•.r:- Sec. SUN
PLANET
f....
,,0
0
I
I
0.0
0.0
0.1
0.1
60
61
62
63
64
1.3
1.3
1.3
1.3
1.3
120
121
122
123
124
2.5
2.5
2.5
2.6
2.6
5
6
7
8
9
0.1
0.1
0.1
0.2
0.2
65
66
67
68
69
1.4
1.4
1.4
1.4
1.4
125
126
127
128
129
2.6
2.6
2.6
2.7
2.7
254.2
254.4
254.7
254.9
255.1
10
11
12
13
14
0.2
0.2
0.3
0.3
0.3
70
71
72
73
74
1.5
1.5
1.5
1.5
1.5
130
131
132
133
134
2.7
2.7
2.8
2.8
2.8
304.3
3 04.5
3 04.8·
305.0
3 05.3
255.4
255.6
255.9
2 56.1
256.3
15
16
17
18
19
0.3
0.3
0.4
0.4
0.4
75
76
77
78
79
1.6
1.6
1.6
1.6
1.6
135
136
137
138
139
05.0
05.3
05.5
05.8
06.0
3 05.5
305.8
3 06.0
3 06.3
306.5.
256.6 20 0.4 80
256.8 21 0.4. 81
257.0 22 0.5 82
257.3 23 0.5 83
257.5 24 0.5 84
1.7
1.7
1.7
1.7
L8
25
26
27
28
29
306.3
3 06.5
3 06.8
3 07.0
3 07.3
306.8
3 07.0
307.3
3 07.5
3 07.8
2 57.8 25 0.5
258.0 26 0.5
2 58.2 27 0.6
258.5 28 0.6
258.7 29 0.6
85
86
87
88
89
30
31
32
33
34
307.5
3 07.8
3 08.0
3 08.3
3 08.5
3 08.0
3 08.3
308.5
3 08.8
3 09.0
259.0
2 59.2
2 59.4
259.7
2 59.9
30
31
32
33
34
0.6
0.6
0.7
0.7
0.7
35
36
37
38
39
3
3
3
3
3
08.8
09.0
09.3
09.5
09.8
3
3
3
3
3
09.3
09.5
09.8
10.0
10.3
3 00.2
3 00.4
300.6
3 00.9
3 01.1
35
36
37
38
39
40
41
42
43
44
3
3
3
3
3
10.0
10.3
10.5
10.8
11.0
3
3
3
3
3
10.5
10.8
11.0
11.3
11.5
3
3
3
3
3
45
46
47
48
49
3
3
3
3
3
11.3
11.5
11.8
12.0
12.3
31L8
3 12.0
3 12.J
3 12.5
3 12.8
50
51
52
53
54
3
3
3
3
3
12.5
12.8
13.0
13.3
13~5
3
3
3
3
3
55
56
57
58
59
3
3
3
3
3
13.8
14.0
14.3
14.5
14.8
3
3
3
3
3
60
3 15.0
,
0 3 15.0
1 3 15.3
2 3 15.5
3 3 15.8
4 3 16.0
'Y'
0
MOON
,
0
,
,jf ~.r:- ;]>0
,,0
Ci
CI
.• .r:f....
,,0
,
I
;]>0
•.r:-
CI
,,0f....
,
3
3
3
3
3
15.5
15.8
16.0
16.3
16.5
306.1
306.4
3 06.6
3 06.8
3 07.1
o 0.0
1
2
3
4
0.0
0.0
0.1
0.1
60
61
62
63
64
1.4
1.4
1.4
1.4
1.4
120
121
122
123
124
2.7
2.7
2.7
2.8
2.8
5
6
7
8
9
0.1
0.1
0.2
0.2
0.2
65
66
67
68
69
1.5
1.5
1.5
1.5
1.6
125
126
127
128
129
2.8
2.8
2.9
2.9
2.9
3
3
3
3
3
16.3
16.5
16.8
17.0
17;3
3
3
3
3
3
16.8
17.0
17.3
17.5
17.8
3 07.3
307.5
3 07.8
308.0
3 08.3
12
13
14
3
3
3
3
3
17.5
17.8
18.0
18.3
18.5
3
3
3
3
3
18.0
18.3
18.5
18.8
19.0
3
3
3
3
3
08.5 10 0.2
08.7 11 0.2
09.0 12 0.3
09.2 13 0.3
09.5 14 0.3
70
71
72
73
74
1.6
1.6
1.6
1.6
1.7
130
131
132
133
134
2.9
2.9
3.0
3.0
3.0
2.8
2.8
2.9
2.9
2.9
15
16
17
18
19
3
3
3
3
3
18.8
19.0
19.3
19.5
19.8
3 19.3
319.5
3 19.8
3 20.0
3 20.3
3
3
3
3
3
09.7
09.9
10.2
10.4
10.7
15
16
17
18
19
0.3
0.4
0.4
0.4
0.4
75
76
77
78
79
1.7
1.7
1.7
1.8
1.8
135
136
137
138
139
3.0
3.1
3.1
3.1
3.1
140
141
142
143
144
2.9
2.9
3.0
3.0
3.0
20
21
22
23
24
3 20.0
3 20.3
3 20.5
3 20.8
32LO
3 20.5
3 20.8
321.0
3 21.3
3 21.6
3
3
3
3
3
10.9
11.1
11.4
11.6
11.8
20
21
22
23
24
0.5
0.5
0.5
0.5
0.5
80
81
82
83
84
1.8
1.8
L8
1.9
1.9
140
141
142
143
144
3.2
3.2
3.2
3.2
3.2
1.8
1.8
1.8
1.8
1.9
145
146
147
148
149
3.0
3.0
3.1
3.1
3.1
25
26
27
28
29
3
3
3
3
3
21.3
21.5
21.8
22.0
22.3
3
3
3
3
3
21.8
22.1
22.3
22.6
22.8
3
3
3
3
3
12.1
12.3
12.6
12.8'
13.0
25
26
27
28
29
0.6
0.6
0.6
0.6
0.7
85
86
87
88
89
L9
1.9
2.0
2.0
2.0
145
146
147
148
149
3.3
3.3
3.3
3.3
3.4
90
91
92
93
94
1.9
1.9
1.9
1.9
2.0
150
151
152
153
154
3.1
3.1
3.2
3.2
3.2
30
31
32
33
34
3
3
3
3
3
22.5
22.8
23.0
23.3
23.5
3
3
3
3
3
23.1
23.3
23.6
23.8
24.1
3
3
3
3
3
13.3
13.5
13.8
14.0
14.2
30
31
32
33
34
0.7
0.7
0.7
0.7
0.8
90
91
92
93
94
2.0
2.0
2.1
2.1
2.1
150
151
152
153
154
3.4
3.4
3.4
3.4
3.5
0.7
0.8
0.8
0.8
0.8
95
96
97
98
99
2.0
2.0
2.0
2.0
2.1
155
156
157
158
159
3.2
3.3
3.3
3.3
3.3
35
36
37
38
39
3
3
3
3
3
23.8
24.0
24.3
24.5
24.8
3
3
3
3
3
24.3
24.6
24.8
25.1
25.3
3
3
3
3
3
14.5
14.7
14.9
15.2
15.4
35
36
37
38
39
0.8
0.8
0.8
0.9
0.9
95
96
97
98
99
2.1
2.2
2.2
2.2
2.2
155
156
157
158
159
3.5
3.5
3.5
3.6
3.6
01.3 40
01.6 41
01.8 42
02.1 43
02.3 44
0.8
0.9
0.9
0.9
0.9
100
101
102
103
104
2.1
2.1
2.1
2.1
2.2
160
161
162
163
164
3.3
3.4
3.4
3..4
3.4
40
41
42
43
44
3
3
3
3
3
25.0
25.3
25.5
25.8
26.0
3
3
3
3
3
25.6
25.8
26.1
26.3
26.6
3
3
3
3
3
15;7
15.9
16.1
16.4
16.6
40
41
42
43
44
0.9
0.9
0.9
1.0
1.0
100
101
102
103
104
2.3
2.3
2.3
2.3
2.3
160
161
162
163
164
3.6
3.6
3.6
3.7
3.7
3
3
3
3
3
02.5
02.8
03.0
03.3
03.5
45
46
47
48
49
0.9
1.0
1.0
LO
LO
105
106
107
108
109
2.2
2.2
2.2
2.3
2.3
165
166
167
168
169
3.4
3.5
3.5
3.5
3.5
45 ·326.3
46 3 26.5
47 3 26.8
48 3 27.0
49 3 27.3
3
3
3
3
3
26.8
27.1
27.3
27.6
27.8
3
3
3
3
3
16.9 45
17.1 46
17.3 47
17.6 48
17.8 49
1.0
1.0
1.1
1.1
1.1
105
106
107
108
109
2.4
2.4
2.4
2.4
2.5
165
166
167
168
169
3.7
3.7
3.8
3.8
3.8
13.0
13.3
13.5
13.8
14.0
3
3
3
3
3
03.7
04.0
04.2
04.4
04.7
50
51
52
53
54
1.0
1.1
1.1
1.1
1.1
110
III
112
113
114
2.3
2.3
2.3
2.4
2.4
170
171
172
173
174
3.5
3.6
3.6
3.6
3.6
50
51
52
53
54
3
3
3
3
3
27.5
27.8
28.0
28.3
28.5
3
3
3
3
3
28.1
28.3
28.6
28.8
29.1
3
3
3
3
3
18.0
18.3
18.5
18.8
19.0
50
51
52
53
54
1.1
1.1
1.2
1.2
1.2
110
III
112
113
114
2.5
2.5
2.5
2.5
2.6
170
171
172
173
174
3.8
3.8
3.9
3.9
3.9
14.3
14.5
14.8
15.0
15.3
304.9
305.2
3 05.4
3 05.6
305.9
55
56
57
58
59
1.1
1.2
1.2
1.2
1.2
115
116
117
118
119
2.4
2.4
2.4
2.5
2.5
175
176
177
178
179
3.6
3.7
3.7
3.7
3.7
55
56
57
58
59
3
3
3
3
3
28.8
29.0
29.3
29.5
29.8
3 29.3
3 29.6
3 29.8
330.1
330.3
3
3
3
3
3
19.2
19.5
19.7
20.0
20.2
55
56
57
58
59
1.2
1.3
1.3
1.3
1.3
115
116
117
U8
119
2.6
2.6
2.6
2.7
2.7
175
176
177
178
179
3.9
4.0
4.0
4.0
4.0
2.5 180
3.8
60
330.0
3 30.6
3 20.4 60 1.4 120
2.7 180
4.1
3 15.5
306.1 60 1.3 120
5
6
7
8
9
10
11
FIGURE
90
5
APPARENT DECLINATIONS OF STAR-PAIRS
For the upper transit, meridian of La Pla.ta
Greenwich mean astronomical dates
Group X
55
,
0
10
11
12
13
14
42.25
42.34
42.41
42.45
42.47
15
16
17
18
19
42.49
42.51
42.52
42.56
42.64
20
21
22
23
24
42.72 11
42.83 12
42.95 12
43.07 9
43.16 7
25
26
27
28
29
43.23
43.29
43.33
43.36
43.39
30
1
2
3
4
5
7
CI
5
8
9
q
9
9
7
4
2
2
5L!·.44
54.49
5455
54.63
54.71
54.79
54.87
54.93
54.96
54.98
55.00
55.00
55.00
55.03
55e09
,
8
8
B
8
6
3
2
2
16.82
16.85
16.91
16.96
17.04
17.11
17.18
17.23
17,26
17.27
17.27
17.26
17.26
17.28
17.32
5
B
7
7
5
3
1
0
41.37
41.39
41.43
41.47
41.53
41.59
41.65
41.70
41.72
41.. 72
7
41.71
41.69
41.67
41.68
41.71
55r16 11
55·27 11
55,38 11
55.49 8
55.57 7
17.39 9
17.48 10
17.58 10
17.68 B
17.76 6
41.76
41.84
41.93
42.02
42.09
55.64
55.69
55.73
55.75
55.78
2
3
3
17.82
17.87
17.90
17.92
17.94
42.15
42.19
42.22
42.23
42.24
43.43 6
43-49 7
43.56 B
43.64 10
43.74 11
55.81 5
55.86 6
55.92 8
56.00 9
56.09 10
17.97
18.01
18.07
18.13
18.22
43.85 11
43.96 11
44.07
56.19 11
56.30 11
56.41
18.31 11
18.42 10
18.52
2
1
4
8
8
6
4
3
:3
4
0
0
3
6
7
5
.,
FIGURE
91
1
0
2
4
5
3
2
2
3
6
.,
6
6
q
q
42.26
42.30
42.34
42.40
42.48
CI
"
2
4
4
6
6
6
5
2
0
1
2
2
1
3
5
8
q
q
7
6
4
:3
1
1
2
26.21
26.21
26.24
26.27
26.32
26.36
26.41
26.45
26.46
26.45
26.43
26.39
26.37
26.35
26.37
26.41
26.47
26.54
26.62
26.69
26.74
26.77
26.79
26.80
26.80
q
26.AO
26.83
26.86
26.90
26.97
42.57 9
42.66 10
42.76
27.. 05
27.13
27.22
(j
4
6
8
I
-3451
-3451
"
3
6
,
0
-3454
"
5
6
,
0
-3457
"
"
9
6
,
60
59
58
57
-'3452
41.85
41.90
4198
42.07
4216
5
6
7
8
7
0
-3451
1947
6
56
"
0
3
3
5
4
5
4
1
1
2
15.35
15.34
15.. 34
15.35
15.38
15.41
15.44
15.47
15.46
15..45
2
4
15.41
15.36
15.32
15.28
15.28
6
7
8
7
5
15.30
15.34
15.40
15.46
15.52
4
2
2
)
2
]
0
0
3
3
4
7
8
8
9
15.56
15.59
15.60
15.60
15.59
15.58
15.59
15.62
15.65
15.70
15.77
15.85
15.93
1
0
1
3
3
3
3
1
1
4
5
4
4
0
2
4
6
6
6
4
,
1
0
1
1
1
3
)
5
7
8
8
OPPOSITION
1947
3 4
h
12
12 12
20 26 12
28 12
4 5 12
13 12
12
20 27
28
4 5
13
h
36.0
- 6 23
29. 7 ~.~ - 5 59
22.5 7'8 - 5 28
14.7 7'7 - 4 52
07.0 7'0 - 4 14
00.0 . - 3 37
24
31
36
38
37
12
12
12
12
12
12
39.1
+ 8
33.8 ~.! + 9
27.4 7'0 + 9
20.4 7'0 +10
13.4 . +10
07.2 6.2 +10
217 Eudora
h
3 4
12
12 12
20 27 12
28 12
4 5 12
13 12
39.7
+ 0 32
34.9 ;.: + 1 27
29.3 6'1 + 226
23.2 6'0 + 3 26
1 7.2 5' 9 + 4 24
11.3 . + 5 17
1486 Marilyn
h
3 12
20
28 28
4 5
13
21
12
12
12
12
12
12
m
55
59
60
58
53
3
4
12
20 28
28
4 5
13
12
12
12
12
12
12
44
50
52
49
40
3 4
0.260*
1
4
12 44.3
12 12 38.6 ~. ~
20 28 12 31.5 8'0
28 12 23.5 8:2
4 5 12 15.3 7 7
·13 12 07.6 •
- 8 53
- 8 30
- 7 56
- 7 13
- 6 24
- 5 34
20
28
5
13
21
12
12
29 12
12
12
12
260
0.514,
-3. 6
4
0.357*
2
12
20 12
28 29 12
5 12
13 12
21 12
34
43
49
50
18
5
0
85
0.510 ,
-37
f!16
0.360*
1
0
-12
-13
-13
-13
-13
-12
50
06
12
09
01
50
M
6
3
8
11
m
h
3 12
0
235
0.376
-6.' 5
20
28 29
4 5
13
21
0.139
12
12
12
12
12
12
3 12
20
28 29
4 5
13
21
12
12
12
12
12
12
Phryne
40.2
35.1 5.1
29.4 ~. ~
23.6 5'5
18.1 •
13.3 4.8
-11 15
-10 33
-
54
56
7 00 54
-
o
3 12
20
28 29
4 5
13
21
0.098*
12
FIGURE
92
7
12
12
12.
12
12
12
0
229
0.507 ,
-4. 3
0.347
13
11."7
0
44.6
-12 02 33
180
37.5 ~.! -11 29 42 0.431,
29.7 7'9 -10 47 49
-6. 5
21.8 . - 9 58 52
14.4 7.4 - 9 06 50
0.232*
07.8 6.6 - 8 16
12
hm
272
0.352
-6.'7
42
9 44 49
- 8 50
- 7 54
m
13~8
o
m
1549 Mikko
0
0.180*
2
13~3
0
1530 1938 SG
h
0
13
0.399,
-9. 8
0
41.8
-10 15 44
307
35.6 ~.~ - 9 31
0.324 ,
28.4 7'4 - 8 36 !~
-107
21.0 . - 7 33
lef.16
14.1 6.9 - 6 28 65
08.4 5.7 - 5 26 62 .0.046*
1
1291
16.m2
23
m
414
36.3 7.1
28.3 :.~
2~3 •
12.8 7.5
06.3 6.5
h
3 12
0
41.4
-16 38 7
220
36.0 ~.: -16 31 17 0.520
,
30.1 6'2 -16 14 26
-47
23.9 6'4 -15 48 34
6'!'6
17.5 . -15 14 38
0.368*
11.6 5.9 -14 36
7
o
55
44
32
929 Algunde
0
m
39.6
+18 45
34.0 ~.: +19 40
28.1 5'9 +20 24
22.2 . +20 56
16.7 5.5 +21 14
12.0 4.7 +21 19
h
3 12
15!"5
o
cr.7
o
m
.385 nmatar
13
1496 1938 SAl
h
~3
15om 6
o
39.7
- 4 19
32.5 ~.! - 3 35
24.6 8'0 - 2 45
16.6 . - 1 53
09.1 7.5 - 1 04
02.7 6.4 - 0 24
m
12
20 12
28 12
4 5 29 12
13 12
21 12
214
0.448 ,
-66
13~8
1485 Isa
h
3 12
0
0
32 39
346
11 36 0.381,
47
-7. 6
17 30
37 2~
0.153*
44
2
o
m
511 Davida
h
12~0
o
m
Misc.
14~5
o
m
278 Paulina
3 4
1947
Misc.
917 Lyka
69
EPHEMERIDES
o
15?'2
0
45.3
+ 5 47 53
85
38.1 7.2 + 6 40 4'1
0.349
30.2 ~.: + 7 29 39
-6.' 8
22.3 7' 1 + 8 08 26
15.2 5'9 + 8 34 12
0.095*
09.3 . + 8 46
12
13
ELEMENTS
No.
m
m
g
Epoch OIlUT
M
0
m
!J
w
i
1950.0
0
0
0
~
0
n
"
a
158,150 170.066 16.101 6.262
41.999 33~805 1~215 13.647
154.127 3.44.283 8.028 9.822
26.630 12.980 4.433 10.421
7.769
1~440 341278 1~666
639.412
649.958
871.097
630.418
682.398
3.1344
3.1003
2.5505
3.1641
319.949 8.683 12.681
286.257 10.088 4.475
294.668 .9.374 6.790
166.265 4.147 2.026
20.990 13.088 14.763
852.425
736.802
674.113
653.676
656.500
2.5877
2.8517
3.0259
3.0886
3.0797
611
612
613
614
A 615
12.3 8,4 1939 II 22 51.990 25~359 19~333 1~407 6.973
9 24.381 116.310 205.785 20.492 15.462
14.6 10.4 1906 X
13.0 9.3 1945 VIII 13 254.689 63.224 355.185 7.670 3.583
13.7 10.4 1919 IX
3 299.931 206.582 218.019 6.996 6.299
0 270.150 242.400 14.650 2.770 6.390
12.6 9.4 1900 I
689.747
636.959
711.389
802.264
831.146
2.9800
3.1424
2.9192
2.6944
2.6316
616
S 617
618
619
620
12.7 9.7
12.6 5.9
12.4 8.2
12.1 9.2
13.6 10.9
3.410
8.130
4.739
4.370
7.660
869.943
299.717
623.700
886.799
933.328
2.5527
5.1943
3.1868
2.5203
2.4358
621
622
623
S 624
625
13.9 9.8 1942 V
1 129.029 30.373 67.485 2.357 7.883
12.8 10.1 1917 IX 15 329.798 254.035 142.956 8.641 14.032
12.8 10.0 1900 I
0 111.540 123.030 309.120 14.170 6.550
13.2 6.4 1940 XII 19 293.458 17~001 34~152 18.267 1.613
12.1 8.9 1950 I
0 168,883 198.759 12~748 1~O93 12.977
641.457
945.316
919.333
304.721
823.989
3.1277
2.4152
2.4604
5.1373
2.6469
859.549
718.676
855.232
641.364
834.894
2.57.33
2.8995
2.5820
3.1280
2.62.37
760.172
816.653
676.596
666.462
637.307
2.7929
2.6626
3.0184
3.0489
3.1413
17 178,170
20 71.960
0 182.506
25 67.547
24 38.853
601
602
603
604
605
12.6 8.5 1943 III
12.1 ao 1938 I
13.9 10.9 1900 I
12.4 8.2 1945 I
12.9 9.0 1937 X
606
607
608
609
610
0 134.016 54.312
12.9 9.8 1900 I
12.6 9.0 1918 II 19 280.374 288.700
14.1 10.2 1942 VII 20 297.605 67.437
31 112.365 122.275
12.8 8,8 1934 V
15.6 11.6 1942 V
1 209.601 350.997
626
627
A 628
629
630
11.4 8.4
·13.1 9.3
12.2 9.2
1.3.8 9.7
13.5 10.3
1900
1940
1945
1900
1900
1922
1933
1900
1946
1950
0 49.570 10~860 356~60 1~000
I
9 353.645 303.410 43.934 22.103
X
VII 20 307.744 244.184 111.473 17.007
I
0 142.400 174.600 188.420 13.740
0 129.740 333.090
0.860 7.770
J
XII 15 3~555 42.387 342.021 25.454 14.041
V
21 293.478 177.613 143.053 6.449 3.383
I
0 294.050 201.720 11~690 11.541 2.453
I
20 23.677 31.984 87.746 9.322 8.819
I
0 35.080 37.780 10~710 13.900 6.500
3~0013
631
632
633
634
S 635
12.3 8.8 1921 IV
14.5 11.3 1900 I
12.9 9.0 1937 II
13.1 9.1 1937 II
12.6 8.5 1944 I
28 68.920
0 97.310
1 153•.372
1 113.825
19 58,589
636
637
638
639
640
12.4 8.7 1937 III
14.0 9.8 ·1941 III
13.5 10.1 1943 III
12.1 8,2 1943 V
13.0 8.8 1936 XI
10 18~976 296.589 35.429 7.939
21 12.612 164.666 357.000 0.324
28 337.124 127.194 103.655 7.708
25 271.534 65.216 280.666 8,559
13 149.283 18.737 236.156 13.374
9.975
6.799
9.182
6.327
3.923
714.847
631.934
784.808
678,516
630.670
2.9098
3.1591
2.7342
3.0127
3.1632
1 30.102 16.416 41.106 1.733
14.5 12.3 1925 I
13.5 9.3 1934 I
13 122.898 110.525
7.776 8.193
13.9 9.4 1945 XII 16 352.206 21~628 255.110 11902
13.1 10.0 1900 I
0 66.090 267.510 108.940 1.040
0.990 7.063
13.5 9.3 1939 V
5 119.839 87.782
7.392
8.481
5.485
8.980
9.920.
1072.666
627.618
581.259
846.504
623.667
2.2200
3.1735
3.3401
2.5997
3.1869
641
642
S 643
644
645
646
647
648
649
650
14.5
13.5
13.1
15.1
14.7
12.1
10.8
8,9
12.1
11.9
1950
1939
1938
1900
1907
I
I
IV
I
X
276.906 225.621 18.8.30 4.811
247.060 358.770 2.262 11.080
188.894 14~937 1~906 4.947
220.407 134.084 12.288 10.494
22~295 18~961 1~979
4.688
0 346.145 36.155 303.047 6.945 12.327 1000.813 2.3251
15 24.619 173.639 254.826 7.299 10.956 928. 740 2.4439
7 72.117 168.389 292.933 9.989 13.102 628.555 3.1704
0 44.768 347.088 357.887 12.680 16.053 87L566 2.5495
3.071 175.990 216.357 2.555 10.770 918.478 2.4620
5
FIGURE
93
8
94
SCIENTIFIC
through the typewriter would permit the superposition of
several curves, if desired. This problem has been solved
in a general way, but some details remain to be worked
out. It will be possible to process data in a Type 602, so
that the resulting punches will operate in the typewriter
to produce the desired plots. Possible uses of such a curve
plotter, with accuracy of this order, may occur to some
of you.
This has been a very brief summary of the machines
and methods in use at the Nautical Almanac Office. As I
have indicated, the scope of our work is continually increasing. There is considerable satisfaction to be derived
from the fact that every astronomical problem to which
the machines have been applied so far has been successfully solved. I have no doubt that this will continue to be
the case.
REFERENCE
1. W. J. ECKERT and R. F. HAUPT, "The Printing of Mathematical
Tables," MTAC, II (1947), pp. 197-202.
DISCUSSION
Dr. Eckert: I would like to point out that copy for the
Air Almanac has been prepared on the typewriter for
COMPUTATION
three and a half years. In that time about ten million figures were computed, typed, proofread as Mr. Hollander
described, and printed in editions running as high as two
hundred thousand copies. Thousands of aviators used
these volumes day after day. And to date there has not
been one single error reported.
Mr. Hollander: Regarding the accuracy of the typewriter, it has never been known to replace one figure with
another. The worst it has ever done is fail to read the hole
in the card and thus leave a blank space; this can be detected by inspection.
Dr. Eckert: This proofreading business may sound laborious, but you must remember you are putting out a
publication that is going to be studied, used, and sweated
over; whether the production of a figure takes ten seconds
or eleven doesn't make too much difference. One girl,
completely inexperienced in technical matters, holding the
lowest grade clerical rating in the Civil Service, can learn
in two weeks to proofread fifteen pages a day by punching
it up. When one error is so important, you certainly can
afford these checks.
Programming and Using the Type 603 -405 Combination
Machine in the Solution of Differential Equations
GEORGE
S.
PENN
Northrop Aircraft, Incorporated
arbitrary input functions, and two of everything like the
f2(Z) in the case just described. We needed something
that would read factors from the accounting machine on
one cycle, multiply, and enter the product in the accounting machine on the succeeding cycle. For example, a Type
603 Electronic Multiplier hooked to the accounting machine could do the job. In addition, we needed a great
many selectors for controlling counter connections; a
means for converting counters for true figure read-in to
the multiplier; and a balance test impulse available for
every cycle for algebraic sign discrimination.
We asked the local IBM representatives for such a
machine. Mr. A. B. Kimball of IBM was called in. He
suggested, on first sight, going after the answer "hammer
and tongs" style. For some possible time saving ideas he
referred us to an account of step-by-step integration done
at the Thomas J. \Vatson Astronomical Computing Bureau
in New York by Dr. Eckert. As the problem would require
more than 200,000 steps with at least six function products
punched, carried to the multiplier, to the collator, and
then back to the accounting machine for each cycle, this
was prohibitive. The problem was then taken to IBM
Headquarters in New York where Mr. J. C. McPherson
went over the problem requirements in detail and agreed
with our conclusion with respect to machine requirements.
The machine was built. In fact, we received much more
than we expected. The request was for a bare minimum of
the items noted above. vVe received a machine comparable
in programming technique to current large computer design-a poor man's ENIAC. Only eighty decimal digits
of memory, only six-by-six multiplication, and only 150
computing cycles per minute, but sequence controlled completely. Our experience indicates this program power is
fundamental.
A description of specific machine features is not amiss.
First, I assume you are familiar with the IBM manuals
on the Type 603 Electronic Multiplier, Type 405 Accounting Machine and Type 517 Summary Punch. Just about
all available extras are included.
FOR THE PAS T two years Northrop Aircraft has
used a small installation of International Business Machines equipment for engineering calculations. While much
of the work was and is routine-stress distribution, reduction of wind tunnel results, and the multiplication and
inversion of matrices-considerable miscellaneous work
in connection with research projects has appeared from
time to time. Last winter two despairing men brought us
a system of differential equations for step-by-step integration.
The differential analyzer at the University of California
was unable to solve the problem because of the presence of
a term a¢
cf> where a is of the order of 107 • The turn
ratio for the ci> and cf> shafts made the time of solution run
into centuries. Large scale digital computers were considered but the earliest schedule times Were remote and the
probability of needing additional work after the first set
of solutions were completed introduced the possibility of
another schedule delay-two years if we were lucky. After
telling us this story they presented us with the system
+
= ax + by - efl (t)
= dy + ef2(z)
Z = fz + g:c + h%
z
y
z>
+~
-~ 0,
v(x,t)~O
as
x~oo
(3)
= l/y;
(4)
(5)
It is desired to obtain solutions of (1) corresponding to
assigned values of the constant m and subject to the
boundary conditions (2) through (4). The usual numerical procedure consists of replacing (1) bya difference
equation and finding the solution of this equation as an
approximation to that of the partial differential equation.
Putting
x
= i6.z,
v(i,j)
t
= j6.t
(i,j
== v(i6.s,j6.t) ,
pj
== p(j6.t),
Letting b
i
CX)
y-; + 2m
Y1r,p(t)
]
.
(8)
where
y-;'
Pf) = --'----
Y1r+2m
'l!(i+l,j)-'l'(i-l,j)
26.x
.
~(I-~m6.xpj) v(i-l,j)
+ ~(I+~m6.xpj)v(i+l,j).
1[
oVt(x,t) dx =2. 1 -
That is,
== 6.t/(6.X)2 = !, this equation may be written
v(i,j+l) =
(7)
In view of the boundary conditions (2), (4), and (3), this
may be written
v( i-l,j) - 2v( i,j) +v( i+ l,j)
(6.X)2
+ tnpj
+ 6v(3,j)]
The solution of (6) with boundary conditions (7), (2),
( 4 ), and (5) leads to a function v (i;j) defined over a
rectangular network of points. The procedure for solving
consists in starting with the function (5) defined for the
row j = 0, and with Po = 0, and proceeding by means of
(6) to the row j = 1. The quantity PI is then determined
by (7) and is used in (6) together with the boundary condition (4) to determine v for j = 2, etc. A check on the
procedure is derived in the following manner: both sides
of (1) are integrated with respect to x between the limits
o and 00, giving
= 0,1,2, ... )
and replacing the partial derivatives in (1) by finite difference approximations, the following difference equation
is obtained:
v( i,j+ 1) -v( i,j)
6.t
- 9v(2,j)
(1)
(2)
-2fJv(0)/fJs = 1 - pet)
v(O,t)
11/y;-
+ [18v(1,j)
,
m = const.)
COMPUTATION
is the steady state value of p. Integrating the two sides of
equation (8) between the limits 0 and t, interchanging the
order of integration on the right, and making use of the
fact that
1
o v(.'r,O)- dx = 2' '
f.
we obtain
(6)
CX)
FORUM
103
PROCEEDINGS
By computing fT(t) in the two ways shown in the foregoing equation, a numerical check is obtained.
The problem described is well suited to a machine of
the type of the IBM Relay Calculator. The machine
method conforms to the hand computation procedure in
that a row of points is computed at a rttn. For any row
j = constant, the input cards in the reproducing feed consist of cards punched with the values v (i,j). In the first
four of these cards are punched the constants appearing
in (7) : 3.6% - l1/y-;, 18, -9, and 6, respectively. At the
beginning of each run the constants
~(1 - ~m.6.r,,"j), ~(1 + ~'m.6%,,"j)
(9)
are set in two rows of switches. A stack of blank cards is
fed into the punch feed. As the input data are fed through
the reproducing feed, the machine uses the values v(i,j)
to compute the values v( i,j + 1) for the next row by (6),
and punches them into the cards going through the punching feed. Using the first four values v(O,O), v( 1,0),
v(2,0), and v(3,0) the machine also CO!:l1putes 3,,"1.6.V by
(7), and the constants (9) to be used in the next run. At
the end of the run the output cards are removed from the
punch stacker, the constants (9) are read and set in the
switches, the cards are placed in the reproducing feed
while another deck of blank cards is placed in the punching
feed, and the process is repeated.
The interval .6 t must be chosen small to begin with, for
,,"(t) behaves like 'Vi for small t, so that ,,"'(0) is infinite.
As the solution progresses and the steady state condition
is approached, .6 t can and should be increased to keep the
computation time within reasonable bounds. It must be remembered that in order to satisfy the condition b = 1/2,
.6% must be changed whenever .6 t is. In the problems carried out on the relay calculators, the interval .6% was
doubled whenever it was found that the results obtained
upon using the larger interval were the same as those obtained by use of the smaller interval. Solutions were carried out for the cases m = 0, 1, 2, 5, and 10.
REFERENCES
1. W. J. ECKERT, "The IBM Pluggable Sequence Relay
lator,') MTAC III (1948), pp. 149-61.
2. J. LYNCH and C. E. JOHNSON, "Programming Principles
IBM Relay Calculators," Ballistic Research Laboratories
705 (1948), Aberdeen Proving Ground, Department
Army.
Calcufor the
Report
of the
DISCUSSION
Dr. Herget: I have found it pretty dangerous to use
b = 1/2. Haven't there been papers showing that the solution may not be stable for the value one-half?
Dr. Levin: The equation we are dealing with is of parabolic type. It may be shown that for such an equation the
condition for stability is .6t/ (.6x) 2 < 1/2. Also a big
reason for using b = 1/2 is that it yields a simpler formula.
Dr. Thomas: May I take two minutes for a statement
about this sort of problem?
Suppose you have fixed ranges and even intervals in .r.
What you are essentially doing in this sort of work in
going to finite differences is exactly equivalent to using a
finite number of terms in a Fourier expansion. You have
terms, for instance, in sin ;r, sin 2x, sin 3x, . . . ; these in
the exact solution should decay with time at rates e- t , e-4t,
e-rJt , • • • The difference formulas, besides only giving a
finite number of terms of the expansion, make the terms
given decay at the wrong rates.
The later terms decay at less and less accurate rates. In
the simplest case, you can show that if you have n points
corresponding to n terms of the expansion, the first few
terms will decay nicely; but as you come to higher terms,
the end terms decay more and more slowly, and it is usually
the case that the last term you keep decays a little less
rapidly than the first. So all you can do is get the limiting
form; you get nothing exact about the details after some
time.
You have to look at this sort of thing fairly closely to
estimate what is happening, especially if you hope to extend it to more complicated cases with curved boundaries,
and so on.
GENERAL DISCUSSION
Dr. Eckert: We have a few minutes to wind up the
Endicott part of our conference. With reference to the
last papers on special equipment, I might say that you will
see two relay calculators of the type Dr. Levin and Dr.
Bramble mentioned, and a prototype combination machine
somewhat like Dr. Fenn's, at the Watson Laboratory. The
card-operated typewriter at the SSEC is also on display .
Are there other comments or questions?
Mr. Kintas: We have been considering solving natural
frequency problems for free crank systems and compound
beams. Mr. Mack suggested that both problems could
probably be handled by a single setup of matrix equations
using the coefficients and unknown terms which describe
the conditions of kinetic and potential energy in the system. Other aircraft men here have investigated that
problem.
Mr. lIarntan: What is the maximum capacity of the
Type 602 in multiplying square matrices of the nth order?
Suppose the elements are two-digit numbers; how many
can be handled at once by the 602?
Dr. BraNtble: If you put one element on a card there is
no limitation on n.
Dr. Grosch: But then multiplying two unsymmetrical
matrices together would require 2n 3 cards to pass through
the 602. I f there were 'an unlimited number of program steps available, one could store perhaps thirty-two
104
two-digit numbers in the machine, reserving only counters
1, 3, 6, 9 and 10 for operational use. I am assuming all
elements are positive. Then every other card passing into'
the 602 would carry up to sixteen elements of matrix A
(one row) and the next would carry up to sixteen elements of B (one column) ; one element of A· B would be
accumulated in LHC 9 and 10 and punched on the B cards.
Only 2n2 cards would pass into the 602, for tz. < 16. But
since fifteen recalculation cycles would be required, the
standard machine would be limited to smaller values of H.
And many selectors would be needed, especially if negative
elements did occur.
I would like to make a statement about the bibliography
on technical applications of punched card methods we put
SCIENTIFIC
COMPUTATION
out at the Watson Laboratory in 1947.1 If any of you have
not seen it, you can secure copies through the IBM Department of Education in Endicott. \Ve want to correct, im~
prove and expand that first attempt. In the foreword we
requested suggestions as to format, and especially addi~
tions to the list of references. I am sorry to have to report
that in the first year not one single item has come in except
through our own efforts. If some of you will help, we
might revise and expand the bibliography annually.
REI"ERENCE
1. International Business Machines Corporation, "Bibliography
[on] the Use of' IBM Machines in Scientific Research, Statistics, and Education," Form 50-3813 (1947).
Simultaneous Linear Equations
FRANCIS
J.
MURRAY
Columbia University
to give about three figures in the result. Eliminating Y
yields an equation in X
ISH 0 U L D L IKE to discuss linear equations and the
solution of linear equations more or less in connection with
the type of machines that we have developed. I am sure
that you are all familiar with many other approaches, and
the mathematical basis of what I would like to say is reasonably well known. The usual theory of linear dependence
such as ,one finds in Bacher, and the relatively simple
vector constructions, are to be considered not absolutely
but to within a certain accuracy. It seems to me that this
has to be done as soon as one considers systems in which
the size of the results is not obvious.
Now let me briefly describe our objectives i11 building
the machine which is now in use at the \Vatson Laboratory.
We set up a machine to solve a system of equations
o.aa 1999 X =
1.000000
where, due to rounding, the last 9 in the coefficient of X
may be off by more than one. The first three digits for
which the matrix of coefficients is singular also corresponds to a loss of accuracy in the answer.
We have a reasonably routine method of handling problems of this type on our three-digit machine. Mathematically this process consists in looking for a linear dependence among the columns of the coefficient matrix and by
a suitable change of variables eliminating them. Notice that
if in our example we let X = (.'r + y) and, Y = ..1," - y
our equations become
(1)
+
1.999000 x
0.001000 y
1.999000 x - 0.001000 y
However at all times the device regards the x/s as well as
the ai/s and b/s as inputs and evaluates by relatively simple
means
(2)
By a change of y scale, i.e., y
=
1.000000
= 0.000000
= 1000 z this becomes
1.999 x + 1.000 z = 1.000
1.000 x - 1.000 z = 0.000
which it represents on a meter. One finds a solution of the
equation by manipulating the x/s so as to minimize p.. This
manipulation process is always convergent, an advantage
enjoyed by no other adjusting device. If there is a solution,
we obtain it.
Another objective was to put in the coefficients digitally,
and we succeeded in doing this. The use of alternating
current raised certain phase difficulties and these Mr.
Walker settled. 1
In attacking a system of equations (1), we begin by
rounding to three figures. Normally 'any accuracy can be
obtained by an iterative process which involves merely a
change of the constants b'i at each step. However, this is
clearly ineffective when dealing with a problem of this
sort:
1.000000 X + 0.999000 Y = 1.000000
0.999000 X + 1.000000 Y = 0.000000 ,
This problem can be readily solved on the machine. The
coefficient of z now has only three significant digits but
these are all usable in the machine and we can obtain all
the accuracy the problem justified.
The advantage of this type of device is that the essential
linear combinations can be found on the machine. Consider
our original system of equations
(~1:l
l ·t 'l
+ ... + a
12 1:l X 12
= b12 t
In the machine, the b/s appear multiplied by a gauge variable t. This has the result that the unknowns Xi are in the
form xi/f and hence not necessarily between -1 and 1.
Suppose that the third column is very nearly a linear combination of the remaining columns. We can readily find
this on the machine as follows: We set t = 0, Xs = 1 and
minimize
p. = ~i (ah ·1,"1 + ai:l x:! + ais + ... + ail!! X 12 ):!
+
since the problem is indeterminate relative to the first three
digits. Nevertheless our machine can be effectively used
here. Before discussing this method, let us point out that
the six figure accuracy of the coefficients is only adequate
relative to
105
.1,"1' X 2 , X,H ••• ,
X 12 •
This yields the coefficients
106
YH
SCI E NT I FIe
Y12' of the linear combination and we can use
these to compute digitally the terms
)'2' •.. ,
In this computation all the digits of the a/sare used. We
then can rewrite our .original system of equations
al
I (Xl -
a121
(%1 -
Y1 X
g)
+a
)'l.t"g)
+a
+ a l 2 (%2 112
Y2-t ·g) + a 1 3 Xg
(X12 - Y12 .t"g) = bIt
+ au 2 (X2 . . .,. Y2%g) +
12 12 (%12 --
a12
+ .,.
g Xg
+ ...
Y12%3) = bi2 t .
vVe now have a system of twelve equations in the twelve
unknowns
:z1 =
%1 -
Yl·t"3' :z2 = %2 - .Y 2%g, :zg =
Z12 = X 12 - Y12 %g •
%g , . . . ,
The a column has been obtained by a minimizing process
relative to the other columns, and if this has been carried
out completely the a column will be orthogonal to the
other columns. The a column will have smaller quantities
in it than the aij and consequently fewer digits. In these
circumstances, however, we can always change the scale
of %g and, in the machine, utiliz~ the values of aia to within
0.1 per cent of the largest aia' The loss of accuracy represented by the fewer digits in the .other a'S corresponds
precisely to loss of accuracy in these equations themselves,
i.e., the loss of accuracy involved in passing from coefficients to unknowns. On the other hand, we have eliminated a linear dependence among the coefficient columns
and thus permitted the solution process to proceed.
When there is only one approximate linear dependence
among the columns, the above process permits us to eliminate it by modifying just one column of the coefficient
matrix, i.e., just changing one card. I f there are a number
of these dependencies, a number of these steps must be
taken. Theoretically, this modification process can be
pushed to an extreme. vVe begin by minimizing the second
c.olumn relative to the first. This yields a second column
orthogonal to the first. Then we can minimize the third
column relative to these two. This process can be repeated
until we have a matrix whose columns are orthogonal, and
this can be inverted immediately with practically no loss
of accuracy. Let A be the original matrix, T the new
orthogonal matrix. We have found a matrix C such that
AC = T and hence A-I = CT-1. Since the columns of T
are mutually orthogonal, the inverse T-1 of T is very
simply connected with its transpose. Thus if T is the
matrix 1tij r where i is the row index, j the column index
and Tj = 'tit'i}' then T-1 is simply 1tj,i/Tjr. Note that the
accuracy of T-1 is essentially that of T and the loss of
COM PUT A T ION
accuracy for the system is represented by the change of
scale factors in C.
These difficulties are present of course in every method
of solving simultaneous linear equations since we are
always limited in the number of digits available. Even in
the elimination process one must always keep in mind
the retention of maximum accuracy. For instance, the
von Neumann-Goldstine estimate of the loss of accuracy
involved in the elimination process is based on the assumption that the largest coefficient available is used as a divisor
in each step.!! Failure to follow this as a policy can lead to
a serious loss of accuracy. If we merely choose a coefficient which is half as big as the largest, in the twelfth
order case at each step we will have lost unnecessarily an
accuracy factor of 211 = 2048. Pyramids have been built
on factors like this.
Finally, one should point out that the type of quasisingularity I have discussed above is by no means uncommon. In the majority of practical problems, full digital
regularity is exceedingly difficult to obtain. Instead some
aspect of high accuracy may be reasonably available. For
instance, t.o locate a point in a plane by linear observation,
two observation points far apart could be used, lines from
each of these points drawn, and their intersection found.
This would yield full digital regularity in the corresponding linear equation problem. However, this is not always
practical. Instead one may have to take the observation
points near together and very accurately measure the direction angles of the lines.
This situation generalizes. Full digital regularity is frequently nonobtainable and we must substitute high accuracy procedures instead. This leads to these quasi-singular
situations.
REFERENCES
1. R. M. WALI
FIGURE 3. PROBLEM 3 :
<
~' AS A FUNCTION OF
~
W
Cl
W
~
U
W
~
~
1.1
W
~
C/)
C/)
0
~
U
<
0
H
~
<
~
W
,~=o.o
1.0
~
::>
C/)
C/)
W
~
Po!
II
IJJ..J)
9
0
4
12
8
16
W
20
24
= ANGLE
28
32
36
40
44
48
52
56
60
64
OF INCIDENCE IN DEGREES
TABLE II
Problem
~
Wt
3
1.0
0.9
0.7
0.5
0.3
0.1
0.0
43.089
44.090
46.442
49.429
53.392
59.009
62.833
0
WI
0
90.000
74.146
65.457
61.945
61.439
64.272
67.792
WL
Problem
~
Wt
1.0
0.9
0.7
0.5
0.3
0.1
0.0
46.911
48.059
50.779
54.295
59.110
66.410
72.009
0
0
69.295
69.295
69.295
69.295
69.295
69.295
69.295
15
WI
WL
0
90.000
74.146
65.457
61.945
61.439
64.272
67.792
unreal
unreal
unreal
unreal
unreal
unreal
unreal
FORUM
111
PROCEEDINGS
100
~=0.9"
95
l; =0.7",
90
" '"
~
85
~ =0.5
80
'" \
"" \
\
\
I
\
\
\
I
\ \
~ =0.3\\ \
~ 75
\~
~
~
~=1.0\
~~
70
\
\
~
Cl 65
Z
..... 60
Z
9f-I 55
U
< 50
P:::
~
45
P:::
~40
~
~
35
II
30
<
FIGURE
"325
4.
PROBLEM
3:
w' AS A FUNCTION .oF w
20
15
10
5
W
10
20
= ANGLE
30
OF INCIDENCE IN DEGREES
40
50
60
70
80
For Problem 3 a reflected rarefaction wave solution exists
at all values of ~ at normal incidence (w = 0) ; whereas
with the gases interchanged (Problem 15) there is always
a reflected shock-wave for n.ormal incidence.
The dependence of the strength f of the reflected wave
on the angle of incidence w for incident shock strengths
~ = 1.0, 0.9, 0.7, 0.5, 0.3, 0.1, 0.0 is plotted in Figure 3.
A pressure ratio greater than 1.0 indicates a reflected
shock-wave; while a value less than 1.0 indicates a reflected rarefaction wave. We see, as we said before, that
for all strengths of the incident shock wave, the type of
c.onfiguration is of the rarefaction variety for normal incidence (w = 0).
A typical curve begins at a value of f less than 1, and
continues upward until it reaches Wt, the transition angle,
90
100
at which point the configuration changes from a rarefaction
to a shock-wave variety. It extends upward until it reaches
an extreme angle WL beyond which no further solutions are
possible. The existence of an extreme angle was first demonstrated by J. Von Neumann f.or the problem of regular
reflection from a rigid wall. The similarity of the two
problems is here illustrated.
Figure 4 is a plot of w' (the angle of refraction) against
w (the angle of incidence) for the same problem. For weak
incident shock-waves Snell's law of refraction holds. This
is the relati.on used for plotting the curve for ~ = 1. It will
be noticed that for any other shock-wave strengths (including infinite strength) the law of refraction does not
differ appreciably from Snell's law. This is an interesting
result.
112
SCIENTIFIC
COMPUTATION
O.960~------~1~O~----~2~O~------~3~O--------4~O--------~5~O--~----~------~7~O--------~80
W =
FIGURE
5.
ANGLE OF INCIDENCE IN DEGREES
PROBLEM
15: f
AS A FUNCTION oF' w
For this problem all configurations are of the reflected shock-wave variety for
normal incidence, which usually change over to the rarefaction type at larger values
of w. For ~ = 0.0 and 0.1, however, no transition takes place, as has been pointed
out previously. In that case, the curves reach an extreme angle WI., beyond which no
further solutions are possible.
FORUM
113
PROCEEDINGS
70
c;=l.O
65
60
U)
~
~
55
0::
t:J
~ 50
q
Z 45
Z
~
0
~
~
l)
~
~
~
0::
~
40
35
30
0
~
~
25
~
20
t:J
Z
II
J
15
10
5
5
10
15
20
W
30
25
= ANGLE
FIGURE
6.
35
40
45
50
55
60
OF INCIDENCE IN DEGREES
PROBLEM
15:
w' AS A FUNCTION OF w
65
70
75
114
SCI E N T I.F ICC 0 M PUT A T ION
We will next consider two other problems (10 and 22)
which exhibit slightly different properties. In both Problem
10 and 22 there are no real solutions for the transition
angle for values of ~ > 0.66. ,This means that for these
values of ~ there cannot possibly be any transition from
one type pattern to the other. In other words, the refraction pattern which begins at normal incidence must continue throughout the problem.
In the case of Problem 10, the starting pattern is always
of the rarefaction type for values of ~ > 0.66, but of the
\\.~=O.l
£=0.7--- \
1.9
\
E-!
U 1.7
0::
\
\
\
\
\
1.6
\
\
0 1.5
\
1.4
E-!
<
0::
~~
\
0::
u
01-1
10:
~=0.3
\
til
til
<
PROBLEM
\
\
'\
0 1.8
7.
e AS A FUNCTION OF
\~
\
\
~
~
~
~
~
FIGURE
1--\ \
~
<
~
~=0.5
\ \
~=O.9 \
2.0
>
shock wave variety for the case ~ < 0.66, or where the
solution for Wt is unreal. For Problem 22 exactly the reverse occurs.
Figure 7 is a plot of e (the ratio of pressures across the
reflected wave) versus (0. In the case of strong shockwaves we start with a rarefaction type of configuration at
W = 0 and go over to the shock wave variety at Wt, while
for weak incident shock-waves we start with the shockwave variety, which persists throughout the entire range
of w.
\ \ \\ ~=O.O
1.3
\
~
~\ \
\ \
\\
0::
::J 1.2
til
til
~
c;=0.3
~=0.7
~=0.1
0:: 1.1
P-r
?'=0.9
"
~ 1.0
O-S =0.0
0.9 0':----."..---1
W
20
= ANGLE
30
40
50
OF INCIDENCE IN DEGREES
TABLE III
Problem
10
~
1.0
0.9
0.7
0.5
0.3
0.1
0.0
Wt
unreal
unreal
unreal
0
11.820
20.739
31.192
37.558
WI
90:000
73.554
65.378
62.929
64.065
70.061
77.690
WL
0
42.706
43.123
44.009
44.972
46.025
47.183
47.807
Problem
22
~
1.0
0.9
0.7
0.5
0.3
0.1
0.0
Wt
unreal
unreal
unreal
0
16.848
29.476
44.916
55.360
WI
0
90.000
74.526
65.566
61.439
60.005
61.294
63.431
WL
unreal
unreal
unreal
unreal
unreal
unreal
unreal
W
FORUM
115
PROCEEDINGS
140
~=I.O ~
~
130
~
~=0.9 -
120
-
~
__
-.......
~
'-.....
110
{;=0.7 - _
r.J)
ril
ril 100
P=:
~=0.5 -
t:)
ril
0
Z
H
90
Z
0H
80
--- -- -......
--- ----
___
£=0.3-- _
___
-- "
"" '"
'---...
"- \
-........
'" \ \ \
"- \
~=O.l "--~ \ ,
\
\
""-....
\~\
\
~
U
<
P=:
\
70
~
ril
P=:
60
~
0
ril
~
50
t:)
Z
<
40
II
-:)
30
20
10
O~-----------------------------------------------------------
o
4
8
12
16
20
24
28
32
36
40
W = ANGLE OF INCIDENCE IN DEGREES
FIGURE
8.
PROBLEM
10:
w' AS A FUNCTION OF w
A refracted angle of 90° indicates a transmitted wave which is normal to the
interface. In acoustics this angle is known as the angle of total reflection, and is the
limiting point beyond which no refraction pattern can take place. For stronger
shock-waves a refraction angle of 90° does occur mathematically; however, it is
on a portion of the curve which we believe to be physically unreal. \'1 e must t1ms
expect other types of limiting angles for shock waves of finite amplitude.
116
SCIENTIFIC
COMPUTATION
1.10
iii
!>
\
2.0
\
::>
en
en
\
~
'"
ct:;
~
"
-..\41)
~=0.5
1.5
5=0.3
g =0.0
5=0.1
M PUT A T 10 N
shock-wave solution always occurs in Problem 11, while a
rarefaction wave solution takes place when the gases are
interchanged (Problem 23).
In Problems 11 and 23 no real solutions for the transition angles Wt exist at all; hence, the type of pattern which
occurs at normal incidence persists throughout the entire
range of permissible values of angle of incidence. Thus a
4.0
ceo
~ =0.9
""
\\
\\\
\
\
"" ~ ~,
',\
~=0.7
'\
\
~,
1.0
0
10
20
30
40
W= ANGLE OF INCIDENCE IN DEGREES
FIGURE 11. PROBLEM 11 : f
AS A
FUNCTION OF
w
The shock-wave solutions do not exist f.or angles of incidence
greater thal1al1 extreme angle WL'
50
FORUM
119
PROCEEDINGS
IV
TABLE
Problem
11
~
Wt
1.0
0.9
0.7
0.5
0.3
0.1
0.0
unreal
unreal
unreal
unreal
unreal
unreal
unreal
Wl
WL
0
0
90.000
73.554
65.378
62.929
64.065
70.061
Problem
Wt
Wl
WL
0
23
50.769
51.323
52.512
53.821
55.275
56.903
57.795
i7.690
~
unreal
unreal
unreal
unreal
unreal
unreal
unreal
1.0
0.9
0.7
0.5
0.3
0.1
0.0
90.000
74.526
65.566
61.439
60.005
61.294
63.431
unreal
unreal
unreal
unreal
unreal
unreal
unreal
140
130
120
~=0.7--_
CI)
~
~
~
~ 100
Cl
Z
1-1
--- --
110
;=0.5 -
_
?'=0.3
-
- -- -- - -
~=0.1- _
90
E:
~
S=O.O
U
...:::::
-
~-
-
-
-
-
_
-
~
-- --- ___
--
........... """
- - "'"""-"'"
- - __
-:::.:..__
--___
" "'- \
\
-................."
,,"-
"
""
"'- \
'~ \ \
--'":::::::...---..........
--~,
Z 80
~
Cl
H
U 70
Z
H
0
-__
"-
~\
\,
60
~
~
d
Z
<
"
3
50
40
30
20
10
0
0
5
10
15
20
25
30
40
35
W = ANGLE OF REFRACTION IN DEGREES
FIGURE
12.
PROBLEM
11 :w'
AS A FUNCTION OF'
W
45
50
I
120
SCIENTIFIC
COMPUTATION
1.00
/
/
/
~=0.7
/
/
/
I
/~soniC.
y
Line
/
/
, f, =0.0
"
--.oJ)
.75~------------------------------------------------------------------------10
20
60
70
o
30
40
50
(J
= ANGLE OF INCIDENCE IN DEGREES
FIGURE
13.
PROBLEM:
e' AS A FUNCTION OF
(0
FORUM
121
PROCEEDINGS
70
60
~A:;::...----~ =0.1
II
:3
10
10
20
30
W
= ANGLE
FIGURE
14.
40
50
0
OF INCIDENCE IN DEGREES
PROBLEM
23:
w' AS A FUNCTION OF w
70
80
122
To sum up, the refraction problem was formulated in
terms of a complex system of algebraic equations. In
order to solve these equations we first had to carry out a
number of preliminary investigations. First, we had five
independent parameters to deal with. To take into account
all possible variations would probably not be feasible even
on a large-scale· machine, and would also fill up several
large volumes which no one would want to read anyway.
We thus selected typical cases of gas combinations; but in
each case we obtained solutions for all strengths of the
incident shock wave and for all possible angles of incidence.
Secondly, we were faced with the problem of finding out
a priori which ones of the many possible mathematical
solutions were physically plausible. We have solved this
difficulty by connecting our soluti.ons with the known solutions for the acoustic case, and for the case of normal
incidence.
\Ve still had several other difficulties. We had to consider the two types of configurations, and we had to predict in each case which type of configuration will occur.
We also wanted to have a good representation of the various types of characteristic curves that may occur. We
wanted to obtain a sufficient number of solutions which
we can later compare with experimental results. In all, we
solved a total of approximately 3000 points. If a single
point were to be solved by hand, it would probably take an
experienced computer at least one day. If the entire job
were done by hand, it would take one computer approximately fi fteen years or a staff of ten computers about a
year and a half. In both cases, I believe, the problem would
get out of the range of feasibility. It is doubtful if the
pr()blem would ever have been solved without the use of a
SCIENTIFIC
COMPUTATION
large-scale machine. The actual running time on the machine· was approximately sixty hours. Preliminary preparati.ons and coding of the problem for the machine took
several weeks. Thus, within a month we actually were able
to obtain a solution to a problem which otherwise probably
would never have been solved.
I would like to make a few general remarks about the
place of machines in the scientific computation field. I do
not believe that machines will ever replace the necessity
for analysis, or for scientific investigation; however, the
automatic calculating machine certainly represents a new
and powerful tool which the scientist will be able to utilize
in the future in the solution of many difficult and otherwise unfeasible problems. There are really only two largescale electronic machines in this country at present in
actual operation. The EN lAC, at Aberdeen, was the. first
machine to prove that large-scale computations of complicated mathematical problems is feasible. The SSEC is the
second machine. When we used it, it had been in operation
for approximately six months, but it had already demonstrated not only that solutions of. complex mathematical
problems on large-scale electronic machines are feasible,
but that these are actually practicable.
In conclusion, I would like to thank the staff of the
IBM Corporation for the excellent spirit of cooperation
they have shown throughout all phases of this problem,
and f.or the valuable assistance they have offered both in
coding the problem and in preparing the problem for the
machine. I would also like to mention that this work has
been carried out in collaboration with my colleague, Dr.
R. J.. Seeger, and that the work was sponsored by the
Mathematics Branch of the Office of Naval Research.
Computation of Statistical Fields for Atoms and Ions
L.
H.THOMAS
Watson Scientific Computing Laboratory
I PRO P 0 S E not to go into the theory of the problem
that I am going to talk about, but to write down the equations with as little explanation as possible, and show how
they were transformed to what seemed to me to be a suitable form for putting on the machines. Then I will show
you what the people who coded the problem and put it on
the machines gave back to me, and how I had to treat that.
I shall be trying to explain the point of view of the person
who brings a problem to the machine, to try to show you
how much that person has to do and how little he has to
do it with.
I t happens that in this case the problem is very small
compared to the capacity of the machine. I do not need
anything like the number of decimal places that the machine has. I do not need anything like the storage that the
machine has. The problem could perhaps have been done
on a smaller sequence computer. It might have been done,
with considerable inconvenience and taking quite a long
time, on the relay calculators. However, it is much easier
to do it on this machine, and much easier-that is, for mesimply because I do not have to worry about capacity at
all, throughout the problem. I just put the equations there
and say, "do it," and it turns the solutions out. There is
everywhere ample capacity.
2"1 ( P
P
- 2/-rr
= S/2-rr
)2
V'J
atr
tain finite distance ro, P shall take the boundary value S/2-rr
while its derivative satisfies the second condition involvinz Z - n, where n is the degree of ionization, so that
Z - n measures the part of the nuclear change unneutralized by the electrons around it.
We then have a problem in which the outer boundary
is not fixed beforehand. The equation is a non-linear ordinary equation, since this is a purely one-dimensional problem. We have one boundary condition at the center; at
r = ro we have the outside movable boundary. It is known
from much work on this and on more complicated problems that it is convenient to adopt a logarithmic scale of
distance, using a finer mesh where r tends to zero. Therefore, we make a number of substitutions directed towards
a logarithmic independent variable in place of r, while
leaving a differential equation of the second order without
a first order term in it.
The following substitutions were made:
Going to a logarithmic scale, if r = e-:C, '" = em/ 2 X removes
the first order term from the resulting equation. It is convenient also to scale the problem so that the movable
boundary is fixed, writing % - %0 = y, at y = O. Writing
X as shown makes the coefficients of the final equation
simpler.
This gives us the equations in the following form:
Z/r
= ro
-r2 d/dr{i-·f 1 IWEf'X·{r
3
~ ~';,Ig~:
925137
925137
11943290
11943290
: ~: ~ ~I~ ~ ~ ~
7522334
7520034
7517453
7514557
7511307
925137
925137
925137
925137
92513'7
11943290
11943290
11943290
11943290
11943290
154809048
454596259
454357620
454089571
453788526
7507659
7503566
7498972
7493816
74R8030
925137
925137
925137
925137
925137
11943290
11943290
11943290
11943290
11943290
453450900
453072225
452647139
452170246
451635280
7481536
7474247
7466066
7456883
7446575
925137
925137
925137
925137
925137
11943290
11943290
11943290
11943290
11942397
451035058
450360601
449603808
44R754723
447801684
74 3!5 004
7422014
7407431
7391059
7372676
925137
925136
925135
925134
925133
11941503
11940609
11939714
11938819
11937027
446732030
445531272
44418312R
442669439
440970187
73512035
7328857
7302829
727.3598
7240770
925131
925128
925124
925118
925110
11935234
11932543
11928951
11925356
11919959
439062590
436921154
434516636
431817792
4287AR093
95932
9SR71
72013900
7162490
7115981
7063747
925098
925082
925059
925026
11913656
11905544
11895619
11RR2975
4253817096
421569000
4172A4260
412477007
481095
473673
9383
9001
~I~ § ~ ~
51186
46942
40757
36878
32397
2R845
467376
464698
462296
4602R9
8643
8476
8315
8163
671286
646746
624402
604265
4S8596
457292
4S6408
4S5151
6428
43265
41
40104
36944
33781
8019
7R85
7762
7649
44098
41238
38360
35461
32539
34322
31826
29386
26955
24577
26568
24395
22278
20260
18294
586549
578428
570868
563872
586499
578378
5.70818
563822
455189
456178
456525
456186
7549
7503
7461
7422
0613
27437
25046
31
23448
21846
29591
26614
24357
22843
21320
2204
19833
21
18100
16914
15727
16377
14508
13140
12228
11315
557442
551581
546292
541578
557392
551531
546242
541528
457560
458252
459134
460096
7386
7353
7325
7299
240
18629
21°
17013
15391
13762
19787
18244
16690
15125
13548
4585
13441
11
12295
11147
9997
110447
9578
8707
7835
6961
537424
535574
533870
532312
537416
535566
533862
532304
461294
461108
462607
463302
7278
7269
7';161
7253
2535
11715
11
10893
10069
9243
12356
11558
10757
9952
9144
9121
8521
7920
7318
6715
6315
5899
5483
5066
4648
530901
529637
528521
527553
530894
529629
528513
527545
464078
464050
465665
466522
7247
7';142
7:>. 37
7233
8415
7585
6753
5918
5081
8332
7517
6698
5875
5048
6111
5506
4900
4292
3683
4230
3811
3391
2970
2548
526743
526074
525555
525187
526726
526057
525538
525170
467456
468429
469442
4704S9
7';131
7230
7230
7230
4241
3398
2553
524970
524905
524993
525234
524953
52488R
524976
525217
471592
472691
473069
475086
7233
7235
7239
7244
4217
3382
2543
1700
0852
3073
2461
1848
1234
p618
2126
1703
1279
0854
10428
%~~~
FIGURE
1
126
these solutions. It turns out that if you try to make a rough
check on results by identical relation between the integrals,
the check is not nearly accurate enough unless you have
added in these extra terms to the sums. I had hoped that
I would not have to do this for the originally assumed
values of y, but it turns out that to get a good check one
must.
vVell, then, these solutions were made originally for
n = 0, for a wide set of values of y. Gamma equals 1, 2,
3, 4, 5, 6-1 think that is as far as it went. For each of
these we computed the corresponding value of Z. Actually,
that was computed as we went along, and we found that
this covered the whole of the periodic table. We then took
the part which included the whole of the periodic table and
took even intervals among these values of y = 2.1, 2.2,
and so on, up to 5.9. And with the integrations of the differential equations at this distance apart in y, we found
the values of Z. That gave a table in which inverse interpolations gave us new values of y which would lead to
assigned values of Z accurate to five decimals. So that
as it turns out, without too much work, it was possible
to get for any values of nand Z-and we have done this
SCIENTIFIC
COMPUTATION
also for n = 1 and 2-by an inverse interpolation, the
corresponding value of y, to put that in as a starting value,
and to compute the whole field for that particular case.
To compute the field by integration, once the sequence
is set up, is a good deal easier than to get the whole field
at each point by interpolation among the solutions that
have been made, although that would be perfectly possible.
I think actually we shall probably tabulate and publish just
certain selected solutions. However, we hope to be able to
do enough so that anyone could interpolate to find any
solution that they might want.
It should again be realized that the statistical fields only
give a rather rough order of approximation to the properties of the atomic ion. You really want to use these as a
starting field to start after the real approximations in any
particular case in which you are interested. I think it is
clear that it would be too much work to work out complete
sets of Hartree approximations for every possible atom.
I f you have once got a statistical field., that is a good start
toward obtaining the Hartree field. You could probably get
the Hartree field with at most two successive approximations of the Hartree type.
'
Source Exif Data:
File Type : PDF File Type Extension : pdf MIME Type : application/pdf PDF Version : 1.3 Linearized : No XMP Toolkit : Adobe XMP Core 4.2.1-c043 52.372728, 2009/01/18-15:56:37 Create Date : 2010:10:08 19:08:59-08:00 Modify Date : 2010:10:08 21:40:29-07:00 Metadata Date : 2010:10:08 21:40:29-07:00 Producer : Adobe Acrobat 9.34 Paper Capture Plug-in Format : application/pdf Document ID : uuid:feda870e-a966-4263-8f33-73c95c5be26e Instance ID : uuid:afc9c23e-2b3d-4c41-82f5-2e367a3a2b42 Page Layout : SinglePage Page Mode : UseNone Page Count : 127EXIF Metadata provided by EXIF.tools