Statistical Engineering Laboratory

National Bureau of Standards

^Washington 25, D. C.

NATIONAL BUREAU OF STANDARDS REPORT

3371

Motivation for working in numerical analysis

by

John Todd

<NBS>

U. S. DEPARTMENT OF COMMERCE

NATIONAL BUREAU OF STANDARDS
U. S. DEPARTMENT OF COMMERCE

Sinclair Vi eeks. Secretary

NATIONAL Bl'REAU OF STANDARDS

A. \ . Astin, Director

THE NATIONAL BUREAU OF STANDARDS

The scope of activities of the National Bureau of Standards is suggested in the following listing of

the divisions and sections engaged in technical work. In general, each section is engaged in specialized research, development, and engineering in the field indicated by its title. A brief description

of the activ ities, and of the resultant reports and publications, appears on the inside of the back

cover of this report.

Electricity. Resistance and Reactance Measurements. Electrical Instruments. Magnetic

Measurements. Electrochemistry.

Optics and Metrology. Photometry and Colorimetry. Optical Instruments. Photographic

Technology. Length. Engineering Metrology.

Heat and Power. Temperature Measurements. Thermodynamics. Cryogenic Physics. Engines

and Lubrication. Engine Fuels. Cryogenic Engineering.

Atomic and Radiation Physics. Spectroscopy. Radiometry. Mass Spectrometry. Solid

State Physics. Electron Physics. Atomic Physics. Neutron Measurements. Infrared Spectroscopy. Nuclear Phvsics. Radioactivity. X-Rav. Betatron. Nucleonic Instrumentation. Radiological Equipment. Atomic Energy Commission Radiation Instruments Branch.

Chemistry. Organic Coatings. Surface Chemistry. Organic Chemistry. Analytical Chemistry.

Inorganic Chemistry. Electrodeposition. Gas Chemistry. Physical Chemistry. Thermochemistry.

Spectrochemistry. Pure Substances.

Mechanics. Sound. Mechanical Instruments. Fluid Mechanics. Engineering Mechanics. Mass

and Scale. Capacity, Density, and Fluid Meters. Combustion Control.

Organic and Fibrous Materials. Rubber. Textiles. Paper. Leather. Testing and Specifications. Polymer Structure. Organic Plastics. Dental Research.

Metallurgy. Thermal Metallurgy. Chemical Metallurgy. Mechanical Metallurgy. Corrosion.

Mineral Products. Porcelain and Pottery. Glass. Refractories. Enameled Metals. Concreting

Materials. Constitution and Microstructure.

Building Technology. Structural Engineering. Fire Protection. Heating and Air Conditioning. Floor, Roof, and ^ all Coverings. Codes and Specifications.

Applied Mathematics. Numerical Analysis. Computation. Statistical Engineering.

Electronics. Engineering Electronics. Electron Tubes. Electronic Computers. Electronic

Instrumentation. Process Technologv. Radio Propagation. Upper Atmosphere Research. Ionospheric Research. Regular Propagation

Services. Frequency Ltilization Research. Tropospheric Propagation Research. High Frequency

Standards. Microwave Standards.

^Office of Basic Instrumentation ^Office of Weights and Measures.
NATIONAL BUREAU OF STANDARDS REPORT

NBS PROJECT NBS REPORT

1 102-10-1 10*+ June 22, 195*+ 3

Motivation for working in numerical analysis

by

John Todd

<NBS>

U. S. DEPARTMENT OF COMMERCE

NATIONAL BUREAU OF STANDARDS

The publication, re Approved for public release by the

unless permission! Director of the National Institute of

25,d,c. Such pei Standards and Technology (NIST)

cal ly prepared if t on October 9, 2015

• in part, is prohibited

tandards, Washington

;port has been specifieport for its own use,
v : . .. .-i • ;

.

*
Motivation for working in numerical analysis

The profession of numerical analysis is not yet so desirable

that it is taken up by choice; indeed, although it is one of

the oldest professions, it is only now becoming respectable.

Most of those who are now working in this field have been, more

or less, drafted into it, either in World War I or in World

War II, or more recently. The question at issue is why have

we stayed in this field and not returned to our earlier interests.

Our answer is that numerical analysis is an attractive subject

where mathematics of practically all sorts can be used significantly, and from which, on the other hand, many of its branches

can benefit. We call attention here to the applications of functional analysis by the Russian school led by Kantorovitch [l]

(A survey of this has been given recently by Collatz [la]). In

the other direction we would recall the developments in analytic

number theory by Lehmer and Rademacher which followed MacMahon’s

computations of p(n) for Hardy and Ramanujan (see [2]) 0 Before proceeding to a discussion of some individual topics

in numerical analysis, some general remarks are in order.

We have, on various occasions, distinguished between classical

and modern numerical analysis, the latter being material required

in connection with the exploitation of high speed automatic

digital computing machines. It now seems desirable to recognize

ultra-modern numerical analysis, which may be specified as adventures with high speed automatic digital computing machine$[See
- 2-

50, 51 ]o There are, of course, no sharp boundaries between these

parts of the subject, and, as we shall see, there is room for

development in the classical phases as well as in the newer.

In distinction from the deliberate explorations contemplated

in ultra-modern numerical analysis, there is much routine work

in numerical analysis which must necessarily be of an experimental

or empirical nature. It is just not feasible to carry out rigorous error estimates for all problems of significant complications it is necessary to place considerable reliance^ on the one

hand, on the experience of those familiar with similar problems

and, on the other, on the good judgment of the setter of the

problem 0 To justify this remark we consider three examples,,

The solution of systems of twenty or more first order differential

equations is being handled regularly* To see the complication

of theoretical error estimates (in whieli the fact that all

numbers handled are finite (binary) 4feciraals is disregarded) we

refer to Bieberbach [3]» The complication of a stability analysis

in a system of fourteen equations is evident from a study carried

out by F a J 0 Murray [4], Again, the extent of a complete error

estimate for the problem's! matrix inversion is familiar from

- i the work of von Neumann and Goldstine [5 ? 5a] and Turing [6] e Finally, there is the analysis of the triple diagonal method

for determining the characteristic roots of a symmetric matrix

by Givens [7, 7a, 7b] 0 What the numerical analyst has to do is to be aware of the
-

3-

precision of results obtained from, for instance, the conformal

mapping of an ellipse on a circle by a certain process, and from

this to extrapolate to cases of regions of comparable shape 0 On

the one hand, he has to examine general error analyses for their

realism by comparison with cases where the explicit exact results are known* On the other hand, he must devote time to the

construction and study of bad examples so as to counteract any

tendency to too much extrapolation.

The main part of this paper is devoted to a discussion of

some topics in numerical analysis which appear attractive. These

have been chosen among those with which we are familiar, to point

out some of the techniques of the subject and to indicate some

of those who have made distinguished contributions in the field.

In addition, our choice has been controlled by our opinion that

separation between theoretical and practical numerical analysis

is undesirable. As an illustration of the practicality of

some of the techniques used, we refer to computations of the

radiation from a simple source which is reflected from a Lambert

plane recently carried out by P. Henrici [8] where the ideas

of ^2 and ^

5a below were used.

1 . Evaluation of polynomials

What is the best way of computing polynomials, f or instance

f (x + a „ x + a

)

= a

Q

x

n +

n-1 n

o o o
-kfor a series of values of x, not equally spaced? (In the case

where the values of f(x) for a series of equally spaced values

of x are required, building up f(x) from its differences

might be the most convenient e ) The usual answer is to suggest

the recurrence scheme;

f = a

o o

f

r+1 = xfr + a r+1

r = 0 » 1 » . . . >n-1 .

which was known to Newton but is usually ascribed to Horner [9].

In this way we get f(x) by n additions and n multiplications

.

Is this the best possible algorithm? Consider an alternative,

in the case of

f(x) = 1 + 2x + 3x2

„

If we proceed as follows?

2x, x

2

, 3x 2, 1 + 2x + 3x2

we need 3 multiplications and 2 additions compared with the

2 multiplications and 2 additions needed in applying the above

algorithm; thus

3x, 3x+2, x(3x+2), x(3x+2) + 1 9 This problem was formulated as one in abstract algebra by

Ostrowski and he showed [9] that the above algorithm was indeed

the best for polynomials of degree not exceeding It was
- 5-

recently taken up again by Motzkin [10] who showed that the result is not true in practice for larger n Q We give a simple

example in the case n = 6* Consider the evaluation of

p = x + Ax^ + Bx + CxJ + Dx + Ex + F c

Introduce the following polynomials

P.j = x + ax

P

2 = ^ P1 + x + ^ P1 +

P

3 = (P2 + d) (P

1 + e)

and determine a,b,c,d,e and f by identifying P and P^+f* This

can be done by the solution of linear equations and a single

quadratic* This evaluation is done once for all and then P

can be evaluated at the expense of three multiplications only,

with a significant economy over the other process if we have

to evaluate P for a sufficiently large number of values of x*

Whether it is possible to improve on the Newton-Horner scheme

by a purely rational process does not seem to have been decided*

The details of the evaluations are as follows* The result

of equating coefficients in P and P^+f iss

(1.1) 3a + 1 = A

(1 .2) 3a^ +2a+b+c+e=B

(1 .3) a^ + a^ + 2ab + 2ac + 2ae + c + e = C

0 0 0 (1» a b+a c + a e + ac + ae + be + be + ce + d — D

(1*5) abc + abe + ace + ad + ce = E

( 1 *6) ebc + de + f = F
-6-

From (1,1) we find a, Hence we can rewrite (1,2) and (1,3)

in the form

(1,2*) b + c + e = B

1

( 1 0 3

T ) 2a ( b + c + e ) + c + e = C

1

(We use primed capitals to indicate new known constants) 0 Using

( 1 .2') in (1 .3* ) we get

(1 .3") c + e = C"

which with (1,2*) gives us b explicitly. Using a,b, c+e, we can

write ( 1 ,4) as

(1 >") d + ce = D

! Using a,b, c+e, d+ce in (1,5) we find

(1 .50 ce = E'

which gives d from (1,4*), From (1,3 M

) and (1,5 ? ) we can find[ - - ii C>e by solving a quadratic equation and then from (1,6) we can

find f,

2, Increasing the speed of convergence of sequences

The construction of processes which increase the speed of

convergence of sequences and series has been a favorite topic

p for many numerical analysts. For instance, there is the h

extrapolation process of L, F, Richardson [11], the convergingfactor method of J, R, Airey [12], [12a], the Euler summation

process [13] and a whole subject associated with the name of
- 7

-

2 Chebyshev [1*+]. We shall discuss the 6-process which

popularized in numerical analysis by A. C. Aitken [l5?e

dates back at least to Kummer [16],

If

x

n

x

and

( 2 . 1 )

then

( 2 . 2)

x - x - A /\ n

, n • /v* m< i

x n+2

- x • x

n+l ~ x

n+1 - x o x - x

n

From (2.2) we find

x

T

X n+2

( x n+2

~ x n+1>

x n+2

" 2x n+1

+ x

n

This suggests that the sequence [ x . ~ j defined by

x n+2

= x n+2

(x n+2

- x

n+

1

2

n_ Q n

, x , 0 - 2x . 1 + x

’ *n+2 n+1 n

converges more rapidly to x than the original sequence

is indeed the case for if

x

n - X = A 7(

n

+ o( A

n

) , IA| <1

then it follows that

x

n " xn = 0( *

n)

*

has been

.g]5 it

This
- 8-

Several remarks are in order. First, this process can be

iterated to remove successively components in the remainder of

the form

A TV > B ^ , C 5000

wherel> 1^1 > lyA| > > • • . . The c ases when there are

equalities such as |AI =|^ can be handled by simple modifications.

Second, it is important to note that this process can make things

worse if the convergence is not geometric as required by (2.1).

Here is a simple example involving two of the standard iterative

processes for determining the reciprocal of a number N. Consider

the sequences

y n+1 = (1 - N) y n +1j zn+1 = z

n (2 - Nzn

)„

In the case N = £ with yQ = 1 , z

Q = 1 , we obtain the following

table:

1 1.5

1.75

1.875

.5

.25

.125

-.125

-.125 2

2

1 1.5 -.125

.375

1.875 -.2578 3.0000

.1172

1 .9922 2. oh-55

The sequence {yn j[ appears to converge more rapidly than |y )• while the sequence

|

z

nJ" appears to converge less rapidly than

^ z

Q

These results. can be easily established. First of all, each
-9-

sequence converges to N MBt if 0 < N < 1 for

y n - N

1 = (1 - N

1 )N

(yQ - N |)

and

z

n

Thus [yn ^ satisfies the condition (2 e l) while ^z

n

^j does not,

converging too rapidly. In the present case we have y n = N

On the other hand it can be shown that

z n fpr 00

z n

Note, however, that to justify the application of this process

it is sufficient to show the existence of an expansion of the

form (2.1) with /A( < 1,

Extensive use of this process was made in experiments in

conformal mapping by G. Blanch and L. K 0 Jackson [17] and by

John Todd and S, E. Warshawski [l8] e For instance, in the

latter the mapping of an ellipse (of axis ratio 5s 1) on a circle, it was found that while about 50 iterations, each

requiring about 30 minutes of computing on SEAC, were required

to secure directly about nine correct decimals in the value of

the boundary function. It was however, possible to obtain

the same accuracy by a double use of the Aitken process on the

first 14 iterants — the extra time required for this being

negligible.
- 10 -

3, Modified Differences

We shall show here how the use of quadratic interpolation

enables the table-maker to cut down on the size of a table at

the expense of some work by the table-user,, We shall then show

how a further saving in space can be accomplished, at no further

expense to the user but at some to the table-maker, by the use

of modified differences. For simplicity and definiteness we

consider the construction of a table of sin x, to four places of

decimals, for x in the range (0, .

a) ly ing a.r_ J.P1S

r

pp laiipji 0 The error involved in linear interpolation, i.e., the

assumption that

f(a+ph) = f(a) +

p

^f ( a+h)-f ( a)j[ , 0 < p < 1 ^

can be estimated as

h

2

| max jf "(x)

j .

If this is to be less than ^xIO an appropriate choice for h

is .02. This requires a table of some 80 entries, part of which

is shown below.

X sin_>

0.00 0.0000

0.02 0.0200

0.04 0.0400 Interpolation, say for x=1.231+, in this table is

0.06 0.0600 carried out as follows:

0.08 0.0799

O sin 1.234 = .9391 +

^ (.9458 -.9391) = .9438

0

CvJ O o!9320

1.22 0.9391

1.24 0.9458

1.26 0.9521

1.28 0.9580

O
- 11 -

b) Quadratic interpolation

We now consider using the Everett interpolation formula

(3.1) f^qfQ+pf^EQ^^Q+E^O2

.

where q = 1-p,E

Q

2 = q(q2-l)/6, E^ 2 = p (

p

2-1 )/6 , . . .

If we retain the first four terms the truncation error can be

estimated as

(V >)

max t

k M

4

< h x .024 x 1

1 — < For this to be less than 2 x 10 we can conveniently take h= 0 2.

The corresponding complete table is given below 0 X sili_x &

•

0.0 0.0000 0 For interpolation we now have either to

0.2

0.4

0.1987

0.3894

- 80

-155 compute the Everett coefficients or to

0.6 0.5646 -224 obtain them from a table; we find

0.8 0.7174 -287

0.8415

p = .17, E

n

^ = -.0430, E^ = -.0275.

1.0 -336 u

1.2

1>

0.9320

0.9854

-371

-392 We then have

1.6 0.9996 -400

1.8 0.9738 -387

sin 1.234 = .9320 + 200 (.9854 -. 9320 ) + c .0430) (.0371)

+ (.0275) (.0392) = .9438
- 12-

c) Comrie 1 s throwback

This device, introduced by Comrie [19] depends essentially on the fact that the ratio

£

1+ V( D ) - - (p±Q(p-3.) K^p; 2 20

E

0

is approximately constant for 0 < p < 1 . Various ways of choosing

a mean value for this have been discussed [20, ]. The

preferred value is k = — . 1 3393 • With this value of k we rewrite

the first four terms of (3,1) as:

(3.2) f

D 7

qf

o

+ P£

1 +

E

0

2

(6

2

f

0 + ko

1

'' f

o

) + E

1 2

(6

2

i

1 + ko^)

Therefore if we define

6

2

f = 6

2

f + k 6^f

m

and use these modified second differences in exactly the same way

as we used the ordinary second differences in the preceding subsection, we can obtain the desired accuracy of interpolation

with a much larger interval. Iq fact the error in (3.2) is

made up of the truncation now bounded by

(3.3) h

6

|(

P

6

3

)| max |f

(6) (x)| < h

6

x .00^9
- 13-

together with the error caused by the modification. It can be

shown that the latter is less than half a unit in the last place

if the fourth differences are less than 1000 units (and the

fifth less than 70 units).

The condition (3.3) gives h < .46, which suggests that h = .5

might be acceptable. For this value of h, a bound for the fourth

difference is

(.5)^ x |f^(x)

| < .0625

which is acceptable, although the bound for the fifth is not.

Nevertheless we shall use h = .5 without carrying out a more precise estimate.

The complete table is given below:

X

0.0

sin x

0.0000

o?m

0

- 0.5 0.4794 -1225 For interpolation we first find

1 .0

1.5

0.8415

0.9975

-2154

-2552 the Everett coefficients

2.0 0.9093 -2326

P = .468: E

0

2 = .0636, E

1 2 = -.0609.

We then have

Sin 1.234 = .8415 + (.9975 - .8415) + (.0636X .2154),

+ (.0609x .2552) = .9145 + .0137 + .0155 = .9437

The discrepancy between the results can be explained either by
the marginal choice of h, or by rounding errors.

More elaborate methods of throw-backs e o g 0 of the sixth difference as well as the fourth into the second have also been

given by Comrie [21], Discussions of the minor disadvantages

of the modified differences are also available [19]*

4 , Characteristic roots of finite matrices

Considerable effort has been expended in problems of

numerical analysis involving matrices [see e.g., 72, 73? and

22, 23] o The two main problems are the inversion of matrices

and the determination of their characteristic values. In both

problems the practical determination of bounds for characteristic

roots is important. In this connection we call attention, to the following lemma of Gerschgorji [2*+] which has many applications.

All the characteristic roots of A = (a. .) lie in the union of

the circular regions

This is proved by use of the fact that a determinant with dominant

main diagonal does not vanish. This last result has been generalized

by many writers: for an account of some of this work (up to 19^7)

see 0. Taussky [25].

One of the preferred and practical methods of getting all the

characteristic roots of a symmetric matrix depends on the reduction of the matrix to pure diagonal forms (Goldstine, von Neumann

and Murray [ 26 ], Gregory [27]) by superposing orthogonal trans-
- 15-

formations involving two variables at a time. Theoretically we

obtain

TAT 1 = diag ( ^ 1

, 7\ 2

, . . . 5 7\ n

)

and then ^ ^ 2

’ 00# ’^n are t *ie exact characteristic values.

In practice we find

TAT’ = ( £ )

where the £ . . are small for i^j. We then ask how near are the

•I ^

to the If we disregard the question of the transformation

not being truly orthogonal — and therefore the characteristic

roots of (/ . .) not being identical with those of (a. .) - the

ij

answer comes at once from the lemma. If the £ . . , i J- j are

^ J sufficiently small then

a, -f 4,1 < 11 i/j ' 1 -> i = 1 , 2 , o o a J n .

Allowance can easily be made in this inequality for round-off

error in the product TAT’.

5 • Quadrature* i ntegral equations

a) We begin with an example to show that there is still scope

for new ideas in classical numerical analysis. A typical quadrature formula is
-16-

b ^

J

f(x) dx = ^

p i

f (x^

)

a and the error

is estimated as a multiple of a (high) derivative f

n

(x) of f(x)

at a point in (a,b). In many cases it is far from convenient to

obtain bounds on f

n (x), or to estimate these by computing the

corresponding differences manually. Recently P, Davis and

P, Rabinowitz [28], [28a] reconsidered this problem in the case

when f(x) is analytic in a region including the segment (a,b) 0 W. F. Eberlein [29] has also contributed to this problem. The

case of ellipses^ with foci at the end-points, which we normalize to (1,0) (-1 ,0) ^can be handled elegantly, in terms of the

Chebyshev polynomials,

(1 - z ; 2 sin( (n+1 ) arccos z)

which are orthogonal over the area of such an ellipse. It can

be shown that

E < II til

where is a constant depending only on the ellipse £ and where

II 1 II = ff/t(z)l

2 dx dy
- 17-

(Note that \\ f \\ increases as ^

expands; however, ^then decreases

and there is a problem of optimal choice of £ ) . The (T^ can be

tabulated once for all, and 1! ffl can be estimated in terms e.g, , of max | f| . As an example consider the evaluation of

J

r (x)dx

3

using a 7~point Gaussian rule. To evaluate and bound the fourteenth derivative of P(z) seems rather out of the question. Simple

estimates can be used in the method just described to find

E < 2.04 x 10 12

. A comparison of this estimate with the usual one [30, 31]*

f (2n) (j) (n!) 4

^ (2nJ) 3 (26+1)3 S<J s< ‘S

where the derivative is now estimated by the use of Cauchy's

f ormula^shows that the new one is somewhat better,

b) Among the basic problems in the numerical analysis of integral

equations is the relation between the eigenvalues of a (symmetric)

kernel and those of an approximating matrix. A satisfactory account

of this was given recently by Wielandt [32] in support of some

experiments on conformal mapping [18] which were being carried out
- 18-

on SEAC, the National Bureau of Standards Eastern Automatic

Computer, The continuous problem is the solution of

J

K(x,j ) y(|) dj = ky(x)

0

We make this discrete by introducing a quadrature

(5.D

J

1 f(j) d$ = I P, f(^)

0

and are therefore led to consider the matrix problem

What are the relations between the finite number of k and the

infinity of the k? We quote a typical result. If we take for

(5.1) the trapezoidal quadrature

S

d

J

= ifK + f

i + —+f „-i

+ K\ » f i =

0

J

then, provided K satisfies

|K(x,y) - KW,6)| < L (|x^l|+ jy -/S J )

where ^

run through the points (Vn, -Vn) and where

& A A |x -ol j< in , jy ~/t>l< i n > we have

k - k < ^

n-1
- 19-

where the constant C = ^

+ V 1/12 is best possible

6. Game theory and related developments

There are here problems in which the intuition of

a geometer can play an essential part; for instance, the theory of

polyhedra and convex bodies and fixed point theory are all

highly relevant. The foundations of a theory of games were laid

down by von Neumann [ 33 ]? [ 39 ? 7^? 75 ] beginning in 1928. The

theory of two-person zero sum games is well developed but the prac

tical problem of finding the value of such a game and the optimal

strategies is difficult and solutions available so far are not

entirely satisfactory. (Among related and essentially equivalent

problems are the solution of systems of linear inequalities and

the solution of linear programs in the sense of G. B.Dantzig.

)

Among the methods of attacking these problems are the Simplex

Method [ 3*+]? the Relaxation Method [ 35 ]? the Bouble Description

Method [3^]. We shall, however, discuss a very simple example

by a natural approach due to G. W. Brown [ 37 ]? the validity of

which was established by J. Robinson [33]# A related continuous

solution of this discrete problem has been given by Brown and

von Neumann [jjfo]; we take up this idea of continuous approach to

discrete problems again

Consider the following game played between two players R and C

each of whom has two strategies, which we may interpret as the
- 20-

choice of a row or a column in the pay-off matrix P = (p — )t

' -

<i !> •

If R chooses the ith row and C chooses the jth column then R gets

p.. from C. This is manifestly an unfair game, and R should pay

to play it.

The value of this game is 2„5 and the optimal strategies are

the following: R should choose 1 and 2 each with probability 25

•1 C should choose 1 with probability ^

and 2 with probability 3/^.

The significance of these statements are the following: if R plays

in this way his expected gain is not less than 2.5 while if C

plays this way his expected loss if not greater than 2.5*

To prove this statement is simple. Let R play 1 with probability

* - r ^ 0 and 2 with probability 1-r ^ 0; let C play 1 with probability

c >0 and 2 with probability 1-c ^ 0. Then the expectation of

R is:

E = 1 x rc + ^(l-rjc + 3r(1-c) + 2 (1-r)(1-c).

We have:

E = -b ( r-*jr) (c- xjO + 5/2.

This shows that when r = 2 then E = 5/2 for any c, and that if

-1 r / jr then c can be chosen to make E < 5/2; similarly, if c = ^

-i then E = 5/2 and if c / ^

then r can be chosen to make E > 5/2.

(ef, McKinsey [39]).
- 21 -

How can we arrive at these results, or approximations to them?

Suppose we consider a sequence of plays of the game in which each

player is trying his best in the following sense. After n plays

suppose C has chosen 1 in c

; x n plays and 2 in the remaining

( n)

c ^ * n plays. If C continues this pattern, the expectation of

R in the next play is

e^ = c^ < 3C

2

if he chooses 1 , 2c^ if he chooses 2.

R therefore chooses 1 or 2 according to whether e^ ^©2 or e

i < e 2°

Similarly, if R has chosen 1 in r^

^ n ^x n plays and 2 in the re-

/ n ) maining * n plays, then C chooses 1 or 2 according to the

expected size of his loss which is

f

.j = r.j + if he chooses 1 , *2 = 3 r

i + ^ r 2

if he chooses 2,

Specifically C chooses 1 if f^ N< f

2 and 2 if f^ > f

2

» This sequence of plays is determined after we make an (arbitrary)

assignment of the initial play e.g,, that each player chooses 1 e The resultant sequence of plays is:

(1,1), (2,1), (2,2), (2,2), (2,2), (2,2), (2,2),

(2,2), (1,2), (1,2), (1,2), (1,2), (1,2), (1,D, (1,D,

(2,1), •N CM

•N CM (2,2), (2,2), (2,2), (1,2), (1,2), (1,2),

(1,1), (i,D, (2,1), (1,1), (2,1), (2,2), (2,2), (2,2), (2,2),

(2,2), (2,2), (2,2), CM CM (2,2), (2,2), (2,2), (2,2), (1,2),

(1,2), (1,2), (1,2), (1 ,2), (1,2), (1,2), (1,2), (1 ,2), (1,2),
- 22-

It has been shown C

3<C3 5 that as n—> 00 , the

sequences^

j , (

r

i converge to the optimal strategy i.e • ?

c r- r C

(n L | ^ IJ? r

-J

1 ( n ) ^ , and that the average pay-off p converges to the value of the game. In our case

(50) _

1 = >8, (50) _1 = .2, (50) _ 2,b,

The structure of the sequence above, consisting of blocks of

identical elements, is typical; this can obviously be used to

speed the computations. For some practical experiments in this

field see [ 40 ],

We shall now discuss an application of the theory of games

to the so-called assignment problem. This problem is to assign

n square pegs to n round holes in such a way as to maximize the

total goodness of fit. In other words, (a. .) being given, we

^ J have to choose a permutation of (1,2,...,n) so as to maximize

(6.1 ) *> iP

i

This is trivial theoretically; we have only to find the largest

of the nl sums of the form (6.1). In practice, however, this may

be out of the question and so we may have to settle for some approximation to the maximum. One way of doing this (suggested

by von Neumann (^ 0),

is to set up an equivalent game theory

problem it turns out to be a sort of hide and seek and

solve this approximately by the method just discussed. The first
-23-

player chooses a pair of indices (i,j) (1 ^ i ^ n, 1 j n);

he has n strategies. The second then elects first to guess the

first or second of these two indices, and then guesses it by

choosing k(1(k<n); he has 2n strategies. In the first case if

k=i, and in the second case if k=j, the first player pays the

_ i second (a. .) , otherwise there is no pay-off.

J Assignment problems for n=12 have been handled by this method.

An up-to-date account of this problem and its generalizations

has been given by Motzkin [42], Among these are the transportation problem, the caterer problem, the problem of contract awards,

the traveling salesman problem [43, 44], Solutions to problems

of this type are now obtained on a routine basis, on high speed

computers, as an aid to management decisions in industrial and

military situations [43], Among other problems of this general

character, which are in the research stage, are those concerned

with organization theory which have been studied by Tompkins

and Marshak [46].

7. Monte Carlo

This is a subject with large areas unsoiled by theorems, as

can be seen by reference to the reports on various symposia held

on the subject [67, 68], For instance, during the last four

years we have been generating millions of pseudo-random numbers

on SEAC using such relations as
-24-

r n

2

-42

or

r n 2

-44

The r behave as if they came from a uniform distribution in the

interval 0.1. The results we obtained were satisfactory in all

cases where we had independent checks. We have, however, no

theorems at all about the "randomness’1 of these sequences or

about the distributions of blocks of the size used in our calculations.

We mention here also the quasi-Monte Carlo processes studied

by Peck and Richtrayer [47,48]. Here high-power algebraic number

theory is used to evaluate the error committed by replacing integrals by sums of the integrands at points determined by certain

algebraic numbers.

We mention here a few areas, with which we are familiar, and

which we have found interesting. This personal selection omits

reference to many areas in which there have been important advances (e.g. meteorology) and to areas which have been discussed

elsewhere in this volume, or its predecessor. .

8 . Recent activity in numerical analysis

a) Ultramodern numerical analysis

One class of experiments may be described as follows
-25-

It has been usual in discussing properties of matter to regard

the medium as continuous, set up differential equations, look

at them for a while, give up and replace them by difference equations. These difference equations were then solved and no attention was paid for their physical significance, if any. An alternative approach is to handle the problem discretely from the

beginning, lumping the "molecules” together in groups as small

as the computing equipment can handle.

Among those who have handled problems in this general way are

Seeger, von Neumann and Polachek (cf. Seeger [49]), who were

concerned with shock wave phenomena. Recently Ulam [50] has

studied the mixing of fluids and the motions of star clusters

in this way. Metropolis and Fermi [5l] have investigated the

equations of state of individually interacting particles forming

an idealized liquid,

b) Biological applications

There has been pioneering work by Turing [52] on the

problem of morphogenesis. Turing constructs a mathematical model

of a growing embryo and shows how well-known physical laws are

sufficient to account for many of the facts about the development of its anatomical structure.

Another application has been the study of the reaction of nervefibers to electrical stimuli. This phenomenon is governed by a system of four non-linear ordinary differential equations
-26-

(Hodgkin-Huxley) . The system has been studied on SEAC (by

K. S. Cole, H. A. Antosiewicz, and P, Rabinowitz [ 53 ]* and in

particular the threshold value of the input current has been

determined. The agreement with the results of millions of experiments indicates the reliability of the model and encourages

further investigation,

c) Combinatorial analysis

recent reports on this topic by S. S e Cairns [5*+] and C. B.

Tompkins [ 55 ]* The numerical analyst, however, soon finds himself out of his depth if he uses straight-forward approaches.

One new idea which was tried is the followings that of a continuous approach to discrete problems, in particular to the search

for perfect difference sets. , A perfect difference set is a set of n+1 integers whose n(n+1) differences

p take on all non-zero values mod n + n+1. For example, the

differences of 1,2,4, are + 1 , + 2, + 3 i,e 0 all non-zero values

mod ^ , and so 1,2,4 form a perfect difference set mod . 2 A perfect difference set S can be specified by N=n+n+1 constants x^ where x = 1 if rf S, x =0 otherwise. In this case

t r

1 r we have

This is an obvious source of problems. There have been

(The subscript v~+s is to be understood mod N); hence
-27-

(8.1) Z. y = (n+1) 2

. r ^

It follows,7 therefore, that such a set x minimizes 7 r

J = (n+1) X. y g

2 - n yQ

2

for, in view of (8.1), J differs by a constant from

(y n - (n+1)) 2

+ (n+1) 21 (y - I)

2

, ^

. S>0

This suggests an attempt to obtain a set of x

p

by minimizing J,

now regarded as a function of the N continuous real variables x,

subject to (8.1), (and perhaps to other relations such as

o <: x r ^ i).

Such an attempt was made on SWAC, by a steepest descent process 0 Although admissible values of y were obtained rapidly, the corresponding values of x were not integers,

d) Number theory and algebra

These subjects are natural sources of problems and there

has been many applications of high speed computers, particularly

in elementary, algebraic and analytic number theory and some in

algebra proper [55a, 56].

Recent work on SWAC (The National Bureau of Standards Western

Automatic Computer), mainly on elementary number theory by D 0 H e and E 0 Lehmer and their collaborators, has been discussed by

E. Lehmer [57].
-28-

Among other work has been a study of the divisibility of

( ( p— 1)1 + 1 )/p by p. This was known to be the case for p=5?

* 4 p=13; K. Goldberg [58] found that it was also the case for p=5&3

and for no other p < 10,000 o Problems in algebraic number theory are more complicated to handle.

A survey of computational problems in this field has been given by

0. Taussky [59]. Since then there has been work by H. Cohn and

S 0 Gorn [ 60 ] on units in cubic fields.

There have been various attempts to study the zeros of the Riemannzeta function; among those is that of Turing [6l] 0

e . Topology

It is clear that approximate computations of quantities known

to be integers serve to define them if the absolute value of the

error is known to be less than This idea is used in the work on

p(n) mentioned earlier. Ulam [ 50 ] has suggested that further applications can be made in essentially topological problems, e 0 g o the structure of the lines of force caused by current in two infinite

straight wires which are skew.

A simple application of this is to the location of the zeros

of a polynomial P(x). We use the fact that
-29-

where n is the number of zeros inside the simple closed rectifiable curve C. It is possible to choose C to be a square so

large as to contain all the zeros of P(z) and then, by process

of quadrisection, to locate the zeros approximately . The quadrature must be accomplished with absolute error less than

if this proves difficult due to vanishing or near vanishing of

P(z) on the boundary, then we know that we are in the neighborhood of a zero and can act accordingly,, A constructive proof

of the fundamental theorem of algebra along these lines had

been given by Rosenbloom [ 62] 0

9, Theory of machines or automata

Among those who have contributed to basic research in this

direction have been Turing [ 63 ], -Shannon [ 6*+] and von Neumann

[65]o There has been some efforts of a character of supporting

research; the use of machines to design circuits for better

machines, the design of self-correcting codes, and improvements

in the use of machine s^ e 0 g„ more automatic coding, Much of

this belongs more to the domain of logicians than to that of the

numerical analyst.
-30-

BIBLIOGRAPHY

References

1* L. Kantorovitch, Functional analysis and applied mathematics,

Uspekhi Mat* Nauk, 3? 19^8, pp„ 89~l85 0 la.L. Collatz, Einige Anwendungen f unktionalanaly tischer Methoden

in der praktischen Analysis, IZeit* ang* Math* u*

Physik, b, 19 53? PPo 327-357*

2. G, H. Hardy, Ramanujan (Cambridge, England 19^0)*

3. L 0 Bieberbach, On the remainder of the Runge-Kutta formula i^

the theory of ordinary differential equations ,2e:J.t e ang. Math* u. Physik 2, 1951? pp. 233-2^8*

b. F* J e Murray, Planning and error consideration for the numerical solution of a system of differential equations on a sequenced calculator,, Math 0 Tables

and Other Aids to Computation, *+, 1950, pp* 1 33“ 1 «,

5o John von Neumann and H 0 H 0 Goldstine, Numerical inverting of

matrices of high order* Bull* Amer 0 Math c Soc* 53?

19^7, PPo 1921-1097.

5a*H» H. Goldstine and J* von Neumann? Numerical inverting of

matrices of high order* Proc c Amer* Math* Soc 0 , 2? 1951, PPo 1 88-202*

6* A. M* Turing, Rounding-off errors in matrix processes* Quarterly

J* Mechanics and Applied Math*; 1, 19^+8, pp*207~308*

7m Wallace Givens, Numerical computation of the characteristic

values of a real symmetric matrix, Oak Ridge National Laboratory, Report 1 57^? 195*+

7a* Wallace Givens, A method of computing eigenvalues and eigenvectors suggested by classical results on symmetric

matrices , pp* 1 1 7“1 22 in ( 72 )*

7b. Wallace Givens, A method of error analysis, pp. in (70)*

8 . P. Henrici, Computation of the intensity of radiation reflected

from an infinite Lambert plane, to appear in

Journal Washington Academy of Sciences*
-31-

9.

A. M. Ostrowski, On two problems in abstract algebra connected

with Horner’s rule, to appear in R. von Mises

Memorial Volume (New York, 195*0 . 10. T. S. Motzkin, unpublished.

11. L. F. Richardson, The deferred approach to the limit, Phil.

Trans. A, 226, 1927, pp. 299-350.

12. J. R, Airey,The converging factor in asymptotic series and the

calculation of Bessel, Laguerre and other functions, Phil. Mag. ser. 7, 25, 1937? pp 0 521-552.

12a.J,C.P, Miller, A method for the determination of converging

factors, applied to the asymptotic expansions for

the parabolic cylinder functions. Proc. Cambridge

Phil. Soc. 48, 1952, pp. 243-254. ^ ' * — . -

; - <1 13. 1. Barkley Rosser, Transformation to speed the convergence

for series, J. Research National Bureau of Standards, 56, 1951? pp 0 56-65.

15.

C. Lanczos, Introduction to Tables of Chebyshev Polynomials.

National Bureau of Standards Applied Mathematics

Series, No. 9 (U 0 S e Government Printing Office,

Washington, D 0 C., 1952).

15. A. C. Aitken, Studies in practical mathematics 6. On the

factorization of polynomials by iterative methods,

Proc. Royal Soc. Edinburgh 63A, 1951? pp 0 175-191.

16. E. E, Kummer, J. reine u. ang. Math. 16, 1837? pp 0 206-215.

17. G. Blanch and L. K. Jackson, Computation of harmonic measure

by L. Ahlf or s’ method , pp. in(7l).

1 8. J. Todd and S. E. Warshawski, On the solution of the Lichtens tein-Ger schgorin integral equations in conformal

mapping? II. Computational Experiments, pp.

in (71).

19. L. J. Comrie, in British Association Mathematical Tables vol.

1(ed. 1, 1931? ed. 3) 1951)) Cambridge, England,

pp. xiii-xiv.

20. M, Abramowitz, in National Bureau Standards, Tables of Bessel

functions of fractional order, vol. 1 (Columbia

Univ. Press, NY) 1957? pp. xxxiii-xxxvi

.
-32-

21. Interpolation and allied tables, reprinted from the Nautical

Almanac (H.M. Stationery Office 1937)? pp. 928-929.

22. G. E. Forsythe, Solving linear equations can be interesting,

Bull. American Math. Soc. 59? 1953? PPo 299“329.

23. M. R. Hestenes, Computational methods, pp. in this

volume [ 66 ],

2^. S.Ger#ehgotio,Uber die Abgrenzung der Eigenwerte einer Matrix.

Izv. Akad. Nauk, SSSR 7? 19*+1? PP. 7^9-75^.

25. 0. Taussky, A recurring theorem on- dete rminants , American

Math. Monthly 56, 1 9j9 , pp. 672-676 . 26. H. H. Goldstine, J. von Neumann and F. J. Murray, unpublished.

27. R. T. Gregory, Computing eigenvalues and eigenvectors of a symmetric matrix on the ILLIAC, Mathematical, Tables

and Othar^Aidi to Computation 7?1953?pp.21 5-220.

28. F. Davis and P. Rabinowitz, On the estimation of quadrature

errors for analytic functions, Math. Tables and

Other Aids to Computation 8 , 195*+? pp.

28*. P. i>avis, Errors of numerical approximation for analytic functions, Journal Rational Mech. and Analysis 2,

1953, pp. 303-313.

29. W. F. Eberlein, , Preliminary report

in Proc. American Math. Soc. 195^? pp.

30. G. Szego, Orthogonal polynomials, American Math. Soc. Colloquium

Publication, 23? 1939.

31. A. N. Lowan, Norman Davids and Arthur Levenson. Table of the

zeros of Legendre polynomials etc., Bull. American

Math. Soc . 48, 1942, pp. 739-763.

32. H. W. Wielandt, Error bounds for eigenvalues of symmetric

matrices, pp. in ( 70 ).

33. J. von Neumann and 0. Morgenstern, Theory of games and economic

behavior, ed. 2(Princeton Univ. Press., Princeton,

N. J.), 1 9^7

o
-33-

34. G. B. Dantzig? Maximization of a linear function of variables subject to linear inequalities, pp 0 339~ 3^7

in Activity analysis of prediction and allocation,

ed. T. G 0 Koopmans (New York, 1951}.

35. S. Agmon, The relaxation method for linear inequalities, Canadian J. Math, , 1954, pp.

36. T # S. Motzkin, H 8 Raiffa, G. L 0 Thompson and R 0 M. Thrall, The

double description method, pp c 5~73 in (75).

37. G. W. Brown, Iterative solution of games by fictitious play in

Activity analysis of prediction and allocation,

ed. T 0 G. Koopmans, New York, 1951* pp 0 374-376*

37$. G. W. Brown and J 0 von Neumann, Solutions of games by differential equations, pp 0 73~79 in (74) 0 33. J. Robinson, An iterative method of solving a game, Annals of

Maths a 54, 1951? PPo 296-

39. J.C # C. McKinsey, Introduction to the theory of games, (New York,

McGraw-Hill, 1 9 52J 0 40. A. J, Hoffman, Me Mannos, 1

0

o Sokolowsky and N 0 Wiegmann, Computational experience in solving linear programs,

Jo Society for Ind 0 and Applied Maths, 1, 1953?

pp. 17-3^.

John von Neumann, A certain zero-sum two-person game equivalent

to the assignment problem 0 pp 0 5“ 12 in (75)

•

42. T. S. Motzkin, The assignment problem, pp 0 in (70) 0 43. H. Kuhn, The traveling salesman problem, pp e in (70) 0 44. G. Dantzig, R. Fulkerson and S c Johnson, Solutions of a large

scale traveling salesman problem,

45. L. Gainen, 0. Honig and E, D. Stanley, Linear programming in

bid-evaluation, Logistics Quarterly, 1, 1954,

PPo

46. J, Marshak and C c B, Tompkins, unpublished 0 by. L. G. Peck, On Uniform distribution of algebraic numbers, Proc 0 American Math, Soc 0 4, 1953? pp* 440-443,
-y+-

48* R 0 D„ Hiichtmyer, The evaluation of definite integrals and a quasi-Monte Carlo method based on the properties

of algebraic numbers 0 49. ft, Jo Seeger, On computational .techniques for certain problems

in fluid dynamics^in Proceedings of a Symposium

on Large Scale Digital Calculating Machinery,

Computation Laboratory, Harvard University Annals

vol. 16 , 1948, pp. 157- 168 . 50. S. M. Ulam, Heuristic studies in problems of mathematical physics on high speed computing machines*

51 * N* Metropolis, Some topics in experimental mathematics, Proceedings of a symposium on large scale digital

computing machines, Argonne National Laboratory,

August 3-5? 1953o Report ANL-5181, pp 0 265-270*

52* A, M* Turing, The chemical basis of morphogenesis, Phil* Trans

Royal Soc * B237, 1952, PPo 37"72*

53. H* A e Antosiewicz, K 0 S* Cole and P* Rabinowitz, unpublished*

9+. s* s* Cairns, American Math* Monthly*

55. C* B. Tompkins, Application of automatic digital computers to

problems with discrete variables, pp,

in ( 70 )*

55a. L* J* Paige and C, B* Tompkins, Systematic generation of permu

tations on an automatic computer and an application to a problem concerning finite groups*

56* K* Goldberg, log e x e^ in a free associative ring, Abstract in

Bull* American Math* Soc* , 1951+o 57. E* Lehmer, Number theory on the SWAG, pp* in (70)*

58* K* Goldberg, A table of Wilson quotients and the third Wilson

prime* Journal London Math* Soc* 28 , 1953?

pp* 252-256*

59. 0* Taussky, Some computational problems in algebraic number

theory, pp* in ( 70 )*

60 . H* Cohn and S* Gorn, unpublished*
-35-

6l e A. M. Turing, Some calculations of the Riemann Zeta-f unctions

Proc 0 London Math 0 Soc 0 {

3

?

1

953 ? PPo 99—1

1

7 0 62. P. C. Rosenbloom, An elementary constructive proof of the fundamental theorem of algebra 0 American Math 0 Monthly,

52,1 9^5r pp. 562-570.

63. A. M. Turing, On computable numbers with an application to the

Ent scheidungsproblem 0 Proc Q London Math. Soc., {2^,

42, 1936 , 230-265 and 43,1937, 544-546.

64. C. Shannon, Computers and automata, Proc. Inst. Radio Eng.,

4-1, 1953, PPo 1234-1241.

65 . von Neumann, The general and logical theory of automata, pp.

1-41 in Cerebral Mechanisms in Behavior, ed. L a A.

Jeff re ss , (New York, Wiley, 1951 9 0 66. eds. , Transactions of the

Symposium on Applied Mathematics, University of Chicago, 1954 (Inter science , New York-London).

6 7. A. S. Householder, G. E. Forsythe, H. H 9 Germond, eds. Monte

Carlo Method, National Bureau of Standards, Applied

Mathematics Series, 12 (U. S. Government Printing

Office) 1 951

o

68. Proceedings of a Symposium on Monte Carlo Method, Gainesville,

Fla., March 1954, to appear.

69 . G. Birkhoff, K. 0. Friedrichs, T. E. Sterne, eds. Transactions

of the Symposium on Fluid Mechanics and Computing,

New York University, 1953 (Inter science , New YorkLondon) . 70. J. H. Curtiss, Proceedings of the Sixth Symposium on Applied

Mathematics* Numerical Analysis, Santa Monica, Calif.

1953 (McGraw-Hill, New York, N. Y.)

71 • J, Todd, ed.. Experiments in the computation of conformal maps,

National Bureau of Standards, Applied Math. Series

42, (U. S. Government Printing Office 1954).

72. L.J. Paige and 0. Taussky, eds. Simultaneous linear equations

and the determination of eigenvalues. National

Bureau of Standards, Applied Mathematic Series 29,

(U, S. Government Printing Of f ice , 1 953)

»
-36-

73. 0. Taussky, ed 0 , Contributions to the solution of linear systems

and the determination of eigenvalues, National Bureau

of Standards, Applied Mathematics Series, 39? (U 0 S.

Government Printing Office, 195*0

7^. H. W. Kuhn and A, W. Tucker, eds. Contributions to the theory

of games, Vol. 1, Annals of Math. , Studies, 2*+, 1950.

75. H. W. Kuhn and A, f, Tucker, eds., Contributions to the theory

of games, Vol. 2, Annals of Math. Studies, 28, 1953.
THE NATIONAL BUREAU OF STANDARDS

Functions and Activities

The functions of the National Bureau of Standards are set forth in the Act of Congress, March

3, 1901, as amended by Congress in Public Law 619, 1950. These include the development and

maintenance of the national standards of measurement and the provision of means and methods

for making measurements consistent with these standards; the determination of physical constants

and properties of materials; the development of methods and instruments for testing materials,

devices, and structures; advisory services to Government Agencies on scientific and technical

problems; invention and development of devices to serve special needs of the Government; and the

development of standard practices, codes, and specifications. The work includes basic and applied

research, development, engineering, instrumentation, testing, evaluation, calibration services, and

various consultation and information services. A major portion of the Bureau’s work is performed

for other Government Agencies, particularly the Department of Defense and the Atomic Energy

Commission. The scope of activities is suggested by the listing of divisions and sections on the

inside of the front cover. Reports and Publications

The results of the Bureau’s work take the form of either actual equipment and devices or published papers and reports. Reports are issued to the sponsoring agency of a particular project

or program. Published papers appear either in the Bureau’s own series of publications or in the

journals of professional and scientific societies. The Bureau itself publishes three monthly periodicals, available from the Government Printing Office: The Journal of Research, which presents

complete papers reporting technical investigations; the Technical News Bulletin, which presents

summary and preliminary reports on work in progress; and Basic Radio Propagation Predictions,

which provides data for determining the best frequencies to use for radio communications throughout

the world. There are also five series of nonperiodical publications: The Applied Mathematics

Series, Circulars, Handbooks, Building Materials and Structures Reports, and Miscellaneous

Publications.

Information on the Bureau’s publications can be found in NBS Circular 460, Publications of

the National Bureau of Standards ($1.00). Information on calibration services and fees can be

found in NBS Circular 483, Testing by the National Bureau of Standards (25 cents). Both are

available from the Government Printing Office. Inquiries regarding the Bureau’s reports and

publications should be addressed to the Office of Scientific Publications, National Bureau of Standards, Washington 25, D. C.