Partial Differential Equations

Download as pdf or txt
Download as pdf or txt
You are on page 1of 64

MAT-51316

Partial Differential Equations


Robert Pich e
Tampere University of Technology
2010
Contents
1 PDE Generalities, Transport Equation, Method of Characteristics 1
1.1 PDE Generalities . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Transport Equation . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Method of Characteristics . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Example: u
t
+ 2u
x
= 0 . . . . . . . . . . . . . . . . . . . . . . . 3
1.5 Example: u
t
+xu
x
= 0 . . . . . . . . . . . . . . . . . . . . . . . 4
1.6 Example: u
t
+ (xu)
x
= 0 . . . . . . . . . . . . . . . . . . . . . . 5
2 Models of Vibration, Diffusion and Heat Conduction; Ill-Posed Prob-
lems 6
2.1 Vibrating String . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 One-dimensional Diffusion . . . . . . . . . . . . . . . . . . . . . 7
2.3 One-dimensional Heat Conduction . . . . . . . . . . . . . . . . . 8
2.4 Ill-posed Problems . . . . . . . . . . . . . . . . . . . . . . . . . 9
3 One Dimensional Wave Equation 10
3.1 General Solution . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.2 Solution of initial value problem . . . . . . . . . . . . . . . . . . 11
3.2.1 Example: Three-nger pluck . . . . . . . . . . . . . . . . 11
3.2.2 Example: Hammer blow . . . . . . . . . . . . . . . . . . 12
3.3 Energy and Uniqueness . . . . . . . . . . . . . . . . . . . . . . . 13
4 One Dimensional Diffusion Equation 14
4.1 Solution of IVP . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
4.2 General Properties of the IVP Solution . . . . . . . . . . . . . . . 15
4.3 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
4.4 Derivation of solution formula using Fourier transforms . . . . . . 17
5 Duhamels Principle; Half-Line Models 19
5.1 Duhamels Principle for ODE . . . . . . . . . . . . . . . . . . . . 19
5.2 Duhamels Principle for Diffusion Equation . . . . . . . . . . . . 19
5.3 Diffusion/Heat on the Half Line: Reection Method . . . . . . . . 20
i
5.4 Diffusion/Heat on the Half Line with Sources . . . . . . . . . . . 21
5.5 Waves on the Half Line . . . . . . . . . . . . . . . . . . . . . . . 22
6 Separation of Variables 24
6.1 Separation of Variables for the Heat Equation . . . . . . . . . . . 24
6.2 Sturm-Liouville Theory . . . . . . . . . . . . . . . . . . . . . . . 25
6.3 Heat IBVP with Constant Coefcients . . . . . . . . . . . . . . . 26
6.4 Wave IBVP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
7 Numerical Solution of PDEs with Matlab 28
7.1 Specifying an IBVP . . . . . . . . . . . . . . . . . . . . . . . . . 28
7.2 Solving an IBVP . . . . . . . . . . . . . . . . . . . . . . . . . . 29
8 Fourier Series 32
8.1 Least-Squares Approximation, Completeness . . . . . . . . . . . 32
8.2 Classical Fourier series . . . . . . . . . . . . . . . . . . . . . . . 34
8.3 Pointwise Convergence . . . . . . . . . . . . . . . . . . . . . . . 35
8.4 Solving PDE initial value problems . . . . . . . . . . . . . . . . . 36
8.5 Heat equation with source term . . . . . . . . . . . . . . . . . . . 38
9 Laplaces Equation 40
9.1 Some Facts from Vector Analysis . . . . . . . . . . . . . . . . . . 40
9.2 Heat Flow in Three Dimensions . . . . . . . . . . . . . . . . . . 40
9.3 Membrane Vibration . . . . . . . . . . . . . . . . . . . . . . . . 41
9.4 Laplaces Equation . . . . . . . . . . . . . . . . . . . . . . . . . 43
10 Solving Two-Dimensional Laplace Equations 48
10.1 Dirichlet Problem in a disk . . . . . . . . . . . . . . . . . . . . . 48
10.2 Dirichlet Problem in a rectangle . . . . . . . . . . . . . . . . . . 50
10.3 Dirichlet-Neumann Problem in a wedge . . . . . . . . . . . . . . 52
10.4 Dirichlet Problem in the region outside a circle . . . . . . . . . . 53
11 Greens Functions 55
11.1 Greens Function for One-Dimensional Equation . . . . . . . . . 55
11.2 Greens Function for Two-Dimensional Poisson Equation . . . . . 57
11.3 Greens functions from eigenfunctions . . . . . . . . . . . . . . . 60
ii
1 PDEGeneralities, Transport Equation, Method of
Characteristics
how to classify PDEs
how to model one dimensional transport phenomena by a rst-order PDE
how to solve initial value problems for this equation using the method of
characteristics
how to compute and plot solutions using Maple function PDEplot
1.1 PDE Generalities
Recall that an ordinary differential equation (ODE) relates a one-variable function
u(x) and its derivatives in an equation of the form
F(x, u, u

, u

, . . . , u
(n)
) = 0.
The order of the ODE is the highest derivative order that appears in the equation.
For example, the Malthus population growth model
u

(t) = ru(t)
is a rst-order ODE with independent variable t (time), dependent variable u (pop-
ulation), and constant parameter r (net growth rate).
An ODE is said to be linear if it has the form
a
0
(x)u(x) +a
1
(x)u

(x) + a
n
(x)u
(n)
(x)
. .
Lu
= g(x)
and a linear ODE is said to be homogeneous if g 0. Linear homogeneous ODEs
have the superposition property: if u
1
and u
2
are solutions then so is the function

1
u
1
+
2
u
2
for any constants
1
,
2
. For example, the Malthus model is a linear
homogeneous ODE.
A solution of an ODE is a function that satises the equation everywhere in
some domain of the dependent variable. General solutions of ODEs generally
contain arbitrary constants. For example, u(t) = Ae
rt
for any constant A is a
solution of the Malthus model ODE.
Apartial differential equation (PDE) relates a multivariable function u(x, y, . . .)
and its derivatives u
x
=
u
x
, u
xy
=

2
u
xy
, . . . in an equation of the form
F(x, y, . . . , u, u
x
, u
y
, . . . , u
xx
, u
xy
, u
yy
, . . .) = 0
The order of the ODE is the highest derivative order that appears. Linear and
homogenous PDEs are dened analogously to ODEs. Here are some examples of
1
two-variable PDEs that are used to model physical phenomena:
1. u
x
+u
y
= 0 (transport; order 1, linear homogeneous)
2. u
x
+yu
y
= 0 (transport; order 1, linear homogenous)
3. u
x
+uu
y
= 0 (shock wave; order 1, nonlinear)
4. u
xx
+u
yy
= 0 (Laplace eqn; order 2, linear homogeneous)
5. u
tt
u
xx
+u
3
= 0 (wave with interaction; order 2, nonlinear)
6. u
t
+uu
x
+u
xxx
= 0 (dispersive wave; order 3, nonlinear)
7. u
tt
+u
xxxx
= 0 (vibrating beam; order 4, linear homog.)
8. u
t
iu
xx
= 0 (quantum mechanics; order 2, linear homog.)
A solution of a PDE is a function that satises the equation everywhere in some
domain of the dependent variables. For example, both u
1
(x, t) = x
2
+ 2t and
u
2
(x, t) = e
t
sin x are solutions of the PDE u
t
u
xx
= 0 for all (x, t). General
solutions of PDEs generally involve arbitrary functions. For example, the general
solution of u
x
= t sin x is u(x, t) = t cos(x) + (t), and the general solution of
u
xy
= 0 is u(x, y) = F(y) +G(x).
In this course we will see how PDEs arise as mathematical models of phe-
nomena, we will present general properties of solutions, and learn some solution
techniques.
1.2 Transport Equation
Consider a substance (e.g. mass or energy) owing in a region of space. Let
u(x, t) denote its density (units: [quantity] [volume]
1
) as a function of position
x and time t, and let (x, t) denote the ux (units: [quantity] [time]
1
[area]
1
).
(Density and ux variations in the y and z directions are assumed to be negligible.)
The amount of substance in an interval a x b of a tube-shaped region of
constant cross section A is
_
b
a
u(x, t)Adx.
A
b
a x
The net ux into the interval is (a, t)A (b, t)A. Let f(x, t, u) denote the
source term, that is, the rate (units: [quantity] [time]
1
[volume]
1
) at which
substance density increases by processes other than ux, for example chemical
reaction. The rate of increase of the total amount of substance in the interval is
then
d
dt
_
b
a
u(x, t)Adx = (a, t)A (b, t)A +
_
b
a
f(x, t, u)Adx,
which can be rearranged to give
_
b
a
(u
t
+
x
f) dx = 0.
Because [a, b] is arbitrary, this implies that the conservation equation
u
t
+
x
= f
2
should hold at every point in the region.
If we know the velocity c(x, t) (units: [length] [time]
1
) then the ux is
= cu. Substituting this constitutive equation into the conservation equation
gives the transport equation
u
t
+ (cu)
x
= f. (1)
In an initial value problem for the transport equation, one seeks the function
u(x, t) that satises (1) and that satises u(x, 0) = u
0
(x) for some given initial
density prole u
0
.
1.3 Method of Characteristics
The transport equation (1) can be written as cu
x
+ u
t
= f c
x
u, that is, as
c u = g where g(x, t, u) = f c
x
u, c =
_
c
1
_
, and u =
_
u
x
u
t
_
. The
transport equation thus has a geometric interpretation: we seek a surface z =
u(x, t) whose directional derivative in the direction of vector c is g(x, t, u). This
geometric interpretation is the basis for the following solution method.
t
x
k
c
Curves x = X(t) in the (x, t) plane that are tangential
to the vector eld
_
c(x, t)
1
_
are called characteristic curves.
From this denition it follows that the characteristic curve that
goes through the point (x, t) = (k, 0) is the graph of the func-
tion X that satises the ODE
dX
dt
= c(X, t) (2)
with initial condition X(0) = k.
Denoting the value of u along a characteristic curve by U(t) = u(X(t), t), we
have
d
dt
U =
u
x
dX
dt
+
u
t
= cu
x
+u
t
= g,
that is, the value of u along the characteristic curve is determined by the ODE
U

= g(X(t), t, U(t)). (3)


The solution of the ODE (2) with initial value U(0) = u
0
(k) determines the value
of u along the characteristic curve that intersects the x-axis at (k, 0), because
U(0) = u(X(0), 0) = u(k, 0) = u
0
(k). The solution surface is the collection (or
envelope) of space curves created as k takes on all real values.
The Maple code PDEplot produces the graph of the solution using numerical
algorithms to solve the ODEs (23) with initial condition
_
X(0)
U(0)
_
=
_
k
u
0
(k)
_
.
In simple enough cases, the ODEs can also be solved analytically by hand.
1.4 Example: u
t
+ 2u
x
= 0
This equation models transport with constant velocity c(x, t) = 2 and no source
term.
3
t
x
k
x 2t k = The characteristic ODE is X

= 2. The solution of
this ODE satisfying the initial condition X(0) = k is the
straight line X = 2t + k. The characteristic curve (in
this case: the line) through a given point (x, t) crosses
the x axis at (k, 0) with k = x 2t.
The ODE describing the value of u along a charac-
teristic line is U

(t) = 0, i.e. the value is constant along the line. The solution
of this ODE satisfying the initial condition U(0) = u
0
(k) is U(t) = u
0
(k). The
solution of the PDE initial value problem is therefore u(x, t) = u
0
(x 2t). In
particular, if u
0
(x) = e
x
2
then the solution is u(x, t) = e
(x2t)
2
.
2
1.5
1
t
0.5
-2
0
2
4
0.2
0
6
x
0.4
u(x,t)
0.6
0.8
The solution of the PDE u
t
+2u
x
= 0 with ini-
tial prole u
0
(x) = e
x
2
can be plotted in Maple
by the commands
> PDE:=diff(u(x,t),t)+2
*
diff(u(x,t),x)=0;
> with(PDEtools):
> PDEplot(PDE,[x,0,exp(-x2)],
x=-3..3,t=0..2);
The plot shows how the initial prole translates
to the right at constant speed without changing
shape.
1.5 Example: u
t
+xu
x
= 0
This equation can also be written as u
t
+ (xu)
x
= u, which is of the form of
the transport equation (1) with source term f(x, t, u) = u. This equation models
transport in a velocity eld c(x, t) = x, that is, the velocity is equal to the distance
from the origin. The source term f(x, t, u) = u models generation of substance
at a rate equal to the amount of substance.
t
x
k
x = ke
t
The characteristic ODE is X

= X. The solution
of this ODE satisfying the initial condition X(0) = k
is X = ke
t
. The characteristic curve through a given
point (x, t) crosses the x axis at (k, 0) with k = xe
t
.
The ODE describing the value of u along a charac-
teristic curve is U

= 0, i.e. the value is constant along


the curve. The solution of this ODE satisfying the ini-
tial condition U(0) = u
0
(k) is U(t) = u
0
(k). The
solution of the PDE initial value problem is therefore
u(x, t) = u
0
(xe
t
). In particular, if the initial prole is u
0
(x) = e
(x3)
2
then the
solution is u(x, y) = e
(xe
t
3)
2
.
2
1.5
t
1
0.5
0
40
30
20
10 x
0
0.2
0.4
u(x,t)
0.6
0.8
The solution of the PDE u
t
+ xu
x
= 0 with
initial prole u
0
(x) = e
(x3)
2
can be plotted in
Maple by the commands
> PDE:=diff(u(x,t),t)
+x
*
diff(u(x,t),x)=0;
> PDEplot(PDE,[x,0,exp(-(x-3)2)],
x=0..6,t=0..2);
4
The PDE solution spreads out as time advances,
and the surface height remains constant along the
characteristic curves, and so the total amount of
substance increases as time advances.
1.6 Example: u
t
+ (xu)
x
= 0
This equation models transport in the same velocity eld as the previous example,
but without the source term. The characteristic curves are the same as in the
previous example.
Rewriting the equation in the form u
t
+ xu
x
= u, we see that the ODE
describing the value of u along a characteristic curve is U

= U.The solution of
this ODE satisfying the initial condition U(0) = u
0
(k) is U(t) = u
0
(k)e
t
. The
solution of the PDE initial value problem is therefore u(x, t) = u
0
(xe
t
)e
t
. In
particular, if the initial prole is u
0
(x) = e
(x3)
2
then the solution is u(x, y) =
e
(xe
t
3)
2
t
.
2
1.5
t
1
0.5
0
40
30
20
x
10
0
0.2
0.4
u(x,t)
0.6
0.8
The solution of the PDE u
t
+ xu
x
= 0 with
initial prole u
0
(x) = e
(x3)
2
can be plotted in
Maple by the commands
> PDE:=diff(u(x,t),t)
+diff(x
*
u(x,t),x)=0;
> PDEplot(PDE,[x,0,exp(-(x-3)2)],
x=0..6,t=0..2);
The PDE solution spreads out as time advances,
and because there is no source term, the solu-
tion also decreases in amplitude, so that the total
amount of substance remains constant (conserva-
tion law).
5
2 Models of Vibration, Diffusion and Heat Conduc-
tion; Ill-Posed Problems
how to derive the PDE for the vibrating string
how to derive the PDE for one-dimensional diffusion or heat conduction
how to model boundary conditions for these PDEs
some examples of ill-posed problems
2.1 Vibrating String
Consider a thin exible string moving in the xz plane. Assume that points of the
string move in the z direction only, and let z = u(x, t) denote the shape of the
string.
Longitudinal force balance: Let T(x, t) be the ten-
sion, assumed to act tangentially along the string. Let
= tan
1
u
x
denote the angle between the string tan-
gent and the x axis. The only longitudinal forces acting
on the part of the string between x = a and x = b are the
x components of the tension force, and because there is
no longitudinal motion these forces must be equal, that
is,
T(b, t) cos (b, t) = T(a, t) cos (a, t).
Because the segment is arbitrary, this implies T cos is constant with respect to
x, say
T(x, t) cos (x, t) = (t).
Mass conservation: Let (x, t) be the strings mass per unit length, which may
vary as the string deforms during the motion, and let
0
(x) be the mass per unit
length when the string is straight. If dx represents an element of length when the
string is straight and dx

=
_
1 +u
2
x
dx is the length element of the deformed
string, then mass conservation requires that dx

=
0
dx.
Transverse force balance: By Newtons law, the net transversal force on a
string segment [a, b] is equal to the time derivative of the momentum:
d
dt
_
b
a
u
t
dx

= T(b, t) sin (b, t) T(a, t) sin (a, t)


= T(b, t) cos (b, t) tan (b, t)
T(a, t) cos (a, t) tan (a, t)
= (t) (u
x
(b, t) u
x
(a, t))
=
_
b
a
u
xx
dx.
Because of mass conservation, the left hand side can be written as
d
dt
_
b
a
u
t

0
dx =
_
b
a

0
u
tt
dx, and so the transverse force balance becomes
_
b
a
(
0
u
tt
u
xx
) dx = 0.
6
Because the interval is arbitrary, this implies

0
u
tt
u
xx
= 0
for all (x, t) in the solution domain. Denoting c =
_
/
0
, this can be written
u
tt
= c
2
u
xx
,
which is the one-dimensional wave equation.
If the string is modelled to have nite length, say x [0, l], then it is necessary
to specify the boundary conditions. If the motion at the ends is known, this can
be modelled by the Dirichlet condition
u(0, t) = h
0
(t), u(l, t) = h
1
(t)
In particular, xed ends are modelled by h
0
0 and h
1
0.
An end support at x = l that is not perfectly rigid can be modelled as a lin-
ear spring where the spring force is ku(l, t), with k being the spring constant.
T
k
The force balance with the transverse component of the
tension
T(l, t) sin (l, t) = T(l, t) cos (l, t)
. .
(t)
tan (l, t)
. .
u
x
(l, t)
.
produces the Robin condition
(t)u
x
(l, t) + ku(l, t) = 0.
Similarly, a exible support at the other end can be modelled by the Robin condi-
tion
(t)u
x
(0, t) +ku(0, t) = 0.
2.2 One-dimensional Diffusion
Recall from1.2 that the one-dimensional conservation equation relating the den-
sity u(x, t) of material moving with ux (x, t) and source term f(x, t, u) is
u
t
+
x
= f.
Diffusion processes, whereby substance ows from areas of high concentration to
areas of low concentration, can be modelled by the constitutive relation (Ficks
law)
= ku
x
,
where k(x) is a material parameter (diffusivity, units: [length]
2
[time]
1
). Com-
bining these equations gives the one-dimensional diffusion equation
u
t
(ku
x
)
x
= f.
7
g
If the domain is of nite length, say x [0, l], then
it is necessary to specify boundary conditions. For ex-
ample, consider a tube whose end x = l is covered by
a thin permeable membrane, beyond which there is a
large well-stirred reservoir with given density g(t). Sup-
posing the ux through the membrane is proportional to the difference in densities
on its two faces, we have
(l, t) = (u(l, t) g(t)),
where is the membrane permeability. Substituting Ficks law gives the Robin
condition
ku
x
(l, t) +u(l, t) = g(t).
In the limiting case /k 0 (impermeable membrane, i.e. the tube end is closed)
this becomes the Neumann condition
u
x
= 0,
while in the limiting case /k (no membrane) we get the Dirichlet condition
u(l, t) = g(t).
2.3 One-dimensional Heat Conduction
The one-dimensional conservation equation is equally valid when u = cT de-
notes density of heat energy. Here c(x) and (x) are material parameters (specic
heat and mass per unit length) and T is the temperature. The constitutive relation
(Fouriers law)
= KT
x
models conduction, whereby heat ows from hot areas to colder areas. The ma-
terial parameter K(x) is called the heat conductivity. Combining the equations
gives the one-dimensional heat equation
c(T)
t
(KT
x
)
x
= f.
This is very similar to the diffusion equation, and is essentially identical to it when
c is constant.
Similarly to the diffusion equation, one can model a thin insulating layer be-
tween the end x = l and a region with given temperature T
1
by a Robin condition
KT
x
(l, t) +T(l, t) = T
1
(t)
with Neumann and Dirichlet conditions obtained as limiting cases corresponding
to perfect insulation and perfect conduction.
8
2.4 Ill-posed Problems
A mathematical model or problem often consists of a set of differential and alge-
braic equations. However, not all such sets of equations are useful models: the set
should have a unique solution, and the solution should be continuously dependent
on the available data. Such models are said to be well posed problems. Here are
some examples of PDE problems that, although they may appear to be all right
at rst glance, are in fact ill-posed. Well-posed PDE problems will be presented
later in the course.
1. The boundary value problem (BVP)
u

(x) = 0 (0 < x < 1), u

(0) = 0, u

(1) = 1
has no solutions. The problem is overdetermined.
2. The BVP
u

(x) = 0 (0 < x < 1), u

(0) = 0, u

(1) = 0
has innitely many solutions, namely, solutions of the form u = constant.
The problem is underdetermined.
3. The two-dimensional Laplace equation
u
xx
+u
yy
= 0
on the domain y > 0 with boundary conditions u(x, 0) = 0 and u
y
(x, 0) =
0 has the solution u 0. If the second boundary condition is changed to
u
y
(x) = e

n
sin(nx), the solution becomes u(x, y) =
1
n
e

n
sin(nx) sinh(ny).
We can choose n large enough to make max
x
|u
y
(x)| as small as we like,
but no matter how small the perturbation, the solution always blows up as
y . Thus, the problem is unstable: the solution does not depend con-
tinuously on the boundary data.
9
3 One Dimensional Wave Equation
general solution of one dimensional wave equation
dAlemberts solution of initial value problem
uniqueness of IVP solution via energy
3.1 General Solution
The one dimensional wave equation
u
tt
c
2
u
xx
= 0 (x R, t > 0)
models the transverse vibration of a long string whose ends are so far away that
they can be neglected. The PDE can be written as the sytem of rst-order PDEs
_

t
c

x
_
v = 0,
_

t
+c

x
_
u = v.
The PDE v
t
cv
x
= 0 has the characteristic equation X

= c, its characteristic
curves are X = ct + k, and the general solution is v(x, t) = h(x + ct) where h
is an arbitrary function.
The PDE u
t
+ cu
x
= v has the characteristic equation X

= c, and its charac-


teristic curves are X = ct + l. The value of u along the characteristic is U(t) =
u(ct +l, t); it satises the ODE U

(t) = v(ct +l, t) = h(ct +l +ct) = h(2ct +l).


Making the change of variables s = 2ct +l and

U(2ct +l) = U(t), we obtain the
ODE 2c

(s) = h(s), which has the solution



U(s) = f(s) + g where f =
1
2c
_
h
and g is constant along the characteristic. Then
u(x, t) = U(t) =

U(2ct +l) = f(2ct +l) +g(l) = f(2ct + (x ct)) +g(x ct)
= f(x +ct) +g(x ct). (4)
The general solution (4) is the sum of a shape f that moves left at speed c and a
shape g that moves right with speed c, as shown here:
t
x
l
x ct l =
t
x
k
x +ct k =
g
x
f
x
10
x
t
domain of dependence
range of influence
(x
0
,t
0
)
The general solution (4) indicates that information
(about the local shape of the string) propagates at a nite
speed c along the characteristics. The displacement at a given
point in time and space (x
0
, t
0
) can be deduced from values
lying in a cone-shaped domain of dependence of previous val-
ues; values outside this domain have no inuence on the value
of u(x
0
, t
0
). Similarly, any point (x
0
, t
0
) has a cone-shaped range of inuence.
3.2 Solution of initial value problem
The initial value problem for the one dimensional wave equation is to nd the
solution given the initial displacement and velocity, that is, we are given the initial
conditions
u(x, 0) = (x), u
t
(x, 0) = (x).
Substituting t = 0 into u(x, t) = f(x + ct) + g(x ct) and u
t
(x, t) = cf

(x +
ct) cg

(x ct) gives
(x) = f(x) +g(x), (x) = cf

(x) cg

(x).
Differentiating the rst equation and solving gives
f

2
+

2c
, g

2


2c
,
which can be integrated to give
f(x) =
1
2
(x) +
1
2c
_
x
0
() d +A
g(x) =
1
2
(x)
1
2c
_
x
0
() d +B.
_
(5)
Adding these together gives f(x) +g(x) = (x) +A+B, so we have A+B = 0.
Then
u(x, t) = f(x +ct) +g(x ct)
=
1
2
(x +ct) +
1
2c
_
x+ct
0
() d +A +
1
2
(x ct)
1
2c
_
xct
0
() d +B
=
1
2
[(x +ct) + (x ct)] +
1
2c
_
x+ct
xct
() d (6)
Formula (6) is known as dAlemberts formula.
3.2.1 Example: Three-nger pluck
Consider the innite-length string with initial displacement given by the hat
function
(x) =
_
1 |x| |x| 1
0 otherwise
11
and initial velocity 0. dAlemberts formula (6) gives the solution as
u(x, t) =
1
2
[(x +ct) +(x ct)].
This could be written more explicitly using a lot of if clauses, but such a formula
would be difcult for a human reader to interpret. A more geometric approach is
to decompose the initial shape according to (5) with A = B = 0, which gives the
left-moving shape f(x) =
1
2
(x) and the right-moving shape g(x) =
1
2
(x). The
solution is thus the sum of two hat functions, one moving to the left and the other
moving to the right, both at speed c (see Figure 1).
x
u
ct = 0
-1 1
x
u
ct = 0.5
-1.5 1.5
x
u
ct = 2
-3
1
-1 3
x
-1 1
Figure 1: Three-nger pluck solution snapshots
3.2.2 Example: Hammer blow
Consider the innite-length string with zero initial displacement and initial veloc-
ity given by the step function
(x) =
_
1 |x| 1
0 otherwise.
x
-1 1
dAlemberts formula (6) gives the solution
u(x, t) =
1
2c
[(x +ct) (x ct)]
where is the ramp function
(x) =
_
x
0
() d =
_
_
_
1 x < 1
x |x| 1
1 x > 1.
x
-1 1
The solution is the sum of two mirror-image ramp functions, one moving to the
left and the other moving to the right, both at speed c.
12
3.3 Energy and Uniqueness
Dene the total energy E(t) of an innitely long string as the sum of kinetic
energy and potential energy
E =
1
2
_

0
u
2
t
dx +
1
2
_

u
2
x
dx.
Differentiating, we have
E

=
1
2
_

0
(2u
t
u
tt
) dx +
1
2
_

(2u
x
u
xt
) dx
=
_

(u
t
u
xx
+u
x
u
xt
)
. .
(u
t
u
x
)
x
dx =

u
t
u
x
.
If we assume that the support of and (i.e. the subset of the domain where they
are nonzero) is contained in a nite interval [a, b], it follows that u(x, t) is zero for
x outside the range of inuence of [a, b]. Then E

is identically zero and the total


energy E remains constant in time.
Similarly, the total energy E(t) of a nite-length string is dened as
E =
1
2
_
L
0

0
u
2
t
dx +
1
2
_
L
0
u
2
x
dx.
If the string is assumed to be clamped at both ends, that is, u is assumed to be
subject to the homogeneous Dirichlet boundary conditions
u(0, t) = 0, u(L, t) = 0,
then u
t
(0, t) = 0 and u
t
(L, t) = 0, so that E

L
0
u
t
u
x
= 0, and the total
energy remains constant.
The above principle of conservation of energy can be used to prove the unique-
ness of the solution of the initial value problem for the nite-length string with
Dirichlet boundary conditions. Let u and v be solutions of the initial-boundary
value problem, that is,
_

_
u
tt
= c
2
u
xx
u(x, 0) = (x)
u
t
(x, 0) = (x)
u(0, t) = h
0
(t)
u(L, t) = h
1
(t)
_

_
and
_

_
v
tt
= c
2
v
xx
v(x, 0) = (x)
v
t
(x, 0) = (x)
v(0, t) = h
0
(t)
v(L, t) = h
1
(t)
_

_
Then the difference w = u v satises the initial-boundary value problem
_

_
w
tt
= c
2
w
xx
w(x, 0) = 0
w
t
(x, 0) = 0
w(0, t) = 0
w(L, t) = 0
_

_
The energy of w at any time is equal to its energy at t = 0, which is zero, and
so w
x
0, which implies that w is constant with respect to x. To satisfy the
boundary conditions, the constant must be zero. Thus w 0, that is, u v.
13
4 One Dimensional Diffusion Equation
Formula for the solution of the initial value problem
General properties of the solution
Example IVPs
Derivation of the formula using Fourier Transforms
4.1 Solution of IVP
The solution of the initial value problem for the diffusion equation u
t
ku
xx
= 0
with u(x, 0) = (x) is given by formula (8) below.
Theorem 1 Let be a bounded and piecewise continuous function on R, and let
S(x, t) =
1

4kt
e

x
2
4kt
(x R, t > 0), (7)
where k is a positive constant. Then
u(x, t) =
_

S(x , t)() d (8)


is a smooth function that satises u
t
= ku
xx
for x R, t > 0, and
lim
t0
u(x, t) =
(x
+
) +(x

)
2
(x R). (9)
PROOF. First, note that S > 0 and that
_

S(x, t) dx =
..
p =
x

4kt
1

e
p
2
dp = 1. (10)
Then
|u(x, t)|
_
max
xR
|(x)|
_
. .
M
_

S(x , t) d = M, (11)
and so the improper integral (8) converges. Also,
u
x
(x, t) =
_

S
x
(x , t)() d
=
1

4kt
_

x
2kt
e

(x)
2
4kt
() d [ p = (x )/

4kt
=
1

kt
_

pe
p
2
(x p

4kt) dp,
and this integral converges because
|u(x, t)|
M

kt
_

|p|e
p
2
dp =
M

kt
.
14
Similarly, it can be shown that u
t
, u
xx
, u
xt
, and derivatives of higher orders all
exist, so u is smooth. It satises the diffusion equation because
(u
t
ku
xx
)(x, t)) =
_

(S
t
kS
xx
. .
=0
)(x , t)() d = 0.
The proof of (9) is omitted.
4.2 General Properties of the IVP Solution
1. The formula (8) is a convolution
1
and can be written u = S (with t
treated as a parameter).
2. The function S is variously called the kernel, source function, or fundamen-
tal solution of the diffusion (or heat) equation. We have lim
t0
S(x, t) = 0
for x = 0 and lim
t0
S(0, t) = . Snapshots of S at various time instants
show an initial narrow peak that spreads out as time advances.
!8 !4 0 4 8
0
0.4
0.8
x
S
kt = 0.1
kt = 1
kt = 10
3. Any jump discontinuities in the initial shape or in its derivatives are in-
stantly smoothed out not like the wave equation.
4. If > 0 on a nite interval [a, b] and is zero elsewhere, we have u(x, t) > 0
for all x (no matter how large) and all t > 0 (no matter how small). Thus, in
this model, information has innite speed of propagation not like the
wave equation.
5. A small change in the initial condition produces a small change in the solu-
tion. That is, if
u
t
= ku
xx
(x R, t > 0)
u(x, 0) = (x)
and
v
t
= kv
xx
(x R, t > 0)
v(x, 0) = (x),
then the difference w = u v satises w
t
= kw
xx
with initial condition
w(x, 0) = (x) (x), and by (11) we have |w(x, t)| max
xR
|(x)
(x)|.
6. The initial value problem has at most one solution. This can be proved by
setting = in the previous argument.
1
The convolution of two functions f and g is denoted f g and is given by (f g)(x) =
_

f(x)g() d. Convolution is semilinear (i.e. (f) g = (f g)), commutative (f g =


gf), associative (f (gh) = (f g)h), and distributes over addition (f (g+h) = f g+f h).
15
7. The identity (10) can be interpreted in terms of the IVP: u 1 is indeed a
solution of the diffusion equation for 1.
4.3 Examples
Initial prole = step Solve the one-dimensional diffusion equation u
t
= ku
xx
with initial condition
u(x, 0) = Heaviside(x) =
_
0 (x < 0)
1 (x > 0)
.
x
1
Write the solution using the standard function erf, which is dened as
erf(u) =
2

_
u
0
e
p
2
dp.
1
!1
!3 3
x
Solution. The solution formula (8) gives
u(x, t) = ( S)(x, t)
=
1

4kt
_

e
z
2
/(4kt)
(x z) dz
=
1

4kt
_
x

e
z
2
/(4kt)
dz [ p = z/

4kt
=
1

_
x/

4kt

e
p
2
dp
=
1

_
0

e
p
2
dp +
1

_
x/

4kt
0
e
p
2
dp
=
1
2
+
1
2
erf(x/

4kt).
!8 !4 0 4 8
0
0.5
1
x
u
kt = 0.1
kt = 1
kt = 10
Initial prole = exponential Solve the one-dimensional diffusion equation u
t
=
ku
xx
with initial condition
u(x, 0) = (x) = e
x
.
16
Solution. The solution formula (8) gives
u(x, t) =
1

4kt
_

e
(x)
2
/(4kt)
e

d
=
1

4kt
_

(x2kt)
2
4kt
+ktx
d
= e
(xkt)
1

4kt
_

(x2kt)
2
4kt
d
. .
=1
= (x kt),
that is, the initial shape translates to the right with speed k.
4.4 Derivation of solution formula using Fourier transforms
The Fourier transform F = Ff and the inverse Fourier transform f = F
1
F are
given by the formulas
F() =
_

f(x)e
ix
dx, f(x) =
1
2
_

F()e
ix
d.
These transforms have many useful properties, including:
Linearity of FT and IFT: If g(x) =
1
f
1
(x)+
2
f
2
(x) then G() =
1
F
1
()+

2
F
2
(). If G() =
1
F
1
() +
2
F
2
() then g(x) =
1
f
1
(x) +
2
f
2
(x).
FT of derivative: (Ff

)() = iF().
FT of convolution: (F[f
1
f
2
])() = F
1
()F
2
().
The Fourier transform can be used to solve the diffusion equation IVP as
follows. Transforming u
t
ku
xx
gives the ordinary differential equation U
t
+
k
2
U = 0, which has the solution U(, t) = Ce

2
kt
. The initial condition gives
C = U(, 0) = (), so we have U(, t) = ()e

2
kt
. Thus u = f , where
f is the inverse Fourier transform of F() = e

2
kt
.
The inverse Fourier transform can be found from standard tables, or it can be
derived as follows. Differentiating f(x) =
1
2
_

2
kt
e
ix
d gives
f

(x) =
1
2
_

2
kt
ie
ix
d
=
1
2

i
2kt
e

2
kt
e
ix

1
2
_

x
2kt
e

2
kt
e
ix
d
=
x
2kt
f(x).
Multiplying by e
x
2
/(4kt)
gives the differential equation
e
x
2
/(4kt)
f

(x) +
x
2kt
e
x
2
/(4kt)
f(x)
. .
_
e
x
2
/(4kt)
f(x)
_

= 0,
17
which has the solution e
x
2
/(4kt)
f(x) = constant. The constant is determined by
the initial condition
f(0) =
1
2
_

2
kt
d =
1

4kt
.
Thus f(x) =
1

4kt
e

x
2
4kt
, which coincides with the formula of the fundamental
solution given in (7).
18
5 Duhamels Principle; Half-Line Models
how to solve linear problems with source terms (Duhamels Principle)
how to solve diffusion models on the half-line x > 0
how to solve wave models on the half-line x > 0
5.1 Duhamels Principle for ODE
Consider the linear homogeneous ODE
u +Au = 0, (12)
where u(t) is a vector and A is a constant square matrix. The source operator
S(t) is a square matrix such that, for any vector , u(t) = S(t) is a solution
of (12) with initial condition u(0) = . It follows that
S(0) = I and

S(t) +AS(t) = 0,
where I is the identity matrix. For example, the scalar IVP u+u = 0, u(0) =
has the solution u(t) = e
t
, that is, the source operator is S(t) = e
t
.
The following theorem shows how the general solution of the homogeneous
ODE can be used to solve the ODE with a source term.
Theorem 2 (Duhamels Principle) The function
u(t) =
_
t
0
S(t )f() d +S(t) (13)
satises the differential equation u +Au = f with initial condition u(0) = .
PROOF. We have u(0) = 0 +S(0) = and
u(t) =
_
t
0

S(t )f() d +S(t t)f(t) +



S(t)
= A
__
t
0
S(t )f() d +S(t)
_
+If(t)
= Au(t) +f(t),
which completes the proof.
Note that the integral in (13) can also be written
_
t
0
S()f(t ) d.
5.2 Duhamels Principle for Diffusion Equation
Consider now the diffusion equation, which can be written as (12) with the dot
representing

t
and A representing the operator k

2
x
2
for functions with spatial
domain x R. The general solution of the one-dimensional diffusion equation
IVP on R was shown earlier to be
u(x, t) =
_

S(x , t)() d, where S(x, t) =


1

4kt
e

x
2
4kt
.
19
From this formula we see that the source operator S(t) for this problem transforms
the function (x) into the function
_

S(x , t)() d. Duhamels principle


thus gives the solution of the IVP with source term
u
t
ku
xx
= f(x, t), u(x, 0) = (x)
as
u(x, t) =
_
t
0
_

S(x , t )f(, ) d d +
_

S(x , t)() d. (14)


Example Solve the diffusion equation u
t
ku
xx
= e
x
on x R, t > 0 with
initial condition u(x, 0) = 0.
SOLUTION. Duhamels formula gives
u(x, t) =
_
t
0
_

1
_
4k(t )
e

(x)
2
4k(t)
e

d d
=
_
t
0
e
x+k(t)
d =
1
k
(e
kt
1)e
x
,
which can be veried (by substitution) to satisfy the equation and initial condition.
5.3 Diffusion/Heat on the Half Line: Reection Method
Consider the initial-boundary value problem
v
t
kv
xx
= 0, (x > 0, t > 0)
v(x, 0) = (x) (x > 0)
v(0, t) = 0 (t > 0).
_
_
_
(15)
This models, for example, the ground temperature at depth x (on a at planet!)
given an initial temperature prole and a surface temperature that is xed at
zero.
To solve this initial-boundary value problem, we exploit the fact that the solu-
tion of the diffusion problem on the whole line is odd whenever the initial prole
is odd. We introduce the odd extension of , that is, the function on R given by

odd
(x) =
_
_
_
(x) (x > 0)
(x) (x < 0)
0 (x = 0)
x
The solution of u
t
ku
xx
= 0 on x R with initial condition u(x, 0) =
odd
(x)
is
u(x, t) =
_

S(x , t)
odd
() d
=
_

0
S(x , t)
odd
() d +
_
0

S(x , t)
odd
() d
=
_

0
S(x , t)() d
_
0

S(x , t)() d
=


_

0
S(x +, t)() d
20
The solution of the original initial-boundary value problem (15) is then the restric-
tion of u(x, t) to the half-line x > 0, that is,
v(x, t) =
_

0
S
haline
(x, , t)() d, (16)
where
S
haline
(x, , t) = S(x , t) S(x +, t)
Example Solve the IBVP (15) with initial prole (x) = 1. This models a
sudden drop in the surface temperature.
SOLUTION. The odd extension of the initial prole can be written as

odd
(x) = 2Heaviside(x)1 x
1
-1
and so, using the result from section 4.3, we have
v(x, t) = 2
_
1
2
+
1
2
erf(x/

4kt)
_
1 = erf(x/

4kt).
0 2 4 6 8
0
0.5
1
x
v
kt = 0.1
kt = 1
kt = 10
5.4 Diffusion/Heat on the Half Line with Sources
The solution formula (16) can be written as S
haline
(t), where the source operator
S
haline
(t) transforms the function (x) into the function
_

0
S
haline
(x, , t)() d.
Then, using Duhamels Principle, the solution of the IBVP with source function
w
t
kw
xx
= f, (x > 0, t > 0)
w(x, 0) = (x) (x > 0)
w(0, t) = 0 (t > 0).
_
_
_
(17)
can be written directly as
w(x, t) =
_
t
0
_

0
S
haline
(x, , t )f(, ) d d +
_

0
S
haline
(x, , t)() d.
Now consider the IBVP with Dirichlet boundary condition
y
t
ky
xx
= 0, (x > 0, t > 0)
y(x, 0) = 0 (x > 0)
y(0, t) = h(t) (t > 0).
21
Setting w(x, t) = y(x, t) h(t) yields the IBVP
w
t
kw
xx
=

h, (x > 0, t > 0)
w(x, 0) = h(0) (x > 0)
w(0, t) = 0 (t > 0),
which is of the same form as (17) so, using its solution and the results of the
example in section 5.3, we have
y(x, t) = h(t)
_
t
0

h()
_

0
S
haline
(x, , t ) d d h(0)
_

0
S
haline
(x, , t) d
= h(t)
_
t
0
erf
_
x
_
4k(t )
_

h() d h(0) erf


_
x

4kt
_
.
5.5 Waves on the Half Line
Consider the wave equation initial-boundary value problem
v
tt
c
2
v
xx
= 0 (x > 0, t > 0)
v(x, 0) = (x), v
t
(x, 0) = (x) (x > 0)
v(0, t) = 0 (t > 0).
_
_
_
(18)
This models a semi-innite string with one end xed.
We can use the reection method here too, because the solution of the wave
equation on the whole line is odd whenever both and are odd. Applying
dAlemberts formula, the solution of the odd-extended problem is
u(x, t) =
1
2
[
odd
(x +ct) +
odd
(x ct)] +
1
2c
_
x+ct
xct

odd
() d,
and so the solution of (18) is
v(x, t) =
_
_
_
1
2
[(x +ct) + (x ct)] +
1
2c
_
x+ct
xct
() d (x ct > 0)
1
2
[(x +ct) (ct x)] +
1
2c
_
x+ct
ctx
() d (x ct < 0).
Example The solution of the IBVP (18) with initial prole (x) = e
(x3)
2
and
zero initial velocity is
v(x, t) =
_
_
_
1
2
[e
(x+ct3)
2
+e
(xct3)
2
] (x ct > 0)
1
2
[e
(x+ct3)
2
e
(ctx3)
2
] (x ct < 0).
22
0 2 4 6 8
!1
0
1
x
ct = 0
ct = 0.75
ct = 1.5
ct = 2.25
ct = 3
ct = 3.75
ct = 4.5
ct = 5.25
23
6 Separation of Variables
rst steps in solving the heat and wave equation on an interval
eigenvalues and eigenfunctions (Sturm-Liouville theory)
6.1 Separation of Variables for the Heat Equation
Consider the IBVP
u
t
(u
x
)
x
= 0 (x (0, l), t > 0)
u(0, t) = 0, u(l, t) = 0 (t 0)
u(x, 0) = (x) (x (0, l))
_
_
_
(19)
where (x) = c(x)(x) > 0 and (x) > 0. This models heat ow in a pipe of
length l with xed temperatures at the ends (Dirichlet boundary conditions).
Substituting a trial solution of the form u(x, t) = X(x)T(t) into the PDE
in (19) and rearranging gives
T

T
=
(X

X
.
The function on the left side of this equation is constant with respect to x and the
function on the right side is constant with respect to t, so they are both equal to a
constant, call it . We then have two ODEs, namely
T

+T = 0,
which has general solution T(t) = T(0)e
t
, and
(X

+X = 0, (20)
which, because of the boundary conditions in (19), has the boundary conditions
X(0) = 0, X(l) = 0. (21)
The solution method is as follows:
1. nd numbers
n
and nonzero functions X
n
that satisfy the BVP (2021);
2. express the initial condition as a linear combination of the X
n
, that is,
(x) =

n
C
n
X
n
(x);
3. then u(x, t) =

n
C
n
X
n
(x)e

n
t
.
In the remainder of this lecture we focus on stage 1 of the method.
24
6.2 Sturm-Liouville Theory
Theorem 3 There are innitely many pairs of numbers
n
(eigenvalues) and nonzero
functions X
n
(eigenfunctions) that are solutions of problem (2021). The eigen-
values are real and positive, the eigenfunctions corresponding to distinct eigen-
values are -orthogonal, that is,

m
=
n

_
l
0
(x)X
m
(x)X
n
(x) dx = 0,
and every eigenvalue has multiplicity 1, that is, the corresponding eigenfunction
is unique up to a multiplicative factor.
PROOF. First, note that for any two eigenfunctions we have
(X

m
X
n
X

n
X
m
)

= (
n

m
)X
m
X
n
;
this can be veried by expanding the left hand side then substituting the ODE.
Now, if (
n
, X
n
) satises (2021) then so does the complex conjugate pair (

n
,

X
n
),
and so
(
n

n
)
_
l
0


X
n
X
n
dx =

l
0
(

X

n
X
n
X

X
n
) = 0,
and dividing this through by
_
l
0
|X
n
|
2
dx (which is > 0) leads us to the result

n
= 0, that is, the eigenvalues are real.
Next, if
m
=
n
,
_
l
0
X
m
X
n
dx =
1

l
0
(X

m
X
n
X

n
X
m
) = 0,
and so X
m
and X
n
are -orthogonal.
Next, multiplying the ODE (20) by X and integrating, we obtain
_
l
0
X(X

dx +
_
l
0
X
2
dx = 0,
from which we can solve for and get
=

_
l
0
X(X

dx
_
l
0
X
2
dx
=
_
l
0
(X

)
2
dx

l
0
XX

_
l
0
X
2
dx
0.
The case = 0 can be ruled out, because
= 0
_
l
0
(X

)
2
dx = 0 X

0 X is constant
and the only constant that is consistent with the boundary conditions (21) is zero.
Finally, if (, X
1
) and (, X
2
) satisfy (2021), we have
(X

1
X
2
X

2
X
1
)

= ( )X
1
X
2
= 0,
and so (X

1
X
2
X

2
X
1
) = constant. The boundary condition (21) implies that
the constant is zero, so X

1
X
2
X

2
X
1
= 0. Then (X
2
/X
1
)

= 0, so X
2
/X
1
is a constant, that is, the eigenfunctions corresponding to are identical up to a
multiplicative factor.
The proof of existence and inniteness of number of eigenvalues is omitted.
25
6.3 Heat IBVP with Constant Coefcients
For the heat equation (19) with constant and , equation (20) has the general
solution
X(x) = Acos(x
_
/k) +Bsin(x
_
/k),
where k = /. The boundary conditions (21) imply that A = 0 and that
sin(l
_
/k) = 0.
This is satised when l
_
/k = n for n Z, so we have the eigenvalues

n
= k(n/l)
2
for n {1, 2, . . .}. (The solution = 0 is rejected because the
general solution of X

= 0 is X(x) = E+Fx, and the boundary conditions imply


E = 0 and F = 0.) The corresponding eigenfunctions are X
n
(x) = sin(
nx
l
).
0 1
!1
0
1
X
1

x/l
0 1
!1
0
1
X
2

0 1
!1
0
1
X
3

0 1
!1
0
1
X
4

0 1
!1
0
1
X
5

0 1
!1
0
1
X
6

For initial conditions of the form
(x) = C
1
sin(
x
l
) +C
2
sin(
2x
l
) + +C
n
sin(
nx
l
),
the solution of (19) is a linear combination of spatial sinusoids with amplitudes
that are exponentially decaying in time:
u(x, t) = C
1
sin(
x
l
)e
k
2
t/l
2
+C
2
sin(
2x
l
)e
4k
2
t/l
2
+ +C
n
sin(
nx
l
)e
kn
2

2
t/l
2
.
The solution tends to zero as time advances: all the heat eventually leaks out of
the ends of the tube. Note that the terms corresponding to larger n have wavier
shape and decay in time faster.
6.4 Wave IBVP
Consider the IBVP

0
u
tt
u
xx
= 0 (x (0, l), t > 0)
u(0, t) = 0, u(l, t) = 0 (t 0)
u(x, 0) = (x), u
t
(x, 0) = (x) (x (0, l))
_
_
_
(22)
where
0
(x) > 0 and > 0. This models the small-amplitude transverse motion
of a taught exible string with xed ends.
26
Proceeding as for the heat equation, we assume a trial solution of the form
u(x, t) = X(x)T(t) and obtain two ODEs,
T

+T = 0,
which has the general solution T(t) = T(0) cos(t

) +
1

(0) sin(t

), and
the eigenvalue problem
X

+
0
X = 0, X(0) = 0, X(l) = 0.
This is a special case of (2021), so the results of Theorem 1 apply here also:
there are real positive eigenvalues 0 <
1
<
2
< with unique eigenfunctions
X
1
(x), X
2
(x), . . .. If the initial conditions are linear combinations of eigenfunc-
tions, that is, if
(x) =

n
A
n
X
n
(x) and (x) =

n
B
n
X
n
(x),
then the solution of (22) is a superposition of shapes whose amplitudes vary sinu-
soidally in time:
u(x, t) =

n
_
A
n
cos(t
_

n
) +
1

n
B
n
sin(t
_

n
)
_
X
n
(x).
The factors

n
are called natural frequencies and have units [radians per time
unit].
In the case where
0
is constant, the eigenvalues are
n
= (nc/l)
2
and the
eigenfunctions are X
n
(x) = sin(nx/l) for n {1, 2, . . .}, where c =
_
/
0
.
In this case, all the natural frequencies are integer multiples of the fundamental
frequency

1
= c/l =

l
_

0
. From this formula we can explain various
musical phenomena associated with guitar or violin strings:
the note rises by one octave (i.e. the frequency is doubled) when the string
is clamped at its midpoint, because the clamping produces two vibrating
strings, each half the length;
the note rises when the string is tightened, because the tightening increases
the value of .
27
7 Numerical Solution of PDEs with Matlab
How to solve IBVPs in one spatial dimension using pdepe
7.1 Specifying an IBVP
The Matlab solver pdepe solves PDEs of the form
(x, t, u, u
x
)u
t
= x
m
(x
m
f(x, t, u, u
x
))
x
+s(x, t, u, u
x
) (x (a, b), t (t
0
, t
nal
])
The constant m may be 0, 1 or 2; if m > 0 then a must be non-negative. The ux
term f must depend on u
x
. The term must be non-negative, and may only be
zero at mesh points. Spatial discontinuities in or the source term s are allowed
but only at mesh points.
The problem has an initial condition of the form
u(x, t
0
) = (x) (a x b).
The boundary conditions are
p
left
(a, t, u(a, t)) +q
left
(a, t)f(a, t, u(a, t), u
x
(a, t)) = 0,
p
right
(b, t, u(b, t)) + q
right
(b, t)f(b, t, u(b, t), u
x
(b, t)) = 0
for t t
0
, where q
left
and q
right
are either identically zero or never zero.
Thus, the mathematical problem is completely dened by the specifying the
values m, a, b, t
0
, t
nal
and by the functions , f, s, p
left
, q
left
, p
right
, q
right
, .
Example 1 Consider the PDE
2
u
t
= u
xx
on 0 < x < 1 and 0 < t 2 with
boundary conditions
u(0, t) = 0, u
x
(1, t) = e
t
and initial condition u(x, 0) = sin(x). This models for example the tempera-
ture in a rod with c =
2
and K = 1 that is insulated along its length, its left
end maintained at constant zero temperature, and ux at the right end given by
Ku
x
(1) = e
t
(the negative sign implies that heat ows out of the rod at this
end).
The specication of the problem for solution by pdepe is
m = 0, a = 0, b = 1, t
0
= 0, t
nal
= 2,
(x, t, u, u
x
) =
2
, f(x, t, u, u
x
) = u
x
, s(x, t, u, u
x
) = 0,
p
left
(a, t, u(a, t)) = u(a, t), q
left
(a, t) = 0,
p
right
(b, t, u(b, t)) = e
t
, q
right
(b, t) = 1,
(x) = sin(x).
The exact solution for this problem can be obtained by the method of separation
of variables:
u(x, t) = e
t
sin(x).
28
Example 2 Consider the PDE
u
t
= x
2
(x
2
f(x, t, u, u
x
))
x
+s(x, t, u, u
x
) (0 < x < 1, 0 < t 1),
where
f(x, t, u, u
x
) =
_
5u
x
x (0, 0.5]
u
x
x (0.5, 1)
, s(x, t, u, u
x
) =
_
1000e
u
x (0, 0.5]
e
u
x (0.5, 1)
The boundary conditions are u
x
(0, t) = 0 and u(1, t) = 1, and the initial condition
is
(x) =
_
0 x (0, 1)
1 x = 1.
The specication of the problem for solution by pdepe is
m = 2, a = 0, b = 1, t
0
= 0, t
nal
= 1,
(x, t, u, u
x
) = 1,
p
left
(a, t, u(a, t)) = 0, q
left
(a, t) = 1
p
right
(b, t, u(b, t)) = u(b, t) 1, q
right
(b, t) = 0.
and with f, s, specied as above.
7.2 Solving an IBVP
The syntax of the Matlab PDE solver for a single PDE is
sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)
where
m is 0, 1 or 2,
pdefun is a handle to a function that computes , f and s, with calling syntax
[mu,f,s] = pdefun(x,t,u,ux)
icfun is a handle to a function that computes the initial condition , with calling
syntax
phi = icfun(x)
bcfun is a handle to a function that computes the boundary condition. Its calling
syntax is
[pleft,qleft,pright,qright] = bcfun(a,ua,b,ub,t)
where ua and ub are the values of u(a, t) and u(b, t). For m > 0 and
a = 0 the solver automatically uses the boundary condition u
x
(0, t) = 0
and ignores the values returned in pleft and qleft.
29
xmesh is a vector of points in [a, b] where the solution is approximated. The
solution interval end points a and b are xmesh(1) and xmesh(end),
the values of xmesh must be monotonically increasing, and the length of
xmesh must be at least 3.
tspan is a vector of time values where the solution is approximated. The start
and end times t
0
and t
nal
are tspan(1) and tspan(end), the values of
tspan must be monotonically increasing, and the length of tspan must
be at least 3.
sol is a three-dimensional array where sol(i,j,1) is the solution value at
time tspan(i) and mesh point xmesh(j). (The third dimension is for
systems of PDEs.)
Example 1 can be coded by the functions
function [mu,f,s] = pdex1pde(x,t,u,ux)
mu = pi2;
f = ux;
s = 0;
function phi = pdex1ic(x)
phi = sin(pi
*
x);
function [pleft,qleft,pright,qright] = pdex1bc(a,ua,b,ub,t)
pleft = ua;
qleft = 0;
pright = pi
*
exp(-t);
qright = 1;
and solved using
x = linspace(0,1,20);
t = linspace(0,2,5);
sol = pdepe(0,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
u = sol(:,:,1);
surf(x,t,u)
30
0
0.2
0.4
0.6
0.8
1
0
0.5
1
1.5
2
0
0.2
0.4
0.6
0.8
1
Distance x
Time t
The code to solve this example is pdex1 which you can run from the Matlab
command line. You can look at it using the command
edit pdex1
See pdex2 for the code that solves Example 2.
31
8 Fourier Series
How to approximate a function by a linear combination of orthogonal func-
tions
Fourier series and its convergence in mean square and pointwise
Solving PDEs using Fourier series
8.1 Least-Squares Approximation, Completeness
Consider the approximation of a function f(x) by a linear combination

N
n=1
c
n
X
n
(x)
of functions X
1
, X
2
, . . . that are -orthogonal on (a, b).
Theorem 4 The coefcients c
1
, c
2
, . . . , c
N
that minimize the mean-square error
of the approximation,
E
N
=
_
b
a
(x)
_
f(x)
N

n=1
c
n
X
n
(x)
_
2
dx,
are the Fourier coefcients
c
n
=
_
b
a
(x)f(x)X
n
(x) dx
_
b
a
(x)(X
n
(x))
2
dx
.
PROOF. We use the notation (f, g) =
_
b
a
(x)f(x)X
n
(x) dx (inner product)
and f =
_
(f, f) (2-norm). Then
E
N
= f

n
c
n
X
n

2
= (f

m
c
m
X
m
, f

n
c
n
X
n
)
= f
2
2

n
c
n
(f, X
n
) +

n
c
m
c
n
(X
m
, X
n
)
. .

n
c
2
n
X
n

2
= f
2
+

n
X
n

2
_
c
n

(f, X
n
)
X
n

2
_
2

n
(f, X
n
)
2
X
n

2
,
which is minimized when c
n
= (f, X
n
)/X
n

2
. 2
Theorem 5 (Bessels inequality) The Fourier coefcients of f satisfy

n=1
c
2
n
X
n

2
f
2
. (23)
32
PROOF. Substituting c
n
= (f, X
n
)/X
n

2
into the last line of the previous proof
gives

N
n=1
c
2
n
X
n

2
f
2
. Because the sequence of partial sums is monotone
and bounded, the series converges and the limit satises (23). 2
We saw earlier that the eigenfunctions of a Sturm-Liouville problem corre-
sponding to distinct eigenvalues are -orthogonal. The eigenfunctions corre-
sponding to a single eigenvalue span a space of dimension at most 2, so one can
nd an orthogonal basis of the space spanned by all the eigenfunctions. The fol-
lowing result, whose proof is omitted, tells about the convergence of the best least
squares approximation that uses this basis.
Theorem 6 The set of orthogonal eigenfunctions of a Sturm-Liouville problem is
complete in the sense that for every function with nite 2-norm f, the best least
squares approximation s
N
(x) =

N
n=1
(f, X
n
)X
n
(x)/X
n

2
converges to f in
the mean square sense: lim
N
f s
N
= 0. Also,

n=1
c
2
n
X
n

2
= f
2
(Parsevals identity).
Example 1 The functions sin(nx/l) are the eigenfunctions of X

+X = 0 on
(0, l) with boundary conditions X(0) = 0, X(l) = 0. The corresponding Fourier
coefcients of the function f(x) = 1 (0 < x < l) are
c
n
=
_
l
0
sin(nx/l) dx
_
l
0
sin
2
(nx/l) dx
=
2(1 (1)
n
)
n
,
and the eigenfunction expansion is
f(x) =
4

_
sin(x/l) +
1
3
sin(3x/l) +
1
5
sin(5x/l) +
_
.
Some partial sums s
N
(x) =

N
n=1
c
n
sin(nx/l) are:
0 1
0
1
s
1
x / l
s
3
s
5
s
21
Parsevals identity for this example gives the interesting series

n=1,3,...
1
n
2
=

2
8
.
33
8.2 Classical Fourier series
The trigonometric basis functions 1, cos(x/l), sin(x/l), cos(2x/l), sin(2x/l),. . .
are 2l-periodic and orthogonal on(l, l), that is,
_
l
l
cos(mx/l) cos(nx/l) dx = 0 (m = n)
_
l
l
sin(mx/l) sin(nx/l) dx = 0 (m = n)
_
l
l
cos(mx/l) sin(nx/l) dx = 0
This can be veried using trigonometric identities; the rst two orthogonality
results can also be derived using the orthogonality of eigenfunctions of distinct
eigenvalues for the Sturm-Liouville problem X

+X = 0 with periodic bound-


ary conditions, as in problem 5 of exercise set 6. Then, using the results
_
l
l
dx = 2l,
_
l
l
cos
2
(nx/l) dx =
_
l
l
sin
2
(nx/l) dx = l (n 1),
we nd that the coefcients that minimize the mean square error of
s
N
(x) =
1
2
a
0
+
N

n=1
a
n
cos(nx/l) +b
n
sin(nx/l)
as an approximation of f are
a
n
=
1
l
_
l
l
f(x) cos(nx/l) dx (n 0),
b
n
=
1
l
_
l
l
f(x) sin(nx/l) dx (n 1).
These are the coefcients of the classical Fourier series, which by Theorem 3
converges in the mean square sense to f provided that f is nite. Parsevals
identity can be written
1
2
a
2
0
+

n=1
_
a
2
n
+b
2
n
_
=
1
l
_
l
l
(f(x))
2
dx.
Example 2 The Fourier series coefcients of the function f(x) = x/l (l <
x < l) are
a
n
=
1
l
_
l
l
x
l
cos(nx/l) dx = 0 (n 0),
b
n
=
1
l
_
l
l
x
l
sin(nx/l) dx =
2
n
(1)
n+1
(n 1).
and so the Fourier series for f is
2

_
sin(x/l)
1
2
sin(2x/l) +
1
3
sin(3x/l) . . .
_
.
34
Some partial sums s
N
(x) are:
!1 0 1
!1
0
1
s
1
x / l
s
2
s
4
s
10
Parsevals identity for this example gives the interesting series

n=1
1
n
2
=

2
6
.
The above example illustrates the fact that the Fourier series of an odd-symmetric
function has only sine terms. Similarly, the Fourier series of an even function has
only cosine terms. These facts can be used to relate Fourier series with sine or co-
sine series. For example, the series in Example 1 extended to the interval (l, l)
is the Fourier series of the step function
f(x) =
_
_
_
1 l < x < 0
0 x = 0
1 0 < x < l.
x
-1
1
-l l
8.3 Pointwise Convergence
A function f is said to be piecewise continuous on the nite interval [a, b] if f(a
+
)
and f(b

) exist and f is continuous on (a, b) except for a nite number of simple


jumps (i.e. points of discontinuity where the left and right limits exist). A function
is piecewise continuous on R if it is piecewise continuous on every nite interval.
A function f is said to be piecewise C
1
on [a, b] if it is piecewise continuous on
[a, b], f

(a
+
) and f

(b

) exist, and f

exists and is continuous in (a, b) except for a


nite number of simple jumps. A function is piecewise C
1
on R if it is piecewise
C
1
on every nite interval.
Theorem 7 (Dirichlets Theorem) The Fourier series of a 2l-periodic function
f that is piecewise C
1
on R converges to
f(x
+
)+f(x

)
2
for all x. In particular, it
converges to f(x) at every point x where f is continuous.
The function in Example 2 can be extended to a 2l-periodic function with
f(nl) = 0 for n Z:
!5 !3 !1 1 3 5
!1
0
1
x / l
This extended function is piecewise C
1
, with a jump of 2 at every odd multiple of
35
l. According to Theorem 4, the Fourier series converges to the extended function
at every x. The oscillation near the jump that is seen in the partial sums (Gibbs
phenomenon) does not spoil the convergence because it becomes innitesimally
narrow as the number of terms is increased.
Similarly, the function in Example 1 can be extended to a 2l-periodic odd
function:
!5 !4 !3 !2 !1 0 1 2 3 4 5
!1
0
1
x / l
This extended function is piecewise C
1
, so the Fourier series converges to it at
every x.
0 1 2
0
1
2
3
s
1
s
2
s
3
x
Pointwise convergence does not imply mean square conver-
gence. For example, the function sequence
s
N
(x) =
N
1 +N
2
x
2
(x > 0)
converges pointwise to the zero function, but does not converge
in the mean-square sense because
s
N

2
>
_
l
0
N
2
(1 + N
2
x
2
)
2
dx = N
_
Nl
0
1
(1 + y
2
)
2
dy,
and
_

0
(1 + y
2
)
2
dy = /4, so s
N
.
Also, mean square convergence does not imply pointwise convergence. For
example, the function sequence
f
N
(x) =
_
N |x| < 1/N
3
0 otherwise
(x R)
x
N
-3
-N
-3
N
converges in mean square to the zero function, because
_

(f
N
(x))
2
dx =
_
1/N
3
1/N
3
N
2
dx =
2
N
0,
but does not converge pointwise at x = 0.
8.4 Solving PDE initial value problems
Heat equation with Dirichlet boundary conditions As we saw in lecture 6,
the solution of the IBVP
u
t
ku
xx
= 0 (x (0, l), t > 0)
u(0, t) = 0, u(l, t) = 0 (t 0)
u(x, 0) = (x) (x (0, l))
36
is
u(x, t) =

n=1
c
n
sin
_
nx
l
_
e
k(n/l)
2
t
,
where the c
n
are the Fourier coefcients of the initial prole,
c
n
=
2
l
_
l
0
(x) sin(nx/l) dx.
0
0.1
0.2 0
1
0
0.5
1
x/l kt/l
2
In particular, if (x) = 1 for 0 < x < l,
we can use the coefcients from Example 1
and obtain the solution
u(x, t) =

n=1
2(1(1)
n
)
n
sin
_
nx
l
_
e
k(n/l)
2
t
.
A plot of the partial sum with 21 terms
shows the rapid decay of the higher fre-
quency terms and the disappearance of the
Gibbs oscillation.
Wave equation with Dirichlet boundary conditions As we also saw in lec-
ture 6, the solution of the IBVP
u
tt
c
2
u
xx
= 0 (x (0, l), t > 0)
u(0, t) = 0, u(l, t) = 0 (t 0)
u(x, 0) = (x), u
t
(x, 0) = (x) (x (0, l))
is
u(x, t) =

n=1
_
A
n
cos(nct/l) +
l
nc
B
n
sin(nct/l)
_
sin(nx/l),
where the A
n
and B
n
are the Fourier coefcients of the initial proles,
A
n
=
2
l
_
l
0
(x) sin(nx/l) dx, B
n
=
2
l
_
l
0
(x) sin(nx/l) dx.
In particular, if (x) = 0 and (x) = 1 for 0 < x < l, we can use the coefcients
from Example 1 and obtain the solution
u(x, t) =

n=1
2(1 (1)
n
)l
n
2

2
c
sin(
nct
l
) sin(
nx
l
)
=

n=1
2(1 (1)
n
)l
n
2

2
c
_
1
2
cos(
n(ct x)
l
)
1
2
cos(
n(ct +x)
l
)
_
.
The second formula has the same formas the general solution of the wave equation
f(x+ct) +g(xct). A plot of a partial sum with a large number of terms shows
the time-periodic response:
37
0
1
2
0
1
!0.1
0
0.1
x/l
ct/l
8.5 Heat equation with source term
Consider the heat equation with Dirichlet boundary conditions
u
t
(u
x
)
x
= 0 (x (0, l), t > 0)
u(0, t) = 0, u(l, t) = 0 (t 0)
u(x, 0) = (x) (x (0, l))
By the method of separation of variables, we have found that the source operator
S(t) for this problem transforms the initial prole (x) into the function

n
(, X
n
)
X
n

2
X
n
(x)e

n
t
,
where X
n
are -orthogonal eigenfunctions corresponding to the eigenvalues
n
.
Then, by Duhamels principle, the solution of the problem with source term, that
is, of
u
t
(u
x
)
x
= (x)f(x, t) (x (0, l), t > 0)
u(0, t) = 0, u(l, t) = 0 (t 0)
u(x, 0) = 0 (x (0, l))
_
_
_
(24)
is
u(x, t) =
_
t
0

n
(f(), X
n
)
X
n

2
X
n
(x)e

n
(t)
d.
Assuming that integration and summation commute, the solution can be written
as
u(x, t) =

n
u
n
(t)X
n
(x),
where
u
n
(t) =
_
t
0
f
n
()e

n
(t)
d (25)
with
f
n
() =
(f(), X
n
)
X
n

2
=
_
l
0
()f(, )X
n
() d
_
l
0
()(X
n
())
2
d
.
38
Here is an alternative derivation of the solution of (24). Substituting an as-
sumed solution of the form u(x, t) =

m
u
m
(t)X
m
(x) into the PDE gives
(x)f(x, t) =

m
(u

m
X
m
u
m
(X

m
)

) =

m
(u

m
+
m
u
m
)X
m
.
Multiplying through by X
n
/X
n

2
and integrating gives the decoupled ODEs
1
X
n

2
(f(t), X
n
)
. .
f
n
(t)
=

m
(u

m
+
m
u
m
)(X
m
, X
n
)/X
n

2
= u

n
(t) +
n
u
n
(t),
each of which (with initial condition u
n
(0) = 0) has the solution (25).
For example, for the constant-coefcient heat equation problem
u
t
ku
xx
= 1 (x (0, l), t > 0)
u(0, t) = 0, u(l, t) = 0 (t 0)
u(x, 0) = 0 (x (0, l))
0
1 0
0.2
0.4
0
0.05
0.1
kt/l
2
x/l
we have
f
n
(t) =
2
l
_
l
0
sin(n/l) d =
2(1 (1)
n
)
n
and so u(x, t) =

n=1
u
n
(t) sin(nx/l)
with
u
n
(t) =
2(1 (1)
n
)
n
_
t
0
e
k(n/l)
2
(t)
d
=
2(1 (1)
n
)
kn
3

3
/l
2
(1 e
k(n/l)
2
t
).
39
9 Laplaces Equation
Vector Analysis facts
Diffusion and Heat ow in three dimensions
Membrane vibration
Laplaces equation
9.1 Some Facts from Vector Analysis
Let D be an open simply connected spatial domain with surface D, let n be the
outward unit normal, and let f be a vector eld. Gausss divergence theorem is
_
D
f dV =
_
D
f ndA,
where the divergence of the vector eld having cartesian coordinates f = f
1
i +
f
2
j +f
3
k is
f =
f
1
x
+
f
2
y
+
f
3
z
.
The gradient of a scalar eld u(x) in cartesian coordinates is the vector
u =
u
x
i +
u
y
j +
u
z
k.
Applying the divergence theorem to f = vu, where u and v are scalar elds,
yields Greens rst identity
_
D
(vu) ndA =
_
D
(v u +vu) dV
where u = u = u
xx
+u
yy
+u
zz
is the laplacian of u.
In two dimensional problems, D denotes a two-dimensional domain and D
denotes the closed curve that is its boundary. The divergence theorem is written
_
D
f dA =
_
D
f ndl,
and Greens rst identity is written
_
D
(vu) ndl =
_
D
(v u +vu) dA.
9.2 Heat Flow in Three Dimensions
Consider a substance (e.g. mass or energy) owing in a region of space. Let
u(x, t) denote its density (units: [quantity] [volume]
1
) as a function of position
x = [x, y, z] and time t, and let

(x, t) denote the ux vector. (units: [quantity]
[time]
1
[area]
1
). The amount of substance in a domain D is given by the
volume integral
_
D
u(x, t) dV .
40
Letting n denote the unit normal vector on the surface D of the domain D, the
net ux out of the domain is given by
_
D

ndA. Let f(x, t, u) denote the


source term, that is, the rate (units: [quantity] [time]
1
[volume]
1
) at which
substance density increases by processes other than ux, for example chemical
reaction. The rate of increase of the total amount of substance in the interval is
then
d
dt
_
D
u(x, t) dV =
_
D

ndA +
_
D
f(x, t, u)dV.
Using the divergence theorem, the surface integral can be replaced by a volume
integral, yielding
_
D
_
u
t
+

f
_
dV = 0.
Because D is arbitrary, this implies that the conservation equation
u
t
+

= f
should hold at every point in the region .
Diffusion processes, whereby substance ows from areas of high concentra-
tion to areas of low concentration, can be modelled by the constitutive relation
(Ficks law)

= ku,
where k(x) is a material parameter (diffusivity, units: [length]
2
[time]
1
). Sub-
stituting this into the conservation equation gives the three-dimensional diffusion
equation
u
t
(ku) = f.
When k is constant, the equation reduces to u
t
ku = f.
Alternatively, we can let u = cT denote density of heat energy, where c(x)
and (x) are material parameters (specic heat and mass per unit length) and T is
the temperature. The constitutive relation (Fouriers law)

= KT
models conduction, whereby heat ows from hot areas to colder areas. The ma-
terial parameter K(x) is called the heat conductivity. Substituting Fouriers law
into the conservation equation gives the three-dimensional heat equation
cT
t
(KT) = f.
This equation is very similar to the diffusion equation, and is essentially identical
to it (except for notation) when c is constant.
9.3 Membrane Vibration
Let u(x, y, t) denote the displacement of a thin membrane that moves in the z di-
rection (vertical) only.
41
Horizontal force balance: Let T(x, y, t) be
the tension (units [force] [length]
1
), as-
sumed to act tangentially along the mem-
brane. Let D be a domain in the xy
(horizontal) plane, let D be its bound-
ary curve, and let n denote the unit out-
ward normal vector (in the xy plane) on
D. Let u
n
= n u denote the direc-
tional derivative of u in the direction n; then
= tan
1
(n u) is the angle between n
and the tension vector. Because there is no horizontal motion, the vector sum of
the horizontal forces acting on the boundary must be zero, that is,
_
D
T cos ndl = 0.
Then, for any constant vector a, we have
a 0 = a
_
D
T cos ndl =
_
D
(aT cos ) ndl
=
_
D
(aT cos ) dA = a
_
D
(T cos ) dA,
and since a is arbitrary, this implies
_
D
(T cos ) dA = 0. Because the domain
is arbitrary, this in turn implies (T cos ) = 0, so that T cos is constant with
respect to x and y, say
T(x, y, t) cos (x, y, t) = (t).
Mass conservation: Let (x, y, t) be the membranes mass per unit area (
may vary as the membrane deforms during the motion), and let
0
(x, y) be the
mass per unit area when the membrane is plane. If dA

represents an area element


of deformed membrane and dA represents the same element when the membrane
is plane, then mass conservation requires that dA

=
0
dA.
Vertical force balance: Consider a membrane piece whose projection onto the
xy plane is D. By Newtons law, the net vertical force on this piece is equal to the
time derivative of the momentum:
d
dt
_
D
u
t
dA

=
_
D
T sin dl =
_
D
T cos tan dl
=
_
D
u ndl =
_
D
(u) dA.
Using mass conservation gives
_
D
(
0
u
tt
u) dA = 0,
which implies

0
u
tt
u = 0.
Denoting c =
_
/
0
, this can be written
u
tt
= c
2
u,
which is the two-dimensional wave equation.
42
9.4 Laplaces Equation
If the source term in the diffusion equation is constant, then the steady state equi-
librium concentration is described by the diffusion equation with the time deriva-
tive terms removed:
(ku) = f.
When k is constant this reduces to the Poisson equation
u = F
where F = f/k. The Poisson equation with no source term is Laplaces equation
u = 0.
Similar equations arise as models of steady state heat ow. A membrane sub-
jected to a transversal static load f(x) is modelled by a two-dimensional Poisson
equation
u = f.
Laplaces and Poissons equations also arise as models of gravitational elds, elec-
trostatic elds, stationary uid ow, brownian motion, and many other phenom-
ena.
Any function that satises Laplaces equation is called a harmonic function.
In one dimension, Laplaces equation is u
xx
= 0, so one-dimensional harmonic
functions are all of the formu(x) = A+Bx. Things get more interesting in higher
dimensions, however! The following results holds in the one, two, and three (and
higher!) dimensional versions of Laplaces equation.
Theorem 8 (Maximum Principle) If u is harmonic in a connected bounded open
set , and continuous in

= , then the maximum value of u is attained on
the boundary .
Proof. Let > 0 and v(x) = u(x) +

N
i=1
x
2
i
. If v has a maximum at a point
x , the hessian matrix v
xx
= [

2
v
x
i
x
j
] is negative semidenite, which implies

2
v
x
2
i
0 for every i. But
v =

2
v
x
2
i
= u
..
=0
+2N > 0,
so v cannot have a maximum inside . Because v is continuous, it has a maximum
somewhere in the compact set

, say at x

. Then, for all x



,
u(x) v(x) v(x

) = u(x

) +

i
(x

i
)
2
u(x
0
) + max
x
|x|
2
where x
0
is a point where u(x
0
) = max
x
u(x). Then taking 0, we
have
u(x) u(x
0
) for all x

,
which completes the proof.
Corollary 1 (Minimum Principle) For u as in Theorem 1, the minimum value is
attained on the boundary .
Proof. Apply the maximum principle to the harmonic function u. 2
43
Example 1 Find the maximum value of f(x, y) = x
2
y
2
in the unit disk
x
2
+y
2
1.
Solution. Because f = 0, the maximum occurs on the disk boundary. Using po-
lar coordinates, f = cos
2
sin
2
= cos(2) on the boundary, and the maximum
value is 1, attained at the points (x, y) = (1, 0) and at (x, y) = (1, 0). 2
Example 2 Prove that the Dirichlet problem
u = f in connected bounded , u = h on .
has at most one solution.
Solution. If u
1
and u
2
are two solutions, their difference w = u
1
u
2
is harmonic
in and w = 0 on . By the maximum/minimum principle,
0 = min
x
w(x) w(x) max
x
w(x) = 0,
for all x

, and so w 0, that is, u
1
u
2
.
An alternative solution is based on Greens rst identity:
_

w
..
=0
w ndA =
_

|w|
2
+w w
..
=0
dV
Thus w 0, so w is constant, and the constant is zero because w is continuous
and is zero on the boundary. 2
The following result indicates that the laplacian is suitable for modelling isotropic
physical phenomena, in which there is no preferred direction. A rotation of the
coordinate axes corresponds to a linear transformation x

= Bx with orthogonal
B (that is, B
T
B = I) and det(B) = 1. (An orthogonal B with det(B) = 1
models a rotation with reection.)
Theorem 9 (Rotational invariance of laplacian) If u

(x

) = u(B
T
x

) with or-
thogonal B then

= u.
Proof. By the chain rule we have
u

i
=

l
u
x
l
x
l
x

i
..
b
il
and

2
u

i
x

j
=

2
u
x
l
x
p
b
il
b
jp
= (Bu
xx
B
T
)
ij
.
Then

= tr(u

x
) = tr(Bu
xx
B
T
) = tr(B
T
Bu
xx
) = tr(u
xx
) = u.
Theorem 10 (Mean value property) If u is harmonic in the ball D = {x : |x
x

| = a} then the value at the centre u(x

) is equal to the mean value on the ball


boundary.
44
Proof. Moving the origin to x

, the mean value of u on the boundary of the three


dimensional ball is given by
m(a) =
_
D
u dA
A
=
1
4a
2
_
2
0
_

0
u(a, , )a
2
sin d d.
From Greens rst identity with v
1
4a
2
and u = 0 we have
0 =
_
|x|=a
vu ndA =
1
4
_
2
0
_

0
u
r
(a, , ) sin d d = m

(a),
so that the mean value on the ball surface is independent of the balls radius.
Taking a 0 gives m = u(0). The proof for one and two dimensions is similar. 2
According to the maximum principle (Theorem 1), harmonic functions attain
their maximum on the boundary. Using the mean value property we can show that
the maximum is not attained inside the region, unless the function is constant.
Corollary 2 If u is harmonic in an open connected region , is continuous in

= , and attains its maximum in , then u is constant in .


Section 2.4, Laplaces Equation 47
x
0
x
1
x
2
x
m
!
"
x
3
Figure 2.1: Diagram for Proof of Maximum Principle.
Proof. Suppose u attains its maximum M := max

u at a point x
0
.
We wish to show that at any other point x
m
we must have u(x
m
) = M.
Let the curve connect x
0
and x
m
, and choose the nite set of points
x
1
, x
2
, . . . x
m1
on to be centers of balls contained , and arranged so that
the point x
i +1
lies on the surface B
i
of the ball B
i
centred at the previous
point x
i
. The values on B
0
are all less than or equal to M. But, by the
mean value property (Theorem 2.17) u(x
0
) must be equal to the average of
the values on the balls surface, and so the surface values must all be equal to
M. In particular, u(x
1
) = M. With similar arguments we obtain u(x
i
) = M
for i = 2, 3, . . . m (Figure 2.1). The proof for the minimum is similar.
From Theorem 2.18 we can obtain the results of section 2.2.1 on con-
tinuous dependence on boundary data and monotonicity of solutions of the
Dirichlet problem.
2.4.4 Existence of Solution
This chapter has given several uniqueness results but has not yet said anything
about the existence of the solution. We close the chapter with a few words
about this.
The Dirichlet problem can in fact fail to have a solution if there are sharp
enough spikes that penetrate into the domain . In the absence of such
spikes, however, a solution will exist; see [9, p.198] for details. Domains
encountered in applications are unlikely to cause trouble in this regard.
An alternative is to replace the PDE by an integral formulation of the
boundary value problem that doesnt require so much smoothness in the
solution. Such variational or weak formulations are the starting point for the
theory of numerical methods such as the Finite Element Method.
Proof. Suppose u attains its maximum
M := max

u at a point x
0
. We wish
to show that at any other point x
m

we must have u(x
m
) = M. Let the curve
connect x
0
and x
m
, and choose
the nite set of points x
1
, x
2
, . . . x
m1
on
to be centers of balls contained , and ar-
ranged so that the point x
i+1
lies on the sur-
face B
i
of the ball B
i
centred at the previ-
ous point x
i
. The values on B
0
are all less than or equal to M. But, by the mean
value property u(x
0
) must be equal to the average of the values on the balls sur-
face, and so the surface values must all be equal to M. In particular, u(x
1
) = M.
With similar arguments we obtain u(x
i
) = M for i = 2, 3, . . . m.
x
y
z
r
!
!
In spherical coordinates the laplacian is
u =
1
r
2
_
r
2
u
r
_
r
+
1
r
2
sin
(u

sin )

+
1
r
2
sin
2

Harmonic functions in three dimensions that depend only on r satisfy


the ODE (r
2
u
r
)
r
= 0, which has the general solution u(r) = c
1
r
1
+
c
2
.
In cylindrical polar coordinates the laplacian is
u =
1
R
_
(Ru
R
)
R
+
1
R
(u

+ (Ru
z
)
z
_
.
Harmonic functions that depend only on distance R from the z-axis satisfy the
ODE (Ru
R
)
R
= 0, which has the general solution u(R) = c
1
ln(R) +c
2
.
45
Example 3 Acircular membrane with a uniformtransverse load and xed bound-
ary is modelled by u = 1 in the domain x
2
+y
2
< a
2
with u = 0 on the boundary.
Find the shape of the membrane.
Solution. Assuming a solution of the formu(R), we have the ODE R
1
(Ru
R
)
R
=
1, which has the solution u(R) = c
1
ln(R) + c
2
+
1
4
R
2
. Taking c
1
= 0 (to ensure
the solution remains bounded on the z axis) and c
2
= a
2
/4 to ensure u(a) = 0,
we have u(R) = (R
2
a
2
)/4.
Theorem 11 (Dirichlets Principle) Among all the functions w that satisfy w =
h on , the one that minimizes the energy E(w) =
1
2
_

|w|
2
dV is har-
monic.
Proof. Let u = w = h on and u = 0 in . Then, with v = u w, we have
E(w) =
1
2
_

|(u v)|
2
dV
=
1
2
_

(u v) (u v) dV
= E(u) + E(v)
_

u v dV
Greens rst identity gives
_

u v dV =
_

v
..
=0
u ndA
_

v u
..
=0
dV = 0,
so E(w) E(u). 2
The Dirichlet principle motivates the Rayleigh-Ritz method for computing an
approximate solution of the Dirichlet problem of Example 2. Choose functions
w
0
, w
1
, . . . , w
n
such that w
0
= h and w
1
= = w
n
= 0 on , and consider
the linear combination w = w
0
+c
1
w
1
+ +c
n
w
n
. Then
E(w) =
1
2
_

_
w
0
+
n

i=1
c
i
w
i
_

_
w
0
+
n

j=1
c
j
w
j
_
dV
=
1
2
a b
T
c +
1
2
c
T
Ac,
where a =
_

|w
0
|
2
dV , b
i
=
_

w
0
w
i
dV , and A
ij
=
_

w
i
w
j
dV .
This energy is minimized when c = A
1
b, because
E(w) =
1
2
a b
T
(c A
1
b +A
1
b) +
1
2
(c A
1
b +A
1
b)
T
A(c A
1
b +A
1
b)
=
1
2
(a b
T
A
1
b) +
1
2
(c A
1
b)
T
A(c A
1
b),
and A is symmetric positive denite.
Example 4 Find an approximate solution to u = 0 on the triangle {(x, y) :
x > 0, y > 0, 3x +y < 3} with the boundary conditions
u(x, 0) = 0, u(0, y) = 3y y
2
, u(x, 3 3x) = 0.
46
0
1
0
1
2
3
0
1
2

x
y
u
Solution. With w
0
= (33xy)y and w
1
= (33xy)xy
we have
b
1
=
_
1
0
_
33x
0
w
0
x
w
1
x
+
w
0
y
w
1
y
dy dx =
9
20
A
11
=
_
1
0
_
33x
0
_
w
1
x
_
2
+
_
w
1
y
_
2
dy dx =
3
2
c
1
= b
1
/A
11
=
3
10
so the approximate solution is w
0

3
10
w
1
= y(3 3x y)(1 0.3x).
47
10 Solving Two-Dimensional Laplace Equations
Laplace equation boundary value problems in a disk, a rectangle, a wedge,
and in a region outside a circle
10.1 Dirichlet Problem in a disk
Consider the two dimensional Laplace equation in the disk x
2
+ y
2
< a
2
with
u = h on the boundary. The Laplace equation in polar coordinates (r, ) is
(ru
r
)
r
+
1
r
(u

= 0
A separation of variables trial solution u(r, ) = R(r)() gives
R

+
1
r
R

+
1
r
2
R

= 0
which can be rearranged to
r
2
R

R
+
rR

R
=

.
Equating both sides to the separation constant , we are left with two ODEs,

+ = 0 (26)
and
r
2
R

+rR

R = 0. (27)
For = 0, the general solution of (26) is () = Acos

+ Bsin

.
Substituting this into the periodic boundary conditions (0) = (2) and

(0) =

(2) gives the homogeneous equations


A(1 + cos(2

)) + Bsin(2

) = 0
A

sin(2

) +B

(1 + cos(2

)) = 0
This system of equations has a nontrivial solution (that is, a solution other than
A = B = 0) if the determinant of the coefcient matrix is zero. The determinant
is 2

(1cos(2

)), and the nonzero values of that give a nontrivial solution


are
n
= n
2
with n = 1, 2, . . .
For = 0, the general solution of (26) is () = A + B, and the only
nontrivial periodic solution is () = A (nonzero constant).
Consider now the equation (27). For = 0, a trial solution of the formR(r) =
r

gives
(
2
)r

= 0,
which implies the solutions =

= n. For = 0, (27) reduces to


r(rR

= 0, which has the general solution R(r) = c


1
ln r +c
2
.
48
Writing the solution as a linear combination of the solutions found above, we
have
u(r, ) = c
1
ln r +c
2
+

n1
(A
n
cos n +B
n
sin n)r
n
+

n1
(

A
n
cos n +

B
n
sin n)r
n
(28)
To ensure continuity and boundedness at the origin, we set c
1
=

A
n
=

B
n
= 0,
leaving us with
u(r, ) =
1
2
A
0
+

n1
(A
n
cos n +B
n
sin n)r
n
where the constant term c
2
has been renamed
1
2
A
0
. At the boundary r = a we
have
u(a, ) =
1
2
A
0
+

n1
(A
n
cos n +B
n
sin n)a
n
= h()
and the coefcents are determined by equating A
n
a
n
and B
n
a
n
with the Fourier
series coefcients of h:
u(r, ) =
1
2
_

h(

) d

n1
__
1

h(

) cos n

_
cos n
+
_
1

h(

) sin n

_
sin n
_
_
r
a
_
n
(29)
Example 1 Find the steady-state temperature in a long cylinder of radius a if
the upper half is kept at u = 100 and the lower half is kept at u = 0.
Solution. The boundary function is h() = 50 + 50f() where
f() =
_
_
_
1 < < 0
0 = 0
1 0 < < .
!
-1
1
-! "
In section 8 we derived the Fourier series
f() =

n1
2(1 (1)
n
)
n
sin(n),
so the PDE solution can be written directly as
u(r, ) = 50+50

n1
2(1 (1)
n
)
n
_
r
a
_
n
sin(n)
Plotting the sum of a large number of terms gives the above gure. 2
49
The series solution (29) for the Dirichlet problem in the disk can be written as
u(r, ) =
1
2
_

h(

)
_
1 + 2

n1
_
r
a
_
n
cos n(

)
_
d

(30)
Letting z =
r
a
e
i(

)
, the term in brackets can be written
1 +

n1
z
n
+ z
n
= 1 +
z
1 z
+
z
1 z
=
1 z z
(1 z)(1 z)
=
1 |z|
2
1 +|z|
2
(z + z)
=
1 (r/a)
2
1 + (r/a)
2
2(r/a) cos(

)
Substituting this into (30) gives the Poisson integral formula
u(r, ) =
1
2
_

h(

)
_
a
2
r
2
a
2
+r
2
2ar cos(

)
_
d

.
One application of this formula is to provide an alternative derivation of the mean
value property: the value at r = 0 is
1
2
_

h(

) d

, which is the average of u on


the circumference r = a.
10.2 Dirichlet Problem in a rectangle
a
b
y
x
u = g
u = 0
u = 0 u = 0 u
xx
+ u
yy
= 0
Consider the two dimensional Laplace equation u = 0 in
the rectangle (0, a) (0, b) with u(x, b) = g(x) and u = 0 on
the rest of the boundary. Substituting a solution of the form
u(x, y) = X(x)Y (y) into u
xx
+u
yy
= 0 gives
X

Y +XY

= 0,
which can be rearranged to
X

X
=
Y

Y
.
Equating both sides to the separation constant , we are left with two ODEs,
X

+X = 0 (31)
and
Y

Y = 0. (32)
For = 0, the general solution of (31) is X = Acos

x + Bsin

x. Substi-
tuting this into the boundary conditions X(0) = 0 and X(a) = 0 gives A = 0
and Bsin(

a) = 0, and the nonzero values of that give a nontrivial solution


are
n
= (n/a)
2
with n = 1, 2, . . .. For = 0, the general solution of (31)
is X(x) = A + Bx, and there is no nontrivial solution satisfying the boundary
conditions.
The general solution of (32) with = (n/a)
2
is Y (y) = Acosh(ny/a) +
Bsinh(ny/a). The boundary condition Y (0) = 0 is satised by setting A = 0.
50
Writing the solution as a linear combination of the solutions found above, we
have
u(x, y) =

n1
B
n
sin(nx/a) sinh(ny/a)
At the boundary y = b, this is
u(x, b) =

n1
B
n
sin(nx/a) sinh(nb/a) = g(x).
The coefcients are obtained by equating B
n
sinh(nb/a) with the Fourier sine
series coefcients of g(x):
u(x, y) =

n1
_
2
a
_
a
0
g(x

) sin(nx

/a) dx

_
sin(nx/a) sinh(ny/a)
sinh(nb/a)
(33)
Example 2 Find the steady-state temperature in a long prismatic tube with square
a a cross-section if the top face is kept at u = 1 and other three faces are kept
at u = 0.
Solution. The boundary function is g(x) = 1, whose Fourier sine series was found
in section 8 to be
g(x) =

n1
2(1 (1)
n
)
n
sin(nx/a).
The Dirichlet problem solution is then
u(x, y) =

n1
2(1 (1)
n
)
n
sin(nx/a) sinh(ny/a)
sinh(n)
Plotting the sum of a large number of terms gives the above gure. 2
The solution for the general Dirichlet problem in the rectangle is found by
superposition of (33) and solutions of similar problems:
For u(x, 0) = h(x), and u = 0 elsewhere
0 0
0
h
, the solution is
u(x, y) =

n1
_
2
a
_
a
0
h(x

) sin(nx

/a) dx

_
sin(nx/a) sinh(n(b y)/a)
sinh(nb/a)
This is found by replacing y by by [the laplacian is invariant to reection]
and g by h in (33).
For u(a, y) = k(y), and u = 0 elsewhere
0
0
0
k
, the solution is
u(x, y) =

n1
_
2
b
_
b
0
k(y

) sin(ny

/b) dy

_
sin(ny/b) sinh(nx/b)
sinh(na/b)
which is found by swapping x y and a b and replacing g by k in (33).
51
for u(0, y) = j(y), and u = 0 elsewhere
0
0
0
j
, the solution is
u(x, y) =

n1
_
2
b
_
b
0
j(y

) sin(ny

/b) dy

_
sin(ny/b) sinh(n(a x)/b)
sinh(na/b)
.
10.3 Dirichlet-Neumann Problem in a wedge
u
r
= h
u = 0
u = 0
!
Consider the two dimensional Laplace equation in the sector
{(r, ) : 0 < < , r < a}, with boundary conditions
u = 0 on the rays = 0 and = and a Neumann condition
u
r
= h on the perimeter r = a. Assuming a solution of the
form u(r, ) = R(r)() and proceeding as for the disk, we
obtain the ODEs (26) and (27). From the boundary conditions (0) = () = 0
we obtain () = sin(n/) with n = 1, 2, . . . and R(r) = r
n/
. Writing the
solution as a linear combination, we have
u(r, ) =

n1
B
n
sin(n/)r
n/
+

B
n
sin(n/)r
n/
To ensure boundedness at the origin, we set

B
n
= 0, leaving
u(r, ) =

n1
B
n
sin(n/)r
n/
Substituting this into the Neumann boundary condition gives
h() = u
r
(a, ) =

n1
n

B
n
sin(n/)a
n

1
and the coefcients are determined by equating
n

B
n
a
n

1
with the Fourier sine
coefcients of h(), leading nally to
u(r, ) =

n1
a
n
_
2

_

0
h(

) sin(n

/) d

_
sin(n/)
_
r
a
_
n/
In a nonconvex sector having > , the term r

1
in the rst term of
the series for u
r
is unbounded as r 0, and so this solution does not, strictly
speaking, satisfy the PDE in fact, the problem does not have a solution (having
second derivatives continuous up to the boundary) in this case.
Example 3 Solve u = 0 in a = /2 sector with u = 0 on the rays and
u
r
= 1 on the perimeter r = a.
Solution. We have
2

_

0
h(

) sin(n

/) d

=
2(1 (1)
n
)
n
and so
u(r, ) =

n1
(1 (1)
n
)a
n
2

sin(2n)
_
r
a
_
2n
.
Notice that Gibbs oscillation is not visible in the
plot of the solution surface. 2
52
10.4 Dirichlet Problem in the region outside a circle
Consider the two dimensional Laplace equation in the region x
2
+ y
2
> a
2
with
u = h on the boundary and u bounded at innity.
We can proceed as in section 10.1 up to formula (28). To ensure the solution
is bounded at innity, we set c
1
, A
n
and B
n
to zero, leaving
u(r, ) =
1
2

A
0
+

n1
(

A
n
cos n +

B
n
sin n)r
n
Imposing the boundary condition on the circle perimeter gives
h() = u(a, ) =
1
2

A
0
+

n1
(

A
n
cos n +

B
n
sin n)a
n
The coefcients are determined by equating a
n

A
n
and and a
n

B
n
to the Fourier
coefcients of h, leading to
u(r, ) =
1
2
_

h(

) d

n1
__
1

h(

) cos n

_
cos n
+
_
1

h(

) sin n

_
sin n
_
_
a
r
_
n
. (34)
An alternative derivation is to use the change of variables r

= a
2
/r and u

(r

, ) =
u(a
2
/r

, ). By the chain rule,


u

r
=
u

=
u
r
r
r

=
a
2
(r

)
2
u
r
and
1
r

(r

r
)
r
=
1
r

(
a
2
r

u
r
)
r
=
r
a
2
(ru
r
)
r
r
r

=
_
r
a
_
4
1
r
(ru
r
)
r
and so

=
1
r

(r

r
)
r
+
1
(r

)
2
(u

=
_
r
a
_
4
u.
Thus, if

= 0 inside the circle then u = 0 outside it. The solution inside the
circle is (29) written with u

and r

, and applying the change of variables to this


solution gives (34). Similarly, the Poisson integral formula for the disk interior is
transformed to
u(r, ) =
1
2
_

h(

)
_
r
2
a
2
a
2
+r
2
2ar cos(

)
_
d

.
for the exterior.
Example 4: Stationary ow past a circular cylinder The steady-state two-
dimensional velocity eld for the irrotational ow of an incompressible inviscid
constant-density uid is given by
y
i
x
j, where is a harmonic function called
the stream function. Because the velocity is orthogonal to , the lines of con-
stant (called streamlines) are tangential to the velocity eld.
53
Find the streamlines for ow past a long circular cylinder of radius a whose
axis is the z axis, assuming the ow far from the cylinder to be constant in the x
direction, that is, = Uy.
Solution. The stream function satises the two dimensional Laplace equation on
the exterior of the circle r = a. The radial component of the velocity is zero on
the cylinder boundary, that is, the cylinder boundary is a streamline, which gives
the boundary condition =constant (say, zero) for r = a. Substituting this into
the general solution (28) we nd
c
1
ln a +c
2
= 0
A
n
a
n
+

A
n
a
n
= 0
B
n
a
n
+

B
n
a
n
= 0
so that
(r, ) = c
1
ln
r
a
+

n1
_
r
n

a
2n
r
n
_
(A
n
cos n +B
n
sin n)
In order to satisfy the condition = Uy = Ur sin at r we set B
1
= U
and the remaining A
n
and B
n
coefcients to zero, leaving
(r, ) = c
1
ln
r
a
+U
_
1
a
2
r
2
_
r sin .
Here are streamlines for various c
1
values, which correspond to different cylinder
clockwise rotation speeds:
c
1
= 0 c
1
= Ua
c
1
= 2Ua c
1
= 2.1Ua
The closely spaced streamlines correspond to regions of low pressure and indicate
the presence of a net lift force in the y direction (Magnus effect).
54
11 Greens Functions
11.1 Greens Function for One-Dimensional Equation
The Greens function provides a complete solution to a boundary value problem in
much the same way that an inverse matrix provides a general solution for systems
of linear equations. In this section the Greens function is introduced in the context
of a simple one-dimensional problem.
Some of the proofs use the identity
_
b
a
(uv

vu

) dx =

b
a
uv

vu

.
This can be obtained by integrating (uv

vu

= uv

vu

.
Asingularity function K(x, ) of the operator Ldened by Lu(x) = u

(x)
c(x)u(x) is characterised by three properties:
1. K is continuous;
2. K
x
is continuous in x < and in x > , and K
x
(x
+
, x) K
x
(x

, x) = 1;
3. K
xx
is continuous and LK = 0 for x = .
Note that the three properties do not dene a singularity function uniquely: if
K is a singularity function then so is K + H, where H(x, ) is any function with
continuous H and H
x
and with LH = 0.
The Greens function G(x, ) for the operator L and the domain (a, b) with
Dirichlet boundary conditions is the singularity function that satises the homo-
geneous Dirichlet conditions G(a, ) = 0 and G(b, ) = 0. The Greens function
provides the solution to the boundary value problem with Dirichlet boundary con-
ditions:
Theorem 12 If u satises the differential equation u

+ cu = f on (a, b) and
the boundary conditions u(a) = 0, u(b) = 0, then
u() =
_
b
a
G(x, )f(x) dx (35)
for all (a, b).
Proof. Letting v(x) = G(x, ), we have
_
b
a
vf dx =
_

a
uv

vu

dx +
_
b

+
uv

vu

dx
=

a
uv

vu

+
uv

vu

b
a
uv

vu

. .
0
u()

. .
1
+

vu

. .
0
,
which completes the proof. 2
55
The load (or source) function f in the differential equation u

cu = f
can be thought of as a superposition of point loads f(x) =

(x )f()d,
where (x ) is concentrated at and has unit magnitude (i.e.
_
dx = 1).
Then formula (35) represents the solution as a weighted sum of Greens functions,
where each G(, x) is the solution to u

(x) = (x ). The Greens function


can thus be thought of as the response to a unit point load.
The Greens function also provides the solution of the boundary value problem
with nonhomogeneous boundary conditions:
Theorem 13 The solution of u

+ cu = 0 with boundary conditions u(a) = h


0
and u(b) = h
1
satises
u() = G
x
(a, )h
0
G
x
(b, )h
1
(36)
for all (a, b).
Proof. Denoting v(x) = G(x, ), we have
0 =
_

a
(uv

vu

) dx +
_
b

+
(uv

vu

) dx
=

a
(uv

vu

) +

+
(uv

vu

)
=

b
a
(uv

vu

(uv

vu

)
=

b
a
uv

+u(),
which completes the proof. 2
The following result tells us that the response at x to a point load applied at
is equal to the response at to a point load applied at x.
Theorem 14 (Reciprocity Principle) G(x, ) = G(, x) for all x, (a, b).
Proof. Let y and be distinct points in (a, b) with y < , let u(x) = G(x, ) and
v(x) = G(x, y). Then
0 =
_
y

a
uv

vu

dx +
_

y
+
uv

vu

dx +
_
b

+
uv

vu

dx
=

a
uv

vu

y
+
uv

vu

+
uv

vu

b
a
uv

vu

. .
0

y
+
y

uv

vu

. .
u(y)

uv

vu

. .
v()
,
leaving u(y) v() = 0, that is, G(y, ) G(, y) = 0. The proof for y > is
similar. 2
As a consequence of Theorem 2, we can rewrite formula (35) as
u(x) =
_
b
a
G(x, )f() d.
56
and formula (36) as
u(x) = G

(x, a)h
0
G

(x, b)h
1
.
A singularity function of
d
2
dx
2
is given by
K(x, ) =
1
2
|x |,
as can readily be veried. The Greens function for the interval (0, 1) can be
0
1
0
1
0
0.25
x
!
G
found by solving H
xx
= 0 with boundary conditions
H(0, ) = K(0, ) and H(1, ) = K(1, ), then
setting G = K +H. This yields
G(x, ) =
1
2
|x | +
1
2
(x +) x
=
_
(1 x) for < x
(1 )x for x < .
The solution of u

= f satisfying u(0) = h
0
and u(1) = h
1
is given by
u(x) =
_
1
0
G(x, )f() d +G

(x, 0)h
0
G

(x, 1)h
1
= (1 x)
_
x
0
f() d +x
_
1
x
(1 )f() d + (1 x)h
0
+xh
1
.
11.2 Greens Function for Two-Dimensional Poisson Equation
Now we go through the same discussion in two dimensions. Some of the proofs
use Greens second identity
_

uv vudA =
_

(uv vu) ndl,


This can be derived by writing Greens rst identity twice, with u and v inter-
changed the second time, and subtracting.
A singularity function K(x, x

) of the operator is characterised by the


three properties
1. For any xed x

, lim
0
_
B

K(x, x

) dl = 0, where B

denotes the radius-


disk centred at x

;
2. For any xed x

, lim
0
_
B

K(x, x

) ndl = 1, where n denotes the


outward unit normal to B

;
3. K is harmonic as a function of x for x = x

.
Note that these three properties do not dene a singularity function uniquely:
if K is a singularity function then so is K + H, where H(x, x

) is any function
that is harmonic as a function of x.
57
A singularity function of is given by
K(x, x

) =
1
2
ln |x x

|. (37)
This assertion can be veried as follows. Without loss of generality we can take
x

= 0. In polar coordinates, we have K =


1
2
ln r with r = |x|, which is
harmonic in R
2
\ {0} (see section 9.4). Also,
_
B

K dl =
1
2
_
2
0
ln d = ln 0
and (because K n = K/r)
_
B

K ndl =
1
2
_
2
0
1

d = 1.
The Greens function for and a domain with Dirichlet boundary con-
ditions is a singularity function that satises G(x, x

) = 0 for x . The
Greens function provides the solution to the Poisson equation with homogeneous
Dirichlet boundary conditions:
Theorem 15 If u = f in and u = 0 on then
u(x

) =
_

G(x, x

)f(x) dA (x

). (38)
Proof. Letting v(x) = G(x, x

), we have
_
\B

vf dA =
_
\B

uv vudA
=
_

(uv vu) ndl


. .
0

_
B

(uv vu) ndl


u(x

)
_
B

v
n
dl
. .
1
+
_
B

vu
n
dl,
and the second term goes to zero because
min
B

u
n

_
B

v dl
. .
0

_
B

vu
n
dl max
B

u
n

_
B

v dl
. .
0
,
which completes the proof. 2
The formula (38) represents the solution as the superposition of Greens func-
tions, which can be thought of as point source responses. The following result
tells us that the response at x to a point source located at x

is equal to the response


at x

to a point source located at x.


Theorem 16 (Reciprocity Principle) G(x, x

) = G(x

, x) for x = x

.
58
Proof. Let y and y

be distinct points in , let B

and B

denote -radius disks


centred at y and y

, and let u(x) = G(x, y

) and v(x) = G(x, y). Then


0 =
_
\B

\B

uv vudA
=
_

(uv vu) ndl


. .
0

_
B


_
B

u(y)
_
B

v ndl
. .
1
+
_
B

vu
n
dl
. .
0

_
B

uv
n
dl
. .
0
+v(y

)
_
B

u
n
dl
. .
1
,
leaving us with u(y) v(y

) = 0, that is, G(y, y

) G(y

, y) = 0. 2
As a consequence of Theorem 4, formula (38) can be written as
u(x) =
_

G(x, x

)f(x

) dA

.
The singularity function (37) is called the free-space Greens function for Pois-
sons equation. It is a Greens function for Poissons equation with the boundary
condition u(x)
1
|x|
as |x| . The free-space Greens function doesnt satisfy
this boundary condition, but the boundary condition does ensure that the term
_

(uv vu) ndl


in the proof of Theorem 3 goes to zero when is taken to be a disk of radius R
centred at x

and R .
The method of images can be used to nd Greens functions for other domains.
For example, using the idea that a unit point source located at x

= (x

, y

) with
y

> 0 and a point source of strength 1 located at (x

, y

) will cancel each other


on the x-axis, we nd the Greens function
G(x, x

) =
1
2
ln
_
(x x

)
2
+ (y y

)
2
_
1
2
+
1
2
ln
_
(x x

)
2
+ (y +y

)
2
_
1
2
for the Poisson equation in the half-plane y > 0.
Using similar ideas, one can derive the formula for Greens function for the
origin-centred disk of radius a as
G(x, x

) =
1
2
ln |x x

| +
1
2
ln

|x|
2
xa
2

a|x|
.
Representing x and x

in polar coordinates as (r, ) and (,

), and using the


cosine law |x x

|
2
= |x|
2
+ |x

|
2
2|x||x

| cos with =

, the disk
Greens function can be written
G(x, x

) =
1
2
ln

a
4
+r
2

2
2a
2
r cos
a
2
(r
2
+
2
2r cos )
. (39)
59
The Greens function G((x, y), (0, 0.5)) for the half-plane and the unit-radius
disk are shown below.
!1
0
1
0
0.5
1
2
0
0.1
0.2
0.3
x
y
G
!1
0
1
!1
0
1
0
0.1
0.2
0.3
0.4
x
y
G
The Greens function also provides the solution of Laplaces equation with
nonhomogeneous Dirichlet boundary conditions:
Theorem 17 If u = 0 in and u = h on then
u(x

) =
_

h(x)G(x, x

) ndl.
Proof. Similar to the proof of Theorem 4. 2
For the half-plane y > 0 the unit normal points in the negative y direction and
we have
G
n

=
G
y

y=0
=
y

/
(x x

)
2
+ (y

)
2
.
The solution of Laplaces equation u = 0 in the half-plane with u(x, 0) = h(x)
is therefore
u(x

, y

) =
1

h(x)
(x x

)
2
+ (y

)
2
dx.
For the origin-centred disk, the unit outward normal points in the radial direc-
tion and we have
G
n

=
G
r

r=a
=
1
2a

a
2

2
a
2
+
2
2a cos
.
The solution of Laplaces equation on the disk with u = h on the boundary r = a
is therefore
u(,

) =
1
2
_
2
0
(a
2

2
)h()
a
2
+
2
2a cos(

)
d,
which is Poissons integral formula (section 10.1).
11.3 Greens functions from eigenfunctions
The Poisson problem can also be solved by the method of eigenfunctions. To
introduce the technique, we start with a one-dimensional problem.
60
Consider the differential equation u

= f with boundary conditions u(0) =


0, u(a) = 0. The associated eigenvalue problem is

+ = 0 with the same


boundary conditions. The eigenvalues and eigenfunctions are

m
=
m
2

2
a
2
,
m
(x) = sin
mx
a
.
Substituting a trial solution of the form

m

1
A
m

m
(x) into the differential
equation, multiplying through by
m
(x), and integrating gives
A
m

m
_
a
0
sin
2
mx
a
dx
. .
a/2
=
_
a
0
f(x) sin
mx
a
dx,
Thus the A
m
are
1

m
the Fourier sine coefcients of f, and the solution at a point
(0, 1) is
u() =
_
a
0
f(x)
_

m1
2a
m
2

2
sin
mx
a
sin
m
a
_
dx.
Comparing this with formula (35), we deduce the Greens function to be the term
in brackets, that is,
G(x, ) =

m1
2a
m
2

2
sin
mx
a
sin
m
a
.
When a = 1, this is the Fourier sine expansion of the Greens function presented
in section 1.
Next, consider the two-dimensional Poisson equation u = f on the do-
main (0, a) (0, b), with u = 0 on the boundary. The associated eigenvalue
problem is u + u = 0, with the same boundary conditions. Assuming a solu-
tion of the form u(x, y) = X(x)Y (y) yields, with
2
as separation constant, the
two eigenvalue problems
X

+
2
X = 0, Y

+ (
2
)Y = 0
with homogeneous boundary conditions X(0) = X(a) = Y (0) = Y (b) = 0. The
eigenfunctions of the one-dimensional problems are
X
m
(x) = sin
mx
a
, Y
n
(y) = sin
ny
b
.
and the eigenvalues are

mn
=
2
_
m
2
a
2
+
n
2
b
2
_
.
Substituting a trial solution of the form

m

1
A
m

,n
X
m
(x)Y
n
(y) into
the Poisson equation, multiplying through by X
m
Y
n
, and integrating gives
A
mn

mn
_
a
0
sin
2
mx
a
dx
. .
a/2
_
b
0
sin
2
ny
b
dy
. .
b/2
=
_
b
0
_
a
0
f(x, y) sin
mx
a
sin
ny
b
dx dy
61
Thus the A
mn
are
1

mn
the two-dimensional Fourier sine coefcients of f, and
the solution at a point x

is
u(x

, y

) =
_
a
0
_
b
0
f(x, y)G(x, y, x

, y

) dx dy
with
G(x, y, x

, y

) =
4ab

m1

n1
1
m
2
b
2
+n
2
a
2
sin
mx
a
sin
ny
b
sin
mx

b
sin
ny

b
.
The Greens function G(x, y, a/2, a/4) for a square domain is plotted below.
0
0.5
1
0
0.25
1
0
0.1
0.2
0.3
x/a
y/a
G
62

You might also like