Explaining the method of characteristics

The basic idea is that we look for parametric curves along which the PDE tells us how the function changes. Suppose you have a smooth parametric curve in the $xy$ plane: $x = X(t)$, $y = Y(t)$. Consider how a smooth function $u(x,y)$ changes as you move along the curve. The chain rule says $$\frac{du}{dt} = \frac{\partial u}{\partial x} \frac{dX}{dt} + \frac{\partial u}{\partial y} \frac{dY}{dt}$$ Now this looks rather like the left side of a first-order PDE $$a(x,y) \frac{\partial u}{\partial x} + b(x,y) \frac{\partial u}{\partial y} = c(x,y)$$ In fact, if you can find $X(t)$ and $Y(t)$ such that $\frac{dX}{dt} = a(X(t),Y(t))$ and $\frac{dY}{dt} = b(X(t),Y(t))$, this tells you that along that curve $\frac{du}{dt} = c(X(t), Y(t))$, so that if you know $u(X(0), Y(0))$ you can get $$u(X(s),Y(s)) = u(X(0), Y(0)) + \int_0^s c(X(t),Y(t))\ dt$$

Now, for any point $(x_1, y_1)$ in the plane, suppose you want to find $u(x_1, y_1)$, where $u(x,y)$ satisfies the PDE $a(x,y) \frac{\partial u}{\partial x} + b(x,y) \frac{\partial u}{\partial y} = c(x,y)$ plus some boundary condition. Then you want to find a curve $x = X(t)$, $y = Y(t)$ satisfying the system of ordinary differential equations $\frac{dX}{dt} = a(X(t),Y(t))$ and $\frac{dY}{dt} = b(X(t),Y(t))$ that passes through the given point $(x_1, y_1)$, follow that curve to some $(x_0, y_0)$ where your boundary condition tells you the value of $u$, and then you can get $u(x_1, y_1)$ by an integral as above.

This is an intuitive take which is intended to be as simple as possible.

Let's say you're measuring some system for two variables $x$ and $t$, and you discover that you can describe the system by $u_t+a(x,t)u_x=0$.

Okay, but what if $x$ and $t$ just superficially look like independent variables, but in fact $x$ is a function of $t$?

If it were so, then $$\frac{d}{dt}u(x(t),t)=u_t+x'(t)u_x$$ which means that if we can find a function $x(t)$ with the property that $x'(t)=a(x,t)$, then $u(x(t),t)$ would have the property that $$\frac{d}{dt}u(x(t),t)=u_t+a(x,t)u_x$$ which is exactly the type of function we are looking for.

So if it really is the case that our measured variable $x$ actually is a function of the underlying variable $t$, we only need to find a function $x(t)$ such that $x'(t)=a(x(t),t)$ with $x(0)=x_0$ and solve the much easier ODE $$\frac{d}{dt}u(x(t),t)=0$$ and then $u(x(t),t)$ satisfies the PDE we started with.

Usually these problems come with some boundary conditions like $u(x,0)=\phi(x)$, so we solve our easy ODE: $$u(x(t),t)=C$$

If we set $C=\phi(x_0)$ we see that $u(x(t),t)=u(x(0),0)=u(x_0,0)=\phi(x_0)$ (indeed, $u$ is constant wrt $t$)

So $$u(x(t),t)=\phi(x_0)$$

is a solution to the PDE provided we can find an $x(t)$ such that $x'(t)=a(x(t),t)$ with $x(0)=x_0$, which we can. When you find $x(t)$, you can solve it for $x_0$ and substitute so you get a solution function consisting of $x$'s and $t$'s, which you can then differentiate to check that it satisfies the PDE.

As an intuitive summary, a surface described by a first order pde inherently has a vector field defined on the surface (see Monge cones). So naturally the vector field provides an associated set of flow lines which cover the surface. Now along a single flow line, as we move down the line, we simply follow the vector field, so the line is described by an ode. The method of characteristics is finding the ode along the flow line.