How to solve (O/P)DE via cellular automaton?

Question

I heard that one can solve differential equations by a cellular automaton. Even though I have searched some papers, which seem too formal to me, I do not get the hang of it. Anyone knows how to get started?

For example, a simple initial value problem: $ y'(x)= y(x), y(0)=1, x\in[0,1] $, how to use the Mathematica built-in cellular automaton capabitlies to find its numerical solution?

Answer 1

Here's an example:

dx = 1./1000;
grid = ConstantArray[0, 1000];
grid[[1]] = 1;

sol = CellularAutomaton[{
     {0., x_, _} :> x,
     {x_, _, _} :> x + dx x
     }, grid, 1000][[-1, 2 ;;]];

sol2 = NDSolveValue[{
    y'[x] == y[x],
    y[0] == 1
    }, y, {x, 0, 1}];

Show[
 ListLinePlot[sol, DataRange -> {0, 1}, PlotStyle -> Thickness[0.02]],
 Plot[sol2[x], {x, 0, 1}, PlotStyle -> {Dashed, White}]
 ]

Solving differential equations with cellular automata is just like simple Euler integration. It can be done for many types of differential equations by discretizing them using e.g. central differences. In this case I'm just applying the simple rule that $$ y(x+\mathrm dx) = y(x) + y(x)\mathrm dx $$ In Mathematica, this rule looks like this:

{x_, _, _} :> x + dx x

The left-hand side in this rule is centered on the element that will be replaced by the right-hand side.

CellularAutomaton assumes the grid to be cyclical, at least I couldn't find a way to turn that off. (Of course, there are other ways of simulating cellular automata in Mathematica, but I'm sticking to this to emphasize that it is a cellular automaton I'm solving the equation with.) This means that the first element in the list, set to 1 because this is the initial condition, will be replaced by zero because the preceding list element is zero. To stop this, I included the rule

{0., x_, _} :> x

The 0 in this rule is the last element in the grid.

Maybe it is arguable if this is truly a cellular automaton because the states are not discrete, but I've also read about "real-valued cellular automata" and "continuous automata". When one talks about solving PDEs using cellular automata, I'm pretty sure that one is talking about a broader class of lattice-based methods.

Next, I will show a problem from an old homework of mine regarding a reaction-diffusion system. This solves a coupled system of differential equations on a grid using methods similar to cellular automata (but definitely using continuous states.) I think that when you heard about PDEs and cellular automata, it might have been something like this that your source was thinking about.

The differential equations are $$ \frac{\partial u}{\partial t} = a - (b+1)u+u^2v+\nabla^2u $$ $$ \frac{\partial v}{\partial t} = bu - u^2v+D\nabla^2 v $$

You can think of $u$ and $v$ as the number of particles of two different types. The particles spread out (diffuse) in space as well as react with each other, which may cause a particle of one type to convert to the other type. One might imagine that this would lead to very simple solutions; that all the particles will eventually be of one type, or that all particles will be spread out uniformly throughout space. It turns out, however, that the solutions can be much more complicated than that.

In order to simulate the system as cellular automata, we first have to discretize the operators in the equations. In this case, we only have the Laplacian. The following can be shown to be a discrete approximation: $$ \Delta f(x,y) \approx \frac{f(x-h, y) + f(x+h,y)+f(x,y-h)+f(x,y+h)-4f(x,y)}{h^2} $$

We could implement this as a convolution with a kernel like this:

ker = {
   {0, 1, 0},
   {1, -4, 1},
   {0, 1, 0}
   };
laplacian = ListConvolve[ker, grid];

But we could also use the built-in Laplacian function:

LaplacianFilter[u, 1, Padding-> "Periodic"]

I specify that the grid should be periodic because that is what we want for this particular PDE.

Doing a (discrete) convolution and evolving a cellular automata is the same thing. They both replace each element in the grid with a new value depending on the neighbors of that grid.

Now on to solving the equations. Since we have two variables, we are going to have to use two lattices, one for each variable. We will evolve them alternately and use them as inputs to each other.

Here it is:

dudt[a_, b_, u_, v_] := a - (b + 1) u + u^2 v + LaplacianFilter[u, 1, Padding -> "Periodic"]
dvdt[b_, u_, v_, d_] := b  u - u^2 v + d LaplacianFilter[v, 1, Padding -> "Periodic"]

step[a_, b_, d_][{uvals_, vvals_}] := {
  uvals + 0.01 dudt[a, b, uvals, vvals],
  vvals + 0.01 dvdt[b, uvals, vvals, d]
  }

simulate[d_, nrOfIterations_] := Module[
  {u, v, a = 3, b = 8, L = 128},
  u = ConstantArray[a, {L, L}] + RandomReal[{-0.1 a, 0.1 a}, {L, L}];
  v = ConstantArray[b/1, {L, L}] + RandomReal[{-0.1 b/a, 0.1 b/a}, {L, L}];
  Nest[step[a, b, d], {u, v}, nrOfIterations]
  ]

res = simulate[2.3, 20000];
MatrixPlot[#, ImageSize -> 250, PlotTheme -> "Monochrome"] & /@ res // Row

The left plot shows the solution for u and the right plot the solution for v. The solutions to the differential equations were found by updating cells in lattices by local rules (only using the cell's value and its four direct neighbors to do the updating.) As I've said, I think this is what people may refer to as solving PDEs with cellular automata, because there are a lot of similarities. A simpler example with only one variable differs less still, mainly in the fact that these partial differential equations are continuous.

The solutions u and v may denote the concentration of two different types of particles, and the interesting thing here is that the distributions of the particles are not uniform. The distributions show quite complicated patterns, in fact, and this type of pattern formation is thought to explain a lot of the complicated patterns that we see everywhere in nature. The property that local rules can create complicated patterns is also what cellular automata with discrete states are known for, so also in this regard there is a relationship.