PDEs

The PDE of diffusion

Just as advection, heat transfer exhibits a familiar dispersal pattern. However, without a flow to carry material around, the pattern lacks a direction:

Instead of being carried around, the heat is exchanged -- with adjacent locations.

The process is also known as diffusion.

The simplest way to model heat propagation is to set up a grid of square rooms and then increment the temperature of each room by a proportion of the average of the temperature of the four adjacent rooms. It is given by the following short formula: $$RC= RC -.25 *h*\bigg( (RC-RC[-1]) + (RC-RC[1]) + (RC-R[-1]C) + (RC-R[1]C) \bigg),$$ Only a proportion $h$, dependent on the presumed length of the time interval, of the total amount is shared.

Below are the setup and the results (after $1,700$ iterations) of such a simulation for a single initial hot spot, $h=.01$, and the temperature at the border fixed at $0$:

On a larger scale, the simulation produces a realistic circular pattern:

A more careful look reveals that to model heat transfer, we need to separately record the exchange of heat with each of the four adjacent rooms.

The process we are to follow is the Newton's Law of Cooling:

• the rate of cooling of an object is proportional to the difference between its temperature and the ambient temperature.

This law is nothing but a version of the ODE of population growth and decay -- with respect to the exterior derivative $d_t$ over time. This time, however, for each cell there are four adjacent cells and four temperature differences, $d_x$, to be taken into account. The result is a partial differential equation (PDE).

Below, we derive a PDE of forms for the diffusion process. The set-up is as follows. A certain material is contained in a grid of rooms and each room (an $n$-cell) exchanges the material with its neighbors through its walls ($(n-1)$-cells). We will assume initially that the space is the cubical complex ${\mathbb R}^n$. The time is the standard complex ${\mathbb R}$. For now, we ignore the geometry of time and space.

Dually, one can think of pipes ($1$-cells) connecting the compartments or ponds ($0$-cells) in the centers of the rooms each of which exchanges the material with its neighbors. We will use both of the two approaches.

Below, the rooms are blue, the walls are green, and the pipes are gray:

We assume that the time increment is long enough

• for the material in each room to spread uniformly; or
• for the material to pass from one end of the pipe to the other.

What makes this different from ODEs is that the forms will have two degrees -- with respect to location $x$ and with respect to time $t$.

The amount of material $U=U(\alpha,t)$ is simply a number assigned to each room $\alpha$ which makes it an $n$-form. It also depends on time which makes it a $0$-form. We call it an $(n,0)$-form, an $n$-form with respect to location and a $0$-form with respect to time.

Definition. In the product of two complexes $K\times L$, an $(k,m)$-form is a real-valued linear operator defined on cells $a\times b$, where $a$ is a $k$-cell in $K$ and $b$ is an $m$-cell in $L$.

The outflow gives the amount of flow across an $(n-1)$-face (from the room to its neighbor) per unit of time. It is also a number assigned to each “pipe” $p$ that goes through each wall from one cell to the next. So, $F=F(p,t)$ is a dual $(n-1,0)$-form.

Notation.

• $d_t$ is the exterior derivative with respect to time (just the difference in ${\mathbb R}$); and
• $d_x$ is the exterior derivative with respect to location.

Below, we routinely suppress $t$ as the second variable.

The “conservation of matter” in cell $\alpha$ gives us the following. The change of the amount of the material in room $\alpha$ over the increment of time is equal to $$d_t U(\alpha)=−\bigg( \text{ sum of the outflow } F(·) \text{ across the walls of } \alpha\bigg).$$ Of course, the walls form the boundary $\partial \alpha$ of $\alpha$. Therefore, we have $$d_t U(\alpha)=- F^\star(\partial \alpha)$$ or, by the Stokes Theorem, $$d_t U(\alpha) = −(d_x F^\star)(\alpha).$$

Now, we need to express $F$ in terms of $U$. The flow $F(a^\star)=F^\star(a)$ through wall $a$ of room $\alpha$ is proportional to the difference of the amounts of material (or, more precisely, the density) in $\alpha$ and the other room adjacent to $a$. So, $$F^\star(a) = - k(a)d_x(\star U)(a^\star).$$ Here, $k(a) \ge 0$ represents the permeability of the wall $a$ at a given time. Since $a$ is an $(n-1)$-cell, $k$ is an $(n-1)$-form with respect to space. It is also a $0$-form with respect to time.

The result of the substitution is a PDE of second degree called the diffusion equation of forms: $$\begin{array}{|c|} \hline \\ \quad d_t U = d_x \star kd_x \star U \quad \\ \\ \hline \end{array}$$

The right-hand side is seen in the Hodge duality diagram below, if we start with $U$ in the right upper corner and make the full circle: $$ \newcommand{\ra}[1]{\!\!\!\!\!\!\!\xrightarrow{\quad#1\quad}\!\!\!\!\!} \newcommand{\da}[1]{\left\downarrow{\scriptstyle#1}\vphantom{\displaystyle\int_0^1}\right.} \newcommand{\la}[1]{\!\!\!\!\!\!\!\xleftarrow{\quad#1\quad}\!\!\!\!\!} \newcommand{\ua}[1]{\left\uparrow{\scriptstyle#1}\vphantom{\displaystyle\int_0^1}\right.} % \begin{array}{cccl} C ^{n-1}(K) & \ra{d_x} & C ^{n}(K) & \ni U \\ \da{\star}\ua{} & \ne & \da{\star}\ua{} \\ C ^{1}(K^\star)& \la{d_x} & C ^{0}(K^\star).\\ \end{array} $$ The multiplication by $k$ is implied.

In general, $k$ cannot be factored out. However, if $k$ happens to be constant with respect to space, the right-hand side is $kU' '$.

This was an outline. In the following, we develop both the mathematics and the simulations for progressively more complex settings.

Diffusion with a spreadsheet

For the simplest Excel simulation, let's “trivialize” the above analysis. We start with dimension $1$.

The material is contained in a row of rooms and each room exchanges the material with its two neighbors through its walls. The amount of material is a $(1,0)$-form.

The outflow, the amount of flow across a wall (from the room to its neighbor) per unit of time, $F=F(p,t)$ is a $(1,1)$-form over the dual complex.

Suppose

• $a=AB$ is one of the rooms;
• $b,c$ are the two adjacent rooms, left and right;
• $A,B$ are the walls of $a$, left and right;
• $p=A^\star,q=B^\star$ are the two pipes from $a$, left and right.

As the material is conserved, we know that $d_t U(a)$ is the negative of the sum of the outflow $F(·)$ across the two walls of $a$. The flow is positive at $A$ if it is from left to right and the opposite for $B$; then: $$d_t U(a)=-\left( F^\star(A)-F^\star(B) \right) = F^\star(A)-F^\star(B),$$ which is the exterior derivative of a $0$-form. The flow through a wall is proportional to the difference of amounts of material in the two adjacent rooms: $$F^\star(A) = - k(A) \left( U(a) - U(b) \right),$$ $$F^\star(B) = - k(B) \left( U(c) - U(a) \right),$$ where $k\ge 0$ is the permeability, a $(0,0)$-form. Then: $$d_t U = \Big(- k(A) \left( U(a) - U(b)\right)\Big) - \Big(- k(B) \left( U(c) - U(a)\right)\Big). $$ The right-hand side becomes the increment in the recursive formula for the simulation: $$U(a, t+1):= U(a,t) + \Big[- k(a) \left( U(a) - U(a-1)\right) + k(a+1) \left( U(a+1) - U(a)\right)\Big].$$ The Excel setup is shown below:

The result after $1,500$ iterations is shown next:

One can clearly see how the density becomes uniform eventually -- but not across an impenetrable wall ($k(A_5)=0$).

Note: When the domain isn't the whole space, the pipes at the border of the region have to be removed. In the spreadsheet, we use boundary conditions to substitute for the missing data.

Link to file: Spreadsheets.

Exercise. For an infinite sequence of rooms, what is the limit of $U$ as $t\to \infty$?

Exercise. Create a simulation for a circular sequence of rooms. What is the limit state?

Exercise. Generalize the formula to the case when the permeability of walls depends on the direction.

Next, dimension $2$.

The amount of material $U=U(\tau,t)$ is a $2$-form with respect to location and a $0$-form with respect to time.

The outflow gives the amount of flow across a $1$-face (from the room to its neighbor) per unit of time. It is simply a number assigned to each “pipe” $p$ that goes through each wall from one cell to the next. So, $F=F(p,t)$ is a $(1,1)$-form, but over the dual grid.

The “conservation of the matter” for cell $\tau$ implies that the change of the material over the increment of time $d_t U(\tau)$ is, as before, the negative sum of the outflow $F(·)$ across the four walls of $\tau$: $a,b,c,d$. Let's rewrite the latter: $$=-\left( F^\star(a)+F^\star(b)+F^\star(c)+F^\star(d)\right) =- F^\star(a+b+c+d).$$

Now, we need to express $F$ in terms of $U$. The flow $F(a^\star)=F^\star(a)$ through face $a$ of cell $\tau$ is proportional to the difference of the amounts of material in $\tau$ and the other $2$-cell adjacent to $a$. So, $$F^\star(a) = - kd_x(\star U)(a^\star).$$

Exercise. Demonstrate how the above analysis leads to the “naive” Excel simulation presented in the beginning.

Note: One needs to use a buffer (an extra sheet); otherwise Excel's sequential -- cell-by-cell in each row and then row-by-row -- manner of evaluation will cause a skewed pattern:

A simulation of heat transfer from single point is shown below:

Link to file: Spreadsheets.

Exercise. (a) What kind of medium would create an oval diffusion pattern? Demonstrate. (b) What about an oval not aligned with the axes?

Exercise. Modify the formula to create a simulation for diffusion on the torus.

Diffusion on geometric complexes

We have used interchangeably “amount of material” and “density of material”. The reason is that we ignore the possibility of cells of different sizes (and shapes) by assuming the simplest geometry.

We have modeled heat transfer with: $$RC= RC - k*\bigg( \left(RC-RC[-1]\right) + \left(RC-RC[1]\right) + \left(RC-R[-1]C\right) + \left(RC-R[1]C\right) \bigg).$$

Next, what if the horizontal walls are longer than the vertical ones? Then more heat should be exchanged across the vertical walls than the horizontal ones.

We can change the coefficients of “vertical” differences in the above formula to reflect that: $$RC= RC - k*\bigg( .5*\left(RC-RC[-1]\right) + .5*\left(RC-RC[1]\right) + 2*\left(RC-R[-1]C\right) + 2*\left(RC-R[1]C\right) \bigg).$$ As a result, the material does travel farther -- as measured by the number of cells -- in the vertical direction (second image) than normal (first image):

However, if we stretch -- justifiably -- the cells in Excel horizontally (third image), the pattern becomes circular again!

In summary, the amount of heat exchanged between two rooms is proportional to:

• the temperature difference,
• the permeability of the wall (dually: the conductance of the pipe),
• the area of the wall that separates them (dually: the cross section of the pipe), and
• inversely, to the distance between the centers of mass of the rooms (dually: the length of the pipe).

Let's split the data into three categories:

• the adjacency of rooms (and the pipes) is the topology,
• the areas of the walls (and the lengths of the pipes) is the geometry, and
• the properties of the material of the walls (and the pipes) is the physics.

They are given by, respectively:

• the cell complex $K$,
• the Hodge star operator $\star^m :C^m(K)\to C^{n-m}(K^\star)$, and
• the $(n-1)$-form $k$ over $K$.

Suppose the geometry of space is supplied by means of the $m$-dimensional volume $|b|$ of each $m$-cell $b$ -- in both primal and dual complexes $K$ and $K^\star$. Now, the $m$th Hodge star is a diagonal matrix whose entries are the ratios of dual and primal volumes (up to a sign) in each dimension $m=0,1,...,n$: $$\star ^m_{ii}=\pm\frac{|a_i^\star|}{|a_i|}.$$

Similarly the geometry of time is supplied by means of the length $|t|$ of each time interval ($1$-cell) $t$ -- in both primal and dual complexes ${\mathbb R}$ and ${\mathbb R}^\star$. The Hodge star is a diagonal matrix whose entries are the reciprocals of the lengths of these intervals: $$\star ^1_{ii}=\pm\frac{1}{|t_i|}.$$

Recall that the right-hand side of our equation is the familiar second derivative (with respect to space) of an $n$-form with an extra factor $k$: $$(kU_x)_x:=(kU')' = d_x \star kd_x \star U.$$ As before, this expression is seen in the Hodge duality diagram, but this time the factors that come from the star operators are also shown: $$ \newcommand{\ra}[1]{\!\!\!\!\!\!\!\xrightarrow{\quad#1\quad}\!\!\!\!\!} \newcommand{\da}[1]{\left\downarrow{\scriptstyle#1}\vphantom{\displaystyle\int_0^1}\right.} \newcommand{\la}[1]{\!\!\!\!\!\!\!\xleftarrow{\quad#1\quad}\!\!\!\!\!} \newcommand{\ua}[1]{\left\uparrow{\scriptstyle#1}\vphantom{\displaystyle\int_0^1}\right.} % \begin{array}{ccccccccccc} & C ^{n-1}(K) & \ra{\quad d_x \quad} & C ^{n}(K) & \ni U \\ &\quad\quad \quad\ua{\times\frac{|a^{n-1}|k(a^{n-1})}{|b^{1}|}} & \ne & \quad\da{\times\frac{|b^0|}{|a^n|}} \\ & C ^{1}(K^\star)& \la{\quad d_x \quad} & C ^{0}(K^\star),\\ \end{array} $$ where $a^m$ is a primal $m$-cell and $b^m$ is a dual $m$-cell.

If we are to use the derivative, instead of the exterior derivative, with respect to time, we need to consider two issues. First, let's recall that when studying ODEs we used the function $q:C_1({\mathbb R})\to C_0({\mathbb R})$ given by $$q([i,i+1])=i,$$ to make the degrees of the forms, over time, match. Second, in addition to the above list, the amount of material exchanged between two rooms is also proportional to the length of the current time interval $|t|$. Then our PDE takes the form: $$d_tUq^{-1}=(kU_x)_x|t|,$$ with both sides $(n,0)$-forms. Consider the first derivative with respect to time: $$U_{t}:=U'=\star d_t U=\frac{1}{|t|}d_tU.$$

Definition. Given a geometric complex $K$ of dimension $n$, the diffusion equation is: $$U_t q^{-1}= (kU_x)_x,$$ where $U$ is an $(n,0)$-form over $K\times {\mathbb R}$.

The abbreviated version is below. $$\begin{array}{|c|} \hline \\ \quad U_t = (kU_x)_x \quad \\ \\ \hline \end{array}$$

Definition. An initial value problem (IVP) for diffusion is a combination of the diffusion equation and an initial condition (IC): $$U_t = (kU_x)_x,\ U(\cdot,A_0)=u_0\in C^n(K).$$ Then a $(n,0)$-form $U$ on $K\times {\mathbb R} \cap \{A\ge A_0\}$ that satisfies both conditions is called a (forward) solution of the IVP.

Since the exterior derivative with respect to time is simply the difference of values, a solution is easy to construct iteratively.

Theorem (Existence and Uniqueness). The solution to the IVP above exists and is given by: $$U(\cdot,A_0):=u_0,\quad U(\cdot,A+1):=U(\cdot,A)+(kU_x)_x(\cdot ,A)|[A,A+1]|,\ \forall A\ge A_0,$$ provided $\frac{1}{|a|}\in R$ for every cell $a$ in $K\sqcup K^\star$.

Exercise. Provide a weaker sufficient condition for existence.

Note: Boundary values problems lie outside the scope of this book.

Example. The choice of $$K={\mathbb R}^2,\quad R={\bf Z}_2,$$ produces this simple “cellular automaton”:

$\square$

Exercise. Derive the dual diffusion equation: “Suppose the material is located in the nodes of a graph in ${\bf R}^n$ and the amount of material is represented by a $0$-form $V$...”

Note: One can let the physics to be absorbed into the geometry. Since the transfer across a primal $(n-1)$-cell $a$ is proportional to its permeability $k(a)$, we can simply replace the volume $|a|$ of $a$ with $k(a)|a|$, provided $k\ne 0$. Then the right-hand side of our equation becomes the Laplacian (with respect to space) $U_{xx}=\Delta U$.

We test the performance of this equation below.

How diffusion patterns depend on the sizes of cells

Let's consider diffusion over a $1$-dimensional geometric cubical complex in ${\mathbb R}$.

The diagonal elements of the Hodge star operator for dimension $n=1$ are: $$\star ^1_{ii}=\frac{|\text{point}|}{|\text{edge}|} = \frac{1}{\text{length}} = \frac{1}{|a_{i}|}.$$ Then $\star ^1$ and $\star ^0$ give us our Hodge duality diagram: $$ \newcommand{\ra}[1]{\!\!\!\!\!\!\!\xrightarrow{\quad#1\quad}\!\!\!\!\!} \newcommand{\da}[1]{\left\downarrow{\scriptstyle#1}\vphantom{\displaystyle\int_0^1}\right.} \newcommand{\la}[1]{\!\!\!\!\!\!\!\xleftarrow{\quad#1\quad}\!\!\!\!\!} \newcommand{\ua}[1]{\left\uparrow{\scriptstyle#1}\vphantom{\displaystyle\int_0^1}\right.} % \begin{array}{ccccccccccc} & C ^{0}(K) & \ra{d} & C ^{1}(K) & \ni U \\ & \ua{\times\frac{|b|}{1}} & \ne & \da{\times\frac{1}{|a|}} \\ & C ^{1}(K^\star)& \la{d} & C ^{0}(K^\star),\\ \end{array} $$ where $a,b$ are primal and dual $1$-cells respectively. The right-hand side of our equation will have extra coefficients: the lengths of the cells appear twice as we make a full circle.

First, it's the length of the cell itself, $\frac{1}{|a|}.$ Then the equation will contain: $$\frac{U(a)}{|a|}.$$ That's the density of the material inside $a$.

Second, it's the lengths of the $1$-cells dual to the end-points of $a=AB$: $$\frac{1}{|A^\star|},\frac{1}{|B^\star|}.$$ The denominators are the lengths of the pipes.

Conclusion: The amount of material exchanged by a primal $1$-cell $a$ with its neighbor $a'$ is

• directly proportional to the difference of density in $a$ and $a'$, and
• inversely proportional to the length of the pipe that leaves $a$ for $a'$.

To confirm these ideas, we run an Excel simulation below:

With a single initial spike in the middle, we see that the amounts of material in the smaller cells on the right: $$|a_5|=|a_6|=|a_7|=|a_8|=|a_9|=1,$$ quickly become uniform, while the larger ones on the left: $$|a_1|=|a_2|=|a_3|=|a_4|=10,$$ develop slower.

Link to file: Spreadsheets.

Exercise. Find the speed of propagation in a uniform grid as a function of the length of the cell.

Exercise. Suppose we have two rods made of two kinds of metal soldered together. The cells will expand at two different rates when heated. Model the change of its geometry and illustrate with Excel. Hint: consider a single rod and then simply compare the geometry of the heated rod and that of the original rod.

For dimension $n=2$, we have the following diagonal entries in the case of a rectangular grid ${\mathbb R}^2$: $$\star ^2_{ii}=\frac{|\text{point|}}{|\text{rectangle}|} = \frac{1}{\text{area}} = \frac{1}{\Delta x \Delta y},$$ $$\star ^1_{ii}=\frac{|\text{edge}|}{|\text{edge}|} = \frac{\text{length}}{\text{length}} = \frac{\Delta y}{\Delta x }.$$ Then $\star ^2$ and $\star ^1$ give us our Hodge duality diagram: $$ \newcommand{\ra}[1]{\!\!\!\!\!\!\!\xrightarrow{\quad#1\quad}\!\!\!\!\!} \newcommand{\da}[1]{\left\downarrow{\scriptstyle#1}\vphantom{\displaystyle\int_0^1}\right.} \newcommand{\la}[1]{\!\!\!\!\!\!\!\xleftarrow{\quad#1\quad}\!\!\!\!\!} \newcommand{\ua}[1]{\left\uparrow{\scriptstyle#1}\vphantom{\displaystyle\int_0^1}\right.} % \begin{array}{ccccccccccc} & C ^{1}(K) & \ra{d} & C ^{2}(K) & \ni U \\ & \ua{\times\frac{|b|}{|a|}} & \ne & \da{\times\frac{1}{|\sigma|}} \\ & C ^{1}(K^\star)& \la{d} & C ^{0}(K^\star),\\ \end{array} $$ where $a,b$ are primal and dual $1$-cells respectively and $\sigma$ is a $2$-cell.

Just as above, we notice that $U(\sigma)/|\sigma|$ is the density inside $\sigma$. The second coefficient is new: $$\frac{|a|}{|a^\star|}.$$ We see here the $1$-cell that represents the wall and its dual that represents the pipe: $$\frac{\text{length of the wall}}{\text{length of the pipe}}.$$

Conclusion: The amount of material exchanged by a primal $2$-cell $\sigma$ with its neighbor $\tau$ is

• directly proportional to the difference of density in $\sigma$ and that of $\tau$,
• inversely proportional to the length of the pipe from $\sigma$ to $\tau$, and
• directly proportional to the length of the wall (thickness of this pipe).

If we use the $2\times 1$ grid, we have the following lengths of $1$-cells:

• horizontal primal: $|a|=2$, dual: $|a^\star|=1$;
• vertical primal: $|a|=1$, dual: $|a^\star|=2$.

Then the Excel formula is the same as the one discussed above: $$RC= RC -.0025*\bigg( .5*(RC-RC[-1]) + .5*(RC-RC[1]) + 2*(RC-R[-1]C) + 2*(RC-R[1]C) \bigg).$$ This formula, with a single source, produces a circular pattern on a spreadsheet with appropriately sized cells, as expected:

Exercise. Show that with the horizontal walls impenetrable, the transfer pattern will be identical to the $1$-dimensional pattern.

To summarize, the material flows from room $\sigma$ to room $\tau$ through what appears to be a pipe of length $|a^\star|$ and width $|a|$:

How diffusion patterns depend on the angles between cells

What if the wall between two rooms is slanted? Does it affect the amount of liquid that crosses to the other room?

If two walls are made of the same material (and of the same thickness), the slanted one will allow less liquid through it. In fact, what matters is the normal component of the flow across the wall.

We arrive to the same idea if, dually, we expand our pipe, $a^\star$, to the largest possible width, $|a|$:

Then $a$ is a (possibly non-orthogonal) cross-section of pipe $a^\star$.

Now, to make this more precise, consider the general case of an $n$-dimensional room, $\sigma$, with an $(n-1)$-dimensional wall, $a$, and still a $1$-dimensional pipe, $a^\star$. It is illustrated below for $n=3$:

As before, we are looking at the angle $\alpha$ between the normal $c=a^\perp$ of the wall $a$ and the pipe $a^\star$. Then

• the normal component of the flow is acquired by multiplying by the cosine of the angle between $a$ and $a^\star$.

Conclusion: The amount of material exchanged by a primal $n$-cell $\sigma$ with its neighbor $\tau$ is

• directly proportional to the difference of density in $\sigma$ and that of $\tau$: $U(\sigma)/|\sigma| - U(\tau)/|\tau|,$
• inversely proportional to the length $|a^\star|$ of the pipe $a^\star$ that leaves $\sigma$ for $\tau$,
• directly proportional to the area $|a|$ of the wall $a$, and
• directly proportional to the cosine of the angle between this pipe and this wall, $\cos \widehat{aa^\star}$.

Next, the iterative formula.

The boundary -- as a chain -- of a cell is the sum of its boundary cells taken with appropriate orientations: $$\partial \sigma=\sum_{a\in \partial \sigma} \pm a.$$

Therefore, the net change of material in cell $\sigma$ is the sum of the amounts exchanged through its boundary cells: $$U(\sigma,t+1)-U(\sigma,t)=\sum_{\partial \sigma} \frac{|a|}{|a^\star|}\cos \widehat{aa^\star} \left[ \frac{U(\sigma)}{|\sigma|} - \frac{U(\tau_a)}{|\tau_a|} \right],$$ where $$\tau_a:=\sigma+(\partial a^\star)^\star$$ is the $n$-cell that shares wall $a$ with $\sigma$.

It is easy to confirm that the right-hand side is nothing but the same expression $\star d \star d U$ as before. What has changed is the coefficients of the matrix of the Hodge star operator of the new geometry of the complex: $$\star^k_{ii}=\frac{|a_i|}{|a_i^\star|}\cos \widehat{a_ia_i^\star},$$ where $a_i$ is the $i$th $k$-cell. Recall that the feature of a non-rectangular geometry is that the angles $\cos\widehat{aa^\star}\ne 1$.

Exercise. Create a spreadsheet for the above formula and confirm the results below.

With this formula, we are able to study diffusion or heat transfer through an anisotropic or even irregular pattern, such as fabric, plant cells, or sunflower seeds:

Consider this trapezoid grid:

It is made from a square grid by turning the vertical edges by the same angle alternating left and right.

Exercise. (a) Compare the speed of the horizontal flow and the vertical to those of the original grid. (b) Verify the conclusion with a simulation.

Example. Consider now the rhomboid grid:

The pattern isn't circular anymore! It's not even aligned with the axis:

This is how this pattern is acquired. Even with slanted walls, the exchange with each of the four adjacent cells is identical. That's why the result of the Excel simulation is still circular (left):

The pattern is shown circular, however, only because Excel displays the shapes of the cells as squares. To see the true pattern then, we skew the resulting image (right).

The result may seem unexpected; after all, the physics appears to be symmetric. It is true that the the pattern of the flow is identical for the two directions along the grid, and therefore, the speed of propagation in either of the two directions is the same, just as before. However, what about the other directions? Below, the red and green segments are equal in length, but to follow the red one, we would have to cross more walls than for the green:

We conclude that the skewed pattern is caused by the “anisotropic” nature of the medium.

$\square$

Exercise. (a) Verify these conclusions with a simulation. (b) What are the directions of the fastest and the slowest propagation? Explain.

Exercise. To study a regular triangular grid one would need a non-cubical cell complex. Instead, use Excel to model diffusion on a triangular grid produced by diagonally cutting the squares.

The PDE of wave propagation

Previously, we studied the motion of an object attached to a wall by a (mass-less) spring. Imagine this time a string of objects connected with springs:

Just as above, we will provide the mathematics to describe the following three parts of the setup:

• the topology of the cell complex $L$ of the objects and springs,
• the geometry given to that complex such as the lengths of the springs, and
• the physics represented by the parameters of the system such as those of the objects and springs.

Let $u(t,x)$ measure the displacement from the equilibrium of the object associated with position $x$ at time $t$ (we will suppress $t$). It is an algebraic quantity: $$u=u(t,x)\in R.$$ As such, it can represent quantities of any nature that may have nothing to do with a system of objects and springs; it could be an oscillating string:

Here, the particles of the string are vertically displaced while waves propagate horizontally. Or, instead, pressure or stress may vary in a solid medium producing sound.

First, we consider the spatial variable, $x\in {\bf Z}$. We think of the array -- at rest -- as the standard $1$-dimensional cubical complex $L={\mathbb R}_x$. The complex may be given a geometry: each object has a (possibly variable) distance $h=\Delta x$ to its neighbor and the distance between the centers of the springs has length $\Delta x^\star$. We think of $u$ as a form of degree $0$ -- with respect to $x$.

According to Hooke's law, the force exerted by the spring is: $$F_{Hooke} = -k df,$$ where $df\in R$ is the displacement of the end of the spring from its equilibrium state and the constant, stiffness, $k\in R$ reflects the physical properties of the spring. Then, if this is the spring that connects locations $x$ and $x+1$, its displacement is the difference of the displacements of the two objects. In other words, we have: $$df=u(x+1) - u(x).$$ Therefore, the force of this spring is $$H_{x,x+1} = k \left[ u(x+1) - u(x) \right].$$ Since $k$ can be location-dependent, it is a $1$-form over $L$.

Now, let $H$ be the force that acted on the object located at $x$. There are two Hooke's forces acting on this object from the two adjacent springs: $H_{x,x-1}$ and $H_{x,x+1}$. Therefore, we have: $$\begin{array}{lll} H &= H_{x,x-1} &+ H_{x,x+1} \\ &= k \left[ u(x-1) - u(x) \right] &+ k\left[ u(x+1) - u(x) \right]\\ &=-(kd_xu)[x-1,x]&+(kd_xu)[x,x+1]. \end{array}$$

Next, we investigate what this means in terms of the Hodge duality. These are the duality relations of the cells involved: $$\begin{array}{llll} &[x-1,x]^\star &= \{ x-1/2 \},\\ &[x,x+1]^\star &= \{ x+1/2 \},\\ &\{x\}^\star &= [x-1/2,x+1/2]. \end{array}$$ Then the computation is straight-forward: $$\begin{array}{llll} H&= (kd_xu)([x+1,x]-[x,x-1])\\ &= (kd_xu)(\{ x+1/2 \}^\star -\{ x-1/2 \}^\star)\\ &=(\star kd_xu)(\{x+1/2\}-\{x-1/2\})\\ &= d_x(\star kd_xu)[x-1/2,x+1/2]\\ &=d_x(\star kd_xu)(x^\star)\\ &=\star d_x\star kd_xu(x)\\ &=\star d_xk^\star \star d_xu(x)\\ &=(k^\star u_x)_x(x). \end{array}$$

Note that for a constant $k$, we are dealing with the second derivative of the $0$-form $u$ with respect to space: $$u' '=\star d_x\star d_xu = \Delta u.$$ Compare it to the second derivative of a $1$-form $U$ with respect to space: $$U' '=d_x\star d_x\star U = \Delta U,$$ which we used to model $1$-dimensional diffusion. The difference is seen in the two different starting points in the same Hodge duality diagram: $$ \newcommand{\ra}[1]{\!\!\!\!\!\!\!\xrightarrow{\quad#1\quad}\!\!\!\!\!} \newcommand{\da}[1]{\left\downarrow{\scriptstyle#1}\vphantom{\displaystyle\int_0^1}\right.} \newcommand{\la}[1]{\!\!\!\!\!\!\!\xleftarrow{\quad#1\quad}\!\!\!\!\!} \newcommand{\ua}[1]{\left\uparrow{\scriptstyle#1}\vphantom{\displaystyle\int_0^1}\right.} % \begin{array}{ccccccccccc} u\in & C ^{0}({\mathbb R}_x) & \ra{d_x} & C ^{1}({\mathbb R}_x) &\ni U \\ & \ua{\star} & \ne & \da{\star} \\ & C ^{1}({\mathbb R}_x^\star)& \la{d_x} & C ^{0}({\mathbb R}_x^\star).\\ \end{array} $$ Then, $$d_xu' ' = (d_xu)' '.$$

Second, we consider the temporal variable, $t\in {\bf Z}$. We think of the time as the standard $1$-dimensional cubical complex ${\mathbb R}_t$. The complex is also given a geometry. It is natural to assume that the geometry has no curvature but each increment of time may have a different duration (and, possibly, $\Delta t \ne \Delta t^\star$). We think of $u$ as a form of degree $0$ with respect to $t$.

Now suppose that each object has mass $m$. Then, by the Second Newton's Law, the total force is: $$F=m \cdot a,$$ where $a$ is the acceleration. It is the second derivative with respect to time, i.e., this $0$-form: $$a=u_{tt}:=\star d_t \star d_t u . $$ The mass $m$ is a $0$-form too and so is $F$. Note that the stiffness $k$ is also a $0$-form with respect to time.

Now, since these two forces are equal, we have the wave equation of differential forms: $$\begin{array}{|c|} \hline \\ \quad m u_{tt} = (k^\star u_x)_x. \quad \\ \\ \hline \end{array}$$

If $k$ and $m$ are constant forms (and $R$ is a field), the wave equation takes a familiar shape: $$u_{tt}=\frac{k}{m}u_{xx}.$$

Exercise. Derive the dual (with respect to space) wave PDE.

Solutions of the wave equation

Now we will derive the recurrence relations.

First, we assume that the geometry of the time is “flat”: $\Delta t =\Delta t^\star.$

Then the left-hand side of our equation is $$m d_{tt} u = m\frac{u(x,t+1)-2u(x,t)+u(x,t-1) }{ (\Delta t)^2}.$$ For the right-hand side, we can use the original expression: $$ \star d_x \star k d_xu = k \left[ u(x-1) - u(x) \right] + k\left[ u(x+1) - u(x) \right].$$

Second, we assume that $k$ and $m$ are constant.

Now just solve for $u(x,t+1)$: $$\begin{array}{ll} u(x,t+1) &= 2u(x,t) - u(x,t-1) + \alpha\Big(u(x+1,t)-2u(x,t)+u(x-1,t)\Big), \end{array}$$ where $$\alpha :=(\Delta t)^2 \frac{k}{m}.$$

To visualize the formula, we arrange the terms in a table: $$\begin{array}{l|cccccc} & x-1 & x & x+1\\ \hline t+1 & & u(x,t+1) \\ t &=\alpha u(x-1,t) & +2(1-\alpha)u(t,x) & +\alpha u(x+1,t)\\ t-1 & & -u(x,t-1) \end{array}$$ Even though the right-hand side is the same, the table is different from that of the (dual) diffusion equation. The presence of the second derivative with respect to time makes it necessary to look two steps back, not just one. That's why we have two initial conditions.

Suppose, for simplicity, that $\alpha=1$.

Example. Choosing the simplified settings allows us to easily solve the following initial value problem: $$\begin{array}{lll} u_{tt}=u_{xx};\\ u(x,t)=\begin{cases} 1 &\text{ if } t=0,x=1;\\ 0 &\text{ if } t=0,x\ne 1;\\ 1 &\text{ if } t=1,x=2;\\ 0 &\text{ if } t=1,x\ne 2. \end{cases} \end{array}$$ Initially, the wave has a single bump and then it takes one step from left to right. The negative values of $x$ are ignored.

Now, setting $k=1$ makes the middle term in the table disappear. Then every new term is computed by taking an alternating sum of the three terms above, as shown below: $$\begin{array}{l|cccccccc} t \backslash x& 1 & 2 & 3 & 4 & 5 & 6 & 7 & ..\\ \hline 0 & 1 & 0 & 0 & 0 & 0 & [0] & 0 & .. \\ 1 & 0 & 1 & 0 & 0 & [0] & 0 & [0] & .. \\ \hline 2 & 0 & 0 & 1 & 0 & 0 & (0) & 0 & .. \\ 3 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & .. \\ 4 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & .. \\ 5 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & .. \\ .. & .. & .. & .. & .. & .. & .. & .. & .. \end{array}$$ The wave still has a single bump and it is running from left to right at speed $1$:

$\square$

Exercise. (a) Solve the two-sided version of the above IVP. (b) Set up and solve an IVP with $2$ bumps, $n$ bumps.

Exercise. Implement an Excel simulation for the case of non-constant $m$. Hint: you will need two buffers.

Next, we consider the case of finite string of springs:

• the array of identical weights ($m$ is constant) consists of $N$ weights,
• the weights are distributed evenly over the length $L := N\Delta x$,
• the total mass is $M := Nm$,
• the springs are identical ($k$ is constant),
• the total spring constant of the array is $K := k/N$,
• the geometry of the space is “flat”: $\Delta x =\Delta x^\star$.

We can re-write the above explicitly: $$\begin{array}{ll} m\frac{u(x,t+1)-2u(x,t)+u(x,t-1) }{ (\Delta t)^2} \\ \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad = k\left[u(x+1,t)-2u(x,t)+u(x-1,t)\right] \end{array}$$ Or, $$\begin{array}{ll} \frac{u(x,t+1)-2u(x,t)+u(x,t-1) }{ (\Delta t)^2} \\ \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad= \frac{k(\Delta x)^2}{m}\frac{u(x+1,t)-2u(x,t)+u(x-1,t)}{(\Delta x)^2}\\ \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad= \frac{k/N(N\Delta x)^2}{Nm}\frac{u(x+1,t)-2u(x,t)+u(x-1,t)}{(\Delta x)^2}\\ \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad= \frac{KL^2}{M}\frac{u(x+1,t)-2u(x,t)+u(x-1,t)}{(\Delta x)^2}. \end{array}$$

Solving for $u(x,t+1)$, we obtain the same formula: $$\begin{array}{ll} u(x,t+1) &= 2u(x,t) - u(x,t-1) + \alpha\left[u(x+1,t)-2u(x,t)+u(x-1,t)\right], \end{array}$$ with a new coefficient: $$\alpha :=\left(\frac{\Delta t}{\Delta x}\right)^2 \frac{KL^2}{M}.$$

Exercise. Set up and solve the IVP for the finite string, for the simplified settings. Hint: mind the ends.

Waves in higher dimensions

We next consider an array of objects connected by springs, both vertical and horizontal:

While the setup for the time is the same, the space is now the complex of the plane. Its geometry is given by the lengths of the edges: $\Delta x$ and $\Delta y$.

The forces exerted on the object at location $x$ are the four forces of the four springs attached to it: $$\begin{array}{llll} H & = H_{x,x-1} + H_{x,x+1} + H_{y,y-1} + H_{y,y+1}\\ & = k[ u(x-1,y) - u(x,y) ] \\ & + k[ u(x+1,y) - u(x,y) ] \\ & + k[ u(x,y-1) - u(x,y) ] \\ & + k[ u(x,y+1) - u(x,y) ] . \end{array}$$

We still consider only $0$- and $1$-forms but on the standard $2$-dimensional cubical complex $L={\mathbb R}^2$. Hodge duality is slightly more complicated here. As an example, these are a $1$-form $\varphi$ and its dual $1$-form $\varphi^*$:

Just as in the $1$-dimensional case, each bracketed term in the expression for $F$ is the exterior derivative: two with respect to $x$ and two with respect to $y$: $$\begin{array}{llll} k[ u(x-1,y) - u(x,y) ] &= kd_xu([x,x-1],y), \\ k[ u(x+1,y) - u(x,y) ] &= kd_xu([x,x+1],y), \\ k[ u(x,y-1) - u(x,y) ] &= kd_yu(x,[y,y-1]),\\ k[ u(x,y+1) - u(x,y) ] &= kd_yu(x,[y,y+1]). \end{array}$$ We can re-arrange these terms as they all start from $(x,y)$: $$H= kdu\left\{\begin{array}{llll} & & +\{x\}\times [y,y+1] \\ &+[x,x-1]\times \{y\} & & +[x,x+1]\times \{y\}\\ & & +\{x\}\times [y,y-1] \end{array}\right\},$$ where $d=(d_x,d_y)$. To get the duals of these edges, we just rotate them counterclockwise. Then, $$H= kd^\star u\left\{\begin{array}{llll} & & +[x^+,x^-]\times \{y^+\} \\ &+\{x^-\}\times[y^+,y^-] & & +\{x^+\}\times[y^-,y^+]\\ & & +[x^-,x^+]\times \{y^-\} \end{array}\right\},$$ where “$^\pm$” stands for “$\pm 1/2$”.

Observe that the dual edges are aligned counterclockwise along the boundary of the square neighborhood $[x^-,x^+]\times [y^-,y^+]$ of the point $(x,y)$. That neighborhood is the Hodge dual of this point. We recognize this expression as a line integral: $$H= \int_{\partial (\star (x,y))} \star kdu.$$

Now by the Stokes Theorem, the integral is equal to: $$H= \star d_x \star kd_xu .$$ Since the left-hand side is the same as before, we have the same wave equation: $$m u_{tt} = (k^\star u_x)_x.$$

In general, the medium may be non-uniform and anisotropic, such as wood:

We model the medium with a graph of springs and objects with a possibly non-trivial geometry:

We are now in the same place as with the diffusion equation. We split the data into three categories:

• 1. the adjacency of the springs (and the objects) is the topology,
• 2. the lengths of the springs (and the angles and the areas of their cross-sections) is the geometry, and
• 3. the Hooke's constants of the springs is the physics.

They are given by, respectively:

• 1. the cell complex $L$,
• 2. the Hodge star operator $\star^m :C^m(L)\to C^{n-m}(L^\star)$, and
• 3. the primal $1$-form $k\in C^1(L)$.

The two complexes dual to each other are shown below:

The total force that affects the object located at vertex $x$ in $L$ is: $$F= \star\int_{\partial (\star x)} \star kdu.$$ Therefore, we end up with the same wave equation as above. Its right-hand side is seen as the full circle in the Hodge duality diagram: $$ \newcommand{\ra}[1]{\!\!\!\!\!\!\!\xrightarrow{\quad#1\quad}\!\!\!\!\!} \newcommand{\da}[1]{\left\downarrow{\scriptstyle#1}\vphantom{\displaystyle\int_0^1}\right.} \newcommand{\la}[1]{\!\!\!\!\!\!\!\xleftarrow{\quad#1\quad}\!\!\!\!\!} \newcommand{\ua}[1]{\left\uparrow{\scriptstyle#1}\vphantom{\displaystyle\int_0^1}\right.} % \begin{array}{ccccccccccc} u\in & C ^{0}(L) & \ra{\quad d_x \quad} & C ^{1}(L) & \\ &\quad\quad \quad\ua{\times\frac{|a^{0}|}{|b^{n}|}} & \ne & \quad\da{\times\frac{|b^{n-1}|}{|a^1|}} \\ & C ^{n}(L^\star)& \la{\quad d_x \quad} & C ^{n-1}(L^\star),\\ \end{array} $$ where $a$ and $b$ are dual.

The geometry of the primal complex $L$ is given first by the lengths of the springs. What geometry should we give to the dual complex $K=L^\star$?

First, the meaning of the coefficients of the Hodge star operator are revealed to be represented by the stiffness $k(a)$ of spring $a$. Indeed, it is known to be directly proportional to its cross-sectional area $|b|$ and inversely proportional to its length $|a|$: $$k(a):=E\frac {|b|} {|a|},$$ where $E$ is the “elastic modulus” of the spring. Of course, when the angles are right, we have simply $b=a^\star$.

Second, we assign the reciprocals of the masses to be the $n$-volumes of the dual $n$-cells: $$|x^\star|:=\frac{1}{m(x)}.$$

Then, we have the simplified -- with the physics taken care of by the geometry -- wave equation, $$u_{tt}=(E^\star u_{x})_x,$$ where $E$ is the $1$-form of elasticity. Its right-hand side is seen as the full circle in the Hodge duality diagram modified accordingly: $$ \newcommand{\ra}[1]{\!\!\!\!\!\!\!\xrightarrow{\quad#1\quad}\!\!\!\!\!} \newcommand{\da}[1]{\left\downarrow{\scriptstyle#1}\vphantom{\displaystyle\int_0^1}\right.} \newcommand{\la}[1]{\!\!\!\!\!\!\!\xleftarrow{\quad#1\quad}\!\!\!\!\!} \newcommand{\ua}[1]{\left\uparrow{\scriptstyle#1}\vphantom{\displaystyle\int_0^1}\right.} % \begin{array}{ccccccccccc} u\in & C ^{0}(L) & \ra{\quad d_x \quad} & C ^{1}(L) & \\ &\quad\quad \quad\ua{\times\frac{1}{m(b^n)}} & \ne & \quad\da{\times k(a^1)} \\ & C ^{n}(L^\star)& \la{\quad d_x \quad} & C ^{n-1}(L^\star),\\ \end{array} $$

Wave propagation with a spreadsheet

The recurrence formula for dimension $1$ wave equation and constant $k$ and $m$ is: $$\begin{array}{ll} u(x,t+1) &= 2u(x,t) - u(x,t-1) + \alpha\left[u(x+1,t)-2u(x,t)+u(x-1,t)\right], \end{array}$$ where $$\alpha :=(\Delta t)^2 \frac{k}{m}.$$ We put these terms in a table to be implemented as a spreadsheet: $$\begin{array}{l|cccccc} & x-1 & x & x+1\\ \hline t-1 & & -u(x,t-1) \\ t &=\alpha u(x-1,t) & +2(1-\alpha)u(t,x) & +\alpha u(x+1,t)\\ t+1 & & u(x,t+1) \end{array}$$

The simplest way to implement the dynamics with the above formula in Excel is to use the first two rows for the initial conditions and then add one row for every moment of time. Then the Excel formula is:

• =R1C5*R[-1]C[-1]+2*(1-R1C5)*R[-1]C+R1C5*R[-1]C[1]-R[-2]C $\\$

Here cell R1C5 contains the value of $\alpha$.

Link to file: Spreadsheets.

Example. The simplest propagation pattern is given by $\alpha=1$. Below we show the propagation of a single bump, a two-cell bump, and two bumps:

$\square$

Exercise. Modify the spreadsheet to introduce walls into the picture:

Exercise. Modify the spreadsheet to accommodate non-constant data by adding the following, consecutively: (a) the stiffness $k$ (as shown below), (b) the masses $m$, (c) the time intervals $|\Delta t|$.

The $2$-dimensional case is treated with the formula given in the last subsection.

Example. The choice of $$K={\mathbb R}^2,\quad R={\bf Z}_2,$$ produces this simple “cellular automaton”:

$\square$