Description

Please solve all of these questions and show all the steps in the Words. If you provide me the photo or picture, make sure it will be clear. Had better send me one answer when you finish one. I provide the materials which is the PDF

28 attachmentsSlide 1 of 28attachment_1attachment_1attachment_2attachment_2attachment_3attachment_3attachment_4attachment_4attachment_5attachment_5attachment_6attachment_6attachment_7attachment_7attachment_8attachment_8attachment_9attachment_9attachment_10attachment_10attachment_11attachment_11attachment_12attachment_12attachment_13attachment_13attachment_14attachment_14attachment_15attachment_15attachment_16attachment_16attachment_17attachment_17attachment_18attachment_18attachment_19attachment_19attachment_20attachment_20attachment_21attachment_21attachment_22attachment_22attachment_23attachment_23attachment_24attachment_24attachment_25attachment_25attachment_26attachment_26attachment_27attachment_27attachment_28attachment_28

Unformatted Attachment Preview

Linear and Nonlinear Waves

1

1.1

Introduction

Wave

A wave is an identifiable (measurable) signal that propagates with a definite

velocity. This signal may be any feature provided that it is identifiable and

its location can be determined at any time. The signal may change its shape,

size and velocity, provided that it remains recognizable.

Examples: Sound waves, electromagnetic waves (light, radio waves, infra-red,

X-rays etc), water waves, waves on strings, elastic waves (earthquakes etc),

chemical waves, biological waves, traffic waves etc.

1.2

Classes of Waves

It is convenient to split waves into two main classes: hyperbolic waves and

dispersive waves.

1.2.1

Hyperbolic Waves

These are governed by hyperbolic partial differential equations. The prototype of equations of this kind is the linear wave equation.

2

∂ 2φ

2∂ φ

=

c

0

∂t2

∂x2

(1.1)

where φ = φ(x, t) is the dependent variable, x, t are the independent variables

and c0 is a constant.

This equation can be factored to give

1

∂

∂

+ c0

∂t

∂x

∂

∂

− c0

∂t

∂x

φ = 0.

(1.2)

The solution to (1.1) is therefore of the form

φ = f (x, t) + g(x, t)

where f, g satisfy

∂

∂

+ c0

∂t

∂x

f = 0 a),

∂

∂

− c0

∂t

∂x

g = 0. b).

(1.3)

Consider equation (1.3a). This is

∂f

∂f

+ c0

= 0.

(1.4)

∂t

∂x

The solution to this is f (x, t) = f (x − c0 t), where f is an arbitrary function.

To see this, we have

∂f

∂f

= −c0 f 0 ,

= f 0,

∂t

∂x

where f 0 is the derivative of f w.r.t. its argument. Substituting this into

(1.4) gives

−c0 f 0 + c0 f 0 = 0.

Similarly, the solution to (1.3b) is g(x, t) = g(x + c0 t). The complete

solution to (1.1) is therefore

φ(x, t) = f (x − c0 t) + g(x + c0 t)

(1.5)

We can see that if we travel along the x-axis with speed c0 , then x − c0 t =

const and so the solution φ = f (x − c0 t) is a wave that travels with speed

c0 without change of shape. Similarly, φ = g(x + c0 t) is wave travelling in

the negative x direction. Since they travel without change of shape, they are

perfect signals.

2

1.2.2

Dispersive Waves

The prototype of a dispersive wave is based on a type of solution rather than

a type of equation. A linear dispersive system is one that admits solutions

of the form

φ = A cos(kx − ωt)

(1.6)

where A is the amplitude), ω the frequency and k the wavenumber.

We have

3

Wavelength λ =

2π

, a)

k

Period T =

2π

, b).

ω

(1.7)

A solution of the form (1.6) is called a monochromatic wave (single frequency).

Example 1.1

The linear wave equation (1.1) has a solution of the form (1.6). Put (1.6) in

(1.1) to get

−ω 2 A cos(kx − ωt) = −c20 k 2 A cos(kx − ωt)

Hence

ω 2 = c20 k 2 → ω = ±c0 k.

(1.8)

As we already know, there are two waves travelling in opposite directions

with speed c0 . The speed, c0 = ω/k is called the phase speed. Note that

waves of all frequencies travel with the same speed i.e. this is not a dispersive

system.

In a dispersive system ω depends upon k

ω = W (k).

(1.9)

This is called the dispersion relation for the system. The phase speed of the

wave is

cp = ω/k = W (k)/k

(1.10)

and is not the same for all k for dispersive waves.

Example 1.2 The beam equation is

4

∂ 2φ

2∂ φ

+

γ

=0

∂t2

∂x4

(1.11)

4

Substituting (1.6) into this gives

−ω 2 A cos(kx − ωt) + γ 2 k 4 A cos(kx − ωt) = 0,

so that

ω 2 = γ 2 k 4 → ω = ±γk 2 .

(1.12)

The phase speed of the wave is now

cp = ω/k = ±γk

(1.13)

i.e. it is not the same for all k.

2

Hyperbolic Waves

Hyperbolic waves are described by hyperbolic equations, such as the linear

wave equation (1.1) discussed in the introduction. That equation is linear

(the wave speed does not depend upon φ) and it has only two waves with

speeds ±c0 . More generally, there can be a number of speeds and if these depend upon the solution, then the equation is nonlinear. The most important

property of such equations is that they define curves called characteristics

which can be used to determine the solutions.

2.1

Characteristics

Consider one of the factors in the linear wave equation (1.1)

∂φ

∂φ

+ c0

= 0.

∂t

∂x

(2.1)

This is called the linear advection equation because it describes the advection

of a passive scalar e.g. a contaminant in a fluid flow that does not diffuse or

react.

Consider straight lines for which

5

dx

= c0 .

dt

(2.2)

Along such a line we have

∂φ

∂φ

∂φ

∂φ

dφ =

dt +

dx =

dt +

c0 dt =

∂t

∂x

∂t

∂x

∂φ

∂φ

+ c0

∂t

∂x

dt = 0

since φ satisfies (2.1). The equation tells us that φ is constant along any line

that satisfies (2.2). These lines are called the characteristics of the equation.

The solution to (2.2) is x − c0 t = constant, which is consistent with the

solution φ = f (x − c0 t) which we obtained in section 1.2.1 (equation 1.5).

2.2

Nonlinear Scalar Hyperbolic Equation

(2.1) is the simplest linear hyperbolic equation. The simplest nonlinear hyperbolic equation is the scalar equation

∂φ

∂φ

+ c(φ)

= 0,

∂t

∂x

(2.3)

where the wave speed now depends upon the solution φ. Here scalar just

means that there is only one dependent variable φ.

Equations of this type arise in the study of traffic flow, flood waves,

chromatography, glacier flows etc.

Even though (2.3) is a nonlinear equation, we can still solve it using

characteristics. Given a solution φ to (2.3), consider the curves defined by

dx

= c[φ(x, t)].

dt

(2.4)

Along such a curve we have

∂φ

∂φ

∂φ

∂φ

dt +

dx =

dt +

cdt =

dφ =

∂t

∂x

∂t

∂x

6

∂φ

∂φ

+c

∂t

∂x

dt = 0

since φ satisfies (2.3). As for the linear equation, we have that φ is constant

along any curve that satisfies (2.4). Since φ is constant along these curves,

so is dx/dt and the curves are actually straight lines, just as in the linear

case. The only difference is that the lines do not all have the same slope. As

for the linear equation, these lines are the characteristics of the equation.

We now show how the characteristics can be used to determine φ(x, t)

given the initial data

φ(x, 0) = f (x),

where f is a given function.

Let C be one of the characteristics (2.4). If C intersects t = 0 at x = ξ

then φ = f (ξ) on the whole of C. Therefore, from

dx

= c[f (ξ)]

dt

we conclude that C = Cξ is given by

Cξ :

x = ξ + F (ξ)t

with F (ξ) = c(f (ξ)) .

(Note that F is known from the initial data.)

As a result, we have

φ(x, t) = f (ξ) ,

x = ξ + F (ξ)t ,

(2.5)

F (ξ) = c(f (ξ)) .

(2.6)

Allowing ξ to vary, we can obtain the solution φ(x, t) by first finding ξ from

(2.6) and then substituting it into (2.5).

Geometrically, to determine the solution at any point (x, t), we look for

the characteristic from (x, t) to the x-axis. Since the characteristic has slope

1/c[φ(x, t)], it intersects the x-axis at

x0 = x − tc[φ(x, t)] ,

at which point we have φ = f (x0 ).

Since φ is constant along the characteristic, we have φ(x, t) = f (x0 ), so

applying f to both sides of the previous formula, we obtain

φ = f (x − t · c[φ]) ,

where f is the initial data φ|t=0 ,

7

(2.7)

and this gives us an implicit equation for determining φ = φ(x, t).

The basic definition of a wave is that some recognizable feature moves

with a finite speed. For hyperbolic equations, characteristics correspond to

this idea. Information about the solution is “carried” along the characteristics.

2.2.1

Example: Rarefaction Wave

Consider the equation (2.3) and suppose that

dc

> 0.

dφ

(An equation for which this derivative always has the same sign is called

convex). Let the initial data be

φ(x, 0) = f (x) =

φr x > 0

φl x < 0,
with φr > φl , so that cr = c(φr ) > cl = c(φl ).

We can see from the diagram that in the region x > cr t we must have

φ = φr since all the characteristics come from x > 0 at t = 0. Similarly, in

the region x < cl t we must have φ = φl . In the region cl t ≤ x ≤ cr t, the
characteristics must have come from the origin, so they must have slope x/t.
9
Denoting by c(x, t) =
therefore have
cl
x/t
c(x, t) =
cr
c[φ(x, t)] the slope of the characteristic at (x, t), we
for x < cl t
for cl t ≤ x ≤ cr t
for x > cr t.

(2.8)

Since we know the c(x, t), we can determine φ(x, t) by solving c(x, t) =

c(φ(x, t)). Such solutions are called rarefaction waves because the propagating wave produces a decrease in c(φ). In this case it is called a centred

rarefaction because there is a “fan” of characteristics emanating from a “centre” (the origin in this case).

2.3

Breaking of the solutions

The formula (2.7) resembles the formula φ(x, t) = f (x − c0 t) that we had in

the linear case of constant c = c0 , and it can be used to visualize the solution.

When c = const, the initial profile φ(x, 0) = f (x) moves with constant speed

c0 , and its shape remains the same. While for nonconstant c(φ), different

points on the graph of f (x) would move horizontally with various speed

c(φ), thus distorting the initial shape.

10

E.g., in the case when c is positive and increasing function, the upper

parts of the graph move faster than the bottom parts. As a result, the

solution φ may break, i.e. become multivalued, as a function of x.

This can also be observed on the (x, t) diagram as the region where the

characteristics overlap. The breaking point correspond to the lowest point

of the overlapping region. To determine its time and location, we notice

that the boundary of the overlapping region serves as the envelope for the

characteristics: the characteristics are tangent to it.

11

Let us consider two neighbouring characteristics Cξ and Cξ0 , with ξ 0 close

to ξ. Their intersection point is determined from the system of equations

x = ξ + F (ξ)t

x = ξ 0 + F (ξ 0 )t

From which we get (by subtracting and dividing by ξ 0 − ξ):

x = ξ + F (ξ)t

F (ξ 0 ) − F (ξ)

t

0 = 1+

ξ0 − ξ

In the limit ξ 0 → ξ, the intersection point of two characteristics tends to a

point on the envelope, so taking the limit we obtain the system describing

the envelope as

x = ξ + F (ξ)t

0 = 1 + F 0 (ξ)t

We can view this as a parametric equation for the envelope, rewriting it as

t=−

1

F 0 (ξ)

,

x=ξ−

F (ξ)

.

F 0 (ξ)

(2.9)

(Recall that F (ξ) = c(f (ξ)).)

The breaking time corresponds to the smallest possible t > 0 in (2.9),

i.e. where F 0 (ξ) takes its most negative value. This gives the formula for the

breaking time:

tb = −

1

,

min F 0 (ξ)

F (ξ) = c(f (ξ)) .

(2.10)

(The breaking location x can then be found from (2.9).)

We thus conclude that breaking does not occur if F 0 (ξ) > 0 for all ξ, i.e.

when c increases on the initial data.

As a final remark, we note that in applications φ usually corresponds to

some physical quantity (e.g. density or pressure) and cannot be multivalued.

Therefore, after t = tb the solution is no longer adequate and discontinuities

(shocks) appear. We will turn to this problem later, after we discuss how

equation (2.3) appears in modelling traffic flow.

12

3

3.1

Application: Traffic Flow.

Hyperbolic conservation law

Consider a very long stretch of road with no entries or exits (e.g. a section

of a motorway). Let ρ(x, t) be the number of cars per unit length and v(x, t)

the speed of the cars at time t, position x. We can then define a flux of cars

q(x, t) = ρv ,

where q is the number of cars passing position x in unit time. Since there

are no entries or exits, the number of cars is conserved. Hence the rate of

change of the total number of cars in any section x1 ≤ x ≤ x2 must be due

to cars entering and leaving the section, i.e.

Z

d x2

ρ(x, t) dx = q(x1 , t) − q(x2 , t) .

(3.1)

dt x1

Since x1 , x2 are fixed, we can write this as

Z x2

∂ρ ∂q

+

dx = 0 ,

∂t ∂x

x1

(3.2)

provided that the partial derivatives of ρ, q exist (smooth solutions). Now

x1 , x2 are arbitrary, so (3.2) can only be true for all such intervals if the

integrand vanishes everywhere. Hence

∂ρ ∂q

+

= 0.

∂t ∂x

(3.3)

As it stands, the problem is incomplete because we have two unknowns,

ρ and q (or v). To close the system, we assume that v and hence q depend

upon the local density of cars, ρ, i.e. q = q(ρ).

Then (3.3) can be written

∂ρ dq ∂ρ

+

= 0,

∂t dρ ∂x

(3.4)

which is just (2.3) with φ = ρ and c(ρ) = dq/dρ.

The equation (3.3) is an example of a hyperbolic conservation law. The

number of cars, ρ, is the conserved quantity and q is the associated flux.

13

Equations of this type also arise in flood waves, chromatography, glacier

flows etc.

Let us make a few remarks on the form of q = q(ρ) for the traffic flow.

Clearly, the speed of the cars should be a decreasing function of the density of

cars: the speed is a maximum when ρ → 0 (empty road) and should be zero

when ρ = ρmax , the value at which the cars are bumper to bumper (traffic

jam). The flow rate q = ρv must be zero at both ρ = 0 and ρ = ρmax and

must have a maximum, q = qm , at some intermediate value of ρ = ρm e.g.

q(ρ) = aρ log

ρmax

ρ

a constant

Observations of traffic for a single lane show that ρmax = 225 vehicles

per mile, ρm = 80 vehicles per mile and qm = 1500 vehicles per hour. The

q

maximum

flow rate occurs for a velocity vm = qm /ρm = 1500/80 = 20 mph.

qm

Slope =

ρm

ρ max

dq

d ρ

=c

ρ

The wavespeed is

c=

dq

d

dv

= (ρv) = v + ρ .

dρ

dρ

dρ

(3.5)

Since dv/dρ < 0, the wave speed is less than the car velocity. This means
that waves propagate backward through the traffic and drivers are warned
of disturbances ahead.
14
For ρ < ρm c(ρ) > 0 Waves move forward relative to road,

For ρ = ρm c(ρ) = 0 Waves are stationary relative to road.

For ρ > ρm c(ρ) < 0 Waves move backward relative to road.
3.2
Example
Now consider a set of traffic lights with a long enough red period for the cars
behind the lights to be stationary with ρ = ρmax i.e. bumper to bumper and
for there to be no cars ahead of the lights, ρ = 0. The initial conditions are
ρ(x, 0) =
ρl = ρmax for x < 0,
ρr = 0
for x > 0,

Note that cl = c(ρmax ) < 0 and cr = c(0) > 0. At t = 0, the lights change to

green. As in Example 2.2.1, we get a rarefaction wave because the wavespeed

is smaller on the left than the right (cl < cr ).
t
x = c lt
(x,t)
x=c r t
On these
lines
ρ = ρl
On these
lines
ρ = ρr
x
Note that the head of the wave, x = cr t, has a positive velocity and the
tail, x = cl t, has a negative velocity with respect to the road. At x = 0, we
must have c = 0, therefore, ρ = ρm ≈ 20 mph. This means that q = qm at
x = 0 i.e. the flow rate attains its maximum at the traffic lights. To
determine the solution, we use (2.8) and find ρ(x, t) from c(x, t) using the
relation c = c(ρ).
15
4
Shock solutions
4.1
Shock relation
In the previous example breaking did not occur, but what if it happens (e.g.,
when the lights turn red again)? Since ρ, the density of traffic, cannot be
multivalued, common sense tells us that there will be a shock (disturbance)
in the traffic, so ρ becomes discontinuous rather that multivalued.
Recall that (3.4) was derived from the conservation law (3.1), which is
more fundamental and has a meaning even if ρ has a discontinuity/jump in
value. Thus, it makes sense to use (3.1) to determine the shock position.
Suppose that we have a shock whose position is x = xs (t). At x = xs we
have a discontinuity in ρ and let ρ = ρl on the left of the shock, x → x−
s,
and ρ = ρr on the right, x → x+
.
We
cannot
use
the
differential
equation
to
s
relate ρl and ρr , but we can use the conservation law, (3.1).
Consider an interval, x1 ≤ x ≤ x2 which contains the shock at time t i.e.
x1 ≤ xs (t) ≤ x2 . We apply (3.1) to this interval and split the integral over x
into the parts on the left and right of the shock.
d
dt
Z
x2
x1
d
ρ(x, t) dx =
dt
"Z
x−
s
Z
ρ dx +
x+
s
x1
#
x2
ρ dx = q(x1 , t) − q(x2 , t). (4.1)
+
Note that limits of these integrals are x1 to x−
s and xs to x2 . We have the
following rules for differentiating an integral with variable limits:
d
dt
Z
x−
s
Z
x−
s
ρ dx =
x1
x1
∂ρ
dxs
dx +
ρl ,
∂t
dt
d
dt
Z
x2
Z
x2
ρ dx =
x+
s
x+
s
dxs
∂ρ
dx −
ρr .
∂t
dt
Since the limits of integration exclude the shock, the differential equation,
(3.4), is valid in the integrals and we can replace ∂ρ/∂t by - ∂q/∂x to get
"Z −
#
Z x2
xs
d
ρ dx +
ρ dx
dt x1
x+
s
Z x−s
Z x2
∂q
dxs
∂q
dxs
=−
dx +
ρl −
dx −
ρr
dt
∂x
dt
x1 ∂x
x+
s
dxs
= q(x1 , t) − ql − q(x2 , t) + qr +
(ρl − ρr ) .
dt
16
+
Here ql = q(x−
s ) = q(ρl ), qr = q(xs ) = q(ρr ). By the conservation law (3.1),
the above expression equals q(x1 , t) − q(x2 , t). We therefore have
ql − qr = q(ρl ) − q(ρr ) = s(ρl − ρr )
(4.2)
where
s=
dxs
dt
(4.3)
is the speed of the shock. These are the shock relations that we need to
determine the solution.
Note that if we let the magnitude of the discontinuity tend to zero, then
dq
q(ρl ) − q(ρr )
= c(ρl )
=
s = lim
ρl →ρr
ρl − ρr
dρ ρ=ρl
i.e. very weak shocks travel at the local wave speed.
4.2
Example
As in Example 2.2.1, let us consider equation (2.3) with dc/dφ > 0. Let the

initial data be

φr x > 0

φ0 (x) =

φl x < 0,
as before, but now with φl > φr ⇒ cl > cr . The (x, t) diagram is now

We see that the characteristics now overlap in the sector cr t < x < cl t,
which means that a shock appears at t = 0 and x = 0. Let x = xs (t) be the
shock position at time t.
Notice that the characteristics intersect at the shock to give the two values, φl and φr on the left and right of the shock. The shock relations, (4.2),
tell us the shock speed must be
s=
q(φl ) − q(φr )
,
φl − φr
which is constant. Therefore, the shock moves with uniform speed and
xs (t) = st.
17
t
On these
lines
φ=φl
On these
lines
φ=φr
x
t
Shock
On these
lines
φ=φl
On these
lines
φ=φr
x
The solution, therefore, is
φr for x > st

φ( x, t) =

φl for x < st,
Note that in the regions to the right and left of the shock it satisfies the
equation trivially.
18
5
Shock fitting
R
For an equation (2.3) , we find q = q(φ) by q = c(φ) dφ. In general, to
determine the shock position xs (t) might be a non-trivial task. Given an initial condition φ(x, 0) = f (x), we can determine the time and location (tb , xb )
where the shock first appears. For t > tb , the analytic solution becomes

multivalued, so we need to replace multivalued part by a discontinuity; this

procedure is known as shock fitting.

Let us consider two approaches, analytic and geometric.

5.1

Analytic method

Suppose that we have found a formula for the multivalued solution φ(x, t)

for t > tb . Typically, φ becomes triple-valued for x between some positions

x1 and x2 (both depend on t). So, φ will have three branches:

for x < x2 (t)
φ1 (x, t)
φ(x, t) = φ2 (x, t)
for x > x1 (t)

φ3 (x, t)

for x1 (t) ≤ x ≤ x2 (t)

with φ3 lying between φ1 and φ2 .

19

The shock location x = xs (t) must be chosen somewhere between x1 and

x2 , and the values of φ immediately to the left and right of the shock are

given by φ1 (xs , t) and φ2 (xs , t). Therefore, we can write the shock relation

(4.2) as follows:

q(φ1 (xs , t)) − q(φ2 (xs , t))

dxs

.

=

dt

φ1 (xs , t) − φ2 (xs , t)

This gives us a first order ODE for xs = xs (t) which should be solved subject

to the initial condition xs (tb ) = xb . The resulting shock solution is then given

for t > tb by

(

φ1 (x, t)

for x < xs (t)
φ(x, t) =
φ2 (x, t)
for x > xs (t)

5.2

Geometric method

Consider equation (3.3) , for which we have a conservation of mass formula

(3.1). Assuming ρ, q → 0 as x → ±∞, the total mass

Z ∞

ρ dx

−∞

20

will, therefore, be constant. This equals the area under the graph of ρ. For

t > tb , ρ becomes multivalued, as a function of x, but the area bounded by

it must remain the same. Recall that when we introduced a shock, we did

this in such a way that (3.1) remained valid. Therefore, the area under the

graph of the shock solution must be the same as the area bounded by the

multivalued solution. As a result, the vertical line x = xs must cut away

lobes of equal area. This is known as the equal areas rule, and it can be

used to determine the location x = xs of the shock.

5.3

Example

Consider the initial value problem

φt + φφx = 0 ,

for x < 0
1
φ(x, 0) = 1 − x for 0 < x < 1
0
for x > 1

We have c(φ) = φ and q =

R

c(φ) dφ = 21 φ2 .

Either from the sketch or by using formulas (2.9), (2.10) , we find that the

shock appears at tb = 1 and xb = 1. From the sketch of the characteristics we

can also see that the shock position x = xs (t) must be within the overlapping

region 1 < x < t and that the two characteristics meeting at the shock carry
21
the values φl = 1 and φr = 0. Using the shock relation, we find that
s=
φ2l − φ2r
1
q(φl ) − q(φr )
=
= .
φl − φr
2(φl − φr )
2
So we obtain that
dxs
1
=
dt
2
with xs (1) = 1 ,
which leads to
xs (t) =
1+t
.
2
As a result, the solution for t > 1 is

(

1 if x < (1 + t)/2
φ(x, t) =
0 if x > (1 + t)/2

22

6

Other boundary conditions

The method of characteristics can also be used in a situation when the boundary condition for equation (2.3) is given along some curve in the x, t plane,

rather than just along the axis t = 0. A standard example is the so-called

signalling problem, when the value of φ is given along the semi-axes:

φ(x, 0) = f (x) for x > 0 ,

and

φ(0, t) = g(t) for t > 0

Here f (x) and g(t) are some given functions.

To solve (2.3) with such boundary data, we construct characteristics

through every point on the boundary. For the points (ξ, 0) on the x-axis,

the characteristics are, as before, the lines

Cξ : x = ξ + c(f (ξ))t ,

(ξ > 0)

with φ = f (ξ) along Cξ .

For a point (0, τ ) on the t-axis, the characteristic that passes through is

given by

Cτ : x = c(g(τ ))(t − τ ) ,

(τ > 0)

and φ = g(τ ) along Cτ .

The rest is the same as before: to determine the solution at point (x, t),

we must determine the characteristic that passes through (x, t) and take the

corresponding value of φ.

Let us illustrate this on one example related to traffic interrupted by

traffic lights.

6.1

Example: Interrupted traffic

Let us suppose that we have a steady flow of cars at a rate qs through a set

of lights and at t = 0 the lights turn red. The cars that have passed through

the lights carry on at the same speed, but those behind must stop. To find

the solution ρ(x, t), we must treat this as two separate signalling problems:

one for the cars ahead of the lights, another – for the cars behind.

23

1. Ahead of the lights: x > 0.

The boundary conditions are ρ = ρs (therefore, q = qs ) along the x-axis

and ρ = 0 (thus, q = 0) along the t-axis since there are no cars passing

through the lights. As a result, the characteristics that originate at the xaxis will have slope cs = c(ρs ), while the characteristics that originate at the

t-axis will have slope c0 = c(0). Note that c = dq/dρ is maximal at ρ = 0

(empty road), thus c0 > cs . It is clear from the sketch that these two families

of characteristics overlap, so there will be a shock formed at t = 0 and x = 0.

This shock will separate the cars which were ahead of the lights from a region

where there are no cars. We will, therefore, have 2 characteristics meeting

at the shock: one with ρ = 0 and another with ρ = ρs . Therefore, the values

on both sides of the shock are:

ρl = 0 , ql = 0 ,

ρr = ρs , qr = qs .

The speed of the shock is, from (4.2),

s1 =

q(ρs )

q(ρs ) − q(0)

=

= v(ρs )

ρs − 0

ρs

i.e. the speed of the last car as we would expect.

ρ

ρs

no cars

x

lights

2. Behind the lights: x < 0.
The boundary conditions here are ρ = ρs (therefore, q = qs ) along the
x-axis and also ρ = ρmax (thus, q = 0) along the t-axis since at x = 0 the
cars are stationary for t > 0. As a result, the characteristics that originate at

the x-axis will have slope cs = c(ρs ), while the characteristics that originate

24

at the t-axis will have slope cmax := c(ρmax ). Note that cmax is the minimum

value of the gradient of q, therefore, cmax < cs . From the sketch, one sees that
these two families of parallel lines overlap, so there will be a shock formed
at t = 0 and x = 0. This shock separates the cars in the queue at the
lights from those still approaching with density ρs . We will, therefore, have
2 characteristics meeting at the shock: one with ρ = ρs and another with
ρ = ρmax . Therefore, the values on both sides of the shock are:
ρl = ρs , ql = qs ,
ρr = ρmax , qr = 0 .
ρ
ρ max
ρs
queue
incoming cars
no cars
lights
x
The speed of this shock is
s3 =
q(ρmax ) − q(ρs )
−q(ρs )
=
0. If ν is small, then our previous assumption that q = Q(ρ) is a
good approximation provided ∂ρ/∂x does not become too large. Note that
no matter how small ν is, the diffusive term will always become important if
the solution requires a shock.
Substituting (7.2) into our basic equation, (3.3), gives
∂ 2ρ
∂ρ ∂Q
+
= ν 2.
∂t
∂x
∂x
The diffusive term causes sharp changes in Q to be smoothed and hence
prevents discontinuities from appearing.
7.1
Travelling wave solution
We can illustrate the effect of the diffusive term by studying solutions of
special type, the so-called travelling wave solutions. First, suppose that
28
ν = 0 in (7.1), and suppose we have a shock solution with φ = φl for x < xs
and φ = φr for x > xs . The shock relations, (4.2), tell us the shock speed.

If φl and φr are constant, then so is s so that the shock moves with uniform

speed. If f (x) is the step function

(

φl for x < 0
f (x) =
,
(7.3)
φr for x > 0

then the shock solution is φ(x, t) = f (x − st), i.e. this is a shock wave

travelling with constant speed s.

In the case ν > 0 we can look for a similar travelling wave solution

φ(x, t) = f (x − st) ,

with a function f is to be determined. As we will see, the presence of the

diffusive term ensures that f is not a discontinuity, but a smooth transition

between φr and φl . To determine f , we think of it as a function of a new

independent variable ξ = x − st. Then

df

∂φ

= −s ,

∂t

dξ

∂φ

df

= ,

∂x

dξ

∂ 2φ

d2 f

=

.

∂x2

dξ 2

Putting this into (7.1) gives

−s

df

dq

d2 f

+

=ν 2 ,

dξ dξ

dξ

q = q(f ) .

This can be integrated to give

−sf + q(f ) = ν

df

+ C,

dξ

(7.4)

where C is a constant of integration. Assume the asymptotic behaviour

ξ → −∞ (far behind) f → φl ,

ξ → +∞ (far ahead) f → φr ,

df /dξ → 0,

df /dξ → 0.

Hence we have two expressions for C:

C = −sφl + ql = −sφr + qr ,

where ql = q(φl ), qr = q(φr ) .

29

(7.5)

It follows that

s=

ql − qr

,

φl − φr

C=

φl qr − φr ql

.

φl − φr

(7.6)

The first relation is precisely the shock relations, (4.2). We could regard this

as another way of deriving the shock relations.

The equation (7.5) is then a first order ODE for f (ξ):

ν

df

= q(f ) − sf − C ,

dξ

(7.7)

with s, C given by (7.6). This can be solved since it is a separable equation:

Z

Z

df

dξ

ξ

=

= .

(7.8)

q(f ) − sf − C

ν

ν

Note that the integration in LHS should be taken for values of f between

φl and φr , to ensure correct asymptotic behaviour. For quadratic q the

integration can be done explicitly.

7.2

Example: Burgers’ equation

The simplest choice is q(φ) = 12 φ2 , in which case the equation (7.1) becomes

φt + φφx = νφxx ,

which is a model used in gas dynamics and called Burgers’ equation.

In this case (7.6) gives

1

1

s = (φl + φr ) C = − φl φr .

2

2

Equation (7.7) for the travelling wave solution then takes the form

ν

df

1

1

= (f − φl )(f − φr ) = − (φl − f )(f − φr ) .

dξ

2

2

(7.9)

df

Note that the derivative dξ

is only zero on the lines f = φl and f = φr .

To achieve both boundary conditions at ξ → ±∞, f must be in the strip

df

φl > f (ξ) > φr , in which case dξ

< 0.
30
To solve (7.9), we separate and use partial fractions to get
Z
df
ξ
1
φl − f
=
=
log
,
2ν
(f − φl )(f − φr )
(φl − φr )
f − φr
(7.10)
which holds for φr < f < φl . We have set the integration constant to zero
as it only introduces an overall shift in ξ. We can express f in terms of ξ by
writing
φl − f
(φl − φr )
= exp
ξ ,
f − φr
2ν
from which
1
1
f = (φl + φr ) − (φl − φr ) tanh
2
2
(φl − φr )
ξ
4ν
.
(7.11)
This describes a smooth transition from φl to φr . For small ν > 0, the

ν

ν

transition mainly takes place over a small interval − φl −φ

< ξ < φl −φ
. In
r
r
the limit ν → 0 we recover the step function (7.3).
8
Hyperbolic Systems
So far we have considered only one dependent variable, φ. We now consider
systems of n equations with n dependent variables u = (u1 , u2 , · · · un ),
∂u
∂u
+A
= 0,
∂t
∂x
(8.1)
where A is an n × n matrix. If A is constant or depends only upon x and t,
then the system is linear. If it also depends on u, but not upon the derivatives
of u, then the system is quasi-linear.
Definition. A system (8.1) is called hyperbolic if A has real eigenvalues and
a full basis of eigenvectors. For instance, it is always true if A is symmetric.
8.1
Linear case: A = const
Let us solve the system (8.1) in the linear case when A is a constant matrix.
In general, the matrix A is not symmetric, so it has right and left eigenvectors
given by
31
li A = λi li ,
Ari = λi ri ,
i = 1, 2, · · · n.
(8.2)
Here li are row vectors and ri are column vectors. We now show that the the
left and right eigenvectors are biorthogonal i.e.
li · rj = 0 for i 6= j .
(8.3)
The product here is the matrix product (of a row by a column), so li · rj is
a 1 × 1 matrix, i.e. a scalar (this coincides with the standard dot product of
two vectors).
To prove (8.3) we multiply Ari = λi ri by lj and lj A = λj lj by ri to get,
lj Ari = λi lj · ri ,
lj Ari = λj lj · ri .
Subtracting gives (λi − λj )lj · ri = 0, which proves (8.3) if the eigenvalues are
all different. Even if they are not, it is possible to find a set of biorthogonal
eigenvectors as long as the ri are linearly independent.
Since the ri are linearly independent, we can write
u=
n
X
vi (x, t)ri .
(8.4)
i=1
This is just a transformation of the dependent variables, u(x, t) = (u1 , . . . , un )
to a new set, v = (v1 , . . . , vn ). Since the eigenvectors are biorthogonal, we
can easily write the new variables in terms of the old ones by multiplying
(8.4) by lj
lj · u =
n
X
vi (x, t)lj · ri = vj lj · rj .
i=1
We therefore have
vi =
li · u
,
li · ri
i = 1 . . . n.
(8.5)
Now substitute for u from (8.4) in (8.1) to get
n
n
n
n
∂ X
∂ X
∂ X
∂ X
vi (x, t)ri +A
vi (x, t)ri =
vi (x, t)ri +
λi vi (x, t)ri = 0 .
∂t i=1
∂x i=1
∂t i=1
∂x i=1
32
Since λi are constant, this can be written as
n
X
∂vi
i=1
∂vi
+ λi
∂t
∂x
ri = 0 ,
and because ri are independent, we get
∂vi
∂vi
+ λi
= 0,
∂t
∂x
i = 1 . . . n.
(8.6)
These are a set of uncoupled equations for the vi , each of which is a linear
advection equation of the form (2.2). As we saw in section 2.1, the solutions
to these are
vi (x, t) = vi (x − λi t, 0)
where vi (x, 0) is the initial data. We can obtain these by substituting the
initial data u = u(x, 0) = u0 (x) in (8.5). Once we have determined v(x, t),
we can get u(x, t) from (8.4). The result is
u(x, t) =
n
X
li · u0 (x − λi t)
i=1
li · ri
ri .
(8.7)
The solution clearly consists of n waves travelling with speeds λi (i = 1 . . . n).
The λi are therefore the wave speeds of the system.
For a nonlinear system A is not constant and we cannot split the system
into uncoupled equations. However, we can always linearise the system locally about the state at some x, t. Small disturbances that vary on a scale
much smaller than the scale of variation of u therefore obey a set of linear
wave equations with wave speeds λi .
33
8.2
Example: a second order equation
Consider the second order linear equation
φtt − φxt − 2φxx = 0 ,
subject to the following initial conditions φ(x, 0) = 0 and φt (x, 0) = f (x),
where f (x) is a given function.
This can be solved by reducing it to a linear system using the variables
u = φt , v = φx . Indeed, we have ut = φtt = φxt + 2φxx from the equation,
hence
ut − ux − 2vx = 0 .
Also, vt = φxt = φtx = ux , therefore,
vt − ux = 0 .
These two equations can be organised into the form (8.1) with
−1 −2
u
φt
.
u=
=
,
A=
−1 0
v
φx
The eigenvalues of A are λ1 = 1, λ2 = −2.
For λ = 1 we have
−2 −2
1
⇒ l1 = (1, −2) , r1 =
.
A − λI =
−1 −1
−1
For λ = −2 we have
1 −2
A − λI =
−1 2
⇒
l2 = (1, 1) , r2 =
2
1
.
Note that li · ri = 3 in both cases.
The initial data for u is
f (x)
u0 (x) =
.
0
From the formula (8.7), the solution is
u(x, t) =
l1 · u0 (x − λ1 t)
l2 · u0 (x − λ2 t)
1
1
r1 +
r2 = f (x−t)r1 + f (x+2t)r2 .
l1 · r1
l2 · r2
3
3
34
Substituting for r1 and r2 gives
2
1
1
1
u = f (x − t) + f (x + 2t) , v = − f (x − t) + f (x + 2t) .
3
3
3
3
Finally, φ can be found by integrating v w.r.t. x, so
1
1
φ(x, t) = − F (x − t) + F (x + 2t) , F 0 = f .
3
3
By the fundamental theorem of calculus, this is the same as
Z
1 x+2t
φ(x, t) =
f (ξ) dξ .
3 x−t
9
Hyperbolic Systems: non-linear case
Let us consider a hyperbolic system
ut + A ux = 0 ,
(9.1)
where A is an n × n matrix. We will assume that A depends on u so the
system is nonlinear. In this case, the eigenvalues and eigenvectors are no
longer constant but also functions of u.
We consider the differential du and use the equation:
du = ut dt + ux dx = −Aux dt + ux dx ,
For any li , we multiply this equation from the left and use the eigenvector
property:
li · du = −(li A)ux dt + li · ux dx = (dx − λi dt)li · ux .
This means that along the characteristics (associated with λi ), we have
dx
= λi (u) ⇒ li · du = 0 .
(9.2)
dt
Notice that the condition li · du = 0 does not mean that du = 0. Since u
is no longer a constant, the first equation in (9.2) has a non-constant right
hand side, so defines a curve, not a straight line.
The second equation in (9.2) defines a nonlinear, first order partial differential equation for u, called the characteristic (or compatibility) condition.
Here i = 1, . . . , n, so we have n such conditions.
35
9.1
Example: Shallow Water Wave Equations
We now consider a simple model to describe a flow of shallow water in a
channel, with h(x, t) being the displacement of the surface above the average
water level and v(x, t) being the downstream velocity,
ht
v h
hx
0
+
=
,
(9.3)
vt
g v
vx
0
where the constant g is acceleration due to gravity.
v h
are
It is easy to see that the eigenvalues of A =
g v
p
λ1 = v − c , λ2 = v + c , where c = gh .
The corresponding left and right eigenvectors are
−c
l1 = (−c, h) , l2 = (c, h) , r1 =
,
g
r2 =
c
g
.
(9.4)
For i = 1, (9.2) takes the form
p
dx
√ dh
= v−c .
−c dh+h dv = 0 ⇒ − g √ +dv = 0 ⇒ d(v−2 gh) = 0 when
dt
h
For i = 2, (9.2) takes the form
c dh + h dv = 0
⇒
d(c + 2
p
dx
= v + c.
gh) = 0 when
dt
We thus have solved the PDEs (9.2) by using integrating factors to reduce
them to exact form, with the result
p
w1 = v − 2 gh = const when
p
w2 = v + 2 gh = const when
dx
= λ1 = v − c ,
dt
dx
= λ2 = v + c .
dt
The functions wi (x, t) are called Riemann invariants.
Since
∂wi
∂wi
∂wi dx ∂wi
dwi =
dt +
dx =
+
dt ,
∂t
∂x
∂t
dt ∂x
36
(9.5)
(9.6)
vanishes whenever
dx
dt
= λi , we must have
∂wi
∂wi
+ λi
= 0.
∂t
∂x
(9.7)
We can invert the formulas (9.5)–(9.6) to write
1
v = (w1 + w2 ) ,
2
c=
p
1
gh = (w2 − w1 ) ,
4
hence
λ1 (w) =
3w1 + w2
,
4
λ2 (w) =
w1 + 3w2
.
4
We can thus rewrite the shallow water wave equations (9.3), purely in terms
of Riemann invariants:
w1t
λ1 (w)
0
w1x
0
+
=
,
(9.8)
w2t
0
λ2 (w)
w2x
0
which has diagonalised the equations. However, because of the form of λi (w),
these equations are not scalar. They form a diagonal, but coupled system.
9.2
Riemann invariants in general
Given a system (9.1) with eigenvalues and eigenvectors λi , li (i = 1, . . . , n),
we say that a function of dependent variables wi = wi (u) is a Riemann
invariant if, for any solution u, we have
dwi = 0 when
dx
= λi ,
dt
(9.9)
or, equivalently,
wi = const when
dx
= λi .
dt
(9.10)
For n = 2 Riemann invariants w1 , w2 always exist. To find w1 , for instance, we write
l1 · du = l11 du1 + l12 du2 = 0
37
and consider this as an ODE by writing it as
l11
du2
=− .
du1
l12
Solving this ODE, we find that its general solution satisfies
(9.11)
w1 (u1 , u2 ) = const
for some function w1 , which will be an i = 1 Riemann invariant. Similarly
for i = 2 and w2 .
E.g., for shallow water equations we have for i = 1:
r
p
dv
g
=
,
l1 · du = −cdh + hdv = 0 ⇒ − ghdh + hdv = 0 ⇒
dh
h
which is a separable equation giving
Z
Z r
p
g
dv =
dh ⇒ v − 2 gh = const = w1 .
h
Similarly for w2 .
In some cases the equation li · du = 0 can be solved by a suitable change
of dependent variables. To illustrate the idea, consider the above example
when
p
l1 · du = −cdh + hdv = 0 with c = gh .
We can use c and v as new dependent variables. Then
c2 = gh
⇒
d(c2 ) = 2cdc = gdh
⇒
dh = (2c/g)dc .
Thus, the equation, written in terms of c, v, becomes
−(2c2 /g)dc + (c2 /g)dv = (c2 /g)(dv − 2dc) = 0
⇒
dv − 2dc = 0 ,
leading to the same conclusion that w1 = v − 2c.
For n > 2 Riemann invariants do not always exist, and if they exist then

the hyperbolic system can be considerably simplified. Namely, similarly to

(9.7) and (9.8), we obtain the system of equations

∂wi

∂wi

+ λi

=0

(i = 1, . . . , n)

∂t

∂x

where wi = wi (u(x, t)). Therefore, a change of dependent variables

(u1 , . . . , un )

−→

(w1 , . . . , wn )

puts the system (9.1) into ‘diagonal form’ (9.12).

38

(9.12)

9.3

Simple wave solutions

Consider n = 2 case, and suppose that we found Riemann invariants w1 , w2

and hence transformed the system (9.1) into a diagonal form

w1t

λ1 (w)

0

w1x

0

+

=

.

(9.13)

w2t

0

λ2 (w)

w2x

0

We cannot solve this in general, but we can construct a large class of solutions

called simple waves.

Namely, set w1 = const for all x, t. Then the first equation in (9.13) is

trivially satisfied. It remains to solve the second equation,

w2t + λ2 (w)w2x = 0 .

Since w1 is constant, this is a scalar hyperbolic equation (2.1) for w2 , so it

can be solved by the method of characteristics. In particular, w2 is constant

along the characteristics (associated with λ2 ). Since w1 is constant everywhere, we conclude that λ2 = λ2 (w1 , w2 ) will be constant, which means that

these characteristics are straight lines. (The characteristics for λ1 will not be

straight lines, in general.)

The solution obtained this way is called an i = 2 simple wave (indicating

that among w1,2 only wi will vary). A general i = 2 simple wave solution

depends on the choice of the constant value of w1 and the choice of initial

data w2 (x, 0) = f (x).

Similarly, we can construct i = 1 simple wave solutions, by setting w2 =

const and solving a scalar hyperbolic equation for w1 . For i = 1 simple waves,

the characteristics for λ1 will be straight lines.

For example, consider the equations (9.3) when w2 = const. Then the

second equation is trivially satisfied, while for w1 we get a scalar hyperbolic

equation

3

1

∂w1

∂w1

+

w1 + w2

= 0,

(9.14)

∂t

4

4

∂x

which can be solved with arbitrary initial data to give w1 = w1 (x, t). Using

this, we can express h and v in terms of w1 , w2 and obtain an i = 1 wave

solution

h=

(w1 − w2 )2

,

16g

v=

w1 + w2

.

2

39

(9.15)

Notice that this method is applicable to (9.8) when the initial data

h0 (x) = h(x, 0) and v0 = v(x, 0)

satisfies the condition

p

w2 = v0 + 2 gh0 is uniform, i.e. constant in x .

More generally, whenever all but one of the Riemann invariants w1 , . . . , wn

for the system (9.1) are uniform on the initial data u0 , the solution to such

initial value problem reduces to a scalar hyperbolic equation.

9.4

Example: The Dam Break Problem

Suppose that for the shallow water equations (9.3) we have the following

initial data at t = 0:

(

(hl , vl ) for x < 0
(h, v) =
(hl , vl ) for x < 0
Find the solution h(x, t), v(x, t), assuming that
p
vl + 2cl = vr + 2cr (c = gh) with vl < vr .
40
(9.16)
Solution: At t = 0 we have
(
wl = vl + 2cl
w2 = v + 2c =
wr = vr + 2cr
for x < 0
for x < 0
Then (9.16) tells us that wl = wr , i.e. w2 = const at t = 0. Therefore, we
can set w2 = const for all x, t and construct the solution in the form of an
i = 1 simple wave.
We know that for i = 1 simple waves, characteristics for λ1 = v − c will
be straight lines, and both w1 , w2 , hence h, v, remain constant along each of
these lines. Note that the i = 2 characteristics are not straight lines because
they cut across i = 1 characteristics with different values of w1 , so that
u 6= const along the i = 2 characteristics.
The characteristic diagram is:
t
i=1 characteristics
i=2 characteristics
x
The i = 2 characteristics are represented by dotted lines. They play no
role in the solution so are not important.
In drawing the above diagram, we used the fact that the slope of the i = 1
characteristics that start at x < 0 equals
λ l = v l − cl ,
while those that start at x > 0 have slope

λr = vr − cr .

41

From (9.16) we get that

vl −vr = 2(cr −cl ) < 0
⇒
cr −cl < 0
⇒
λl −λr = vl −cl −vr +cr < 0 ,
so the i = 1 characteristics do not overlap and there will be no shock in the
solution. Between the lines x = λl t and x = λr t we must have a rarefaction
fan of characteristics.
As a result, the solution is
(hl , vl ) x < λl t ,
(hf , vf ) λl t ≤ x < λr t ,
(h, v) =
(hr , vr ) x > λr t ,

where (hf , vf ) denotes the solution in the rarefaction fan. To determine this

we note that, inside the fan, all the i = 1 characteristics come from the origin.

The slope of the characteristic through (x, t) is therefore x/t, so we have

λ1 = v − c =

x

.

t

We also have

w2 = v + 2c = const = vl + 2cl (= vr + 2cr ) .

These two equations can be

1

c(x, t) =

vl + 2cl −

3

solved to give c and v inside the fan as follows:

x

1

x

, v(x, t) =

vl + 2cl + 2

.

t

3

t

Finally, h(x, t) can be obtained from c(x, t)

h(x, t) = c2 /g =

x 2

1

vl + 2cl −

.

9g

t

42

10

Conservation laws and shock relations

The notion of a local conservation law is unchanged in the context of the

systems. We just require that there exist scalar functions ρ and q of u1 , . . . , un

satisfying

ρ t + qx = 0 ,

where ρt =

n

X

∂ρ ∂ui

,

∂u

∂t

i

i=1

qx =

n

X

∂q ∂ui

.

∂u

∂x

i

i=1

A given system (9.1) may admit more than one local conservation law. Corresponding to the local conservation law with density ρ and flux q, we have

the integral conservation law:

d

dt

Z

x2

ρ(x, t) dx = q(x1 , t) − q(x2 , t) ,

(10.1)

x1

where ρ, q are evaluated on an arbitrary solution u of (9.1). Again, in physical

models, the integral conservation laws are usually more fundamental, with

partial differential equations being a consequence.

10.1

Example: shallow water wave equations

The shallow water wave equations (9.3) are themselves in conservation form:

ht + (vh)x = 0

1

vt + (gh + v 2 )x = 0

2

⇒

ρ = h , q = vh ,

⇒

1

ρ = v , q = gh + v 2 .

2

The first of these represents the conservation of mass.

With a little manipulation, it is possible to find more. For example,

(hv)t = ht v + hvt = −(vhx + vx h)v − h(ghx + vvx )

1 2

2

= −(v hx + 2vvx h + ghhx ) = − hv + gh

,

2

x

2

giving

ρ = hv ,

1

q = hv 2 + gh2 ,

2

(10.2)

43

which represents the conservation of momentum.

The conservation of energy corresponds to

1

ρ = (gh2 + hv 2 ) ,

2

1

q = gh2 v + hv 3 .

2

There are further conservation laws for the shallow water wave equations,

but only three of them have a direct physical meaning.

Shock Relations: We can again use the calculations of Chapter 2 to derive

the differential equation describing the path of a shock wave. Each conservation law, with density and flux (ρ, q) will give an equation for the speed of

the shock:

ql − qr

dxs

=s=

.

dt

ρl − ρr

(10.3)

The “correct” choice of conservation laws depends upon physical arguments,

specific for each case. This will not be discussed here.

10.2

Bore Formation

We now describe a situation with a steady flow of water (height h0 , velocity

v0 ) in a channel (coming from the left) and meeting a wall. At the wall, the

water must stop moving, whilst the steady flow keeps arriving. The result is

that near the wall the height rises to hr and the step between the two levels

hl = h0 and hr , moves at some speed s to the left. The profile at t > 0 is

We use the shallow water wave equations in conservation form:

ht + (vh)x = 0 ,

1

(hv)t + (hv 2 + gh2 )x = 0 ,

2

44

(10.4)

with the channel in the region x < 0 and the wall at x = 0.
The speed of the bore is s < 0, which must satisfy two conditions (10.3),
given by the mass conservation (the first of (10.4)) and conservation of momentum (10.2), which give
1
1
(h0 v0 −hr vr )s = (h0 v02 + gh20 )−(hr vr2 + gh2r ) .
2
2
(h0 −hr )s = v0 h0 −vr hr ,
With vr = 0, these give
(h0 − hr )s = v0 h0 ,
1
g(h20 − h2r ) + h0 v02 = h0 v0 s .
2
(10.5)
In these equations, we take h0 and v0 as given, and use the equations to find
both s and hr . The first equation of (10.5) gives
v0
hr = h0 1 −
s
⇒
1 2 1 h20 (s − v0 )2
gh − g
+ h0 v02 = h0 v0 s ,
2 0 2
s2
giving a cubic equation for s:
1
s3 − v0 s2 − gh0 s + gh0 v0 = 0 .
2
It is easy to see that when s = 0, the cubic is positive, while at s = v0 it
is negative. Since the cubic tends to ±∞ as s → ±∞, it clearly has two
positive roots and one negative root. Taking this negative root, the formula
for hr gives hr > h0 , as expected.

We have thus determined hr and s, as required, giving

(

(h0 , v0 )

for x < st ,
(h, v) =
(hr , 0)
for st < x < 0 .
11
Applications to Gas Dynamics
We now consider one-dimensional gas flow, such as the flow along a pipe.
Since gas is a compressible fluid (think of a bicycle pump!) its density ρ is
variable. Other variables are the pressure p and the velocity (along the tube)
v. An important parameter is γ = cp /cv , where cp is the specific heat at fixed
pressure and cv is the specific heat at fixed volume. For a monatomic gas,
45
γ = 5/3 and for a diatomic gas, γ = 7/5 (this is a good approximation for
air).
Gas dynamics therefore gives us an example of a 3-dimensional system,
for the variables p = (ρ, v, p)T , with an additional parameter γ. Before
stating the equations, it is convenient to define the energy e:
e=
1
p
+ ρv 2 .
γ−1 2
The first term is the thermal energy and the second the kinetic energy per
unit volume. We first state three conservation laws
ρt + (ρv)x = 0 ,
(ρv)t + (p + ρv 2 )x = 0 ,
et + (v(e + p))x = 0 , (11.1)
representing, respectively, the conservation of mass, momentum and energy.
These, in turn, imply the system of equations
0
v ρ
0
ρx
ρt
vt + 0 v 1/ρ vx = 0 .
(11.2)
0
0 γp v
px
pt
The eigenvalue equation is
v−λ
ρ
0
0
v − λ 1/ρ
0
γp v − λ
γp
2
= (v − λ) (v − λ) −
= 0,
ρ
giving
r
λ1 = v − c , λ2 = v , λ3 = v + c ,
where c =
γp
(the sound speed) .
ρ
The corresponding eigenvectors are
l1 = (0, −cρ, 1), l2 = (−c2 , 0, 1), l3 = (0, cρ, 1),
ρ
1
ρ
r1 = −c , r2 = 0 , r3 = c .
c2 ρ
0
c2 ρ
Do Riemann invariants exist in general? It has already been mentioned
that for 3 or more components there is no guarantee that Riemann invariants
exist. We require
li · dp = 0 along the characteristic
46
dx
= λi .
dt
For i = 1 we have
−cρdv + dp = 0
⇒
√
− γpρdv + dp = 0 ,
(11.3)
which can never be transformed (by use of integrating factor) into an exact
differential, since there is no dρ term, and yet the coefficients depend upon
ρ. Thus, the Riemann invariant w1 does not exist. A similar calculation for
i = 3 gives the non-existence of the Riemann invaraint w3 .
However, for i = 2 we have
−c2 dρ + dp = 0
⇒
γp
dp
= c2 =
.
dρ
ρ
This equation is separable and can readily be integrated to give
w2 =
11.1
p
= const on characteristics
ργ
dx
= λ2 = v .
dt
(11.4)
Isentropic Flow
Equation (11.4) just says that along each λ2 characteristic , the Riemann
invariant w2 remains constant. In general w2 takes different constant values
along different characteristics. By definition, a gas flow is called isentropic if
S=
p
= const
ργ
everywhere in the gas, not just along the λ2 characteristics. (The term“isentropic”
is a combination of the Greek word “iso” (same) and entropy. Isentropic flows
occur when the change in flow variables isp
small and gradual.) Without loss
γ
of generality, we may take p = ρ , so c = γργ−1 .
In fact, we can rewrite (11.2) as
ρt
v
ρ
0
ρx
0
vt + γργ−2 S v ργ−1 vx = 0 ,
(11.5)
St
0
0
v
Sx
0
so it is consistent to set S = 1 and reduce the 3-component system to the
following 2-component system:
ρt
v
ρ
ρx
0
+
=
.
(11.6)
γ−2
vt
γρ
v
vx
0
47
This system has eigenvalues and lef eigenvectors as follows:
p
λ± = v ± c , l± = (±c, ρ) , where c = γργ−1 .
(11.7)
This c is just the isentropic form of the previous sound speed and λ± are just
the old λ3 and λ1 , respectively.
Along λ+ characteristics, we have
2
c = 0,
cdρ + ρdv = 0 ⇒ d v +
γ−1
whilst along λ− characteristics, we have
2
−cdρ + ρdv = 0 ⇒ d v −
c = 0,
γ−1
giving the two Riemann invariants
w+ = v +
11.2
2
c,
γ−1
w− = v −
2
c.
γ−1
Example: Piston Driven Expansion
Consider an infinite tube with a piston placed at x = 0 and stationary gas
with constant density and pressure (i.e. ρ = ρ0 , v = 0, p = p0 ) in the region
x > 0. At t = 0 the piston starts to move to the left, with constant speed

vp < 0, which is assumed to be slow enough to make the process isentropic.
The gas will expand to fill the region x > vp t to the right of the piston. The

gas particles at x = vp t will be travelling with the speed vp , same as that of

the piston.

Since the process is isentropic, we can use the 2-component system (11.6)

to find ρ, v and then calculate p = ργ .

Along the characteristics with eigenvalues λ+ = v + c, we have that the

2

Riemann invariant w+ = v + γ−1

c takes a constant value. Along the characteristics with eigenvalues λ− = v − c, we have that the Riemann invariant

2

w− = v − γ−1

c takes a constant value.

Through every point x, t we have two characteristics C± . They have slopes

λ− = v − c < v < λ+ = v + c .
The piston path in the x, t plane is a straight line x = vp t. The fluid velocity
at the piston is v = vp , i.e. it is the same as that of the piston. Thus, if a
48
point (x, t) lies on the piston path x = vp t then C− and C+ look outwards and
inwards, respectively. The conclusion is that on the piston path we have the
characteristics for λ+ = v+c entering the region x > vp t, while characteristics

for λ= v − c can only end at the piston path so they have to start at x ≥ 0

semi-axis.

In particular, this means that the region x > vp t is completely covered

by the characteristics C− that start at the semi-axis x ≥ 0. Since the initial

state is uniform, w− is constant at t = 0 and, therefore, it remains constant

in the whole of the region x > vp t, i.e. we have a simple wave solution. This

means that the solution is constant along the C+ characteristics which are

straight lines.

At the piston we have

v = vp , vp −

2

γ−1

2

cp = w− = const = −

c0 ⇒ cp =

vp + c0

γ−1

γ−1

2

where c = cp at the piston.

We can now sketch the (x, t)-diagram for the C+ characteristics. Some

of them start

qat x > 0, t = 0; these have slope λ+ = v + c = c0 > 0,

where c0 = γρ0γ−1 . Other start at the piston x = vp t; those have slope

v + c = vp + cp . Using the expression for cp , we see that vp + cp < c0 .
Thus, we must also have a fan of C+ characteristics coming out of the origin.
Therefore, the characteristic diagram of C+ characteristics is
We see from the diagram that the disturbance travels into the tube at speed
49
c0 . Within the undisturbed region,
ρ = ρ0 , v = 0 for x > c0 t.

At the piston we have

v = vp = const , c = cp =

γ−1

vp + c0 = const ,

2

thus, these constant values are carried along the C+ characteristics. This

means that the flow is uniform with v = vp , c = cp between the piston,

x = vp t and x = (vp + cp )t.

The C+ characteristics in the region (vp + cp )t ≤ x ≤ c0 t all come from

the origin and therefore have slope x/t. Thus, inside the fan we have

λ+ = v + c =

x

,

t

w− = v −

2

2

c = const = −

c0 .

γ−1

γ−1

From this

i

x

1 h

(γ − 1) + 2c0 ,

c=

γ+1

t

2 x

v=

− c0 .

γ+1 t

The complete solution is therefore

v = vp , c = cp =

γ−1

vp

2

+ c0

for vp t < x < (vp + cp )t,
i
1 h
x
2 x
− c0 , c =
(γ − 1) + 2c0 for (vp + cp )t ≤ x ≤ c0 t,
v=
γ+1 t
γ+1
t
v = 0, c = c0
for x > c0 t.

Note that both the density ρ(x, t) and the pressure p(x, t) are determined

from c(x, t):

c2 = γργ−1 ,

p = ργ =

c2

γ

γ

γ−1

.

The solution looks like

50

v

p

v=0

v=vp

p=p0

p=pp

x

x

51

11.3

Shock relations for gas dynamics

Recall that for a conservation law ρt + qx = 0, the shock relations are given

by (10.3):

ql − qr

dxs

=s=

.

dt

ρl − ρr

For gas dynamics (compressible flow) we have three conservation laws (11.1):

ρt + (ρv)x = 0 ,

(ρv)t + (p + ρv 2 )x = 0 ,

et + (v(e + p))x = 0 ,

where

p

1

e=

+ ρv 2

γ−1 2

⇒

v(e + p) = v

γp

1 2

+ ρv .

γ−1 2

The shock relations are therefore

vl

γpl

γ−1

+

ρl vl − ρr vr

2

2

pl + ρlvl − pr − ρr vr

γpr

1

ρ v 2 − vr γ−1

+ 12 ρr vr2

2 l l

= (ρl − ρr )s ,

= (ρ

l vl − ρr vr )s ,

pl −pr

1

2

2

=

+

(ρ

v

−

ρ

v

)

s.

l

r

r

l

γ−1

2

The first relation can be rewritten as

ρl (vl − s) = ρr (vr − s)

⇒

ρl ul = ρr ur ,

where ul = vl − s and ur = vr − s represent velocities relative to the shock.

With a bit of algebra, one can rewrite all three shock relations in the frame

of shock to get

ρl ul = ρr ur ,

pl + ρl u2l = pr + ρr u2r ,

1 2

γpr

1

γpl

2

ul

+ ρ l u = ur

+ ρr ur .

γ−1 2 l

γ−1 2

(11.8a)

(11.8b)

(11.8c)

The first thing to notice is that a possible solution is

ul = ur = 0, pl = pr , ρl , ρr arbitrary,

which means that s = vl = vr . This type of shock in compressible flow is

called a contact discontinuity, and the associated wave speed is the same on

both sides.

52

It is useful to obtain expressions for the state on one side of the shock in

terms of that on the other. First define

τ=

ul

ρr

= .

ρl

ur

Assuming τ 6= 1 (discontinuity), one then finds from (11.8a)–(11.8c) that

τ=

ρr

ul

γ−1

2

=

=

+

,

ρl

ur

γ + 1 (γ + 1)Mr2

where Mr = ur /cr .

(11.9)

(See Q5 on Examples 3 for the derivation of (11.9).) This relates the density

and velocity jumps across the shock. The quantity Mr is the so-called Mach

number of the shock.

To get the pressure jump, we have from (11.8b)

ul

2

2

2

, by using (11.8a) .

pl = pr + ρr ur − ρl ul = pr + ρr ur 1 −

ur

Therefore

pl

γ−1

2

2γ

2

= 1 + γMr 1 −

−

(M 2 − 1), (11.10)

=

1

+

pr

γ + 1 (γ + 1)Mr2

γ+1 r

which gives the pressure jump across the shock. This shows that if pl > pr

then Mr2 > 1, and vice versa.

By symmetry, there are similar relations with “lef” and “right” swapped.

For example,

2γ

pr

=1+

(M 2 − 1) .

pl

γ+1 l

(11.11)

Hence, Mr2 > 1 if and only if Ml2 < 1. In other words, the flow must be
supersonic on one side of the shock and subsonic on the other.
We can also relate the Mach numbers on both sides of the shock by
multiplying (11.10) and (11.11):
2γ
2γ
2
2
1+
1+
(M − 1)
(M − 1) = 1 .
γ+1 r
γ+1 l
53
11.4
Example: Reflection of a shock wave
Consider a shock wave travelling in the region x > 0 towards a vertical wall

positioned at x = 0. Let us assume that the state behind the shock front is

uniform with constant ρ = ρs , v = vs < 0 and p = ps . The state between the
shock wave and the wall is undisturbed with ρ = ρ0 , v = 0 and p = p0 . We
want to find the pressure at the wall immediately after the shock is reflected.
This is the approximate situation when a shock wave caused by an explosion
hits a solid structure.
Before the reflection, the speed of the shock is, say, s1 < 0. From (11.9)
we have
ul
γ−1
2 c2r
=
+
.
ur
γ + 1 γ + 1 u2r
(11.12)
Since
ul = 0 − s1 , ur = vs − s1
⇒
ul = ur − vs ,
the substitution into (11.12) gives
ur − vs
γ−1
2 c2r
=
+
.
ur
γ + 1 γ + 1 u2r
54
This can be rewritten as a quadratic equation for ur :
2 2
γ−1 2
ur + (vs − ur )ur +
c =0
γ+1
γ+1 r
⇒
(ur )2 −
γ+1
vs ur − c2r = 0 .
2
(Note that ur > vs since s1 < 0.)
After the reflection, the shock wave travels backwards with some speed
s2 > 0. The state between the wall and the shock is unknown, but we still

have v = 0 between the shock and the wall. That is, ρ = ρ1 , v = 0 and

p = p1 . The state on the other side of the shock remains the same, i.e.

ρ = ρs , v = vs and p = ps . Let us put tilde to distinguish the quantities after

the reflection from those before, i.e. we have

uel = 0 − s2 , uer = vs − s2

⇒

uel = uer − vs .

Thus, we obtain the same quadratic equation for uer :

γ+1

vs uer − c2r = 0 ,

2

p

where we used that cer = cr = γps /ρs . Note that now uer < vs since s2 > 0.

The conclusion is that z1,2 = ur , uer are two distinct roots of the equation

(uer )2 −

z2 −

γ+1

vs z − c2r = 0 .

2

The product of the roots is therefore

z1 z2 = ur uer = −c2r

⇒

fr = ur uer = −1 .

Mr M

cr cr

(11.13)

Next, from relation (11.10) we obtain that

pl

γ−1

2γ

+

=

Mr2 .

pr γ + 1

γ+1

Applying this before and after the reflection gives

p0 γ − 1

2γ

+

=

M2

ps γ + 1

γ+1 r

and

p1 γ − 1

2γ f 2

+

=

Mr .

ps γ + 1

γ+1

Multiplying these and using (11.13), we get

p0 γ − 1

p1 γ − 1

4γ 2

+

+

=

.

ps γ + 1

ps γ + 1

(γ + 1)2

55

(11.14)

For a detonation wave we expect the pressure behind the shock wave to be

much greater than that in front, i.e. the term pp0s will be small and can be

neglected. As a result, (11.14) simplifies to

p1

3γ − 1

=

.

ps

γ−1

E.g. for air γ ≈ 1.4 which gives p1 ≈ 8ps . Thus, for a solid structure to

protect from a blast, it should be able to withstand the pressure 8 times

greater than that in the blast!

56

12

Dispersive Waves

In previous chapters the emphasis has been on hyperbolic wave equations,

and mainly on nonlinear equations.

In the present chapter we introduce the concept of dispersion, which is

a property of linear equations which support “sinusoidal solutions” (plane

wave solutions). Dispersion is the phenomenon when waves of different wave

lengths travel at different speeds. Such equations can be solved by Fourier

methods, and dispersion means that a localised wave packet at t = 0 can lose

its coherent structure.

12.1

Dispersion Relations

A linear dispersive equation is one which admits solutions of the form of a

plane wave

φ = A cos(kx − ωt),

(12.1)

where

Wavelength λ =

2π

, a)

k

Period T =

2π

, b).

ω

(12.2)

It is convenient to write this as

φ = A cos(kx − ωt) = Re(Aei(kx−ωt) ),

(12.3)

where Re means the real part. Since we are discussing linear equations, we

can remove the amplitude A. For convenience, we mainly consider solutions

in complex form

φ = ei(kx−ωt) ,

(12.4)

since exponentials are easier to manipulate. The parameter k is called the

wave number and ω the frequency.

In order to satisfy the equation, ω must be a function of k, i.e. we must

have the dispersion relation

ω = ω(k) .

(12.5)

56

In general there will a number of such relations and these are referred to as

different “modes”. Since

φx = ikφ,

φt = −iωφ,

substitution into a linear PDE (with constant coefficients) converts this differential equation into an algebraic relation between ω and k, from which the

dispersion relation (12.5) is determined.

Examples:

1. The linear wave equation. The equation (in one spatial dimension) is

φtt − c20 φxx = 0,

so substitution of (12.4) leads to

(ω 2 − c20 k 2 )φ = 0

ω 2 − c20 k 2 = 0

⇒

⇒

ω = ±c0 k .

We therefore have two modes, φ = eik(x−c0 t) and φ = eik(x+c0 t) . These modes

travel to right and left with speed c0 , which is independent of k, so all wavelengths travel at the same speed.

2. The Klein-Gordon equation. This equation appears in relativistic quantum theory and is a slight modification of the wave equation:

φtt − c20 φxx + m2 φ = 0.

Substitution of (12.4) leads to

(ω 2 − c20 k 2 − m2 )φ = 0

ω 2 − c20 k 2 − m2 = 0,

⇒

p

so ω = ± c20 k 2 + m2 . Again, we have two modes

φ+ = e

ik(x−cp t)

ik(x+cp t)

and φ− = e

,

p

c20 k 2 + m2

where cp =

,

k

which travel to right and left, but now with speed cp , which is a function

of k. This means that different wave lengths travel at different speeds, so a

wave packet (built out of many Fourier components) will change its shape in

time. This is called dispersion.

57

3. The heat equation: φt − µφxx = 0. Substitution of (12.4) leads to

(−iω + µk 2 )φ = 0

⇒

ω + iµk 2 = 0,

so we just have the one mode

φ = ei(kx+iµk

2 t)

2

= e−µk t eikx ,

which represents a damped (spatial) oscillation. Thus, the heat equation is

dissipative rather than dispersive.

For a pane wave solution (12.4), the quantity

θ = kx − ωt,

(12.6)

is called the “phase”. It measures the position in the cycle between a “crest”

where φ is a maximum (cos θ = 1) and a “trough” where φ is a minimum

(cos θ = −1). Note that

∂θ

= −ω ,

∂t

∂θ

= k.

∂x

(12.7)

Note also that particular values of θ move with velocity

cp =

ω

,

k

(12.8)

called the phase velocity. Generally, cp is a function of k, so plane waves with

different wavenumbers move with different speeds.

Note that the dispersion relation is given by a polynomial whose degree in

ω is that of the highest time derivative in the equation. It follows that there

are as many solutions (modes) as the order of the highest time derivative.

Also note that, since we require a real dispersion relation, the factors of i

must cancel and hence only equations which contain only even or only odd

derivatives have dispersive wave solutions. For example, this is not true for

the heat equation, as it contains first order time derivative and second order

x-derivative and so is not dispersive.

58

12.2

Fourier Decomposition

In general a disturbance will consist of many different wavenumbers. Since we

are considering linear equations, we can find the solution for each wavenumber and then add them together to get the complete solution (principle of

superposition). If there is only one mode, then solution can be written as a

Fourier integral

∞

Z

Φ(k)eikx−iω(k)t dk.

φ(x, t) ==

(12.9)

−∞

Since there is only one mode, the equation must be first order in time and

an appropriate initial condition is

Z

∞

φ(x, 0) = φ0 (x) =

Φ(k)eikx dk.

(12.10)

−∞

A standard result from the theory of Fourier transforms tells that this relation

holds true if Φ(k) is the Fourier transform of the left-hand side, namely,

∞

Z

1

Φ(k) =

2π

φ0 (x)e−ikx dx.

(12.11)

−∞

One can view (12.10) is a Fourier decomposition of the initial data into a

superposition of sinusoidal waves with varying wavenumbers. Given Φ(k),

the solution φ(x, t) is given by (12.9).

If there are two modes, ω = ω1 (k), ω = ω2 (k), then (12.9) becomes

Z

∞

[Φ1 (k)eikx−iω1 (k)t + Φ2 (k)eikx−iω2 (k)t ] dk.

φ(x, t) =

(12.12)

−∞

Since the equation must be second order in time, appropriate initial conditions are now

φ = φ0 (x),

∂φ

= ψ0 (x) at t = 0.

∂t

Imposing the first condition gives

59

1

Φ1 (k) + Φ2 (k) =

2π

Z

∞

φ0 (x)e−ikx dx,

(12.13a)

−∞

and the second condition gives

i

ω1 (k)Φ1 (k) + ω2 (k)Φ2 (k) =

2π

Z

∞

ψ0 (x)e−ikx dx.

(12.13b)

−∞

These two equations determine Φ1 and Φ2 . The extension to n modes is

obvious.

12.3

Group Velocity

Although the Fourier transform described in the previous section gives us

the solution for arbitrary initial data, it is generally not possible to obtain

an explicit expression for the solution unless both the initial data and the

dispersion relation are very simple. However, we are usually not interested in

the complete solution from t = 0 to t = ∞, we merely want to have some idea

of how the waves behave at large distances from the location of the original

disturbance.

Let us start with a simple example. Consider the linear superposition of

two plane waves with nearby values of (k, ω):

cos((k + ∆k)x − (ω + ∆ω)t) + cos((k − ∆k)x − (ω − ∆ω)t)

= 2 cos(∆k x − ∆ω t) cos(kx − ωt),

with ∆k k, which represents a plane wave with wavelength λ = 2π/k, with

a modulated amplitude, where the modulation has much longer wavelength

λ̄ = 2π/(∆k), as shown in the figure:

2

1

20

40

60

-1

-2

60

80

100

This can be heard as beats when 2 tuning forks with nearby frequencies

are played together.

We can see that, whereas the original (carrier) wave has phase velocity

∆ω

dω

ω

, the modulation has velocity

≈

as ∆k ≈ 0. This leads us to define

k

∆k

dk

the group velocity of plane waves as

cg =

dω

.

dk

(12.14)

p

Example: For the Klein-Gordon equation, we have ω = c20 k 2 + m2 , so

p

c20 k 2 + m2

ω

dω

c2 k

cp = =

,

cg =

=p 2 0

.

k

k

dk

c0 k 2 + m2

We see that

cp

m2

= 1 + 2 2 > 1,

cg

c0 k

so cp > cg . If we travelled with the envelope as in the figure above (at speed

cg ) we would notice wave crests overtaking us.

The definition (12.14) of the group velocity can also be justified from the

Fourier representation (12.9). Indeed, consider the real part of the integrand

(as a function of k)

Re Φ(k)ei(kx−ωt) = Φ(k) cos(kx − ωt),

(12.15)

which, for large t, will rapidly oscillate. As an illustration, here is the graph

1 2

of (12.15) for Φ(k) = e− 2 k and ω = k − k 3 . :

1.0

0.5

k

0.2

0.4

0.6

0.8

1.0

1.2

1.4

-0.5

-1.0

Integration over the rapidly oscillating part of this graph gives approximately zero. The main contribution is coming from the central part, which

is around the stationary point of the phase

θ = kx − ωt

⇒

dθ

dω

=x−

t = 0 when

dk

dk

61

x

= cg ,

t

(12.16)

so is determined by the group velocity. In the case of the above graph, and

for large time asymptotics, we have

dθ

= x − (1 − 3k 2 )t

dk

⇒

1

k ≈ √ = 0.577,

3

when

x

1,

t

corresponding to the less oscillating part of the graph.

This illustrates the general principle (known as the method of stationary

phase) which says that the main contribution of Fourier integrals, such as

(12.9), corresponds to the stationary points of the phase, given by (12.16). It

implies that for large (fixed) x and t, the plane waves that make a significant

contribution to φ(x, t) are those for which x/t = cg (k).

12.4

General Wave Packets

One of the disadvantages of studying plane wave solutions is that only constant coefficient equations have such solutions. If we wish to include nonhomogeneities of the medium, then we must allow some coefficients to depend

upon the spatial variables. Furthermore, nonlinear equations usually don’t

have such solutions.

Instead of an exact plane wave solution, we may consider a solution of

the form

φ(x, t) = A(x, t)eiθ(x,t) ,

(12.17)

with slowly varying amplitude A(x, t). Such functions can be solutions of the

above types of equation and can also be used in perturbation theory, which

is beyond the scope of this course. Instead, our goal here is just to give a

little more insight into phase and group velocities.

We may consider the phase function θ(x, t) to be more general than the

usual θ = kx − ωt. We generalise k and ω to be functions of (x, t), defined

by

k(x, t) =

∂θ

,

∂x

ω(x, t) = −

∂θ

,

∂t

(12.18)

considered as local wave number and frequency. We can follow curves x(t),

along which either θ(x, t) or k(x, t) are constant. For example

dθ

∂θ dx ∂θ

dx

=

+

= −ω + k

= 0 when

dt

∂t

dt ∂x

dt

62

dx

ω

= = cp .

dt

k

Since cp (k(x, t)) is a function of (x, t), this defines a curve.

Note that the crests/troughs of the wave (12.17) correspond to θ = πn

with integer n. Therefore, to follow crests/troughs, we must travel at the

phase velocity.

On the other hand, the consistency condition (θxt = θtx ) of definitions

(12.18) is

kt + ωx = 0,

which is in the form of a local conservation law. Furthermore, the dispersion

relation, giving ω(k), leads to

∂k

∂k

+ cg (k)

= 0 since

∂t

∂x

dω

= cg .

dk

This implies that

dk

∂k dx ∂k

=

+

= 0 when

dt

∂t

dt ∂x

dx

= cg (k).

dt

Since k is constant along these characteristics, cg (k) is also constant, so we

have

dx

k = const on straight lines for which

= cg (k).

dt

Hence, if we wish to follow waves of a given wave length (fixed k) then we

must travel at the group velocity for that particular wave length.

This shows the difference between the phase velocity and the group velocity.

If we travel at the phase velocity, we follow crests/troughs but the wavelengths around us change (i.e. distance to nearby crests/troughs change).

On the other hand, if we travel at the group velocity, the wavelength and

frequency of the nearby waves remain constant, but crests/troughs will keep

passing us.

We can use the above formulas to obtain the functions k(x, t), ω(x, t) and

θ(x, t) for a given dispersion relation.

If we imagine a disturbance created at time t = 0, close to the origin

x = 0, then along the line x/t = cg , we have a constant wave number k. For

a given dispersion relation, we know the function cg (k), so can rearrange, to

get the function k(x, t). The dispersion relation then gives ω(x, t) and then

we solve the system (12.18) to obtain θ(x, t).

63

Example: Consider the beam equation

φtt + γ 2 φxxxx = 0

⇒

ω = ±γk 2 ,

giving

cg = 2γk =

x

t

⇒

k(x, t) =

x

2γt

⇒

ω(x, t) =

x2

4γt2

⇒

θ(x, t) =

x2

+δ.

4γt

We therefore have that k(x, t) is constant on straight lines and θ(x, t) is

constant on parabolas.

64

13

Water Waves

Consider a body of water with undisturbed depth h and suppose that the

surface is disturbed so that the depth is h + φ(x, t), where φ is small. Then

the theory of water waves tells us that there are solutions of the form

φ = A cos(kx − ωt),

(13.1)

provided ω satisfies the dispersion relation for water waves,

σk 2

.

ω = k tanh(kh) g +

ρ

2

(13.2)

(see King & Billingham, Chapter 4). Here g is the acceleration due to gravity,

ρ is the water density and σ is the surface tension. We have g = 981 cm

s−2 , ρ = 1 g cm−3 , σ = 74 g s−2 . Note that this system has two modes with

opposite speeds.

As explained below, for certain ranges of values of k and h, the relation

(13.2) can be simplified. It will be convenient to introduce the wavenumber

r

km =

13.1

ρg

σ

r

⇒

λm = 2π

σ

= 1.73 cm.

ρg

(13.3)

Gravity Waves

If k km , λ λm , then

2 !

σk 2

σk 2

k

g+

=g 1+

=g 1+

≈g

ρ

ρg

km

and so the surface tension term can be neglected. The dispersion relation is

then

ω 2 = gk tanh(kh).

(13.4)

These are called gravity waves. The phase and group velocities are

65

r

cp (k) =

1

2kh

cg = cp (k) 1 +

.

2

sinh(2kh)

g tanh(kh)

,

k

(13.5)

Consider the limiting case kh = 2πh/λ 1, i.e. wavelengths that are

short compared to the depth. In this limit, tanh(kh) → 1 and so

p

ω = gk,

r

cp =

g

,

k

1

cg =

2

r

g

,

k

1

i.e. cg = cp .

2

(13.6a)

Note that, since we are considering gravity waves, this requires λm λ h.

Now consider the case kh = 2πh/λ → 0, i.e. wavelengths that are very

long compared to the depth h. Then tanh(kh) is well approximated by kh

and so

ω=

p

ghk,

cp =

p

gh,

cg =

p

gh ,

i.e. cg = cp ,

(13.6b)

which is exactly what we had for the shallow water equations, Sec. 9.1, with

v = 0.

13.2

Capillary Waves

If k km or λ λm , then

σk 2

=

ρg

k

km

2

σk 2

g 1+

ρg

1

⇒

≈g

σk 2

σk 2

=

ρg

ρ

so that surface tension now dominates over gravity. The dispersion relation

reduces in this limit to

ω2 =

σk 3

tanh(kh).

ρ

(13.7)

Note that this does not depend on g. These waves are called capillary waves.

The phase and group velocities are

66

r

cp =

σ

k tanh(kh) ,

ρ

3

2kh

cg = cp 1 +

.

2

3 sinh(2kh)

(13.8)

In the short wavelength limit, kh → ∞, we have

s

cp =

σk

,

ρ

3

cg = cp ,

2

(13.9)

while in the long wavelength limit, kh → 0, we have

s

cp =

13.3

σh

k,

ρ

cg = 2cp .

(13.10)

Combined Gravity and Capillary Waves

Unless we are dealing with waves in a puddle, we usually have kh = 2πh/λ

1 when λ ‘ λm = 1.73 cm. In this limit we have

σk 3

,

ω 2 = gk +

ρ

s

cp =

g σk

+

,

k

ρ

1

cg = cp

2

ρg + 3σk 2

ρg + σk 2

. (13.11)

The phase velocity has a minimum

cp = cg = cm = 23.2 cm s

We have

λ > λm

λ < λm
cg < cp
cg > cp

−1

r

at k = km =

2π

ρg

, λm =

= 1.73 cm.

σ

km

(13.12)

(gravity wave branch),

(capillary wave branch).

The minimum group velocity is

cg = 17.9 cm s−1 at λ = 2.54λm = 4.39 cm.

For the derivation of all these facts, see Q4 on Examples 4.

This is illustrated by the next plot.

67

(13.13)

70

60

50

cg

40

30

cp

20

0.5

14

14.1

1.0

1.5

2.0

2.5

k/k m

Application to Water Patterns

Point Source

Consider an object dropped into deep water. The disturbance introduces

waves with a range of wavenumbers which propagate out from the point

of impact. Since wavenumbers travel at the group velocity, at time t, a

wavenumber k0 will be found at

x = cg (k0 )t .

For a stone of size l ‘ λm = 1.73 cm, thrown into a pond, the combined

gravity and capillary limit from Sec. 13.3 is appropriate. Since the group

velocity has a minimum at cg = 17.9 cm s−1 , there is a quiescent (calm)

circular region of radius x = 17.9t cm after t seconds.

The amplitude of the different wavenumbers is determined by the initial

disturbance. The main waves (with the largest amplitude) have k ‘ π/l,

where l is the size of the stone. The main rings will therefore be found at

68

x = cg (π/l)t.

For a larger stone of size l λm in a deep pond, the limit (13.6a) is more

appropriate. Hence,

r

1

1 g

= cp

cg =

2 k

2

which is a decreasing function of k. The longest (and more prominent) waves

are therefore to be found on the outside of the expanding pattern.

Remark: Since cg < cp in the latter case, an observer travelling with the
wave crests will see the waves suddenly disappear. This is because the crests
travel at the phase velocity, but the energy (amplitude) travels at the group
velocity.
14.2
Obstacle
Consider an obstacle in a steady stream of water flowing with speed V past
a stationary obstacle that spans the stream e.g. a log tied to the banks, or a
step in the stream bed:
Surface
Surface
V
or
Bed
Bed
In some cases the obstacle generates a steady wave pattern i.e. the position of the crests and troughs does not change with time. For this to happen,
the crests must move with speed V relative to the stream and must therefore
have wavenumbers such that
69
cp (k) = V.
Since the phase velocity has a minimum at cp = 23.2 cm s−1 , this can only
happen if V ≥ 23.2 cm s−1 . There are then two values of k such that
cp (k) = V , one on the gravity branch, k = kg < km and one on the capillary
branch, k = kc > km .

We have

cg (kg ) < cp (kg ) = V on the gravity branch,
cg (kc ) > cp (kc ) = V on the capillary branch.

The longer, gravity waves (with k < km ) therefore cannot move upstream, so
these waves are found on the downstream side of the obstacle. The shorter,
capillary waves (with k = kc ) move faster than the stream, so these are found
on the upstream side.
14.3
Ship Waves
Consider a ship with size L moving in water with depth h L (i.e. deep
water). The ship will produce waves with wavelength ' L, so that we have
kh '
h
1,
L
i.e. (13.6a) applies. This means that
1
cg = cp .
2
We assume that the ship acts as a point source moving with constant
speed V in the positive x direction. We are looking for a wave pattern that
is stationary in the frame of the ship. Consider a wave crest in this pattern
at some point C, and assume that it was emitted by the ship at some angle φ
when the ship was at a point A. We assume that the ship moved to a point
B since. Suppose that the ship moved from A to B in time t i.e. AB = V t.
A wave crest moves in the direction of φ with phase velocity, cp (k). To
be stationary with respect to the ship, such a crest must have
cp (k) = V cos φ .
(13.14)
70
This is illustrated in the picture.
However, wavenumbers travel with the group speed, therefore, we must
have
1
AC = cg (k)t = cp (k)t .
2
Using (13.14), we get
1
1
AC = V t cos φ = AB cos φ .
2
2
This means that the point C must lie on a small semicircle, as shown in the
picture.
D
T
C
φ
A
θ
B
E
71
We conclude that among all waves emitted from A, those that contribute
to the steady pattern (when the ship is at B) must be located on the small
circle. Keeping B fixed and varying A, we get a family of circles that fill a
wedge with semi-angle θ. As a result, the waves that are stationary (relative
to the ship) fill that wedge. To determine the angle θ, we draw a tangent
from B to the small semi-circle in the picture. If E is its centre and T the
tangency point, then its radius is R = AE = ET and AB = 4R. Clearly we
have
sin θ =
R
1
ET
=
= .
EB
3R
3
In particular, we see that the angle does not depend on the ship velocity,
which is a surprising result.
72
University of Leeds
MATH 5373M
Advanced Linear and Nonlinear Waves
The Korteweg-deVries Equation:
Special Solutions and the Bäcklund Transformation
1
Introduction
We have seen in the historical introduction, that the KdV equation is a remarkable equation.
Kruskal and his coworkers intensively studied the KdV equation and published their results
in a series of papers between 1967 and 1974. This was the birth of a new subject (now
referred to as Integrable Systems Theory), which had (and continues to have) enormous
impact on both pure and applied mathematics.
Here we just discuss a few very simple properties.
1.1
The Effect of Dispersion
We already saw in Chapter 2 (of the joint Level 3/5 part of the module) that the equation
ut + uux = 0
can lead to wave breaking and shocks. Also in Chapter 2, we saw that the shock could be
avoided (smoothed out) by adding small dissipation, in which case Burgers’ equation was
shown to have a travelling wave solution with a “tanh” type profile. In place of dissipation,
Zabusky and Kruskal added dispersion, with a small parameter δ:
ut + uux + δ 2 uxxx = 0.
Looking carefully at the dashed plot of Figure 4 of the historical introduction, we see that
at breaking time tB (of the dispersionless equation with δ = 0) we no longer have a vertical
profile, but that a small peak has formed. This develops into the array of 8 solitons shown
at time t = 36tB .
2
If the same equation (with δ = 0.06) is used to evolve the initial data u(x, 0) = 2e−x ,
then we get Figure 1. Again it looks like a number of solitons is emerging.
2.0
3.0
2.0
2.5
1.5
1.5
1.0
Out[53]=
2.0
1.5
1.0
1.0
0.5
0.5
0.5
-4
-2
2
4
-4
2
-2
4
-4
-2
2
4
2
Figure 1: The solution of the KdV equation with initial data u(x, 0) = 2e−x .
1
In the next section we calculate the travelling wave solution of the KdV equation. The
existence of such a solution is not at all remarkable, since many equations have travelling wave solutions. What is remarkable is the existence of multi-soliton solutions, which
are exact solutions which represent a number of travelling wave solutions which interact
elastically like particles (no radiation emitted). The 2−soliton solution is presented below.
The mathematical machinery used for constructing this solution is called a Bäcklund
transformation.
2
Solitary Wave Solutions of the KdV Equation
We start with the KdV equation of the specific form
ut = uxxx + 6uux = (uxx + 3u2 )x ,
(1)
so waves will move to the left.
A travelling wave solution, is a solution of the form:
u(x, t) = q(ξ) ,
ξ = x + ct,
so that ut = cqξ , ux = qξ , uxx = qξξ , etc. If u(x, t) satisfies the equation (1), then the
function q(ξ) satisfies the ordinary differential equation:
cqξ = (qξξ + 3q 2 )ξ
⇒
qξξ + 3q 2 − cq = α,
where α is a constant. This is a Newtonian equation with energy:
E=
1 2
qξ − cq 2 + q 3 − αq,
2
(2)
with potential (when α = 0) being graphed (with c = 1) in Figure 2(a). The phase portrait
V(q)
qξ
0.4
0.02
0.2
0.3
-0.2
0.5
q
0.2
-0.2
0.4
0.6
q
-0.02
-0.2
(a) The potential V (q).
(b) The phase portrait.
Figure 2: The KdV Travelling Wave as a Dynamical System.
(q, qξ ) is obtained by plotting ‘level curves’ E=constant. A series of these is numerically
plotted in Figure 2(b). The equilibrium point (at the minimum of the graph) q = 31 c (with
2
1 3
qξ = 0) has E = − 54
c . Close to this we have small amplitude periodic solutions which
‘look like’ sinusoidal oscillations. As we increase the amplitude (keeping E < 0) the solution
remains periodic, but loses its sinusoidal nature. The curve which circulates this equilibrium
point, but with E = 0, is the separatrix. In this diagram we can see a negative energy
‘particle’ coming from the left and being ‘scattered’ (unbounded, non-periodic solution). A
positive energy ‘particle’ comes from the left, is temporarily captured by the minimum (at
the equilibrium point), and then continues on its unbounded orbit.
If we plot q(ξ) we find a function which is periodic (for q(ξ) > 0, E < 0), but whose period
increases as E → 0 from below (unlike SHM, whose period is independent of amplitude).
As the amplitude approaches 0.5 ( 12 c), the period approaches ∞ and the solution takes the
“bell shaped” form seen by John Scott Russell.
-5
0.5
0.5
0.25
0.25
0
5
(a) Small periodic orbit.
-15
-10
-5
0.5
0.25
0
5
10
15
(b) Large periodic orbit.
-10
0
-5
5
10
(c) Infinite period orbit.
Figure 3: The soliton as a limit of periodic orbits.
Actually, the periodic solutions can be explicitly written down in terms of elliptic functions. For a given energy E, (2) is just a separable, first order ODE for q(ξ), giving
Z
dq
p
= ξ − δ,
2E + cq 2 − 2q 3
where δ is an arbitrary constant, called the phase (indicating the initial point on the periodic
orbit). For general E we can express q(ξ) in terms of the Weierstrass elliptic function.
The separatrix (infinite period solution) corresponds to E = 0, so has a very simple
formula. The above integral now takes the form
Z
dq
√
= ξ − δ,
q c − 2q
We can explicitly evaluate the integral on the left by substituting q =
c sech2 η tanhη dη. We then have:
Z
Z
dq
2
√
dη,
= √
c
q c − 2q
so η =
√
c
2 (ξ
− δ), giving:
q=
c
sech2
2
√
c
(ξ − δ) ,
2
c
2
sech2 η, so dq =
(3)
which is an analytic expression for the separatrix in Figure 2(b).
Again we emphasise that the amplitude (the maximum height of the graph) is directly
proportional to the speed c, so that bigger waves travel at greater speed. The phase δ just
gives the position of the maximum (at ξ = δ).
3
Since the equation is nonlinear, a “linear superposition” of two or more solitary waves
is no longer a solution. However, there do exist exact solutions of the KdV equation, which
represent multi-soliton solutions (meaning, the interaction of many solitary waves).
3
Symmetries and Similarity Solutions
The KdV equation has 4 Lie point symmetries:
x−translation X = x + τ, T = t, U = u,
t−translation X = x, T = t + τ , U = u,
scaling symmetry X = λx, T = λ3 t, U = λ−2 u, where λ = eτ .
Galilean symmetry X = x − ct , T = t , U = u + 61 c, where c = τ .
In terms of the parameter τ , each of these is close to the identity transformation, meaning
that it reduces to the identity transformation when τ → 0.
In each case:
UT = UXXX + 6U UX ⇔ ut = uxxx + 6uux.
For instance, in the case of the scaling symmetry,
UT = λ−5 ut ,
UX = λ−3 ux ,
UXXX = λ−5 uxxx,
leading to:
UT − UXXX − 6U UX = λ−5 (ut − uxxx − 6uux).
This means that, if U (X, T ) is a solution of the KdV equation, then so is u(x, t) =
λ2 U (λx, λ3 t). For instance, if we had started with the travelling wave solution with c = 1:
1
1
2
q = sech
(x + t − δ) ,
2
2
√
the general formula just corresponds to rescaling, with λ = c.
3.1
Similarity Solutions
For each symmetry, we can build invariants and invariant solutions.
Example 1 (The Scaling Symmetry) With respect to the scaling symmetry, we have
invariants
z = x/t1/3 and ϕ = u(x, t)t2/3 ,
since
X
λx
x
= 3 1/3 = 1/3
T 1/3
(λ t)
t
and U T 2/3 = λ−2 u(λ3 t)2/3 = ut2/3 .
Writing u(x, t) = t−2/3 ϕ(z), the KdV equation implies
2
1
ϕzzz + 6ϕϕz + zϕz + ϕ = 0.
3
3
This equation belongs to a very special class of ODE classified by Paul Painlevé and his
collaborators around 1900. Its general solution cannot be written in terms of elementary
functions, but needs the introduction of the second Painlevé transcendent.
4
Example 2 (The Galilean Symmetry) With respect to the Galilean symmetry, we eliminate c to obtain
X
1
x − ct
x
U+
=u+ c+
=u+ .
6T
6
6t
6t
Since t is itself an invariant, we can consider a solution of the KdV equation of the form
u(x, t) = −
x
+ ϕ(t)
6t
⇒
ϕt = −
ϕ
t
⇒
u(x, t) = −
x+α
.
6t
The freedom x → x + α is just the x translation symmetry, so the Galilean symmetry gives
x
.
us a simple rational solution u(x, t) = − 6t
4
The Bäcklund Transformation
Consider the variable w(x, t), defined by u(x, t) = wx . If u(x, t) satisfies the KdV equation,
then w(x, t) satisfies
wt = wxxx + 3wx2 ,
(4)
which is called the potential KdV equation (PKdV).
Let two functions w1 (x, t), w2 (x, t) satisfy
1
w1x + w2x = 2k 2 − (w1 − w2 )2 .
2
(5)
If w1 (x, t) is a solution of (4) then so is w2 (x, t). The proof of this is to check the compatibility
of (4) and (5) by differentiating (5) with respect to t and then using (4) to eliminate w1xt .
If we then use (5) and its derivatives to eliminate all x−derivatives of w1 , then we find that
w2 satisfies (4).
The calculation is rather laborious, but the result is very useful. The Bäcklund transformation (5) enables us to construct a new solution w2 (x, t) out of a known solution w1 (x, t).
The 1−Soliton Solution. The function w1 (x, t) ≡ 0 is clearly a solution of (4). The
Bäcklund transformation (5) then gives
Z
dw2
1
1
= (x − δ),
(6)
w2x + w22 = 2k 2 ⇒
2
2
2
4k − w2
2
which gives one of the following
w2 (x) = s1 (x) = 2k tanh(k(x − δ))
or w2 (x) = s2 (x) = 2k coth(k(x − δ)).
To find the t−dependence, we differentiate (6) with respect to x, to obtain
w2xx + w2 w2x = 0,
2
w2xxx + 3w2x
− 4k 2 w2x = 0
⇒
w2t − 4k 2 w2x = 0,
so
w2 (x, t) = 2k tanh(k(x + 4k 2 t − δ))
For the KdV equation, this gives
u = 2k 2 sech2 (k(x + 4k 2 t − δ))
or w2 (x, t) = 2k coth(k(x + 4k 2 t − δ)).
(7)
or u = −2k 2 cosech2 (k(x + 4k 2 t − δ)),
the first of which is the travelling wave solution (3) (with c = 4k 2 ) previously found. The
second solution is singular when x+4k 2 t−δ = 0, but is needed for constructing the 2−soliton
solution below.
5
4.1
Bianchi Permutability
Suppose we start with a solution w (of the PKdV equation) and use (5) with parameters k1
and k2 to respectively build solutions w1 and w2 . We then use (5) with parameters k2 and
k1 to respectively build new solutions from w1 and w2 and require that these are identical
(named w12 ). This is depicted in Figure 4. Each arrow corresponds to a Bäcklund transw1
w
H
H k
HH2
HH
j w
12
*
k1
*
k1
HH
k2HHH
j
w2
Figure 4: A commutative Bianchi diagram
formation. This commutativity is called “Bianchi permutability” and leads to an algebraic
relation between the functions at the four vertices of the diagram.
To calculate this formula, consider the Bäcklund relations (5) corresponding to each of
the 4 arrows and eliminate the derivatives. The remaining terms are then simplified to give
(w12 − w)(w1 − w2 ) = 4(k12 − k22 ).
(8)
The 2−Soliton Solution. If we choose w = 0, w1 = 2k1 tanh(k1 (x + 4k12 t)), w2 =
2k2 coth(k2 (x + 4k22 t)) then
w12 =
2(k12 − k22 )
.
k1 tanh(k1 (x + 4k12 t)) − k2 coth(k2 (x + 4k22 t))
With k1 = 1, k2 = 2, we have
u(x, t) =
12(3 + 4 cosh(2(x + 4t)) + cosh(4(x + 16t)))
∂w12
.
=
∂x
(3 cosh(x + 28t) + cosh(3(x + 12t)))2
(9)
This represents a pair of solitons (one of height 8 to the right of another of height 2) coming
from the right, overtaking at x = 0, with them emerging (unscathed!) to the left (see Figure
5). The only change is that the fast soliton shifts its position relative to the slow one (called
a phase shift), as can be seen in the density plot of Figure 6.
4.2
The Wronskian Representation
Equation (6) is an example of a Riccati equation, which can be linearised by the transformation
fxx
fx2
w2 = 2(log f )x ⇒ w2x = 2
⇒ fxx − k 2 f = 0.
− 2
f
f
In particular, the solutions (7) correspond to
f1 (k) = cosh(k(x + 4k 2 t − δ)) and f2 (k) = sinh(k(x + 4k 2 t − δ)).
6
(10)
10
u
10
8
8
6
6
6
4
4
4
2
0
-5
Out[72]=
10
5
10
x
-10
10
8
8
6
6
4
4
-5
u
2
0
-5
u
2
-10
10
8
2
-10
u
5
10
x
5
10
x
-10
-5
0
5
10
x
u
2
0
5
10
x
-10
0
-5
Figure 5: The 2−soliton solution of the KdV equation
0.4
0.2
0.0
-0.2
-0.4
-10
0
-5
5
10
Figure 6: Density plot, showing phase shift
The 1−soliton solution of the KdV equation is then written
u(x, t) = 2k 2 sech2 (k(x + 4k 2 t − δ)) = 2(log(τ ))xx ,
(11)
where τ = f1 (k).
Remarkably, the 2−soliton solution (9) has a compact form of the same type:
u(x, t) = 2(log(τ ))xx ,
(12)
where τ = 12 (3 cosh(x + 28t) + cosh(3(x + 12t)). Furthermore, this τ −function has a remarkably simple construction as the Wronskian of the two functions f1 and f2 of (10):
τ=
1
(3 cosh(x + 28t) + cosh(3(x + 12t)) =
2
f1 (1) f2 (2)
.
f1x (1) f2x (2)
Of course, the overall constant multiple is unnecessary, since we are going to differentiate
log(τ ).
7
The general N −soliton solution is built in a similar way, with an N × N Wronskian:
τN = W [f1 (k1 ), f2 (k2 ), f1 (k3 ), . . .],
(13)
where we alternately choose f1 and f2 , as shown, and where the parameters ki are distinct.
We then have the corresponding N −soliton solution:
uN (x, t) = 2(log(τN ))xx
5
Rational Solutions of the KdV Equation
As well as the multi-soliton solutions, we can also build explicit formulae for rational solutions of the KdV equation:
P (x, t)
u(x, t) =
,
Q(x, t)
where P (x, t) and Q(x, t) are polynomial functions of x and t. For instance, we have already
seen that the similarity solution, corresponding to the Galilean symmetry, is u = −x/(6t).
We now show how to construct a family of rational solutions as singular limits of the soliton
solutions.
We consider the two functions f1 and f2 (of (10)), but slightly changing the definition
of δ:
f1 (k, δ) = cosh(kx + 4k 3 t + δ) and f2 (k, δ) = sinh(kx + 4k 3 t + δ).
(14)
We consider the τ −function (13), with ki = εκi , where ε ≪ 1, and choose some specific
values of δi (which we now allow to be complex). We expand τ as a series in ε, the coefficients
being polynomial in the variables x and t. With appropriate choices of δi , we retain the
lowest order part of the surviving series and use this as a polynomial τ −function, which
then gives a rational solution of the KdV equation.
5.1
Limit of the 1−Soliton Solution
Corresponding to the 1−soliton solution, we have
1
τ1 = cosh(κ1 εx + 4κ31 ε3 t + δ1 ) = cosh δ1 + εκ1 x sinh δ1 + ε2 κ21 x2 cosh δ1 + · · · .
2
If we choose δ1 = iπ/2, so that cosh δ1 = 0, sinh δ1 = i, then
τ1 = εiκ1 x + O(ε3 ) = εiκ1 (x + O(ε2 )).
Dropping the constant multiplier (since we are taking the logarithm and then differentiating), and then taking the limit ε → 0, we have:
τ1r = x
⇒
ur1 = 2(log τ1r )xx = −
which is a rational solution of the KdV equation.
8
2
,
x2
(15)
5.2
Limit of the 2−Soliton Solution
Corresponding to the 2−soliton solution, we have
τ2
= εκ2 cosh(εκ1 x + 4ε3 κ31 t + δ1 ) cosh(εκ2 x + 4ε3 κ32 t + δ2 )
−εκ1 sinh(εκ1 x + 4ε3 κ31 t + δ1 ) sinh(εκ2 x + 4ε3 κ32 t + δ2 )
= (κ2 cosh(δ1 ) cosh(δ2 ) − κ1 sinh(δ1 ) sinh(δ2 )ε + (κ22 − κ21 ) cosh(δ1 ) sinh(δ2 )ε2 x + · · ·
Choosing δ1 = iπ/2, δ2 = 0, the first surviving term is ε4 , so we write:
τ2 =
1
iκ1 κ2 (κ22 − κ21 )ε4 (x3 − 12t) + O(ε6 )).
3
We define
τ2r = x3 − 12t
5.3
⇒
ur2 = 2(log τ2r )xx = −
6x(x3 + 24t)
.
(x3 − 12t)2
(16)
A Direct Calculation of Polynomial τ Functions
The function w(x, t) = 2(log(τ ))x is a solution of the PKdV equation (4), if and only of
τ (x, t) is a solution of
2
τ τxt − τx τt = τ τxxxx − 4τx τxxx + 3τxx
,
(17)
and then u(x, t) = 2(log(τ ))xx is a solution of the KdV equation (1).
The functions τ1r and τ2r are just polynomial solutions of this equation. The next polynomial solution is
τ3r = x6 − 60x3 t − 720t2
⇒
ur3 = 2(log τ3r )xx =
9
12x(43200t3 − 5400t2 x3 − x9 )
.
(x6 − 60tx3 − 720t2 )2
(18)
University of Leeds
MATH3374/5373M
Examples 1
1. Consider the initial value problem
φt + φφx = 0 ,
0 for x < 0
φ(x, 0) = f (x) = x for 0 < x < 1
1 for x > 1

(a) Sketch the characteristics and hence find the solution φ(x, t).

(b) Repeat the same for the step-like initial data

(

0 for x < 0
φ(x, 0) =
1 for x > 0 .

(c) How your answer in part (b) would change if the equation was

replaced by φt + φ2 φx = 0?

2. For the following initial value problems, (1) write down the parametric equation of the envelope of the characteristics, and (2) determine

the breaking time and location:

(a)

φt + φφx = 0 ,

φ(x, 0) = 1 − tanh x ,

(b)

2

φ(x, 0) = e−x .

φt + φ2 φx = 0 ,

3. Consider the initial value problem

for x < −1
1
ut + uux = 0 ,
u(x, 0) = −x for − 1 < x < 0
0
for x > 0

(a) Sketch the characteristics and hence find the solution u(x, t) for

t < 1. Show, using the sketch or otherwise, that for t > 1 the solution u(x, t) becomes triple-valued, and determine all possible values of

u(x, t) for t > 1.

(b) Conclude that a shock forms at t = 1. Use the shock relation to

find the speed, s, of the shock and hence determine the solution u(x, t)

for t > 1.

1

4. Consider the initial value problem for the equation ut + uux = 0

with the following initial data (tent function):

for 0 ≤ x ≤ 12

x

u(x, 0) = 1 − x for 12 ≤ x ≤ 1

0

otherwise .

(a) Sketch the characteristics and hence find the solution u(x, t) for

t < 1. Show, using the sketch or otherwise, that for t > 1 the solution u(x, t) becomes triple-valued, and determine all possible values of

u(x, t) for t > 1.

(b) Conclude that a shock forms at t = 1. Use the shock relation to

find the speed, s, of the shock and hence determine the solution u(x, t)

for t > 1.

5. Consider the equation

∂φ 1 ∂φ2

+

=0

∂t

2 ∂x

with the initial data

x 1.

(a) Show that a shock forms at x = 1, t = 1. Determine the solution

for t < 1. (Hint: For determining φ, use the implicit equation φ =
f (x − c(φ)t).)
(b) Show that, for t > 1, the shock position xs (t) is determined by the

first order ODE

p

dxs

1

= 2 [2xs t − 1 − t + 1 − 4xs t + 2t + t2 ] ,

dt

4t

with the initial condition xs = 1 at t = 1.

6. A simplified model of traffic flow without diffusion is governed by

the equation

∂ρ ∂Q

+

= 0,

∂t

∂x

where Q(ρ) = ρ(1 − ρ). Find the shock speed for a shock with ρl , ρr

on the left and right.

(a) Obtain ρ(x, t) for the initial data

ρr x > 0

,

ρl x < 0
when: (i) c(ρl ) < c(ρr ) (rarefaction); (ii) c(ρl ) > c(ρr ) (shock).

ρ(x, 0) =

(b) Consider a stationary queue of traffic at the traffic lights at t = 0:

0 x < −L
1 −L < x < 0
ρ(x, 0) =
0 x>0

At t = 0, the lights turn green. Find the density of traffic ρ(x, t) for

t > 0.

University of Leeds

MATH3374/5373M

Solutions 1

1. (a) Done in lecture on Jan 25.

(b) See the video on Minerva.

(c) See the video on Minerva.

2. (a) See the video on Minerva.

(b) We have c(φ) = φ2 , so

2

F (ξ) = c(f (ξ)) = e−2ξ .

The breaking time corresponds to the minimum negative value of F 0 (ξ). Since

F 0 (ξ) = −4ξe−2ξ

2

is smooth, its minimum is attained at points with F 00 (ξ) = 0. Computing F 00

gives

2

F 00 (ξ) = e−2ξ (16ξ 2 − 4) ,

which vanishes at ξ = ±1/2. From a sketch of F we see that the minimum

of F 0 is attained at ξb = 1/2.

From that

2

tb = −

√

e2ξb

1

=

=

e/2 ,

F 0 (ξb )

4ξb

and

xb = ξb −

F (ξb )

1

1

= +

= 1.

0

F (ξb )

2 4ξb

The characteristics for this problem are given by

2

Cξ : x = ξ + e−2ξ t .

To determine their envelope, we apply formula (2.9) from the notes and

obtain:

1

,

x=ξ+

4ξ

2

e2ξ

t=

.

4ξ

The envelope has two branches: one for ξ > 0 and another for ξ < 0. We
need ξ > 0 only, since we are interested in the branch that lies in the upper

half-plane t > 0. One can further express ξ in terms of x from the first

formula,

√

x ± x2 − 1

ξ=

,

x ≥ 1,

2

and substitute this to find t = t(x). This represents two curves, each defined

for x ≥ 1. The result is illustrated in the picture.

3. (a) Done in lecture on Jan 26.

(b) From part (a) we know that the solution breaks at tb = 1 at xb = 0.

After that, it becomes triple-valued, with

1

if x < t − 1
x
φ(x, ...
Purchase answer to see full
attachment
Tags:
Exercise Example
Arbitrary function
Linear waves
Nonlinear waves
linear dispersive system
User generated content is uploaded by users for the purposes of learning and should be used following Studypool's honor code & terms of service.

## Reviews, comments, and love from our customers and community:

This page is having a slideshow that uses Javascript. Your browser either doesn't support Javascript or you have it turned off. To see this page as it is meant to appear please use a Javascript enabled browser.