No good ending can be expected in the absence of the right beginning.
— I Ching
Before we can make our predictions better, we need to be able to make predictions in the first place. For physical systems, this means we need a mathematical model of how the system works. Once we have a set of equations that describe the system, we can get approximate solutions to them and make predictions about how the system will behave. It doesn't matter if we're looking at how fluids flow, how electromagnetic waves propagate, or how a building will flex in the wind; the first step is turning it into math.
Our Principal Equation
Fortunately, a lot of problems can be described by one family of equations. Unfortunately, they're partial differential equations (PDEs), which means its rare that analytic solutions exists for all but the simplest problems. This family of equations has a few names but we will call it the unsteady non-linear advection-diffusion-reaction equation,
$$
\frac{\partial u}{\partial t} + \nabla\cdot\left[\vec{F}\left(u\right) - \kappa\left(u\right)\nabla u\right] + S\left(u, \grad u\right) - f\left(\vec{x}, t\right) = 0.
$$
This is an equation for $u\left(\vec{x}, t\right)$, where $\vec{x}$ is any point in the $domain$ of the PDE i.e. the space where the PDE is defined. For example, if we were looking at how water in a river was flowing, the domain would be the river we're interested in. It wouldn't make sense to try and solve "how the river is flowing" in any impermeable rock nearby. We will refer to the domain of the PDE with $\Omega$, and the boundaries of the domain as $\partial\Omega$.
By changing the definitions of $\vec{F}$, $\kappa$, $S$, and $f$ we can change the problem we're modeling. Any of these terms can be zero. Whether a term is present of not can have significant effects on the behavior of the solutions.
Let's take a look at this equation term-by-term:
- $\frac{\partial u}{\partial t}$ — This term represents that the value of $u$ at a specific point in our domain may change in time. It is the presence of this term which makes the equation unsteady.
- $\nabla\cdot\vec{F}\left(u\right)$ — This term is called the advective term, and usually represents the movement of $u$ about the domain by some bulk motion. Dye being carried downstream in a river would be represented by this term.
- $\nabla\cdot\left[\kappa\left(u\right)\nabla u\right]$ — This term is the diffusion term. It represents the change in $u$ due to gradients. The change in temperature as you're defrosting a Thanksgiving turkey would be represented by this term.
- $S\left(u, \grad u\right)$ — This is term is reaction term (although we will call it the source term). It represents changes in $u$ that are a function of $u$ (and its gradient) at a specific point in the domain. Chemical reactions would be represented by this term.
- $f\left(\vec{x}, t\right)$ — This is the forcing term, and represents changes to $u$ that are independent of its value.
We can also think of $u$ as a vector of multiple variables, or as a single variable. Thinking of $u$ representing a vector is convenient if there are multiple PDEs of interest on the same domain, especially if they're coupled. In this case, we can think of each variable as a component of $u$, and each equation contributing to the components of $\vec{F}$, $\kappa\left(u\right)$ etc. Since we are dealing with spatial vectors and solution vectors, for the dot product we will use $a \cdot b$ if $a$ and $b$ are spatial vectors, and $a^{T} b$ if they are solution vectors.
The unsteady non-linear advection-diffusion-reaction equation is so powerful because it can describe (almost) any conservation law. That is, if you have a system that can be thought of as balancing the net transport of a quantity into a region with the change in that quantity with with time, and any creation or destruction if it, then you can probably use this equation to describe it. This includes:
- Fluid flows
- The temperature of solids, or fluids
- The stress and strain in solids
- Pollutant concentrations
- Oil moving through porous rock
and many more.
The Strong Form and its Limitations
PDEs can generally be written in two forms: strong form, and weak form. We wrote the unsteady non-linear advection-diffusion-reaction equation above in strong form. The problem with using the strong form is that it actually prevents you from accessing certain solutions. To demonstrate this, lets think about the shape of a string under tension with fixed ends, and a distributed load applied to it. If the string is stationary, we can describe the deflection of the string ($u$) as,
$$
-\frac{d^{2} u}{dx^{2}} = \frac{p}{T},
$$
where $T$ is the tension of the string and $p$ is the distributed load (i.e. force per unit length) on the string. We will assume that the string is fixed at $x = 0$ and $x = L$. Let's imagine we apply a distributed load (say $P$) over the middle half of the string only. Anywhere where we're not applying $P$, the string will be straight because $\frac{d^{2} u}{dx^{2}} = 0$, and anywhere we are applying $P$ the string will curve. Now let's double the distributed load ($2P$), but only apply it over the middle quarter of the string — more of the string will be straight, and the curvature of the middle quarter will be larger. We can repeat this process, doubling the distributed load but halving the portion of the string we apply it over. Each time, more of the string will be straight, and the curvature of the string in the middle will increase.
What happens if we repeat this process and infinite number of times? We have an infinitely large distributed load over an infinitely small length of string. We now have $\frac{d^{2} u}{dx^{2}} = \infty$ at the point at the center of the string, and $\frac{d^{2} u}{dx^{2}} = 0$ everywhere else. Intuitively, we expect the string to have two straight section with a "kink" in the middle. So, how do we solve $\frac{d^{2} u}{dx^{2}} = \infty$ at that center point where the kink is? The answer is, we can't. To solve this PDE, we can only get solutions where $\frac{d^{2} u}{dx^{2}}$ is finite... so no strings with kinks. We call the set of functions where $\frac{d^{2} u}{dx^{2}}$ is finite $\mathcal{C}^{2}$ — $\mathcal{C}$ for continuous, and $2$ for second derivative. In this case, solutions to the strong form can only be functions in $\mathcal{C}^{2}$.
From Strong Form to Weak Form
It turns out, there is a way to solve these kind of problems, but we cannot use the strong form of the PDE. To see this, let's take our equation for the deflection of the string, and multiply it by a weight function $w$,
$$ - w \frac{d^{2} u}{dx^{2}} = w \frac{p}{T}, $$
where $w$ is any function that fulfills some properties. We don't know what these properties are yet, so we'll use $\mathcal{X}$ as short hand for "all functions that fulfill some properties we'll define later". Out of convenience, we can also use the following set notation: $w \in \mathcal{X}$, which is short hand for "$w$ is any function from $\mathcal{X}$".
Our next steps are to integrate this equation over the entire domain of the PDE,
$$
\int_{0}^{L} -w \frac{d^{2} u}{dx^{2}} dx = \int_{0}^{L} w\frac{p}{T} dx,
$$
and apply integration-by-parts,
$$
\int_{0}^{L} \frac{d w}{d x} \frac{d u}{dx} dx - \left[w\left(L\right)\frac{d u\left(L\right)}{d x} - w\left(0\right)\frac{d u\left(0\right)}{d x}\right] = \int_{0}^{L} w\frac{p}{T} dx.
$$
We've got rid of the second derivative, and replaced it with a derivative on $w$ instead — this is the weak form!
Lets go one step further and say that $u$ and $w$ should be functions from the same set of functions, $\mathcal{X}$. In that case, we can use the weak form to find any solution where $\int_{0}^{L} \left(\frac{d u}{dx}\right)^{2} dx$ is finite. It turns out many more functions fulfill this criteria than have finite $\frac{d^{2} u}{dx^{2}}$. The space of functions with finite $\int_{0}^{L} \left(\frac{d u}{dx}\right)^{2} dx$ is called $\mathcal{H}^{1}$ — $\mathcal{H}$ means Hilbert (after the mathematician David Hilbert), and $1$ for the first derivative. For completeness, it is worth noting that in the problem described above $P$ becomes the delta function, so the integral on the right-hand-side of our weak form is also well defined.
Closing Remarks
The strong form of the unsteady non-linear advection-diffusion-reaction equation is given by,
$$
\mathcal{L}\left(u\right) \equiv \frac{\partial u}{\partial t} + \nabla\cdot\left[\vec{F}\left(u\right) - \kappa\left(u\right)\nabla u\right] + S\left(u, \grad u\right) - f\left(\vec{x}, t\right).
$$
Solutions to the strong form are found by solving,
$$
\mathcal{L}\left(u\right) = 0
$$
for $u$. This form of the PDE requires that $u$ has finite second derivitives over the entire domain.
The weak form of the unsteady non-linear advection-diffusion-reaction equation is given by,
$$
\begin{eqnarray*}
\mathcal{R}\left(u, w\right) \equiv \int_{\Omega} \Bigg[-\left(\frac{\partial w}{\partial t}\right)^{T} u - \left(\nabla w\right)^{T} \cdot \left[\vec{F}\left(u\right) - \kappa\left(u\right)\nabla u\right] + w^{T} S\left(u, \grad u\right) \\\ - w^{T} f\left(\vec{x}, t\right)\Bigg] dx + \int_{\partial\Omega} w \left[\vec{F}\left(u\right) - \kappa\left(u\right)\nabla u\right] ds.
\end{eqnarray*}
$$
Solutions to the weak form involves finding $u \in \mathcal{X}$, such that
$$
\mathcal{R}\left(u, w\right) = 0
$$
$\forall w \in \mathcal{X}$ (that is, for all $w$ in $\mathcal{X}$). Here, we've reduced (weakened as it were) the continuity requirements for solutions when compared to the strong form.
The weak form has advantages beyond the reduced continuity requirements. It is also much more amenable to mathematical analysis — results from functional analysis can be used to prove several important results for weak solutions.
We will be making significant use of the weak form of the unsteady non-linear advection-diffusion-reaction equation over the rest of this series.
Any questions or comments? Feel free to reach out to me on LinkedIn or Twitter and let me know!