First-Order Equations and Characteristics

Book 6A was about second-order equations, where diffusion smooths everything. First-order equations are the opposite extreme. A first-order equation carries information along curves, the characteristics, without spreading it out: whatever the initial data look like, the solution is made of the same values, moved. That makes first-order equations solvable by ordinary differential equations, which is the content of the method of characteristics, and it also makes them fragile: when characteristics cross, the solution is asked to take two values at one point, and a smooth solution stops existing.

Both facts matter far beyond first-order equations. Characteristics are how waves carry information (6B.3 The Wave Equation), how conservation laws form shocks (6B.2 Conservation Laws and Shocks), how light rays and seismic rays travel, and, through the Hamilton–Jacobi equation, how classical mechanics and optimal control are written as PDE (6C.9 Hamilton–Jacobi Equations and Viscosity Solutions, 6D.8 Stochastic Control and the HJB Equation). In geometry, the distance function from a point satisfies a first-order equation whose characteristics are geodesics; where geodesics from a point stop minimising, its characteristics cross, and that set, the cut locus, is exactly where the distance function is not smooth (9B.1 Laplacian Comparison). Perelman's reduced distance satisfies a Hamilton–Jacobi equation of the same kind (12A.5 Reduced Distance and Reduced Volume).

By the end of this chapter you will be able to:

• solve the transport equation, with and without sources, and with variable coefficients, by following characteristics;
• solve linear and quasilinear first-order equations in two variables, and recognise when Cauchy data are characteristic;
• derive the characteristic ODEs for a fully nonlinear equation $F(Du, u, x) = 0$, and recognise Hamilton's equations as the special case of a Hamilton–Jacobi equation;
• prove the local existence theorem for noncharacteristic data, by the inverse function theorem;
• compute when characteristics cross for Burgers' equation and for the eikonal equation, and explain what the crossing means.

1.1Travel times and rays

1.2The transport equation

The simplest first-order equation is

with a constant vector $b \in \mathbb{R}^n$. Its left side is the derivative of $u$ in the space-time direction $(b, 1)$. So fix a point $(x, t)$ and look along the line through it in that direction: $z(s) = u(x + sb, t + s)$. By the chain rule,

So $u$ is constant along each such line. The line through $(x, t)$ meets the plane $t = 0$ at $s = -t$, at the point $x - tb$, where $u = g(x - tb)$.

The solution is the initial profile moving with velocity $b$. If $g$ is not differentiable, $u = g(x - tb)$ is still the obvious answer, though it does not satisfy the equation pointwise; it is a weak solution in the sense of distributions (Exercise 1.12). Unlike the heat equation, transport does not smooth: a corner in $g$ stays a corner forever, carried along the characteristic through it.

Sources. For $u_t + b\cdot\nabla u = f(x, t)$, the same computation gives $\dot z(s) = f(x + sb, t + s)$, and integrating from $s = -t$ to $s = 0$,

This is Duhamel's principle (6A.3 The Heat Equation on ℝⁿ) for transport: each instant's source is carried along from where it was injected.

Variable coefficients. For $u_t + b(x, t)\cdot\nabla u = 0$, with $b$ a $C^1$ vector field, the same idea works with lines replaced by the trajectories of the ODE

Along such a trajectory, $\frac{d}{ds}u(X(s), s) = \nabla u\cdot b + u_t = 0$: $u$ is constant. So $u(x, t) = g(X_{x,t}(0))$, where $X_{x,t}$ is the trajectory with $X_{x,t}(t) = x$, followed backward to time $0$. Existence and uniqueness of trajectories, and smooth dependence on the starting point, come from the ODE theory of 2B.10 Ordinary Differential Equations. The function $u$ is a passive scalar: a dye carried by the flow $b$, its value at each particle unchanged.

Transport versus conservation. A related equation describes a density carried by a flow, such as a concentration of particles:

This is the continuity equation, the integral law "the amount in a region changes only by flux through its boundary", made differential exactly as in 6A.1 What a PDE Is. Expanding, $\rho_t + b\cdot\nabla\rho = -\rho\operatorname{div}b$: along trajectories $\rho$ is not constant but changes at the rate $-\rho\operatorname{div}b$, increasing where the flow converges. When $\operatorname{div}b = 0$ (an incompressible flow, such as the rotation above), the two equations coincide. The continuity equation is the first-order part of the Fokker–Planck equation of 6D.3 Stochastic Differential Equations and Feynman–Kac, and the conservation laws of 6B.2 Conservation Laws and Shocks are its nonlinear relatives.

1.3Linear and quasilinear equations in two variables

Consider a first-order equation for $u(x, y)$ of the form

quasilinear because it is linear in the derivatives, with coefficients that may depend on $u$ (6A.1 What a PDE Is). Geometrically: the graph $z = u(x, y)$ is a surface in $(x, y, z)$-space, with upward normal proportional to $(u_x, u_y, -1)$. The equation says

the vector field $V = (a, b, c)$ is tangent to the graph. So the graph is a union of integral curves of $V$, the characteristic curves

Their projections to the $(x, y)$-plane are the characteristic projections, often also just called characteristics.

To pick out one solution we prescribe $u$ along a curve $\Gamma$ in the plane, say $(x, y) = (\alpha(r), \beta(r))$ with $u = \varphi(r)$ there: the Cauchy problem. The lifted curve $(\alpha(r), \beta(r), \varphi(r))$ lies in $(x, y, z)$-space; from each of its points start a characteristic curve, and the surface they sweep out is the graph of the solution. This works as long as the characteristic projections cross $\Gamma$ transversally:

When the data are noncharacteristic, the map $(r, s) \mapsto (x(r, s), y(r, s))$, sending a point $r$ of $\Gamma$ and a time $s$ along the characteristic from it to the point reached, has invertible Jacobian at $(r_0, 0)$: its columns are the tangent $(\alpha', \beta')$ to $\Gamma$ and the characteristic direction $(a, b)$, which are independent. By the inverse function theorem (2B.9 The Inverse and Implicit Function Theorems), it is a diffeomorphism near $(r_0, 0)$, so each point $(x, y)$ near $\Gamma$ lies on exactly one characteristic, at a definite time $s$; define $u(x, y) = z(r, s)$ there. This is the local existence theorem in two variables, and the general case below is the same argument.

If the data are characteristic, $\Gamma$ itself is (the projection of) a characteristic curve, and the data on it cannot be arbitrary: they must satisfy the ODE $\dot z = c$ along it. If they do, there are infinitely many solutions; if not, none. Either way the problem is not well posed.

1.3.1Burgers' equation and the first crossing

The most important quasilinear example is the inviscid Burgers equation

Here $(t, x)$ play the roles of $(x, y)$, and the characteristic equations are $\dot t = 1$, $\dot x = z$, $\dot z = 0$. So $z$ is constant along each characteristic, and the projected characteristics are straight lines whose speed is the value carried:

The solution is given implicitly by $u(x, t) = g(\xi)$ where $x = \xi + tg(\xi)$, that is,

Points where $u$ is large move fast. If $g$ is increasing, faster points start to the right of slower ones, characteristics fan out, and all is well. If $g$ decreases anywhere, faster points start behind slower ones, and eventually overtake them.

The quantity $u_x$ satisfies its own ODE along characteristics, $\frac{d}{dt}u_x = -u_x^2$ (Exercise 1.13), the Riccati equation $y' = -y^2$ whose solutions blow up in finite time from negative data. This is the same mechanism as the blow-up of $y' = y^2$ in 6A.8 Nonlinear Parabolic Equations, now forced by geometry rather than by a reaction term. What happens after $T^*$ is the subject of 6B.2 Conservation Laws and Shocks: the solution develops a jump, a shock, which moves at a speed fixed by conservation of $\int u$.

1.4Fully nonlinear equations

Now let $U \subseteq \mathbb{R}^n$ be open and consider the general first-order equation

where $F = F(p, z, x)$ is a smooth function of $p \in \mathbb{R}^n$, $z \in \mathbb{R}$, $x \in U$. (Here $Du = \nabla u$, and we follow Evans's letters: $p$ for the gradient slot, $z$ for the value slot.) The eikonal equation is $F = |p|^2 - c(x)^{-2}$; a Hamilton–Jacobi equation $u_t + H(Du, x) = 0$ fits with $(x, t)$ as the variables.

The idea is the same as before: find curves $x(s)$ along which the values $z(s) = u(x(s))$ and the gradients $p(s) = Du(x(s))$ can be computed by solving ODEs. The difficulty is that the gradient is now part of the unknown.

Derivation. Suppose $u$ is a $C^2$ solution and $x(s)$ is a curve, and set $z(s) = u(x(s))$, $p(s) = Du(x(s))$. Then

This involves second derivatives, which we do not know. But differentiate the equation $F(Du, u, x) = 0$ with respect to $x_i$:

Since $u_{x_jx_i} = u_{x_ix_j}$, the second derivatives appear in the same combination in both formulas if we choose the curve so that $\dot x^j = F_{p_j}$. Then

with no second derivatives left. And $\dot z = Du\cdot\dot x = p\cdot D_pF$. So:

We have shown: if $u$ is a $C^2$ solution and $x(s)$ solves $\dot x = D_pF(Du(x), u(x), x)$, then $(x(s), z(s), p(s))$ solves the characteristic equations. The point of the method is the converse: solve the ODEs from the boundary data, and assemble $u$ from the solutions. Note that $F$ itself is constant along solutions of the characteristic equations (Exercise 1.15), so a curve that starts on $\{F = 0\}$ stays there.

Special cases.

• Linear equations $b(x)\cdot Du + c(x)u = 0$: $F = b(x)\cdot p + c(x)z$, $D_pF = b$, so $\dot x = b(x)$ and $\dot z = b\cdot p = -c(x)z$ (using $F = 0$). The equations for $x$ and $z$ decouple from $p$: the method for the transport equation.
• Quasilinear equations $b(x, u)\cdot Du + c(x, u) = 0$: $\dot x = b(x, z)$, $\dot z = -c(x, z)$. Again $p$ is not needed, as in the two-variable section.
• Hamilton–Jacobi equations $u_t + H(Du, x) = 0$. Take the variables $(x, t)$, the gradient slots $(p, q)$ with $q = u_t$, and $F = q + H(p, x)$. Then $\dot t = F_q = 1$, so $t$ itself is the parameter, and

using $q = -H$ on solutions. The first two are Hamilton's equations of classical mechanics, and $\dot z$ is the Lagrangian along the trajectory (by the Legendre transform, 6C.9 Hamilton–Jacobi Equations and Viscosity Solutions). So $u$ is the action, and the characteristics of the Hamilton–Jacobi equation are the trajectories of the mechanical system. This is the link between PDE and mechanics that Hamilton and Jacobi found in the 1830s, and it reappears as the link between optimal control and the HJB equation (6C.9 Hamilton–Jacobi Equations and Viscosity Solutions, 6D.8 Stochastic Control and the HJB Equation).

• The eikonal equation $|Du| = 1$. Take $F = \frac12(|p|^2 - 1)$, so $D_pF = p$, $F_z = 0$, $D_xF = 0$:

The characteristics are straight lines, traversed at unit speed in the direction of the gradient, and $u$ increases at unit rate along them: $u$ grows like distance along straight rays. With variable speed $c(x)$ the rays bend, and they are the geodesics of the metric $c^{-2}|dx|^2$, which is Fermat's principle.

1.5The local existence theorem

Now the converse in general. To keep the notation simple, let the boundary be flat near the point of interest: $U \subseteq \{x_n > 0\}$, with data $u = g$ on a piece $\Gamma$ of the plane $\{x_n = 0\}$ (any $C^1$ boundary can be flattened by a change of coordinates, 6A.1 What a PDE Is). Write points of $\Gamma$ as $y = (y_1, \ldots, y_{n-1}, 0)$.

Initial values for the ODEs. At a point $y \in \Gamma$, the characteristic must start with $x(0) = y$ and $z(0) = g(y)$. What about $p(0)$? The tangential derivatives of $u$ along $\Gamma$ are those of $g$, so $p^i(0) = g_{x_i}(y)$ for $i < n$, and the last component is constrained only by the equation: $F(p(0), g(y), y) = 0$. A triple $(p_0, z_0, x_0)$ with $x_0 \in \Gamma$, $z_0 = g(x_0)$, $p_0^i = g_{x_i}(x_0)$ for $i < n$, and $F(p_0, z_0, x_0) = 0$ is called admissible. (For the eikonal equation in Example 1.7 there were two admissible $p_0$ at each point, pointing in and out; we chose the inward one.)

The proof needed exactly one non-degeneracy condition, $F_{p_n} \neq 0$: the characteristics must leave the boundary. Everything else is ODE theory. And it gives only a local solution: near the boundary, for $|s|$ small. How far a solution extends is decided by whether characteristics cross, and that is a global question that the theorem does not address.

1.6Where classical solutions stop

Two different things can stop the method.

Characteristics cross. The map $(y, s) \mapsto x(y, s)$ of Step 2 is a local diffeomorphism near the boundary, but globally it may fail to be injective: two characteristics, carrying different values or different gradients, reach the same point. Burgers' equation and the distance function show both cases:

• In Burgers' equation the crossing characteristics carry different values of $u$. A single-valued solution has to choose, the natural choice is a jump, and the rule for where the jump goes comes from the integral form of the conservation law (6B.2 Conservation Laws and Shocks).
• For the distance function the crossing characteristics carry the same value but different gradients. The natural continuation is $u = \min$ over the competing characteristics, the first arrival, which is continuous but has a kink. The theory of viscosity solutions (6C.9 Hamilton–Jacobi Equations and Viscosity Solutions) makes "first arrival" a definition, and proves it is the unique such continuation.

The Jacobian degenerates without crossing. The map $(y, s) \mapsto x$ can have a singular Jacobian at a point without two characteristics meeting there exactly, where a family of rays is tangent to a curve, its envelope or caustic. In optics, caustics are the bright curves of light on the bottom of a swimming pool or inside a coffee cup, where rays concentrate. Taylor I, §6.7, analyses how they form. At a caustic the travel-time function stops being smooth, and the wave field is not described by rays alone.

1.7History

First-order PDE were studied by Euler, Lagrange and Monge in the eighteenth century; Lagrange and Paul Charpit (around 1784) found the characteristic method for general first-order equations in two variables, and the general theory of characteristics was completed by Augustin Cauchy in the 1820s. William Rowan Hamilton developed his formulation of optics and mechanics in 1827–1835, and Carl Jacobi turned it into a method for integrating the equations of motion through the Hamilton–Jacobi equation in the 1830s and 1840s. Burgers' equation is named after Johannes Burgers, who used it in the 1940s as a model of turbulence; the formation of shocks in compressible gas flow had been analysed by Bernhard Riemann in 1860. Sethian's fast marching method appeared in 1996.

1.8Exercises

Solution

Along $x = \xi + 2t$, $\frac{du}{dt} = e^{-t}$, so $u = \sin\xi + 1 - e^{-t}$. Hence $u(x, t) = \sin(x - 2t) + 1 - e^{-t}$. Check: $u_t = -2\cos(x - 2t) + e^{-t}$ and $2u_x = 2\cos(x - 2t)$, which add to $e^{-t}$.

Solution

In the variables $(\xi, \tau)$, the function $\phi(\xi + \tau, \tau) = \psi(\xi, \tau)$ has $\psi_\tau = \phi_t + \phi_x$, the Jacobian of the change of variables is $1$, and $u = |\xi|$. So the integral is $\int_{\mathbb{R}}|\xi|\int_0^\infty\psi_\tau\,d\tau\,d\xi = 0$, because $\psi$ has compact support in $\tau > 0$.

Solution

(a) $u_{tx} + u_xu_x + uu_{xx} = 0$. (b) Along $\dot x = u$, $\frac{dw}{dt} = w_t + uw_x = -w^2$. With $w(0) = g'(\xi)$, $\frac{1}{w} = \frac{1}{g'(\xi)} + t$, so $w = \frac{g'(\xi)}{1 + tg'(\xi)}$.

Hint

Maximise $-g'(x) = \frac{2x}{(1 + x^2)^2}$.

Solution

$\frac{d}{dx}\frac{2x}{(1 + x^2)^2} = \frac{2(1 + x^2) - 8x^2}{(1 + x^2)^3} = 0$ at $x^2 = \frac13$, and the maximum (at $x = \frac{1}{\sqrt3}$) is $\frac{2/\sqrt3}{(4/3)^2} = \frac{9}{8\sqrt3} \approx 0.6495$. So $T^* = \frac{8\sqrt3}{9} \approx 1.540$. The characteristic from $\xi = \frac{1}{\sqrt3}$ carries $u = g(\xi) = \frac34$, so the crossing begins at $x = \xi + T^*g(\xi) = \frac{1}{\sqrt3} + \frac{8\sqrt3}{9}\cdot\frac34 = \frac{1}{\sqrt3} + \frac{2}{\sqrt3} = \sqrt3 \approx 1.732$. The front of the bump steepens: the hump moves right, faster at its top.

Solution

The computation is Step 4 of Theorem 1.10. For $F = q + H(p, x)$ with $\dot q = -F_t - F_zq = 0$ (no explicit $t$ or $z$), $q$ is constant, hence so is $H = -q$; directly, $\frac{d}{dt}H(p, x) = D_pH\cdot\dot p + D_xH\cdot\dot x = -D_pH\cdot D_xH + D_xH\cdot D_pH = 0$.

Solution

(a) $H = \frac12|p|^2$, so $\dot x = p$, $\dot p = 0$, $\dot z = |p|^2 - \frac12|p|^2 = \frac12|p|^2$. (b) From $y$, $p = ay$ and $x = (1 + at)y$, $z = \frac12a|y|^2 + \frac12a^2|y|^2t = \frac12a|y|^2(1 + at)$. With $y = \frac{x}{1 + at}$, $u = \frac{a|x|^2}{2(1 + at)}$. Check: $u_t = -\frac{a^2|x|^2}{2(1 + at)^2}$ and $\frac12|Du|^2 = \frac12\frac{a^2|x|^2}{(1 + at)^2}$. (c) If $a < 0$, $x = (1 + at)y = 0$ for every $y$ at $t = -\frac1a$: all characteristics pass through the origin then, and $u$ blows up there. This is the PDE picture of a converging lens.

Solution

The distance from $(x, y)$ to the side $x = 1$ is $1 - x$, and similarly for the other sides; the distance to $\partial U$ is the minimum over the four sides, which is $\min(1 - |x|, 1 - |y|)$. Off the diagonals one term is strictly smaller and $u$ is locally $1 - |x|$ or $1 - |y|$ (or the same with signs), with gradient of length $1$. On a diagonal the two terms agree but have different gradients, $(\mp1, 0)$ and $(0, \mp1)$, so $u$ has a ridge. The horizontal characteristics from the vertical sides meet the vertical ones from the horizontal sides there: the diagonals are the medial axis of the square.

Solution

(a) $\sin\theta_c = \frac{c_1}{c_2} = \frac23$, $\theta_c \approx 41.8°$. (b) $T_y^2 = c_2^{-2} - \frac{\sin^2\theta_1}{c_1^2}$, negative when $\sin\theta_1 > \frac{c_1}{c_2}$. No admissible gradient exists, so no transmitted ray: the light is reflected.

Solution

(a) $\ell_\tau = -\frac{|x|^2}{4\tau^2}$ and $\nabla\ell = \frac{x}{2\tau}$, so $|\nabla\ell|^2 = \frac{|x|^2}{4\tau^2}$. (b) $\dot x = D_pH = 2p$, $\dot p = 0$: straight lines traversed at constant velocity $2p$. Along $x = \tau v$, $\nabla\ell = \frac v2 = p$, and $2p = v = \dot x$ as required; $\dot z = 2|p|^2 - |p|^2 = \frac{|v|^2}{4}$, and indeed $\frac{d}{d\tau}\frac{|\tau v|^2}{4\tau} = \frac{|v|^2}{4}$. (c) $\Delta\ell = \frac{n}{2\tau}$, so $2\Delta\ell - |\nabla\ell|^2 + \frac{\ell - n}{\tau} = \frac n\tau - \frac\ell\tau + \frac\ell\tau - \frac n\tau = 0$.