This site is being phased out.

Diffusion

From Mathematics Is A Science
Jump to navigationJump to search

The PDE of diffusion

A naive way to model heat propagation is to set up a grid of square rooms and then increment the temperature of each room by the average of the temperature of the four adjacent rooms. This simple model is implemented with an Excel simulation with the following short formula: $$RC= RC -.25 *\bigg( (RC-RC[-1]) + (RC-RC[1]) + (RC-R[-1]C) + (RC-R[1]C) \bigg),$$ Normally, only a proportion $k$ of this amount is shared.

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

Diffusion with Excel.png

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

Diffusion large.png

A more careful look reveals that to model heat transfer, we need to establish that a room exchanges heat with the four adjacent rooms. Moreover, the process has to follow 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.

Below we derive a PDE for discrete differential forms. Then, we will gradually develop a broad approach to the problem, which will be applicable to a variety of processes, and use the forms as simulations.

The set-up is as follows. A material (fluid, gas, heat, etc.) is contained in a grid of rooms and each room (an $n$-cell) exchanges the material with its neighbors through its walls.

Alternatively, one can think of pipes connecting the compartments or sinks in the centers of the rooms that average the amount of the material in the adjacent locations.

Duality for diffusion.png

Here, the rooms are blue, the walls are green, and the pipes are gray.

Of course, we recognize the two set-ups as Hodge dual!

We will the two approaches:

  • primary: rooms and walls; and
  • dual: sinks and pipes.

We can also think of pipes as passing through these walls.

We assume that the time increment is long enough

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

The forms will have two degrees -- with respect to location $x$ and with respect to time $t$.

The amount of material $U=U(a,t)$ is simply a number assigned to each room $a$ which makes it a $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.

More generally, 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 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 $(n,1)$-form, but over the dual complex.

Let

  • $d_t$ is the exterior derivative with respect to time (just the difference since the dimension is $1$); and
  • $d_x$ is the exterior derivative with respect to location.

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

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

Now, we need to express $F$ in terms of $U$. The flow $F(A^\star)=F^\star(A)$ through wall $A$ of room $a$ is proportional to the difference of amounts of material (or, more precisely, the density) in $a$ and the other room adjacent to $A$. So, $$F^\star(A) = - kd_x(\star U)(A^\star).$$ Here $k = k(A)$ is a $(n-1,0)$-form that represents the permeability of the wall $A$ at a given time. For realistic results choose: $$0\le k \le 1.$$ 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}$$

If $k$ happens to be constant, we are dealing with the Laplacian of $U$ with its second term missing because $C^{n+1}(K)=0$. Then the equation takes a more familiar form: $$d_t U = k\Delta U.$$

The right-hand side is given in this (non-commutative) diagram, 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} & C ^{n}(K) & \ni U \\ \da{\star} & \ne & \da{\star} \\ C ^{1}(K^\star)& \la{d} & C ^{0}(K^\star).\\ \end{array} $$ The multiplication by $k$ is implied. The second term of the Laplacian is missing because the next square on the right is trivial.

This was an outline. In the following, we develop both the mathematics (the generalized Hodge star) and the simulations for progressively more complex situations.

Simulation in dimensions 1 and 2

For the simplest Excel simulation, let's “trivialize” the above analysis.

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.
Duality for diffusion dim1.png

Since 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$. Positive flow at $A$ is from left to right and vice versa for $B$, then: $$d_t U(a)=-\left( F^\star(A)-F^\star(B) \right) = F^\star(A)-F^\star(B),$$ which is simply the exterior derivative in dimension $0$. 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$ is the permeability $(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). $$

If $k$ is constant, we have: $$d_t U =-k(U(b)+U(c)).$$ The result is a recursive formula for the simulation: $$U(a, t+1)= U(a,t) -k \Big[ U(a-1,t) + U(a+1,t) \Big].$$ The Excel setup is shown below:

Diffusion dim 1.png

The result after $1500$ iterations is shown:

Diffusion dim 1 iterated.png

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”. Here we use boundary conditions to substitute for the missing data.

Files:

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 forms will have two degrees -- with respect to location $x$ and with respect to time $t$. The amount of material $U=U(\tau,t)$ is a $2$-form with respect to location and a $0$-form with respect to time.

Closed 2cell.png

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. It can be written simply as $Qdt$, where $Q$ is a dual $1$-form with respect to space.

The “conservation of the matter” in 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$. 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 amounts of material in $\tau$ and the other $2$-cell adjacent to $a$. So, $$F^\star(a) = - kd_x(\star U)(a^\star).$$ Here $k = k(a)$ is a $(1,0)$-form that represents the permeability of the wall $a$ at a given time.

To see how the above analysis leads to the “naive” Excel simulation presented in the beginning, let's rewrite it as follows: $$U_{n+1}^{(i,j)}-U_{n}^{(i,j)} = -\frac{k}{4} \Big[ (U_{n}^{(i,j)}-U_{n}^{(i,j-1)})+(U_{n}^{(i,j)}-U_{n}^{(i,j+1)})+(U_{n}^{(i,j)}-U_{n}^{(i-1,j)})+(U_{n}^{(i,j)}-U_{n}^{(i+1,j)})\Big].$$ Now compare it to these two formulas from above: $$d_t U(\tau)=- F^\star(\partial\tau),$$ before the Stokes Theorem is applied, and $$F^\star(a) = - kd_x(\star U)(a^\star).$$ Then, the expression in the brackets is the integral of the exterior derivative of the flow over boundary of the cell. The $\frac{1}{4}$ coefficient is explained by the choice of geometry of the cell with $4$ edges of equal length, discussed later.

A simulation of heat transfer from single point on a larger space is shown below:

Heat transfer from single point.png

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

We simplified the general diffusion formula for the two specific cases. The general recursive formula: $$U(a, t+1)= U(a,t) + \left( d_x \star k(\cdot,t)d_x \star U \right)(a, t),$$ is implemented via $d_x$ and $\star$.

The dual equation with inflow

In the general cell complex case the material is still contained in $n$-dimensional cells (we assume that they form an $n$-dimensional combinatorial manifold) and each cell exchanges the material with its neighbors through its boundary that consists of its $(n-1)$-dimensional faces. The amount of material in a cell $U$ is a $(n,0)$-form, an $n$-form with respect to location (in the $n$-dimensional space) and a $0$-form with respect to time.

An alternative (dual) approach is to think of the nodes, i.e., $0$-cells, as the containers instead of the $n$-cells. The fluid is then exchanges via the pipes, i.e., edges or $1$-cells, that connect the nodes.

Let's derive the equation in this dual setting.

The amount of material $M$ is a dual $(1,0)$-form, a $1$-form with respect to location (on the grid) and a $0$-form with respect to time.

The outflow $F$ is a dual $(1,1)$-form: the amount of flow across an $(n-1)$-face (from the $n$-cell to its neighbor) per unit of time.

We will also include this time the inflow $S$ of the material to a given cell from the outside of the system which is given by a dual $(1,1)$-form. It can be written simply as $S=Pdt$, where $P$ is an $1$-form with respect to space.

Then the preservation of the material in cell $\sigma=A^\star$, where $A$ is the dual vertex, is given by $$d_t M(A,t) = −\int_{\partial A^\star} \star F(·,t) + S(A,t).$$

By the Stokes Theorem, we rewrite: $$d_t M(A,t) = −(d_x \star F)(A^\star,t) + S(A,t).$$ And, based on $\phi (\sigma ^\star)=\phi ^\star(\sigma)$, we re-write: $$d_t M = − \star d_x \star F + S.$$

Now, the flow $F$ through pipe $b$ is proportional to the difference of amounts of material at the end-points of $b$. So, $$F(b,t) = - k(b,t)d_x M(b,t).$$ Here $k$ is a $1$-form that may represent the thickness of the pipe. Substitute: $$d_t M= \star d_x \star kd_x M + S.$$

Exercise. Derive the primary equation with inflow.

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.

Cells with arbitrary geometry

We have modeled heat transfer with: $$RC= RC -.25 *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 -.25 *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):

Diffusion w rectangles.png

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

To sum up, the amount of heat exchanged between two rooms is proportional to:

  1. the temperature difference,
  2. the conductance of the wall,
  3. and, inversely, to the conductance of the pipe,.
  4. the area of the wall that separates them,
  5. and, inversely, to the distance between the centers of mass of the rooms (the length of the pipe).

On a deeper level, we split the data into three categories:

  1. the adjacency of rooms (and the pipes) is the topology,
  2. the areas of the walls (and the lengths of the pipes) is the geometry, and
  3. the conductance of the walls (and the pipes) is the physics.

They are given by, respectively, by:

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

Suppose the geometry is supplied by means of the $m$-dimensional volume $|b|$ of each $m$-cell $b$ -- in both primal and dual complexes. 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|}.$$

We realize that, since the geometry is taken care of by the Hodge star and, therefore, we don't have to change anything about the diffusion equation: $$d_t U = -d_x \star kd_x \star U.$$ As before, the right-hand side of the equation is seen in the Hodge duality diagram, but this time we indicate the multiples of the star operators: $$ \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{d} & C ^{n}(K) & \ni U \\ & \ua{\times\frac{|b^{1}|}{|a^{n-1}|}} & \ne & \da{\times\frac{|b^0|}{|a^n|}} \\ & C ^{1}(K^\star)& \la{d} & C ^{0}(K^\star),\\ \end{array} $$ where $a^m$ is a primary $m$-cell and $b^m$ is a dual $m$-cell.

We test this approach below.

Dimension 1: transfer depends on lengths

Let's consider diffusion over a $1$-dimensional geometric cubical complex such as this combination of primary and dual complexes:

Geometric complex dim 1.png

From the last section, the diagonal elements of the Hodge star operator for dimension $n=1$ are: $$\star ^1_{ii}=\frac{|point|}{|edge|} = \frac{1}{length} = \frac{1}{|a_{i}|}.$$

In the diffusion equation, we have $\star ^1$ and $\star ^0$ in our 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 primary and dual $1$-cells respectively.

It follows then that the right-hand side of our equation will have extra coefficients. The lengths will appear twice.

First, it's the length of the cell itself: $$\frac{1}{|a|},$$ which gives us $U(a)/|a|,$ the density of the material inside $a$.

Second, it's the lengths of the $1$-cells dual to the end-points of $a$: $$\frac{1}{|A^\star|},\frac{1}{|B^\star|},$$ which gives us the lengths of the pipes.

Conclusion:. The amount of material exchanged by a primary $1$-cell $a$ with its neighbors is

  • directly proportional to the difference of density in $a$ and those of its two neighbors, and
  • inversely proportional to the lengths of the two pipes that leave $a$.

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

Diffusion dim 1 w geometry.png

With a single initial spike in the middle, the result is that the amount 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.

Source file: diffusion_dim_1_w_geometry.xlsx

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

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 to the heated rod and that of the original rod.

Dimension 2: transfer depends on lengths

For dimension $n=2$, we have the following diagonal entries in the case of a rectangular grid: $$\star ^2_{ii}=\frac{|point|}{|rectangle|} = \frac{1}{area} = \frac{1}{\Delta x \Delta y},$$ $$\star ^1_{ii}=\frac{|edge|}{|edge|} = \frac{length}{length} = \frac{\Delta y}{\Delta x }.$$

In the diffusion equation, we have $\star ^2$ and $\star ^1$ in our 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 primary and dual $1$-cells respectively and $\sigma$ is a $2$-cell.

It follows then that its right-hand side of our equation will have these coefficients: $$\frac{1}{|\sigma|},$$ the area of the cell itself, and $$\frac{|a|}{|a^\star|},...$$ the length of the $1$-cell over the length of its dual. Again we observe that $U(\sigma)/|\sigma|$ is the density inside $\sigma$.

Conclusion: The amount of material exchanged by a primary $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 thickness of this pipe.

To illustrate the last point, let's consider our grid:

Duality and pipes on grid.png

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

  • horizontal primary $|a|=2$ and its dual $|a^\star|=1$;
  • vertical primary $|a|=1$ and its 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, as expected:

Diffusion w rectangles 2.png

So, the material flows left to right through what appears to be a pipe of length $|a^\star|$ and width $|a|$:

Duality and pipes on grid 2.png

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

Dimension 2: transfer depends on angles

What if the wall is sloped? Does it affect the amount of material that crosses to the other room?

Duality and pipes on grid 3.png

The answer is, of course, Yes. If we assume that the liquid flows along the line (the pipe $a^\star$) that connects the centers of the rooms, then what matters is the normal component of the flow.

Diffusion through sloped wall.png

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

Duality and pipes 2.png

Suppose $a$ is the wall (edge) and $b=a^\star$ is the pipe between the centers. If we think of them as vectors, the flow is proportional to the sine of the angle between $a$ and $b$. Therefore, we need an adjustment coefficient $\sin \widehat{ab}$. Note that when the metric tensor is Euclidean, $a$ and $b$ are perpendicular; then the coefficient is $1$.

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

  • directly proportional to the difference of density in $\sigma$ and $\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 sine of the angle between this pipe and this wall.

We assume that the cells in this cell complex have no self-intersections. Then the boundary -- as a chain -- of a cell is the sum of its boundary cells taken with appropriate orientations: $$\partial x=\sum_{a\in \partial x} \pm a.$$ In particular, if $x$ is an edge, its boundary is the difference of its end-points: $$\partial x=B-A.$$

Therefore, the net change of material in cell $\sigma$ is the sum of the amounts exchanged with its boundary cells: $$U(\sigma,t+1)-U(\sigma,t)=\sum_{(\partial \sigma)^\star} k\frac{|a|}{|a^\star|}\sin \widehat{aa^\star} \left[ \frac{U(\sigma)}{|\sigma|} - \frac{U(\tau)}{|\tau|} \right],$$ where the summation is over: $$(\partial \sigma)^\star =\sum a,\ a\in C_1(K^\star),$$ and $$\tau:=\sigma-\partial a^\star.$$

Exercise. Confirm that, up to the coefficients, the right-hand side is nothing but $\star d \star d U$.

Exercise. Show that when $a⊥a^\star$, we have the original formula.

Consider this trapezoid grid:

Trapezoid grid.png

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 will to those of the original grid. (b) Verify the conclusion with a simulation.

Example. Consider now the rhomboid grid:

Rhomboid grid.png

Whatever the pattern the flow has at these edges, it is the identical for both of the directions along the grid. Therefore, the speed of propagation in any of the two directions along the grid is the same, just as on the square grid. This is confirmed by a simulation of a single source diffusion, which produces a skewed, non-circular pattern:

Square vs rhomboid diffusion.png

The anisotropy of the pattern is explained by the anisotropy of the setup. Below, the red and green segments are equal in length, but to follow the red one on would have to cross $9$ walls vs. only $4$ for the green:

Diffusion on rhombus anisotropy.png

$\square$

Exercise. Consider the triangular grid.

The Hodge star operator

Next, we consider the general case: $n$-dimensional rooms, $\sigma$, with $(n-1)$-dimensional walls, $a$, but still $1$-dimensional pipes, $a^\star$.

We use the known fact that the amount of liquid that crosses the surface of a hyperplane is proportional to:

  • the component of the flow orthogonal to the surface; and
  • the area of the this part of the surface.
Duality and pipes.png

To be able to compute these numbers, we need to capture the common geometry of primary and dual complexes.

It is insufficient to have the two metric tensors, for a $K$ and $K^\star$: $$\begin{array}{ccc} <\cdot,\cdot>(\cdot): T(K) \to {\bf R},\\ <\cdot,\cdot>(\cdot): T(K^\star) \to {\bf R}, \end{array}$$ as they only give the angles between cells of the same complex.

It is also insufficient to have the Hodge star operators, as defined previously, provided for the primary and dual complexes: $$\begin{array}{ccc} \star^k : C^k(K)\to C^{n-k}(K^\star),\\ \star^k : C^k(K^\star)\to C^{n-k}(K), \end{array}$$ as they only keep record of the volumes of the cells.

The angle between $1$- and $(n-1)$-dimensional subspaces is always understood as the angle between the direction vector of the former and the normal vector of the letter: $$<a^\star, a>:=<a^\star, a^\perp>.$$

Exercise. In an inner product space, what is the angle between two subspaces of complementary dimensions?

We can then extend the metric tensor to include the inner product of any pair of dual cells for all dimensions since $\dim a +\dim a^\star =n$.

To accommodate these changes, how should we adjust the coefficients of the star operator? The star operator is still just a diagonal matrix with positive entries; however, we will provide a new interpretation.

We have realized that the formula for the Hodge star operator (the ratio of the dual and primal volume) that we have used is only applicable to rectangular grids. The necessary generalization is below.

Definition. The Hodge star operator is a diagonal matrix whose $i$th entry is the ratio of the dual volume and the projection of the primal volume of the $i$th cell. In other words, it is $$\star ^m_{ii}:=\frac{|a^\star_i|^2}{<a_i^\star, a_i>},$$ where $\dim a_i =m$.

Exercise. Compute the Hodge star operator for the regular triangular grid.

Exercise. What is the geometry of a complex the Hodge star of which is the identity?


One-dimensional wave

The wave equation in the one-dimensional case is derived below.

We studied previously the motion of a weight attached to a wall by a (mass-less) spring. Imagine this time a string of weights connected with springs:

Array of masses.png

As before, we will separate the following from each other:

  • the topology of the cell complex,
  • the geometry given to that complex, and
  • the physics represented by the parameters of the system.

Let $u(t,x)$ measure the displacement from the equilibrium of the weight 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 weights and springs:

Wave wo springs.png

First, we consider the spacial variable, $x\in {\bf Z}$. We think of the array -- at rest -- as the standard $1$-dimensional cubical complex ${\mathbb R}_x$. The complex may also be given a geometry: each weight 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$.

R complex.png

According to Hooke's law, the force exerted by the spring is: $$F_{Hooke} = -k df,$$ where $df\in R$ is the the displacement of the spring from its equilibrium state and the constant $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 weights. 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.

As we see, only the adjacency of the weights matters as the lengths $h$ of the edges of the complex don't come into play. Then, we can think of the geometry of ${\mathbb R}$ as trivial: $$\Delta x=\Delta x^\star=1.$$

Now, let $H$ be the force that acted on the weight located at $x$. There are two Hooke's forces affected this weight 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}{ll} [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). \end{array}$$

Note that for a constant $k$, we are dealing with the second derivative of $0$-form $u$ with respect to space: $$u^0_{xx}:=\star d_x\star d_xu.$$ Compare it to the second derivative of a $1$-form $U$ with respect to space: $$U^1_{xx}:=d_x\star d_x\star U = \Delta U,$$ which we used to model $1$-dimensional diffusion. The difference is seen in 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^0_{xx}=(d_xu)^1_{xx}.$$

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 weight 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 = (\star d_t)^2u. $$ The mass $m$ is a $0$-form too and so is $F$.

Now, since these two forces are equal, we have the wave equation of differential forms: $$\begin{array}{|c|} \hline \\ \quad m (\star d_t)^2 u = \star d_x\star kd_xu. \quad \\ \\ \hline \end{array}$$

If $k$ is a constant form, the wave equation takes a more familiar shape: $$u_{tt}=\frac{k}{m}u_{xx}.$$

Solutions of the wave equation

Now we will derive the recurrence equations.

First, we assume that the geometry of the time is trivial: $$\Delta t =\Delta t^\star.$$ Then the left-hand side of our equation is $$m (\star d_t)^2 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].$$

Now just solve for $u(x,t+1)$: $$\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}.$$

Or, we can 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}$$

The table is different from that of the 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 will need two buffers for our simulation.

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 it is 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$:

Traveling wave.png

$\square$

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

In a more general setting, the mass $m$ and the Hooke's constant $k$ are determined by the physics of the problem to be solved. Most commonly this data comes in the digital format, as a table of numbers. These numbers give you the values of $m$ and $k$ as they vary with location.

Exercise. Implement an Excel simulation for the case of non-constant $m$.

Next, we add geometry to the system and 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 trivial: $\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}$$

If solve 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 consdier an array of weights with springs both vertical and horizontal:

Array of masses dim 2.png

While the setup for the time is the same, the space is now given by the complex ${\mathbb R}^2$. Its geometry is given by the lengths of the edges: $\Delta x$ and $\Delta y$. We will see, as before, that the geometry can be trivial $\Delta x=\Delta y=1$.

The forces exerted on the weight at location $x$ are the four forces of the four springs attached to the weight: $$\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 a $2$-dimensional cubical complex. Hodge duality is slightly more complicated here. As an example, these is a $1$-form $\varphi$ and its dual $1$-form $\varphi^*$:

Hodge duality of forms.png

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, this is: $$H= \star d \star kdu .$$

Since the left-hand side is the same as before, we have the same equation: $$m (\star d_t)^2 u = \star d \star kdu.$$

In the general case, we deal with a graph:

Graph of masses and springs.png

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 weights) is the topology,
  2. the the lengths of the springs (and areas of the walls) is the geometry, and
  3. the Hooke's constant of the springs (and the area of the walls) is the physics.

They are given by, respectively, by:

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

Then the force looks like this $$F= \star\int_{\partial (\star x)} \star kdu.$$ Therefore we end up with the same equation.

For the right-hand side, it is, with the dimensions specified, $$\star_0\partial _{0}\star _1 k^1\cdot\partial _1^T u^0 .$$ All the parts are available. Recall that $\star d$ is equal to and $d$ is the transpose of the boundary operator $\partial$.

Exercise. Show that this matrix has $1$s on the main diagonal and $-1$s above it. Same for $d^\star_t$ and $d_t$.