THE WAVE EQUATION

1. Vibration of an elastic spring

We consider 1-D compressive vibration of a linear elastic spring. Suppose

that the spring is uniform with mass density ρ and Young’s modulus E.
Suppose there is only a small displacement u(x) in x direction, and dis-
placements in other directions are negligible. Then, the spring section origi-
nally takes the place of [x1 , x2 ] now has a length

s(x2 ) − s(x1 ) = [x2 + u(x2 )] − [x1 + u(x1 )] = (x2 − x1 ) + (u(x2 ) − u(x1 )).

Thus, the relative prolongation of the spring at a point x is given by

ds
(x) − 1 = ux (x).
dx
Applying the Newton’s law to the section [x1 , x2 ], we have
Z x2
ρ utt (x, t) dx = Eux (x2 , t) − Eux (x1 , t).
x1

Since
Z x2
ux (x2 , t) − ux (x1 , t) = uxx (x, t) dx,
x1

we have
Z x2 Z x2
ρ utt (x, t) dx = E uxx (x, t) dx.
x1 x1

This leads to the wave equation

utt = c2 uxx , (1.1)

p
with c = E/ρ > 0. We will see that the constant c is the wave speed, i.e.
the waves characterized by the wave equation (1.1) travel at the speed c.

1
 2. Solutions consist of traveling waves

Suppose that u(x, t) satisfies the wave equation (1.1). We introduce new
variables
ξ = x + c t and η = x − c t,
and define the function
v(ξ, η) = u(x, t). (2.1)
Then, by the chain rule, we have
ux = vξ ξx + vη ηx = vξ + vη
uxx = vξξ + 2vξη + vηη .
Similarly, we have
utt = c2 (vξξ − 2vξη + vηη ).
Thus, as a consequence of the wave equation (1.1), we have

0 = utt − c2 uxx = −4 c2 vξη .

Hence, the function v defined by (2.1) satisfies the following equation
vξη = 0. (2.2)
This implies that the function v(ξ, η) is of the form
v(ξ, η) = f (ξ) + g(η) (2.3)
for some functions f and g, which need to be determined by the initial and
boundary conditions. On the other hand, it is easily verified that for any
smooth functions f and g, the function
u(x, t) = f (x + c t) + g(x − c t) (2.4)
solves the wave equation (1.1). In fact, we have

utt = c2 (f 00 + g 00 ) and uxx = f 00 + g 00 ,

which immediately gives utt = c2 uxx .
Noticing that the function g(x − c t) represents a traveling wave moving to
the right with a speed of c, while f (x + c t) represents a traveling wave moving
to the left with a speed of c, we see that a solution to the wave equation (1.1)
essentially consists of two traveling waves moving in the opposite directions
with the same speed c.

2
 3. Traveling waves solutions of the Cauchy problems

Consider the Cauchy problem of the wave equation (1.1) with initial data

u(x, 0) = ϕ(x), (3.1)

ut (x, 0) = ψ(x). (3.2)

Suppose the solution u(x, t) is of the form (2.4), then, we have

ut (x, t) = c(f 0 (x + c t) − g 0 (x − c t)).

By setting t = 0, and by the initial data, we obtain the relations

ϕ(x) = f (x) + g(x), (3.3)

ψ(x) = c(f 0 (x) − g 0 (x)). (3.4)

On the other hand, it follows from (3.3) that

ϕ0 (x) = f 0 (x) + g 0 (x). (3.5)

Combining (3.4) and (3.5) yields

1 1
f0 = (c ϕ0 + ψ) and g 0 = (c ϕ0 − ψ),
2c 2c
and thus, by integration, we have
Z x
1 1
f (x) = c1 + ϕ(x) + ψ(s) ds, (3.6)
2 2c 0

Z x
1 1
g(x) = c2 + ϕ(x) − ψ(s) ds, (3.7)
2 2c 0

where c1 and c2 are integral constants. Since ϕ(x) = f (x) + g(x) (see (3.3)),
by (3.6) and (3.7), we have
c1 + c2 = 0.

As a conclusion of the above analysis, we obtain the so called d’Alembert

solution to the Cauchy problem of the wave equation
Z
1 1 x+c t
u(x, t) = (ϕ(x + c t) + ϕ(x − c t)) + ψ(s)ds. (3.8)
2 2c x−c t

3
Example 1. Consider the Cauchy problem

 2
utt = c uxx , x ∈ R, t > 0,
u(x, 0) = 0, x ∈ R, (3.9)

u (x, 0) = cos(x), x ∈ R.
t

Notice that here ϕ(x) = 0 and ψ(x) = cos(x), it follows from (3.8) that the
solution to (3.9) is given by
Z x+c t
1 1
u(x, t) = cos(s)ds = (sin(x + c t) − sin(x − c t))
2c x−c t 2c
1
= cos(x) sin(c t).
c
It is interesting to see that the solution separates the variables, i.e. the solution
can be written in the form of separation of variables
u(x, t) = X(x) T (t). (3.10)

4. Initial-boundary value problem and separation of variables

In the rest of this chapter, we consider the following initial-boundary value

problem of the wave equation
utt = uxx , ∀x ∈ (0, 1), ∀t > 0, (4.1a)
u(0, t) = u(1, t) = 0, ∀t > 0, (4.1b)
u(x, 0) = f (x), ut (x, 0) = g(x), ∀x ∈ (0, 1). (4.1c)
Here, for simplicity but without loss of generality, we have assumed that the
wave speed c = 1. In fact, in the general case, by a change of the time scale
t̃ = c t, we have
utt = c2 ut̃t̃ ,
and thus, in the new time variable t̃, the wave equation utt = c2 uxx is expressed
as
ut̃t̃ = uxx .
For the initial-boundary value problems, the solutions in the form of the
oppositely moving traveling waves are not as easily found as for the initial value

4
problems, since the boundary effect must be taken into account. For example,
the waves can not simply travel because of the existence of the boundaries,
they may be absorbed, reflected by the boundaries, and there may be incom-
ing waves traveling into the domain through the boundaries. In general, the
problem is much more complicated because of the boundary effects introduced
by the boundary conditions.
Let us try to find solutions of the form (3.10) to the initial boundary value
problem (4.1). Inserting the ansatz u(x, t) = X(x)T (t) into the wave equation
(4.1a), we obtain
T 00 (t) X 00 (x)
= . (4.2)
T (t) X(x)
Since the left hand side is independent of x, while the right hand side is inde-
pendent of t, both expressions must be independent of x and t, and therefore,
be a constant. That is
T 00 (t) X 00 (x)
= = −λ (4.3)
T (t) X(x)
for a suitable λ ∈ R. In particular, this and the boundary conditions (4.1b)
imply that the functions X(x) are the solutions to the eigenvalue problem
(
−X 00 (x) = λX(x), x ∈ (0, 1),
(4.4)
X(0) = X(1) = 0.

Hence, we have a family of solutions

Xk (x) = sin(kπx), k = 1, 2, . . . , (4.5)

with respect to the eigenvalues

λ = λk = (kπ)2 , k = 1, 2, . . . . (4.6)

For each eigenfunction Xk (x), the corresponding Tk (t) must satisfy

−Tk00 (t) = λk Tk (t) = (kπ)2 Tk (t), (4.7)

which has two linearly independent solutions given by

Tk (t) = eikπt and Tk (t) = e−ikπt .

5
Thus, the general real solution to the problem (4.7) is of the form

Tk (t) = ak cos(kπt) + bk sin(kπt), (4.8)

where ak , bk ∈ R are arbitrary constants. By combining (4.5) and (4.8), we

conclude that the functions
uk (x, t) = sin(kπx)(ak cos(kπt) + bk sin(kπt)) (4.9)

satisfy the wave equation (4.1a) and the boundary conditions (4.1b), and so
do the functions of their linear combinations
N
X
u(x, t) = sin(kπx)(ak cos(kπt) + bk sin(kπt)). (4.10)
k=1

In particular, u(x, t) given by (4.10) is a solution to the initial-boundary value

problem (4.1) with the initial data
N
X N
X
u(x, 0) = ak sin(kπx) and ut (x, 0) = kπbk sin(kπx). (4.11)
k=1 k=1

Hence, if the initial data can be expressed in the form

N
X N
X
f (x) = ak sin(kπx) and g(x) = b̃k sin(kπx), (4.12)
k=1 k=1

then, at least formally, the solution of the initial-boundary value problem (4.1)
is given by
∞
à !
X b̃k
u(x, t) = sin(kπx) ak cos(kπt) + sin(kπt) . (4.13)
k=1
kπ

Example 2. Consider the initial-boundary value problem (4.1) with

f (x) = x(1 − x) and g(x) = 0.

Recall that the Fourier coefficients are given by

Z 1
ak = 2 f (x) sin(kπx) dx, k = 1, 2, . . . ,
0

6
we have
∞
X 8
x(1 − x) = sin((2k − 1)πx),
k=1
π 3 (2k− 1)3

by (4.13), the formal solution of the problem (4.1) is given by

∞
X 8
u(x, t) = sin((2k − 1)πx) cos((2k − 1)πt).
k=1
π 3 (2k − 1)3

5. Energy arguments and uniqueness of solutions

Suppose that u ∈ C2 ([0, 1] × [0, ∞) is a solution to the initial-boundary

value problem (4.1). Then, the mechanical energy E(t) of the elastic string at
time t is given by
Z 1
¡ 2 ¢
E(t) = ux (x, t) + u2t (x, t) dx, (5.1)
0

where the first part is the elastic energy, and the second part is the kinetic en-
ergy. We are going to see how the energy evolves with time. By differentiating
E(t) with respect to t, we obtain
Z 1
0
E (t) = 2 (ux (x, t)uxt (x, t) + ut (x, t)utt (x, t)) dx. (5.2)
0

Using integration by parts for the first term, we obtain

Z 1 Z 1
ux (x, t)uxt (x, t) dx = ux (x, t)utx (x, t) dx
0 0
Z 1 Z 1
= ux (x, t)ut (x, t)|x=1
x=0 − uxx (x, t)ut (x, t) dx = − uxx (x, t)ut (x, t) dx,
0 0

where the last equality follows from the homogeneous boundary conditions at
x = 0, 1. Inserting this into (5.2) and noticing that u solves the wave equation
(4.1a), we otain
E 0 (t) = 0.
This implies
E(t) = E(0), ∀t ≥ 0. (5.3)

7
Hence, for the wave equation, the mechanical energy E(t) is preserved for all
time. Especially, if at the initial time t = 0, E(0) = 0, then E(t) ≡ 0.
Let u1 and u2 be two solutions of (4.1) with initial data (f1 , g1 ) and (f2 , g2 )
respectively, and let w = u1 − u2 . Then, it is easily verified that w is a solution
of (4.1) with initial data (f1 − f2 , g1 − g2 ). Hence, it follows from (5.3) that
Z 1 Z 1
£ 2 ¤ £ ¤
wx (x, t) + wt2 (x, t) dx = (u1 − u2 )2x (x, t) + (u1 − u2 )2t (x, t) dx
0 0
Z 1 £ ¤
= (f1 − f2 )2x (x) + (g1 − g2 )2 (x) dx. (5.4)
0

This equality implies the stability of the system with respec to the initial data,
more precisely, it states that if the initial data of the solutions are close to each
other, then the solutions will remain close for all time. In particular, it impleas
the following uniqueness result:

Theorem 5.1. If u1 , u2 ∈ C2 ([0, 1] × [0, ∞) are two solutions to the initial-

boundary value problem (4.1) with the same initial data, then u1 ≡ u2 .

Proof. Since the initial data are the same, it follows from (5.4) that u1,x =
u2,x and u1,t = u2,t . Hence, the two solutions can only differ by a constant.
However, since they have the same initial and boundary data, the constant
must be zero, and thus we have u1 ≡ u2 . ¤

6. Initial value problems of the convection equation

Consider certain solute carried by some incompressible fluid flowing in a

pipe with invariant cross section. Let u(x, t) be the density of the solute mea-
sured in mass per unit length, let c be the speed of the flow, which is assumed
to be a constant. If we neglect the effect of diffusion, then the difference of
the solute mass contained in the pipe section [x1 , x2 ] at time t2 and at time
t1 < t2 is the net silute mass carried into this section by the fluid through the
two end sides of the section, that is
Z x2 Z x2 Z t2 Z t2
u(x, t2 ) dx − u(x, t1 ) dx = c u(x1 , t) dt − c u(x2 , t) dt. (6.1)
x1 x1 t1 t1

8
Suppose that u(x, t) is suffciently smooth, then, we have
Z x2 Z t2
[ut (x, t) + cux (x, t)] dt dx, ∀x1 < x2 , t1 < t2 . (6.2)
x1 t1

This leads to to the convection equation

ut + cux = 0. (6.3)

The convection equation, also known as the advection equation, may be re-
garded as a first order wave equation. To see this, we introduce a new variable

η = x − c t,

and define
u(x, t) , u0 (x − c t) = u0 (η). (6.4)
Suppose that the function u0 is smooth, then, by the chain rule, we have

ut (x, t) = u00 (η)ηt = −c u00 (η) and ux (x, t) = u00 (η)ηx = u00 (η).

This implies that u(x, t) given by (6.4) satisfies the convection equation (6.3).
Notice that u(x, t) = u0 (x − c t) is a traveling wave moving at speed c, if c > 0,
it moves to the right, and if c < 0, it moves to the left, in particular, at the
initial time t = 0, u(x, 0) = u0 (x). Hence the solution to the initial value
problem of the convection equation
(
ut (x, t) + c ux (x, t) = 0, ∀x ∈ R, ∀t > 0,
(6.5)
u(x, 0) = u0 (x),

is simply given by u(x, t) = u0 (x − c t).

We see that, for a traveling wave solution u(x, t) = u0 (x − c t), the value of
the function is a constant along any line of the form X(t) = x0 +c t. This is the
characteristic property of the convection equation, and we call the equation

dX
(t) = c (6.6)
dt
the characteristic equation of the convection equation (6.3), and we call the
lines satisfying the characteristic equation the characteristic lines.

9
Example 3. For u0 (x) = sin(πx), the solution of the initial value problem of
the convection equation (6.5) is

1 ¡ iπ(x−c t) ¢
u(x, t) = sin(π(x − c t)) = e − e−iπ(x−c t)
2i
= sin(πx) cos(πc t) − cos(πx) sin(πc t).

Notice that the above solution can also be written in the form of a linear
combination of functions with separated variables. This inspired us to try to
solve the convection equation (6.3) by the method of separation of variables,
i.e. to find particular solutions of the form u(x, t) = X(x)T (t). Inserting
u(x, t) = X(x)T (t) into the equation, we obtain

X 0 (x) −1 T 0 (t)
= = iλ,
X(x) c T (t)

where λ ∈ R is an arbitrary constant, which lead to the solution of the form

u(x, t) = eiλx e−iλc t = eiλ(x−c t) . (6.7)

Remark 6.1. Of cause, we may generally obtain

X 0 (x) −1 T 0 (t)
= = µ + iλ,
X(x) c T (t)

which lead to the solution of the form

u(x, t) = eµ (x−c t) eiλ(x−c t) .

However, if µ 6= 0, the function would be unbounded in space and time. On

the other hand, a physical solution to the convection equation will remain
bounded for all time. That is why we only consider particular solutions of the
form (6.7).

Remark 6.2. There is no bounded real solutions in the form of separation of

variables. However, we have
1 ¡ iλ(x−c t) ¢
e − e−iλ(x−c t) = sin(λx) cos(λc t) − cos(λx) sin(λc t) = sin λ(x − c t),
2i
1 ¡ iλ(x−c t) ¢
e + e−iλ(x−c t) = cos(λx) cos(λc t) + sin(λx) sin(λc t) = cos λ(x − c t).
2

10
It follows from the suer-position principle, the Fourier series method can be
easily applied to solving the initial value problem of the convection equation.

Example 4. Consider the initial value problem of the convection equation

(6.5) with

∞
X
u0 (x) = {ak cos(kπx) + bk sin(kπx)} .
k=1

The traveling wave solution is given by

∞
X
u(x, t) = {ak cos(kπ(x − c t)) + bk sin(kπ(x − c t))}
k=1
∞
X
= {ak [cos(kπx) cos(kπc t) + sin(kπx) sin(kπc t)]
k=1

+bk [sin(kπx) cos(kπc t) − cos(kπx) sin(kπc t)]} .

7. Initial-boundary value problems of the convection equation

Since the solution of the convection equation (6.3) is a traveling wave

solution of the form u(x, t) = f (x − c t), we see that, for c > 0, the solution
at (x, t) always picks up its information from the left, while for c < 0, the
information comes from the right. In other words, the solution information
comes from the upwind direction. Hence, for a convection equation defined
on a finite interval, say (0, 1), only on the boundary in the upwind direction
should a boundary condition be imposed. For c > 0, a proper initial-boundary
value problem of the convection equation is given as



ut + c ux = 0, ∀x ∈ (0, 1), ∀t > 0,
u(x, 0) = u0 (x), ∀x ∈ (0, 1), (7.1)

u(0, t) = g(t), ∀t > 0.

11
While for c < 0, a proper initial-boundary value problem of the convection
equation is given as


ut + c ux = 0, ∀x ∈ (0, 1), ∀t > 0,
u(x, 0) = u0 (x), ∀x ∈ (0, 1), (7.2)

u(1, t) = g(t), ∀t > 0.

It is easily verified that the solution to the problem (7.1) is given by

(
u0 (x − c t), if x − c t > 0,
u(x, t) = (7.3)
g(− x−c
c
t
), if x − c t ≤ 0,

and the solution to the problem (7.2) is given by

(
u0 (x − c t), if x − c t < 1,
u(x, t) = (7.4)
g( 1−(x−c
c
t)
), if x − c t ≥ 1,

as long as u0 and g are smooth and satisfy certain consistent conditions (for
example, if u0 (0) = g(0) for c > 0 and u0 (1) = g(0) for c < 0, then the above
defined functions are globally continuous).

8. Finite difference schemes for the convection equation

The simplest finite difference scheme is the so called upwind scheme:

Ujn+1 − Ujn Ujn − Uj−1

n
+c = 0, if c > 0, (8.1)
4t 4x

and
Ujn+1 − Ujn n
Uj+1 − Ujn
+c = 0, if c < 0, (8.2)
4t 4x

where ut (xj , tn ) is approximated by the first order forward finite difference in

the temporal variable t with the time step 4t

Ujn+1 − Ujn
ut (xj , tn ) ≈ ,
4t

12
and ux (xj , tn ) is approximated by the first order backward (c > 0) or forward
(c < 0) finite difference in the space variable x with the space step 4x
Ujn − Uj−1
n n
Uj+1 − Ujn
ux (xj , tn ) ≈ , if c > 0; ux (xj , tn ) ≈ , if c < 0.
4x 4x
The key point here is that the information needed to compute the unknown
Ujn+1 is picked up from the upwind direction.
To see how accurate the convection equation (6.3) is approximated by the
upwind scheme, substitute a smooth solution of the equation into the scheme,
and by applying the Taylor expasions at (xj , tn )
1
un+1
j = unj + 4t ut (xj , tn ) + (4t)2 utt (xj , tn ) + O((4t)3 ),
2
1
unj−1 = unj − 4x ux (xj , tn ) + (4x)2 uxx (xj , tn ) + O((4x)3 ),
2
1
unj+1 = unj + 4x ux (xj , tn ) + (4x)2 uxx (xj , tn ) + O((4x)3 ),
2
we obtain

un+1
j − unj unj − unj−1 1
+c = [ut + c ux ]nj + [4tutt − 4xc uxx ]nj
4t 4x 2
+ O((4t)2 + (4x)2 ), if c > 0, (8.3)

un+1
j − unj unj+1 − unj 1
+c = [ut + c ux ]nj + [4tutt + 4xc uxx ]nj
4t 4x 2
+ O((4t)2 + (4x)2 ), if c < 0, (8.4)

In fact, for any smooth function u(x, t), donote r = 4t/4x, denote 4+x and
4−x the first order forward and backward finite difference operators for the
space variable, denote 4+t the first order forward finite difference operator for
the time, we have
£ ¤
((4t)−1 4+t + c (4x)−1 4−x ) − (∂t + c ∂x ) unj
1
= − c 4x (1 − c r) [uxx ]nj + O((4x)2 ), if c > 0. (8.5)
2

13
 £ ¤
((4t)−1 4+t + c (4x)−1 4+x ) − (∂t + c ∂x ) unj
1
= c 4x (1 + c r) [uxx ]nj + O((4x)2 ), if c < 0, (8.6)
2

It follows from the definition that the truncation error of the upwind scheme
of the convection equation is O(4x), and the leading term of the truncation
error is
1
− 4x|c|(1 − |c| r)uxx . (8.7)
2

To obtain an accurate numerical approximation to the solution of the con-

vection equation, it is necessary to have a good approximation of the equation,
in other words, the truncation error must be small. However, a smaller trun-
cation error does not necessarily imply a better numerical approximation to
the solution, since, as shown by many numerical experiments, in some circum-
stances a small initial error may propagate and grow extremely fast, and may
eventually dominate the numerical results, while in some other circumstances
the growth of the error will be well under control.
A key concept involved here is the so called the stability of a scheme.
A simple but very important stability condition named as CFL condition
(Courant-Friedrichs-Lewy) comes from a direct observation, which states that
the dependent domain of the numerical solution must cover, at least in the
sense of limit, that of the analytical solution. Since the analytical solution of
the convection equation is completely determined by the characteristic lines,
the CFL condition for the upwind scheme can simply be expressed as

4t
0 ≤ |c| = |c| r ≤ 1. (8.8)
4x

Recall that the convection equation has traveling wave solutions eikπ(x−c t) ,
k = 0, 1, 2, · · · , as a sequence of particular solutions (see (6.7)). We may as
well work out a sequence of particular solutions that a finite difference scheme
has, and compare them with the corresponding traveling wave solutions. We
seek solutions of separation of variables of the form
arg(λk )
Ujn = λnk eikπj4x = |λk |n eikπ(j4x+ c kπ 4t c n4t) .

14
Substitute it into the upwind scheme (8.1) and (8.2), we see that the amplifi-
cation factor λk must satisfy

λk = (1 − |c| r) + |c| re−i sign(c) kπ4x . (8.9)

This yields

|λk |2 = [(1 − |c| r) + |c| r cos kπ4x]2 + [|c| r sin kπ4x]2

= (1 − |c| r)2 + (|c| r)2 + 2|c| r(1 − |c| r) cos kπ4x
= 1 − 2|c| r(1 − |c| r)(1 − cos kπ4x)
1
= 1 − 4|c| r(1 − |c| r) sin2 kπ4x. (8.10)
2
Hence, we have
|λk | ≤ 1, ∀k ⇔ |c| r ≤ 1, (8.11)
which implies that the upwind scheme is stable if and only if the CFL condition
is satisfied.
π
In fact, if |c| r > 1, then, for k̂π = 4x
, we have

λk̂ = 1 + 4|c| r|(1 − |c| r)| > 1,

hence the corresponding particular solution Ujn = λnk̂ eik̂πj4x grows exponen-
tially fast. Noticing that this is also the fastest growing frequency, we con-
clude that in such a case the numerical solutions are unstable and will soon
be dominated by the exponentially growing mesh size oscillation, since errors
with all frequencies will inevitably appear in numerical computations.
On the other hand, if |c| r < 1, then, it follows from (8.10) that there exists
damping effect for all frequencies, and the higher the frequency the stronger
1
the damping effect (for k ≤ k̂ = 4x
). For low frequencies with k4x ¿ 1, since
sin kπ4x ≈ kπ4x in such a case, we have
1
|λk | ≈ 1 − |c| r(1 − |c| r)k 2 π 2 (4x)2 ,
2
which means a drop of the amplitude by O((4x)2 ) in each time step. Compare
with the corresponding traveling wave solutions, we see that the amplitude of
the numerical solutions has a global error of O(4x).

15
 To see how accurate the numerical solution recovers the speed of the travel-
ing waves, we need to calculate arg(λk ), which represents the phase increment
in a time step. By (8.9), we have
µ ¶
− sign(c) |c| r sin kπ4x
arg(λk ) = arctan . (8.12)
(1 − |c| r) + |c| r cos kπ4x
For low frequencies with k4x ¿ 1, by Taylor expansion, we have
· ¸
1 2
arg(λk ) = − sign(c) |c| rkπ4x 1 − (1 − |c| r)(1 − 2|c| r)(kπ4x) + · · ·
6
· ¸
1 2
= −c kπ4t 1 − (1 − |c| r)(1 − 2|c| r)(kπ4x) + · · · .
6
Hence, the corresponding wave travels at speed
· ¸
1 2
ck = c 1 − (1 − |c| r)(1 − 2|c| r)(kπ4x) + · · · . (8.13)
6

This shows that the error of the speed (or the phase angle) is O((k4x)2 ).
Compare with the error of the amplitude, which is O(k4x), this is almost
negligible. We also notice that, for ν = c r < 21 , ck < c, this correspond to a
phase lag, while for ν = c r > 12 , ck > c, this correspond to a phase adavance.

Remark 8.1. The upwind scheme can be modified to obtain higher order finite
difference schemes. For example, we may use the second order central finite
difference operator to replace the second order differential operator in (8.5) or
(8.6). Since
unj+1 − 2unj + unj−1
2
= [uxx ]nj + O((4x)2 ),
(4x)
the induced finite difference scheme
1 1
Ujn+1 = c r(1 + c r)Uj−1
n
+ (1 − c2 r2 )Ujn − c r(1 − c r)Uj+1
n
, (8.14)
2 2
known as the Lax-Wendroff Scheme, is of order ((4x)2 ). It can be shown that,
for the Lax-Wendroff scheme, the global error of the amplitude is O((k4x)3 ),
while the error of the speed (or the phase angle) is O((k4x)2 ). Obviously, in
this case, compared with the error of the speed, the error of the amplitude is

16
almost negligible. In other words, for the Lax-Wendroff scheme, the error of
the speed is of main concern. In fact, we may observe some artificial oscilla-
tions in numerical solutions, which is caused by the numerical dispersion (a
phenomenon which reflects the fact that waves of different frequencies travel
in different speed).

Remark 8.2. Let a linear numerical scheme of the convection equation be given
as
U n+1 = N (U n )
with the truncation error
un+1 − N (un )
τn = .
4t
Let en = un −U n be the error of the numerical solution. Suppose, the numerical
scheme is stable in the sense that
kN (V )k ≤ kV k, ∀4x > 0 and for all mesh function V ,
4t
whenever the ratio of the mesh size r = 4x
satisfies a certain specified stability
condition. Then, we have

ken+1 k = kun+1 − U n+1 k = kun+1 − N (un ) + N (un ) − U n+1 k

n
X
n n 0
≤ kN (e )k+4tkτ k ≤ ke k+4t kτ i k ≤ ke0 k+T τmax , ∀(n+1)4t ≤ T.
i=0

This implies that, for a stable linear numerical scheme, the error of the numer-
ical solution is of the same order as its truncation error, provided its initial
data is approximated accordingly.
In the stability analysis using Fourier modes, the norm kU n k = kU n k2 ,
P
( j |Ujn |2 4x)1/2 is used. The idea can be easily applied to other equations and
P
other norms. For example, L1 -norm kU n k1 , ( j |Ujn |4x)1/2 and L∞ -norm
kU n k∞ , supj |Ujn |.

Remark 8.3. CFL condition is a necessary condition for the stability. However,
in genergal, it is not a sufficient condtion. For example, the scheme using the
central difference instead of the one-sided upwind difference is always unstable.

17
 9. Finite difference method for the wave equation

Naturally, using the second order central finite difference operator to ap-
proximate the second order differential operators for both the time and the
space variables, we obtain the following second order finite difference scheme
for the wave equation

Ujn+1 − 2Ujn + Ujn−1 U n − 2Ujn + Uj−1

2 j+1
n
− c = 0, (9.1)
(4t)2 (4x)2

or
Ujn+1 = ν 2 Uj−1
n
+ 2(1 − ν 2 )Ujn + ν 2 Uj+1
n
− Ujn−1 , (9.2)

where ν = c r = c 4t/4x. By substituting a smooth solution of the wave

equation into the finite difference scheme (9.1), and by the Taylor expansion,
it is easily verified that the truncation error of the scheme is O((4t)2 +(4x)2 ).
To start the scheme, we need to provide initial data on the first two time
levels n = 0, 1, which may be obtained by certain proper manipulations of the
initial data (4.1c). A simplest set of discrete initial data is given by

Uj0 = f (xj ), and Uj1 = f (xj ) + 4t g(xj ). (9.3)

A better discrete initial data can be obtained by using

1 1
u(xj , s) = u(xj , 0) + s ut (xj , 0) + s2 utt (xj , 0) + s3 uttt (xj , 0) + O(s4 )
2 6
1 1
= f (xj ) + s g(xj ) + s2 c2 uxx (xj , 0) + s3 c2 uxxt (xj , 0) + O(s4 )
2 6
1 1
= f (xj ) + s g(xj ) + s2 c2 fxx (xj ) + s3 c2 gxx (xj ) + O(s4 ).
2 6

Recall that there are two traveling waves with speed c and −c, so there
are two characteristic lines (X 0 = ±c) passing through each point (x, t) (see
also the form of the d’Alembert solution (3.8)), hence the CFL condition for
the scheme (9.1) is
c 4t
ν , cr = ≤ 1. (9.4)
4x

18
 We may also seek solutions of separation of variables of the form
arg(λk )
Ujn = λnk eikπj4x = |λk |n eikπ(j4x+ c kπ 4t c n4t) .

Substitute it into the scheme (9.2), we see that the amplification factor λk
must satisfy

λ2k − (ν 2 e−i kπ4x + 2(1 − ν 2 ) + ν 2 ei kπ4x )λk + 1 = 0,

which may be simplified as

λ2k − (2 − sk )λk + 1 = 0, (9.5)

where sk = 4ν 2 sin2 ( kπ4x

2
) ∈ (0, 4ν 2 ]. Hence, we obtain
p
2 − sk ± sk (sk − 4)
λk = . (9.6)
2
Lemma 9.1. Let sk ≥ 0 be given. Let λ±
k be the roots of (9.5). Then

|λ+ −
k | ≤ 1 and |λk | ≤ 1 ⇔ sk ≤ 4. (9.7)

Proof. If sk = 0, then λ± ±
k = 1. If sk = 4, then λk = −1.

If sk ∈ (0, 4), then λ+ −

k and λk are conjugate complex numbers, that is

λ+
k = ρk e
i θk
and λ−
k = ρk e
−i θk

for some ρk > 0 and θk ∈ (0, π). Furthermore, it follows from (9.5) that the
product of the roots is 1, i.e.

λ+ − 2
k λk = ρk = 1.

Hence, |λ+ −
k | = |λk | = ρk = 1.

If sk > 4, then λ+ − + −
k and λk are two distinct real roots satisfying λk λk = 1.
Hence, one of them must have absolute value greater than 1. ¤

As a consequence of the lemma, we conclude that the scheme (9.1) is stable

if and only if the CFL condition is satisfied. It is interesting to notice that in
this case there are two traveling waves moving in the opposite direction with
the same speed, and in particular there is no damping and thus no error to the

19
amplitude of the traveling waves (unlike the upwind scheme for the convection
equation). On the other hand, we have
 q 
2ν sin 21 kπ4x 1 − ν 2 sin2 21 kπ4x
arg λ±
k = ± arctan
 
1 − 2ν 2 sin2 12 kπ4x
µ ¶
1
= ±c kπ4t 1 − (1 − ν 2 )(kπ4x)2 + · · · . (9.8)
24

Hence, in general, the error to the speed (the relative phase error) is O((k4x)2 ),
and the numerical waves always travel at speeds slower than that of the an-
alytical waves. Since there is no damping, the phase lag effect of the scheme
should be apparently observable in certain cases.

10. Implicit finite difference schemes of the wave equation

We may also consider the implicit finite difference schemes of the wave
equation (1.1) in the following form
"
Ujn+1 − 2Ujn + Ujn−1 2
n+1
Uj+1 − 2Ujn+1 + Uj−1
n+1
= c θ
(4t)2 (4x)2
#
n−1
n
Uj+1 n
− 2Ujn + Uj−1 Uj+1 − 2Ujn−1 + Uj−1
n−1
+(1 − 2θ) +θ , (10.1)
(4x)2 (4x)2

where θ ∈ (0, 1] is a parameter. It is obvious, for θ = 0, the scheme (10.1)

reduces to the explicit scheme (9.1).
It is not difficult to work out that the truncation error of the scheme (10.1)
is of order O((4t)2 + (4x)2 ).
Substitute the Fourier modes
arg(λk )
Ujn = λnk eikπj4x = |λk |n eikπ(j4x+ c kπ 4t c n4t) .

into the scheme (10.1), we obtain the equation for the amplification factor
¡ ¢
λ2k − 2λk + 1 = (θν 2 λ2k + (1 − 2θ)λk + θν 2 ) eikπ4x − 2 + e−ikπ4x ,

20
or equivalently
µ ¶
4ν 2 sin2 12 kπ4x
λ2k − 2− λk + 1 = 0. (10.2)
1 + 4θν 2 sin2 12 kπ4x

4ν 2 sin2 21 kπ4x
Denote sk = 1+4θν 2 sin2 21 kπ4x
, and by Lemma 9.1, we conclude that the scheme
(10.1) is stable if
(
(1 − 4θ)ν 2 ≤ 1, ∀θ ∈ (0, 14 );
(10.3)
unconditional, ∀θ ∈ [ 41 , 1].

On the other hand, we have

 q 
2ν sin 12 kπ4x 1 + ν 2 (4θ − 1) sin2 21 kπ4x
arg λ±
k = ± arctan
 
1 + 2ν 2 (2θ − 1) sin2 12 kπ4x
µ ¶
1
= ±c kπ4t 1 − (1 + (12θ − 1)ν 2 )(kπ4x)2 + · · · . (10.4)
24

Hence, in general, there is a phase lag of order O((k4x)2 ), and it is obvious

that a smaller θ produces a smaller relative phase error.
1
We may conclude from the above analysis that θ = 4
should be the most
practical choice.

11. Matrix method for stability analysis

If the boundary condition is nonhomogeneous and nonperiodic, then the

stability analysis using the Fourier modes does not work, since the Fourier
modes are no longer the particular solutions. However, a two level linear
scheme U n+1 = N (U n ) can always be written in the form

U n+1 = An U n + bn .
To show the stability, we only need to check that, for some matrix norm k · k,
Q
there exists a constant C such that k ni=0 Ai k ≤ C, for all n. For example, if
Ai = A for all i, then, we only need to show that

kAn k ≤ C, ∀n.

21
In particular, if
kAk ≤ 1,
then, the scheme is stable.
Notice that for different mesh size 4x, the size of the matrix A is also
different, and its norm kAk is in general a function of 4x and 4t.

22