Chapter 2
First-Order Partial Differential Equations

2.1 The Method of Characteristics

Figure 2.2

This figure in the book was generated by the command which is now known as PlotVectorField.  The appropriate command in the current version of Mathematica is

<<Graphics`PlotField`

r=Sqrt[x^2+y^2];
PlotVectorField[{x,y}/r,{x,-1,1},{y,-1,1}]

This plots the vector field in the x-y plane.  A three-dimensional representation of this vector field can be obtained with PlotVectorField3D.  The corresponding commands are

<<Graphics`PlotField3D`

r=Sqrt[x^2+y^2+z^2];
PlotVectorField3D[{x,y,z}/r,{x,-1,1},{y,-1,1},
     {z,-1,1},VectorHeads->True,PlotPoints->5]

Example 2.1

Mathematica can be used to show directly that the expression in Equation (2.37) solves the differential equation in (2.30) for any differentiable function F.  We first define the function f according to the general solution in (2.37):

f=F[x/y,x/u];

The general solution is now substituted into the differential equation in (2.30) and Simplify is used to show that an identity results for any function f:

Simplify[x D[f,x] + y D[f,y] + u D[f,u] == 0]

Analogous steps can be carried out to verify that the expression u=xG[x/y] solves Equation (2.27) for any function f that can be differentiated.

The transformation of (2.27) carried out after Equation (2.38) can also be performed within Mathematica.  The transformed function u is defined as
u[x,y]=xG[x/y,y] and substituted into the original differential equation in (2.27).  Simplify is then used to combine terms to obtain the transformed differential equation.

u=x G[x/y,y];

Simplify[x D[u,x] + y D[u,y] - u == 0]

Clear[u]

Example 2.3

The time development of the solution in (2.60) from a particular initial condition can be constructed with remarkable ease using Mathematica.  Suppose the initial condition in (2.54) is given by f[x]=1-x^2, if |x| is less than or equal to 1, and f[x]=0, if |x| is greater than 1.  This initial function has been chosen because the solution in (2.60) can be solved for u as a function of x and t (Problem 18).  However, the procedure described below for constructing the solution using Mathematica can be applied to initial functions that do not permit such an explicit solution.

The evolution of the point x=x0, u=f[x0] on the initial function at t=0 can be determined from (2.56) and (2.57).  At t=0, we have c2=x0, so at time t the
point initially at x0 has moved to

x = x0 + f[x0] t

These are the characteristics of Equation (2.53).

To construct the characteristics and the evolution of the solution in Mathematica, the initial profile (2.61) and the characteristics (2.62) are first entered as

f[x_]:=If[Abs[x]<1,1-x^2,0]

x[t_,x0_]:=x0+f[x0]t

Both the characteristics and the evolution of the solution can now be displayed
with ParametricPlot.  The trajectory of the characteristic emanating from the
point x0=0 up to time t=1 is displayed with

ParametricPlot[{x[t,0],t},{t,0,1}]

Similarly, the solution at time t over the region -2<x0<2 is displayed with the command

u[t_]:=ParametricPlot[{x[t,x0],f[x0]},{x0,-2,2}]

Notice that in generating the solution with ParametricPlot the quantity being
plotted is the pair of points {x[x0,t],f[x0]}, which corresponds to the position at time t of the point on the profile u[x,t] that was initially located at x0.  Thus, this construction represents a point-by-point assembly of the solution at time t in terms of the initial profile and the characteristics.

The trajectories of several characteristics projected on to the x-t plane and the early stages of evolution of the solution from the initial profile (2.61) are shown in Figure 2.7.  For short times the initial solution is deformed because the characteristics move with different speeds for different values of u.  At the point where two characteristics first cross, the solution becomes double-valued, and the slope in the u-x plane becomes infinite.  From this point onward, the solution is no longer a continuous function of position, and corresponds to the onset of a 'shock' wave.  This type of behavior is due to the form of the characteristics and not to a specific choice of initial condition

Figure 2.7

The top panel of Figure 2.7 can be generated by using the command Show:

u[t_]:=ParametricPlot[{x[t,x0],f[x0]},{x0,-2,2},
     DisplayFunction->Identity]

Show[{u[0],u[1],u[2]},AxesLabel->{x,u},
     DisplayFunction->$DisplayFunction]

The bottom panel is similarly generated by first defining the trajectory of the characteristic that emanates from the point x0 at time t=0 and ends at the point r at time t:

r[x0_,t_]:=ParametricPlot[{x[s,x0],s},{s,0,1},
     DisplayFunction->Identity]

A collection of characteristics is now assembled into Trajectories, which is then displayed using Show, as above:

Trajectories[t_]:=Table[r[x0,t],{x0,-1.5,1.5,.25}]

Show[Trajectories[1],AxesLabel->{x,t},
     DisplayFunction->$DisplayFunction]