Diffusion

The heat equation and diffusion processes.

In this chapter 4

Tick a unit to add it to your PDF basket.

The heat equation

The heat equation

Theorem

The heat kernel

Adapted

Add to PDF

Convolving bounded continuous initial data with the Gaussian kernel (4 pi k t)^{-1/2} exp(-x^2/4kt) solves the heat equation on the line, reproduces the data as t tends to 0, makes the solution positive everywhere at once, and bounds its m-th derivative by t^{-m/2}.

Prerequisites: The convolution theorem, Gaussian integrals, Differentiation under the integral sign

Definition(The heat equation and its kernel)AdaptedAdapted from Strauss (2008), §2.4, Eq. (1) p. 46 and Eq. (7) p. 49; Evans (2010), §2.3.1, p. 44

Fix a constant $k > 0$ , the diffusivity. The heat equation (or diffusion equation) on the line is

u_t = k\,u_{xx}, \qquad x \in \mathbb{R},\ t > 0 .

Its heat kernel is the function

S(x,t) = \frac{1}{\sqrt{4\pi k t}}\, e^{-x^2/4kt}, \qquad x \in \mathbb{R},\ t > 0 ,

extended by $S(x,t) = 0$ for $t < 0$ . For each fixed $t > 0$ , $S(\cdot,t)$ is the Gaussian probability density of mean $0$ and variance $2kt$ .

The same function is called the source function, the propagator, the diffusion kernel, the Gaussian kernel, and the fundamental solution of the heat equation; (Strauss (2008), p. 50) lists these names together. Evans (Evans (2010), p. 44) takes $k = 1$ and works in $\mathbb{R}^n$ , where the kernel is $(4\pi t)^{-n/2}e^{-|x|^2/4t}$ ; the $n = 1$ case of that formula is $S$ with $k = 1$ . The normalising constant is fixed by the requirement $\int_{\mathbb{R}} S(x,t)\,dx = 1$ .

Theorem(Solution of the initial-value problem)AdaptedAdapted from Evans (2010), §2.3.1, Thm 1, p. 45 and pp. 46–47; Strauss (2008), §2.4, Eq. (8), p. 49

Let $k > 0$ and let $g : \mathbb{R} \to \mathbb{R}$ be continuous and bounded. Define

u(x,t) = \int_{-\infty}^{\infty} S(x - y, t)\, g(y)\,dy = \frac{1}{\sqrt{4\pi k t}} \int_{-\infty}^{\infty} e^{-(x-y)^2/4kt}\, g(y)\,dy , \qquad x \in \mathbb{R},\ t > 0 .

Then:

Smoothness. $u \in C^{\infty}\bigl(\mathbb{R} \times (0,\infty)\bigr)$ .
Equation. $u_t(x,t) = k\,u_{xx}(x,t)$ for all $x \in \mathbb{R}$ and all $t > 0$ .
Initial data. For every $x_0 \in \mathbb{R}$ , $\lim_{\substack{(x,t) \to (x_0, 0)\\ x \in \mathbb{R},\ t > 0}} u(x,t) = g(x_0) .$ Hence setting $u(x,0) = g(x)$ makes $u$ continuous on $\mathbb{R} \times [0,\infty)$ .
Averaging and infinite propagation speed. $\inf_{\mathbb{R}} g \le u(x,t) \le \sup_{\mathbb{R}} g$ for all $x$ and all $t > 0$ . If moreover $g \ge 0$ and $g \not\equiv 0$ , then $u(x,t) > 0$ for every $x \in \mathbb{R}$ and every $t > 0$ , however small.

Part 4 is the statement that the heat equation propagates disturbances at infinite speed: data supported in a small interval is felt everywhere immediately.

Idea of the proofOriginal

Everything comes from three properties of the Gaussian: it solves the equation, it has total integral $1$ , and for small $t$ almost all of that integral sits in a tiny interval around the origin. The first two make $u$ a weighted average of the initial values which satisfies the equation, once differentiation under the integral sign is justified; the third makes the weight concentrate at $x_0$ , so the average tends to $g(x_0)$ . The one quantitative step is to split the integral into the part where $g$ is within $\varepsilon$ of $g(x_0)$ and the rest, and to show that the Gaussian gives the rest almost no weight.

ProofAdaptedAdapted from Evans (2010), Thm 1, pp. 45–46

Throughout, $\|g\|_\infty = \sup_{\mathbb{R}}|g|$ , which is finite because $g$ is bounded.

Step 1: $S$ solves the equation for $t > 0$ . Differentiating the explicit formula,

S_t = \Bigl(\frac{x^2}{4kt^2} - \frac{1}{2t}\Bigr) S, \qquad S_x = -\frac{x}{2kt}\,S, \qquad S_{xx} = \Bigl(\frac{x^2}{4k^2t^2} - \frac{1}{2kt}\Bigr) S ,

so $k\,S_{xx} = S_t$ . By translation invariance of the equation in $x$ , the function $(x,t) \mapsto S(x - y, t)$ also solves it for each fixed $y$ .

Step 2: unit mass. Substituting $x = z\sqrt{4kt}$ for fixed $t > 0$ ,

\int_{-\infty}^{\infty} S(x,t)\,dx = \frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty} e^{-z^2}\,dz = 1 ,

by the Gaussian integral $\int e^{-z^2}dz = \sqrt{\pi}$ . In particular the integral defining $u$ converges absolutely, with $|u| \le \|g\|_\infty$ .

Step 3: smoothness and the equation (parts 1 and 2). Fix $R > 0$ and $0 < \delta < T$ , and let $Q = [-R, R] \times [\delta, T]$ . Each partial derivative $\partial_t^{\,j}\partial_x^{\,m} S(x,t)$ is $P_{j,m}(x, 1/t)\,S(x,t)$ for a polynomial $P_{j,m}$ , as Step 1 illustrates; hence on $Q$ there is a constant $C$ with

\bigl|\partial_t^{\,j}\partial_x^{\,m} S(x - y, t)\bigr| \le C\,\bigl(1 + |x-y|^{2j+m}\bigr)\, e^{-(x-y)^2/4kT} , \qquad (x,t) \in Q,\ y \in \mathbb{R} ,

because $t \ge \delta$ bounds the negative powers of $t$ and $t \le T$ bounds the Gaussian from above. For $|y| \ge 2R$ we have $|x - y| \ge |y| - R \ge |y|/2$ , so the right-hand side is at most $C(1 + |y|^{2j+m})e^{-y^2/16kT}$ ; for $|y| \le 2R$ it is at most a constant. In both cases the bound is independent of $(x,t) \in Q$ and integrable in $y$ , and so is its product with the bounded function $g$ . By differentiation under the integral sign, $u$ has partial derivatives of every order on the interior of $Q$ , obtained by differentiating $S$ inside the integral. Since $R$ , $\delta$ , $T$ were arbitrary, $u \in C^\infty(\mathbb{R}\times(0,\infty))$ , which is part 1, and

u_t - k\,u_{xx} = \int_{-\infty}^{\infty} \bigl(S_t - k\,S_{xx}\bigr)(x-y,t)\,g(y)\,dy = 0

by Step 1, which is part 2.

Step 4: the initial data (part 3). Fix $x_0$ and $\varepsilon > 0$ . Since $g$ is continuous at $x_0$ , choose $\rho > 0$ with

|g(y) - g(x_0)| < \varepsilon \qquad \text{whenever } |y - x_0| < \rho .

Let $|x - x_0| \le \rho/2$ and $t > 0$ . By Step 2 we may subtract the constant $g(x_0)$ inside the integral:

u(x,t) - g(x_0) = \int_{-\infty}^{\infty} S(x-y,t)\,\bigl(g(y) - g(x_0)\bigr)\,dy .

Split the integral at $|y - x_0| = \rho$ and bound the two pieces separately.

Near piece. On $|y - x_0| < \rho$ the integrand is at most $\varepsilon\,S(x-y,t)$ in absolute value, by the choice of $\rho$ , and $S > 0$ , so by Step 2 this piece is at most $\varepsilon$ in absolute value.
Far piece. On $|y - x_0| \ge \rho$ we use boundedness of $g$ : the integrand is at most $2\|g\|_\infty S(x-y,t)$ . Moreover $|y - x| \ge |y - x_0| - \rho/2 \ge \tfrac12|y-x_0|$ , so $S(x-y,t) \le (4\pi kt)^{-1/2}e^{-(y-x_0)^2/16kt}$ . Substituting $y - x_0 = z\sqrt{16kt}$ , $\int_{|y-x_0|\ge\rho} S(x-y,t)\,dy \;\le\; \frac{2}{\sqrt{\pi}} \int_{|z| \ge \rho/(4\sqrt{kt})} e^{-z^2}\,dz ,$ and the right-hand side tends to $0$ as $t \to 0^{+}$ , because the domain of integration shrinks to the empty set while the integrand is fixed and integrable. So there is $\tau > 0$ , depending on $\rho$ and $\varepsilon$ but not on $x$ , with far piece $< \varepsilon$ for $0 < t < \tau$ .

Hence $|u(x,t) - g(x_0)| < 2\varepsilon$ whenever $|x - x_0| \le \rho/2$ and $0 < t < \tau$ , which is part 3.

Step 5: averaging and positivity (part 4). Since $S(\cdot,t) > 0$ and integrates to $1$ , the number $u(x,t)$ is a weighted average of the values of $g$ , so it lies between $\inf g$ and $\sup g$ . Suppose now $g \ge 0$ and $g(y_0) > 0$ for some $y_0$ . By continuity there are $c > 0$ and an interval $[a,b] \ni y_0$ of positive length with $g \ge c$ on $[a,b]$ . Then, for any $x$ and any $t > 0$ ,

u(x,t) \ \ge\ c \int_a^b S(x - y, t)\,dy \ > \ 0 ,

because $S(x-y,t) > 0$ for every $y$ . $\square$

Theorem(Smoothing)Original

Let $g$ be continuous and bounded and let $u$ be as above. For every integer $m \ge 1$ and every $t > 0$ ,

\sup_{x \in \mathbb{R}} \bigl|\partial_x^m u(x,t)\bigr| \ \le\ \gamma_m\, \|g\|_\infty\, (2kt)^{-m/2}, \qquad \gamma_m = \int_{-\infty}^{\infty} \bigl|\varphi^{(m)}(z)\bigr|\,dz ,

where $\varphi(z) = (2\pi)^{-1/2}e^{-z^2/2}$ is the standard Gaussian density. The first two constants are

\gamma_1 = \sqrt{\tfrac{2}{\pi}} = 0.797884\ldots, \qquad \gamma_2 = 4\varphi(1) = \sqrt{\tfrac{8}{\pi e}} = 0.967882\ldots ,

so that $|u_x| \le \|g\|_\infty/\sqrt{\pi k t}$ . No smoothness whatever is assumed of $g$ : the bound involves only its supremum. It degrades as $t \to 0^+$ and improves without limit as $t$ grows.

ProofOriginal

Write $\sigma = \sqrt{2kt}$ , so that $S(x,t) = \sigma^{-1}\varphi(x/\sigma)$ by the definition of $S$ . Differentiating $m$ times in $x$ ,

\partial_x^m S(x,t) = \sigma^{-(m+1)}\,\varphi^{(m)}(x/\sigma) ,

and substituting $x = \sigma z$ gives $\int |\partial_x^m S(x,t)|\,dx = \sigma^{-m}\gamma_m$ , which is finite because $\varphi^{(m)}$ is a polynomial times $\varphi$ . Step 3 of the previous proof allows the derivative to be taken inside the integral, so

\bigl|\partial_x^m u(x,t)\bigr| = \Bigl|\int \partial_x^m S(x-y,t)\,g(y)\,dy\Bigr| \le \|g\|_\infty \int \bigl|\partial_x^m S(x-y,t)\bigr|\,dy = \gamma_m\,\|g\|_\infty\,\sigma^{-m} ,

which is the claim. For the constants: $\varphi'(z) = -z\varphi(z)$ , so $\gamma_1 = \int |z|\varphi(z)\,dz = 2\int_0^\infty z\varphi(z)\,dz = 2\varphi(0) = \sqrt{2/\pi}$ . Next $\varphi''(z) = (z^2-1)\varphi(z)$ , so $|z^2-1|\varphi$ equals $\varphi''$ on $|z| > 1$ and $-\varphi''$ on $|z| < 1$ ; integrating each of the three pieces with $\varphi'' = (\varphi')'$ ,

\gamma_2 = \int_{-\infty}^{\infty}|z^2-1|\varphi(z)\,dz = \bigl[\varphi'\bigr]_{-\infty}^{-1} - \bigl[\varphi'\bigr]_{-1}^{1} + \bigl[\varphi'\bigr]_{1}^{\infty} = 4\varphi(1) ,

using $\varphi'(\pm1) = \mp\varphi(1)$ and $\varphi'(\pm\infty) = 0$ . $\square$

Remark(Three ways to find the kernel)AdaptedAdapted from Strauss (2008), §2.4, pp. 47–49; Evans (2010), §2.3.1, pp. 43–44; Stein & Shakarchi (2003), Ch. 5, §2.1, p. 146; Stakgold & Holst (2011), §2.5, Example 12, pp. 177–178

The theorem says the kernel works; it does not say how to find it. The books take three routes.

Dilation invariance. If $u$ solves $u_t = k u_{xx}$ , so does $u(\sqrt{a}\,x, at)$ for every $a > 0$ (Strauss (2008), property (e), p. 47). Strauss therefore first solves the equation with the step initial datum $Q(x,0) = 1$ for $x>0$ , $0$ for $x<0$ , which is itself dilation invariant, so $Q$ can depend on $x$ and $t$ only through $p = x/\sqrt{4kt}$ ; the equation becomes the ordinary differential equation $g'' + 2pg' = 0$ , whose solution is an error function, and $S = \partial Q/\partial x$ (Strauss (2008), Steps 1–4, pp. 47–49). Evans instead looks for a solution invariant under $u(x,t) \mapsto \lambda^{\alpha}u(\lambda^{\beta}x, \lambda t)$ , finds $\beta = 1/2$ and $\alpha = n/2$ forced, and reduces to a radial ordinary differential equation (Evans (2010), pp. 43–44).
Fourier transform. With the convention $\hat f(\xi) = \int f(x)e^{-2\pi i x\xi}\,dx$ used by Stein and Shakarchi (Stein & Shakarchi (2003), p. 146), transforming the equation in $x$ turns it into $\partial_t \hat u = -4\pi^2 k\xi^2 \hat u$ for each fixed $\xi$ , so $\hat u(\xi,t) = \hat g(\xi)\,e^{-4\pi^2 k \xi^2 t}$ . The transform of a Gaussian is a Gaussian, so $e^{-4\pi^2 k\xi^2 t}$ is itself the transform of a Gaussian in $x$ ; inverting it gives exactly $S(\cdot,t)$ , whence $u = g * S$ by the convolution theorem. Stakgold and Holst (Stakgold & Holst (2011), Example 12, pp. 177–178) run the same computation with the convention $\hat f(\omega)=\int f(x)e^{i\omega x}dx$ and with a delta source in time, obtaining $H(t)(4\pi t)^{-1/2}e^{-x^2/4t}$ for $k=1$ as the causal fundamental solution, the one vanishing for $t<0$ ; this is why $S$ is defined to be $0$ for $t < 0$ .
As a limit of data. Equivalently, $S(\cdot,t)$ is the solution whose initial datum is a unit point mass: $S \to \delta_0$ as $t \to 0^+$ . Both Evans and Stakgold–Holst state the problem in this form, the first by writing $\Phi = \delta_0$ at $t=0$ (Evans (2010), p. 46), the second by noting that $E$ solves the homogeneous equation with $\lim_{t\to0}E(x,t) = \delta(x)$ (Stakgold & Holst (2011), p. 178).

Remark(Where the hypotheses are used)Original

Boundedness of $g$ is used twice in Step 4: to make the far piece small, and (through Step 2) to make the integral converge at all. It can be relaxed, but not freely: against the Gaussian weight $e^{-(x-y)^2/4kt}$ , the datum $g(y) = e^{y^2}$ gives an integrand whose exponent has $y^2$ -coefficient $1 - 1/(4kt)$ , so the integral converges for $0 < t < 1/4k$ and diverges for every $t \ge 1/4k$ . A formula that stops existing after a finite time is no solution, which is why boundedness (or at least sub-Gaussian growth) is assumed from the start.
Continuity of $g$ only at $x_0$ is what Step 4 actually uses; boundedness is needed globally. So part 3 holds at every point of continuity of a bounded, piecewise-continuous datum, as the second example below shows. Parts 1 and 2, and the averaging bound $\inf g \le u \le \sup g$ of part 4, need no continuity at all — only boundedness and measurability. The strict positivity in part 4 does use continuity, at one point where $g > 0$ : for the bounded $g \ge 0$ that is $1$ at $y = 0$ and $0$ elsewhere, $g \not\equiv 0$ but $u \equiv 0$ .
$t > 0$ and $k > 0$ . $S$ is defined only for $t > 0$ : at $t = 0$ the formula is meaningless, as (Strauss (2008), p. 50) points out, and the solution is recovered there only as a limit. The sign of $k$ matters: replacing $t$ by $-t$ turns the Gaussian into $e^{+x^2/4kt}$ , which is not integrable, so the formula gives nothing backwards in time.
Positivity in part 4 cannot be upgraded to a bound on the speed of a signal. The value $u(x,t)$ for $|x|$ large is positive but exponentially small, of size $e^{-x^2/4kt}$ ; see the numbers in the second example. The wave equation behaves differently and does have a finite propagation speed (d'Alembert's formula).

NoteAdaptedAdapted from Stein & Shakarchi (2003), Ch. 5, Thm 2.3, p. 148

The theorem produces a solution; it does not say that there is only one. Uniqueness needs a restriction on the growth of $u$ . Stein and Shakarchi prove that if $u$ is continuous on $\mathbb{R}\times[0,\infty)$ , solves the heat equation for $t>0$ , vanishes at $t=0$ , and $u(\cdot,t)$ lies in the Schwartz class uniformly for $0<t<T$ in the sense that $\sup_{x,\,0<t<T}|x|^{j}|\partial_x^{\ell}u(x,t)| < \infty$ for all $j,\ell \ge 0$ , then $u \equiv 0$ ; the proof differentiates the energy $\int |u|^2dx$ and shows that it decreases. For merely bounded solutions, uniqueness follows instead from the maximum principle.

Example(Gaussian initial data spread as a Gaussian)Original

Take $g(x) = \theta_0\,e^{-x^2/2a^2}$ , a bell of height $\theta_0$ and width $a > 0$ . The convolution of two Gaussians is Gaussian, and the variances add: the datum has variance $a^2$ , the kernel variance $2kt$ , so

u(x,t) = \theta_0\,\frac{a}{\sqrt{a^2 + 2kt}}\;\exp\!\Bigl(-\frac{x^2}{2(a^2 + 2kt)}\Bigr) .

Check that this solves the problem. Write $w(t) = a^2 + 2kt$ , so $u = \theta_0 a\,w^{-1/2}e^{-x^2/2w}$ . Then $u_t = \theta_0 a\,e^{-x^2/2w}\bigl(k x^2 w^{-5/2} - k w^{-3/2}\bigr)$ and $u_{xx} = \bigl(x^2/w^2 - 1/w\bigr)u$ , so $u_t = k\,u_{xx}$ ; and $w(0) = a^2$ gives $u(x,0) = g(x)$ .

Numbers. With $\theta_0 = 1$ , $a = 1$ and $k = 1$ (any consistent units), the width grows as $\sqrt{1 + 2t}$ and the peak falls as $1/\sqrt{1+2t}$ . At $t = 1.5$ the variance is $1 + 3 = 4$ , so the peak has halved exactly: $u(0,1.5) = 1/2$ . The product (peak) $\times$ (width) is $\theta_0 a$ for all $t$ , which is Step 2 again: $\int u\,dx = \int g\,dx = \theta_0 a\sqrt{2\pi}$ is conserved.

Smoothing bound. Here $\|g\|_\infty = 1$ and $\sup_x|u_x| = \theta_0 a\,w^{-1}e^{-1/2}$ , attained at $x = \pm\sqrt{w}$ . At $t=1.5$ , $k=1$ this is $e^{-1/2}/4 = 0.1516$ , comfortably under the general bound $1/\sqrt{\pi k t} = 0.4607$ of the smoothing statement.

Example(A box of heat, and infinite propagation speed in numbers)Original

Let $g = \theta_0$ on $[-L, L]$ and $g = 0$ outside: a slab of hot material, initially at temperature $\theta_0$ , placed in a cold line. This $g$ is bounded and continuous except at $x = \pm L$ , so parts 1, 2 and 4 of the theorem apply as stated and part 3 applies at every $x_0 \ne \pm L$ , by the second item of the remark on hypotheses.

Substituting $z = (y-x)/\sqrt{4kt}$ in the solution formula and using $\operatorname{erf}(s) = \frac{2}{\sqrt\pi}\int_0^s e^{-z^2}dz$ ,

u(x,t) = \frac{\theta_0}{2}\left[ \operatorname{erf}\!\Bigl(\frac{L-x}{\sqrt{4kt}}\Bigr) + \operatorname{erf}\!\Bigl(\frac{L+x}{\sqrt{4kt}}\Bigr) \right] .

Take $\theta_0 = 1$ , $L = 1$ , $k = 1$ .

Centre. $u(0,1) = \operatorname{erf}(1/2) = 0.520500$ : after a time $L^2/k$ the centre has cooled by nearly half.
Outside the slab. $u(2,1) = \frac12[\operatorname{erf}(-1/2) + \operatorname{erf}(3/2)] = 0.222803$ .
Far outside. $u(10,1) = \frac12[\operatorname{erf}(-9/2)+\operatorname{erf}(11/2)] = 9.83 \times 10^{-11}$ . It is positive, as part 4 demands, and utterly negligible: this is what infinite propagation speed amounts to in practice.
At the jump. $u(L,t) = \frac12\bigl[\operatorname{erf}(0) + \operatorname{erf}(L/\sqrt{kt})\bigr] \to \tfrac12$ as $t \to 0^+$ , the average of the two one-sided limits of $g$ — not $g(L)$ , whatever value we assign there. Part 3 genuinely fails at a jump.
Smoothing. Differentiating, $u_x(x,t) = (4\pi kt)^{-1/2}\bigl[e^{-(L+x)^2/4kt} - e^{-(L-x)^2/4kt}\bigr]$ . At $t = 1$ its largest absolute value is $0.20604$ , at $x = \pm 1.54340$ , against the bound $\|g\|_\infty/\sqrt{\pi k t} = 0.56419$ . The initial datum was not even continuous; a moment later the solution is $C^\infty$ with a derivative bounded by $0.57$ .

ApplicationOriginal

The diffusion length. Every use of the kernel in the laboratory is really a use of its width: after a time $t$ , matter or heat initially at a point has spread over a distance of order $\sqrt{2kt}$ (one standard deviation) or $\sqrt{4kt}$ (the length appearing in the exponent). Turned around, the time to spread a distance $\ell$ is of order $\ell^2/k$ , quadratic in the distance, which is why diffusion is fast over micrometres and hopeless over metres.

Molecules in water. A small molecule in water has $D$ of order $10^{-5}\,\mathrm{cm^2/s}$ . Spreading one standard deviation of $1\,\mathrm{cm}$ takes $t = \ell^2/2D \approx 5\times10^4\,\mathrm{s}$ , about $14$ hours: unstirred solutions do not mix. Across a $10\,\mu\mathrm{m} = 10^{-3}\,\mathrm{cm}$ cell the same estimate gives $10^{-6}/(2\times10^{-5}) = 5\times10^{-2}\,\mathrm{s}$ , which is why diffusion suffices inside a cell and circulation is needed between cells.
Heat in metal. Copper has thermal diffusivity of order $1.1\,\mathrm{cm^2/s}$ . Heat applied to one face of a $1\,\mathrm{cm}$ plate is felt on the other in $\ell^2/4k \approx 0.23\,\mathrm{s}$ ; for a stainless steel of diffusivity $0.04\,\mathrm{cm^2/s}$ the same thickness takes about $6\,\mathrm{s}$ . This single number governs quench rates, soldering times and the thermal lag of furnace walls.
Probability. The kernel is the transition density of Brownian motion. Strauss states this as: the probability that a particle starting at $x$ lies in $(a,b)$ at time $t$ is $\int_a^b S(x-y,t)\,dy$ for some constant $k$ , so that a concentration profile and a probability density obey the same equation (Strauss (2008), p. 51). He does not fix the constant; matching the kernel's variance $2kt$ against the convention in which the motion has variance $\sigma^2 t$ gives $k = \sigma^2/2$ , half the diffusion coefficient in that normalisation. See Brownian motion and the Fokker-Planck equation.
Smoothing, and its reverse. The derivative bound $|\partial_x^m u| \le \gamma_m\|g\|_\infty(2kt)^{-m/2}$ says that diffusion erases detail at a rate set by $t^{-m/2}$ : a feature of width $\lambda$ in the data is invisible once $\sqrt{2kt} \gg \lambda$ . Recovering the data from a later profile therefore amplifies small errors without bound, which is the reason inverse heat-conduction problems (inferring a surface temperature history from an interior measurement) must be regularised.

The kernel is also the building block for the inhomogeneous problem and for bounded domains: convolving it in time against a source gives Duhamel's formula (Green's functions for the wave and heat equations), and reflecting it in a boundary gives the kernel of a half-line or an interval (diffusion in bounded domains).

Theorem

Maximum principle for the heat equation

Adapted

Add to PDF

A solution of the heat equation in a bounded space-time cylinder takes its largest and smallest values on the initial slice or the lateral boundary, which gives uniqueness and a quantitative stability bound for the initial-boundary-value problem without solving it.

Prerequisites: Classification of second-order equations, Types of boundary and initial condition

Definition(Parabolic cylinder and parabolic boundary)AdaptedAdapted from Evans (2010), §2.3.2, pp. 49–50; Stakgold & Holst (2011), §8.2, Fig. 8.8, p. 483

Let $U \subset \mathbb{R}^n$ be open and bounded and let $T > 0$ . The parabolic cylinder is

\Omega = U \times (0, T] ,

which contains its top face $U \times \{t = T\}$ but no part of the lateral surface. Its closure is $\overline{\Omega} = \overline{U}\times[0,T]$ , and the parabolic boundary is what the closure adds:

\Gamma = \overline{\Omega} \setminus \Omega = \bigl(\overline{U}\times\{t=0\}\bigr) \cup \bigl(\partial U \times (0,T]\bigr) ,

the bottom together with the lateral sides. Evans (Evans (2010), pp. 49–50) writes $U_T$ and $\Gamma_T$ for these and stresses that $U_T$ includes the top. For $n = 1$ and $U = (0,\ell)$ the picture is a rectangle in the $(x,t)$ plane whose parabolic boundary is its bottom edge and its two vertical edges, the top edge excluded; this is the setting of Strauss (Strauss (2008), §2.3, p. 42).

A function $u$ is a subsolution of the heat equation with diffusivity $k > 0$ on $\Omega$ if $u$ is continuous on $\overline{\Omega}$ , has continuous $u_t$ and second spatial derivatives in $\Omega$ (the derivative in $t$ at $t = T$ being the left derivative), and satisfies

u_t - k\,\Delta u \le 0 \qquad\text{in } \Omega .

Physically $u_t - k\Delta u$ is the rate at which heat is being supplied, so a subsolution is one with no interior source.

Theorem(Weak maximum principle)AdaptedAdapted from Stakgold & Holst (2011), §8.2, Thm 8.2.1, p. 484; Strauss (2008), §2.3, p. 42

Let $U \subset \mathbb{R}^n$ be open and bounded, $T > 0$ , $k > 0$ , and let $u$ be a subsolution of the heat equation on $\Omega = U\times(0,T]$ . If $u \le M$ on the parabolic boundary $\Gamma$ , then $u \le M$ on all of $\overline{\Omega}$ . Equivalently,

\max_{\overline{\Omega}}\, u = \max_{\Gamma}\, u .

Applying this to $-u$ : if $u_t - k\Delta u \ge 0$ in $\Omega$ and $u \ge m$ on $\Gamma$ , then $u \ge m$ on $\overline{\Omega}$ .

NoteAdaptedAdapted from Stakgold & Holst (2011), §8.2, Thm 8.2.1, p. 484; Strauss (2008), §2.3, p. 42; Evans (2010), §2.3.3, Thm 4, pp. 52–53

The three sources state this result with different generality, and the version above is the one followed from Stakgold and Holst (Stakgold & Holst (2011), Thm 8.2.1, p. 484): a differential inequality, in $n$ space dimensions, on a parabolic domain. Their constant is $k = 1$ ; the diffusivity is carried along here, which changes nothing in the argument. Strauss (Strauss (2008), p. 42) states the equation case in one dimension on a rectangle, with the diffusivity $k$ . Evans (Evans (2010), Thm 4, pp. 52–53) proves the sharper strong maximum principle, also for the equation: if a solution attains its maximum over $\overline{U_T}$ at a point $(x_0,t_0)$ of $U_T$ and $U$ is connected, then $u$ is constant on $\overline{U_{t_0}}$ , that is, at all earlier times. Strauss remarks that the strong version is considerably harder to prove and refers the reader to Protter and Weinberger (Strauss (2008), p. 43); Evans derives it from a mean-value property over "heat balls" (Evans (2010), Thm 3, p. 50). Nothing below uses the strong form.

Idea of the proofOriginal

At an interior maximum of a function of $x$ alone, calculus gives $\Delta u \le 0$ ; at an interior maximum in $t$ , $u_t = 0$ . Together these say $u_t - k\Delta u \ge 0$ , which just fails to contradict $u_t - k\Delta u \le 0$ : the inequality is not strict, so the maximum could sit inside. The repair is to add a small bump that makes the inequality strict, $v = u + \varepsilon|x|^2$ , whose Laplacian is larger than $u$ 's by the fixed amount $2n\varepsilon$ . For $v$ the contradiction is real, so $v$ has no interior maximum; then let $\varepsilon \to 0$ . The top face needs separate attention, because a maximum there is one-sided in $t$ and gives only $u_t \ge 0$ — which is still the wrong sign for the contradiction, so the argument survives.

ProofAdaptedAdapted from Stakgold & Holst (2011), §8.2, Lemma and Thm 8.2.1, pp. 484–485; Strauss (2008), §2.3, pp. 43–44

Step 1: a strict subsolution has no maximum in $\Omega$ . Suppose $v$ is continuous on $\overline{\Omega}$ , has the required derivatives in $\Omega$ , and satisfies the strict inequality $v_t - k\Delta v < 0$ there. Assume, for contradiction, that $v$ attains its maximum over the compact set $\overline{\Omega}$ at a point $P = (x_0, t_0) \in \Omega$ , so that $x_0 \in U$ and $0 < t_0 \le T$ .

Space. Since $U$ is open, $x_0$ is an interior maximum of the function $x \mapsto v(x, t_0)$ on $U$ . Hence for each $i$ the function of one variable $s \mapsto v(x_0 + s e_i, t_0)$ has an interior maximum at $s = 0$ , so its second derivative there is $\le 0$ : that is, $\partial^2 v/\partial x_i^2 \le 0$ at $P$ . Summing, $\Delta v(P) \le 0$ .
Time. If $t_0 < T$ , then $t \mapsto v(x_0, t)$ has an interior maximum at $t_0$ , so $v_t(P) = 0$ . If $t_0 = T$ , the maximum is one-sided: $v(x_0, T) \ge v(x_0, T - \delta)$ for small $\delta > 0$ , so the difference quotients $\bigl(v(x_0,T) - v(x_0, T-\delta)\bigr)/\delta$ are $\ge 0$ and their limit, the left derivative $v_t(P)$ , satisfies $v_t(P) \ge 0$ .

In both cases $v_t(P) \ge 0 \ge k\,\Delta v(P)$ , using $k > 0$ in the second inequality, so $v_t(P) - k\Delta v(P) \ge 0$ . This contradicts the strict inequality. Hence the maximum of $v$ over $\overline{\Omega}$ is attained only on $\Gamma$ .

Step 2: from $\le$ to $<$ . Let $u$ be as in the statement and let $\varepsilon > 0$ . Because $U$ is bounded, there is $a > 0$ with $|x| \le a$ for all $x \in \overline{U}$ . Put

v(x,t) = u(x,t) + \varepsilon |x|^2 .

Then $v_t = u_t$ , while $\Delta|x|^2 = 2n$ , so $\Delta v = \Delta u + 2n\varepsilon$ and

v_t - k\,\Delta v = \bigl(u_t - k\,\Delta u\bigr) - 2nk\varepsilon \le -2nk\varepsilon < 0 ,

where the first inequality is the hypothesis on $u$ and the second needs $k > 0$ and $\varepsilon > 0$ . So Step 1 applies to $v$ .

Step 3: conclusion. For any $(x,t) \in \overline{\Omega}$ ,

u(x,t) \ \le\ v(x,t) \ \le\ \max_{\Gamma} v \ \le\ \max_{\Gamma} u + \varepsilon a^2 \ \le\ M + \varepsilon a^2 ,

the first step because $\varepsilon|x|^2 \ge 0$ , the second by Step 2 and Step 1, the third because $|x| \le a$ on $\Gamma$ , the fourth by hypothesis. The left-hand side does not depend on $\varepsilon$ , so letting $\varepsilon \to 0^+$ gives $u(x,t) \le M$ . Since $\Gamma \subset \overline\Omega$ , taking $M = \max_\Gamma u$ gives the displayed equality of maxima. $\square$

Theorem(Comparison, uniqueness and stability)AdaptedAdapted from Stakgold & Holst (2011), §8.2, Thms 8.2.2–8.2.5, pp. 485–486; Strauss (2008), §2.3, uniqueness p. 44 and Eq. (6) p. 45

Let $U$ , $T$ , $k$ , $\Omega$ , $\Gamma$ be as above, and let all functions named be continuous on $\overline{\Omega}$ with continuous $u_t$ and second spatial derivatives in $\Omega$ .

Two-sided bound. If $u_t = k\Delta u$ in $\Omega$ and $m \le u \le M$ on $\Gamma$ , then $m \le u \le M$ throughout $\overline{\Omega}$ .
Comparison. If $u_t - k\Delta u \le v_t - k\Delta v$ in $\Omega$ and $u \le v$ on $\Gamma$ , then $u \le v$ throughout $\overline{\Omega}$ .
Uniqueness. For given $f$ on $\Omega$ and $g$ on $\Gamma$ there is at most one $u$ with $u_t - k\,\Delta u = f \ \text{ in } \Omega, \qquad u = g \ \text{ on } \Gamma .$
Stability. Let $u_1, u_2$ satisfy $\partial_t u_i - k\Delta u_i = f_i$ in $\Omega$ , and suppose $|f_1 - f_2| \le \alpha \ \text{ in } \Omega, \qquad |u_1 - u_2| \le \beta \ \text{ on } \Gamma .$ Then $|u_1 - u_2| \le \alpha T + \beta$ throughout $\overline{\Omega}$ . In particular, for the same source ( $\alpha = 0$ ) the solution never departs from its data by more than the data departed: $\max_{\overline\Omega}|u_1 - u_2| \le \max_{\Gamma}|u_1-u_2|$ .

ProofAdaptedAdapted from Stakgold & Holst (2011), §8.2, proofs of Thms 8.2.3–8.2.5, pp. 485–486

Each part reduces to the maximum principle applied to a suitable combination.

The principle gives $u \le M$ ; applied to $-u$ , which satisfies $(-u)_t - k\Delta(-u) = 0 \le 0$ and $-u \le -m$ on $\Gamma$ , it gives $-u \le -m$ , that is, $u \ge m$ .
Let $w = u - v$ . Then $w_t - k\Delta w \le 0$ in $\Omega$ by hypothesis, so $w$ is a subsolution, and $w \le 0$ on $\Gamma$ . The principle with $M = 0$ gives $w \le 0$ on $\overline\Omega$ .
If $u_1$ and $u_2$ both solve the problem, then each of the hypotheses of part 2 holds with $(u,v) = (u_1,u_2)$ and also with $(u,v) = (u_2,u_1)$ , since the two differential expressions are equal and the two boundary functions are equal. Hence $u_1 \le u_2$ and $u_2 \le u_1$ .
Let $w = u_1 - u_2$ , so that $w_t - k\Delta w = f_1 - f_2$ in $\Omega$ and $|w| \le \beta$ on $\Gamma$ . Put $z(x,t) = w(x,t) - \alpha t - \beta .$ Then $z_t - k\Delta z = (f_1 - f_2) - \alpha \le 0$ in $\Omega$ , using the bound on $f_1 - f_2$ ; and on $\Gamma$ we have $z \le \beta - \alpha t - \beta \le 0$ , using $|w|\le\beta$ and $\alpha, t \ge 0$ . By the maximum principle $z \le 0$ on $\overline\Omega$ , that is, $w \le \alpha t + \beta \le \alpha T + \beta$ . Repeating the argument with $-w$ in place of $w$ (which satisfies the same two bounds) gives $-w \le \alpha T + \beta$ . $\square$

Remark(Where the hypotheses are used)Original

$U$ bounded is used twice: to know that the continuous function $u$ attains a maximum on the compact set $\overline\Omega$ at all, and in Step 3 to bound $\varepsilon|x|^2$ by $\varepsilon a^2$ uniformly. On an unbounded domain the conclusion is false without a growth restriction; see the note below.
$k > 0$ enters in Step 1 (to turn $\Delta v \le 0$ into $k\Delta v \le 0$ ) and in Step 2 (to make $-2nk\varepsilon$ negative). It cannot be dropped: for the backward equation $u_t = -k u_{xx}$ on $0 \le x \le 1$ , the function $u(x,t) = e^{\pi^2 k t}\sin \pi x$ is a solution which vanishes on the two lateral sides and has $u \le 1$ at $t=0$ , yet $u(\tfrac12, T) = e^{\pi^2 k T}$ , which exceeds every bound as $T$ grows. The example also shows why $\Gamma$ must contain the bottom and not the top: for the backward equation the natural data sit on the top.
The top face belongs to $\Omega$ , not to $\Gamma$ . This is the only delicate point of Step 1, and it works because a one-sided maximum in $t$ gives $v_t \ge 0$ , which has the same sign as the two-sided conclusion $v_t = 0$ needs. Consequently the principle, applied for every $T$ , bounds the solution for all time by its initial and lateral data.
Only an inequality is assumed. Nothing in the proof uses $u_t = k\Delta u$ . This is what makes the comparison and stability statements cheap: each is obtained by feeding an explicitly built subsolution, such as $w - \alpha t - \beta$ , into the same theorem.
The conclusion does not say the maximum is attained only on $\Gamma$ . A solution with constant data is constant, and its maximum is attained everywhere. Stakgold and Holst (Stakgold & Holst (2011), Remark, p. 485) note that constant initial and boundary values up to the time of an interior extremum are the only way this can happen; that sharper statement is the strong maximum principle.

NoteAdaptedAdapted from Evans (2010), §2.3.3, Thm 6, pp. 55–56

On the whole space the principle needs a growth hypothesis. Evans proves that if $u$ is continuous on $\mathbb{R}^n\times[0,T]$ , solves the heat equation (with $k=1$ ) in $\mathbb{R}^n\times(0,T)$ , equals $g$ at $t=0$ , and satisfies

u(x,t) \le A\,e^{a|x|^2} \qquad (x\in\mathbb{R}^n,\ 0 \le t \le T)

for some constants $A, a > 0$ , then $\sup_{\mathbb{R}^n\times[0,T]} u = \sup_{\mathbb{R}^n} g$ ; uniqueness for the Cauchy problem follows for solutions obeying the two-sided version of that bound (Evans (2010), Thms 6 and 7, pp. 55–56). The proof compares $u$ with the explicit solutions $\mu(T+\varepsilon-t)^{-n/2}e^{|x-y|^2/4(T+\varepsilon-t)}$ on large balls, and the restriction $4aT < 1$ it needs is removed by stepping through short time intervals. The solution built in the heat kernel is bounded, hence covered.

Example(A decaying mode on top of a steady state)Original

Take $n = 1$ , $U = (0,1)$ , $k = 1$ , and

u(x,t) = x + \tfrac12\, e^{-\pi^2 t}\sin(\pi x) .

It is a solution. $u_t = -\tfrac{\pi^2}{2}e^{-\pi^2t}\sin\pi x$ and $u_{xx} = -\tfrac{\pi^2}{2}e^{-\pi^2 t}\sin \pi x$ , since $x$ is annihilated by $\partial_{xx}$ ; so $u_t = u_{xx}$ everywhere, and $u$ is continuous on $\overline\Omega$ for any $T$ .

Its data. On the lateral sides $u(0,t) = 0$ and $u(1,t) = 1$ ; on the bottom $u(x,0) = x + \tfrac12\sin\pi x$ .

The maximum over $\Gamma$ . The bottom profile has derivative $1 + \tfrac{\pi}{2}\cos\pi x$ , which vanishes where $\cos \pi x = -2/\pi$ , that is at

x_* = \tfrac1\pi\arccos(-2/\pi) = 0.7196679\ldots,

and there $\sin \pi x_* = \sqrt{1 - 4/\pi^2}$ , so the bottom maximum is

M = x_* + \tfrac12\sqrt{1 - 4/\pi^2} = 1.1052568\ldots

The lateral values $0$ and $1$ are smaller, so $\max_\Gamma u = M = 1.10526$ and $\min_\Gamma u = 0$ .

The theorem's prediction, checked. Part 1 says $0 \le u \le 1.10526$ in the whole rectangle, for every $T$ . Directly: $u \ge 0$ because both terms are non-negative on $[0,1]$ , and $\max_x u(x,t)$ decreases from $1.10526$ at $t = 0$ to $1.04163$ at $t = 0.02$ . Beyond $t_* = \pi^{-2}\ln(\pi/2) = 0.045755\ldots$ it equals exactly $1$ , attained at the end $x = 1$ : indeed $u_x = 1 + \tfrac{\pi}{2}e^{-\pi^2t}\cos \pi x \ge 1 - \tfrac{\pi}{2}e^{-\pi^2 t}$ , which is non-negative precisely for $t \ge t_*$ , so from then on $u$ increases in $x$ . The maximum over the closed rectangle is therefore attained only at $t = 0$ , as the strong principle would also predict. Note that $u$ is not monotone in $x$ at early times, and its maximum over $\Gamma$ is attained at an interior point of the bottom edge: the principle locates the maximum on the parabolic boundary, not at a corner.

Stability in numbers. Let $\tilde u$ solve the same equation with the initial profile changed by at most $\beta = 10^{-2}$ and the same lateral data. Part 4 with $\alpha = 0$ gives $|u - \tilde u| \le 10^{-2}$ for all $x$ and all $t$ , with no computation of $\tilde u$ . If instead a source of size at most $\alpha = 0.1$ is switched on, then over $T = 2$ the two solutions can differ by at most $\alpha T + \beta = 0.21$ .

What fails for the backward equation. The function $v(x,t) = e^{\pi^2 t}\sin \pi x$ solves $v_t = -v_{xx}$ with $v = 0$ on both lateral sides and $v \le 1$ at $t = 0$ ; but $v(\tfrac12, 0.1) = e^{0.1\pi^2} = 2.6831$ , nearly three times the bound its parabolic data would suggest. Running diffusion backwards is a different business, and the principle genuinely needs $k > 0$ .

ApplicationOriginal

Bounds before solutions. The principle is the cheapest tool in the subject: it bounds a solution using only the data, with no formula, no series and no mesh. If a reactor wall starts between $20^\circ$ and $300^\circ\mathrm{C}$ and its two faces are held in that range, the interior stays in that range for ever, whatever the geometry. If a diffusing concentration is non-negative initially and at the boundary, and there are no sinks, it stays non-negative — which is what makes a concentration a concentration. For a model with a source $f$ , part 4 converts an uncertainty $\alpha$ in the source and $\beta$ in the data into the guaranteed error bar $\alpha T + \beta$ on the answer, a model-validation statement that holds before any numerics.

Comparison with explicit barriers. Part 2 is used by choosing $v$ to be something one can write down. To show that a temperature in a cooled body decays at least exponentially, compare it with a multiple of the slowest decaying mode; to bound a reaction-diffusion solution, compare it with the solution of the ordinary differential equation obtained by deleting the Laplacian, which is a subsolution or supersolution according to the sign of the reaction term (see advection, diffusion and reaction).

Numerical schemes that cannot overshoot. Discretise $u_t = k u_{xx}$ on a grid $x_j = j\,\Delta x$ , $t_n = n\,\Delta t$ by the explicit scheme

U_j^{n+1} = \lambda\,U_{j-1}^{n} + (1 - 2\lambda)\,U_j^{n} + \lambda\,U_{j+1}^{n}, \qquad \lambda = \frac{k\,\Delta t}{\Delta x^2} .

The three weights sum to $1$ , and they are all non-negative exactly when $\lambda \le \tfrac12$ . In that case each new value is a weighted average of three old ones, so $\min_j U^n_j \le U^{n+1}_j \le \max_j U^n_j$ : the scheme obeys a discrete maximum principle and can never manufacture a temperature outside the range of the data. For $\lambda > \tfrac12$ the middle weight is negative and the computed solution oscillates with growing amplitude. This is the same threshold that von Neumann's Fourier analysis of the scheme produces, reached here without any Fourier analysis. For copper, $k \approx 1.1\,\mathrm{cm^2/s}$ , a $1\,\mathrm{mm}$ grid forces $\Delta t \le \Delta x^2/2k = 4.5\times10^{-3}\,\mathrm{s}$ . See finite differences and von Neumann stability.

Uniqueness as a licence. Part 3 is what allows the separation-of-variables series, the Green's-function integral and the finite-element approximation to be called "the" solution and compared with each other: any two of them that satisfy the same initial and boundary conditions must agree (see diffusion in bounded domains).

Method

Diffusion in bounded domains

Adapted

Add to PDF

How to solve the heat equation on an interval by separating variables: the eigenvalue problem set by the boundary condition, the exponentially decaying modes it produces, the rate at which the solution settles to its equilibrium, and two ways to handle boundary data that are not zero.

Prerequisites: Eigenfunction expansions, Maximum principle for the heat equation, Uniform convergence, Integrating and differentiating series term by term

Definition(The initial-boundary-value problem and its eigenvalue problem)AdaptedAdapted from Strauss (2008), §4.1, Eqs. (13)–(16), p. 87; Strauss (2008), §4.2, Eq. (1), p. 89

On an interval $0 < x < \ell$ with diffusivity $k > 0$ we consider

u_t = k\,u_{xx} \quad (0 < x < \ell,\ t > 0), \qquad u(x,0) = \varphi(x),

together with one boundary condition at each end. The three standard homogeneous choices are

Dirichlet (ends held at the reference value): $u(0,t) = u(\ell,t) = 0$ ;
Neumann (ends insulated, no flux): $u_x(0,t) = u_x(\ell,t) = 0$ ;
mixed: one of each, for instance $u(0,t) = u_x(\ell,t) = 0$ .

A separated solution is one of the form $u = X(x)\,T(t)$ . Substituting and dividing by $kXT$ gives

\frac{T'}{kT} = \frac{X''}{X} = -\lambda ,

a quantity that depends on neither $x$ nor $t$ and so is a constant. Hence $T' = -\lambda k T$ , so $T(t) = e^{-\lambda k t}$ up to a factor, and $X$ solves the eigenvalue problem

-X'' = \lambda X \quad \text{on } (0,\ell), \qquad \text{plus the boundary condition.}

A number $\lambda$ for which this has a solution $X \not\equiv 0$ is an eigenvalue and such an $X$ an eigenfunction; Strauss (Strauss (2008), p. 88) notes that physicists and engineers call the eigenfunctions normal modes, since they are the shapes that persist for all time, decaying as a whole. We write $\lambda_n$ , $X_n$ for them.

Method(The modes of an interval)AdaptedAdapted from Strauss (2008), §4.1, Eq. (8), p. 85 and the eigenvalue terminology, p. 88; Strauss (2008), §4.2, Eqs. (2)–(4), p. 90 and Eq. (10), p. 91

For $-X'' = \lambda X$ on $(0,\ell)$ :

condition	eigenvalues $\lambda_n$	eigenfunctions $X_n(x)$
$X(0)=X(\ell)=0$	$(n\pi/\ell)^2$ , $n = 1,2,\dots$	$\sin(n\pi x/\ell)$
$X'(0)=X'(\ell)=0$	$(n\pi/\ell)^2$ , $n = 0,1,2,\dots$	$\cos(n\pi x/\ell)$
$X(0)=X'(\ell)=0$	$\bigl((n+\tfrac12)\pi/\ell\bigr)^2$ , $n = 0,1,2,\dots$	$\sin\bigl((n+\tfrac12)\pi x/\ell\bigr)$

In each case these are all the eigenvalues, including complex ones. The Dirichlet eigenvalues are strictly positive; the Neumann list contains $\lambda_0 = 0$ , whose eigenfunction is any constant, and its smallest positive eigenvalue is $(\pi/\ell)^2$ . There are infinitely many of them, in contrast with an $N\times N$ matrix, which has at most $N$ (Strauss (2008), p. 88).

NoteAdaptedAdapted from Strauss (2008), §4.1, p. 88; Strauss (2008), §4.2, p. 90 and Eq. (10), p. 91 with Exercises 1–2, p. 92

Strauss proves the Dirichlet claim by hand: $\lambda = 0$ would force $X = C + Dx$ with $C = D = 0$ ; $\lambda = -\gamma^2 < 0$ would force $D\sinh\gamma\ell = 0$ with $\sinh\gamma\ell \ne 0$ ; and for complex $\lambda$ , writing $X = Ce^{\gamma x}+De^{-\gamma x}$ with $\gamma^2 = -\lambda$ , the two boundary conditions give $e^{2\gamma\ell}=1$ , hence $\gamma = n\pi i/\ell$ and $\lambda = n^2\pi^2/\ell^2$ (Strauss (2008), p. 88). For the Neumann problem he verifies the positive eigenvalues and $\lambda = 0$ directly and states that the exclusion of negative and complex $\lambda$ goes the same way (Strauss (2008), p. 90). For the mixed case he gives the eigenvalues $(n+\tfrac12)^2\pi^2/\ell^2$ and the eigenfunctions outright at Eq. (10) (Strauss (2008), p. 91), leaving only their derivation to Exercises 1–2 (Strauss (2008), p. 92). The general reason behind all three is that $-d^2/dx^2$ with any of these conditions is self-adjoint and non-negative, which is the subject of Sturm-Liouville theory.

Method(Separation of variables, homogeneous boundary data)AdaptedAdapted from Strauss (2008), §4.1, pp. 87–88; Strauss (2008), §4.2, Eqs. (5)–(6), p. 90

Find the modes. Solve $-X'' = \lambda X$ with the given boundary condition, obtaining $\lambda_n$ and $X_n$ from the table.
Attach the time factors. Each $X_n(x)\,e^{-\lambda_n k t}$ solves the equation and the boundary conditions. Write the candidate solution as the superposition $u(x,t) = \sum_n A_n\, e^{-\lambda_n k t}\, X_n(x) .$
Match the initial data. Setting $t = 0$ requires $\varphi = \sum_n A_n X_n$ , so the $A_n$ are the eigenfunction coefficients of $\varphi$ (eigenfunction expansions). For the interval these are Fourier sine or cosine coefficients: $A_n = \frac{2}{\ell}\int_0^{\ell}\varphi(x)\,\sin\frac{n\pi x}{\ell}\,dx \quad\text{(Dirichlet)}, \qquad A_n = \frac{2}{\ell}\int_0^{\ell}\varphi(x)\,\cos\frac{n\pi x}{\ell}\,dx \quad\text{(Neumann)},$ the Neumann series being written $u = \tfrac12 A_0 + \sum_{n\ge1}A_n e^{-\lambda_n kt}\cos(n\pi x/\ell)$ so that the same formula gives $A_0$ .
Read off the long-time behaviour. The slowest-decaying mode dominates. Under Dirichlet conditions every $\lambda_n > 0$ , so $u \to 0$ at the rate $e^{-\lambda_1 k t}$ with $\lambda_1 = (\pi/\ell)^2$ . Under Neumann conditions the $n=0$ mode does not decay, so $u \to \tfrac12 A_0 = \frac{1}{\ell}\int_0^\ell \varphi$ , the average of the initial data: nothing escapes through an insulated end, so the substance spreads out uniformly.
Check. Verify the series term by term against the equation and the boundary conditions, and bound the answer with the maximum principle, which also guarantees that no second solution continuous on the closed rectangle $[0,\ell]\times[0,T]$ takes the same initial and boundary values.

Method(The series is a solution, and how fast it decays)Original

Let $k>0$ , $\ell>0$ , let $\varphi$ be square-integrable on $(0,\ell)$ , let $A_n$ be its sine coefficients and $\lambda_n = (n\pi/\ell)^2$ , and put

u(x,t) = \sum_{n=1}^{\infty} A_n\,e^{-\lambda_n k t}\,\sin\frac{n\pi x}{\ell}, \qquad t > 0 .

Then:

The series, and the series obtained by differentiating it term by term any number of times in $x$ and $t$ , converge absolutely and uniformly on $[0,\ell]\times[\delta,\infty)$ for every $\delta > 0$ . Hence $u$ is infinitely differentiable on $[0,\ell]\times(0,\infty)$ .
$u_t = k\,u_{xx}$ for $0<x<\ell$ , $t>0$ , and $u(0,t) = u(\ell,t) = 0$ for every $t>0$ .
$\|u(\cdot,t) - \varphi\|_{L^2(0,\ell)} \to 0$ as $t \to 0^{+}$ .
$\|u(\cdot,t)\|_{L^2(0,\ell)} \le e^{-\lambda_1 k t}\,\|\varphi\|_{L^2(0,\ell)}$ for all $t>0$ , with $\lambda_1 = (\pi/\ell)^2$ . The relaxation time is therefore $\tau = \frac{1}{\lambda_1 k} = \frac{\ell^2}{\pi^2 k} .$

The same statements hold for the cosine series under Neumann conditions, with part 4 applied to $u - \tfrac12 A_0$ and $\varphi - \tfrac12 A_0$ , the first positive Neumann eigenvalue again being $(\pi/\ell)^2$ .

Idea of the proofOriginal

Everything rests on the factor $e^{-\lambda_n k t} = e^{-(n\pi/\ell)^2 k t}$ . For any fixed positive time it beats every power of $n$ , so differentiating the series as often as one likes, which multiplies the $n$ -th term by a power of $n$ , cannot stop it converging uniformly. That makes term-by-term differentiation legal for $t>0$ , and then each term does the work. At $t = 0$ the factor is $1$ and this argument collapses, which is why the initial condition is recovered only in the mean-square sense, from Parseval's identity rather than from any smoothness. The decay rate is Parseval again: every coefficient is multiplied by at most $e^{-\lambda_1 kt}$ .

ProofOriginal

Write $S_n(x) = \sin(n\pi x/\ell)$ and $c = k\delta$ .

Step 1 (the coefficients are bounded). By Parseval's identity for the sine system on $(0,\ell)$ (eigenfunction expansions),

\frac{\ell}{2}\sum_{n\ge1} A_n^2 = \int_0^\ell \varphi^2\,dx < \infty ,

so $A_n \to 0$ ; in particular $|A_n| \le C$ for all $n$ , with $C = \sup_n |A_n|$ .

Step 2 (uniform convergence of every derived series). Differentiating the $n$ -th term $j$ times in $t$ and $m$ times in $x$ multiplies it by $(-\lambda_n k)^{j}$ and by $(n\pi/\ell)^{m}$ up to replacing the sine by a cosine, so for $t \ge \delta$ its absolute value is at most

C\,(k\pi^2/\ell^2)^{j}\,(\pi/\ell)^{m}\; n^{2j+m}\, e^{-c\,n^2\pi^2/\ell^2} ,

using $|A_n| \le C$ from Step 1 and $e^{-\lambda_n k t} \le e^{-\lambda_n c}$ from $t \ge \delta$ . The numbers $n^{p}e^{-cn^2\pi^2/\ell^2}$ are summable for every fixed $p$ , because $n^{p}e^{-an^2} \le n^p e^{-an}$ for $n \ge 1$ and $a>0$ . So each derived series is dominated by a convergent series of constants, hence converges absolutely and uniformly on $[0,\ell]\times[\delta,\infty)$ (uniform convergence). The three hypotheses of the term-by-term differentiation theorem are therefore met in each variable separately on every closed interval inside that region — the series converges uniformly, each term is continuously differentiable, and the differentiated series converges uniformly — so the series may be differentiated term by term, once in each variable and then, by the same argument applied to the differentiated series, as often as one likes. Hence $u$ is $C^\infty$ on $[0,\ell]\times(\delta,\infty)$ ; as $\delta>0$ was arbitrary, part 1 follows.

Step 3 (equation and boundary conditions). By Step 2 we may compute

u_t - k\,u_{xx} = \sum_{n\ge1} A_n\Bigl[(-\lambda_n k) + k\lambda_n\Bigr] e^{-\lambda_n kt} S_n(x) = 0 ,

because $-S_n'' = \lambda_n S_n$ by the choice of $\lambda_n$ . Each $S_n$ vanishes at $x = 0$ and $x = \ell$ , and the series converges uniformly up to those endpoints for $t \ge \delta$ , so $u(0,t) = u(\ell,t) = 0$ . This is part 2.

Step 4 (the initial data, in mean square). By Parseval again, applied to the function $u(\cdot,t)-\varphi$ whose $n$ -th sine coefficient is $A_n(e^{-\lambda_n kt}-1)$ ,

\|u(\cdot,t)-\varphi\|_{L^2}^2 = \frac{\ell}{2}\sum_{n\ge1}A_n^2\bigl(1 - e^{-\lambda_n k t}\bigr)^2 .

Let $\varepsilon>0$ . Choose $N$ with $\frac{\ell}{2}\sum_{n>N}A_n^2 < \varepsilon/2$ , possible by Step 1; the tail of the sum above is then $< \varepsilon/2$ as well, since $0 \le 1-e^{-\lambda_n kt} \le 1$ . The remaining finite sum has at most $N$ terms, each tending to $0$ as $t\to0^+$ , so it is $<\varepsilon/2$ for $t$ small. This is part 3.

Step 5 (decay). Since $\lambda_1 \le \lambda_n$ for all $n$ and $k, t > 0$ , we have $e^{-\lambda_n kt} \le e^{-\lambda_1 kt}$ , so

\|u(\cdot,t)\|^2_{L^2} = \frac{\ell}{2}\sum_{n\ge1}A_n^2 e^{-2\lambda_n kt} \le e^{-2\lambda_1 kt}\,\frac{\ell}{2}\sum_{n\ge1}A_n^2 = e^{-2\lambda_1 kt}\,\|\varphi\|^2_{L^2} ,

which is part 4. For the Neumann case, replace the sine system by the cosine system; Parseval reads $\int_0^\ell f^2 = \ell(A_0/2)^2 + \frac{\ell}{2}\sum_{n\ge1}A_n^2$ , the $n = 0$ term is unaffected by the time factor $e^{-\lambda_0 kt} = 1$ , and the estimate above applies to the remaining terms, whose eigenvalues are all at least $(\pi/\ell)^2$ . $\square$

Remark(Where the hypotheses are used, and two warnings)Original

Only square-integrability of $\varphi$ is used, through Step 1. No continuity, and certainly no compatibility with the boundary conditions, is needed: the solution is smooth for every $t>0$ regardless, which is the smoothing already seen in the heat kernel.
$t \ge \delta > 0$ is essential in Step 2, and so part 1 says nothing at $t = 0$ . If $\varphi$ is continuous with $\varphi(0)=\varphi(\ell)=0$ and piecewise $C^1$ , its sine series converges uniformly and then $u \to \varphi$ uniformly as well; without such a hypothesis part 3 is all that is available.
$k>0$ again. Reversing the sign of $k$ turns every factor $e^{-\lambda_n kt}$ into $e^{+\lambda_n|k|t}$ , and the series diverges for every $t>0$ unless the coefficients decay faster than every exponential. Backwards diffusion on an interval is as ill-posed as on the line.
Uniqueness is not part of the construction. It comes from the maximum principle, which also gives the bound $\min\{0,\inf\varphi\} \le u \le \max\{0,\sup\varphi\}$ under Dirichlet conditions — a useful check on a truncated series, which can oscillate. Both need more than this statement supplies: the principle is about functions continuous on the closed rectangle $[0,\ell]\times[0,T]$ with finite data, so it applies when $\varphi$ is bounded with $\varphi(0)=\varphi(\ell)=0$ and $u$ extends continuously to $t=0$ , as in the previous bullet, and not to a merely square-integrable $\varphi$ , for which $\inf\varphi$ and $\sup\varphi$ may be infinite and $u$ need have no limit at $t=0$ at each point.
Term-by-term differentiation at the boundary. Strauss (Strauss (2008), §5.6, p. 147) makes this a cautionary tale: with inhomogeneous boundary data, inserting the sine series into the equation and differentiating twice term by term yields $u_n' + k\lambda_n u_n = 0$ , which cannot meet the boundary data at all, and he draws the moral that the differentiation step is the illegal one. The correct version is in the next method.

Method(Inhomogeneous boundary data)AdaptedAdapted from Strauss (2008), §5.6, Method of shifting the data, pp. 149–150; Strauss (2008), §5.6, Expansion method, Eqs. (2)–(10), pp. 147–148

Now let the ends be driven:

u_t = k\,u_{xx}, \qquad u(0,t) = h(t), \quad u(\ell,t) = j(t), \qquad u(x,0) = \varphi(x) .

No separated solution can satisfy these boundary conditions, so one of two routes is needed.

Route A: shift the data. Subtract any known function that meets the boundary conditions. With

\bar u(x,t) = \Bigl(1 - \frac{x}{\ell}\Bigr)h(t) + \frac{x}{\ell}\,j(t) ,

the difference $v = u - \bar u$ has homogeneous Dirichlet conditions, initial data $\varphi - \bar u(\cdot,0)$ , and the source term $-(\bar u_t - k\bar u_{xx}) = -\bar u_t$ on the right-hand side. When $h$ and $j$ are constants this is the best possible choice: $\bar u$ is then the steady state, $\bar u_t = 0$ , and $v$ solves the homogeneous problem of the previous method exactly. The physical reading is that the solution is the equilibrium profile plus a transient that dies at the rate $e^{-t/\tau}$ , $\tau = \ell^2/\pi^2 k$ .

Route B: expand, and keep the boundary terms. Expand the unknown itself, $u(x,t) = \sum_{n\ge1}u_n(t)\sin(n\pi x/\ell)$ with $u_n(t) = \frac2\ell\int_0^\ell u(x,t)\sin(n\pi x/\ell)dx$ , and differentiate under the integral sign rather than term by term. Then $u_n' = \frac2\ell\int_0^\ell u_t\, S_n\,dx = \frac{2k}{\ell}\int_0^\ell u_{xx}S_n\,dx$ , and integrating by parts twice, with $S_n(0)=S_n(\ell)=0$ and $S_n'(0) = n\pi/\ell$ , $S_n'(\ell) = (-1)^n n\pi/\ell$ :

\int_0^\ell u_{xx}S_n\,dx = -\bigl[u\,S_n'\bigr]_0^{\ell} - \lambda_n\int_0^\ell u\,S_n\,dx = \frac{n\pi}{\ell}\Bigl(h(t) - (-1)^n j(t)\Bigr) - \frac{\ell\lambda_n}{2}u_n(t) .

The boundary data enter here, and nowhere else. So each coefficient obeys a forced linear ordinary differential equation,

u_n'(t) + k\lambda_n u_n(t) = \frac{2n\pi k}{\ell^2}\Bigl(h(t) - (-1)^n j(t)\Bigr), \qquad u_n(0) = \frac{2}{\ell}\int_0^\ell \varphi\, S_n\,dx ,

whose solution by the integrating factor $e^{k\lambda_n t}$ is

u_n(t) = u_n(0)\,e^{-k\lambda_n t} + \frac{2n\pi k}{\ell^2}\int_0^t e^{-k\lambda_n(t-s)}\Bigl(h(s) - (-1)^n j(s)\Bigr)\,ds .

Route B needs no guess and handles time-dependent data directly; Route A is quicker whenever the data are steady or share a single time dependence, for instance $h = H\cos\omega t$ , $j = J\cos\omega t$ , when one looks for $\bar u = u_0(x)\cos\omega t$ (Strauss (2008), p. 150).

Example(A Dirichlet transient, with numbers)Original

Take $\ell = 1$ , $k = 1$ , Dirichlet ends, and $\varphi(x) = x(1-x)$ .

Coefficients. $A_n = 2\int_0^1 x(1-x)\sin(n\pi x)\,dx = \dfrac{4\bigl(1-(-1)^n\bigr)}{n^3\pi^3}$ , which is $8/(n^3\pi^3)$ for odd $n$ and $0$ for even $n$ — as it must be, since $\varphi$ is symmetric about $x = \tfrac12$ and the even modes are antisymmetric.
Solution. $u(x,t) = \displaystyle\sum_{n \text{ odd}} \frac{8}{n^3\pi^3}\,e^{-n^2\pi^2 t}\sin(n\pi x)$ .
Numbers at the midpoint. $u(\tfrac12,0) = \tfrac14$ . At $t = 0.1$ the series gives $u(\tfrac12,0.1) = 0.0961619$ , of which the first mode alone contributes $\frac{8}{\pi^3}e^{-\pi^2/10} = 0.0961632$ : after one relaxation time the higher modes have gone, the error from truncating after one term being $1.3\times10^{-6}$ . At $t = 0.5$ , $u(\tfrac12,0.5) = 0.00185559$ .
Decay bound, checked. Here $\|\varphi\|_{L^2} = 1/\sqrt{30} = 0.1825742$ and $\lambda_1 = \pi^2$ . Part 4 of the statement predicts $\|u(\cdot,0.1)\|_{L^2} \le e^{-\pi^2/10}/\sqrt{30} = 0.0680468$ ; the actual value from Parseval is $0.0679976$ , just inside, the small slack being the contribution of the third mode.
Relaxation time. $\tau = 1/\pi^2 = 0.101321$ . Doubling the length to $\ell = 2$ multiplies $\tau$ by four: a bar twice as long takes four times as long to cool, the same quadratic law as the diffusion length.
Maximum principle, checked. $0 \le u \le \tfrac14$ everywhere, consistent with the bound from the data.

Example(Insulated ends settle to the mean)Original

Same interval and diffusivity, but with insulated ends $u_x(0,t) = u_x(1,t) = 0$ and the initial ramp $\varphi(x) = x$ .

Cosine coefficients: $A_0 = 2\int_0^1 x\,dx = 1$ and $A_n = 2\int_0^1 x\cos(n\pi x)\,dx = 2\bigl((-1)^n-1\bigr)/(n^2\pi^2)$ , which is $-4/(n^2\pi^2)$ for odd $n$ and $0$ for even $n$ . Hence

u(x,t) = \frac12 - \sum_{n\text{ odd}}\frac{4}{n^2\pi^2}\,e^{-n^2\pi^2 t}\cos(n\pi x) .

The $n=0$ term is $\tfrac12 A_0 = \tfrac12 = \int_0^1 x\,dx$ : the average of the data, which is conserved because no flux leaves through an insulated end. Indeed $\frac{d}{dt}\int_0^1 u\,dx = k\int_0^1 u_{xx}dx = k[u_x]_0^1 = 0$ .

Numbers at the cold end. $u(0,t)$ rises from $0$ towards $\tfrac12$ : $u(0,0.05) = 0.252044$ , $u(0,0.5) = 0.497085$ . At the hot end $u(1,0.05) = 0.747956$ , and $u(0,t)+u(1,t) = 1$ exactly for every $t$ , by the antisymmetry of $\varphi - \tfrac12$ about $x = \tfrac12$ . The deviation from the mean decays like $e^{-\pi^2 t}$ , so after $t = 3\tau = 3/\pi^2 = 0.304$ it has fallen by a factor $e^{-3} \approx 0.05$ , which is what "settled" means in practice.

Example(A driven end: the steady state plus a transient)Original

Take $\ell = 1$ , $k = 1$ , $u(0,t) = 0$ , $u(1,t) = 1$ (a bar with one end suddenly clamped at temperature $1$ ), and $u(x,0) = 0$ . The boundary data are inhomogeneous but constant, so Route A applies with the steady state $\bar u(x) = x$ .

Shift. $v = u - x$ has $v(0,t)=v(1,t)=0$ and $v(x,0) = -x$ .
Expand. The sine coefficients of $-x$ on $(0,1)$ are $-2\int_0^1 x\sin(n\pi x)dx = 2(-1)^{n}/(n\pi)$ .
Solution. $u(x,t) = x + \frac{2}{\pi}\sum_{n=1}^{\infty}\frac{(-1)^n}{n}\,e^{-n^2\pi^2 t}\,\sin(n\pi x) .$ At $t=0$ the series is the sine series of $-x$ , so $u(x,0)=0$ , as required; for $t>0$ it decays and $u \to x$ .
Numbers at the midpoint. $u(\tfrac12,0.01) = 0.000407$ , $u(\tfrac12,0.05) = 0.113844$ , $u(\tfrac12,0.2) = 0.4115664$ . The last value is $\tfrac12 - \tfrac{2}{\pi}e^{-\pi^2/5} = 0.5 - 0.0884336$ to seven digits, again because only the first mode survives. The midpoint reaches $99\%$ of its final value $\tfrac12$ when $\frac{2}{\pi}e^{-\pi^2 t} < 0.005$ , that is $t > 0.4911$ , about five relaxation times.
Where the modes come from. Note that the $n$ -th coefficient decays only like $1/n$ , because $-x$ does not vanish at $x=1$ : the shifted initial data are incompatible with the boundary condition there. The series therefore converges slowly at $t = 0$ — and converges beautifully for every $t>0$ , which is exactly the content of part 1 of the statement.

ApplicationOriginal

Which number actually matters. For a bounded body the whole of the long-time behaviour is in the smallest eigenvalue: the solution relaxes as $e^{-t/\tau}$ with $\tau = 1/\lambda_1 k$ , and $\lambda_1$ is set by the geometry alone. For an interval of length $\ell$ with both ends held fixed, $\tau = \ell^2/\pi^2 k$ . A $2\,\mathrm{cm}$ copper rod ( $k \approx 1.1\,\mathrm{cm^2/s}$ ) has $\tau = 4/(\pi^2\cdot1.1) = 0.37\,\mathrm{s}$ ; a $2\,\mathrm{cm}$ layer of dry soil ( $k \approx 2\times10^{-3}\,\mathrm{cm^2/s}$ ) has $\tau = 200\,\mathrm{s}$ . Designers of thermal insulation, annealing schedules and reactor shields quote this one number.

Mode truncation is a numerical method. Because mode $n$ decays as $e^{-n^2 t/\tau}$ , the series is its own error estimate: after one relaxation time the second mode is down by $e^{-3} = 0.05$ relative to the first and the third by $e^{-8} = 3\times10^{-4}$ . Two or three terms then give several digits, as the examples show, which is why spectral methods are so efficient for diffusion and so delicate at $t = 0$ (see spectral methods).

Insulated versus open ends. The difference between the Dirichlet and Neumann tables is a single eigenvalue, $\lambda_0 = 0$ , and it changes the physics entirely: with open ends the contents of the interval drain away to zero; with insulated ends the total amount is conserved and the profile flattens to its average. A calorimeter is designed to realise the second case, a diffusion cell with large reservoirs at both ends the first.

The same answer as the Green's function. Stakgold and Holst (Stakgold & Holst (2011), Eq. (8.2.36a), p. 482) write the Dirichlet heat Green's function for $(0,\ell)$ both as an image sum of Gaussians and as the eigenfunction series $\frac{2}{\ell}\sum_{n\ge1}\sin\frac{n\pi x}{\ell}\sin\frac{n\pi \xi}{\ell}e^{-k\lambda_n t}$ , the two being linked by the theta-function identity. Convolving either against $\varphi$ reproduces the series of this unit, so the modal method and the Green's-function method (Green's functions for the wave and heat equations) are two faces of one formula; the Gaussian form is the practical one at very early times, the modal form at late times.

Beyond the interval. Steps 1 to 5 are unchanged in higher dimensions and in other coordinate systems: only the eigenvalue problem changes, giving Bessel functions on a disc and spherical harmonics with spherical Bessel functions on a ball (separation of variables). The decay law $e^{-\lambda_1 k t}$ and the identification of the limit with the $\lambda = 0$ mode survive verbatim.

Application

Advection, diffusion and reaction

Adapted

Add to PDF

Adding transport and reaction to the diffusion equation: the Peclet number that says which of advection and diffusion wins, the change to a moving frame that solves the linear problem exactly, and the travelling fronts that a bistable reaction term supports.

Prerequisites: The heat kernel, First-order equations and characteristics, Travelling-wave solutions

Definition(Advection, diffusion, reaction)AdaptedAdapted from Stakgold & Holst (2011), §0.1, Remark 1, Eq. (0.1.9), p. 6; Stakgold & Holst (2011), §0.3, Eq. (0.3.1), p. 12, Eq. (0.3.11), p. 14 and Eqs. (0.3.13)–(0.3.18), pp. 15–16; Strauss (2008), §1.3, p. 15

Let $u(x,t)$ be the concentration of a substance carried along a pipe. Conservation of mass on an arbitrary interval says $u_t = -J_x + p$ , where $J$ is the flux and $p$ the volume rate of production. Two mechanisms contribute to $J$ :

diffusion, by Fick's law $J_{\mathrm{diff}} = -D\,u_x$ with $D>0$ the diffusivity, as in (Strauss (2008), p. 15);
advection (also called convection), the bodily transport of the substance by a fluid moving with velocity $V$ : $J_{\mathrm{adv}} = V u$ . Stakgold and Holst (Stakgold & Holst (2011), Remark 1, p. 6) introduce exactly this term by replacing Fourier's law $J = -k\,\mathrm{grad}\,u$ with $J = -k\,\mathrm{grad}\,u + Cu\mathbf{v}$ for a medium moving with velocity $\mathbf{v}$ , and note that for an incompressible fluid, where $\mathrm{div}\,\mathbf{v} = 0$ , this adds a term proportional to $\mathbf{v}\cdot\mathrm{grad}\,u$ to the heat equation.

With $D$ and $V$ constant and $p = 0$ this gives the advection-diffusion equation

u_t + V\,u_x = D\,u_{xx} .

If instead the production term depends on the concentration itself, $p = f(u)$ , the result is a reaction-diffusion equation

u_t = D\,u_{xx} + f(u) .

Stakgold and Holst (Stakgold & Holst (2011), §0.3, Eq. (0.3.11), p. 14) derive such systems for a first-order irreversible reaction in a catalyst pellet, where $f$ carries the Arrhenius factor $e^{-B/u}$ of their Eq. (0.3.1) (Stakgold & Holst (2011), p. 12); later in the same section they build an energy-balance climate model of the same shape, with $f(u) = QS(x)a(u) - (A + Bu)$ (Stakgold & Holst (2011), Eqs. (0.3.13)–(0.3.18), pp. 15–16). Dropping $D$ leaves the first-order equation $u_t + Vu_x = 0$ of characteristics; dropping $V$ and $f$ leaves the heat equation of the heat kernel.

Definition(The Peclet number)Original

Fix a length $L$ characteristic of the problem — a pipe length, the width of an initial release, the size of a grid cell. Advection moves the substance across that length in the time $L/V$ ; diffusion spreads it across that length in the time $L^2/D$ . The ratio of these times is the dimensionless Peclet number

\mathrm{Pe} = \frac{L^2/D}{L/V} = \frac{V L}{D} .

Equivalently, scaling $x = L\xi$ and $t = (L^2/D)\tau$ in the advection-diffusion equation gives

u_\tau + \mathrm{Pe}\,u_\xi = u_{\xi\xi} ,

so $\mathrm{Pe}$ is the only parameter left: all advection-diffusion problems with the same Peclet number and the same shaped data are the same problem. Large $\mathrm{Pe}$ means transport dominates and the equation is nearly the first-order one; small $\mathrm{Pe}$ means diffusion dominates and the velocity is a small correction.

Application(Solution in a moving frame)Original

Let $D > 0$ and $V \in \mathbb{R}$ be constants, let $g$ be continuous and bounded on $\mathbb{R}$ , and let

S_D(x,t) = \frac{1}{\sqrt{4\pi D t}}\,e^{-x^2/4Dt}

be the heat kernel with diffusivity $D$ . Then

u(x,t) = \int_{-\infty}^{\infty} S_D\bigl(x - Vt - y,\, t\bigr)\, g(y)\,dy , \qquad t>0 ,

is infinitely differentiable for $t>0$ , solves $u_t + Vu_x = D u_{xx}$ , and tends to $g(x_0)$ as $(x,t)\to(x_0,0)$ with $t>0$ . In words: the profile diffuses exactly as if $V$ were zero, while the whole pattern is carried along at speed $V$ . Three consequences:

Mean and width. Suppose in addition that $g \ge 0$ , that its total mass $M = \int g$ is finite and positive, and that its second moment is finite, $\int y^2 g(y)\,dy < \infty$ . Then $\int u(x,t)\,dx = M$ for all $t$ , the centre of mass is at $\bar x(t) = \bar x(0) + Vt$ , and the variance about it is $\mathrm{var}(t) = \mathrm{var}(0) + 2Dt$ . (A finite first moment need not be assumed: Cauchy-Schwarz applied to $|y|\sqrt{g}\cdot\sqrt{g}$ gives $\int |y|\,g \le \bigl(\int y^2 g\bigr)^{1/2}\bigl(\int g\bigr)^{1/2} < \infty$ . Some decay hypothesis is needed, and the second moment is the stronger of the two: $g(y) = (1+y^2)^{-1}$ is bounded, continuous and non-negative with finite mass but with neither a mean nor a variance, while $g(y) = (1+|y|^3)^{-1}$ has a finite mass and a finite mean, the integrands there decaying like $|y|^{-3}$ and $|y|^{-2}$ , yet $\int y^2 g = \infty$ because there the integrand decays only like $|y|^{-1}$ .)
Relative spread. Keep the hypotheses of part 1 and let $V > 0$ . By the time the cloud has travelled a distance $L$ , that is at $t = L/V$ , diffusion has added $2DL/V$ to the variance, so the standard deviation is $\sqrt{\mathrm{var}(0) + 2DL/V} .$ For a release that starts concentrated, meaning $\mathrm{var}(0)$ negligible beside $2DL/V$ , this is close to $\sqrt{2DL/V} = L\sqrt{2/\mathrm{Pe}}$ with $\mathrm{Pe} = VL/D$ . (The hypotheses of part 1 force $\mathrm{var}(0) > 0$ strictly, since a continuous $g \ge 0$ of positive mass cannot be concentrated at a point; the limiting value $\mathrm{var}(0)=0$ belongs to the point release of the example below, which the statement does not cover but the kernel itself handles.) The cloud therefore stays sharp relative to the distance travelled exactly when $\mathrm{Pe}$ is large and the initial width is small compared with $L$ .
A linear decay term. If $u_t + Vu_x = Du_{xx} - r u$ with $r$ constant, then $u = e^{-rt}w$ where $w$ solves the equation without the decay term, so the formula above acquires the factor $e^{-rt}$ .

Idea of the proofOriginal

The equation treats $x$ and $t$ asymmetrically only through the drift, and a drift is exactly what a change of observer removes. Sitting on the fluid, that is substituting $y = x - Vt$ , the advection term and the time derivative combine into the material derivative, and what is left is the heat equation. Everything then follows from the heat kernel, read in the new coordinate.

ProofOriginal

Step 1 (change of frame). Suppose $w$ is a solution of $w_t = D w_{yy}$ and set $u(x,t) = w(x - Vt, t)$ . By the chain rule, with $y = x - Vt$ ,

u_t(x,t) = w_t(y,t) - V\,w_y(y,t), \qquad u_x = w_y, \qquad u_{xx} = w_{yy} ,

so $u_t + V u_x = w_t = D w_{yy} = D u_{xx}$ : the change of variable converts solutions of the heat equation into solutions of the advection-diffusion equation. The computation is reversible, so it is a bijection between the two solution sets.

Step 2 (apply the heat kernel). Let $w(y,t) = \int S_D(y-z,t)g(z)\,dz$ . By the heat-kernel theorem, applied with diffusivity $D$ , $w$ is $C^\infty$ for $t>0$ , solves $w_t = Dw_{yy}$ , and $w(y,t)\to g(y_0)$ as $(y,t)\to(y_0,0)$ . Step 1 turns this into the stated $u$ , which is therefore $C^\infty$ and solves the equation; and since $(x,t)\to(x_0,0)$ forces $(x-Vt,t)\to(x_0,0)$ , the initial condition carries over.

Step 3 (mass, mean, variance). Each is a property of $S_D$ transported by Step 1. $\int S_D(\cdot,t) = 1$ gives conservation of mass by Tonelli's theorem, the integrand being non-negative. Treating $u(\cdot,t)/M$ as a probability density — legitimate because $M$ is finite and positive by hypothesis — it is the convolution of the two densities $g/M$ and $S_D(\cdot,t)$ , translated by $Vt$ . The mean and variance of $g/M$ exist because its second moment was assumed finite, which as noted in part 1 also makes its first moment finite, and those of $S_D(\cdot,t)$ are $0$ and $2Dt$ by the definition of $S_D$ . Convolution of two densities with finite second moments adds their means and adds their variances, and a translation adds to the mean only. This is part 1. Substituting $t = L/V$ into $\mathrm{var}(t) = \mathrm{var}(0) + 2Dt$ gives part 2, with $2DL/V = L^2\cdot 2D/(VL) = 2L^2/\mathrm{Pe}$ , so that the term diffusion contributes to the standard deviation is $L\sqrt{2/\mathrm{Pe}}$ .

Step 4 (decay term). If $u = e^{-rt}w$ then $u_t + Vu_x - Du_{xx} + ru = e^{-rt}\bigl(w_t - rw + Vw_x - Dw_{xx} + rw\bigr) = e^{-rt}(w_t + Vw_x - Dw_{xx})$ , which vanishes precisely when $w$ solves the undecayed equation. $\square$

NoteAdaptedAdapted from Stakgold & Holst (2011), §8.2, Exercise 8.2.1, p. 486; Strauss (2008), §2.4, Exercises 16 and 18, p. 54

Both reductions used above are the ones the books themselves point to. Strauss (Strauss (2008), Exercises 16 and 18, p. 54) sets the moving frame $y = x - Vt$ for $u_t - ku_{xx} + Vu_x = 0$ and the factor $e^{-bt}$ for $u_t - ku_{xx} + bu = 0$ as exercises with exactly those hints. Stakgold and Holst (Stakgold & Holst (2011), Exercise 8.2.1, p. 486) pose the same problem for a pollutant of concentration $C$ with constant diffusivity $D$ and a time-dependent wind $U(t)$ , asking for the transformation $t' = t$ , $x' = x - \mathcal{V}(t)$ that removes the convection term, and they call the curves $x - \mathcal{V}(t) = \text{constant}$ drift curves. For $\mathcal{V}(t) = \int_0^t V(s)\,ds$ Step 1 goes through unchanged with $Vt$ replaced by $\mathcal{V}(t)$ , since $\partial_t \mathcal{V} = V(t)$ is all that the chain rule uses: the cloud follows the drift curve through its release point while spreading as $\sqrt{2Dt}$ .

Application(Travelling fronts of a bistable reaction-diffusion equation)AdaptedAdapted from Evans (2010), §4.2.1c, Eqs. (12)–(16) and (26), pp. 170–175

Consider the scalar reaction-diffusion equation

u_t - u_{xx} = f(u) \qquad \text{in } \mathbb{R}\times(0,\infty) ,

where $f : \mathbb{R}\to\mathbb{R}$ is smooth and, for some point $0 < a < 1$ , satisfies

(a) $f(0) = f(a) = f(1) = 0$ ;
(b) $f < 0$ on $(0,a)$ and $f > 0$ on $(a,1)$ ;
(c) $f'(0) < 0$ and $f'(1) < 0$ ;
(d) $\int_0^1 f(z)\,dz > 0$ .

A travelling wave is a solution of the form $u(x,t) = v(x - \sigma t)$ . Substituting into the equation gives the ordinary differential equation

v'' + \sigma v' + f(v) = 0 ,

and to it one imposes the conditions that make the wave a front between the two stable states, Evans' (16):

\lim_{s\to+\infty}v(s) = 1,\qquad \lim_{s\to-\infty}v(s) = 0,\qquad \lim_{s\to\pm\infty}v'(s)=0 .

These limits do not follow from the substitution; they are part of what is being asked for. Evans shows by a phase-plane argument that there exists a velocity $\sigma < 0$ for which this problem has a solution, so the equation has a travelling wave interpolating between the two stable states $0$ and $1$ ; and he states that a more refined analysis shows this velocity is unique (Evans (2010), p. 175). Hypothesis (c) is what makes $0$ and $1$ stable as constant solutions, with the intermediate constant $a$ unstable; the front is therefore an invasion of one stable state by the other. This is the existence theorem of travelling waves, restated here because the identity below and the examples need its hypotheses in front of them; that unit works Evans' argument through in detail.

NoteAdaptedAdapted from Evans (2010), §4.2.1c, pp. 171–175

No proof is given here, for two reasons. The phase-plane route — writing $w = v'$ , reading a front as a trajectory joining the saddle at $(0,0)$ to the saddle at $(1,0)$ , and moving $\sigma$ until the unstable and stable curves meet — is set out step by step in travelling waves, and repeating it here would duplicate that unit. And the source itself stops short of a complete argument: Evans says he outlines the analysis "without complete proofs" (Evans (2010), p. 171), and the uniqueness of $\sigma$ is asserted rather than demonstrated (Evans (2010), p. 175). What this unit adds is what can be settled in full: the identity below, which fixes the sign of $\sigma$ from $f$ alone, and the exact front in the second example.

Application(The velocity is fixed by the areas)AdaptedAdapted from Evans (2010), §4.2.1c, the function $E$ and its derivative, p. 173

Let $f$ be continuous and let $v$ solve $v'' + \sigma v' + f(v) = 0$ on $\mathbb{R}$ with $v(-\infty) = 0$ , $v(+\infty) = 1$ and $v'(\pm\infty) = 0$ , with $v' \in L^2(\mathbb{R})$ . Then

\sigma \int_{-\infty}^{\infty} v'(s)^2\,ds \;=\; -\int_0^1 f(z)\,dz .

In particular $\sigma$ has the opposite sign to $\int_0^1 f$ , so Evans' hypothesis (d) forces $\sigma<0$ : the state $1$ , the one with the greater "area" under $f$ , advances. If $\int_0^1 f = 0$ the front is stationary.

ProofAdaptedAdapted from Evans (2010), p. 173

Evans' energy is $E(v,w) = \tfrac12 w^2 + \int_0^{v} f(z)\,dz$ evaluated along a trajectory, with $w = v'$ . Along a solution,

\frac{d}{ds}E\bigl(v(s),v'(s)\bigr) = v'v'' + f(v)\,v' = v'\bigl(v'' + f(v)\bigr) = -\sigma\, v'(s)^2 ,

the last step by the differential equation. The right-hand side is integrable on $\mathbb{R}$ because $v'\in L^2$ , so $s \mapsto E(v(s),v'(s))$ is an antiderivative of an integrable function. It is $E$ itself, not its derivative, that has limits at the two ends: $E$ is continuous, and $v\to 0$ , $v'\to 0$ as $s\to-\infty$ while $v\to1$ , $v'\to0$ as $s\to+\infty$ , so $E(v(s),v'(s))\to E(0,0)=0$ and $\to E(1,0)=\int_0^1 f$ respectively. Hence, by the fundamental theorem of calculus,

\int_{-\infty}^{\infty}\frac{d}{ds}E\,ds = E(1,0) - E(0,0) = \int_0^1 f(z)\,dz ,

the two end values being the limits just identified. Equating this with the integral of $-\sigma (v')^2$ gives $\int_0^1 f = -\sigma\int (v')^2$ . Since $\int (v')^2 > 0$ (otherwise $v$ is constant and cannot have two different limits), the sign statement follows. $\square$

Remark(Where the hypotheses are used)Original

$D>0$ in the moving-frame theorem. Step 2 is the heat-kernel theorem, which needs $D>0$ . With $D = 0$ the solution is $g(x-Vt)$ , pure translation with no spreading, and the problem is the first-order one; the moving frame is then the characteristic (characteristics).
$V$ constant is used only in Step 1, and only through $\partial_t(Vt) = V$ ; a $V$ depending on $t$ alone is as easy, as the note above says. A $V$ depending on $x$ is a different problem: the frame change is no longer a rigid shift, the flux $Vu$ contributes $V_x u$ as well, and no reduction to the heat equation is available in general.
Bistability, that is Evans' (c). If $f'(0)>0$ the state $0$ is unstable and the structure of the travelling-wave problem changes: $(0,0)$ is then not a saddle, and the conclusion of the theorem — one wave with one velocity — need not hold. Hypotheses (a)-(d) are stated exactly as Evans gives them, and the statement above claims nothing outside them.
$\int_0^1 f \ne 0$ is what makes the front move, by the identity. It is a statement about areas, not about shapes: two nonlinearities with the same $\int_0^1 f$ but very different profiles give fronts of the same sign of velocity, though not the same speed, because $\int(v')^2$ differs.
The Peclet number needs a length. $\mathrm{Pe}$ is not a property of the equation alone; different choices of $L$ give different numbers for the same physical set-up, and a statement such as "advection dominates" is meaningless until $L$ is named.

Example(A release into a river)Original

A mass $M$ of tracer is released at $x=0$ , $t=0$ into a channel with mean speed $V = 0.5\,\mathrm{m/s}$ and longitudinal dispersion coefficient $D = 5\,\mathrm{m^2/s}$ .

A point release is not one of the continuous bounded data the statement above admits, so we do not use its formula. What we use instead is the drifting kernel itself: by Step 1 of the proof, $u(x,t) = M\,S_D(x - Vt,\,t)$ , that is

u(x,t) = \frac{M}{\sqrt{4\pi D t}}\,\exp\!\Bigl(-\frac{(x-Vt)^2}{4Dt}\Bigr) ,

solves $u_t + Vu_x = Du_{xx}$ for every $t>0$ , carries mass $M$ at every time, and concentrates at $x = 0$ as $t \to 0^+$ — the point-source reading of the kernel described in the heat kernel. It is a Gaussian of mean $Vt$ and variance $2Dt$ , which is what parts 1 and 2 of the statement would say if $\mathrm{var}(0)=0$ were allowed in them. It is not — those parts need a continuous $g\ge0$ , whose variance is positive — so read the point release as the limiting case they approach, and the formulas below as read off the kernel directly rather than quoted from the statement.

Take $L = 1000\,\mathrm{m}$ , the distance to a water intake downstream.

Peclet number. $\mathrm{Pe} = VL/D = (0.5)(1000)/5 = 100$ . Advection dominates: the tracer arrives as a recognisable pulse, not a slow seepage.
Arrival. The maximum of the travelling profile reaches the intake at $t = L/V = 2000\,\mathrm{s} = 33.3\,\mathrm{min}$ . The concentration measured at the intake peaks a little earlier, at $t = 1980.1\,\mathrm{s}$ : the $t^{-1/2}$ prefactor is still falling while the cloud arrives, and it pulls the maximum about $1\%$ early.
Width. At $t = L/V$ the standard deviation of the cloud is $\sqrt{2Dt} = \sqrt{2\cdot5\cdot2000} = 141.4\,\mathrm{m}$ , which is $L\sqrt{2/\mathrm{Pe}} = 1000\sqrt{0.02} = 141.4\,\mathrm{m}$ : the formula of part 2 in its $\mathrm{var}(0)=0$ limit.
Peak concentration. $u_{\max} = M/\sqrt{4\pi D t} = M/354.5$ , in units of mass per metre: a release of $10\,\mathrm{kg}$ gives $28.2\,\mathrm{g/m}$ at the peak.
Duration. Two standard deviations either side of the centre the Gaussian factor is $e^{-2} = 0.135$ , so the intake sees concentrations above about $14\%$ of the peak for a time $4\sqrt{2Dt}/V = 4(141.4)/0.5 = 1131\,\mathrm{s}$ , about $19\,\mathrm{min}$ . (The letter $\sigma$ is reserved in this unit for the velocity of a travelling front.)
What would change at low Peclet number. Were the channel ten times slower, $\mathrm{Pe} = 10$ and the relative spread $\sqrt{2/\mathrm{Pe}} = 0.45$ : the cloud would be nearly as long as its journey, and tracer would have reached the intake long before the peak did, since the Gaussian is positive everywhere: the leading edge of the cloud runs ahead of its centre, downstream of the peak.

Example(An exact bistable front)Original

Take the cubic $f(u) = u(1-u)(u-a)$ with $0 < a < \tfrac12$ , so that Evans' hypotheses hold: $f$ vanishes at $0,a,1$ ; $f<0$ on $(0,a)$ and $f>0$ on $(a,1)$ ; $f'(0) = -a < 0$ and $f'(1) = a-1 < 0$ ; and

\int_0^1 f(z)\,dz = \int_0^1\bigl(z^2 - z^3 - az + az^2\bigr)dz = \frac{1}{12} - \frac{a}{6} = \frac{1-2a}{12} > 0 .

Find the front. Try the logistic profile $v(s) = \bigl(1+e^{-bs}\bigr)^{-1}$ , which satisfies $v' = b\,v(1-v)$ and therefore $v'' = b\,v'(1-2v) = b^2 v(1-v)(1-2v)$ . Substituting into $v''+\sigma v'+f(v)=0$ and dividing by $v(1-v)$ , which is positive,

b^2(1-2v) + \sigma b + (v - a) = 0 \qquad\Longleftrightarrow\qquad (1 - 2b^2)\,v + \bigl(b^2 + \sigma b - a\bigr) = 0 .

This holds for all $v$ exactly when both brackets vanish:

b = \frac{1}{\sqrt2}, \qquad \sigma = \frac{a - b^2}{b} = \sqrt2\Bigl(a - \tfrac12\Bigr) .

So the front and its velocity are

u(x,t) = \frac{1}{1 + \exp\bigl(-(x-\sigma t)/\sqrt2\bigr)}, \qquad \sigma = \sqrt{2}\,\bigl(a-\tfrac12\bigr) < 0 ,

with $v(-\infty)=0$ , $v(+\infty)=1$ , $v'(\pm\infty)=0$ , exactly as the statement requires — and $\sigma<0$ , as the area identity demands.

Check the identity. With $a = \tfrac14$ : $\sigma = \sqrt2(-\tfrac14) = -0.353553$ and $\int_0^1 f = \tfrac{1}{24} = 0.0416667$ . For the integral of $(v')^2$ , use $v' = \tfrac{1}{\sqrt2}v(1-v)$ and $v(1-v) = \tfrac14\operatorname{sech}^2(s/2\sqrt2)$ , together with $\int\operatorname{sech}^4\theta\,d\theta = \tfrac43$ :

\int_{-\infty}^{\infty}(v')^2 ds = \frac12\cdot\frac{1}{16}\int \operatorname{sech}^4\!\Bigl(\frac{s}{2\sqrt2}\Bigr)ds = \frac{1}{32}\cdot 2\sqrt2 \cdot \frac43 = \frac{\sqrt2}{12} = 0.117851 .

Then $\sigma\int(v')^2 = (-0.353553)(0.117851) = -0.0416667 = -\int_0^1 f$ . The identity is satisfied to every digit, for every $a$ . (The front itself, and the way the guess $v' = b\,v(1-v)$ selects $b$ and $\sigma$ , is worked out at more length in travelling waves; what is new here is the check against the area identity.)

Numbers with dimensions. For $u_t = D u_{xx} + r\,u(1-u)(u-a)$ , scaling $x = \sqrt{D/r}\,\xi$ and $t = \tau/r$ returns the equation above, so the physical speed is $\sqrt{Dr}\,|\sigma|$ and the front thickness is $\sqrt{D/r}$ times the dimensionless one. The profile rises from $0.1$ to $0.9$ over $2\sqrt2\ln 9 = 6.215$ dimensionless units. For a gas-phase front with $D = 0.2\,\mathrm{cm^2/s}$ and $r = 10^4\,\mathrm{s^{-1}}$ and $a = \tfrac14$ , the speed is $\sqrt{2000}\cdot0.353553 = 15.8\,\mathrm{cm/s}$ and the front is $6.215\sqrt{2\times10^{-5}} = 0.028\,\mathrm{cm}$ thick: thin and fast, which is what a flame looks like.

ApplicationOriginal

Which term to keep. The Peclet number decides how to model, how to measure and how to compute. In a river or a ventilation duct $\mathrm{Pe}$ is in the hundreds, so a pollutant arrives as a translating pulse and the measurement that matters is its arrival time $L/V$ ; in a stagnant pore, a gel, or across a cell membrane $\mathrm{Pe}$ is far below one and the velocity can be dropped. Chromatography lives in between: the retention time is set by $V$ and the peak width, which sets the resolution, by $D$ — and by part 2 of the moving-frame statement the peak width relative to the distance travelled is $\sqrt{2/\mathrm{Pe}}$ , so the resolving power of a column improves only as $\sqrt{\mathrm{Pe}}$ , that is as the square root of its length.

Numerical warning. When the advection-diffusion equation is discretised, the relevant Peclet number uses the cell size: $\mathrm{Pe}_{\text{cell}} = V\Delta x/D$ . Centred differences for $Vu_x$ combined with centred differences for $Du_{xx}$ produce a scheme whose coefficients lose their sign once $\mathrm{Pe}_{\text{cell}} > 2$ , and the computed solution then oscillates around steep fronts — the same loss of a discrete maximum principle seen in the maximum principle. The classical responses are to refine until $\mathrm{Pe}_{\text{cell}} \le 2$ or to use one-sided (upwind) differences for the advection term; see finite differences.

The probabilistic reading. Grimmett and Stirzaker (Grimmett & Stirzaker (2020), Example (10), p. 584) derive, for a one-dimensional Brownian particle with constant drift rate $m$ and instantaneous variance $\beta^2$ , the forward equation $\partial_t f = -m\,\partial_y f + \tfrac12\beta^2\partial_{yy} f$ and observe that the corresponding process is $\beta W(t) + mt$ . (They write $\sigma$ for what is called $\beta$ here; in this unit $\sigma$ is the front velocity.) That is precisely the moving-frame theorem, in probabilistic clothing: $V = m$ , $D = \tfrac12\beta^2$ , and "carry the diffusing profile along at speed $V$ " is "add a deterministic drift $mt$ to a Wiener path". The same correspondence makes the advection-diffusion equation the simplest Fokker-Planck equation.

Fronts. The bistable front is the mechanism behind a long list of propagating transitions: the nerve impulse in the Nagumo model, where $u$ is a membrane potential and the two stable states are rest and excitation; flame propagation, where $u$ is a scaled temperature and the Arrhenius factor of (Stakgold & Holst (2011), §0.3, Eq. (0.3.1), p. 12) supplies the nonlinearity; the advance of an advantageous gene through a population; and phase transitions in which one phase invades another. In every case the two numbers the model has to deliver are the ones computed above: the speed, which the area identity ties to $\int_0^1 f$ , and the thickness $\sqrt{D/r}$ , which says how fine a mesh or how fast a camera is needed to see it. When a reaction is coupled to two diffusing species with different diffusivities, uniform states that are stable without diffusion can become unstable with it — the Turing instability, which Evans (Evans (2010), §4.1, pp. 165–166) analyses as a condition on the eigenvalues, and which is the standard account of biological pattern formation.