Interim CAPOW Project Report
November 10, 1994.
Report on a project, involving CA modelling of
the wave equation, being carried out at SJSU Math & CS CAMCOS, (San Jose
State University Department of Mathematics and Computer Science Center for
     Rudy Rucker

     Department of Mathematics and Computer Science

     San Jose State University, San Jose CA 95192

Cellular automata are an intrinsically parallel kind of computation. A single-processor machine can be thought of as a point-like structure; a computing point, also known as a zero-dimensional cellular automaton. A one-dimensional cellular automaton is like a string or wire of computing points. A two-dimensional CA is a scarf of computation, and a threedimensional CA is a computational jelly.

A CA is a manifold of computing cells which repeatedly update their internal states. The update process is characterized by parallelism, homogeneity, and locality. Parallelism means that all the cells are updated at the same time, in lock-step synchronization. Homogeneity says all cells use the same update rule; and locality says cells only gather information from their nearby neighbors.

The cellular automata studied by mathematicians and computer scientists are closely related to the finite element methods used by engineers. A goal of the CAPOW project is to apply cellular automata to the specific engineering problem of simulating alternating current in a power grid. It is hoped that we can develop technology that will be of use for studying power grids, and that the lessons learned from this endeavor will advance the theory and practice of cellular automata.

Much of the mathematical theory of cellular automata focus on the large-scale dynamical behaviors of classes of CAs; this involves looking at many different kinds of CAs. As a practical matter, it is hard to get large CA computations to run rapidly, so mathematicians have mostly looked at CAs in which the cells have only a small amount of information --- the famous Game of Life CA uses only one bit of internal state per cell. Computer scientists have done a lot of work on getting cellular automata to run rapidly on desktop machines, but still most of the CAs investigated have fairly simple internal state; eight bits might be typical.

But in engineering finite element simulations, the cells typically have real-number values that occupy 64 or more bits of information. So in the CAPOW project, we want to look at CAs with states based on real numbers.

We began our project by constructing an object-oriented Windows program that allow us to work with cellular automata. This graphically interfaced program, called CAPOW.EXE, allows the user to view multiple CAs, to save and load them, to alter their parameters, to seed them, to breed them, to measure their entropy, and so on.

As a starter we modeled CAs with small numbers of states; some look at only two neighbors, while some look at up to six neighbors. Then we began adding real-number states.

The idea is to first model the one-dimensional wave equation corresponding to, e.g., a vibrating string, a slightly perturbed water surface, or a sound wave, then to add a driving element to the system, and then to try and simulate alternating electrical current, and finally to simulate power transmission and distribution networks, complete with generators and loads. I will now describe the wave-simulation methods that we are currently investigating.

A Cellular Automaton Version of the Wave Equation:

Think of a particular one-dimensional cellular automaton with discrete cells located next to each other. Each cell has internal structure, and each cell can sense the states of its nearest neighbors. How many neighbors? To begin with, let's have one of our cells only look at its two nearest neighbors; the one on the left and the one on the right.

As we want our array of cells to represent a continuous medium, we think of the cell array as being laid on top of a real-number-line.

We use DX to stand for the continuous spatial step corresponding the distance across a single cell. That is, if there are M cells and the line under the cells has length LengthX, then DX is LengthX/M.

As concerns internal structure, let's suppose that each cell has two floating-point real-number values which we call intensity and velocity. Sometimes we refer to these values as u and v for short.

The CA is initialized by placing some values in the intensity and velocity fields. One simple initialization would be to set all values to 0.0, except for setting the intensity of one cell to some non-zero value. Another way to seed the CA might be to put in random (bounded) values for each cell's intensity and velocity.

As time goes by, the CA is repeatedly updated. The parallelism requirement means that the update happens all at once, in every cell at the same time. (In computational practice, we use buffers to simulate parallelism on our serial 486DX machine.)

In displaying the CA we represent a cell by a pixel, with the color of the pixel derived from the cell's state. The entire state of the CA is represented by a horizontal scan line with one pixel for each cell. As time goes by and the CA is updated, we show the updates as successive scan lines further and further down the screen.

We use DT to stand for the continuous temporal step corresponding the gap between two scan-lines. That is, if you perform N updates and the elapsed time is to represent a duration DurationT, then DT is DurationT/N.

The wave equation says:

The second derivative of intensity with respect to time equals

c-squared times the second derivative of intensity with respect

to space. If we write u for intensity, this takes the form: [Note that in this ASCII printout, some of the 2's should be superscripts!]

  1. d2u / dt2 = c2 * d2u / dx2

    The constant c represents the wave velocity, that is, the maximum velocity for propagating a disturbance in the medium. In other words, c can be thought of as a "speed of light" or a "speed of sound."

    In studying one-dimensional CAs, we often speak of "the speed of light" as being one, because if a particular CA looks only at its nearest neighbor on either side, than no effect can move across the CA faster than one cell per update. If, however, a CA happens to look at its two nearest neighbors on either side, then the "speed of light" will be two cells per update, and so on. As a final possibility it may be that a cell looks at its nearest neighbors, but it takes two time steps for the effects to fully percolate into the cell. This is in fact the case with our CA's cells. They look at their nearest neighbors, but effects take two steps to propagate. This means our wave velocity will one-half cell per update. In terms of DX (the width of a cell) and DT (the time step per update), this makes c equal to DX divided by 2*DT. If we were looking at a CA which the effects were immediate, then c would be DX/DT; while if we had immediate effects and looked at two neighbors on each side, then c would be 2*DX/DT, and so on. Replacing C by DX/2*DT changes the wave equation (1) to:

  2. d2u / dt2 = (DX/2*DT)2 * d2u / dx2

    Recall that each of our cells has an intensity u and a velocity v. We plan to maintain our CA so that a cell's velocity represents the rate at which the cell's intensity is changing. That is, we require that our CA settle into a state such that for each cell du/dt = v. Plugging this into

    the wave equation (2), we get
  3. dv/dt = (DX/2*DT)2 * d2u / dx2

    How can we get a handle on d2u / dx2? This is the Laplacian. For a first-order approximation to it, we can think of a 1D CA cell with intensity C, with left neighbor of intensity L and right neighbor of intensity R. The derivative of intensity with respect to time can be approximated as (C - L)/DX or as (R - C)/DX, according to which end of the cell one looks at. The second derivative of intensity with respect to time

    must then be (((C-L)/DX) - ((R-C)/DX))/DX), or

(4) (L -2C + R)/DX2.

Suppose we call the quantity on the right-hand side of the wave equation the acceleration. Then the acceleration is approximated by:
  1. (DX2)/(4*DT2) * (1/(DX2))(L -2C + R), or (5a) (1/(4*DT2))(L -2C + R).
The update rule should be:
  1. newintensity = oldintensity + DT*velocity.
  2. newvelocity = oldvelocity + DT*acceleration;
The first equation can be rewritten: (7a) newvelocity = oldvelocity + DT*(1/(4*DT2))(-2C + L + R), or (7b) newvelocity = oldvelocity + (1/4*DT)(L -2C + R).

So the rules for doing a CA version of the wave equations are:

  1. newintensity = oldintensity + DT*velocity.
  2. newvelocity = oldvelocity + (1/4*DT)(L -2C + R);
This approach was in fact developed second, and we are still experimenting with it. At this point it seems to give the kinds of effects we expect: wave fronts, reflection, pulsing from a generator, and orgainzation into domains when started with a random seed.

Higher Order Neighborhood:

We can use a higher order technique to approximate the Laplacian d2u / dx2. If we think of five successive cells with intensities LL, L, C, R, and RR, with an x step of DX per cell, then the second-order Laplacian at position C has the following form instead of equation (4):

(10) (-LL + 16*L - 30*C + 16*R - RR) / (12*DX2)

This Laplacian was computed by finding the formula for the degree four Bezier curve through the five values, taking the formula's second derivative, and evaluating the second derivative at the center cell.

Before applying (10), note one other consideration. The "speed of light" c value should be DX/DT in a situation where a wave CA looks at two neighbors on each side. This is because an effect can move two cells per update but, as before, it takes two updates for one intensity to affect a neighboring intensity. Replacing c by DX/DT and using (1), we get this expression for the right-hand side of the wave equation, in analogy to (5):

(11) (DX2 /DT2) * (1/(12*DX2))*(-LL + 16*L - 30*C + 16*R - RR)

(11a)          (1/(12*DT2))*(-LL + 16*L - 30*C + 16*R - RR)
And now the update equations are
  1. newintensity = oldintensity + DT*velocity. (13)newvelocity = oldvelocity + (1/(12*DT))*(-LL + 16*L - 30*C + 16*R - RR) I implemented this rule recently, and it works fine, reproducing the
same results as the order-one Laplacian.

A Non-Physical False Path:

     The approach in this section was the first technique we tried.  It is
physically incorrect, but perhaps of mathematical interest. At the least, it makes some some interesting patterns. We spent several months looking at this approach; and in the process of trying to figure out what was happening we designed and built a lot of software tools that are proving useful for investigating the correct wave equation CAs described above.
     Once again we give each cell two real number values, intensity and
velocity. The new array is updated from the old array by four commands involving four parameters belonging to the CA: acceleration_multiplier, velocity_multiplier, max_velocity, and max_intensity.

new[i].velocity = old[i].velocity + acceleration_multiplier *

     (old[i-1].intensity + old[i+1].intensity - 2.0 * old[i].intensity);
new[i].intensity = old[i].intensity + velocity_multiplier *
Clamp new[i].velocity to lie between -max_velocity and max_velocity; Clamp new[i].intensity to lie between -max_intensity and max_intensity.

acceleration_multiplier, velocity_multiplier, max_velocity, and max_intensity are tweakable parameters for the CAs, these are in fact the genome that we use for our evolutionary CA experiments in finding a better wave. The ranges of values investigated for the _multiplier variables are from -0.5 to 2.0. The investigated ranges for the max_ variables are from 1 to 10,000.

We seeded in any of three ways: set all cells' velocity and intensity to zero except for one cell which has its intensity set to max_intensity, randomize all intensity and velocity values, or put the intensities and velocities into a sine/cosine pattern with arbitrary tweakable values for the pattern's periodicity.

We observed beautiful waves, but they seem non-physical as there is a seemingly unremoveable numerical chop in them, meaning that within the largescale oscillations there is a typical maximum to minimum oscillation of the intensity. We carried out by hand some numerical simulations that confirm this is a real effect for the rule as formulated above. This phenomenon seems to be related to the physics of driven oscillators and of beats in waves.

How can I simulate wind in a CA?
From: Bruce Boghosian <>
From: (Kenneth Tolman)
From: Rudy Rucker <rucker@jupiter.SJSU.EDU>
For a lattice gas treatment of incompressible, convective flow in the Boussinesq approximation, see C. Burges and S. Zaleski, "Buoyant mixtures of cellular automaton gases," Complex Systems 1 (1987) 31.

{\bf Rucker:} A simple way to make wind is to just have each cell copy the cell to its right with each update. This makes a wind that blows to the left. In the RC module of the CA LAB package, I ran this rule alternating with the heat equation (rug rule). I also allowed there to be some stable blocks of cells in the pattern which were *not* shifted to the left with each update --- these were copied in as masks with each update. The effect was of turbulence in the wind behind the blocks, vortices, slipstreams, the whole shmear. Try it you'll like it.

What are Cellular Neural Nets?

Frank Puffer

A CNN is practically
a CA with continuous states that may be disrete or continuous in time. (A good overview article on CNNs is available via ftp at in /pub/cells/

What are the commonly used rules of interaction between the cells in a lattice gas? What sort of fluid behavior can and cannot be modelled by them?

From: (Bruce Boghosian)
From: (Christoph L Mikkelson)
From: (Dave Stack)
From: Mohamed <>
There's a huge literature on this. Briefly, examples include two- and three-dimensional Navier-Stokes flow, magnetohydrodynamics, immiscible fluids with a surface tension interface, convection, two-phase liquid-gas flow, Burgers' equation, and reaction-diffusion equations. For a summary of recent works on the subject, see any of the reviews listed below; also there should be a very recent (or imminent) Physics Reports devoted to lattice gases by Rothman and Zaleski.
\bibitem{doolena} ``Lattice Gas Methods for Partial
   Differential Equations,'' Proceedings of the Workshop on Large
   Nonlinear Systems, held in Los Alamos, New Mexico, August, 1987,
   Proceedings Volume IV of the Santa Fe Institute Studies in the
   Sciences of Complexity, Addison-Wesley (1990), Doolen, G.D., ed.
\bibitem{doolenb} ``Lattice Gas Methods for PDE's, Theory,
   Applications and Hardware,'' Proceedings of the NATO Advanced
   Research Workshop, held at Los Alamos National Laboratory, September
   6-8, 1989, North-Holland (1991), Doolen, G.D., ed.
\bibitem{monaco} Proceedings of the Workshop on ``Discrete
   Kinetic Theory, Lattice Gas Dynamics and Foundations of
   Hydrodynamics,'' held in Torino, Italy, September 20-24, 1988, World
   Scientific (1989), Monaco, R., ed.
\bibitem{manneville} ``Cellular Automata and Modeling of Complex
   Physical Systems,'' Proceedings of the Winter School, Les Houches,
   France, February 21-28, 1989, Springer (1989), Manneville, P.,
   Boccara, N., Vichniac, G., Bidaux, R., eds.
\bibitem{alves} Proceedings of the Colloquium Euromech No.
   267 on ``Discrete Models of Fluid Dynamics,'' held in Figueira da
   Foz, Portugal, September 19-22, 1990, World Scientific (1991), Alves,
   A.S. ed.
\bibitem{nicemtng} Proceedings of the NATO Advanced Research Workshop on
   Lattice Gas Automata:  Theory, Implementation and Simulation, held at
   l'Observatoire de la C\^{o}te d'Azur, Nice, France, June 25-28, 1991,
   to appear in {\sl J. Stat. Phys.}, Boon, J.-P., Lebowitz, J.L., eds.
One of the seminal papers is "Lattice Gas Hydrodynamics in Two and Three Dimensions" by Frisch et. al. in Complex Systems Volumne 1 (1987) 649-707.

There is an excellent 3 part review paper by Brest Hasslacher in a special issue of Los Almos Science in 1987. The title is "discrete fluids" and the the title of part I is "background for lattice gas automata". Also see Computer in physics, Nov. 1991, page 585 for an article by B.M. Boghosian.

Viscosity in LGA?

From: (Bruce Boghosian)
From: (Dollinger Juergen)

Relationship between the shear viscosity
and the collision probabilities.
Stephen Wolfram,
``Cellular Automaton Fluids 1: Basic Theory'', {\em Journal of Statistical Physics} {\bf 45} Nos.~3/4 (1986) 471--526.
(equation 4.6.8 )
U.~Frisch, D.~d'Humi\`{e}res, B.~Hasslacher, Y.~Pomeau, J.~Rivet,
``Lattice Gas Hydrodynamics in Two and Three Dimensions'',
{\em Complex Systems} {\bf 1} (1987) 649--707.
(equation 8.25)
From general kinetic theoretical arguments, the viscosity of a fluid goes as the product of the thermal velocity and the mean free path. Now the thermal velocity of the particles of a lattice gas automaton is fixed at one lattice spacing per time step. It follows that the viscosity goes as the mean free path. Thus, decreasing the collision probability clearly increases the mean free path, and hence the viscosity.