Download Characteristics & Riemann Invariants in First-Order PDE Systems and more Study notes Mathematics in PDF only on Docsity!
B5b Applied Partial Differential Equations 2–
2 First-order systems
2.1 Introduction
In this chapter we consider first-order systems of PDEs of the form
A(x, y, u)
∂u
∂x
+ B(x, y, u)
∂u
∂y
= c(x, y, u), (2.1)
where now the unknown u is an n-dimensional vector function of x and y, c is also an
n-dimensional vector, A and B are n × n matrices. Again, we only consider quasilinear
equations, so that (2.1) is linear in the first derivatives of u. Equations of the form (2.1) arise
frequently in physical problems (e.g. gas dynamics), in which it is often the case that one of
the independent variables represents time, so we will also write (2.1) in the form
A(x, t, u)
∂u
∂t
+ B(x, t, u)
∂u
∂x
= c(x, t, u). (2.2)
The other motivation for studying (2.1) is that many higher-order scalar equations may be
written in this form.
Example 2.1 The shallow-water equations describe the flow of a thin layer of liquid over a flat
surface under gravity. If x measures horizontal distance and t is time, then the dimensionless height
h(x, t) and velocity u(x, t) satisfy
∂h
∂t
+ u
∂h
∂x
+ h
∂u
∂x
∂u
∂t
+ u
∂u
∂x
∂h
∂x
which may be written in the form (2.2), with
u =
h
u
, A =
, B =
u h
1 u
, c =
Example 2.2 Laplace’s equation
∂x
∂y
may be written as a first-order system by setting u = ∂φ/∂x, v = −∂φ/∂y. Then u and v satisfy the
Cauchy–Riemann equations
∂u
∂x
∂v
∂y
∂u
∂y
∂v
∂x
which are in the form of (2.1), with
u =
u
v
, A =
, B =
, c =
(Once u and v are known, φ is found by solving either ∂φ/∂x = u or ∂φ/∂y = −v.)
2–2 OCIAM Mathematical Institute University of Oxford
2.2 Cauchy data
Definition
For the PDE system (2.1), Cauchy data is to specify u on a curve Γ in the (x, y) plane, i.e.
x = x
0
(s), y = y
0
(s), u = u
0
(s), s
1
≤ s ≤ s
2
, (2.8)
and the PDE (2.1) together with the data (2.8) is known as the Cauchy problem. As for the
scalar case, we now ask whether the Cauchy problem determines the first derivatives of u.
Differentiating (2.8) with respect to s, we find
du
0
ds
=
dx
0
ds
∂u
∂x
dy
0
ds
∂u
∂y
. (2.9)
Now we consider (2.1) and (2.9) as a (2n) × (2n) matrix system for the (2n)-dimensional
vector (∂u/∂x, ∂u/∂y)
T
:
(
A B
x
′
0
I y
′
0
I
) (
∂u/∂x
∂u/∂y
)
=
(
c
u
′
0
)
, (2.10)
where
′
is shorthand for d/ds and I is the n × n identity matrix. The system (2.10) may be
solved uniquely for the first derivatives of u provided the determinant of the matrix on the
left-hand side is nonzero, and this determinant may be rearranged to give
det
(
x
′
0
B − y
′
0
A
)
6 = 0. (2.11)
This is the condition on the initial data for the the first derivatives of u to be locally deter-
mined. It clearly reduces to the condition found for scalar equations when n = 1.
Cauchy–Kowalevski theorem
Now we state a generalisation of the theorem previously introduced for scalar PDEs. For
simplicity we suppose that a coordinate transformation is used to shift the boundary Γ onto
the y-axis, where we specify u:
u = u
0
(y) on x = 0, y
1
≤ y ≤ y
2
. (2.12)
Clearly, so long as u 0
is differentiable, we can calculate ∂u/∂y directly:
∂u
∂y
=
du 0
dy
on x = 0, y
1
≤ y ≤ y
2
. (2.13)
We can then use the PDE (2.1) to solve for ∂u/∂x,
∂u
∂x
= A
− 1
c − A
− 1
B
∂u
∂y
= f
(
x, y, u,
∂u
∂y
)
, say, (2.14)
so long as A is invertible. Thus the PDE may be written in the form (2.14) provided |A| 6 = 0,
which is the same as condition (2.11), with x 0
= 0, y
0
= s.
Now suppose that u 0
(y) is analytic at a point y 0
∈ (y 1
, y 2
) and that f is analytic in all
its arguments at the point
(
0 , y
0
, u
0
(y
0
),
du 0
dy
(y
0
)
)
.
2–4 OCIAM Mathematical Institute University of Oxford
The homogeneous system (2.17, 2.18) implies that the jumps in ∂u/∂x and ∂u/∂y must be
zero unless the determinant of the system is zero, which implies that
det
(
dx
dξ
B −
dy
dξ
A
)
= 0, (2.19)
i.e. that C is a characteristic.
Classification
The slopes of the characteristics satisfy the eigenvalue problem
dy
dx
= λ where det (B − λA) = 0. (2.20)
Thus λ satisfies an nth-order polynomial equation, whose roots may be complex in general.
A 2 × 2 system may be classified as follows, depending on the eigenvalues λ.
- If there are two distinct real eigenvalues, then the system is said to be hyperbolic.
- If there is one repeated real eigenvalue, then the system is parabolic.
- If the eigenvalues are complex, then the system is called elliptic.
Example 2.3 Consider the quasilinear second-order PDE
a
∂x
+ 2b
∂x∂y
+ c
∂y
= f, (2.21)
where a, b and f are in general functions of x, y, φ, ∂φ/∂x and ∂φ/∂y. We can write (2.21) as the
first-order system
a b
∂u
∂x
b c
∂u
∂y
f
, where u =
∂φ/∂x
∂φ/∂y
If a, b, c and f are independent of φ, then we can ignore the uncoupled equation ∂φ/∂x = u, and the
characteristic slopes satisfy
b − λa c − λb
⇒ aλ
− 2 bλ + c = 0. (2.23)
The system is thus hyperbolic if b
> ac, parabolic if b
= ac, elliptic if b
< ac.
In dimensions higher than two, there are clearly many possible combinations of real,
complex, distinct and repeated roots of the polynomial equation (2.20), and there is no such
simple classification. However, we still define an equation as hyperbolic if (2.20) has n
distinct real roots λ. Since the matrices A and B depend on x, y and, in general, also on the
solution u, the type of the equation (i.e., hyperbolic, elliptic, parabolic or some hybrid) may
also vary with position.
Now, according to the Cauchy–Kowalevski theorem, provided all our coefficients and initial
data are analytic and the condition (2.11) is satisfied, there is a unique solution for u in a
neighbourhood of Γ. Nevertheless, unless the PDE is hyperbolic, the Cauchy problem is in
general ill posed. This may manifest itself in several ways. For example, the unique local
solution may blow up arbitrarily close to Γ or may be pathalogically sensitive to the initial
data.
B5b Applied Partial Differential Equations 2–
Example 2.4 For the Cauchy–Riemann equations, the characteristic slopes satisfy
= 0 ⇒ λ = ±i, (2.24)
so the system is elliptic.
Suppose we try to solve (2.6) in y > 0 subject to the Cauchy data
u = u
(x) = 0, v = v
(x) =
x
on y = 0, (2.25)
where and δ are positive constants. Notice that the initial data are analytic and bounded: v
(x) is in
the range (0, ] for all x. It may be shown
that the solution of this problem is
u =
xy
[x
+ (y − δ)
] [x
+ (y + δ)
]
, v =
[
x
− y
]
[x
+ (y − δ)
] [x
+ (y + δ)
]
Thus, however small is, both u and v blow up at the point (0, δ): by choosing δ small, we can make
the solution break down arbitrarily close to the initial curve y = 0.
For the remainder of this chapter we restrict our attention to hyperbolic systems, for which
the Cauchy problem is generally well posed, and for which characteristic methods analogous
to those used for scalar equations can be applied. So, at each point in the (x, y) plane, we
assume that (2.15) defines n distinct real eigenvalues λ. Thus, by solving dy/dx = λ for each
of these n characteristic slopes, we can obtain in principle n families of characteristics for an
n-dimensional hyperbolic system.
2.4 Integration along characteristics
Reduction to an ODE
Suppose λ is a real eigenvalue of (2.20); recall that λ is in general a function of x, y and u,
since A and B are. Now the matrix (B − λA) is singular, so there exists a left eigenvector l
T
,
such that
l
T
(B − λA) = 0
T
, that is l
T
B = λl
T
A. (2.27)
Multiplying the PDE (2.1) on the left by l
T
, we obtain
l
T
A
∂u
∂x
T
B
∂u
∂y
= l
T
c
⇒ l
T
A
(
∂u
∂x
∂u
∂y
)
= l
T
c. (2.28)
Along characteristics, whose slope is dy/dx = λ, we have
l
T
A
du
dx
= l
T
c. (2.29)
This is the equivalent of the ODE satisfied by u along characteristics in the scalar case. There
is one ODE of the form (2.29) satisfied along each of the n families of characteristics.
the Cauchy–Riemann equations imply that u + iv is a function of z = x + iy; here, u and v are the real
and imaginary parts of the complex function
f (z) =
δ
i
z
.
B5b Applied Partial Differential Equations 2–
For linear PDEs with c = 0 , as in Example 2.5, we can always find a complete set of n
Riemann invariants. Furthermore, for linear PDEs, the characteristics may be found indepen-
dently of the solution. We thus obtain a system of n algebraic equations for the components of
u in terms of arbitrary functions that are constant along each family of characteristics. This
suggests a plausible method for solving hyperbolic systems numerically. If A, B and c are
approximated as being locally constant near Γ, then the resulting autonomous linear system
has a complete set of Riemann invariants. Thus the solution u a small distance from Γ may
be found by solving the resulting system of algebraic equations. By repeating this process,
the solution may be continued further still from the initial data. This proposed procedure
suggests that the Cauchy problem should usually be well posed for hyperbolic systems.
The following two examples show that, when c is nonzero, even linear PDEs have no
Riemann invariants in general.
Example 2.6 The system
∂u
∂x
∂v
∂y
= u + v,
∂u
∂y
∂v
∂x
may be written as
∂u
∂x
∂u
∂y
u + v
where u = (u, v)
T
. As in Example 2.5, the characteristic directions are dy/dx = λ = ± 1 and the
corresponding left eigenvectors are l
T
= (1, ±1). Now, the ODEs satisfied along the characteristics are
d
dx
(u ∓ v) = u + v on
dy
dx
only one of which is integrable:
d
dx
(u + v) = u + v ⇒ e
−x
(u + v) = const on
dy
dx
Nevertheless, we can use this single Riemann invariant to find the general solution in this case.
Since e
−x
(u + v) is constant when (y + x) is constant, we may write
u + v = e
x
f (y + x), (2.43)
where f is an arbitrary function. Along the other family of characteristics we have (y − x) = const = k,
say, so
d
dx
(u − v) = u + v = e
x
f (y + x) = e
x
f (k + 2x) on y − x = k. (2.44)
This may be integrated with respect to x to give
u − v =
e
(x−y)/ 2
y+x
e
s/ 2
f (s) ds + g(y − x), (2.45)
where g is a second arbitrary function. The general solution is, therefore,
u =
e
x
f (y + x) +
g(y − x) +
e
(x−y)/ 2
y+x
e
s/ 2
f (s) ds, (2.46a)
v =
e
x
f (y + x) −
g(y − x) −
e
(x−y)/ 2
y+x
e
s/ 2
f (s) ds. (2.46b)
2–8 OCIAM Mathematical Institute University of Oxford
Example 2.7 The system
∂u
∂x
∂v
∂y
= u,
∂u
∂y
∂v
∂x
has characteristic equations
d
dx
(u ∓ v) = u on
dy
dx
which cannot be integrated. There are no Riemann invariants and no way to find the general solution
explicitly.
The situation is even worse for fully nonlinear systems, where the characteristics depend
on the solution u. Even when such systems have a complete set of n Riemann invariants,
since we do not know in advance the curves along which each is conserved, we cannot in
general find an explicit solution.
Example 2.8 The shallow-water equations (2.3) have characteristic slopes given by
det (B − λA) =
u − λ h
1 u − λ
= (u − λ)
− h = 0
dx
dt
= λ = u ±
h, (2.49)
and the corresponding left eigenvectors are (1, ±
h). The characteristic ODEs are
dh
dt
h
du
dt
= 0 on
dx
dt
= u ±
h
d
dt
u ± 2
h
= 0 on
dx
dt
= u ±
h, (2.50)
so the Riemann invariants u ± 2
h are preserved along the characteristics dx/dt = u ±
h. Although
the system has two Riemann invariants, we cannot infer the general solution.
Regions of influence
Recall that, for scalar PDEs, where Cauchy data are only given on a finite initial curve Γ, the
solution is only determined in the so-called domain of definition, penetrated by characteristics
emanating from Γ. In n dimensions, the domain of definition is the region penetrated by all
n families of characteristics originating at Γ. Where there are at least one but fewer than n
families of characteristics, the solution is influenced but not fully determined by the Cauchy
data on Γ. The region swept out by all the characteristics intersecting Γ is therefore called
the region of influence.
Example 2.9 We return to the system considered in Example 2.5, which has the general solution
u =
f (y + x) +
g(y − x), v =
f (y + x) −
g(y − x). (2.51)
If we are given the Cauchy data u = u
(x), v = v
(x) on y = 0, 0 < x < 1 , then we can determine the
two arbitrary functions f and g:
f (s) = u
(s) + v
(s), 0 < s < 1 , (2.52a)
g(s) = u
(−s) − v
(−s), − 1 < s < 0. (2.52b)
2–10 OCIAM Mathematical Institute University of Oxford
ppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppp
pppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppp
D
y
x
Γ
γ
Figure 2: Schematic showing the boundary curve Γ, closed by a curve γ to enclose a region
D.
for which ∂u/∂x becomes unbounded in finite time if u
′
0
is ever negative. To continue such
solutions, it is necessary to allow u to be discontinuous across curves in the (x, y) plane, again
referred to as shocks. Since the PDE (2.1) does not make sense on such a curve, we have to
use a weak formulation of the problem. The theory is very similar to the scalar case, so we
omit most of the details.
The first step is to write the system in conservation form
∂P
∂x
∂Q
∂y
= R, (2.56)
where P , Q and R are vector-valued functions of x, y and u.
2
Now, as illustrated in Figure 2,
we form a closed region D by closing Γ with a second curve γ, then multiply (2.56) through
by a test function ψ, assumed to be suitably differentiable and to vanish on γ. Then we
integrate over D and, just as for the scalar case, Green’s theorem leads to the following weak
formulation of (2.1):
∫
Γ
ψ (P dy − Q dx) =
∫ ∫
D
P
∂ψ
∂x
∂ψ
∂y
A function u(x, y) that satisfies (2.57) for all suitable test functions ψ is called a weak
solution of (2.1). If u is continuously differentiable and satisfies (2.57), then it is also a
classical solution of (2.1). However, (2.57) also makes sense if u is discontinuous.
Shocks
Now we look for a weak solution in which u is smooth everywhere except a curve C, across
which it is discontinuous. As shown in Figure 3, C divides the region D into two sub-regions
D 1
and D
2
. The integral on the right-hand side of (2.57) may be split up into two integrals
over D 1
and D
2
respectively. Since u is smooth within D
1
and within D
2
, Green’s theorem
In fact, an arbitrary PDE system cannot always be written in this form, but it is usually possible for
physically-motivated problems that are based on conservation laws.
B5b Applied Partial Differential Equations 2–
ppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppp
pppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppp
ppppppppppppppppppppppppppppppppppppp pppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppp
C
D
1
D
2
C
1
C
2
Figure 3: Schematic showing the shock C dividing D into two regions D 1
and D
2
. The
integration paths on either side of C are denoted C 1
and C
2
.
may then be used, along with the fact that u satisfies the PDE (2.1), giving
∫ ∫
D
P
∂ψ
∂x
∂ψ
∂y
∫ ∫
D
P
∂ψ
∂x
∂ψ
∂y
∫ ∫
D
P
∂ψ
∂x
∂ψ
∂y
=
∮
∂D
ψ (P dy − Q dx) +
∮
∂D
ψ (P dy − Q dx). (2.58)
Then, since ψ is assumed to be zero on γ, (2.57) reduces to
∫
C
ψ([P ]
−
dy − [Q]
−
dx) = 0, (2.59)
where [ ]
−
denotes the jump across the shock. This holds for all test functions ψ, and the
slope of the shock must therefore satisfy the Rankine–Hugoniot condition
[P ]
−
dy
dx
= [Q]
−
. (2.60)
The scalar Rankine–Hugoniot condition is clearly reproduced if n = 1 but, in higher
dimensions, (2.60) gives us n relations between dy/dx and the jumps in the n components of
u. For semilinear equations, we have
P = Au, Q = Bu, R = c +
∂A
∂x
u +
∂B
∂y
u, (2.61)
so the Rankine–Hugoniot condition is
[Au]
−
dy
dx
= [Bu]
−
⇒
(
B −
dy
dx
A
)
[u]
−
= 0. (2.62)
Thus u can only be discontinuous if the determinant of the matrix on the left-hand side is
zero, which implies that dy/dx is equal to a characteristic slope. In other words, shocks occur
on characteristics for semilinear equations. This is not true for general quasilinear systems,
though.
B5b Applied Partial Differential Equations 2–
in which the two equations represent conservation of mass and energy respectively. The corresponding
Rankine–Hugoniot relations,
v =
[uh]
[h]
[
hu
+ uh
]
[
hu
h
]
give shock conditions that are quite different from (2.65). We have a choice of conserving either mass
and momentum or mass and energy. In fact, two different kinds of bores are observed in practice:
turbulent bores that conserve momentum but lose energy, and undular bores that conserve energy but
not momentum.
Causality
Once we have chosen a particular conservation form and, thus, a particular weak formulation,
there is still the possibility of more than one solution existing if we allow shocks. As for scalar
PDEs, there are some shock solutions that, although they satisfy the Rankine–Hugoniot
conditions, are unphysical and should be eliminated. There are several methods for doing
this, of which we concentrate on causality: making sure that information propagates into a
shock, rather than out of it.
An n-dimensional hyperbolic system has n families of characteristics, so a shock intersects
2 n of them: n from either side. If there are k outgoing characteristics, then there are (2n − k)
characteristics going in. We also have the n Rankine–Hugoniot relations, giving a total of
(3n − k) pieces of information on the shock. The unknowns are the n components of u, on
either side of the shock, and the shock slope dy/dx, giving a total of (2n + 1). For the number
of equations to equal the number of unknowns, we require (3n − k) = (2n + 1), that is
k = (n − 1). (2.69)
This is the condition for a shock to be causal: there must be (n − 1) characteristics leaving
the shock (and, therefore, (n + 1) going in).
For scalar equations, n = 1 so there should be no characteristics leaving the shock, 2
going in, as imposed previously. For two-dimensional systems, n = 2 so we need one family
of characteristics going out of a shock and three going in. Schematics of the characteristics
for two alternative shock solutions are shown in Figure 5. In diagram (a), three families of
characteristics enter the shock and one leaves: this solution is causal. In diagram (b), only
one family of characteristics propagates into the shock, so this solution is non-causal and
should be discarded.
Example 2.13 Consider a shock solution of the shallow-water equations satisfying the momentum-
conserving Rankine–Hugoniot conditions (2.65). Recall that the characteristic velocities are given by
dx/dt = u ±
h. Assume the shock is moving in the positive x-direction. Then, for the solution to
be causal, as shown in Figure 5, there should be one set of characteristics entering the shock from
behind and two from the front. In other words, the shock should be moving faster than one of the
characteristic speeds behind the shock and faster than both the characteristic speeds ahead of the shock.
Thus we obtain the following four inequalities
u
h
> v, u
h
< v, u
h
< v, u
h
< v, (2.70)
where, as before, v is the speed of the shock. From these it follows that
(u
− v)
< h
and (u
− v)
> h
2–14 OCIAM Mathematical Institute University of Oxford
-4 -2 0 2 4 -4 -2 0 2 4
0 0
1 1
2 2
3 3
4 4
5 5
t t
(a) (b)
x x
Figure 5: Schematics of the characteristics for two alternative shock solutions; the shock is
shown as a heavy dashed line.
and the Rankine–Hugoniot condition (2.65) then leads to the condition
h
< h
Thus the height behind a bore must be greater than that ahead; otherwise the bore is not causal and
the discontinuity cannot be maintained. That all the inequalities in (2.70) follow from (2.72), may be
verified from the identities
(u
− v)
h
(h
+ h
2 h
, (u
− v)
h
(h
+ h
2 h
which follow from (2.65).
