# 13.5: Simulation of Continuous Field Models

\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)

\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)

Simulation of continuous ﬁeld models written in PDEs is not an easy task, because it easily involves an enormous amount of computation if you want to obtain fairly accurate simulation results, and also because certain kinds of spatial dynamics may sometimes cause highly sensitive behaviors. In fact, most “supercomputers” built for scientiﬁc computation^{2} are used to solve large-scale PDE problems, e.g., to simulate global climate systems, complex material properties, etc.

Covering all the advanced techniques of numerical simulations of PDEs is way beyond the scope of this introductory textbook. Instead, we will try a much simpler approach. As we already discussed cellular automata (CA), we will convert PDE-based models into CA by discretizing space and time, and then simulate them simply as CA models with continuous state variables.

Discretizing time in a PDE is nothing different from what we discussed in Section 6.3. By replacing the time derivative by its discretized analog, the original PDE

\[\dfrac{\partial{f}}{\partial{t}}= F(f, \dfrac{\partial{f}}{\partial{x}}, \dfrac{\partial^{2}{f}}{\partial{x^{2}}} , \dots, x, t)\label{(13.29)}\]

becomes

\[\dfrac{\Delta{f}}{\Delta{t}} =\dfrac{f(x, t +\Delta{t}) -f(x, t)}{\Delta{t}}=F(f, \dfrac{\partial{f}}{\partial{x}}, \dfrac{\partial^{2}{f}}{\partial{x^{2}}} , \dots, x, t), \label{(13.30)}\]

\[f(x, t+ \Delta{t}) =f(x,t) +F(f, \dfrac{\partial{f}}{\partial{x}}, \dfrac{\partial^{2}{f}}{\partial{x^{2}} }, \dots, x, t)\Delta{t} \label{(13.31)}\]

which is now a difference equation over time (note that ∆ above is not a Laplacian!).

We are not done yet, because the right hand side may still contain spatial derivatives (\(∂f/∂x\), \(∂2f/∂x2\), ...). We need to discretize them over space too, in order to develop a CA version of this PDE. Here is one way to discretize a ﬁrst-order spatial derivative:

\[\dfrac{\partial{f}}{\partial{x}} \approx \dfrac{\Delta{f}}{\Delta{x}} =\dfrac{f(x, t +\Delta{t}) -f(x, t)}{\Delta{x}} \label{(13.32)}\]

This is not incorrect mathematically, but there is one practical issue. This discretization makes the space asymmetric along the x-axis, because the spatial derivative at \(x\) depends on \(f(x + ∆x,t)\), but not on \(f(x−∆x,t)\). Unlike time that has a clear asymmetric direction from past to future, the spatial axis should better be modeled symmetrically between the left and the right (or up and down). Therefore, a better way of performing the spatial discretization is as follows:

\[\dfrac{\partial{f}}{\partial{x}} \approx \dfrac{2\Delta{f}}{2\Delta{x}}= \dfrac{f(x, t +\Delta{t}) -f(x, t)}{2\Delta{x}}\label{(13.33)}\]

This version treats the left and the right symmetrically.

Similarly, we can derive the discretized version of a second-order spatial derivative, as follows:

\[ \dfrac{\partial^{2}{f}}{\partial{x^{2}}} \approx \dfrac{\dfrac{\partial{f}}{\partial{x}}|_{x+\Delta{x}} -\dfrac{\partial{f}}{\partial{x}}|_{x-\Delta{x}}}{2\Delta{x}}\]

\[f(x + ∆x + ∆x,t ) \qquad f(x−∆x + ∆x,t) \]

\[ \approx \dfrac{\dfrac{-f(x +\Delta{x} -\Delta{x,t})}{ 2\Delta{x}} - \dfrac{-f(x-\Delta{x} -\Delta{x,t})}{2\Delta{x}}} {2\Delta{x}} \]

\[ =\dfrac{f(x+2\Delta{x,t}) +f(x-2\Delta{x,t}) -2f(x,t)}{2\Delta{x^{2}}} \]

\[=\dfrac{f(x+\Delta{x, t}) +f(x-\Delta{x, t}) -2f(x,t)}{\Delta{x^{2}}} \]

\[ =\dfrac{f(x +\Delta{x, t}) +f(x-\Delta{x, t}) -2f(x,t)}{\Delta{x^{2}}} \text{by replacing 2∆x → ∆x } \label{(13.34)}\]

Moreover, we can use this result to discretize a Laplacian as follows:

\[\nabla^{2} f =\dfrac{\partial^{2} {f}}{\partial{x}^{2}_{1}} + \dfrac{\partial^{2} {f}}{\partial{x}^{2}_{2}} +.....+ \dfrac{\partial^{2} {f}}{\partial{x}^{2}_{n}} \]

\[\approx \dfrac{f(x_1 + ∆x_1,x_2,...,x_n,t) + f(x_1 −∆x_1,x_2,...,x_n,t)−2f(x_1,x_2,...,x_n,t) }{∆x^2_1 } \]

\[ + \dfrac{f(∆x_1,x_2,...,x_n,t) + f(∆x_1,x_2 -\Delta{x_2},...,x_n,t)−2f(x_1,x_2,...,x_n,t) }{∆x^2 _2} \]

\[+\dots\]

\[ + f(x_1,x_2,...x_n +\Delta{x_n},t) + f(x_1,x_2,...,x_n -\Delta{x_n},t)\]

\[= ( f(x_1+ \Delta{x}, x_{2},\dots, x_{n}, t ) +f(x_1 -\Delta{x}, x_2, ...x_n, t) -2f(x_{1}, x_{2},...., x_{n}, t) \]

\[+ f(x_1, x_2, +\Delta{x}, ...., x_n , t)+f(x_1, x_2, -\Delta{x}, ...,x_n, t)\]

\[+....\]

\[ + f(x_1, x_2, ...., x_n +\Delta{x, t}) +f(x_1, x_2,..., x_n -\Delta{x, t})\]

\[-2nf(x_1, x_2,...,x_n, t)/ {\Delta{x^{2}} } \label{(13.35)} \]

\[ \text {(by replacing ∆xi → ∆x for all i)}\]

The numerator of the formula above has an intuitive interpretation: “Just add the values of \(f\) in all of your nearest neighbors (e.g., top, bottom, left, and right for a 2-D space) and then subtract your own value \(f(x,t)\) as many times as the number of those neighbors.” Or, equivalently: “Measure the difference between your neighbor and yourself, and then sum up all those differences.” It is interesting to see that a higher-order operator like Laplacians can be approximated by such a simple set of arithmetic operations once the space is discretized.

You can also derive higher-order spatial derivatives in a similar manner, but I don’t think we need to cover those cases here. Anyway, now we have a basic set of mathematical tools (Eqs. \ref{(13.31)}, \ref{(13.33)}, \ref{(13.34)}, \ref{(13.35)}) to discretize PDE-based models so we can simulate their behavior as CA models.

Let’s work on an example. Consider discretizing a simple transport equation without source or sink in a 2-D space, where the velocity of particles is given by a constant vector:

\[ \dfrac{\partial{c}}{\partial{t}} =-\nabla \cdot(cw) \label{(13.36)} \]

\[ w = \binom{w_x}{w_y} \label{(13.37)} \]

We can discretize this equation in both time and space, as follows:

\[c(x, y, t +\delta{t}) \approx c(x, y,t) -\nabla \cdot (c(x,y,t)w)\delta{t} \label{(13.38)} \]

\[=c(x, y, t ) - {\begin {pmatrix} \dfrac{\partial}{\partial{x}} \\ \dfrac{\partial}{\partial{y}} \end{pmatrix}}^{T} (c(x,y,t) \begin{pmatrix} w_x \\ w_y \end{pmatrix})\Delta{t} \label{(13.39)}\]

\[=c(x, y, t ) - (w_x \dfrac{\partial{c}}{\partial{x}} + w_y \dfrac{\partial{c}}{\partial{y}}) \Delta{ t} \label{(13.40)}\]

\[\approx c(x, y, t ) - (w_x \dfrac{c(x +\Delta{h, y, t}) - c(x-\Delta{h,y, t})}{2\Delta{h}} + w_y \dfrac{c(x,y +\Delta{h, t}) - c(x, y-\Delta{h, t})}{2\Delta{h}}) \Delta{t} \label{(13.41)} \]

\[= c(x, y, t ) -(w_{x}c(x +\Delta{h, y, t}) - w_{x}c(x -\Delta{h, y, t}) + w_{y}c(x,y +\Delta{h, t}) - w_{y}c(x, y-\Delta{h, t})) \dfrac{\Delta{t}}{2\Delta{h}} \label{(13.42)} \]

This is now fully discretized for both time and space, and thus it can be used as a state-transition function of CA. We can implement this CA model in Python, using Code 11.5 as a template. Here is an example of implementation:

In this example, we simulate the transport equation with \((w_x,w_y) = (0.01,0.03)\) in a 2-D `[0,1]×[0,1]`

space with periodic boundary conditions, starting from an initial conﬁguration that has a peak in the middle of the space:

\[c(x,y, 0) =exp(-\dfrac{(x-0.5)^{2} +(y+ 0.5)^{2}}{0.2^{2}}) \label{(13.43)}\]

The temporal and spatial resolutions are set to \(∆t = 0.01\) and \(∆h = 0.01\), respectively (so that the size of the grid is `100×100`

). This is a very slow simulation process, so the` stepSize `

is set to 50 to skip unnecessary visualizations (see the last line). The surface plot is used to visualize the conﬁguration in 3-D axes. Note that you need to call the `show`

command at the end of the `observe `

function, because changes made to 3-D axes are not automatically reﬂected to the visualization (at least in the current implementation of matplotlib).

The result of this simulation is shown in Fig.\ref{13.13}, where you can clearly see the peak being transported to a direction given by \((w_x,w_y)\). Moreover, thanks to the interactive nature of pycxsimulator.py, this 3-D visualization is now interactively manipulatable even if you are not running it from Anaconda Spyder, Enthought Canopy, or other interactive environments. You can click and drag on it to rotate the surface to observe its 3-D structure from various angles.

If you don’t need a 3-D surface plot, yo ucan use` imshow`

instead. This makes the code a bit simpler, and yet the result can be as good as before (Fig. 13.14):

Figure 13.13: Visual output of Code 13.5. Time ﬂows from left to right.

One thing I should warn you about is that the simulation method discussed above is still based on the Euler forward method, so it easily accumulates numerical errors. This is the same issue of the stability of solutions and the possibility of “artifacts” that we discussed in Section 6.4, now arising from discretization of both space and time. In the example of the transport equation above, such issues can occur if you choose a value that is too large for \(∆t\) or \(∆h\), relative to the spatial/temporal scale of the system’s behavior (which is determined by w in this case). For example, if you increase \(∆t\) to 0.1, you will get the

Figure 13.14: Visual output of Code 13.6. Time ﬂows from left to right.

result shown in Fig. 13.15. This may look cool, but it is actually an invalid result. You can tell it is invalid if you understand the nature of the transport equation; it should only transport the surface and not change its shape. But even if you are unsure if the result you obtained is valid or not, you can try increasing the temporal/spatial resolutions (i.e., making \(∆t\) and \(∆h\) smaller) to see if the result varies. If the result remains unaffected by this, then you know that the original spatial/temporal resolutions were already sufﬁcient and that the original result was thus valid. But if not, the original result probably contained numerical errors, so you need to keep increasing the resolutions until the result becomes consistent.

Figure 13.15: Accumulation of numerical errors caused by increasing \(∆t (Dt)\) to 0.1 in Code 13.5. Time ﬂows from left to right.

Before moving on to the next topic, I would like to point out that the discretization/simulation methods discussed above can be extended to PDEs that involve multiple state variables. Here is such an example: interactions between two spatially distributed populations, a escaping from diffusing \(b\):

\[\dfrac{\partial{a}}{\partial{t}} =- \mu_{a} \nabla \cdot(a(-\nabla{b})) \label{(13.44)}\]

\[ \dfrac{\partial{b}}{\partial{t}} =\mu_{b} \nabla^{2} b \label{(13.45)}\]

Here \(µ_a\) and \(µ_b\) are the parameters that determine the mobility of the two species. These equations can be expanded and then discretized for a 2-D space, as follows (you should also try doing this all by yourself!):

\[\dfrac{\partial{a}}{\partial{t}} =-\mu_{a} (\begin{pmatrix} \dfrac{\partial}{\partial{x}} \\ \dfrac{\partial}{\partial{y}}) \end{pmatrix}^{T} \begin{pmatrix} -a\dfrac{\partial{b}}{\partial{x}} \\ -a\dfrac{∂b}{∂y}\end{pmatrix} \label{(13.46)}\]

\[=\mu_{a} (\dfrac{\partial}{\partial{x}} (a\dfrac{∂b}{∂x}) + \dfrac{\partial}{\partial{y}}(a\dfrac{∂b}{∂y})) \label{(13.47)}\]

\[= \mu_{a} ( \dfrac{\partial{a}}{\partial{x}} \dfrac{∂b}{∂x} + a\dfrac{∂^{2}b}{∂x^{2}} + \dfrac{\partial{a}}{\partial{y}} \dfrac{∂b}{∂y} +a\dfrac{∂^{2}b}{∂y^{2}} \label{(13.48)}\]

\[=\mu_{a} (\dfrac{\partial{a}}{\partial{x}} \dfrac{∂b}{∂x} +\dfrac{\partial{a}}{\partial{y}} \dfrac{∂b}{∂y} a\nabla^{2}b) \\label{(13.49)} \]

\[ a(x,y, t +\Delta{t}) \approx a(x, y, t) + µ_a(\dfrac{a(x + ∆h,y,t)−a(x−∆h,y,t)}{ 2∆h}

\dfrac{b(x + ∆h,y,t)−b(x−∆h,y,t)}{ 2∆h} + \dfrac{a(x,y + ∆h,t)−a(x,y−∆h,t)}{ 2∆h}

\dfrac{b(x,y + ∆h,t)−b(x,y−∆h,t)}{ 2∆h}

+ a(x,y,t)× \dfrac{b(x + ∆h,y,t) + b(x−∆h,y,t) + b(x,y + ∆h,t) + b(x,y−∆h,t)−4b(x,y,t)}{\Delta{h^{2}}})\Delta{t} \label{(13.50)}\]

\[a'c \approx a_c +\mu_{a} (\dfrac{a_R −a_L}{2∆h} \dfrac{b_R −b_L}{2∆h} + \dfrac{a_V −a_D}{2∆h} \dfrac{b_V −b_D}{2∆h} + a_C \dfrac{b_R + b_L + b_U + b_D −4b_C} {∆h^{2}} )∆t \label{(13.51)} \]

Note that I used simpler notations in the last equation, where subscripts \((C, R, L, U, D)\) represent states of the central cell as well as its four neighbors, and \(a'C\) is the next value of \(a_C\). In the meantime, the discretized equation for b is simply given by

\[b'_C \approx b_C +\mu_b \dfrac{b_R + b_L + b_U + b_D −4b_C}{\Delta{h^{2}}} \Delta{t}. \label{(13.52)}\]

Below is a sample code for simulating this PDE model starting with an initial conﬁguration with two peaks, one for a and the other for \(b\):

Note that the `observe`

function now uses` subplot`

to show two scalar ﬁelds for \(a\) and \(b\). The result is shown in Fig. 13.16.

Figure 13.16: Visual output of Code 13.7. Time ﬂows from left to right.

Exercise 13.15

Modify Code 13.7 so that both species diffuse and try to escape from each other. Run simulations to see how the system behaves.

__ __

^{2}If you are not clear on what I am talking about here, see `http://top500.org/`

.