4—
Heuristic Studies in Problems of Mathematical Physics On High Speed Computing Machines:
With John Pasta (LA-1557, 1953)

This report introduces some elementary but basic methods of a specific version of "brute force" approaches in problems of hydrodynamics which do not yield to analytical methods.

The ideas in this report were further developed, generalized, and modified by Frank Harlow, and continue to play an important role in modern computer calculations in continuum mechanics. (Author's note.)

This paper, the first of a series, is intended primarily to illustrate some possible uses of electronic computers as a means of performing "mental experiments" on mathematical theories and methods of calculation, for a variety of physical phenomena.

There is no unifying principle in the few problems selected for this first report; on the contrary, the aim is to illustrate a certain variety in the problems which one can consider by performing model calculations. Mathematically, these problems are mainly attempts to solve special cases of various partial differential equations; in some cases one prefers to have direct recourse to the physical formulation and perform, as it were, the model experiments on paper rather than solving, by standard conversion of differential to difference expressions, the equations themselves. One general tendency might be notedwe stress an attempt to "observe" a few functionals of the unknown

functions, rather than to put credence in the solutions themselves "at each individual point."

The calculations were performed on the Los Alamos electronic computer. It is a pleasure to express our thanks to N. Metropolis for making it available and for his generous help in general. We are indebted to Miss Mary Tsingou who performed much of the work of "coding" the problems, running the machine and checking the procedure and results.

One purpose of this work was to gain a feeling for: first, the time necessary to formulate and prepare the problems for such machine study; second, the time of computing as a function of the size and complexity of the formulation. In general, we restrict ourselves in this first paper to problems where meaningful results could be obtained in a few hours of actual computing.

All the problems discussed in this paper were run on our machine. We append a small selection of the results which we intend to discuss at greater length in our second paper. We want to thank Mrs. Connie Snowden for preparation of the graphs.

1—
Hydrodynamical Problems; Heuristic Considerations

In general for problems involving two or more spatial variables, there is at present little hope for obtaining solutions analytically-in closed form. In the several sections that follow, heuristic considerations are set forth with exploratory calculations on the high speed computing machines performed in several cases. The purpose of these was primarily to establish the feasibility of such calculations, to estimate the sizes of problems which can be handled in a reasonable time, and in general to gain experience in the new methods and new fields which, in our opinion, are now open for investigation.

Our approach to the problem of dynamics of continua can be called perhaps "kinetic"--the continuum is treated, in an approximation, as a collection of a finite number of elements of "points;" these "points" can represent actual points of the fluid, or centers of mass of zones, i.e., globules of the fluid, or, more abstractly, coefficients of functions, representing the fluid, developed into a series, e.g., Fourier or Rademacher series. This corresponds to the use of general Lagrangian coordinates in classical mechanics; their use can be rigorously justified in problems where entropy is constant-i.e., holonomic systems.

They are always functions of time, which proceeds by discrete intervals.

We thus replace partial differential equations or integer-differential equations by systems of total differential equations. The number of

elements which we can at present handle is always less than 1,000, the limitation being primarily in the "memory" of the machines, the number of time intervals ("cycles") of the order of 100.

The problems which we study are characterized by lack of symmetry. The positions of our points become, as a rule, more and more irregular as time goes on. This has the consequence that the meaningful results of the calculations are not so much the precise positions of our elements themselves as the behavior, in time, of a few functionals of the motion of the continuum.

Thus in the problem relating to the mixing of two fluids, it is not the exact position of each globule that is of interest but quantities such as the degree of mixing (suitably defined); in problems of turbulence, not the shapes of each portion of the fluid, but the overall rate at which energy goes from simple modes of motion to higher frequencies; in problems involving dynamics of a star cluster, not the individual positions, but quantities like angular momenta of subsystems of the whole system, of size smaller but still comparable to the whole system, etc.

Needless to say, our investigations are in a most preliminary and rudimentary state; we have not made rigorous estimates of the disturbing or smoothing action of the roundoff errors which accumulate, nor the effects of finiteness of the time interval (i.e., the errors due to replacement of differential by differences expressions). In the individual cases that follow the reader will be able to judge for himself how far elementary common sense permits to estimate these effects.

We hope in the future to multiply considerably the number of calculations performed for each of our problems to assay the influence of changes in initial conditions on our conclusions. We repeat that so far the value, if any, of this work is only heuristic.

The first problem considered concerns the behavior of a gas confined in a vessel, expanding into vacuum under its own pressure and weight. The surface is not plane but has an irregularity of a finite size. In other words we consider the problem of instability in the compressible case.

2—
Instability and Mixing

The conditions at time t = 0 are the following: a vessel (twodimensional) contains a gas filling it partly, with a boundary against a vacuum. This boundary is not flat but has an irregularity in it in form of a triangular prominence jutting out with dimensions comparable to the diameter of the vessel (about one-fourth of its width).

We want to follow the behavior of the expansion of this gas, under its own pressure, in the vacuum below it. Two problems were considered:

1) the gas was assumed weightless, i.e., there was no external force acting on it;

2) the gas, in addition to its own pressure, was acted upon by a constant gravitational field.

The hydrodynamical setup used was the following: the gas was represented by 256 material points. These represent centers of mass of regions in the gas. The treatment is Lagrangian; that is, the calculation follows, in time, the position of these masses. The pressure gradients are represented by forces which our points exert on each other. These are repulsive forces, depending only on the distance between points (and thus having a potential "pressure") and are of the form Fij = a/rjI; where the exponent x depends on the adiabatic equation of the state of the gas considered.

We chose oc = I for our first problem. Given a point Pi (xi,yi), one considers all forces exerted on it by other points (Xj, yj) and computes their resultant vector.

Actually, we limit the pj in computing the forces to the "neighbors" of pi; these are defined as the nearest eight points to pi. This is done for two reasons: the economy of the computation-we have to calculate only eight instead of 256 of the total number of forces for each point under consideration; the second, more fundamental reason is, of course, that in the gradient of pressure only the local configuration matters. The 1/r force law gives divergence at infinity. We might mention here parenthetically that in the search or scanning of points for the nearest eight to the given one, the following was adapted for purely practical, economy, reasons. The points really scanned were 50 "candidates" for the nearest neighbor. They were the 50 originally nearest to pi; the problem was not run long enough so that we had to relocate the original candidates but it is, of course, possible to recheck this periodically. In addition, in order to avoid the use of multiplications and operate

only with the much shorter addition times, we used, in the search for the nearest eight points, the non-Euclidean metric p(vi, vj) =I(xi xj) I + I (Yi- yj) Once the eight points were found, however, the true Euclidean distances were computed in order to find correctly the resulting force.

We shall not describe here the special treatment which has to be accorded to points which adjoin the walls of the vessel and the points on the boundary of the gas (with the vacuum).

Among the quantities that were printed as the results of the calculation we shall mention here only these: an interesting functional of the motion is the kinetic energy of the gas divided into two parts: the energy of the motion in the x direction and in the y direction. We study E1 = 1/2 Emi.² and E2 = 1/2 Emi]^2.

The ratio of the two is a function of time and can serve, in a way, as a measure of mixing or irregularity of the motion. One expects due to the initial irregularity of the boundary this ratio to be positive. From the beginning, sidewise motions ensue. Later on one would expect, of course, the motion to be predominantly downward; as the irregularity increases, the ratio of the two quantities should, in Problem I, increase again and approach a constant less than unity.

It is perhaps remarkable that the time behavior found for this ratio was extremely regular; a graph (Fig. 1) for the first 36 cycles is appended.

One word to explain the need to resort to the rather unorthodox procedures outlined above:

It was found impractical to use a "classical" method of calculation for this hydrodynamical problem, involving two independent spatial variables in an essential way (since the gas interface had an irregularity assumed from the beginning). This "classical" procedure, correct for infinitesimal steps in time and space, breaks down for any reasonable (i.e., practical) finite length of step in time. The reason is, of course, that the computation of Jacobians which define the compression assumes that "neighboring" points, determining a "small" area, stay as neighbors for a considerable number of cycles. It is clear that in problems which involve mixing specifically this is not true. Calculations were made just in order to observe the rapidity of change in the "neighbor" relations on the classical pattern and have shown just what was expected to happen: the proximity relations change radically, for points near the boundary, just when the mixing to be studied is starting and the neighboring relation of our points has to be redefined; i.e., the classical way of computing by referring to initial (at time t = 0) ordering of points becomes meaningless.

The problem can be treated, of course, using the Eulerian variables (of a set-up due to von Neumann to be discussed in our next report)* where this difficulty does not occur. This Eulerian treatment

* There was no subsequent report. (Eds.)

is not suitable, however, for the study of the shape of the boundary-a fictitious, i.e., purely calculational, diffusion and mixing obscures the very phenomenon one wants to study. We shall return to this question in our second report.

In the "classical" set-up involving calculation of Jacobians for determining pressure gradients the cycle was approximately two minutes.

In the problem involving the calculation of forces from the eight "neighbors," the cycle time was three minutes. Meaningful results are obtained in about 100 cycles. Appended are graphs showing the results of a few dozen cycles.

3—
Billowing Transformations

The aim of the model calculation here is to exhibit transformations of space which, when applied iteratively to a sphere (or, actually, so far, to a circular region in a plane) will show a sequence of regions, imitating the familiar phenomenon of a ball of smoke billowing outwards. The parameter of iteration is, of course, the time again increasing by discrete intervals.

We assume that the billowing is due to initial irregularities or deviations from a spherical form.

The gas is contained in a region R, and is under internal pressures. We assume first, for simplicity, that it expands into vacuum and we want to study the motion in a highly stylized form-keeping, we hope, the quantitative phenomena essentially correct. The situation then is the following: R is the region occupied by the gas under pressure in the initial position. The following assumptions are made and made plausible qualitatively:

1. The accelerations of points on the boundary are in the direction of the outer normal to the boundary.

2. The motion is computed only for the points on the boundary above; the accelerations at each point depend on the shape of the boundary which is approximately correct in the following approximations:

The density distribution and the pressure well inside the ball is essentially uniform-if the time elapsed since the origin of our ball is long enough-i.e., the instantaneous changes in the position of the particles on the boundary are small compared to the dimensions of our object. "The pressures inside have had time to equalize."

3. On the boundary itself-the accelerations along the normal to each point depend on the curvature at each point in the following way:

There is a positive term outward in n with an additive term whose sign is that of the curvature c of the boundary at the point which is considered. This is due to the convergence (or divergence) effects of the streamlines.

E.g., in case 2 the local density at P is probably greater than in the nearby points because of converging motions of points in a band near the surface; in case 1 there is correspondingly a local rarefaction. So we set ntt=p + ac (1) where c is the curvature of the boundary. The terms p and a are not constant but vary, say, with the volume enclosed by the boundary-the pressure diminishing, e.g., inversely with the volume.

A qualitatively similar formulation, (in a discussion with John von Neumann) in the small, would be this: consider the part of the curve given by y = f(x)

f is positive and the direction of expansion is upward. We write Yttt = Yxx (2)

with the initial conditions at t = 0, ytt = 1, yt = 0, say, and y = f(x). The tendency is again the same in convex points. The accelerations

upwards are decreased, while in concave ones (from which one expects "jetting") they are increased.

Setting y =- (t) X (x) and say X (x)= cos (nx - /), one gets f"'(t) = -n ²). Setting ) = e"(t), a is the cubic root of -n2^. It follows, since the dominant part will be due to the root with positive real and imaginary parts that the motion is highly unstable. A kind of billowing or swirling will take place--the concave parts will puff out, giving rise later to at least two concavities on its sides. This will repeat the same phenomenon later, etc. It would seem then that a multiplication of irregularities takes place-their number will increase rapidly-accompanied with a continuous increase in size.

In three dimensions the problem is much harder to treat. The additional term in the acceleration must be set up in terms of the two principal radii of curvature at each point.

Our computation program started in two dimensions as follows. Initially the boundary is taken to be composed of sixty-four equal segments defined by points (xi, Yi) on the surface. The magnitude and sense of the curvature at (xi, yi) is derived from the cross product of the two vectors (xi - xi,_, i - yi-1i) and (xi+ ₁- xi, Yi+1 - yi). From the coordinates (x, y) at two consecutive time levels the positions at the next time level are found by integrating ntt = p +ac.

Because of the difficulty of plotting many points every time and in order to observe the motion, it was decided to display the points on one of the cathode ray tubes of the computer memory section. Each such tube has a 32 x 32 array of points, each representing a binary digit in one of the machine's 1024 forty bigit numbers. There are forty such tubes. The most convenient tube to use is tube number one, the "sign bigit" tube. For the proper display it is convenient that all constants of the problem be taken as positive, which may be done without loss of generality. In the computer used here, machine orders are independent of the bigit in the sign position. Thus, the abscissa and ordinate of a point are transformed into an instruction for changing the sign of the appropriate memory position, which, in turn, lights that spot on the tube face. In this way a picture of the boundary is "painted" on the tube face. The picture is displayed for a fixed time, is erased by dropping the sign of all memory positions, and a picture of the new position is then plotted and displayed. A scaling routine keeps the

surface within the confines of the tube matrix.

This procedure should be useful in many problems of motion of gases or liquids.

It is certainly of value in obtaining a quick overall check on the correctness of the code; reasonableness of the time intervals, etc.

The cycle time on this problem is of the order of ten seconds. No printing is involved. Meaningful results are obtained in a matter of minutes.

4—
Problems on Rotational Motions in Gravitating Systems

An interesting set of questions in statics concerns the properties of moments of forces exerted on each other by randomly distributed points forming a system S. In dynamics the questions concern angular momenta of subsystems a contained in S; as a function of the time.

Let us imagine the following situation: E consists of a number of mass points mi, m2, ... mn located at t = 0 at positions rl,... rn given at random, say, in a unit sphere (let us assume, for example, a uniform probability distribution for the position of each point). We assume further Newtonian attractive forces Fij acting between any two points mi, mj. Denote by Gi the sum over all j of forces acting upon the point mi and let us imagine the vector Gi applied at the point ri. It is clear that the sum over all i of Gi is equivalent to zero. What we propose to study is, at first, the statistical behavior of the forces Gi if summed over subsystems a of the whole set with the following questions to be investigated: let p be any number and let us consider subsystems located in a circle with radius p and an arbitrary center ro. Let us form the sum of all Gi located in such and we obtain, in general, a single force $ and a couple P referred to ro with magnitude which we shall call 0. Both ( and 0 thus computed are functions of p and ro but we can integrate these quantities over all initial positions ro and will obtain )p and Op, a single force and moment, which will be now functions of p alone. It is our aim to obtain these functions for a random dynamical system, that is to say, the expected values in a random distribution. This can be obtained in practice by computing, on the machine, these quantities for a large number of systems each chosen by a random process. Our statistical computations will be confined at first to plane

cases, the three dimensional systems requiring too great a memory at present.

The next, more interesting, thing to study is the following. In a situation as described above at time t = 0 let t increase. Motions of our

points will ensue and we intend to investigate the angular momenta of subsystems o as functions of the size of a and of $q in time.

Such a situation is perhaps exemplified by star clusters. What we want to study are dynamical systems with many particles, but not gases; that is to say, by mean free path for "collision" we mean an appreciable change in the velocities due to gravitational forces acting between just two of our mass points. It is known that clusters or galaxies possess a rotational motion as a whole. These could, perhaps, originate as follows: the original distribution of matter now separated in galaxies was more or less uniform and random-like. Our system E can be imagined infinite. Finite subsystems have angular momenta as a result of fluctuations in the distribution and then, if fluctuations of density occurred also, some subsystems would isolate themselves, stay together due to gravitation (cf. the work of Jeans) and if the whole space expanded with time these condensations would have kept all or most of their angular momenta due to the original fluctuations in the system of vector forces and as the condensations receded from each other, their non-zero angular momenta would have stayed constant in time. Another way to look upon the problem is to study the distribution of vorticity of finite subsystem a in a very large or infinite system E whose points exert forces on each other.

Our proposed calculations consist then of producing a sizeable number of randomly chosen systems E and following the behavior of subsystems in time, i.e., using a discrete series of time intervals or "cycles" on the machine and computing the following set of averages: L(ro, p), the angular momentum as a function of the radius of the subsystem p and position ro. For a system of 100 points four positions along a radius with four values of p at each point are adequate. The running time per cycle for 100 mass points is less than five minutes with this machine. The running time increases as the square of the number of mass points, but statistics can be improved by running many problems with few mass points, in preference to the less economical method of increasing the number of mass points.

The total kinetic and potential energy is calculated. The system should stabilize at some given size. In this steady state the number of double and triple "stars" is of interest. The total energy serves as a check on the problem.

The moment of inertia of the complete system is calculated, and finally the number of particles in each of the sixteen subsystems mentioned above. The latter numbers permit an approximate plot of the density distribution of the system. The largest value of p is chosen to be of the order of the dimension of the system.

Several hundred cycles will be necessary for results on this problem.

Appended are graphs (Figs. 2 to 5) showing the results of 100 cycles run so far.

5—
Magnetic Lines of Force

It appears that our computing machines are especially well adapted to the study of properties in the large of the system of lines of magnetic force, due to given currents in space. The renewed interest in the qualitative, ergodic, or even just topological behavior of such families of lines is due to studies in magneto-hydrodynamics, applications in astronomy, questions of origin of the cosmic rays (Fermi), not to mention the importance of such knowledge for applications in the construction of high energy accelerating machines (cyclotrons, synchrotons, etc).

In order to study this subject systematically, it is best to consider first steady currents following through given lines (wires) in space. If there is only one current through a straight line extending infinitely, the system of lines of magnetic force is, of course, very well known. They form circles linking the line; the same is, or course, true for a current in a single closed circle.

However, the topology of the system of lines of force seems to be very complicated and the ergodic behavior of single lines of force unknown for the case where the single closed curve, through which the current flows, forms a knotted loop, say in the simplest case a clover leaf knot. Some single lines of force will probably be dense on two dimensional surfaces, probably some singularities in the field of lines exist, independently of the metric appearance of the knot, but are present necessarily in every topological knot of this sort. There seems to be little hope of obtaining analytically closed expressions describing the system of lines of force.

The situation is complicated in case of two given currents. It is easily seen, in the case given below, that almost all lines of force will be bounded, not closed and each dense on a surface! Let the two currents flow as follows: current 1 on a straight line, say the z axis, current 2 on a circle x² + y² = 1. In general, except at points where the ratio of the two current strengths is rational, a line of force will exhibit an "ergodic" behavior on a surface of a torus.

We propose to investigate, on the computing machines, the properties of lines of force due to two currents-each on a straight line extending to infinity. The two lines are skew. T₁ flows on the y axis, T2 on the line z = d, y = 0.

One is interested, among other things, in the following questions: Do there exist lines which, although not closed, cross a surface of a fixed sphere infinitely many times? Are there lines going arbitrarily far from a fixed point and returning to its neighborhood? Do there exist lines braiding or linking both given wires any number of times?

The computations of such lines of force do not involve much of the "memory" of the machine. The procedure is this: starting at a point (xo, yo, zo) we compute the direction of the magnetic fields, simply adding the two field strengths, given elementarily from each wire; we perform "short" step (Ax, Ay, Az) in the direction of the (constant in time!) force. The computation of this step is done as follows: we calculate a provisional set of increments (Ax)', (Ay)', (Az)' of the variables, in the new position we calculate the new set (Ax)", (Ay)", (Az)"; we then take Ax = ((Ax)' + (Ax)")/2, Ay, = ((Ay)' + (Ay)")/2, Az = ((Az' + (Az)")/2 and proceed anew. In general this way of solving a system of equations dx/X = dy/Y = dz/Z works well. In our case it is seen that each step is computed in the order of 50 milliseconds; to perform, say, 1,000 steps will take of the order of a minute. The idea is now that with the order of 10⁴ steps we shall be able to get

some qualitative information about a single line of force as follows: it seems practical to take each step long enough so that in, say, 50 to 100 steps one complete "loop" can be described around a wire (in positions where it is expected that the lines of force surround the current). One would expect then to obtain a number of "loops" of the order of a few hundred.

The quantities printed as a result of each such calculations could be for instance:

1. The number of "returns" of a line to a given sphere. One simply has to record on the machine the number of times our line crosses the surface of a given sphere.

2. The number of times and the sense in which a line loops the two wires separately and the number of loops surrounding both together. This can be done simply by computing the work done by moving on the line of force, calculating the loops around each wire, as if it alone had current flowing through it. The Gaussian looping coefficients for any two given curves in space can be quickly computed on the machine.

It is convenient to take the length of the "step" along the magnetic field vector to be inversely proportional to the magnetic field strength at that point. In this way the step is proportional to the distance from the wire since the magnetic field is inversely proportional to that distance. The number of steps per looping of the wire is then constant and the step is appropriately shorter at points where the curvature is greater.

It is worthwhile to point out that some integer valued topological invariants may be computed exactly even though we use difference expressions instead of differential ones, and have also round-off errors. This is due to "e-invariance" theorems on simplicial approximations in topology¹ . Also, so to say, the field of "error vectors" is in general "curl-free."

Reference

1. K. Borsuk and S. Ulam, "Uber Gewisse Invarianten der Abbildungen," Math. Ann. 108, 311-318, 1933.

blank

Problem Of The Attracting Points Fig. 2.

Fig. 3.

Fig. 4.

Fig. 5.