GWU Linear & Nonlinear Waves Arbitrary Functions & Linear Dispersive System Exercises

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.