The Method of Characteristics for First-Order PDEs
This post aims to introduce the method of characteristics for first-order partial differential equations as clearly and concisely as possible. Personally, I think this is one of the relatively simple yet surprisingly confusion-prone parts of PDE theory, so I'm going to try to explain it in my own words. Of course, more accurately, this is really just my own personal memo on the subject.
The Quasilinear Case
General procedure
Consider the partial differential equation
$$\begin{equation}\boldsymbol{\alpha}(\boldsymbol{x},u) \cdot \frac{\partial}{\partial \boldsymbol{x}} u = \beta(\boldsymbol{x},u)\end{equation}$$
where $\boldsymbol{\alpha}$ is an $n$-dimensional vector function, $\beta$ is a scalar function, $\cdot$ denotes the dot product of vectors, $u\equiv u(\boldsymbol{x})$ is a function of $n$ variables, and $\boldsymbol{x}$ are its independent variables. more
The idea behind the method of characteristics is to treat $\boldsymbol{x}$ as a function of some parameter $s$, so that $\boldsymbol{x}(s)$ effectively becomes the parametric equation of some high-dimensional curve — the so-called characteristic curve. In this way, $u$ also becomes a function of the parameter $s$. We then have
$$\begin{equation}\frac{du}{ds}=\frac{\partial u}{\partial \boldsymbol{x}}\cdot\frac{d\boldsymbol{x}}{ds}\end{equation}$$
Comparing this with the original PDE $(1)$, we notice that we can set
$$\begin{equation}\frac{d\boldsymbol{x}}{ds}=\boldsymbol{\alpha}(\boldsymbol{x},u)\end{equation}$$
which gives us
$$\begin{equation}\frac{du}{ds}=\beta(\boldsymbol{x},u)\end{equation}$$
Combining $(3)$ and $(4)$, we obtain a system of ordinary differential equations consisting of $n+1$ equations.
$$\begin{equation}\left\{\begin{aligned}&\frac{d\boldsymbol{x}}{ds}=\boldsymbol{\alpha}(\boldsymbol{x},u)\\ &\frac{du}{ds}=\beta(\boldsymbol{x},u)\end{aligned}\right.\end{equation}$$
Since $s$ is just an auxiliary variable that we introduced, in principle we can solve for results that don't depend on $s$:
$$\begin{equation}\boldsymbol{c}=\boldsymbol{f}(\boldsymbol{x},u)\end{equation}$$
where $\boldsymbol{c}$ is an $n$-dimensional vector representing the integration constants of this ODE system, and $\boldsymbol{f}$ is an $n$-dimensional vector function. What remains is to determine the relations among the various integration constants based on the initial conditions. Of course, if we want to obtain a general-solution expression, it would be
$$\begin{equation}G(\boldsymbol{f}(\boldsymbol{x},u))=0\end{equation}$$
where $G$ is an arbitrary function of $n$ variables, from which we can in principle solve for $u$ as a function of $\boldsymbol{x}$.
A simple example
The procedure above is fairly abstract; in actual practice we can be more flexible. Let's now solve
$$\begin{equation}\frac{\partial u}{\partial t} + x \frac{\partial u}{\partial x} = u^2,\quad u(x,0)=f(x)\end{equation}$$
We obtain the characteristic equations
$$\begin{equation}dt = \frac{dx}{x}=\frac{du}{u^2}\end{equation}$$
Solving these gives
$$\begin{equation}x=C_1 e^t, u = \frac{1}{C_2-t}\end{equation}$$
When $t=0$, we have $x=C_1,u=\frac{1}{C_2}=f(C_1)$, from which we solve $C_2 = \frac{1}{f(C_1)}$. Since we also have $C_1 = xe^{-t}, C_2 = u^{-1} + t$, substituting gives
$$\begin{equation}u^{-1} + t = \frac{1}{f(xe^{-t})}\end{equation}$$
that is
$$\begin{equation}u = \frac{f(xe^{-t})}{1-t\times f(xe^{-t})} \end{equation}$$
A bit of discussion
What does a characteristic curve actually mean? For beginners, the process above might look like a magic trick — first solving for constants, then eliminating them — without a clear sense of what's really going on. This was exactly my own confusion when I first learned the method of characteristics.
In fact, one useful way to think about it is this: a characteristic curve is already a solution of the PDE, except that it's just one line lying on the solution, whereas the full solution should really be a higher-dimensional surface. Clearly, a moving point traces out a line, and a moving line sweeps out a surface — so if we can find a way to let these lines "move," we obtain the equation of that surface. In other words, the integration constants $\boldsymbol{c}$ need to be allowed to "move."
Of course, they can't move around with total freedom — if they did, they might sweep out the entire space. How they're allowed to move is determined by the initial conditions, so we use the initial conditions to pin down the constraints among the integration constants. Once that's done, we've effectively obtained the parametric equations of this surface. From a purely "solving" standpoint, there's no strict need to eliminate the constants, but since we often prefer explicit solutions, we go ahead and eliminate them anyway.
That, roughly, is the logic of the whole process.
The General Case
Most textbooks limit their treatment of the method of characteristics to quasilinear PDEs. In fact, however, for a general first-order PDE
$$\begin{equation}F\left(\boldsymbol{x}, u, \frac{\partial u}{\partial\boldsymbol{x}}\right)=0\end{equation}$$
the method of characteristics still applies, where $F$ can be an arbitrary function of multiple variables.
This part of the discussion is mainly based on the English Wikipedia:
https://en.wikipedia.org/wiki/Method_of_characteristics
Derivation
To this end, let us first write
$$\begin{equation}\boldsymbol{p} = \frac{\partial u}{\partial\boldsymbol{x}}\end{equation}$$
Then, differentiating both sides of $F\left(\boldsymbol{x}, u, \boldsymbol{p}\right)=0$, we obtain
$$\begin{equation}\begin{aligned}0 =& \frac{\partial F}{\partial\boldsymbol{x}}\cdot\frac{d\boldsymbol{x}}{ds}+\frac{\partial F}{\partial u}\frac{\partial u}{\partial \boldsymbol{x}}\cdot\frac{d\boldsymbol{x}}{ds}+\frac{\partial F}{\partial\boldsymbol{p}}\cdot\frac{d\boldsymbol{p}}{ds}\\ &=\left(\frac{\partial F}{\partial\boldsymbol{x}}+\frac{\partial F}{\partial u}\boldsymbol{p}\right)\cdot\frac{d\boldsymbol{x}}{ds}+\frac{\partial F}{\partial\boldsymbol{p}}\cdot\frac{d\boldsymbol{p}}{ds} \end{aligned}\end{equation}$$
We notice that the above is the sum of two dot products of vectors that add up to zero. An interesting way to resolve this is to require
$$\begin{equation}\frac{d\boldsymbol{x}}{ds}=\frac{\partial F}{\partial\boldsymbol{p}},\quad \frac{d\boldsymbol{p}}{ds}=-\frac{\partial F}{\partial\boldsymbol{x}}-\frac{\partial F}{\partial u}\boldsymbol{p}\end{equation}$$
and additionally we have
$$\begin{equation}\frac{du}{ds}=\frac{\partial u}{\partial\boldsymbol{x}}\cdot\frac{d\boldsymbol{x}}{ds}=\boldsymbol{p}\cdot\frac{\partial F}{\partial\boldsymbol{p}}\end{equation}$$
Putting these together gives us the ODE system
$$\begin{equation}\left\{\begin{aligned}&\frac{d\boldsymbol{x}}{ds}=\frac{\partial F}{\partial\boldsymbol{p}}\\ &\frac{d\boldsymbol{p}}{ds}=-\frac{\partial F}{\partial\boldsymbol{x}}-\frac{\partial F}{\partial u}\boldsymbol{p}\\ &\frac{du}{ds}=\boldsymbol{p}\cdot\frac{\partial F}{\partial\boldsymbol{p}}\\ &F\left(\boldsymbol{x}, u, \boldsymbol{p}\right)=0\end{aligned}\right.\end{equation}$$
From here, the remaining steps are essentially the same as in the quasilinear case, except that we've introduced $n$ additional variables $\boldsymbol{p}$. Solving this system, we obtain integration constants together with $2n$ constants associated with $s$, and we determine the relations among the constants using the initial conditions. The difference is that, because of the extra $n$ variables $\boldsymbol{p}$, we also need to account for partial derivatives of the initial conditions, which makes the solving process more involved. See the example below.
Yet another example
$$\begin{equation}\frac{\partial u}{\partial t} = \left(\frac{\partial u}{\partial x}\right)^2,\quad u(x,0)=f(x)\end{equation}$$
that is, $0=F(\boldsymbol{p})=p_t - p_x^2$, so from equation $(18)$ we obtain the characteristic equations
$$\begin{equation}\left\{\begin{aligned}&\frac{dt}{ds}=1,\frac{dx}{ds}=-2p_x\\ &\frac{dp_t}{ds}=0,\frac{dp_x}{ds}=0\\ &\frac{du}{ds}=p_t - 2p_x^2 = -p_x^2\end{aligned}\right.\end{equation}$$
Taking $s=t$ directly, we find that $p_t=C_1,p_x=C_2$ are all constants, and hence $x=-2C_2 t + C_3, u = -C_2^2 t + C_4$. Next, using the initial condition, when $t=0$ we have $x=C_3,u=C_4$, which implies $C_4=f(C_3)$. Substituting this in gives
$$\begin{equation}u = -C_2^2 t + f(x + 2C_2 t)\end{equation}$$
Then, taking the partial derivative of the initial condition with respect to the variable $x$, we find that when $t=0$,
$$\begin{equation}C_2 = p_x = f'(x)=f'(C_3)\end{equation}$$
Note that $f$ is a predetermined initial-value function, so the equation above is in fact just an algebraic equation, giving us
$$\begin{equation}\begin{aligned}u =& -t\times f'(C_3)^2 + f\Big(x + 2t\times f'(C_3)\Big)\\ =& -t\times f'(C_3)^2 + f(C_3) \end{aligned}\end{equation}$$
and
$$\begin{equation}x=-2C_2 t + C_3=-2t\times f'(C_3) + C_3\end{equation}$$
Solving these for $C_3$ and substituting back into the previous expression yields the complete solution. For instance, when $f(x)=x^2$, we solve to get
$$\begin{equation}C_3=\frac{x}{1-4t}\end{equation}$$
and substituting gives
$$\begin{equation}u=\frac{x^2}{1-4t}\end{equation}$$
The top-down view
On the surface, the method of characteristics for the general case looks quite different from the quasilinear case. So can the general-case method be reduced back to the quasilinear one?
Indeed, substituting $F=\boldsymbol{\alpha}\cdot\boldsymbol{p}-\beta$ into equation $(18)$ gives
$$\begin{equation}\left\{\begin{aligned}&\frac{d\boldsymbol{x}}{ds}=\boldsymbol{\alpha}\\ &\frac{du}{ds}=\boldsymbol{p}\cdot\boldsymbol{\alpha}\end{aligned}\right.\end{equation}$$
We don't need to write out the equation for $\boldsymbol{p}$ separately, because we already have $\boldsymbol{p}\cdot\boldsymbol{\alpha}=\beta$, and substituting this into the equation above closes the system. This is exactly a reduction to equation $(5)$.
There's also a particularly simple special case where $F$ depends only on $\boldsymbol{p}$, in which case $\frac{\partial F}{\partial\boldsymbol{x}}$ and $\frac{\partial F}{\partial u}$ are both zero. Then, in the characteristic equations, $\boldsymbol{p}$ is constant, so $\boldsymbol{x}$ and $u$ are both linear functions of $s$, and the whole system becomes completely solvable. After that, it turns into a purely algebraic problem.
The Case of Systems of Equations
Can the technique of characteristics above be extended to systems of first-order PDEs? Generally, no — because solving general systems of first-order PDEs would be equivalent to solving PDEs of arbitrary order, and we clearly haven't seen such a general method exist (if it were possible, someone surely would have done it already).
That said, if the differential-operator part of the system is shared across all equations, the method of characteristics can still be used. Specifically, consider the PDE
$$\begin{equation}\left(\boldsymbol{\alpha}(\boldsymbol{x},\boldsymbol{u}) \cdot \frac{\partial}{\partial \boldsymbol{x}}\right) \boldsymbol{u} = \boldsymbol{\beta}(\boldsymbol{x},\boldsymbol{u})\end{equation}$$
Here $\boldsymbol{u}$ is also a vector, but the differential operator on the left-hand side is shared across components, while the components of $\boldsymbol{\beta}$ on the right-hand side can differ. In this case, the method of characteristics can still be applied to obtain
$$\begin{equation}\left\{\begin{aligned}&\frac{d\boldsymbol{x}}{ds}=\boldsymbol{\alpha}(\boldsymbol{x},u)\\ &\frac{d\boldsymbol{u}}{ds}=\boldsymbol{\beta}(\boldsymbol{x},u)\end{aligned}\right.\end{equation}$$
Of course, this is just a straightforward extension of the original method.
Translated automatically with claude-sonnet-5; all equations are reproduced verbatim from the source. Copyright remains with the original author.