This site is devoted to mathematics and its applications. Created and run by Peter Saveliev.

# More ODEs

### From Mathematics Is A Science

## Contents

## 1 Vector fields and systems of ODEs

The IVP we considered in chapter 1 was: $$df=Pf,\ f(A_0)=x_0\in R.$$ Here a $0$-form $f$, which represents the variable quantity to be found, is defined on the ray ${\mathbb R} \cap \{A\ge 0\}$ and satisfies the following: $$df\Big(A,[A,A+1] \Big)=P(f(A)),\ \forall A \in {\bf Z}^+.$$ A solution is given by $$f(0):=x_0,\ f(A+1):=f(A)+P(f(A)),\ \forall A\ge 0.$$

What if we have *two* variable quantities?

**Example.** The simplest example is as follows:

- $x$ is the horizontal location, and
- $y$ is the vertical location. $\\$

Both are $0$-forms on ${\mathbb R}$. Their respective displacements are $dx$ and $dy$. If either depends only on its own component of the location, the motion is described by this pair of ODEs of forms: $$dx=g(x),\ dy=h(y).$$ The solution consists of two solutions to the two, unrelated, ODEs either found with the formula above. $\square$

The solution is easy because the two components are independent of each other.

**Example (predator-prey model).** As an example of quantities that do interact, let $x$ be the number of the prey and $y$ be the number of the predators in the forest. Let $R:={\bf R}$.

The prey have an unlimited food supply and reproduces geometrically as described above -- when there is no predator. Therefore, the gain of the prey population per unit of time is $\alpha x$ for some $\alpha\in {\bf R}^+$. The rate of predation upon the prey is assumed to be proportional to the rate at which the predators and the prey meet, which, in turn, is proportional to $xy$. Therefore, the loss per unit of time is $\beta x y$ for some $\beta\in {\bf R}^+$. Then the prey's ODE is: $$dx = \alpha x - \beta x y.$$

The predators have only a limited food supply, i.e., the prey. We use, again, the rate of predation upon the prey proportional to $xy$. Then the gain of predator population per unit of time is $\delta xy$ for some $\delta\in {\bf R}^+$. Without food, the predators die out geometrically as described previously. Therefore, the loss per unit of time is $\gamma y$ for some $\gamma\in {\bf R}^+$. The predator's ODE is: $$dy = \delta xy - \gamma y.$$

This system of these two, dependent, ODEs cannot be solved by reusing the formula for the $1$-dimensional case. $\square$

**Exercise.** (a) Set up and solve the IVP for this system. (b) Find its equilibria.

**Example.** We also take as a model a fluid flow. The “phase space” ${\bf R}^2$ is the space of all possible locations. Then the position of a given particle is a function $u:{\bf Z}^+\to {\bf R}^2$ of time $t\ge 0$. Meanwhile, the dynamics of the particle is governed by the velocity of the flow, at each moment of time. Let $V:{\bf R}^2\to {\bf R}^2$ be a function. Then the velocity of a particle if it happens to be at point $u$ is $V(u)$ and the next position may be seen as $u+V(u)$.

$\square$

A *vector field* is given when there is a vector attached to each point of the plane:
$${\rm point} \mapsto {\rm vector}.$$
Then it is just a function
$$V : {\bf R}^2 \to {\bf R}^2 .$$

Further, one can think of a vector field as a system of time-independent ODE: $$\begin{cases} dx=g(x,y),\\ dy=h(x,y), \end{cases}$$ where $g,h:{\bf R}^2 \to {\bf R}$ are functions and $x,y:{\bf Z} \to {\bf R}$ are $0$-forms. The corresponding IVP adds an initial condition: $$\begin{cases} x(t_0)=x_0,\\ y(t_0)=y_0. \end{cases}$$ As before, a (forward) solution to this IVP is a $0$-form on the ray ${\mathbb R} \cap \{A\ge 0\}$ that satisfies, $\forall A \ge 0$: $$\begin{array}{llll} dx\Big(A,[A,A+1] \Big)=g(x(A),y(A)),\\ dy\Big(A,[A,A+1] \Big)=h(x(A),y(A)). \end{array}$$

Next, we combine the two variables to form a vector: $$\begin{array}{lll} u:=(x,y)\in {\bf R}^2,\\ du:=(dx,dy)\in {\bf R}^2. \end{array}$$ Note that the new vector variables still form a ring!

Then the setup given at the top of the subsection: $$df=Pf,\ f(A_0)=x_0\in R,$$ is now applicable to vector fields: in all the three examples, we set $$R:={\bf R}^2.$$

More generally, a vector field supplies a *direction to every location*, i.e., a vector to every point. Thus, there is one vector at each point picked from a whole vector space:

In a cell complex $X$, this is a function that give every vertex an adjacent edge. The adjacent edges of vertex $A$ form a $1$-dimensional *tangent space*
$$T_A:=\{AB\in X: B\in X\}.$$

**Definition.** A *vector field* over complex $X$ is any function
$$V:X\to C_1(X)$$
such that
$$V(A)=AB\in X,$$
where $B=B(A)$ is a vertex adjacent to $A$ that depends on $A$.

Then the motion produced by the vector field takes a particle from vertex to vertex along an edge, or simply: $$A \mapsto B.$$ Its solution can be written as an iteration: $$A_{n+1}:=A_n+\partial_1 V(A_n),$$ where $\partial_1:C_1(X)\to C_0(X)$ is the boundary operator of $X$.

**Exercise.** Verify the last statement. What kind of function is this?

## 2 Simulating a flow with a spreadsheet

Spreadsheets, such as Excel, can be used to model various dynamics and some of the processes can also be visualized.

We represent a flow with a vector field starting with *dimension* $1$; i.e., $R:={\bf Z}$. This vector field is the right-hand side function $P:{\bf Z}\to {\bf Z}$ of our ODE:
$$df=Pf.$$
The number in each cell tells us the velocity of a particle if it is located in this cell. The values of the velocities are provided in the spreadsheet. They were chosen in such a way that the particle can't leave this segment of the spreadsheet. They can be changed manually.

Then one can trace the motion of a single particle given by a $0$-form $f:{\bf Z}^+\to {\bf Z}$, shown as a red square:

Here, the velocities are seen in the elongated cells (edges), while the locations of the particles are the small square cells (vertices). Also,

- “x” stands for the current position,
- “o” for the previous position. $\\$

The formula is as follows: $$\texttt{ = IF( R[-1]C=R9C3,' 'x' ', IF(R[-1]C=R9C2, ' 'o' ',' ' ' ') )}$$ The data, as it changes, is shown on the left.

For motion simulation, we carry out one iteration for every two pushes of the button. This is an example of such dynamics with three iterations:

Link to file: Spreadsheets.

For the file to work as intended, the settings should be set, under “Formulas”:

- Workbook Calculation: Manual,
- Maximum Iterations: $1$.

To see the motion of the particle, simply push “Calculate now” as many times as necessary.

**Exercise.** Modify the code so that the spreadsheet models motion on the circle, no matter what velocities are chosen. Hint: choose an appropriate ring $R$.

**Exercise.** Modify the code so that, no matter what velocities are chosen, the object stays within, say, interval $[0, ...,7]$, and bounces off its ends.

Next, we simulate a flow with a vector field in *dimension* $2$; i.e., $R:={\bf Z}^2$. This vector field is the right-hand side function $P:{\bf Z}^2\to {\bf Z}^2$ of our ODE. Each pair (right and down) of elongated cells starting at every square cell contains a pair of numbers that represents the velocity of a particle located in this cell.

The vector field of dimension $2$ is made of two vector fields of dimension $1$ (on the same square) independently evaluated for the $x$ and $y$ variables. The components of the velocities are chosen to be only $-1$, $0$, or $1$ (illustrated with the little arrows) so that the jump is always to an adjacent location. They were also chosen in such a way that the particle can't leave this rectangle.

We trace the motion of a single particle given by a $0$-form $f:{\bf Z}^+\to {\bf Z}^2$, shown as a red “x”.

Several iterations are shown in a single picture below:

**Exercise.** Modify the code to model motion on (a) the cylinder and (b) the torus.

**Exercise.** Modify the code so that, no matter what velocities are chosen, the object stays within, say, rectangle $[0,8]\times [0,7]$, and bounces off its walls.

## 3 ODEs of cell maps

So far, we have represented motion by means of a $0$-form, as a function
$$f:{\mathbb R}^+\to R,$$
where $R$ is our ring of coefficients. Now, is this motion *continuous*?

Here, $R$ is our space of locations and it is purely algebraic. Even if we extend $f$ from ${\mathbb R}^+$ to ${\bf R}^+$, without a topology on $R$, the continuity of such an $f$ has no meaning! In contrast, motion is typically represented by a *parametric curve*, i.e., a continuous map
$$p:{\bf R}^+\to U\subset {\bf R}^n,$$
with ${\bf R}^+$ serving as the time and $U$ as the space.

If we choose to recognize $R={\bf R}^n$ as a topological space with the Euclidean topology, then, on the flip side, its subsets, such as $U={\bf R}^n\setminus \{0\}$, may lack an algebraic structure. For example, these rings aren't *rings*:

Then we can't use forms anymore.

To match the standard approach, below we represent motion in cell complex $K$ by a cell map, $$f:{\mathbb R} \supset I \to K,$$ where $I$ is a cell complex, typically an interval.

We can see forms if we look at the curve one point at a time:

Next, we will demonstrate how these curves are generated by vector fields. The analysis largely follows the study of ODEs of forms discussed previously, in spite of the differences. Compare their derivatives in the simplest case:

- for a $0$-form $f:K \to R$, we have

$$df\Big([A,A+1] \Big):=f(A+1)-f(A) \in R;$$

- for a cell map $f:{\mathbb R} \to {\mathbb R}$, we have

$$f'\Big([A,A+1] \Big):= f_1\Big([A,A+1] \Big) \in C_1({\mathbb R};R).$$

**Definition.** Suppose $P$ is a vector field on complex $K$. Then an *ODE of cell maps* generated by $P$ on complex $K$ is:
$$f'=Pf,$$
and its *solution* is a cell map $f:{\mathbb R}\supset I \to K$ that satisfies:
$$f'\Big([A,A+1] \Big)=P(f(A)),\ \forall A\in I\cap {\bf Z}.$$

For $K={\mathbb R}$, we can simply take the examples of $1$st order ODEs of forms given in this section and use them to illustrate ODEs of maps. We just need to set $R:={\bf Z}$. The difference is, this time the values of the right-hand side $P$ can only be $1$, $0$, or $-1$.

For $K={\mathbb R}^2$, the cubical case illustrated with a spreadsheet above doesn't apply anymore because a cell map has to follow the edges and can't jump diagonally. Then one of the components of the vector field has to be $0$, as shown below:

**Example.** We change the values of the vector field in the spreadsheet above as described above in order to model an ODE of cubical maps. Several iterations can be shown in a single picture, below:

$\square$

**Exercise.** Can the original spreadsheet be understood as a model of an ODE of cell maps?

Thus, if we represent motion with an ODE of cell maps, *both time and space are topological*.

**Definition.** An *initial value problem (IVP)* on complex $K$ is a combination of an ODE and an *initial condition (IC)*:
$$f'=Pf,\ f(A_0)=x_0\in K,\ A_0\in {\bf Z},$$
and a solution $f:{\mathbb R}\cap\{t\ge A_0\} \to K$ of the ODE is called a *solution of the IVP* if it satisfies the initial condition.

**Theorem (Uniqueness).** The solution to the IVP above, if exists, is given iteratively by
$$f(A_0):=x_0,\ f(A+1):=f(A)+\partial P(f(A)),\ \forall A\ge A_0.$$

**Proof.** Suppose $f$ is a solution, $f'=Pf$. Then we compute:
$$\begin{array}{llllllll}
&f_0(A+1)-f_0(A)&=f_0((A+1)-(A))\\
&&=f_0\Big(\partial [A,A+1] \Big)\\
&&=\partial f_1[A,A+1]\\
\hspace{.25in}&&=\partial P(f(A)).&\hspace{.25in}\blacksquare
\end{array}$$

Just as before, only the right-sided uniqueness is guaranteed: when two solutions meet, they stay together.

What about existence? We have the formula, but is this a cell map? It can only be guaranteed under special restrictions.

**Theorem (Existence).** If the values of the vector field $P$ are *edges* of $K$, a solution to the IVP above always exists.

**Proof.** For $f$ to be a cell map, $f(A)$ and $f(A+1)$ have to be the endpoints of an edge in $K$. Clearly, they are, because $f(A+1)-f(A) =\partial P(f(A))$. $\blacksquare$

**Definition.** For a given ODE, the *forward propagation map* of depth $c\in {\bf Z}^+$ at $A_0$ is a map
$$Q_c:K \to K$$
defined by
$$Q_c(x_0):=f(A_0+c),$$
where $f$ is the solution with the IC:
$$f(A_0)=x_0.$$

**Exercise.** Prove that if the vector field $P$ is time-independent, $Q_c$ is independent of $A_0$.

**Exercise.** Compute the forward propagation map for the vector field illustrated above.

**Exercise.** Suppose $P$ is edge-valued. Is the resulting forward propagation map $Q_c:K\to K$ a cell map?

## 4 ODEs of chain maps

A cell map can't model jumping diagonally across a square.

That's why the existence theorem above requires the values of the vector field to be edges. Then, in the case of ${\mathbb R}^n$, all but one of the components must be $0$.

This is a significant limitation of ODEs of cell maps, even in comparison to ODEs of forms.

The issue is related to one previously discussed: cell map extensions vs. chain map extensions (subsection V.3.10). Recall that in the former case, extensions may require subdivisions of the cell complex. The situation when the domain is $1$-dimensional is transparent:

In the former case, we can create a cell map:
$$g(AB):=XY,$$
by extending its values from vertices to edges. In the latter case, an attempt of cell extension (without subdivisions) fails as there is no *single* edge connecting the two vertices. However, there is a *chain* of edges:
$$g(AB):=XY+YZ.$$

Now, ODEs.

Once again, a vector field at vertex $X$ only takes values in the tangent space $T_X(K)$ at $X$. We need to incorporate the possibility of moving in the direction $XY+YZ$ (while $YZ\not\in T_X(K)$). From this point on, the development is purely algebraic.

Suppose we are given a (short) chain complex $C$ to represent *space*:
$$\partial: C_1 \to C_0.$$
Here, of the two modules, $C_0$ is meant to represent the chains of locations and $C_1$ the chains of directions. Even though all the information is contained in the complex, we still define an analog of the tangent bundle to demonstrate the parallels with ODEs of forms and cell maps.

**Definition.** A *chain field* over chain complex $C$ is any linear map
$$W:C_0\to C_1.$$

Here is an illustration of possible values of $W$ at vertex $X$ (even though the chains don't have to start at $X$):

The dynamics produced by the chain field is a sequence of $0$-chains given by this iteration: $$X_{n+1}:=X_n+\partial W(X_n).$$

The following is obvious.

**Proposition.** A vector field
$$V:K^{(0)}\to C_1(K)$$
over a cell complex generates -- by linear extension -- a chain field:
$$W:C_0(K)\to C_1(K).$$

With the new definition, the motion produced by such a chain field is as before. For example, if $W(X)=XY$, then $$X_1= X+\partial W(X) = X+(Y-X) =Y.$$

**Definition.** Suppose $W(\cdot)$ is a chain field on a chain complex $C$. Then an *ODE of chain maps* generated by $W$ on complex $C$ is:
$$g_1=Wg_0,$$
and its *solution* is a chain map
$$g=\{g_0,g_1\}:C(I) \to C,$$
for some subcomplex $I \subset {\mathbb R}$, that satisfies the equation
$$g_1\Big([A,A+1]\Big)=W(g_0(A)),\ \forall A\in I\cap {\bf Z}.$$

**Example.** Under this definition, a step of the moving particle on a square grid can be diagonal, without explicitly following the vertical and horizontal edges. For example, to get from $X=(0,0)$ to $Z=(1,1)$ without stopping at $Y=(0,1)$, set
$$W(X):=(X,XY)+(Y,YZ) \in T(C).$$
Indeed, we have
$$\hspace{.25in} X\mapsto X+\partial W(X) = X+ (Y-X)+(Z-Y)=Z. \hspace{.25in}\square$$

Since every cell map generates its chain map, we have the following:

**Proposition.** If a map $f$ is a solution of the ODE of cell maps generated by a vector field $V$, then $g:=f_*$ is a solution of the ODE of chain maps generated by the chain field $W=V_*$.

**Definition.** An *initial value problem* is a combination of an ODE and an *initial condition*:
$$g_1=Wg_0,\ g_0(A_0)=X_0\in C,\ A_0\in {\bf Z},$$
and a solution $g:C({\mathbb R}\cap\{t\ge A_0\}) \to C$ of the ODE is called a (forward) *solution of the IVP* if it satisfies the initial condition.

**Theorem (Existence and Uniqueness).** The only (forward) solution to the IVP above is given iteratively by
$$\begin{array}{lllllll}
g_0(A_0):=X_0,\quad \text{ and for all } A\ge A_0,\\
g_0(A+1):=g_0(A)+\partial W(g_0(A)),\\
g_1\Big([A,A+1]\Big):=W(g_0(A)).
\end{array}$$

**Proof.** The linearity of $W$ on the second argument implies that $g_0,g_1$ so defined form a chain map. The proof that the ODE holds is identical to the proof of the uniqueness theorem in the last subsection. $\blacksquare$

The uniqueness holds only in terms of chains. When applied to cell complexes, the result may be that the initial location is a vertex $X_n$ but the next “location” $X_{n+1}$ is not.

**Example.** Where do we go from $X=(0,0)$ if
$$W(X):=(X,XY)+(X,XZ),$$
with $Y=(0,1),Z=(1,1)$? This is the result:
$$\hspace{.2in} X\mapsto X+\partial W(X) = X+ (Y-X)+(Z-X)=-X+Y+Z. \hspace{.2in}\square$$

This isn't a cell map because the output isn't a single vertex! The latter can be guaranteed though under certain constraints on the equation.

## 5 Cochains are chain maps

Since every cell map generates its chain map, the existence theorem from the last subsection is still applicable: if a vector field $V$ is edge-valued, a solution to the IVP for maps always exists. Below is a stronger statement.

**Theorem (Cell maps as solutions).** If a chain field over the chain complex $C(K)$ of a cell complex $K$ satisfies the condition:
$$\partial W(A,X)=Y-X,\text{ for some } Y\in K^{(0)},$$
then the solution to the IVP for chain maps on $K$ is a cell map.

The difference between a cell map as a solution of

- an ODE of cell maps and
- an ODE of chain maps $\\$

is illustrated below:

We have thus demonstrated that ODEs of *cell maps* are included in the new setting of ODEs of chain maps. What about ODEs of *forms*?

Suppose we are given $f$, a $0$-form on ${\mathbb R}$. Then we would like to interpret the pair $g=\{f,df\}$ as some chain map defined on $C({\mathbb R})$, the chain complex of time. What is the other chain complex $C$, the chain complex of space? Since these two forms take their values in ring $R$, we can choose $C$ to be the trivial combination of two copies of $R$: $$\partial=\operatorname{Id}:R \to R.$$ Below, we consider a more general setting of $k$-forms.

**Theorem.** Cochains are chain maps, in the following sense: for every $k$-cochain $f$ on $K$, there is a chain map from $C(K)$ to the chain complex $C$ with only one non-zero part, $\operatorname{Id}:C_{k+1}=R \to C_k=R$, as shown in the following commutative diagram:
$$
\newcommand{\ra}[1]{\!\!\!\!\!\!\!\xrightarrow{\quad#1\quad}\!\!\!\!\!}
\newcommand{\da}[1]{\left\downarrow{\scriptstyle#1}\vphantom{\displaystyle\int_0^1}\right.}
%
\begin{array}{r|ccccccccccc}
C(K):& ... & C_{k+2}(K) & \ra{\partial} & C_{k+1}(K) & \ra{\partial} & C_k(K) & \ra{\partial} & C_{k-1}(K) & ... \\
f: & \ & \ \da{0} & & \ \da{df} & & \ \da{f} & & \ \da{0}&\\
C: & ... & 0 & \ra{\partial=0} & R & \ra{\partial=\operatorname{Id}} & R & \ra{\partial=0} &0&...
\end{array}
$$

**Proof.** We need to prove the commutativity of each of these squares. We go diagonally in two ways and demonstrate that the result is the same. We use the duality $d=\partial^*$.

For the first square: $$df \partial =(\operatorname{Id}^{-1}f\partial)\partial =\operatorname{Id}^{-1}f0=0.$$ For the second square: $$f\partial =df=\operatorname{Id}df.$$ The third square (and the rest) is zero. $\blacksquare$

Let's review how do ODEs of forms generate ODEs of chain maps. Suppose $df=Pf$ is an ODE of forms with $P:R\to R$ some function. Then, for $K={\mathbb R}$ and $k=0$, the above diagram becomes: $$ \newcommand{\ra}[1]{\!\!\!\!\!\!\!\xrightarrow{\quad#1\quad}\!\!\!\!\!} \newcommand{\da}[1]{\left\downarrow{\scriptstyle#1}\vphantom{\displaystyle\int_0^1}\right.} % \begin{array}{r|cccccccccc} C({\mathbb R}):& ... &\ \ra{0} & C_1({\mathbb R}) & \ra{\partial} & C_0({\mathbb R}) & \ra{0} \\ f: & \ & & \ \da{df} & & \ \da{f} & & \ \\ C: & ... &\ \ra{0} & R & \ra{\operatorname{Id}} & R & \ra{0} \end{array} $$ Then $P$ becomes our chain field: $$W=P:C_0=R\to C_1=R.$$

**Proposition.** A solution of an ODE of forms is a solution of the corresponding ODE of chain maps.

**Exercise.** Prove the proposition.

Thus, the examples at the beginning of this section can be understood as ODEs of chain maps.

**Exercise.** Provide chain fields for the examples of ODEs in this section.

**Exercise.** (a) Define the forward propagation map for ODEs of chain maps. (b) Explain its continuity.

## 6 Simulating advection with a spreadsheet

We have shown that we can model a flow with a *discrete vector field* that provides the velocities of the flow:

We have also modeled motion with cell maps and chain maps:

The last example of an ODE of chain maps, however, gives us something new.

**Example.** Suppose our chain field on ${\mathbb R}^2$ is given by
$$P(X):=\tfrac{1}{2}(X,XY)+\tfrac{1}{2}(X,XZ).$$
with $Y=(1,0),Z=(0,1)$. Where do we go from $X=(0,0)$? This is the result:
$$X\mapsto X+\partial P(X) = X+ \tfrac{1}{2}(Y-X)+\tfrac{1}{2}(Z-X)=\tfrac{1}{2}Y+\tfrac{1}{2}Z.$$
We can interpret this outcome as if some material initially located at $X$ has been split between the two adjacent edges and these halves ended up in $Y$ and $Z$. $\square$

A process that transfers material along a vector field is called *advection*.

**Example.** To model advection, we can use the vector field on the spreadsheet given above. It is understood differently though. If the values at the two adjacent to $X$ cells are $1$ and $1$, it used to mean that the particle follows vector $(1,1)$. Now it means that the material is split and each half goes along one of the two directions. Several iterations of this process are shown with the new cells with non-zero values shown in blue:

That's what a forward propagation map of a chain ODE looks like!

The coloring in the illustration does not reflect the fact that the material is spread thinner and thinner with time. That's why the simulation also resembles spreading of an *infection* -- along preset directions -- with blue squares indicating newly infected individuals, as below:

$\square$

For a more direct interpretation of the spreadsheet, imagine that a liquid moves along the square grid which is thought of as a *system of pipes*:

Then, instead of the velocity vector attached to each cell, we think of our vector field as two numbers each of which represents the amount of liquid that follows that pipe. (Then this is just a cubical $1$-form.)

In the picture above, the numbers represent “flows” through these “pipes”, with the direction determined by the direction of the $x,y$-axes. In particular,

- “$1$” means “$1$ cubic foot per second from left to right”.
- “$2$” means “$2$ cubic feet per second upward”, etc. $\\$

For the sake of *conservation of matter*, these numbers have to be normalized. Then $1/3$ of the amount goes right and $2/3$ goes up.

To confirm that this model makes sense, let's see what happens if a drop of dye be taken by such a flow. For computation, we follow the rule that

- the amount of dye in each cell is split and passed to the adjacent cells along the vectors of the vector field proportional to the values attached to the edges. $\\$

For instance, we can choose $1$'s on the horizontal edges and $0$s on the vertical edges, i.e., $dx$. Then the flow will be purely horizontal. If we choose $dy$, it will be purely vertical.

**Example.** Let's choose the $1$-form with $1$'s on both vertical and horizontal edges, i.e., $dx+dy$:
$$
\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}{cccccccccccccccccc}
& &\ \da{1} & & \ \da{1} & & \ \da{1} \\
& \ra{1} & \bullet & \ra{1} & \bullet & \ra{1} & \bullet & \ra{1}\\
& &\ \da{1} & &\ \da{1} & & \ \da{1} \\
& \ra{1} & \bullet & \ra{1} & \bullet & \ra{1} & \bullet & \ra{1}\\
& &\ \da{1} & & \ \da{1} & & \ \da{1}
\end{array}
$$
It is simple then: the amount of dye in each cell is split in half and passed to the two adjacent cells along the vectors: down and right. How well does it model the flow? Even though dispersal is inevitable, the predominantly diagonal direction of the spreading of the dye is evident in this spreadsheet simulation:

$\square$

**Exercise.** Find the formula for the diagonal elements.

**Exercise.** Create such a spreadsheet and confirm the pattern. What happens if the flow takes horizontally twice as much material as vertically?

**Exercise.** Devise an advection simulation with a circular dispersal pattern.

**Exercise.** (a) Explain how an arbitrary directed graph gives rise to advection and present its ODE. (b) What if, in addition to carrying material around, the flow also deposits a proportion of it in every location?