Interim CAPOW Project Report

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

     rucker@jupiter.sjsu.edu

Introduction:

     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!]

 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:
 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

 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:

 (DX2)/(4*DT2) * (1/(DX2))(L -2C + R), or
(5a) (1/(4*DT2))(L -2C + R).

The update rule should be:

 newintensity = oldintensity + DT*velocity.
 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:

 newintensity = oldintensity + DT*velocity.
 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


 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 *
     old[i].velocity;

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.


 This project is being carried out at SJSU Math & CS CAMCOS, (San Jose
State University Department of Mathematics and Computer Science Center for
Applied Mathematics and Computer Science) with the support of a grant from
EPRI (the Electric Power Research Institute) in Palo Alto.  Student team
leaders: Juyong Lee and Alan Borecky; team members: Tuyen Ly, Robert
Westergaard, Howard Lin, Andy Chu, and Charles Miller.

How can I simulate wind in a CA?
From: Bruce Boghosian <bruceb@conx.bu.edu>
From: tolman@asylum.cs.utah.edu (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 volterra.science.unitn.it in /pub/cells/cimagalli-review.ps.Z.)

  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: bruceb@conx.bu.edu (Bruce Boghosian)
From: mikk0022@maroon.tc.umn.edu (Christoph L Mikkelson)
From: david@alpha1.csd.uwm.edu (Dave Stack)
From: Mohamed <osman@eecs.wsu.edu>

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: bruceb@conx.bu.edu (Bruce Boghosian)
From: S_DOLLINGE@main01.rz.uni-ulm.de (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.