Method of Characteristics

The method is the backbone of every code that resolves shocks. Modern computational fluid dynamics does not integrate the Euler equations blindly; it solves a Riemann problem — Riemann's two-constant-state jump — at each cell interface thousands of times per time step, and Godunov-type finite-volume schemes, ENO and WENO reconstructions, and approximate Riemann solvers (Roe, HLLC) are all industrial descendants of the characteristic analysis in this article. The same mathematics designs supersonic inlets and predicts sonic-boom footprints for quiet-supersonic aircraft, where compression characteristics coalescing into a shock is exactly the calculation an aerodynamicist runs. In traffic engineering the Lighthill–Whitham–Richards model, $\rho_t+q'(\rho)\rho_x=0$, is solved by characteristics to place ramp meters and predict the upstream march of a jam's tail at the Rankine–Hugoniot speed. In seismology and non-destructive testing, ray tracing through a varying medium is the characteristic (eikonal) solution of the high-frequency wave equation, and fast-marching and fast-sweeping algorithms compute first-arrival times by advancing along characteristics. Even general relativity uses the idea: the characteristic (null) surfaces of Einstein's equations are light cones, and characteristic and null-cone formulations extract gravitational waveforms at future null infinity. Wherever information propagates at finite speed — and especially wherever it piles up into a discontinuity — the characteristic curve is the natural coordinate.

Picture the solution as a landscape $z=u(x,y)$ and the PDE as a steady wind $(a,b,c)$ blowing through three-dimensional space. The integral curves of that wind are the characteristics, and the solution surface is simply ruled out of them, like a fabric woven from threads all pinned to one initial seam — the Cauchy-data curve $\Gamma$. Because each thread obeys an ordinary differential equation, building the surface costs only integration, not the full apparatus of PDE theory. The one place the weave can fail is where two threads would have to pass through the same base point carrying different heights: there the smooth surface tears and a shock is born. The first figure fixes this geometry — data transported along characteristics, unique wherever $\Gamma$ cuts the flow transversally.

Now let the characteristic speed depend on $u$ itself. Each point of the initial profile moves with its own velocity, so a compressive region — where taller, faster data sit behind shorter, slower data — folds the profile forward until it goes vertical. The instant of vertical tangent is the moment neighbouring characteristics first intersect, and past it the naive formula returns three values of $u$ at one point. The physics forbids that for a conserved quantity, and the resolution is a shock: a moving discontinuity that swallows the overturned region. The second figure computes this for a concrete ramp, where every characteristic of the compressive fan focuses on a single point and the breaking time comes out to an exact number.

For second-order equations the same word — characteristic — names the curves along which the highest derivatives are not fixed by the data. Two real families mean signals propagate (hyperbolic, the wave equation and its light cones); one repeated family means diffusion (parabolic, heat); no real family means no propagation and instantaneous global influence (elliptic, Laplace). The transport picture of this article applies exactly to the first case and degenerates gracefully into the classification of the second.

The natural home of the method is the quasilinear first-order equation. Its left-hand side is a directional derivative, which is what makes the whole construction geometric rather than merely formal.

where: $u(x,y)$ — the unknown scalar field; $u_x,\,u_y$ — its first partial derivatives $\partial u/\partial x,\ \partial u/\partial y$; $a,b,c$ — given coefficient functions, allowed to depend on $u$ (which is what "quasilinear" means). Reading it: the combination $a\,u_x+b\,u_y$ is the rate of change of $u$ in the direction $(a,b)$, so the equation prescribes that rate to be $c$ and constrains nothing transverse. Demanding that $u$ satisfy (1) along a curve tangent to $(a,b)$ turns the PDE into three coupled ODEs — the characteristic system that the whole method rests on.

where: $s$ — a parameter running along one characteristic; $(x(s),y(s))$ — the projection of the characteristic onto the base plane; $u(s)$ — the value of the solution carried along it. Reading it: the first two equations trace the curve, the third transports the solution, and Derivation 1 shows they follow from the chain rule alone. Eliminating the parameter $s$ gives Lagrange's symmetric form, the shape in which the characteristic equations are usually integrated by inspection.

where: each ratio equals $\mathrm{d}s$; setting any two equal removes $s$ and yields a first integral, and two independent first integrals give the general solution. Reading it: solving a first-order PDE is reduced to finding two conserved combinations of $x,y,u$ along the flow. Whether the construction succeeds at all is a single determinant — the transversality, or non-characteristic, condition on the data curve $\Gamma$.

where: $(\dot{x}_0,\dot{y}_0)$ — the tangent to $\Gamma$; $(a,b)$ — the characteristic direction evaluated on $\Gamma$; $J$ — their Jacobian. Reading it: $J\neq0$ says $\Gamma$ is nowhere parallel to a characteristic, so exactly one characteristic issues from each data point and the local solution is unique. Where $J=0$, the problem is over- or under-determined. For the archetypal nonlinear case, the inviscid Burgers equation $u_t+u\,u_x=0$, integrating the characteristic system gives the solution in implicit form.

where: $f$ — the initial profile $u(x,0)=f(x)$; $x-ut$ — the foot of the characteristic on the initial line, a quantity conserved along it. Reading it: $u$ is constant on straight characteristics of $u$-dependent slope, so the profile advects itself; because the slope varies, the implicit relation can fail to be single-valued. The failure has an exact time. Differentiating (5) shows the gradient blows up when neighbouring characteristics first meet.

where: $f'(x_0)$ — the slope of the initial profile; $t_{\text{shock}}$ — the earliest time at which $u_x\to\infty$. Reading it: a shock forms only where $f'<0$ (a compression), and the first break happens at the steepest such slope. If $f'\ge0$ everywhere the solution stays smooth for all time. Past the break, the discontinuity's speed is fixed not by the differential equation but by the underlying conservation law — the Rankine–Hugoniot condition.

where: $f(u)$ — the flux (for Burgers $f=\tfrac12u^2$, giving $s=\tfrac12(u_L+u_R)$); $u_L,u_R$ — the states on the left and right of the jump; $[\cdot]$ — the jump across it. Reading it: the shock moves at the average of the two states' velocities, and the entropy condition $u_L>u_R$ picks the physical shock over an unphysical rarefaction-shock. Fully nonlinear equations $F(x,y,u,p,q)=0$ need the gradient $(p,q)=(u_x,u_y)$ tracked as well, giving the Charpit–Lagrange system on the five-dimensional space of one-jets.

where: $F(x,y,u,p,q)=0$ — the general first-order equation with $p=u_x,\ q=u_y$; $F_p,F_q,\dots$ — partial derivatives of $F$. Reading it: the characteristics now live in $(x,y,u,p,q)$ space; projected to the base plane they are the light rays of geometrical optics when $F$ is the eikonal equation. Finally, the classification of second-order operators falls out of the same idea: the characteristic curves solve a quadratic.

where: $a,b,c$ — the principal-part coefficients of $a\,u_{xx}+b\,u_{xy}+c\,u_{yy}=\cdots$; the discriminant $b^2-4ac$ decides everything — $>0$ hyperbolic (two real characteristics), $=0$ parabolic (one), $<0$ elliptic (none). Reading it: real characteristics are exactly the curves along which second derivatives are indeterminate, so "how many real roots" is "how many propagation directions".

The single method specializes cleanly across the whole zoo of first- and second-order equations; the table reads left to right as "which equation, what the characteristics do, what it means physically".

1 · Breaking time and shock speed for a compressive ramp (fully worked)

Solve $u_t+u\,u_x=0$ with $u(x,0)=f(x)$ equal to $1$ for $x\le 0$, to $1-x$ for $0\le x\le 1$, and to $0$ for $x\ge 1$. Find the breaking time, the point where the shock first appears, and the speed at which it then travels.

Step 1. Characteristics carry constant $u$: from $x_0$ the line is $x=x_0+f(x_0)\,t$. On the ramp ($0\le x_0\le1$) this is $x=x_0+(1-x_0)\,t$. Set $t=1$: $x=x_0+(1-x_0)=1$ for every $x_0$ in the ramp — the entire compressive fan focuses on the single point $x=1$. So the classical solution cannot survive $t=1$.

Step 2. Confirm with the formula. On the ramp $f'(x_0)=-1$, and this is the minimum slope, so $\min f'=-1$ and $t_{\text{shock}}=-1/\min f'=-1/(-1)=\mathbf{1}$. The break happens at $t=1$, and tracking the focusing lines puts it at $\mathbf{x=1}$. Everywhere $f'\ge0$ there is no crossing, consistent with the shock coming only from the ramp.

Step 3. Determine the states across the nascent shock. Just left of $x=1$ the data came from the top of the ramp, $u_L=1$; just right, from the flat region, $u_R=0$. The entropy condition $u_L>u_R$ ($1>0$) is satisfied, so this is an admissible shock, not a spurious one.

Step 4. Apply Rankine–Hugoniot with flux $f=\tfrac12u^2$: $s=\dfrac{f(u_L)-f(u_R)}{u_L-u_R}=\dfrac{\tfrac12(1)^2-\tfrac12(0)^2}{1-0}=\tfrac12(u_L+u_R)=\mathbf{0.5}$. The shock is born at $(x,t)=(1,1)$ and thereafter follows $x_s=1+0.5\,(t-1)$. Breaking at t = 1, x = 1; the shock then advances at speed s = 0.5.

2 · The full characteristic recipe on a variable-coefficient equation

Solve $x\,u_x+y\,u_y=u$ subject to $u=g(x)$ on the line $y=1$, then evaluate the special case $g(x)=x^2$.

Step 1. Read off $(a,b,c)=(x,y,u)$ and write the characteristic system $\dfrac{\mathrm{d}x}{\mathrm{d}s}=x,\ \dfrac{\mathrm{d}y}{\mathrm{d}s}=y,\ \dfrac{\mathrm{d}u}{\mathrm{d}s}=u$. Each is a simple exponential ODE.

Step 2. Integrate with data at $s=0$ on $y=1$, parametrized by $x_0=\xi$, so $u_0=g(\xi)$: $x=\xi\,e^{s},\ y=e^{s},\ u=g(\xi)\,e^{s}$. The characteristics are rays through the origin.

Step 3. Eliminate the parameters. From $y=e^{s}$ we get $e^{s}=y$; then $x=\xi y$ gives $\xi=x/y$. Substituting, $u=g(\xi)\,e^{s}=y\,g(x/y)$, so the general solution is $\boxed{u(x,y)=y\,g(x/y)}$.

Step 4. Check transversality: on $y=1$ the tangent is $(\dot x_0,\dot y_0)=(1,0)$ and $(a,b)=(x,y)=(x,1)$, so $J=\det\begin{pmatrix}x&1\\1&0\end{pmatrix}=-1\neq0$ — the data curve is non-characteristic and the solution is unique near it.

Step 5. Specialize to $g(x)=x^2$: $u=y\,(x/y)^2=\mathbf{x^2/y}$. Verify directly: $u_x=2x/y$, $u_y=-x^2/y^2$, so $x\,u_x+y\,u_y=2x^2/y-x^2/y=x^2/y=u$, and $u(x,1)=x^2$. The recipe reproduces the exact closed form u = x²/y.

The method of characteristics is a nineteenth-century idea perfected by many hands, but its dramatic climax — the discovery that smooth data can generate a discontinuity — belongs to Bernhard Riemann, who found it in the physics of sound.

Texts. Stone & Goldbart, Mathematics for Physics, Ch. 8 (Partial Differential Equations) · Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., Ch. 9 · Hassani, Mathematical Physics, Ch. 22 (first-order PDEs and characteristics) · Riley, Hobson & Bence, Mathematical Methods for Physics and Engineering, 3rd ed., Ch. 20 · Courant & Hilbert, Methods of Mathematical Physics, Vol. II, Ch. II · Whitham, Linear and Nonlinear Waves, Ch. 2 · Evans, Partial Differential Equations, 2nd ed., §3.2. Key checks across the article: the chain-rule reduction yields $\mathrm{d}u/\mathrm{d}s=c$; re-substituting $u=y\,g(x/y)$ (and its special case $u=x^2/y$) into $x u_x+y u_y=u$ confirms the solution; the ramp gives $t_{\text{shock}}=-1/\min f'=1$ at focusing point $x=1$ with shock speed $s=\tfrac12(u_L+u_R)=0.5$; the coordinate change $\xi=x-ct,\ \eta=x+ct$ reduces the wave operator to $u_{tt}-c^2u_{xx}=-4c^2u_{\xi\eta}$; and the discriminant $b^2-4ac$ governs the hyperbolic/parabolic/elliptic classification of equation (9).

Builds on. Classification of Second-Order Linear PDEs · Ordinary differential equations and existence-uniqueness · Boundary Conditions and Well-Posedness. Leads to. d'Alembert's Solution and Characteristics of the Wave Equation · The Wave Equation · Fundamental Solutions and Green's Functions for PDEs.

Ask the tutor: Want to see a specific initial profile pushed to its breaking time, or the eikonal equation solved as a fully nonlinear characteristic system? Ask and we will trace the curves together.