Continuous Field-(or PDE-)Models¶

Models that acount for change in more than one dimension consist of Partial Differential Equations (PDEs). A PDE is a differential equation that contains functions of multiple variables and their partial derivatives. Partial derivatives simply mean that the function being differentiated has more than one independent variable (e.g., $x, t$ indicating space and time) and that the differentiation is being done while other independent variables are kept as constants. Note that in this case the variable $x$ is not just a state variable of the system, but could represent a position in a continuous space. E.g. it could stand for the Cartesian coordinates $(u, v, w)$.

Here the state of the system is represented by a function, or a field $f(x, t)$, which is defined over a continuous space as well as over time. The value of function $f$ could be scalar $(f(u, v, w)$ or a vector. A scalar field is given if the quality of the represented points is undirected, like temperature, pressure, density (e.g. of a population) etc. A vector field is given if the quality of the represented points is directed, like flows (e.g. of ocean currents), gravitation, magnetic or electrical fields etc.

A corresponding model can include spatial derivatives of the field, which gives information about how the system’s state is shaped spatially.

$\frac{\partial f}{\partial t}=F(f, \frac{\partial f}{\partial x}, \frac{\partial^2 f}{\partial x^2}, ..., x, t)$

An example of a spatial function (a field), which could indicate the density of a population, gives the function $f(x, y,) = e^{-(x^2 + y^2)}$ which is depicted below, including the Python-code that generates these plots.

In [2]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from mpl_toolkits.mplot3d import Axes3D

# preparing a meshgrid as base
X, Y = np.meshgrid(np.arange(-2, 2, 0.02), np.arange(-2, 2, 0.02))
# function for z-values
Z = np.exp(-(X**2 + Y**2))
#Z = X*np.exp(-(X**2 + Y**2))
#Z = np.exp(-((X - 0.5)**2 + (Y - 0.5)**2) / 0.04)

fig = plt.figure(figsize=(15,6))
plt.plot_surface(X, Y, Z, color = 'gold', rstride=10, cstride=10)

plt.plot_wireframe(X, Y, Z, rstride=10, cstride=10)

Out[2]:
<mpl_toolkits.mplot3d.art3d.Line3DCollection at 0xb32c7b8>

Contour¶

A contour is a set of spatial positions $x$ that satisfy $f(x) = C$ for a scalar field $f$, where $C$ is a constant. Contours are also called level curves or surfaces.

The following shows the contour $f(x, y) = \frac{1}{2}$ (blue solid circle) of the above spatial function, again including the Python-code that generates it.

In [11]:
import matplotlib.pyplot as plt
%matplotlib inline
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(15,6))

plt.plot_wireframe(X, Y, Z, rstride=10, cstride=10, colors = 'r')
# this partitions the wireframe into equal parts and draws n contours
plt.contour(X, Y, Z, 1, linewidths=10, colors = 'b')

CS = plt.contour(X, Y, Z, 9,
linewidths=(.5, .5, .5, .5, 4, .5, .5, .5, .5),
colors=('r', 'r', 'r', 'r', 'b', 'r', 'r', 'r', 'r')
)
plt.clabel(CS, inline=1, fontsize=10)

Out[11]:
<a list of 9 text.Text objects>

The gradient of a scalar field $f$ is a vector locally defined as $\nabla f = \vec(\frac{\partial f}{\partial x_1}, \frac{\partial f}{\partial x_2}, ..., \frac{\partial f}{\partial x_n})$ with $n$ indicating the number of dimensions of the space. The symbol $\nabla$ is called "del" or "nabla" and indicates a formal vector consisting of the partial differential operators $\nabla = \vec(\frac{\partial}{\partial x_1}, \frac{\partial}{\partial x_2}, ..., \frac{\partial}{\partial x_n})$.

The gradient indicates the largest change-rate in each point in the scalar-field.

For example, if $f(u, v)$ is a scalar field, then the gradient of $f$ is a vector which points into the direction of the largest change of $f$ at point $P(u, v)$. The length of the vector represents how steep the slope is. Consider the map of a territory which indicates height above sea level at any point $(u, v)$ with the function $f(u, v)$. The gradient of $f$ then is a vector pointing in the direction of the steepest slope at that point. The steepness of the slope at that point is given by the magnitude of the gradient vector.

The gradiant in an $n$-dimensional space consists of the $n$ partial derivatives of the function that determines whatever is depicted in this space. Similarly to the derivative in ODEs, it represents the slope of the tangent of the graph of the function.

The gradiant field of the spatial function $f(x, y)=e^{-(x^2+y^2)}$ looks like the following:

In [13]:
import matplotlib.pyplot as plt
%matplotlib inline
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(16,6))

# using the quiver-plot function of matplotlib needs information on the slope
# this is won with simple finite-difference, as opposed to analytically calculating the
# derivative of the function
dy, dx = np.gradient(Z, np.diff(Y[:2, 0]), np.diff(X[0, :2]))

# consider only every 10th point in each direction
skip = (slice(None, None, 10), slice(None, None, 10))
plt.quiver(X[skip], Y[skip], dx[skip], dy[skip], Z[skip])
plt.set_title('Gradient field of $f(x, y)=e^{-(x^2+y^2)}$')

plt.streamplot(X, Y, dx, dy, color=Z, density=0.5, cmap='gist_earth')
CS = plt.contour(X, Y, Z, cmap='gist_earth')
plt.clabel(CS)

Out[13]:
<a list of 6 text.Text objects>

Divergence¶

The divergence $\nabla \cdot$ of a vector field, $\nabla \cdot v$, is a scalar field indicating the rate at which "density" exits a given region of space. The density within this region can change only by stuff (e.g. particles) flowing into or out of the region. By measuring the net flux of this stuff passing through a surface surrounding the region, it is therefore possible to say how the density of the interior changes. If divergence is positive, it means that stuff is escaping from this local region. If it is negative, the stuff is flowing into the local area. Divergence thus indicates the temporal change of concentration of the stuff in this region. (While air is heated in a region, it expands in all directions, and thus the velocity field points outward from that region. The divergence of the velocity field in that region thus has a positive value. While the air is cooled and thus contracting, the divergence of the velocity has a negative value.)

The individual scalars in divergence (in the scalar field) thus indicate the slope of the differential of the flux at the given points.

The sign for divergence is $\nabla \cdot$. The divergence of a vector field, $\nabla \cdot v$, is seen as a dot product, indicated by the tiny dot between $\nabla$ and $v$.

$$\nabla \cdot = {\frac{\partial}{\partial x}\choose \frac{\partial}{\partial y}}^T$$

Divergence can be used to calculate the transport of stuff in space, for example the transport of a concentration or of a density in space. This is done with

The Transport equation¶

$$\frac{\partial c}{\partial t}=- \nabla \cdot J + s$$

where $c$ is the system’s state defined as a spatio-temporal function that represents the concentration of the stuff moving in space. $J$ is a vector field called the flux of $c$. $\nabla \cdot J$ hence is a scalar field that indicates the rate at which $c$ exits a given region of space.

With $J = -\alpha \nabla c$ we can plug this into $\nabla \cdot J + s$ obtaining $\frac{\partial c}{\partial t}= - \nabla \cdot (-\alpha \nabla c) + s = \nabla \cdot (\alpha \nabla c) + s$. If $\alpha$ doesn't depend on spatial locations, one can take it out of the parenthesis and obtains

$$\frac{\partial c}{\partial t}= \alpha \nabla^2 c + s$$

$\nabla^2$ is called the Lapace-operator, which is discussed below. The term $s$ is often called a source/sink term, which is a scalar field that represents any increase or decrease of $c$ taking place locally - for example if $c$ would also be subject to diffusion.

Simulating the Transport equation (by way of discretizing)¶

The simulation of continuous field- or PDE-models can involve an enormous amount of computation. Simulating global climate systems for instance takes real powerful computers and elaborate techniques for numerical simulation. However, just like differential equations can be converted into difference equations, PDE-models can be converted into Cellulare Automata-like models by discretizing space and time. Be aware though that similar accuracy-problems may arise if differences are taken too small.

Discretizing time:¶

Like in the case of ODEs, $\frac{\partial f}{\partial t}=F(f, \frac{\partial f}{\partial x}, \frac{\partial^2 f}{\partial x^2}, ..., x, t)$ becomes $\frac{\Delta f}{\Delta t}= \frac{f(x, t+\Delta t)-f(x,t)}{\Delta t}$

(Note the difference between $\Delta$ and $\nabla$)

Discretizing Space¶

Here, one has to consider that the spatial dimension - other than time - should be modeled symmetrically. So the actual calculation point should not be in the beginning, but in the middle of the difference step considered. Therfore the spatial derivatives of $\frac{\partial f}{\partial x}$ become $\frac{2\Delta f}{2\Delta x}=\frac{f(x+\Delta x,t)-f(x-\Delta x, t)}{2 \Delta x}$

With this, the transport of a concentration of stuff in space and in time can be modelled like this (in a 2-D space without source or sink):

$$\frac{\partial c}{\partial t}=- \nabla \cdot (c w)$$

with $w$ being a constant vector indicating the velocity of particles $w = {w_x \choose w_y}$

This can be discretized in time $(\Delta t)$ and space $(\Delta s)$ as follows

$$c(x,y,t+\Delta t) \approx c(x,y,t) - \nabla \cdot (c(x,y,t)w)\Delta t$$$$= c(x,y,t) - {\frac{\partial}{\partial x}\choose \frac{\partial}{\partial y}}^T \left( c(x,y,t){w_x \choose w_y} \right) \Delta t$$$$= c(x,y,t) - (w_x c(x+\Delta s, y, t) - w_x c(x - \Delta s, y, t)+ w_y c(x,y + \Delta s, t)-w_y c (x,y-\Delta s, t))\frac{\Delta t}{2\Delta s}$$

This gives the base for the numerical simulation, as used in the following code

In [14]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from mpl_toolkits.mplot3d import Axes3D

n = 100 # size of grid: n * n
ds = 1. / n # spatial resolution, assuming space is [0,1] * [0,1]
dt = 0.01 # temporal resolution

# preparing a meshgrid as base
x, y = np.meshgrid(np.arange(0, 1, ds), np.arange(0, 1, ds))
global z, nextz

# function for z-values
z = np.exp(-((x - 0.5)**2 + (y - 0.5)**2) / 0.04)
# prepare empty container for next step
nextz = np.zeros([n, n])

# plot initial state
fig = plt.figure(figsize=(20,6))
plt.plot_surface(x, y, z, color = 'gold', rstride = 5, cstride = 5)
plt.set_title('time = 0')

# constant velocity of movement, here only into y-direction
wx, wy = 0., 0.03

def move():
global z, nextz
for u in xrange(n):
for v in xrange(n):
# the state-transition function as described in the explanation above
nextz[u, v] = z[u, v] - ( wx * z[(u + 1)%n, v] - wx * z[(u - 1)%n, v]
+ wy * z[u, (v + 1)%n] - wy * z[u, (v - 1)%n]) * dt/(2*ds)
z, nextz = nextz, z

c = 2
for i in range(1, 1001):
# draw only every 500th result
if i%500 == 0:
plt.plot_surface(x, y, nextz, color = 'gold', rstride = 5, cstride = 5)
plt.set_title('time = %0.f' %i)
plt.set_zlim(0,1)
c += 1
move()


Divergence is also needed to comprehend the Laplace operator.

The Laplace operator¶

The Laplace operator (or Laplacian) is a differential operator given by the divergence of the gradient of a function. It is denoted as $\nabla \cdot \nabla$ or $\nabla^2$ (in mathematics sometimes also $\Delta$).

$$\nabla^2 = \sum_{i=1}^n \frac{\partial^2}{\partial x_i^2}$$

in $n$ dimensions.

The Laplace operator represents the flux density of the gradient flow of a function. For instance, the net rate at which a chemical dissolved in a fluid moves toward or away from some point is proportional to the Laplace operator of the chemical concentration at that point. The resulting equation is the diffusion equation (see below).

The Laplace operator is a second order differential operator in the n-dimensional Euclidean space, defined as the divergence $\nabla \cdot$ of the gradient. Thus if $ƒ$ is a twice-differentiable real-valued function, then the Laplace operator of $ƒ$ is defined by $\nabla \cdot \nabla$ or $\nabla^2$, with the partial differential operator $\nabla = \vec(\frac{\partial}{\partial x_1}, \frac{\partial}{\partial x_2}, ..., \frac{\partial}{\partial x_n})$.

The Laplace operator can be used to fomulate

The Diffusion equation¶

as derived from the Transport equation

$\frac{\partial c}{\partial t}= \alpha \nabla^2 c$ , also called the heat equation or Fick’s second law of diffusion. $\alpha$ is called the diffusion constant, which determines how fast the diffusion takes place.

Expressing the Laplacian, the diffusion equation (in 2 dimensions) can be written $\frac{\partial c}{\partial t} = \alpha (\frac{\partial^2 c}{\partial x^2} + \frac{\partial^2 c}{\partial y^2})$

Simulating the Diffusion equation (by way of discretizing)¶

Space and time discretizing is done in the same way as above. The state transition however, here uses the abbreviations $zC, zR, zL, zU, zD$ for the five cells in a Von-Neumann-neighborhood (that is, the C = center cell, the R = right cell, the L = left cell, the U = up cell and the D = down cell, for the cells to the East, the West, the North and the South of a given cell in two dimensions), corresponding here to $(u,v), (u+\Delta s,v), (u-\Delta s,v),(u, v+\Delta s)$ and $(u, v-\Delta s)$.

The Laplacian then is obtained by summing the values of $f$ in all nearest neighbor cells (that is, the left, right, top and bottom cell in a 2-D space) and then subtracting the value of the center cell $f(x; t)$ as many times as the number of neighbors (i.e. here $zLap = zR + zL + zU + zD - 4zC$).

In [117]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from mpl_toolkits.mplot3d import Axes3D

n = 100 # size of grid: n * n
ds = 1. / n # spatial resolution, assuming space is [0,1] * [0,1]
dt = 0.02 # temporal resolution

alpha = 0.001 # diffusion constant of c

# preparing a meshgrid as base
x, y = np.meshgrid(np.arange(0, 1, ds), np.arange(0, 1, ds))
global z, nextz

# function for z-values
z = np.exp(-((x - 0.5)**2 + (y - 0.5)**2) / 0.04)
# prepare empty container for next step
nextz = np.zeros([n, n])

# plot initial state
fig = plt.figure(figsize=(20,6))
plt.plot_surface(x, y, z, color = 'gold', rstride = 5, cstride = 5)
plt.set_title('time = 0')

def diffuse():
global z, nextz
for u in xrange(n):
for v in xrange(n):
# state-transition function
# center, right, left, up, down
# %n as a simple way to consider periodic boundery conditions
zC, zR, zL, zU, zD = z[u,v], z[(u+1)%n,v], z[(u-1)%n,v], z[u,(v+1)%n], z[u,(v-1)%n]
# Laplacian - Just add the values of f in all of your nearest neighbors
# (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.
zLap = zR + zL + zU + zD - 4 * zC
# the actual diffusion equation (discretized)
nextz[u,v] = zC + (alpha * zLap) * dt / (ds**2)
z, nextz = nextz, z

r = 2
for i in range(1, 1001):
# draw only every 500th result
if i%500 == 0: