First order PDEs

In this section, we describe a general technique for solving first-order equations. We start our discussion with solutions to the transport equation to demonstrate the geometric concept behind PDE solution methods. We then move on to semilinear, first-order PDEs to introduce the idea behind the method of characteristics. These ideas can then be extended to the more complicated PDEs given by quasilinear and fully nonlinear cases. Before we continue however, we expand on the classification of the PDEs discussed in Sec. 2.1 giving a broader definition of the various classes of PDEs we can encounter.

1. A PDE of order $n$ is quasilinear if it is linear in the derivatives of order $n$ with coefficients that depend on the independent variables and derivatives of the unknown function of order strictly less than $n$.
2. A quasilinear PDE where the coefficients of derivatives of order $n$ are functions of the independent variables alone is called semilinear.
3. A PDE which is linear in the unknown function and all its derivatives with coefficients depending on the independent variables alone is called linear.
4. A PDE which is not quasilinear is fully nonlinear.

The transport equation

The simplest transport equation takes the form,

where $c$ is a constant and $u=u(x, t)$. We seek functions $u(x, t)$ that satisfy the transport equation. We approach this by looking a solution method via the directional derivative. The directional derivative of $u$ in the direction of a vector $\boldsymbol{v}$ is given by:

where $\hat{\boldsymbol{v}}$ is the unit vector in the direction of $\boldsymbol{v}$. In particular, if the directional derivative in Eq. (2.19) is 0 at all points then $u$ is constant along all lines that are tangential to $\boldsymbol{v}$. This can be thought of as a generalisation of a partial derivative. Now, by rewriting the transport equation as a dot product, we have

Notice that the RHS of Eq. (2.20) is almost the directional derivative of $u$ in the direction of $\boldsymbol{v}$. To make it look exactly like the directional derivative, we can simply divide by $\sqrt{c^{2}+1}$, i.e.:

which is equivalent to Eq. (2.19). Hence, every solution $u$ must be constant along lines that are tangential to $(c, 1)$. In this case, the lines are also parallel to the vector $(c, 1)$ (see Fig. 2.1). These lines, since they are parallel to $(c, 1)$, they have slope,

Upon integrating, we obtain:

where $x_{0}$ is a constant. Rearranging to give

we can see that the equation describes infinitely many lines which are tangential and parallel to the vector $(c, 1)$. Along each of these lines $u$ is constant and, so, $u$ depends only on $x-c$. The general solution to the PDE therefore takes the form,

where $f$ is some arbitrary, differentiable function.

Method of characteristics

Consider the PDE:

where $u=u(x, y)$. Equation (2.24) is semilinear and it includes the linear form,

Figure 2.1: Characteristic curves; the PDE solution is constant along each curve in $(x, t)$ space and each curve is equal to a different constant value, $x_{0}$. The initial data is propagated along these characteristic curves.

with $f(x, y, u)=g(x, y)-c(x, y) u$. We also consider the initial condition given by:

Suppose we have found a solution $u(x, y)$ to Eq. (2.24). Plotting the relationship of the solution as it varies with the independent variables $x$ and $y$, gives a surface, $z=u(x, y)$. Let us denote this surface by $S=u(x, y)-z$, or equivalently,

for all real variables $x$ and $y$. A surface $u=u(x, y)$ is a solution to the PDE and is known as the integral surface. In Eq. (2.27), $u(x, y)$ is a function of $x$ and $y$ where $u$ is a variable. For instance, if $u(x, y)=x-y$ then $S=x-y-u$. From vector calculus, we know that the gradient, $\nabla S$ gives a normal to the surface $S=0$, i.e.:

This is a downward-pointing vector (because of the -1 term in the $z$-direction), as shown in Fig. 2.2. Note that the PDE (2.24) may be rewritten in the following dot product form:

Equation (2.29) is the scalar product of two vectors; a zero dot product indicates that the two vectors are at right angles. We can write Eq.

(2.29) as:

The vector,

is normal to $\boldsymbol{N}$. Now, since $\boldsymbol{N}$ is normal to the surface $S=0$, this means that $\boldsymbol{v}$ is tangential to the surface $S=0$ which, in turn, implies that $\boldsymbol{v}$ lies in the tangent plane to $S$ (see Fig. 2.2).

Figure 2.2: The surface $S$ showing the normal to the surface, $N(x, y)$ at a point $(x, y, u(x, y))$ and the vector $\boldsymbol{v}$ [given by Eq. (2.31)] which lies in the tangent plane to $S$.

As a result, the PDE dictates that any integral surface through a given point must be tangent to $\boldsymbol{v}$. This means that we start from an initial condition (which lies in the integral surface) and we move in the direction of $\boldsymbol{v}$; since $\boldsymbol{v}$ lies in the tangent plane, we move along a curve which lies within $S=0$. This curve is referred to as the characteristic curve. The question is how do we go about constructing such a surface? We start by describing points on a characteristic curve through the parametrisation,

Differentiating $\boldsymbol{r}$ with respect to $t$ gives the tangent vector

where the primes denote differentiation with respect to $t$. This tangent vector, $\boldsymbol{r}^{\prime}(t)$ then, belongs to the tangent plane to the surface $S=0$ at a given point. Recall that the characteristic curve moves in the direction of $\boldsymbol{v}$ and as such, $\boldsymbol{v}$ and $\boldsymbol{r}^{\prime}$ must be proportional. Therefore, we have

for some $\lambda$. Using Eqs. (2.31), (2.33) and (2.34), we can write down the following:

Or, in differential form:

From Eq. (2.36), we can form two pairs of ODEs:

or

The solutions to Eq. (2.37) [or Eq. (2.38)] determine the characteristics of the PDE. Example 2.1 Find a solution to the following semilinear PDE,

Solution We start off by writing the ODEs as in Eq. (2.36):

The objective here is to find a pair of differential equations which are easy to solve with basic ODE techniques. As mentioned previously, we have the two options given by Eqs. (2.37) and (2.38). We see that if we choose the first pair [Eq. (2.37)], we have:

which we can easily solve to obtain,

The second ODE in the pair yields:

which, upon integration, gives:

Now, remember that we are solving these ODEs along characteristic curves. Equation (2.44) gives us $u$ through $c_{2}$ (note that $c_{2}$ varies between characteristics but it remains constant along the same characteristic). The particular characteristic curve we are considering is given in terms of the independent variables $x$ and $y$ whose relationship is dependent on the constant $c_{1}$. For this reason, $c_{1}$ and $c_{2}$ are related. At this point we don't exactly know how they are related so we introduce some arbitrary (but differentiable) function for their relationship:

Hence, using Eqs. (2.42), (2.44) and (2.45), we have:

which is an implicit, general solution to Eq. (2.55). In Example 2.1, the solution to the PDE is given in a general form. To determine the functionality of $g$, we need to apply the initial condition given by $u(x, 0)=u(x)$ where $u$ is a function of a single variable.

Summary of the method

Suppose we have a semilinear PDE of the general form given by Eq. (2.24). To solve using the method of characteristics:

• Consider the system:

The equations in (2.47) are equivalent to the PDE.

• Solve a pair of ODEs to obtain the functions $h(x, y)=c_{1}$ and $j(x, y, u)=$ $c_{2}$.
• In particular, solve the following ODE to obtain the function $h(x, y)=c_{1}$,

• And the following for $j(x, y, u)=c_{2}$

• Relate the constants $c_{1}$ and $c_{2}$ through an arbitrary but differentiable function, $g$, i.e. $j(x, y, u)=g(h(x, y))$.
• The function $g$ can be chosen to satisfy an initial condition of the form $u(x, 0)=u_{0}(x)$ where $u_{0}$ is a function of a single variable.

Quasilinear equations

Next, we look at quasilinear PDEs. Note that by Definition 2.2, a quasilinear PDE may also be linear or nonlinear since the definition places no restriction on the function $u$ itself. The method outlined above can be applied to quasilinear equations as well. Next, we walk through an example as the method can be a little trickier to use in more complicated cases.

Consider the quasilinear PDE of the following form,

Following Eq. (2.47), we have:

From Eqs. (2.51), if we can set up two first-order ODEs [i.e. analogous to the pairs in Eqs. (2.48) or (2.49)], then it is possible to obtain two equations of the form:

Just like in the case of semilinear PDEs (and as seen in Example 2.1) the constants $c_{1}$ and $c_{2}$ are related and hence the solution to the PDE (2.50), is given by:

where $F$ is an arbitrary (but differentiable) function. Example 2.2 Solve the following quasilinear PDE,

Solution Using Eq. (2.36), we write:

Recall the objective is to form first-order ODEs in the form given by Eqs. (2.48) and (2.49). The algebra in more complicated first-order PDEs tends to be more tedious. Here, we make use of the following:

Then, we have:

We now have the following ODE:

Rearranging gives:

which is equivalent to:

Similarly, we have:

which is:

Finally, we relate the two constants through Eq. (2.53), i.e. $c_{2}=F\left(c_{1}\right)$, to obtain the PDE solution as:

where $F$ is an arbitrary, differentiable function.

Traffic flow

We now consider an application, namely one-dimensional and one-directional traffic flow. We consider cars on a single-lane road with no entry or exit ramps such that we have conservation of cars. We take $x$ to denote a unique point on the road at time $t$. While cars are discrete objects, we model this in a continuum sense using the notion of traffic density. To ensure that the model is valid, we assume that we are dealing on a large length scale. Then, since car flow is conserved, we can model the traffic flow with the following continuity equation,

The traffic density is denoted by $\rho(x, t)$ and it represents the number of cars per unit length at position $x$ and time $t$. The traffic flow or flux is given by $q(x, t)$ at position $x$ and time $t$ and it represents the number of cars passing a fixed point in $x$ per unit time. Before we proceed, we need a constitutive relation between $\rho$ and $q$; we will get one by considering the speed at which the cars move. Suppose this is given by $c$ but this quantity dpes not represent a constant speed. Of course the speed at which cars move on the freeway is a factor of many things but, keeping things simple, let us say that this is predominantly be affected by the traffic density $\rho$ such that $c=c(\rho)$. If the traffic density were very low (i.e. very close to zero) then, the cars would be allowed to move at a maximum (hopefully, legal!) speed, say $c_{0}$. The traffic density cannot reach infinity as this would be physically impossible and must therefore be bounded by a maximum value; this may be represented by how many cars can fit on an arbitrary length of the freeway when they are closely packed (i.e. bumper-to-bumper) together. As the traffic density approaches this maximum value, the speed $c$ approaches zero. Suppose this maximum value is denoted by $\rho_{\max }$. It follows that the relation we are looking for (i.e. $c \rightarrow c_{0}$ as $\rho \rightarrow 0$ and $c \rightarrow 0$ as $\rho \rightarrow \rho_{\max }$, may be of the following form,

Now, since $q$ represents the flux, it represents the amount of cars per unit time per unit area. This is therefore given by:

We take $c_{0}=1$ and $\rho_{\max }=1$ for simplicity and, by differentiating Eq. (2.66) with respect to $x$, we have:

Hence, Eq. (2.64) becomes,

This is a quasilinear PDE in $\rho(x, t)$ so we may apply the method of characteristics. Remember these are curves in the $x t$-plane on which solutions to Eq. (2.68) are constant. If we now apply Eqs. (2.37) [note: we could of course apply Eq. (2.38) instead], we have the following ODEs:

Recall that $\rho=\rho(x, t)$ and so the first equation in (2.69) does not help us significantly in finding the characteristics; however, the second equation saves the day. The second equation, $\frac{d \rho}{d t}=0$ tells us that the value of $\rho$ along the characteristics curves is constant, say $\rho_{0}$ which depends on the value of $x_{0}$ at time $t=0$. Thus, the first equation in Eq. (2.69) becomes,

with $\rho_{0}$ constant which integrates to:

Equation (2.71) implies that the characteristics are straight lines with $x_{0}$ being the starting point. Different initial conditions therefore yield different characteristic curves (i.e. different slope). Since we typically plot $t$ against $x$, it is handy to have an equation that relates $t$ to $x$; from Eq. (2.71), we have:

Figure 2.3 shows various characteristic curves; on the leftmost plot, we show a single curve passing through $x_{0}$. Now, $\rho\left(x_{0}, 0\right)$ represents the initial distribution of the traffic density which varies according to the value of $x$ at $t=0$. Let us call this function $f(x)$ i.e. $\rho(x, 0)=f(x)$. At two distinct points, say $x_{0}$ and $x_{1}$, the value of the density function is $f\left(x_{0}\right)$ and $f\left(x_{1}\right)$, respectively. Suppose that $f\left(x_{0}\right)>f\left(x_{1}\right)$. Graphically, this implies that the slope of the characteristic at $x_{0}$ is less steep than that of the characteristic at $x_{1} 1$ (see the middle plot in Fig. 2.3). If the value of $f$ decreases continuously as we move in the direction of increasing $x$, the slopes (i.e. $d x / d t$ ) increase producing a somewhat 'fan'-like structure as shown in Fig. 2.3. Since the traffic density is less at $x_{1}$ than at $x_{0}$, the traffic is moving faster at $x_{1}$. It follows that there is a gradual separation between the traffic that started at $x_{0}$ and at $x_{1}$.

Figure 2.3: Characteristic curves in the $x t$-plane. The curve emanates from the $x$-axis satisfying a starting point $x_{0}$ (left). Two characteristic curves starting from initial points $x_{0}$ and $x_{1}$ (middle). Gradual separation between the points $x_{0}$ and $x_{1}$, showing a 'fan'-like pattern (right).

A question that arises is the following: what if the traffic density starts low at $x_{0}$ and increases as $x$ increases? Then, the respective characteristic curves will be pointing inward and toward each other. The fate of these curves will ultimately be that they intersect. Now, recall that the model involves a single-lane road and therefore cars cannot pass other cars. If the speed is less at $x_{1}$ compared to the speed at $x_{0}$, then cars are approaching the cars in front of them and are forced to reduce their speed to the lower speed the cars in the front maintain. We will see the effect this has on the solutions in the next subsection.

Finally, to conclude this subsection, we note that in order to compute the solution $\rho(x, t)$ at any point, we find the characteristics using Eq. (2.72), follow it back to the initial point $\left(x_{0}, 0\right)$ and determine $\rho(x, t)=f\left(x_{0}\right)$. For the model we looked at in this section where we have heavier traffic to the left than to the right, a possible initial condition given by $\rho(x, 0)=f(x)$ may be:

We make the following observation: recall that the maximum density was scaled by $\rho_{\max }=1$, therefore the initial constant density for $x<0$ is merely a quarter of the maximum. It is possible that the characteristic curves intersect at some point in time even if light traffic heads into heavier traffic.

${ }^{1}$ Note that the slope is defined as the change of $x$ with $t$; however, we are plotting $t$ against $x$ and so in Fig. 2.3 the slope appears to be less at $x_{1}$ compared to $x_{0}$.

Shocks

As mentioned in the previous subsection, it is possible that more than one characteristic passes through a given point in the $x t$-plane (note: this is a consequence of the nonlinear term given by $\rho \rho_{x}$ ). This is in fact a very likely scenario leading to the overlap of the characteristics for sufficiently large time. The point at which characteristics first touch is called a shock wave or, simply, a shock. The term comes from gas dynamics where this phenomenon was first encountered.

Consider the traffic flow equation [given by Eq. (2.68)] with the following initial condition,

This represents a case where there is fast traffic coming from behind, there's a relatively slow moving traffic at the very front and a transition region where the traffic is slowing down - this is in $0<x<1$. For such a case, at some point in time, the characteristics overlap (see Fig. 2.4). Once the characteristics intersect, the density function is multi-valued at the intersection point since the density function is composed of points where $\rho=f\left(x_{0}\right)$ and $\rho=f\left(x_{1}\right)$. Now,

Figure 2.4: Characteristic curves in the $x t$-plane. The characteristics are pointing inward as faster-moving traffic from behind is catching up to slower-moving traffic in the front.

what happens in terms of the traffic density solution? As Eq. (2.183) indicates the traffic density starts low and gradually approaches the value of 0.25 within $0<x<1$; thereafter, the density remains constant at 0.25 at time $t=0$. As time progresses, the solution gets steeper and steeper until it becomes vertical which coincides with when the shock wave develops. After the shock develops, the solution in the form we derived in in the previous section becomes invalid: it becomes a multi-valued function which poses a problem as we don't know which value to assign at the solution where the curves intersect. So what can be done? Two things: either we modify the model (which will consequently give a different solution) such that shocks are prevented or, modify the solution procedure. We proceed with the latter.

We give an alternative definition for the 'solution' to the PDE. More specifically, we allow the solution to become vertical at some point in time. This implies that at that point in time, the solution has a jump discontinuity (see notes from Topic B1 for definition). The solution at that point is no longer differentiable and hence cannot be a solution to the PDE in the classical sense. We can find ways around this though. We can think of a solution with a jump discontinuity as a shock, travelling in the direction of increasing $x$ as time increases (see Fig. 2.5). This 'solution' should be valid if we can show that the speed of the shock does not violate the principle of conservation (recall that in deriving the model we assume that the number of cars entering and leaving a certain length of the freeway, is constant). We can therefore divide the domain in two regions: upstream and downstream of the shock. A solution to a PDE that does not need to be continuous is referred to as a weak solution. A strong or classical solution is one that it is defined by a continuous function.

Figure 2.5: Shock formation. At $t=t_{1}$ a shock forms with a jump discontiniuity at $x_{s}\left(t_{1}\right)$. At a later time, $t_{2}$, the jump discontinuity is at $x_{s}\left(t_{2}\right)$.

Suppose that a shock develops and moves with position $x_{s}(t)$ and now, the solution $\rho(x, t)$ has a jump discontinuity at $x_{s}$ resulting in the values of the traffic density at $x_{s}-\epsilon$ and at $x_{s}+\epsilon$, (where $\epsilon>0$ and small) to be finite but not equal. We introduce the following notation: $x_{s_{-}}=x_{s}-\epsilon$ and $x_{s_{+}}=x_{s}+\epsilon$. In particular, we have:

Away from the shock, the solution is described by a smooth function, $\rho(x, t)$. To determine the position and hence the speed of the shock wave, we demand the following:

1. The original PDE describing the model, i.e. Eq. (2.68) needs to be satisfied on either side of the shock: at $x<x_{s}$ and $x>x_{s}$ for all time, $t$. So, using the characteristics, we construct information on both sides of the shock.
2. As previously mentioned, the total flow, $q=(1-\rho) \rho$ needs to be conserved, relative to the moving shock. The flow from the left into the shock must be equal to the flow from right away from the shock. The following theorem ensures that this is satisfied if $\rho(x, t)$ is to be a weak solution to the PDE.

Theorem: Rankine-Hugoniot condition

If $\rho(x, t)$ is a weak solution to the quasilinear traffic flow PDE (2.68), such that $\rho$ is discontinuous across the curve $x=x_{s}(t)$ but it is smooth on either side of $x_{s}$, then $\rho$ must satisfy the following condition:

where $x_{s}^{\prime}$ is the speed of the shock wave and $q$ evaluated at $x_{s_{-}}$and $x_{s_{+}}$ denotes the flux into and out of the moving shock, respectively. Equation (2.76) is known as the Rankine-Hugoniot formula for the shock speed.

For the traffic flow model described in this set of notes, the flux is given by $q=(1-\rho) \rho$. Now, let $\rho_{-}=\rho\left(x_{s_{-}}, t\right)$ and $\rho_{+}=\rho\left(x_{s_{+}}, t\right)$. Then, Eq. (2.76) implies that the shocks move at speeds which obey:

which implies that the shock speed for the traffic equation is determined by the density value on either side of the shock.