1B Methods 98

1B METHODS
LECTURE NOTES
Richard Jozsa, DAMTP Cambridge
[email protected]
October 2013

PART IV:
PDEs on
unbounded domains
1B Methods 99

9 CHARACTERISTICS

9.1 Well-posed problems

In chapters 3,4,5 our development of PDEs on bounded domains was based largely on
physical considerations, not just in the origin of the equations themselves, but also in
motivating the various kinds of BCs and ICs that are necessary and sufficient to guarantee
a unique solution.
• For the 1D wave equation on a finite string we specified the boundary values y(0, t)
and y(L, t) for all t as well as the initial displacement y(x, 0) and initial velocity yt (x, 0).
• For the heat equation in a finite bar for temperature distribution θ(x, t) we specified
the boundary values θ(0, t) and θ(L, t) and the initial distribution θ(x, 0) but not the
initial time derivative.
• For the Laplace equation (e.g. for steady heat conduction) all that was required was the
boundary value. More generally for the Laplace equation on a domain D with boundary
∂D (in 1D, 2D or 3D) the following is standard terminology: solve for φ throughout D
with –
Dirichlet problem: φ being specified on the boundary ∂D;
Neumann problem: the normal derivative n.∇φ being specified on the boundary.
It can be shown that the Dirichlet problem has a unique solution whereas the solution
of the Neumann problem is unique up to just an additive constant.
This illustrates that the exact nature of data that is necessary and sufficient for a unique
solution depends non-trivially on the kind of equation being considered. Also if the
equation is just viewed mathematically in terms of all its independent variables (without
identifying space or time) then an IC is just another kind of BC (for the variable t). If we
have a PDE for φ with auxiliary data (values of φ and its derivatives) specified on some
surface (e.g. along a line in 2D or on a surface in 3D) then this data is called Cauchy
data and the problem is called a Cauchy problem.
In generality, a problem comprises a PDE plus some auxiliary data. The problem is
said to be well-posed (in the sense of Hadamard) if three conditions hold:
(WP1): a solution exists;
(WP2): the solution is unique;
(WP3): the solution depends continuously on the auxiliary data. (Intuitively speaking a
small change in the data should result in a small change in the full solution; but a rigorous
statement of this condition requires a notion of nearness (topology) and continuity in
spaces of functions, which we don’t develop in this course).
So far we have always considered well-posed problems. Note that (WP1) can be violated if
we attempt to impose too many auxiliary conditions and (WP2), if we don’t make enough
demands on the solution. (WP3) is an interesting additional requirement, motivated
again by physical considerations – since physical systems cannot be measured or set up
with infinite precision, we would want the solution to be suitably stable under small
perturbations of the data, if the equation is to usefully model the physical situation.
(But the mathematical subject of chaos theory studies precisely the opposite behaviour!)
1B Methods 100

Example. An intuitive example of an ill-posed problem that satisfies (WP1) and (WP2)
but not (WP3) is the “backwards-in-time” heat equation. Consider the heat equation
for u(x, t) on 0 ≤ x ≤ L with BCs specified at x = 0, L. Usually we specify an IC of
u at t = 0. Diffusion is a “smoothing” process so we expect that small perturbations
in u(x, 0) will not grow, and the problem is indeed well-posed. But suppose instead we
specify u(x, T ) at some (late) time T and ask for u(x, 0) at an earlier time t = 0 say. This
problem can be shown to violate (WP3) – intuitively if u(x, 0) has a complex detailed
structure e.g. with high spatial gradients, they will be quickly smoothed out by heat flow
gradients and tend to leave only a small imprint on u(x, T ) at the later time T i.e. small
perturbations in the latter will back-evolve into large changes in the former. Note also
that if (as physically usual) we have an asymptotic steady state, all initial conditions will
evolve to become very close to it, and hence remain encoded only in arbitrarily small
perturbations of the steady state.
Example. (An ill-posed Laplace equation problem).
Consider the following problem – we have Laplace’s equation uxx + uyy = 0 on the upper
half plane y ≥ 0 and −∞ < x < ∞, with BCs

u(x, 0) = f (x) and uy (x, 0) = g(x)

where f and g are specified functions. If f (x) = f1 (x) = 0 and g(x) = g1 (x) = 0 then we
have the solution u1 (x, y) = 0 which can be shown to be unique.
If now we set f (x) = f2 (x) = 0 and g(x) = g2 (x) = sinAAx then again we have a unique
solution (which can be found by looking for a separated variable solution u(x, y) =
X(x)Y (y) etc.) given by
sin Ax sinh Ay
u2 (x, y) = .
A2
Now if we consider A → ∞ we have f2 → f1 (actually they’re equal) and g2 → g1
but |u1 − u2 | can become arbitrarily large e.g. at x = π/(2A) we have u2 (x, y) =
(sinh Ay)/A2 ≈ eAy /A2 → ∞ as A → ∞. Hence the problem is ill-posed.

9.2 Characteristics for first order PDEs

We consider the general 1st order linear (actually so-called quasi-linear because of RHS
dependence on u) PDE in 2D :

α(x, y)ux + β(x, y)uy = γ(x, y, u) (1)

(where the subscripts denote partial derivatives), together with Cauchy data specifying
u(x, y) along a specified “initial” curve B in the xy plane. B will be described paramet-
rically (parameter t):
x = xB (t) y = y B (t)
and u along this curve is a specified function h(t):

u(xB (t), y B (t)) = h(t).

1B Methods 101

The homogeneous case γ = 0

Consider first the homogeneous case γ = 0 in eq. (1)
αux + βuy = 0. (2)
Introducing the vector field A(x, y) = (α(x, y), β(x, y)) (depending on the PDE but not
the Cauchy data) we see that eq. (2) is
A.∇u = 0. (3)
Interlude: facts about parameterised curves and vector fields
(a) if (x = xC (s), y = y C (s)) is any parameterised curve C, its tangent vector at point
labelled by s is T = (dxC /ds, dy C /ds).
(b) if u(x, y) is a function on the plane, its restriction to C is u(xC (s), y C (s)) and the
directional derivative along C is (by the chain rule)
du dxC dy C
= ux + uy = T .∇u.
ds ds ds
Hence if T .∇u = 0 (as a function of s) then u is constant along the curve C.
(c) if A(x, y) = (α(x, y), β(x, y)) is any (suitably regular) vector field on (a part of) R2 ,
its integral curves area 1-parameter family of non-intersecting curves filling (that part
of) R2 , defined by the requirement that at any point (x, y), the curve through (x, y) has
tangent A(x, y) (e.g. if we think of A as the velocity field of a fluid then the integral curves
are the flow lines of the fluid elements). More explicitly let x = xB (t) y = y B (t) be some
specified curve B that is transverse to A i.e. along B the tangent vector (dxB /dt, dy B /dt)
is nowhere parallel to A at the same point. We then label the integral curves of A by
t, the point where they intersect B, and use parameter s along the integral curves, with
s = 0 at B. Thus the tth curve is x = x(t, s) y = y(t, s) (here t is a curve-label and s is
the parameter along the curve) satisfying
dx dy
= α(x, y) = β(x, y) (4)
ds ds
subject to ODE initial conditions x(t, 0) = xB (t) and y(t, 0) = y B (t) (stating that the
curves pass through B when the parameter s is 0).
Now back to our PDE
For αux + βuy = 0 the integral curves of the vector field (α, β) are called the charac-
teristic curves of the PDE. Then eq. (3) just says that u(x, y) is constant along the
characteristic curves C. Now taking B of interlude item (c) above, to be the Cauchy
data curve we see that the Cauchy data values h(t) along the curve B are propagated
constantly along the characteristic curves to define u(x, y) in the plane (or at least in a
neighbourhood of the curve B) i.e. u(xC (t, s), y C (t, s)) = h(t) is constant in s. Finally
to get an explicit expression for u(x, y) as a function of our original co-ordinates (x, y),
we view the characteristic curve equations (for the tth curve parameterised by s)
x = xC (t, s) y = y C (t, s)
1B Methods 102

as a co-ordinate transformation (x, y) ↔ (t, s). If the Jacobian J = xC C C C

t ys − xs yt is
non-zero then the transformation can be inverted, solving for (t, s) in terms of (x, y):

t = t(x, y) s = s(x, y)

and then
u(x, y) = h(t(x, y))
gives the solution of our Cauchy problem.
Note the following features of the above construction:
(a) If any characteristic curve intersects the initial curve B more than once then the
problem is over-determined – the value of h(t) must then be constrained to be the same
at all such multiple intersection points for a solution to exist (cf (WP1)).
(b) If the initial curve B itself is a characteristic curve then we must have h(t) = const
for a solution to exist, and in that case the solution is not unique (cf (WP2)) as it can
be freely specified on all the other characteristics.
(c) If the initial curve is transversal to all characteristics and intersects them once only,
then the problem is well-posed for any h(t) and has a unique solution u(x, y) (at least
in a neighbourhood of B). Note that the initial data cannot be propagated from one
characteristic to another, so, for example, discontinuities in the initial data propagate
along the corresponding characteristic curve.
In summary,
to solve αux + βuy = 0 with u(x, y) = h(t) on initial curve (xB (t), y B (t)):
(1) Write down the equations eq. (4) for the characteristics with the s = 0 parameter’s
initial condition given by the initial curve B i.e. for each t we have the system of two
ODEs:
dx dy
=α = β with x(t, 0) = xB (t), y(t, 0) = y B (t).
ds ds
(2) Solve for the characteristics x = xC (t, s) y = y C (t, s).
(3) Algebraically invert these relations to obtain t = t(x, y) and s = s(x, y).
(4) Using the Cauchy data h(t) construct the solution as u(x, y) = h(t(x, y)), which
represents the initial data propagated constantly along the characteristics off B.
Example. (the simplest possible example!)
Consider ux (x, y) = 0 with Cauchy data u(0, y) = h(y) on the y axis. Without any fancy
theory of characteristics etc., the solution to ux = 0 is obviously u(x, y) = f (y) i.e. an
arbitrary function of y only, and then the Cauchy data immediately gives u(x, y) = h(y).
But it’s instructive to identify the characteristics in this example. We have (α, β) =
(1, 0), a field of constant horizontal unit vectors, and initial curve B (the y axis) can be
parameterised as (xB (t) = 0, y B (t) = t) with u = h(t) along B. So the characteristics are
clearly horizontal lines parallel to the x axis. Formally the tth characteristic curve goes
through the point t on B viz. (x, y) = (0, t) and has
dx dy
=α=1 =β=0 x(0) = 0, y(0) = t
ds ds
1B Methods 103

so x(t, s) = s + c1 and y(t, s) = c2 , and the initial s = 0 conditions give c1 = 0, c2 = t i.e.

x(t, s) = s and y(t, s) = t as expected (the tth line being the horizontal line at height t).
Inverting these relations gives s(x, y) = x and t(x, y) = y so u(x, y) = h(t(x, y)) = h(y)
as expected.
Example. Consider
ex ux + uy = 0 with u(x, 0) = cosh x.
The initial curve B is the x axis (x = t, y = 0) and the initial data along B is h(t) = cosh t.
We have (α, β) = (ex , 1) so the characteristics, labelled by t, satisfy

dx dy
= ex = 1, x(0) = t y(0) = 0
ds ds
so e−x = −s + c1 and y = s + c2 where the integration constants depend on t via the
initial conditions, which give c1 = e−t and c2 = 0 so e−x = e−t − s and y = s. Inverting
these expressions we get
s=y t = − log (y + e−x )
and the solution to the Cauchy problem is

u(x, y) = h(t(x, y)) = cosh − log(y + e−x )

 

whose validity you can verify directly (by substitution).

9.3 Inhomogeneous problems and characteristics

For the general quasi-linear case eq. (1) viz. α(x, y)ux + β(x, y)uy = γ(x, y, u) together
with u(x, y) = h(t) on curve B: x = xB (t), y = y B (t), we still have characteristic curves
exactly as before (not depending on h or γ) and writing u(s) = u(x(s), y(s)) for u along
any characteristic curve, eq. (1) states precisely that (for each curve labelled by t and
parameterised by s):
du
= γ(x, y, u) u(0) = h(t). (5)
ds
Thus the only difference from our previous (homogeneous) situation is that now, u is not
propagated as a constant along characteristic curves off from B, but instead we have to
solve the ODE eq. (5) to obtain u(s) for each curve t i.e. to obtain u(t, s). Then just as
before, we invert the relations x = x(t, s) and y = y(t, s) that define the characteristics,
to get s = s(x, y) and t = t(x, y) and finally our solution is u(x, y) = u(t(x, y), s(x, y)).
Example. Consider the Cauchy problem

ux + 2uy = yex with u = sin x when y = x.

Thus α = 1, β = 2 γ = yex and the initial curve B is (x = t, y = t), with initial data
h(t) = sin t. We first solve for the characteristic curves:

dx dy
=1 = 2, x(0) = t y(0) = t
ds ds
1B Methods 104

giving x = s + t and y = 2s + t. Thus along the tth characteristic curve we have

du
= γ(x, y, u) = yex = (2s + t)es+t with u(0) = h(t) = sin t.
ds
The solution is (exercise, but remember that t is just a parameter!)

u(s, t) = (2 − t)et + sin t + es+t (t + 2s − 2). (6)

Next inverting the characteristic curve equations gives

s=y−x t = 2x − y

and substituting into eq. (6) gives the solution

u(x, y) = (2 − 2x + y)e2x−y + sin(2x − y) + ex (y − 2).

9.4 Classification of second order linear PDEs

We can usefully generalise the idea of characteristics to apply to some classes of second
order PDEs. Consider the general second-order linear PDE (in two variables)

a(x, y)uxx + 2b(x, y)uxy + c(x, y)uyy + d(x, y)ux + e(x, y)uy + f (x, y)u = g(x, y) (7)

with Cauchy data u = h(t), ux = m(t), uy = n(t) specified along some curve B given
parametrically as x = xB (t), y = y B (t). Note that if we differentiate u along t we get

dxB dy B dxB dy B
h0 (t) = ux + uy = m(t) + n(t)
dt dt dt dt
giving a relation between h, m and n. Hence no more than two of these functions can be
freely specified.
In terms of the coefficient functions for the second derivatives (and note the conventionally
used extra factor of 2 in the uxy term), we introduce the following classification of linear
2nd order PDEs. This will be significant for the behaviour of solutions, and especially for
their relation to the Cauchy data, via an associated notion of characteristic curves. The
equation (7) is called:
• hyperbolic if b2 − ac > 0.
• parabolic if b2 − ac = 0,
• elliptic if b2 − ac < 0.
Note that the wave equation is hyperbolic (c = 1, b = 0, a = -wave speed2 ), the diffusion
equation is parabolic (a = 0, b = 0, c = −D) and the Laplace equation is elliptic
(a = c = 1, b = 0).
Note that in general, a, b and c are functions of (x, y) so a single equation can be hyper-
bolic, parabolic and elliptic in different parts of the plane.
1B Methods 105

Now consider introducing new independent variables ξ, η by the substitution

ξ = φ(x, y) η = ψ(x, y)

transforming the equation into one of the same type:

A(ξ, η)uξξ + 2B(ξ, η)uξη + C(ξ, η)uηη + · · · = G(ξ, η).

The substitution is straightforward but cumbersome to calculate. For example we have

ux = uξ φx + uη ψx so

uxx = [uξξ φx + uξη ψx ] φx + uξ φxx + [uηξ φx + uηη ψx ] ψx + uη ψxx

etc. and in particular we get the new 2nd derivative coefficients to be

A = aφ2x + 2bφx φy + cφ2y

B = aφx ψx + b(φx ψy + φy ψx ) + cφy ψy
C = aψx2 + 2bψx ψy + cψy2 .

Note that A (and/or C) can be made zero if φ (and/or ψ) can be chosen to satisfy

aφ2x + 2bφx φy + cφ2y = 0

i.e. aM 2 + 2bM + c = 0 for M = φx /φy

√
b2 − ac
−b ±
i.e. M= .
a
In that case the corresponding co-ordinate curves ξ = φ(x, y) = const (and/or η =
ψ(x, y) = const) are called characteristic curves of the original PDE. Note that if
φ(x, y) = const defines a curve y = y(x) then
dy dy φx
φx + φy =0 so =− .
dx dx φy
Hence in terms of a, b and c the characteristic curves are given by solving
 √ 
dy −b ± b2 − ac
=−
dx a
(obtaining the solution in the form f (x, y) = const and then using f for φ or ψ.)
For elliptic equations there are no real characteristics.
dy
For parabolic equations we get a single family of characteristics dx
= ab .
For hyperbolic equations we get two families of characteristics corresponding to the two
roots √
−b ± b2 − ac
M± =
a
and using the corresponding change of variables

ξ = φ(x, y) η = ψ(x, y)
1B Methods 106

with co-ordinate lines being the characteristic curves, the hyperbolic equation takes the
form
uξη + D(ξ, η)uξ + E(ξ, η)uη + F (ξ, η)u = G(ξ, η)
having no uξξ or uηη term, which is called the canonical form of the hyperbolic
equation.
Example. The equation uyy −xyuxx = 0 has a = −xy, b = 0 and c = 1 so b2 −ac = xy so
the equation is hyperbolic in the first (x > 0, y > 0) and third (x < 0, y < 0) quadrants,
elliptic in the second and fourth quadrants and parabolic along the axes x = 0 or y = 0.
In the hyperbolic region we have
√
−b ± b2 − ac 1
= ±√
a xy

so the two families of characteristics are given by

dy 1 y 3/2
= ±√ i.e. ± x1/2 = c
dx xy 3

and the substitution

y 3/2 1/2 y 3/2
ξ= +x η= − x1/2
3 3
(after a lengthy but straightforward calculation, that we omit!) will reduce the equation
to canonical form in the hyperbolic region.

9.5 D’Alembert’s general solution of the wave equation

An especially important example of a hyperbolic equation and its characteristics is the

(1D) wave equation

∂ 2u 2
2∂ u
= c with u(x, 0) = φ(x), ut (x, 0) = ψ(x)
∂t2 ∂x2
which is hyperbolic in the whole (x, t) plane. The characteristics are easily calculated to
be x ± ct = const so we introduce new variables

ξ = x + ct η = x − ct

and the equation takes an especially simple canonical form

∂ 2u
= 0.
∂ξ∂η
The general solution is easily obtained by first integrating w.r.t. ξ giving ∂u/∂η = F (η)
and then Z η
u= F (y)dy + g(ξ) = f (η) + g(ξ) (8)

for arbitrary functions f and g (g(ξ) being the η-integration constant for each ξ).
1B Methods 107

Thus the solution is the sum of two terms that are constant respectively on the two
families of characteristics. Let us now impose our initial conditions on the general solution
eq. (8). At t = 0 (recalling that ξ, η = x ± ct) we get

u(x, 0) = f (x) + g(x) = φ(x)

ut (x, 0) = −cf 0 (x) + cg 0 (x) = ψ(x).

Differentiating the first equation then gives with the second equation:
 
0 1 0 1
g (x) = φ (x) + ψ(x)
2 c
so
1 x
Z
1
g(x) = (φ(x) − φ(0)) + ψ(y) dy
2 2c 0
1 x
Z
1
f (x) = (φ(x) + φ(0)) − ψ(y) dy
2 2c 0
and we obtain the very elegant d’Alembert’s solution of the wave equation:

u(x, t) = f (x − ct) + g(x + ct)

1 x+ct
Z
1
= [φ(x + ct) + φ(x − ct)] + ψ(y) dy. (9)
2 2c x−ct

We see that u(x, t) is determined fully by the values of the initial functions φ, ψ in the
interval [x − ct, x + ct] of the x-axis, whose endpoints are cut out by the characteristics
through the point (x, t). This interval is called the domain of dependence for the
solution at (x, t).
Conversely the initial data at the point (ξ, 0) of the x-axis at time t = 0 influences u(x, t)
at points (x, t) in the wedge-shaped region bounded by the characteristics x ± ct = ξ
through (ξ, 0) i.e. the region ξ − ct < x < ξ + ct. Thus disturbances or signals travel
only with speed c.
In particular, discontinuities in the initial data φ propagate along characteristics. For
example consider φ(x) = H(x) and ψ(x) = 0. Then u(x, t) = 21 [H(x − ct) + H(x + ct)]
which is an initial (t = 0) unit step discontinuity at x = 0 simply propagating to the left
and right at speeds ±c, each with half unit heights.
1B Methods 108

10 Green’s functions for PDEs

In this final chapter we will apply the idea of Green’s functions to PDEs, enabling us to
solve the wave equation, diffusion equation and Laplace equation in unbounded domains,
and also with forcing terms (i.e. inhomogeneous PDEs). In some of these developments,
the Fourier transform will play a key role (via the “differentiation becomes multiplica-
tion” rule) and to each of our equations we will associate a fundamental FT pair and a
corresponding so-called fundamental solution.

10.1 FTs for the diffusion equation

Consider the Cauchy problem for the 1D diffusion equation

∂θ ∂ 2θ
=D 2 θ(x, 0) = f (x)
∂t ∂x
and θ → 0 as |x| → ∞ (e.g. diffusion of heat in an infinitely long bar with initial
temperature f (x)). Taking FTs of the diffusion equation w.r.t. x and writing the FT of
θ(x, t) as θ̃(k, t) we have

∂
θ̃(k, t) = −Dk 2 θ̃(k, t) θ̃(k, 0) = f˜(k)
∂t
and a simple t integration (using the initial condition) gives
2
θ̃(k, t) = f˜(k) e−Dtk . (10)

To apply the convolution theorem (to invert the RHS product in eq. (10)) we need to
identify the function g(x, t) such that
2
g̃(k, t) = e−Dtk

i.e. we want the inverse FT of a Gaussian. On exercise sheet 3 question 9 you’ve already
derived the FT pair for a general Gaussian viz.
√
−a2 x2 π − k22
φ(x) = e ↔ φ̃(k) = e 4a .
a
√
Hence identifying 1/4a2 as Dt and rescaling by π/a we get the fundamental Fourier
transform pair for the diffusion equation
1 −x2 2
g(x, t) = √ e 4Dt ↔ g̃(k, t) = e−Dk t , (11)
4πDt
and so, using the convolution theorem, the general solution to our Cauchy problem is
Z ∞ Z ∞
1 (x−u)2
− 4Dt
θ(x, t) = √ f (u) e du = f (u) Sd (x − u, t) du, (12)
4πDt −∞ −∞
1B Methods 109

where the function Sd (x − u, t) is called the fundamental solution or the source func-
tion for the diffusion equation. Here it is associated with the unforced (homogeneous)
equation with inhomogeneous ICs, but below we will see that it arises also in the solution
of the forced (inhomogeneous) equation with homogeneous ICs.
The Gaussian pulse as IC
Suppose the IC for the heat equation is the Gaussian
r
a −ax2
f (x) = θ0 e
π
normalised so that Z ∞
f (x) dx = θ0 .
−∞

Substituting into eq. (12) we obtain

Z ∞
θ0 a1/2 (x − u)2
 
2
θ(x, t) = √ exp −au − du
4π 2 Dt −∞ 4Dt
Z ∞
θ0 a1/2 −([1 + 4aDt]u2 − 2xu + x2 )
 
= √ exp du
4π 2 Dt −∞ 4Dt
h i
ax2
θ0 a1/2 exp − 1+4aDt
= √
4π 2 Dt
"
Z ∞  2 #
−(1 + 4aDt) x
× exp u− du,
−∞ 4Dt 1 + 4aDt
r Z ∞
ax2
 
1 4aDt 2
= θ0 exp − √ e−v dv
1 + 4aDt 4π 2 Dt 1 + 4aDt −∞
where we have used the substitution
r   r
1 + 4aDt x 4Dt
v= u− du = dv.
4Dt 1 + 4aDt 1 + 4aDt
Thus finally
ax2
r  
a
θ(x, t) = θ0 exp − .
π(1 + 4aDt) 1 + 4aDt
Thus an initial Gaussian retains a Gaussian form, with its squared width (1 + 4aDt)/a
spreading linearly with t (recall linear growth of variance for diffusing probabilistic pro-
cesses) while the total area remains constant (cf conservation law of diffusing material
that is built in to the diffusion equation) and the peak at x = 0 drops as t−1/2 . These
features are illustrated in the figure.
The delta function pulse as IC
If the IC for the diffusion equation is the delta function

f (x) = θ0 δ(x)
1B Methods 110

0.5

0.45

0.4

0.35

0.3

0.25

0.2

0.15

0.1

0.05

0
−10 −8 −6 −4 −2 0 2 4 6 8 10

Figure 1: Plots (with a = D = 1 = θ0 ) of the Gaussian pulse solution for: t = 0.1 (solid
line); t = 1 (dashed); t = 10 (dotted); and t = 100 (dot-dashed).
1B Methods 111

then substitution into eq. (12) and the sampling property of δ(x) gives
Z ∞
−(x − u)2
 
θ0
θ = √ δ(u) exp du (13)
4πDt −∞ 4Dt
 2
θ0 −x
= √ exp = θ0 Sd (x, t), (14)
4πDt 4Dt
i.e. the solution is the source function itself. Thus for all t > 0 the δ-pulse spreads as a
Gaussian and as t → 0 from above, we regain the δ function as a Gaussian in the limit
of zero width while keeping the area constant (and hence unbounded height).
x2
Note also the ubiquitous appearance of the dimensionless group η 2 = 4Dt
that we previ-
ously introduced (in chapter 4) as the so-called similarity variable.

10.2 The forced heat equation

Consider the forced 1D heat equation on an infinite domain, with homogeneous BCs
∂ ∂2
θf (x, t) − D 2 θf (x, t) = f (x, t) θf (x, 0) = 0, (15)
∂t ∂x
(with subscript ‘f’ for ‘forced’).
Remark. In the previous section we solved the homogeneous equation with inhomo-
geneous BCs and here we are considering the forced (inhomogeneous) equation with
homogeneous BCs. For the general case of the forced equation with inhomogeneous BCs
we can reduce the problem to these two cases, writing the solution as y = yh + yf where
yh satisfies the homogeneous equation with the given inhomogeneous BCs and yf satisfies
the forced equation with homogeneous BCs. This decomposition will clearly apply to our
other equations as well. 
To solve the problem of eq. (15) we seek a Green’s function G(x, t; ξ, τ ) such that
∂ ∂2
G(x, t; ξ, τ ) − D 2 G(x, t; ξ, τ ) = δ(x − ξ)δ(t − τ ), G(x, 0; ξ, τ ) = 0, (16)
∂t ∂x
because then (formally at least)
Z ∞ Z ∞
θf (x, t) = G(x, t; ξ, τ )f (ξ, τ ) dξ dτ,
0 −∞

when substituted into the heat equation (using a double application of the sampling
property of delta functions) will satisfy eq. (15).
Taking FTs of eq. (16) w.r.t. x (and recalling that the FT of δ(x − ξ) is e−ikξ ) we obtain
2
(after multiplying through by eDk t )
∂ h Dk2 t i 2
e G̃(k, t; ξ, τ ) = e−ikξ+Dk t δ(t − τ ), G̃(k, 0; ξ, τ ) = 0,
∂t Z t
2 2
eDk t G̃(k, t; ξ, τ ) = e−ikξ eDk u δ(u − τ ) du.
0
1B Methods 112

2τ
so by the sampling property, the remaining integral is zero if t < τ , and is equal to eDk
if t > τ . Thus
2
G̃(k, t; ξ, τ ) = H(t − τ )e−ikξ e−Dk (t−τ )
H(t − τ ) ∞ ik(x−ξ) −Dk2 (t−τ )
Z
G(x, t; ξ, τ ) = e e dk
2π −∞
H(t0 ) ∞ ikx0 −Dk2 t0
Z
= e e dk
2π −∞

where we have written x0 = x − ξ and t0 = t − τ . But now we have just recovered

the Gaussian function written in the FT pair eq. (11) (with variables x0 , t0 ) and so we
immediately get
h i
−(x−ξ)2
H(t − τ ) exp 4D(t−τ )
G(x, t; ξ, τ ) = p (17)
4πD(t − τ )
= H(t − τ )Sd (x − ξ, t − τ ). (18)

Note that this involves the same source function

x2
 
1
Sd (x, t) = √ exp −
4πDt 4Dt

that we encountered in the solution of the homogeneous equation with inhomogeneous

ICs viz. eq.(12).
Duhamel’s principle
Duhamel’s principle asserts a relationship between (a) the solution of a forced equation
with homogeneous ICs, and (b) the solution of the unforced (homogeneous) equation
with non-homogeneous ICs. Let us recall our solution for (b) as given in eq. (12) with
a slight generalisation – let us take our initial data at time t = τ (instead of t = 0) and
solve for t > τ . A simple time translation of eq. (12) gives the solution of

θt − Dθxx = 0 with θ(x, τ ) = g(x)

to be Z ∞
θ(x, t) = g(u)Sd (x − u, t − τ ) du t>τ (19)
−∞

i.e. the initial data g(x) at t = τ has been propagated for a time t − τ up to time t.
Now consider our solution of the forced problem obtained above viz. the problem

θt − Dθxx = f (x, t) with θ(x, 0) = 0

having solution
Z t Z ∞ 
θf (x, t) = f (u, τ ) Sd (x − u, t − τ ) du dτ. (20)
0 −∞
1B Methods 113

If we view the forcing term f (x, t) for each fixed time t = τ as an IC(!) at initial time
t = τ then the integral in the square brackets represents its effect propagated to time
t via eq. (19). Then the time integral in eq. (20) expresses the forced solution θf as
the accumulation (superposition) of all these IC effects from all times τ earlier than t,
propagated for time interval t − τ up to time t. The upper limit t of this integral is a
causality effect (expressed in the Green’s function by the H(t − τ ) factor) – an IC applied
at a time τ cannot influence the past t < τ , but contributes to the solution only for all
later times t > τ .

10.3 The forced wave equation

Consider the forced 1D wave equation with homogeneous ICs:

∂ 2 yf 2
2 ∂ yf
− c = f (x, t) (21)
∂t2 ∂x2
∂yf
yf (x, 0) = 0, (x, 0) = 0
∂t
with −∞ < x < ∞ and t ≥ 0. As before we seek a Green’s function G(x, t; ξ, τ ) satisfying

∂2 2 ∂
2
G(x, t; ξ, τ ) − c G(x, t; ξ, τ ) = δ(x − ξ)δ(t − τ ),
∂t2 ∂x2
∂
G(x, 0; ξ, τ ) = 0, G(x, 0; ξ, τ ) = 0 (22)
∂t
so then Z ∞ Z ∞
yf (x, t) = f (ξ, τ )G(x, t; ξ, τ ) dξ dτ (23)
0 −∞

will be a solution of eq. (21).

Taking FTs with respect to x (and wlog taking c > 0) we get

∂2
G̃(k, t; ξ, τ ) + k 2 c2 G̃(k, t; ξ, τ ) = e−ikξ δ(t − τ ).
∂t2
with both G̃ and ∂ G̃/∂t being zero at t = 0.
Now we recognise this as an ODE Green’s function problem for an IVP, essentially the
same as the example given previously at the end of section 7.3 – the latter example just
needs to be slightly generalised by introducing the constant k 2 c2 and the scale factor
e−ikξ to multiply δ(t − τ ) on RHS. Thus the solution is:

e−ikξ sin[kc(t − τ )]H(t − τ )

G̃(k, t; ξ, τ ) = .
kc
To invert this FT we use Dirichlet’s discontinuous formula (in the last line below) to get:

H(t − τ ) ∞ eik(x−ξ) sin[kc(t − τ )]

Z
G(x, t; ξ, τ ) = dk,
2πc −∞ k
1B Methods 114

H(t − τ ) ∞ cos[k(x − ξ)] sin[kc(t − τ )]

Z
= dk,
πc 0 k
H(t − τ ) ∞ sin[k(x − ξ + c[t − τ ])]
Z
= dk
2πc 0 k
H(t − τ ) ∞ sin[k(x − ξ − c[t − τ ])]
Z
− dk
2πc 0 k
H(t − τ )
= { sgn(x − ξ + c[t − τ ]) − sgn(x − ξ − c[t − τ ]) } .
4c
Note that the H(t − τ ) factor imposes t − τ > 0 and we can simplify the above formula
using the following
Exercise: Suppose that B > 0. Then

sgn(A + B) − sgn(A − B) = 2H(B − |A|).

(Just consider the cases A > B, A < −B and −B < A < B separately.) 
Thus we get the causal fundamental solution of the wave equation

H(c(t − τ ) − |x − ξ|)
G(x, t; ξ, τ ) = . (24)
2c
(Note that this expression already requires t − τ > 0 for G to be non-zero so we need not
include the H(t − τ ) factor, which is redundant.)
We see that G is non-zero only when |x − ξ| < c(t − τ ) so a forcing disturbance f (ξ, τ )
at position ξ and time τ can affect a point at x only for times t > τ + |x − ξ|/c i.e. we
have a a finite speed c of propagation of disturbances.
Using our Green’s function, the solution (21) is now
Z tZ x+c(t−τ )
1
yf (x, t) = f (ξ, τ ) dξ dτ. (25)
2c 0 x−c(t−τ )

Notice that the order of integration in this case is important, correctly capturing the
τ -dependent domain of influence of the forcing i.e. the ξ limits depending on τ .
Exercise. (Duhamel’s principle for the wave equation).
Recalling D’Alembert’s solution eq. (9) of the unforced (homogeneous) wave equation
with non-homogeneous ICs, note that
Z x+c(t−τ )
1
I(x, t) = f (ξ, τ ) dξ t>τ
2c x−c(t−τ )

is the D’Alembert solution with initial data given at time τ : u(x, τ ) = 0 ut (x, τ ) = f (x, τ ).
Hence the solution eq. (25) of the forced equation with homogeneous initial data can
again be viewed as a superposition of influences from the forcing term f (x, τ ) used as an
IC (here an IC for the time derivative) for each τ < t.
1B Methods 115

10.4 Poisson’s equation

Poisson’s equation is the forced Laplace equation:

∇2 u = −f (x) (26)
and we will consider it in domains D (possibly infinite) in 2 and 3 space dimensions x
with suitable BCs (cf. later). We will approach its solution by establishing a notion
of Green’s function for it. Unlike the diffusion and wave equations, we have no time
co-ordinate and correspondingly no causality constraint. We will also encounter other
differences, making its solution via Green’s functions less straightforward than that of
the latter two equations. We’ll use a generalisation of the delta function to more than one
dimension. The natural generalisation, denoted δ(r − r0 ) is deemed to have the following
properties:
Z 
1 r0 ∈ D
δ(r − r0 ) = 0, ∀r 6= r0 , δ(r − r0 )dr = (27)
D 0 otherwise,
and the sampling property
Z
f (r) δ(r − r0 ) dV = f (r0 ) if r0 ∈ D.
D

The integration is over either a part of R2 or R3 and then dV denotes respectively the
surface of volume element.
The free-space Green’s function
The fundamental solution to Poisson’s equation is defined to be the solution to the
problem
∇2 G(r; r0 ) = δ(r − r0 ).
Since the problem rotationally symmetric (in 2D or 3D) about the special point r0 , the
fundamental solution can only depend on the scalar distance from that point:
G(r; r0 ) = G(|r − r0 |) = G(r).
Let us write G3 and G2 for this function in 3D and 2D respectively.
Integrating over a sphere (in 3D) or circle (in 2D) of radius r centred at r0 we obtain
Z Z
2
∇ G3 (r) dV = 1 = ∇G3 (r).n̂ dS
V S3
Z I
2
∇ G2 (r) dS = 1 = ∇G2 (r).n̂ dl
S C

where we have used the divergence theorem and its 2D analogue (usually called Green’s
theorem). Here S3 is the surface of the sphere and C is the circumference of the circle.
In both cases ∇G.n̂ = dG
dr
and we get
dG3 G2
4πr2 =1 2πr =1
dr dr
1B Methods 116

and so
1
G3 (r) = − + c3 (28)
4πr
log r
G2 (r) = + c2 (29)
2π
and r = |r − r0 |. In 3D we often apply the far field BC G3 (r) → 0 as r → ∞ setting
c3 = 0 to get the so-called free-space Green’s function:
1
G3 (r) = − .
4π|r − r0 |
Hereafter, G3 will denote this function. In 2D we cannot apply an analogous far field
condition. Note that these fundamental solutions satisfy the Laplace equation for all
r 6= r0 but at r0 they are singular so we need to exercise special care (cf. next section)
in using them as Green’s functions for the Poisson equation throughout a domain.
Green’s identities
Green’s identities (sometimes called Green’s theorems, but not to be confused with the
2D version of the divergence theorem!) establish a relationship between the fundamental
solutions G2 and G3 and solutions of Poisson’s equation. We will discuss here only the
3D case. The 2D case is similar (exercise and cf. also exercise sheet IV).
For two scalar functions φ and ψ in some volume V with surface S and outward normal
n̂, the divergence theorem and product rule gives us that
Z Z Z
2
∇.(φ∇ψ) dV = [φ∇ ψ + (∇φ).(∇ψ)] dV = φ∇ψ.n̂ dS,
V V S

which is Green’s first identity. By simply interchanging the roles of φ and ψ we also
have
Z Z
ψ∇φ.n̂ dS = [ψ∇2 φ + (∇φ).(∇ψ)] dV,
S V

and subtracting this from the previous form we get Green’s second identity:
Z   Z
∂ψ ∂φ
φ −ψ dS = [φ∇2 ψ − ψ∇2 φ] dV (30)
S ∂n ∂n V

(where we have written ∂φ/∂n for ∇φ.n̂ etc.)

Now we’ll want to use ψ = G3 in Green’s second identity but G3 is singular at the point
r0 in V so it is not a priori clear that eq. (30) will still hold (since the divergence theorem
is usually stated as requiring functions that are regular throughout V ). However below
we will show that in fact, Green’s second identity does hold with ψ = G3 . Then using
∇2 G3 = δ(r − r0 ) the first volume integral on RHS becomed φ(r0 ) and we obtain
Z Z  
2 ∂G3 (r; r0 ) ∂φ(r)
φ(r0 ) = G3 (r; r0 ) ∇ φ dV + φ(r) − G3 (r; r0 ) dS (31)
V S ∂n ∂n
1B Methods 117

(the integrations being over the space of variable r).

—————————————————————–
To see that eq. (31) is actually valid, for simplicity let us choose the point r0 to be the
origin. Then we define a new volume, V , which is V with a little ball B centred at
the origin, of radius  (and with surface S ), cut out of V . Now, in V , G3 is perfectly
well-behaved, so we can legitimately apply Green’s second identity in V . Furthermore
we have ∇2 G3 = 0. Therefore eq. (30) gives
Z Z
2 2
[φ∇ G3 − G3 ∇ φ] dV = − G3 ∇2 φ dV
V V
Z   Z  
∂G3 ∂φ ∂G3 ∂φ
= φ − G3 dS + φ − G3 dS. (32)
S ∂n ∂n S ∂n ∂n
Now consider the behaviour of the surface integral over S (second integral above), as 
becomes arbitrarily small. Since G3 = −1/(4π) on S the second term of the S integral
is
4π2
Z Z
∂φ 1 ∂φ
− G3 dS = dS = A
S ∂n 4π S ∂n 4π
∂φ
where A is the average value of ∂n on S (i.e. the integral over the surface divided by
∂φ
the total surface area). Since φ is regular at the origin, A → ∂n |0 (a finite value) as
 → 0, and so the above integral tends to zero as  → 0.
On the other hand, for the first term of the S integral, the r derivative acts on G3 =
−1/(4πr). Thus (remembering that the outward normal from V on the surface S points
in the negative r direction) we get
Z Z
∂G3 1
φ dS = − φ dS = −φ
S ∂n 4π2 S

where φ is the average value of φ over S , which tends to φ(0) as  → 0. Hence

Z
∂G3
φ dS → −φ(0) as  → 0
S ∂n
Putting all this together, (and wlog removing the reliance on the special point r0 being
at the origin) we get precisely eq. (31) as required.
—————————————————————-
Now returning to Poisson’s equation ∇2 φ = −f we substitute φ = u into eq. (31) to get
Green’s third identity:
Z
u(r0 ) = (−f (r))G3 (r; r0 ) dV
V
Z  
∂G3 (r; r0 ) ∂u
+ u − G3 (r; r0 ) dS. (33)
S ∂n ∂n
This is a remarkable formula, as it describes the solution throughout the interior of the
domain in terms of of the properties of the solution on the boundary (and the free-
space Green’s function) Also notice that (unlike the previous cases we have considered
1B Methods 118

in the course) the Green’s function is here providing an expression for the solution with
inhomogeneous BCs.
Exercise: Green’s third identity in 2D (cf. also exercise sheet IV).
Show that the equivalent result in two dimensions in a domain S with perimeter C (with
arc length element dl) is
Z I  
∂G2 ∂u
u(r0 ) = (−f )G2 dS + u − G2 dl,
S C ∂n ∂n
where G2 as defined by eq. (29) with c2 = 0. 
We make two remarks about eq. (33).
∂u
Remark. In the infinite domain of full 3D space, if u and ∂n tend to zero asymptotically
suitably fast, then eq. (33) shows that G3 is the free space Green’s function for the
Poisson
R equation. This also follows formally from ∇2 G3 = δ(r − r0 ), since if u(r0 ) =
V
−f (r)G3 (r − r0 ) dV then (taking ∇2 w.r.t. variable r0 ):
Z Z
2 2
∇ u(r0 ) = −f (r)∇ G3 dV = −f (r)δ(r − r0 ) dV = −f (r0 ).
V V

Remark. By setting f = 0, eq. (33) also provides an expression for the solution u of the
∂u
Laplace equation in the interior of a domain, in terms of the values of u and ∂n on the
boundary. But if D is say, a closed bounded domain, then specification of u alone on the
boundary already determines u (satisfying Laplace’s equation) uniquely (the Dirichlet
∂u
problem) and the normal derivative ∂n on the boundary cannot be freely independently
∂u
specified. Similarly (Neumann problem) specifying ∂n on the boundary determines u up
to an additive constant, so u cannot be freely specified on the boundary. Nevertheless
eq. (33) is a valid expression for any function that is harmonic (i.e. a solution of the
∂u
Laplace equation) in terms of its values of u and ∂n on the boundary. As such it is not
a useful expression to solve for u in D (since the surface data cannot be freely specified)
but nevertheless it is a useful mathematical relation for proving further properties of
harmonic functions.
In view of this curious “over-determined” property of eq. (33) we will now finally ask
how can the Dirichlet problem, for the Laplace and Poisson equations, be solved using
Green’s function techniques, in domains with boundaries. One may also ask about the
Neumann problem (where the normal derivative is specified on the boundary) but in
this course we will discuss only the Dirichlet problem (in 3D here, and a 2D example is
given on exercise sheet IV). One brief remark though, about the Neumann problem –
it is somewhat more complicated than the Dirichlet problem since (unlike the Dirichlet
case) a consistency condition must be satisfied by the boundary data if a solution is to
∂u
exist at all. Indeed suppose ∂n = h(r) is specified on the boundary ∂D. Then by the
divergence theorem we must have
Z Z Z Z
2
h dS = ∇u.n dS = ∇ u dV = − f dV
∂D ∂D D D

(where f is the RHS forcing function in the Poisson equation) and this consistency
condition then complicates the construction of a Green’s function.
1B Methods 119

10.5 Dirichlet Green’s functions

The Dirichlet Green’s function for the Laplacian operator on some domain D (con-
taining both r and r0 ) is defined be the function G(r; r0 ) such that:
(GR1): G(r; r0 ) = G3 (r; r0 ) + H(r, r0 ) where H is finite throughout D (including at
r = r0 ) and H satisfies the Laplace equation throughout D; and
(GR2): G(r; r0 ) = 0 on the boundary of D i.e. G is G3 modified by a harmonic function
H chosen to make G zero on the boundary.
Note that G also satisfies the Laplace equation at all r 6= r0 .
Assuming for the moment that we can find such a Green’s function, we are then able
to find the solution to Poisson’s equation on the domain D with Dirichlet boundary
conditions, as follows.
Let ∇2 u = −f in D with u = h(r) given on ∂D. Substituting G3 = G − H into Green’s
third identity eq. (33) gives
Z  
∂(G − H)
Z
∂u
u(r0 ) = u − (G − H) dS + (−f )(G − H)dV. (34)
∂D ∂n ∂n D

Now H is harmonic throughout D so Green’s second identity (with u and H) gives

Z   Z
∂H ∂u
u −H dS = − (−f )HdV.
∂D ∂n ∂n D

Hence all the H terms in eq. (34) cancel and we get

Z   Z
∂G(r; r0 ) ∂u(r)
u(r0 ) = u(r) − G(r; r0 ) dS + [−f (r)]G(r; r0 )dV. (35)
∂D ∂n ∂n D

Finally having chosen H so that G is zero on ∂D (and u = h(r) on ∂D) we get

Z Z
∂G(r; r0 )
u(r0 ) = h(r) dS + [−f (r)]G(r; r0 )dV. (36)
∂D ∂n D

This expression is now constructive, in the sense that the solution throughout the domain
is given in terms of the known boundary conditions, and the Green’s function.
Remark: if instead we had a Neumann problem with ∂u/∂n = k(r) being specified on
the boundary ∂D then we would instead seek H in (GR1) to make ∂G/∂n = 0 on ∂D
in (GR2) (so H is defined up to an additive constant). Then from eq. (35) we get the
solution to the Neumann problem as
Z Z
u(r0 ) = −k(r)G(r; r0 ) dS + [−f (r)]G(r; r0 ) dV.
∂D D

Exercise: Symmetry of the Green’s function

Using Green’s second identity we can prove (cf. exercise sheet IV) that the Green’s func-
tion is always symmetric, i.e. G(r; r0 ) = G(r0 ; r), for all r 6= r0 . This is the mathematical
statement of the principle of reciprocity in electrostatics; a source at x has the same
effect at x0 as a source at x0 would have at x.
1B Methods 120

10.6 Method of images

So how do we find the Green’s function (G above, i.e a suitable H) in a domain with
boundaries? In general this is difficult (i.e. we need to solve the Dirichlet problem for
the Laplace equation with non-trivial BCs to get H!), but for domains with sufficient
symmetry we can sometimes use the method of images (also called the reflection
method) to construct the required Green’s function. Indeed, this method can also be
used for the Green’s functions of the forced heat and wave equations too, as well as for
homogeneous equations – the key concept is match the boundary conditions. We will
illustrate it with some examples (cf. exercise sheet IV for more examples).
Example: Laplace equation – Dirichlet Green’s function for the half-space
Consider a domain D = {(x, y, z) : z > 0}. We want to find the solution to the following
problem defined in D:

∇2 u = 0, u(x) → 0 as |x| → ∞, u(x, y, 0) = h(x, y).

To use the formula eq. (36), we need to construct a Green’s function which satisfies the
conditions (GR1) and (GR2), interpreting the zero boundary condition in the far field
in the natural way of requiring G → 0 as |x| → ∞. Since the problem has a natural
Cartesian geometry, let us write r = x = [x, y, z], and r0 = x+
0 = [x0 , y0 , z0 ].

We know that the free space Green’s function G3 satisfies all the conditions except the
homogeneous BC G = 0 boundary z = 0. We need to “cancel” its nonzero values there
while retaining all our other desired properties. To achieve this, the elegant method of
images proceeds as follows:
Start from the free space Green’s function and then try to add some other solution of
Laplace’s equation to it to get the required boundary condition. To do this, imagine that
there is an image of the special point outside of the domain, exactly the same vertical
distance away from the boundary as the special point i.e. we reflect the special point in
the boundary as a mirror. To cancel the values of G3 on the boundary we take the image
as a point source with opposite sign, and postulate that
−1 1
G(x; x0 ) = + + ,
4π|x − x0 | 4π|x − x− 0|
1  −1/2
= − (x − x0 )2 + (y − y0 )2 + (z − z0 )2
4π
1  −1/2
+ (x − x0 )2 + (y − y0 )2 + (z + z0 )2 ,
4π
where x− −
0 = [x0 , y0 , −z0 ]. Now since x0 is definitely out of the domain, the second term
satisfies Laplace’s equation everywhere in D, so we have (GR1). H and G3 are both
asymptotically zero at the far boundary |x| → ∞, and finally, for any point xb on the
boundary z = 0, the two terms cancel perfectly; so we have (GR2). Therefore we have
found the Dirichlet Green’s function!
To apply the formula (36), note also that f = 0 (for Laplace’s equation here) and also
that there is no contribution from the far field since u → 0 is the far field BC of the
1B Methods 121

problem. The outward normal from the domain at z = 0 is in the negative z-direction,
and so the only contribution to the expression comes from the lower boundary, and is
 
∂G  ∂G 
= −
∂n z=0 ∂z z=0
 
1 z + z0 z − z0
= −
4π |x − x− 0|
3 |x − x+0|
3
z=0
z0  2 2 2 −3/2

= (x − x0 ) + (y − y0 ) + z0 .
2π
Therefore we get
z0 ∞ ∞ 
Z Z
−3/2
u(x0 , y0 , z0 ) = (x − x0 )2 + (y − y0 )2 + z02 h(x, y) dx dy,
2π −∞ −∞
as the solution of the Dirichlet problem for the Laplace equation in the half space.
The method of images can also be applied to some more complicated domains – e.g. for
the domain a ≤ z ≤ b between two planar boundaries, we can view the two boundaries
as mirrors and then any special point r0 in the domain gives rise to a double infinity of
multiply reflected images. If these are taken with alternating plus and minus signs, we
can construct a Green’s function that vanishes on both boundaries.
Example: Images for wave problems
We can also apply the method of images to the Green’s functions that we developed
for the forced wave and heat equations (with ICs all being zero). For example consider
the semi-infinite domain 0 ≤ x < ∞ and suppose we wish to impose the homogeneous
Dirichlet condition G(0, t; ξ, τ ) = 0 on the causal Green’s function. Exactly the same
idea applies as above: add an equal amplitude Green’s function with opposite sign and
having ‘special point’ not at ξ but at −ξ. For example, the appropriate Dirichlet Green’s
function for the wave equation is
H(c[t − τ ] − |x − ξ|) H(c[t − τ ] − |x + ξ|)
G(x, t; ξ, τ ) = − ,
2c 2c
i.e. an odd function, which automatically imposes the zero Dirichlet boundary condition.
Similarly for a homogeneous Neumann BC at x = 0 viz.

∂G(x, t; ξ, τ ) 
=0 for all t,
∂x 
x=0
we make the derivative an odd function in x, which we achieve by an even extension of
the Green’s function. For example, the appropriate Neumann Green’s function for the
wave equation is
H(c[t − τ ] − |x − ξ|) H(c[t − τ ] − |x + ξ|)
G(x, t; ξ, τ ) = + ,
2c 2c
and the image has the same sign. Indeed, for sufficiently small, yet still positive x,
|x − ξ| = ξ − x, and |x + ξ| = ξ + x, so far all t

∂G  1
= (δ[c(t − τ ) − |x − ξ|][1] + δ[c(t − τ ) − |x + ξ|][−1])x=0
∂x 
x=0 2c
= 0 for all t
1B Methods 122

as desired.
Example: images for homogeneous problems
For homogeneous diffusion or wave equation problems we do not use a Green’s function,
but the method of images can still be applied to the initial data itself to adapt our previous
solutions to be applicable in the presence of extra boundaries. For example consider the
wave equation for y(x, t) on the half-line x ≥ 0 with Neumann BC yx (0, t) = 0 at the
boundary x = 0, and with initial conditions

y(x, 0) = b(x) yt (x, 0) = 0 for x ≥ 0

where b(x) is the ‘top-hat’ function having value 1 in [x0 − a, x0 + a] and zero elsewhere
(and a < x0 ) i.e. a box of width 2a, and height 1 centred on x0 .
According to d’Alembert’s general solution eq. (9), the solution for the whole line will be
a pair of boxes, each of height half and width 2a, one going to the left and one going to
the right, with speed c. This solution clearly does not satisfy the Neumann BC at x = 0
(when the left moving box passes over the origin). Thus we consider an initial image box
(on the negative x axis) by reflecting the given ICs in x = 0 to extend the problem to
the whole line (and take the image with a plus sign, for the Neumann condition).
Applying d’Alembert’s solution we now have two pairs of moving boxes (of heights half).
The solution satisfies the wave equation everywhere, satisfies the given ICs for x ≥ 0 and
satisfies the Neumann BC at x = 0 (as the solution is always an even function of x for
any fixed t. Thus taking this solution restricted to just x ≥ 0, we have solved our original
problem!
It is instructive to picture how this solution evolves in time (especially in the region
x ≥ 0). Let BR and BL denote the right and left moving boxes for the original IC, and
IR and IL the boxes from the image IC. IL just moves off to the left and never appears as
part of the solution for x ≥ 0. BR just moves off to the right and is always present in the
x ≥ 0 region. At time t = (x0 − a)/c both BL and IR reach x = 0 (with their respective
leading edges). As they pass through each other (for a time 2a/c) the solution piles up
to double height where they overlap and then BL goes off into x < 0 while IR emerges
fully into x ≥ 0. Thereafter, the solution in x ≥ 0 is just two boxes of heights half, with
centres separated by constant distance 2x0 (the initial separation of IL and BR ) moving
off to the right at speed c.
If we view the solution near x = 0 just on the positive side, we see BL “hit the wall”,
piling up to double height (and enclosing the same total area at all times) as it appears to
be reflected back to the right into x > 0 with final shape unchanged, giving the trailing
second right-moving box.
The Neumann condition at x = 0 physically models small amplitude waves in the vicinity
of a frictionless wall e.g. an elastic string for x ≥ 0 (with zero gravity) that is attached
to a vertical pole at x = 0 by a frictionless light (i.e. massless) ring, so the string can
sustain no spatial gradient at x = 0 i.e. cannot have any vertical force component, and
is thus always horizontal there.
1B Methods 123

==== End of Notes! ====