Chapter 1

Variational Methods

1.1 Stationary Values of Functions

Recall Taylor’s Theorem for a function f(x) in three dimensions with a displacement δx = (δx, δy, δz):

f(x + δx) = f(x) + \frac{∂f}{∂x} δx + \frac{∂f}{∂y} δy + \frac{∂f}{∂z} δz + higher order terms  (1.1)

so that

δf = f(x + δx) − f(x) = \frac{∂f}{∂x} δx + \frac{∂f}{∂y} δy + \frac{∂f}{∂z} δz + ⋯
= ∇f . δx + · · · .

In the limit |δx| → 0 we write

df = ∇f . dx.  (1.2)

This result is true in any number n of dimensions.

At an extremum (a maximum or minimum) f must be stationary, i.e. the first variation df must vanish for all possible directions of dx. This can only happen if ∇f = 0 there.

Note that if we try to find the extrema of f by solving ∇f = 0, we may also find other stationary points of f which are neither maxima nor minima, for instance saddle points. (This is the same difficulty as in one dimension, where a stationary point may be a point of inflection rather than a maximum or minimum.)

If we need to find the extrema of f in a bounded region – for instance, within a two-dimensional unit square – then not only must we solve ∇f = 0 but we must also compare the resulting values of f with those on the boundary of the square. It is quite possible for the maximum value to occur on the boundary without that point being a stationary one.

Example: Consider the Gaussian ‘bump’:

f(x, y, z) = exp(−x^{2} − y^{2} − z^{2}).

Then ∇f = −2(x, y, z)f. Setting this equal to zero gives (x, y, z) = (0, 0, 0), i.e. maximum at the origin.

Constrained stationary values

Suppose that we wish to find the extrema of f(x) subject to a constraint of the form g(x) = c, where c is some constant. In this case, the first variation df must still vanish, but now not all possible directions for dx are allowed: only those which lie in the surface defined by g(x) = c. Hence, since df = ∇f.dx, the vector ∇f must lie perpendicular to the surface.

But recall that the normal to a surface of the form g(x) = c is in the direction ∇g. Hence ∇f must be parallel to ∇g, i.e., ∇f = λ∇g for some scalar λ.

This gives us the method of Lagrange’s undetermined multiplier: solve the n equations

∇(f − λg) = 0  (1.3)

for x together with the single constraint equation

g(x) = c.  (1.4)

The resulting values of x give the stationary points of f subject to the constraint. Note that while solving the total of n+1 equations it is usually possible to eliminate λ without ever finding its value; hence the description “undetermined”.

If there are two constraints g(x) = c and h(x) = k, then we need a multiplier for each constraint, and we solve

∇(f − λg − μh) = 0  (1.5)

together with the two constraints. The extension to higher numbers of constraints is straightforward.

Worked Example 1: Constrained Maximization

A cuboid is inscribed in an ellipsoid with semi-axes a, b and c. What is its maximum volume?

We must find values of x, y and z which maximize the cuboid’s volume f(x, y, z) ≡ 8xyz subject to the constraint

g(x, y, z) ≡ \frac{x^{2}}{a^{2}} + \frac{y^{2}}{b^{2}} + \frac{z^{2}}{c^{2}} = 1

which ensures that the vertices of the cuboid lie on the surface of the ellipse. We introduce an undetermined multiplier λ and consider the three components of the equation ∇(f − λg) = 0:

8yz − \frac{2λx}{a^{2}} = 0,
8xz − \frac{2λy}{b^{2}} = 0,
8xy − \frac{2λz}{c^{2}} = 0.

Multiplying these equations by x, y and z respectively, we see that

8xyz = 2λ \frac{x^{2}}{a^{2}} = 2λ \frac{y^{2}}{b^{2}} = 2λ \frac{z^{2}}{c^{2}} .

Hence either λ = 0, or

\frac{x^{2}}{a^{2}} = \frac{y^{2}}{b^{2}} = \frac{z^{2}}{c^{2}} .

The first possibility would imply that the volume 8xyz is zero, which is clearly a minimum rather than the maximum which we seek, so we exclude this from now on. Remembering the constraint

\frac{x^{2}}{a^{2}} + \frac{y^{2}}{b^{2}} + \frac{z^{2}}{c^{2}} = 1

we conclude that

\frac{x^{2}}{a^{2}} = \frac{y^{2}}{b^{2}} = \frac{z^{2}}{c^{2}} = \frac{1}{3},

i.e., (x, y, z) = \frac{1}{\sqrt{3}}(a, b, c). The required maximum volume is therefore 8abc/3\sqrt{3}.

Note that we can also, if we wish, deduce that λ = 4abc/\sqrt{3}, but this is of no practical relevance.

1.2 Functionals

Now introduce functionals, real-valued mappings whose arguments are functions, i.e. F : one or more functions → ℝ.

Let y(x) be a function of x in some interval a < x < b, and consider the definite integral

F = ∫_{a}^{b}({y(x)}^{2} + y′(x)y′′(x)) dx.  (1.6)

F is clearly independent of x; instead it depends only on the function y(x). F is a simple example of a functional, and to show the dependence on y we normally denote it F [y].

We can also think of functionals as an extension of the concept of a function of many variables – e.g. g(x_{1}, x_{2}, …, x_{n}), a function of n variables – to a function of an infinite number of variables, because F depends on every single value that y takes in the range a < x < b.

We shall be concerned in this chapter with functionals of the form

F [y] = ∫_{a}^{b} f(x, y, y′) dx  (1.7)

where f depends only on x and the value of y and its first derivative at x. However, the theory can be extended to more general functionals (for example, with functions f(x, y, y′, y′′, y′′′, …) which depend on higher derivatives, or double integrals with two independent variables x_{1} and x_{2} instead of just x).

1.3 Variational Principles

Functionals are useful because many laws of physics and of physical chemistry can be recast as statements that some functional F [y] is minimised.

For example, a heavy chain suspended between two fixed points hangs in equilibrium in such a way that its total gravitational potential energy (which can be expressed as a functional) is minimised. A mechanical system of heavy elastic strings minimises the total potential energy, both elastic and gravitational. Similar principles apply when electric fields and charged particles are present (we include the electrostatic potential energy) and when chemical reactions take place (we include the chemical potential energy).

Two fundamental examples of such variational principles are due to Fermat and Hamilton.

Fermat’s Principle

Consider a light ray passing through a medium of variable refractive index μ(r). The path it takes between two fixed points A and B is such as to minimise the optical path length

∫_{A}^{B} μ(r) dl,  (1.8)

where dl is the length of a path element.

Strictly speaking, Fermat’s principle only applies in the geometrical optics approximation; i.e., when the wavelength of the light is small compared with the physical dimensions of the optical system, so that light may be regarded as rays. This is true for a telescope, but not for Young’s slits: when the geometrical optics approximation fails to hold, diffraction occurs.

For example, consider air above a hot surface, say a tarmac road on a hot day. The air is hotter near the road and cooler above, so that μ is smaller closer to the road surface. A light ray travelling from a car to an observer minimises the optical path length by staying close to the road, and so bends appropriately. The light seems to the observer to come from a low angle, leading to a virtual image (and hence to the “mirage” effect).

Hamilton’s Principle of Least Action

Consider a mechanical system with kinetic energy T and potential energy V which is in some given configuration at time t_{1} and some other configuration at time t_{2}. Define the Lagrangian of the system by

ℒ = T − V,  (1.9)

NST IB Mathematical Methods II 6 and define the action to be

𝒮 = ∫_{t_{1}}^{t_{2}} ℒ dt  (1.10)

(a functional which depends on the way the system moves). Hamilton’s principle states that the actual motion of the system is such as to minimise the action.

1.4 The Calculus of Variations

How do we find the function y(x) which minimises, or more generally makes stationary, our archetypal functional

F [y] = ∫_{a}^{b} f(x, y, y′) dx,  (1.11)

with fixed values of y at the end-points (viz. fixed y(a) and y(b))?

We consider changing y to some “nearby” function y(x) + δy(x), and calculate the corresponding change δF in F (to first order in δy). Then F is stationary when δF = 0 for all possible small variations δy.

Note that a more “natural” notation would be to write dF rather than δF, since we will consider only the first-order change and ignore terms which are second order in δy. However, the notation δ is traditional in this context.

Now

δF = F [y + δy] − F [y]

= ∫_{a}^{b} f(x, y + δy, y′ + (δy)′) dx − ∫_{a}^{b} f(x, y, y′) dx

= ∫_{a}^{b} \{f(x, y, y′) + \frac{∂f}{∂y} δy + \frac{∂f}{∂y′} (δy)′\} dx − ∫_{a}^{b} f(x, y, y′) dx

[using a Taylor expansion to first order]

= ∫_{a}^{b} \{\frac{∂f}{∂y} δy + \frac{∂f}{∂y′} (δy)′\} dx

= [\frac{∂f}{∂y′} δy]_{a}^{b} + ∫_{a}^{b} \{\frac{∂f}{∂y} δy − \frac{d}{dx} (\frac{∂f}{∂y′} ) δy\} dx

[integrating by parts]

= ∫_{a}^{b} \{\frac{∂f}{∂y} − \frac{d}{dx} (\frac{∂f}{∂y′} )\} δy dx

since δy = 0 at x = a, b (because y(x) is fixed there). It is clear that δF = 0 for all possible small variations δy(x) if and only if

\frac{d}{dx} (\frac{∂f}{∂y′} ) = \frac{∂f}{∂y} .  (1.12)

This is Euler’s equation.

Notation

∂f/∂y′ means “differentiate with respect to y′, keeping x and y constant”, and it seems impossible for y′ to change if y does not. But ∂/∂y and ∂/∂y′ in Euler’s equation are just formal derivatives (as though y and y′ were unconnected) and in practice it is easy to do straightforward “ordinary” partial differentiation.

Example: if f(x, y, y′) = x(y′^{2} − y^{2}) then

\frac{∂f}{∂y} = −2xy,  \frac{∂f}{∂y′} = 2xy′,  (1.13)

\frac{∂f}{∂x} = (y^{′2} − y^{2}).

Note however that d/dx and ∂/∂x mean very different things: ∂/∂x means “keep y and y′ constant” whereas d/dx is a so-called “full derivative”, so that y and y′ are differentiated with respect to x as well.

Continuing with the above example (1.13),

\frac{∂}{∂x} (\frac{∂f}{∂y′} ) = 2y′,  (1.14)

but

\frac{d}{dx} (\frac{∂f}{∂y′} ) = \frac{d}{dx}(2xy′) = 2y′ + 2xy′′.  (1.15)

Hence Euler’s equation for this example is

2y′ + 2xy′′ = −2xy  (1.16)

or

y′′ + \frac{1}{x} y′ + y = 0  (1.17)

(Bessel’s equation of order 0).

Several Dependent Variables

Suppose, instead of just one dependent variable y(x), we have n dependent variables y_{1}(x), y_{2}(x), …, y_{n}(x), so that our functional is

F [y_{1}, …, y_{n}] = ∫_{a}^{b} f(x, y_{1}, …, y_{n}, y′_{1}, …, y′_{n}) dx?  (1.18)

NST IB Mathematical Methods II 9 In this case, Euler’s equation applies to each y_{i}(x) independently, so that

\frac{d}{dx} (\frac{∂f}{∂y′_{i}} ) = \frac{∂f}{∂y_{i}}  (1.19)

for i = 1, …, n.

The proof is very similar to before:

δF = ∫_{a}^{b} \{\frac{∂f}{∂y_{1}} δy_{1} + ⋯ + \frac{∂f}{∂y_{n}} δy_{n} + \frac{∂f}{∂y′_{1}} (δy_{1})′ + ⋯ + \frac{∂f}{∂y′_{n}} (δy_{n})′\} dx

= ∫_{a}^{b} ∑_{i=1}^{n} \{\frac{∂f}{∂y_{i}} δy_{i} + \frac{∂f}{∂y′_{i}} (δy_{i})′\} dx

= ∑_{i=1}^{n} ∫_{a}^{b} \{\frac{∂f}{∂y_{i}} − \frac{d}{dx} (\frac{∂f}{∂y′_{i}} )\} δy_{i} dx

using the same manipulations (Taylor expansion and integration by parts). It is now clear that we can only have δF = 0 for all possible variations of all the y_{i}(x) if Euler’s equation applies simultaneously to each and every one of the y_{i}.

1.5 A First Integral

In some cases, it is possible to find a first integral (i.e., a constant of motion) of Euler’s equation. Consider

\frac{df}{dx} = \frac{∂f}{∂x} + y′ \frac{∂f}{∂y} + y′′ \frac{∂f}{∂y′}  (1.20)

(calculating \frac{d}{dx}f(x, y(x), y′(x)) using the chain rule). Using Euler’s equation,

\frac{df}{dx} = \frac{∂f}{∂x} + y′ \frac{d}{dx} (\frac{∂f}{∂y′} ) + y′′ \frac{∂f}{∂y′}

= \frac{∂f}{∂x} + \frac{d}{dx} (y′ \frac{∂f}{∂y′} )

[product rule]

so that

\frac{d}{dx} (f − y′ \frac{∂f}{∂y′} ) = \frac{∂f}{∂x}.  (1.21)

Now, if f has no explicit x-dependence, so that ∂f/∂x = 0, we immediately deduce that

f − y′ \frac{∂f}{∂y′} = constant.

(1.22)

(Note that “f has no explicit x-dependence” means that x does not itself appear in the expression for f, even though y and y′ implicitly depend on x; so f = y′^{2} − y^{2} has no explicit x-dependence while f = x(y′^{2} − y^{2}) does.)

If there are n dependent variables y_{1}(x), …, y_{n}(x), then the first integral above is easily generalised to

f − ∑_{i=1}^{n} y′_{i} \frac{∂f}{∂y′_{i}} = constant  (1.23)

if f has no explicit x-dependence.

1.6 Applications of Euler’s Equation

Geodesics

A geodesic is the shortest path on a given surface between two specified points A and B. We will illustrate the use of Euler’s equation with a trivial example: geodesics on the Euclidean plane.

The total length of a path from (x_{1}, y_{1}) to (x_{2}, y_{2}) along the path y(x) is given by

L = ∫_{A}^{B} dl = ∫_{A}^{B} \sqrt{dx^{2} + dy^{2}}

= ∫_{A}^{B} \sqrt{1 + (\frac{dy}{dx})^{2}} dx = ∫_{x_{1}}^{x_{2}} \sqrt{1 + y′^{2}} dx.

Note that we assume that y(x) is single-valued, i.e., the path does not curve back on itself.

We wish to minimise L over all possible paths y(x) with the end-points held fixed, so that y(x_{1}) = y_{1} and y(x_{2}) = y_{2} for all paths. This is precisely our archetypal variational problem with

f(x, y, y′) = \sqrt{1 + y′^{2}},  (1.24)

NST IB Mathematical Methods II 11

and hence

\frac{∂f}{∂y} = 0,  \frac{∂f}{∂y′} = \frac{y′}{\sqrt{1 + y′^{2}}} .  (1.25)

The Euler equation is therefore

\frac{d}{dx} (\frac{y′}{\sqrt{1 + y′^{2}}} ) = 0  ⟹  \frac{y′}{\sqrt{1 + y′^{2}}} = k,  a constant.  (1.26)

So y′^{2} = k^{2}/(1 − k^{2}). It is clear that k ≠ ±1, so y′ is a constant, m say. Hence the solutions of Euler’s equation are the functions

y = mx + c  (1.27)

(where m and c are constants) – i.e., straight lines! To find the particular values of m and c required in this case we now substitute in the boundary conditions y(x_{1}) = y_{1}, y(x_{2}) = y_{2}.

It is important to note two similarities with the technique of minimising a function f(x) by solving ∇f = 0.

Firstly, we have not shown that this straight line does indeed produce a minimum of L: we have shown only that L is stationary for this choice, so it might be a maximum or even some kind of “point of inflection”. It is usually easy to confirm that we have the correct solution by inspection – in this case it is obviously a minimum. (There is no equivalent of the one-dimensional test f′′(x) > 0 for functionals, or at least not one which is simple enough to be of any use.)

Secondly, assuming that we have indeed found a minimum, we have shown only that it is a local minimum, not a global one. That is, we have shown only that “nearby” paths have greater length. Once again, however, we usually confirm that we have the correct solution by inspection. Compare this difficulty with the equivalent problem for functions, illustrated by the graph below.

An alternative method of solution for this simple geodesic problem is to note that f(x, y, y′) = \sqrt{1 + y′^{2}} has no explicit xdependence, so we can use the first integral:

const. = f − y′ \frac{∂f}{∂y′} = \sqrt{1 + y′^{2}} − y′ \frac{y′}{\sqrt{1 + y′^{2}}}

= \frac{1}{\sqrt{1 + y′^{2}}},

i.e., y′ is constant (as before).

Worked Example 2: Geodesics on the Surface of a Sphere

Recall that in orthogonal curvilinear coordinates (q_{1}, q_{2}, q_{3}),

dr = h_{1} dq_{1} e_{1} + h_{2} dq_{2} e_{2} + h_{3} dq_{3} e_{3}.

In spherical polar coordinates,

dr = dr e_{r} + r dθ e_{θ} + r sin θ dϕ e_{ϕ}.

Without loss of generality, we may take the sphere to be of unit radius: the length of a path from A to B is then

L  = ∫_{A}^{B} |dr|
= ∫_{A}^{B} \sqrt{dθ^{2} + sin^{2} θ dϕ^{2}} [since dr = 0]
= ∫_{θ_{A}}^{θ_{B}} \sqrt{1 + sin^{2} θ ϕ^{′2}} dθ

where the path is described by the function ϕ(θ). Using Euler’s equation,

\frac{d}{dθ} (\frac{∂}{∂ϕ′} \sqrt{1 + sin^{2} θ ϕ^{′2}}) = \frac{∂}{∂ϕ}\sqrt{1 + sin^{2} θ ϕ^{′2}} = 0

so that

\frac{sin^{2} θ ϕ′}{\sqrt{1 + sin^{2} θ ϕ^{′2}}}

is a constant, c say. Hence

ϕ′ = \frac{c}{sin θ\sqrt{sin^{2} θ − c^{2}}}

and the problem reduces to integrating this with respect to θ.

Substitute u = cot θ so that du = − cosec^{2} θ dθ. Then

ϕ = ∫ \frac{−c du}{\sqrt{1−c^{2} cosec^{2} θ}}
= ∫ \frac{−c du}{\sqrt{1−c^{2}(1+u^{2})}}
= ∫ \frac{−du}{\sqrt{a^{2}−u^{2}}}  where a = \frac{\sqrt{1 − c^{2}}}{c}
= cos^{−1}(u/a) + ϕ_{0}

where ϕ_{0} is a constant of integration. Hence the geodesic path is given by

cot θ = a cos(ϕ − ϕ_{0})

and the arbitrary constants a and ϕ_{0} must be found using the end-points. This is a great circle path.

The Brachistochrone

A bead slides down a frictionless wire, starting from rest at a point A. What shape must the wire have for the bead to reach some lower point B in the shortest time? (A similar device was used in some early clock mechanisms.)

Using conservation of energy, \frac{1}{2}mv^{2} = mgy, i.e., v = \sqrt{2gy}. Also dl = v dt, so

dt = \begin{matrix}√ & \overline{dx^{2} + dy^{2}} \\ √ & \overline{2gy} \end{matrix} = \frac{1}{\sqrt{2g}} \begin{matrix}√ & \overline{1 + y′^{2}} \\ √ & \overline{y} \end{matrix} dx.  (1.28)

The time taken to reach B is therefore

T [y] = \frac{1}{\sqrt{2g}} ∫_{0}^{x_{B}} \sqrt{\frac{1 + y′^{2}}{y}} dx  (1.29)

and we wish to minimise this, subject to y(0) = 0, y(x_{B}) = y_{B}. We note that the integrand has no explicit x-dependence, so we use the first integral

const. = \sqrt{\frac{1 + y′^{2}}{y}} − y′ \frac{∂}{∂y′} \sqrt{\frac{1 + y′^{2}}{y}}

= \sqrt{\frac{1 + y′^{2}}{y}} − \frac{y′^{2}}{\sqrt{y}\sqrt{1 + y′^{2}}}

= \frac{1}{\sqrt{y}\sqrt{1 + y′^{2}}} .

Hence y(1 + y′^{2}) = c, say, a constant, so that

y′ = \sqrt{\frac{c − y}{y}}  or \sqrt{\frac{y}{c − y}} dy = dx.  (1.30)

Substitute y = c sin^{2} θ; then

dx = 2c \sqrt{\frac{sin^{2} θ}{1 − sin^{2} θ}} sin θ cos θ dθ
= 2c sin^{2} θ dθ
= c(1 − cos 2θ) dθ.

Using the initial condition that when y = 0 (i.e., θ = 0), x = 0, we obtain

x = c(θ − \frac{1}{2} sin 2θ),
y = c sin^{2} θ

which is an inverted cycloid. The constant c is found by applying the other condition, y = y_{B} when x = x_{B}.

Note that strictly speaking we should have said that y′ = ±\sqrt{(c − y)/y} above. Taking the negative root instead of the positive one would have lead to

x = −c(θ − \frac{1}{2} sin 2θ),
y = c sin^{2} θ,

which is exactly the same curve but parameterised in the opposite direction.

Light and Sound

Consider light rays travelling through a medium with refractive index inversely proportional to \sqrt{z} where z is the height. By Fermat’s principle, we must minimise

∫ \frac{dl}{\sqrt{z}} .  (1.31)

This is exactly the same variational problem as for the Brachistochrone, so we conclude that light rays will follow the path of a cycloid.

More realistically, consider sound waves in air. Sound waves obey a principle similar to Fermat’s: except at very long wavelengths, they travel in such a way as to minimise the time taken to travel from A to B,

∫_{A}^{B} \frac{dl}{c},  (1.32)

where c is the (variable) speed of sound (comparable to 1/μ for light). Consider a situation where the absolute temperature T of the air is linearly related to the height z, so that T = αz + T_{0} for some temperature T_{0} at ground level. Since the speed of sound is proportional to the square root of the absolute temperature, we have c ∝ \sqrt{αz + T_{0}} = \sqrt{Z} say. This leads once again to the Brachistochrone problem (for Z rather than z), and we conclude that sound waves follow paths z(x) which are parts of cycloids, scaled vertically by a factor 1/α (check this as an exercise).

1.7 Hamilton’s Principle in Mechanical Problems

Hamilton’s principle can be used to solve many complicated problems in rigid-body mechanics. Consider a mechanical system whose configuration can be described by a number of so-called generalised coordinates q_{1}, q_{2}, …, q_{n}. Examples:

• A particle with position vector r = (x_{1}, x_{2}, x_{3}) moving through space. Here we can simply let q_{1} = x_{1}, q_{2} = x_{2} and q_{3} = x_{3}: there are three generalised coordinates.

• A pendulum swinging in a vertical plane: here there is only one generalised coordinate, q_{1} = θ, the angle to the vertical.

• A rigid body (say a top) spinning on its axis on a smooth plane. This requires five generalised coordinates: two to describe the position of the point of contact on the plane, one for the angle of the axis to the vertical, one for the rotation of the axis about the vertical, and one for the rotation of the top about its own axis.

The Lagrangian ℒ = T − V is a function of t, q_{1}, …, q_{n} and q̇_{1}, …, q̇_{n}, so

𝒮 = ∫ ℒ (t, q_{1}(t), …, q_{n}(t), q̇_{1}(t), …, q̇_{n}(t)) dt.  (1.33)

This is a functional with n dependent variables q_{i}(t), so we can use Euler’s equation (with t playing the role of x, and q_{i}(t) playing the role of y_{i}(x)) for each of the q_{i} independently:

\frac{d}{dt} (\frac{∂ℒ}{∂q̇_{i}} ) = \frac{∂ℒ}{∂q_{i}}  (1.34)

for each i. In this context these equations are known as the Euler–Lagrange equations.

In the case when ℒ has no explicit time-dependence, the first integral (from §1.5) gives us that

ℒ − ∑_{i=1}^{n} q̇_{i} \frac{∂ℒ}{∂q̇_{i}} = constant. (1.35)

It is frequently the case that T is a homogeneous quadratic in the q̇_{i}, i.e., it is of the form

∑_{i=1}^{n} ∑_{j=1}^{n} a_{ij}(q_{1}, …, q_{n}) q̇_{i}q̇_{j}  (1.36)

where the coefficients a_{ij} do not depend on any of the “generalised velocities” q̇_{i} or on t, and V also does not depend on the velocities or time so that V = V (q_{1}, …, q_{n}). Then it can be shown that

ℒ − ∑_{i=1}^{n} q̇_{i} \frac{∂ℒ}{∂q̇_{i}} = (T − V) − 2T = −(T + V),  (1.37)

i.e., the total energy E = T + V is conserved when there is no explicit time-dependence. This fails however when the external forces vary with time or when the potential is velocity-dependent, e.g., for motion in a magnetic field.

A Particle in a Conservative Force Field

Here

ℒ = \frac{1}{2}m(ẋ_{1}^{2} + ẋ_{2}^{2} + ẋ_{3}^{2}) − V (x_{1}, x_{2}, x_{3});  (1.38)

hence the Euler–Lagrange equations are

\frac{d}{dt}(mẋ_{1}) = −\frac{∂V}{∂x_{1}}, \frac{d}{dt}(mẋ_{2}) = −\frac{∂V}{∂x_{2}}, \frac{d}{dt}(mẋ_{3}) = −\frac{∂V}{∂x_{3}},

or in vector notation

\frac{d}{dt}(mṙ) = −∇V,  (1.39)

i.e., F = ma where F = −∇V is the force and a = r̈ is the acceleration.

Two Interacting Particles

Consider a Lagrangian

ℒ = \frac{1}{2}m_{1}|ṙ_{1}|^{2} + \frac{1}{2}m_{2}|ṙ_{2}|^{2} − V (r_{1} − r_{2}),  (1.40)

where the only force is a conservative one between two particles with masses m_{1} and m_{2} at r_{1} and r_{2} respectively, and depends only on their (vector) separation.

We could use the six Cartesian coordinates of the particles as generalised coordinates; but instead define

r = r_{1} − r_{2},  (1.41)

the relative position vector, and

R = \frac{m_{1}r_{1} + m_{2}r_{2}}{M} ,  (1.42)

the position vector of the centre of mass, where M = m_{1} + m_{2} is the total mass. Now

|ṙ_{1}|^{2} = |Ṙ + \frac{m_{2}}{M} ṙ|^{2} = (Ṙ + \frac{m_{2}}{M} ṙ). (Ṙ + \frac{m_{2}}{M} ṙ) = |Ṙ|^{2} + \frac{m_{2}^{2}}{M^{2}} |ṙ|^{2} + \frac{2m_{2}}{M} Ṙ . ṙ

and similarly

|ṙ_{2}|^{2} = |Ṙ|^{2} + \frac{m_{1}^{2}}{M^{2}} |ṙ|^{2} − \frac{2m_{1}}{M} Ṙ . ṙ.

Let r = (x_{1}, x_{2}, x_{3}), R = (X_{1}, X_{2}, X_{3}), and use these as generalised coordinates. Then

ℒ = \frac{1}{2}M|Ṙ|^{2} + \frac{m_{1}m_{2}}{2M} |ṙ|^{2} − V (r)

= \frac{1}{2}M(Ẋ_{1}^{2} + Ẋ_{2}^{2} + Ẋ_{3}^{2}) + \frac{m_{1}m_{2}}{2M} (ẋ_{1}^{2} + ẋ_{2}^{2} + ẋ_{3}^{2}) − V (x_{1}, x_{2}, x_{3}).

The Euler–Lagrange equation for X_{i} is therefore

\frac{d}{dt}(MẊ_{i}) = 0,  (1.43)

i.e., R̈ = 0 (the centre of mass moves with constant velocity); and for x_{i} is

\frac{d}{dt} (\frac{m_{1}m_{2}}{M} ẋ_{i}) = −\frac{∂V}{∂x_{i}} ,  (1.44)

i.e., μr̈ = −∇V where μ is the reduced mass m_{1}m_{2}/(m_{1} + m_{2}) (the relative position vector behaves like a particle of mass μ).

Note that the kinetic energy T is a homogeneous quadratic in the Ẋ_{i} and ẋ_{i}; that V does not depend on the Ẋ_{i} and ẋ_{i}; and that ℒ has no explicit t-dependence. We can deduce immediately that the total energy E = T + V is conserved.

1.8 The Calculus of Variations with Constraint

In §1.1 we studied constrained variation of functions of several variables. The extension of this method to functionals (i.e., functions of an infinite number of variables) is straightforward: to find the stationary values of a functional F [y] subject to G[y] = c, we instead find the stationary values of F [y] − λG[y], i.e., find the function y which solves δ(F − λG) = 0, and then eliminate λ using G[y] = c.

Worked Example 3: The Catenary

Consider a uniform chain of length L, with mass per unit length ρ, hanging under gravity between the points (−1, 1) and (1, 1). It adopts a form of minimum potential energy, that is it minimises

∫_{−1}^{1} ρgy dl ∝ F [y] ≡ ∫_{−1}^{1} y\sqrt{1 + y^{′2}} dx

subject to the prescribed length,

L = G[y] ≡ ∫_{−1}^{1} \sqrt{1 + y^{′2}} dx.

This is equivalent to minimising F − λG, i.e., to solving

δ ∫_{−1}^{1}(y − λ)\sqrt{1 + y^{′2}} dx = 0.

The integrand has no explicit x-dependence, so we use the first integral

c = (y − λ)\sqrt{1 + y^{′2}} − y′(y − λ)\frac{y′}{\sqrt{1+y^{′2}}}
= \frac{y−λ}{\sqrt{1+y^{′2}}} ,

where c is a constant, whence

x = ∫ \frac{c dy}{\sqrt{(y − λ)^{2} − c^{2}}} .

Making the substitution y = λ + c cosh θ we obtain

x = c cosh^{−1} (\frac{y − λ}{c}) + x_{0}

where x_{0} is an arbitrary constant of integration. Hence the solution is

y = λ + c cosh (\frac{x − x_{0}}{c} ) ,

which is a catenary.

We have three unknown constants, to be found using the equation for y at each of the two end-points, together with the constraint equation. We immediately obtain x_{0} = 0 by symmetry (or by solving the end-point equations for x_{0}). Now y′ = sinh(x/c) and hence \sqrt{1 + y^{′2}} = cosh(x/c); so

L = ∫_{−1}^{1} cosh \frac{x}{c} dx
= 2c sinh \frac{1}{c}.

This equation must, in general, be solved numerically for c given L; then λ can be found using the end-point at (1, 1),

1 = λ + c cosh \frac{1}{c}.

This completes the solution.

1.9 The Variational Principle for Sturm–Liouville Equations

We will now show how variational calculus can be used in the solution of second order linear ordinary differential equations of fairly general type. Under appropriate conditions, any second order ODE can be written as a Sturm-Liouville problem with weight function w. We shall show in this section that the following three problems are equivalent:

(i) Find the eigenvalues λ and eigenfunctions y(x) which solve the Sturm–Liouville problem

−\frac{d}{dx}(p(x)y′) + q(x)y = λw(x)y  (1.45)

in a < x < b, where neither p nor w vanish in the interval.

(ii) Find the functions y(x) for which

F [y] = ∫_{a}^{b}(py′^{2} + qy^{2}) dx  (1.46)

is stationary subject to G[y] = 1 where

G[y] = ∫_{a}^{b} wy^{2} dx.  (1.47)

The eigenvalues of the equivalent Sturm–Liouville problem in (i) are then given by the values of F [y].

(iii) Find the functions y(x) for which

Λ[y] = \frac{F [y]}{G[y]}  (1.48)

is stationary; the eigenvalues of the equivalent Sturm–Liouville problem are then given by the values of Λ[y].

Hence Sturm–Liouville problems can be reformulated as variational problems.

Note the similarity between (iii) and the stationary property of the eigenvalues of a symmetric matrix (recall that it is possible to find the eigenvalues of a symmetric matrix A by finding the stationary values of a^{T}Aa/a^{T}a over all possible vectors a). The two facts are in fact closely related.

To show that (ii) is equivalent to (i), consider

δ(F − λG) = δ ∫_{a}^{b}(py′^{2} + qy^{2} − λwy^{2}) dx.  (1.49)

Using Euler’s equation, F − λG is stationary when

\frac{d}{dx}(2py′) = 2qy − 2λwy,  (1.50)

i.e.,

−\frac{d}{dx}(py′) + qy = λwy,  (1.51)

which is the required Sturm–Liouville problem: note that the Lagrange multiplier of the variational problem is the same as the eigenvalue of the Sturm–Liouville problem.

Furthermore, multiplying the Sturm–Liouville equation by y and integrating, we obtain

∫_{a}^{b}(−y \frac{d}{dx}(py′) + qy^{2}) dx = λ ∫_{a}^{b} wy^{2} dx = λG[y] = λ  (1.52)

using the constraint. Hence

λ = ∫_{a}^{b}(−y \frac{d}{dx}(py′) + qy^{2}) dx
= [−ypy′]_{a}^{b} + ∫_{a}^{b}(py′^{2} + qy^{2}) dx

[integrating by parts]

= ∫_{a}^{b}(py′^{2} + qy^{2}) dx = F [y],

using “appropriate” boundary conditions. This proves that the stationary values of F [y] give the eigenvalues.

There are two ways of showing that (ii) is equivalent to (iii).

The first, informal way is to note that multiplying y by some constant α say does not in fact change the value of Λ[y]. This implies that when finding the stationary values of Λ we can choose to normalise y so that G[y] = 1, in which case Λ is just equal to F [y]. So finding the stationary values of Λ is equivalent to finding the stationary values of F subject to G = 1.

The second, formal way is to calculate

δΛ = \frac{F + δF}{G + δG} − \frac{F}{G}

= \frac{F + δF}{G} (1 − \frac{δG}{G}) − \frac{F}{G}

[using a Taylor expansion for (1 + δG/G)^{−1} to first order]

= \frac{δF}{G} − \frac{F δG}{G^{2}}

(again to first order). Hence δΛ = 0 if and only if δF = (F/G) δG; that is, Λ is stationary if and only if

δF − Λ δG = 0.  (1.53)

But this is just the same problem as (ii); so finding the stationary values of Λ is the same as finding the stationary values of F subject to G = 1.

In the usual case that p(x), q(x) and w(x) are all positive, we have that Λ[y] ≥ 0. Hence all the eigenvalues must be non-negative, and there must be a smallest eigenvalue λ_{0}; Λ takes the value λ_{0} when y = y_{0}, the corresponding eigenfunction. But what is the absolute minimum value of Λ over all functions y(x)? If it were some value μ < λ_{0}, then μ would be a stationary (minimum) value of Λ and would therefore be an eigenvalue, contradicting the statement that λ_{0} is the smallest eigenvalue. Hence Λ[y] ≥ λ_{0} for any function y(x).

As an example, consider the simple harmonic oscillator

−y′′ + x^{2}y = λy  (1.54)

subject to y → 0 as |x| → ∞. This is an important example as it is a good model for many physical oscillating systems. For instance, the Schrödinger equation for a diatomic molecule has approximately this form, where λ is proportional to the quantum mechanical energy level E; we would like to know the ground state energy, i.e., the eigenfunction with the lowest eigenvalue λ.

Here p(x) = 1, q(x) = x^{2} and w(x) = 1, so

Λ[y] = \frac{∫_{−∞}^{∞}(y′^{2} + x^{2}y^{2}) dx}{∫_{−∞}^{∞} y^{2} dx}.  (1.55)

We can solve this Sturm–Liouville problem exactly: the lowest eigenvalue turns out to be λ_{0} = 1 with corresponding eigenfunction y_{0} = exp(−\frac{1}{2}x^{2}). But suppose instead that we didn’t know this; we can use the above facts about Λ to try to guess at the value of λ_{0}. Let us use a trial function

y_{trial} = exp(−\frac{1}{2}αx^{2}),  (1.56)

where α is a positive constant (in order to satisfy the boundary conditions). Then

Λ[y_{trial}] = \frac{(α^{2} + 1) ∫_{−∞}^{∞} x^{2} exp(−αx^{2}) dx}{∫_{−∞}^{∞} exp(−αx^{2}) dx}.  (1.57)

We recall that ∫_{−∞}^{∞} exp(−αx^{2}) dx = \sqrt{π/α} and ∫_{−∞}^{∞} x^{2} exp(−αx^{2}) dx = \frac{1}{2} \sqrt{π/α^{3}} (by integration by parts). Hence Λ[y_{trial}] = (α^{2} + 1)/2α.

We know that Λ[y_{trial}], for any α, cannot be less than λ_{0}. The smallest value of (α^{2} + 1)/2α is 1, when α = 1; we conclude that λ_{0} ≤ 1, which gives us an upper bound on the lowest eigenvalue.

In fact this method has given us the exact eigenvalue and eigenfunction; but that is an accident caused by the fact that this is a particularly simple example!

The Rayleigh–Ritz Method

The Rayleigh–Ritz method is a systematic way of estimating the eigenvalues, and in particular the lowest eigenvalue, of a Sturm–Liouville problem. The first step is to reformulate the problem as the variational principle that Λ[y], the Rayleigh quotient, is stationary. Secondly, using whatever clues are available (for example, symmetry considerations or general theorems such as “the ground state wavefunction has no nodes”) we make an “educated guess” y_{trial}(x) at the true eigenfunction y_{0}(x) with lowest eigenvalue λ_{0}. It is preferable for y_{trial} to contain a number of adjustable parameters (e.g., α in the example above).

We can now find Λ[y_{trial}], which will depend on these adjustable parameters. We calculate the minimum value Λ_{min} of Λ with respect to all the adjustable parameters; we can then state that the lowest eigenvalue λ_{0} ≤ Λ_{min}. If the trial function was a reasonable guess then Λ_{min} should actually be a good approximation to λ_{0}.

If we wish, we can improve the approximation by introducing more adjustable parameters.

The fact that Λ[y] is stationary with respect to variations in the function y means that if y_{trial} is close to the true eigenfunction y_{0} (say within O(ε) of it) then the final calculated value Λ_{min} will be a very good approximation to λ_{0} (within O(ε^{2}) in fact). If the inclusion of further adjustable parameters fails to significantly improve the approximation then we can be reasonably sure that the approximation is a good one.

Note that if the trial function happens to include the exact solution y_{0} as a special case of the adjustable parameters, then the Rayleigh–Ritz method will find both y_{0} and λ_{0} exactly. This is what happened in the example above.

An alternative to calculating Λ[y_{trial}] and minimising it with respect to the adjustable parameters is to calculate F [y_{trial}] and G[y_{trial}], and to minimise F subject to G = 1. These procedures are equivalent, as we showed at the start of this section.

Worked Example 4: The Rayleigh–Ritz Method

The oscillations of a drum (e.g., a timpani, or more generally any circular membrane under tension and fixed at its boundary) obey Bessel’s equation of order zero,

y′′ + \frac{1}{x} y′ + λy = 0,

in 0 ≤ x ≤ 1, with boundary conditions that y should be non-singular at x = 0 and that y(1) = 0. Here λ = ω^{2}/c^{2} where ω is the frequency of oscillation of the drum and c is the wave speed. (This equation may be derived by converting the two-dimensional wave equation on the surface of the drum into plane polar coordinates and assuming a radially symmetric solution with a fixed frequency ω of oscillation.)

A drum can oscillate at many different frequencies, corresponding to the different eigenvalues of this Sturm–Liouville problem; but the fundamental (i.e., lowest) frequency is of the greatest interest since this is the one a listener will hear. (When a drum is struck, all of the possible frequencies are produced to varying extents, but the harmonics, i.e., the higher frequencies, usually decay rapidly leaving only the fundamental.) It is therefore natural to use the Rayleigh–Ritz method to estimate the lowest eigenvalue of Bessel’s equation (and thereby estimate the fundamental frequency).

Before we can proceed, we must put the equation into standard Sturm–Liouville self-adjoint form. By inspection we see that the appropriate equation is

−\frac{d}{dx} (x\frac{dy}{dx}) = λxy.

The equivalent variational problem is therefore that F [y] = ∫_{0}^{1} xy^{′2} dx is stationary subject to G[y] = ∫_{0}^{1} xy^{2} dx = 1. We shall use a trial solution of the form

y_{trial} = a + bx^{2} + cx^{4}

(chosen because we anticipate that the lowest eigenvalue corresponds to a solution which is even in x). This trial solution trivially satisfies the boundary condition at x = 0, and satisfies the condition at x = 1 so long as a + b + c = 0.

We now calculate

F [y_{trial}] = ∫_{0}^{1} xy_{trial}^{′2} dx = ∫_{0}^{1} x(2bx + 4cx^{3})^{2} dx
= b^{2} + \frac{8}{3}bc + 2c^{2}

and

G[y_{trial}] = ∫_{0}^{1} xy_{trial}^{2} dx  = ∫_{0}^{1} x(a + bx^{2} + cx^{4})^{2} dx
= \frac{1}{2}a^{2} + \frac{1}{2}ab + \frac{1}{6}(b^{2} + 2ac) + \frac{1}{4}bc + \frac{1}{10}c^{2}
= \frac{1}{6}b^{2} + \frac{5}{12}bc + \frac{4}{15}c^{2},

using a = −b − c.

We must either minimise F/G – which turns out to be rather messy algebraically – or minimise F subject to G = 1. We choose the latter; hence we minimise F − λG with respect to both b and c. So

0 = \frac{∂}{∂b}(F − λG) = (2 − \frac{1}{3}λ)b + (\frac{8}{3} − \frac{5}{12}λ)c (1)

and

0 = \frac{∂}{∂c}(F − λG) = (\frac{8}{3} − \frac{5}{12}λ)b + (4 − \frac{8}{15}λ)c. (2)

Eliminating b and c from these equations, and rearranging, we find that

3λ^{2} − 128λ + 640 = 0,

which has two solutions

λ = \frac{1}{3}(64 ± \sqrt{2176} ) = 5.7841 … or 36.8825 . . ..

We recall that the eigenvalues of the Sturm–Liouville equation are given by the values of the Lagrange multiplier λ. Therefore the lowest eigenvalue of this problem is approximately 5.7841 (and certainly no larger). We could find the corresponding values of b and c (and hence a) by substituting this value of λ into either equation (1) or (2) (both give the same result); note that we find only the ratio a : b : c because the normalisation of y_{trial} is not important.

In fact, the true value of the lowest eigenvalue is 5.7832 …, so the Rayleigh–Ritz method has produced an extremely good estimate.

Higher eigenvalues

Once we have found a good approximation y_{0 trial} to y_{0}, we can proceed to find approximations to the higher eigenvalues λ_{1}, λ_{2}, . . .. Just as λ_{0} is the absolute minimum of Λ[y] over all possible functions y, so λ_{1} is the absolute minimum of Λ[y] over functions which are constrained to be orthogonal to y_{0}. (Recall that y_{1} is orthogonal to y_{0} in the sense that ∫_{a}^{b} wy_{0}y_{1} dx = 0.) Hence, to estimate λ_{1} we proceed as before but choose our new trial function y_{1 trial} in such a way that it is orthogonal to our previous best approximation y_{0 trial}. This process can be continued to higher and higher eigenvalues but of course becomes less and less accurate.

Poisson’s Equation

2.1 Physical Origins

Poisson’s equation,

∇^{2}Φ = σ(x),  (2.1)

arises in many varied physical situations. Here σ(x) is the “source term”, and is often zero, either everywhere or everywhere bar some specific region (maybe only specific points). In this case, Laplace’s equation,

∇^{2}Φ = 0,  (2.2)

results.

The Diffusion Equation

Consider some quantity Φ(x) which diffuses. (This might be say the concentration of some (dilute) chemical solute, as a function of position x, or the temperature T in some heat conducting medium, which behaves in an entirely analogous way.) There is a corresponding flux, F, of Φ – that is, the amount crossing an (imaginary) unit area per unit time. Experimentally, it is known that, in the case of a solute, the flux is given by F = −k∇Φ where k is the diffusivity; in the case of temperature, the flux of heat is given by F = −k∇T where k is the coefficient of heat conductivity. (Note that the minus sign occurs because the flux is directed towards regions of lower concentration.)

The governing equation for this diffusion process is

\frac{∂Φ}{∂t} = k∇^{2}Φ  (2.3)

where k is referred to, generically, as the diffusion constant. If we are interested in the steady-state distribution of solute or of temperature, then ∂Φ/∂t = 0 and Laplace’s equation, ∇^{2}Φ = 0, follows.

When there are sources S(x) of solute (for example, where solute is piped in or where the solute is generated by a chemical reaction), or of heat (e.g., an exothermic reaction), the steady-state diffusion is governed by Poisson’s equation in the form

∇^{2}Φ = −\frac{S(x)}{k}.  (2.4)

The diffusion equation for a solute can be derived as follows. Let Φ(x) be the concentration of solute at the point x, and F(x) = −k∇Φ be the corresponding flux. (We assume here that there is no advection of Φ by the underlying medium.)

Let V be a fixed volume of space enclosed by an (imaginary) surface S. In a small time δt, the quantity of solute leaving V is given by

∫∫_{S} Fδt. n dS.  (2.5)

Hence

[∫∫∫_{V} Φ dV ]_{t}^{t+δt} = −∫∫_{S} F . n dS δt.  (2.6)

Dividing by δt and taking the limit as δt → 0,

\frac{d}{dt} ∫∫∫_{V} Φ dV = −∫∫_{S} F . n dS = ∫∫_{S} k∇Φ . n dS,  (2.7)

and hence by the Divergence Theorem,

∫∫∫_{V} \frac{∂Φ}{∂t} dV = ∫∫∫_{V} ∇. (k∇Φ) dV.  (2.8)

As this is true for any fixed volume V, we must have

\frac{∂Φ}{∂t} = ∇. (k∇Φ)  (2.9)

everywhere. Assuming that k is constant, we obtain the diffusion equation

\frac{∂Φ}{∂t} = k∇^{2}Φ.  (2.10)

If there are also sources (or sinks) of solute, then an additional source term results:

\frac{∂Φ}{∂t} = k∇^{2}Φ + S(x)  (2.11)

where S(x) is the quantity of solute (per unit volume and time) being added to the solution at the location x. Poisson’s equation for steady-state diffusion with sources, as given above, follows immediately.

The heat diffusion equation is derived similarly. Let T (x) be the temperature field in some substance (not necessarily a solid), and H(x) the corresponding heat field. We have the relation H = ρcT where ρ is the density of the material and c its specific heat. The corresponding heat flux is −k∇T. A similar argument to the above applies again, resulting in

\frac{∂H}{∂t} = k∇^{2}T + S(x)  (2.12)

where S represents possible sources of heat. Hence

\frac{∂T}{∂t} = κ∇^{2}T + (ρc)^{−1}S(x)  (2.13)

where κ = k/ρc is the coefficient of thermal diffusivity. The equation for steady-state heat diffusion with sources is as before.

Electrostatics

The laws of electrostatics are

∇. E = ρ/ϵ_{0}  ∇ × E = 0
∇. B = 0  ∇ × B = μ_{0}J  (2.13)

where ρ and J are the electric charge and current fields respectively. Since ∇ × E = 0, there is an electric potential Φ such that E = −∇Φ; hence ∇. E = ρ/ϵ_{0} gives Poisson’s equation

∇^{2}Φ = −ρ/ϵ_{0}.  (2.14)

In a region where there are no charges or currents, ρ and J vanish. Hence we obtain Laplace’s equation

∇^{2}Φ = 0.  (2.15)

Also ∇ × B = 0 so there exists a magnetostatic potential ψ such that B = −μ_{0}∇ψ; and ∇^{2}ψ = 0.

Gravitation

Consider a mass distribution with density ρ(x). There is a corresponding gravitational field F(x) which we may express in terms of

NST IB Mathematical Methods II 29

a gravitational potential Φ(x). Consider an arbitrary fixed volume V with surface S containing a total mass M_{V} = ∫∫∫_{V} ρ(x) dV. Gauss showed that the flux of the gravitational field through S is equal to −4πGM_{V}. Hence

∫∫_{S} F . n dS = −4πGM_{V}
⟹  − ∫∫_{S} ∇Φ . n dS = −4πG ∫∫∫_{V} ρ(x) dV
⟹  ∫∫∫_{V} ∇. (∇Φ) dV = 4πG ∫∫∫_{V} ρ(x) dV.

This is true for all volumes V, so we must have

∇^{2}Φ = ∇. (∇Φ) = 4πGρ.  (2.16)

Other applications

These include the motion of an inviscid fluid; Schrödinger’s equation in Quantum Mechanics; and the motion of biological organisms in a solution.

2.2 Separation of Variables for Laplace’s Equation

Plane Polar Coordinates

We shall solve Laplace’s equation ∇^{2}Φ = 0 in plane polar coordinates (r, θ) where the equation becomes

\frac{1}{r} \frac{∂}{∂r} (r\frac{∂Φ}{∂r}) + \frac{1}{r^{2}} \frac{∂^{2}Φ}{∂θ^{2}} = 0.  (2.17)

Consider solutions of the form Φ(r, θ) = R(r)Θ(θ) where each function R, Θ is a function of one variable only. Then

\frac{1}{r} \frac{∂}{∂r} (r\frac{∂Φ}{∂r}) = \frac{Θ(θ)}{r} \frac{d}{dr} (r\frac{dR}{dr})  (2.18)

and

\frac{1}{r^{2}} \frac{∂^{2}Φ}{∂θ^{2}} = \frac{R(r)}{r^{2}} \frac{d^{2}Θ}{dθ^{2}}.  (2.19)

Hence after rearrangement,

\frac{r}{R} \frac{d}{dr} (r\frac{dR}{dr}) = −\frac{Θ′′}{Θ} .  (2.20)

The LHS is a function of r only, and the RHS of θ only; hence both must be constant, λ say. Then

Θ′′ = −λΘ

⟹ Θ = \begin{cases}A + Bθ & λ = 0 \\ A cos \sqrt{λ} θ + B sin \sqrt{λ} θ & λ ≠ 0 \end{cases}

To obtain a sensible physical solution, replacing θ by θ + 2π should give the same value of ∇Φ (see later). This is true only if Θ′(θ + 2π) = Θ′(θ) ∀ θ; i.e., either λ = 0 or

cos 2π\sqrt{λ} = 1 and sin 2π\sqrt{λ} = 0  (2.21)

which implies 2π\sqrt{λ} = 2nπ for some integer n. (Note that the possibility that λ < 0 is ruled out at this stage.) Hence

Θ = \begin{cases}A + Bθ & n = 0 \\ A cos nθ + B sin nθ & n ≠ 0 \end{cases}  (2.22)

Returning to (2.18),

\frac{r}{R} \frac{d}{dr} (r\frac{dR}{dr}) = λ = n^{2}

⟹ r^{2}R′′ + rR′ − n^{2}R = 0.

It is easily shown (either by direct verification or by making the substitution u = ln r) that the solutions of this equation are

R = \begin{cases}C + D ln r & n = 0 \\ Cr^{n} + Dr^{−n} & n ≠ 0 \end{cases}  (2.23)

Hence, we obtain possible solutions to (2.17) as

Φ = RΘ = \begin{cases}(C + D ln r)(A + Bθ) & n = 0 \\ (Cr^{n} + Dr^{−n})(A cos nθ + B sin nθ) & n ≠ 0 \end{cases}  (2.24)

We note that the combination θ ln r does not satisfy the requirement above for 2πperiodicity of ∇Φ, and so we exclude it. Equation (2.17) is linear and so we may form a superposition of the above solutions; in fact the general solution is an arbitrary linear combination of all the possible solutions obtained above, that is

Φ = A_{0} + B_{0}θ + C_{0} ln r + ∑_{n=1}^{∞}(A_{n}r^{n} + C_{n}r^{−n}) cos nθ + ∑_{n=1}^{∞}(B_{n}r^{n} + D_{n}r^{−n}) sin nθ

(2.25)

where we have relabelled all the arbitrary constants, e.g., AC has become A_{n} and BD has become D_{n}. We can make this expression more compact by defining A_{−n} = C_{n} and B_{−n} = D_{n} for n > 0; then

Φ = A_{0} + B_{0}θ + C_{0} ln r + ∑_{{n=−∞ \atop n≠0}}^{∞} r^{n}(A_{n} cos nθ + B_{n} sin nθ).  (2.26)

Although this is more compact, the first expression is often easier to use.

Notes:

(i) Why did we require that ∇Φ, rather than Φ itself, be periodic? In many cases (e.g. temperature, diffusion), Φ must clearly be periodic and so we shall further need B_{0} = 0. But in other cases (e.g. electrostatics, gravitation), Φ is not itself a physical quantity, only a potential; it is ∇Φ which has a physical significance (e.g., the force). For example, consider the magnetostatic potential around a wire carrying a current I; here ψ = −(I/2π)θ, which is multi-valued, but B = −μ_{0}∇ψ (the quantity of physical interest) is of magnitude μ_{0}I/2πr and is single valued.

(ii) A common mistake made during separation of variables is to retain too many arbitrary constants; e.g. to write

∑ C_{n}r^{n}(A_{n} cos nθ + B_{n} sin nθ).  (2.27)

For each n, this looks like 3 arbitrary constants (A_{n}, B_{n}, C_{n}); but of course there are really only two arbitrary quantities (C_{n}A_{n} and C_{n}B_{n}, which we have relabelled as A_{n} and B_{n} above).

(iii) The above derivation also applies to 3D cylindrical polar coordinates in the case when Φ is independent of z.

Worked Example 5: Steady-State Temperature Distribution in a Cylinder

An infinitely long cylinder of radius a is heated on its boundary as shown. The steady-state temperature T (r, θ) (note no dependence on z) satisfies

∇^{2}T = 0 in r < a

subject to

T (a, θ) = \begin{cases}+T_{0} & 0 ≤ θ < π, \\ −T_{0} & π ≤ θ < 2π. \end{cases}

The general solution for plane polar coordinates applies; we choose to use it in its second form as given in the lecture notes. We require that the temperature be finite at r = 0 for a physically realistic solution: so C_{0} = 0, and also, for all negative n, A_{n} = B_{n} = 0 (since they are the coefficients of r^{n} \left\{{cos \atop sin}\right\} nθ). Finally, T must be periodic in θ (i.e., not multi-valued), so B_{0} = 0. Hence

T (r, θ) = A_{0} + ∑_{n=1}^{∞} r^{n}(A_{n} cos nθ + B_{n} sin nθ).

On r = a this gives

T (a, θ) = A_{0} + ∑_{n=1}^{∞}(A_{n}a^{n} cos nθ + B_{n}a^{n} sin nθ).

This is a standard Fourier series, so we may calculate the Fourier coefficients using the standard formulae:

A_{0} = \frac{1}{2π} ∫_{0}^{2π} T (a, θ) dθ = 0 (by anti-symmetry of T (a, θ))
A_{n}a^{n} = \frac{1}{π} ∫_{0}^{2π} T (a, θ) cos nθ dθ = 0
B_{n}a^{n} = \frac{1}{π} ∫_{0}^{2π} T (a, θ) sin nθ dθ
= \frac{1}{π} ∫_{0}^{π} T_{0} sin nθ dθ − \frac{1}{π} ∫_{π}^{2π} T_{0} sin nθ dθ

= \begin{cases}4T_{0}/nπ & n odd, \\ 0 & n even. \end{cases}

Hence the final solution for all r and θ is

T = \frac{4T_{0}}{π} ∑_{n odd} \frac{r^{n}}{na^{n}} sin nθ.

Spherical Polar Coordinates: Axisymmetric Case

In spherical polars (r, θ, ϕ), in the case when we know Φ to be axisymmetric (i.e., independent of ϕ, so that ∂Φ/∂ϕ = 0), Laplace’s equation becomes

\frac{1}{r^{2}} \frac{∂}{∂r} (r^{2} \frac{∂Φ}{∂r}) + \frac{1}{r^{2} sin θ} \frac{∂}{∂θ} (sin θ \frac{∂Φ}{∂θ}) = 0. (2.28)

Seek solutions of the form Φ(r, θ) = R(r)Θ(θ). Then

\frac{1}{R}(r^{2}R′)′ = −\frac{1}{Θ sin θ}(Θ′ sin θ)′  (2.29)

and both sides must be constant, λ say. So

(Θ′ sin θ)′ = −λΘ sin θ.  (2.30)

Let ζ = cos θ, and use the chain rule to replace d/dθ by d/dζ:

\frac{d}{dθ} = \frac{dζ}{dθ} \frac{d}{dζ} = − sin θ \frac{d}{dζ} .  (2.31)

So

− sin θ \frac{d}{dζ} (− sin^{2} θ \frac{dΘ}{dζ}) = −λΘ sin θ

⟹ \frac{d}{dζ} ((1 − ζ^{2})\frac{dΘ}{dζ}) + λΘ = 0.

This is Legendre’s equation; for well-behaved solutions at ζ = ±1 (i.e., θ = 0, π) we need λ = n(n + 1) for some non-negative integer n, in which case

Θ = CP_{n}(ζ) = CP_{n}(cos θ)  (2.32)

where C is an arbitrary constant.

Returning to (2.21),

(r^{2}R′)′ = λR
⟹ r^{2}R′′ + 2rR′ − n(n + 1)R = 0,

to which the solution is

R = Ar^{n} + Br^{−n−1}.  (2.33)

The general solution to Laplace’s equation in the axisymmetric case is therefore (absorbing the constant C into A and B)

Φ(r, θ) = ∑_{n=0}^{∞}(A_{n}r^{n} + B_{n}r^{−n−1})P_{n}(cos θ).  (2.34)

Non-axisymmetric Case

A similar analysis when Φ may depend on ϕ shows that the general solution is

Φ(r, θ, ϕ) = ∑_{n=0}^{∞} ∑_{m=−n}^{n}(A_{mn}r^{n} + B_{mn}r^{−n−1})P_{n}^{m}(cos θ)e^{imϕ}  (2.35)

where P_{n}^{m}(ζ) are the associated Legendre functions which satisfy the associated Legendre equation

\frac{d}{dζ} ((1 − ζ^{2})\frac{dΘ}{dζ}) + (n(n + 1) + \frac{m}{1 − ζ^{2}} ) Θ = 0)  (2.36)

when m and n are integers, n ≥ 0, −n ≤ m ≤ n.

Worked Example 6: Diffusion of a Solute past a Solid Sphere

Consider fluid at rest surrounding a fixed solid sphere of radius a at the origin. The fluid contains a solute which diffuses through the fluid, and we are interested in the steady state. At large distances from the sphere (where the sphere has negligible effect) we assume that there is a constant flux of solute parallel to the z-axis of magnitude F (possibly due, for example, to an externally imposed concentration gradient).

The flux is −k∇Φ where Φ is the concentration. There can be no flux across r = a, so ê_{r} . ∇Φ = 0 on r = a, or equivalently \frac{∂Φ}{∂r}(a, θ) = 0 for all θ.

Far from the sphere, we must have ∇Φ ∼ −\frac{F}{k} ê_{z}, i.e., Φ ∼ −\frac{F}{k} z; so we require that as r → ∞, Φ ∼ −\frac{F}{k} r cos θ = −\frac{F}{k} rP_{1}(cos θ).

We use the general axisymmetric solution, and must choose the arbitrary constants to ensure the correct behaviour as r → ∞. This can only occur if A_{1} = −\frac{F}{k} and A_{n} = 0 for all n ≥ 2. Thus

Φ = A_{0} − \frac{F}{k} rP_{1}(cos θ) + ∑_{n=0}^{∞} B_{n}r^{−n−1}P_{n}(cos θ).

On r = a, we must have

\frac{∂Φ}{∂r}|_{r=a}= −\frac{F}{k} P_{1}(cos θ) − ∑_{n=0}^{∞}(n + 1)B_{n}a^{−n−2}P_{n}(cos θ) = 0

for all θ. Using the orthogonality of Legendre polynomials (multiply by P_{m}(cos θ), substitute ζ = cos θ, and integrate from ζ = −1 to 1), or by inspection, we find that B_{0} = 0, B_{1} = −F a^{3}/2k and B_{n} = 0 for all n ≥ 2. So the solution is

Φ = A_{0} − \frac{F}{k} (r + \frac{a^{3}}{2r^{2}} ) cos θ,

and A_{0} remains an arbitrary constant (it measures, in some sense, the average of the concentrations far up and downstream).

Note that the boundary conditions involved only P_{1}(cos θ) and no other P_{n}; and so does the solution. This is usual: boundary conditions can often be expressed in terms of just a few P_{n}, and only those terms need be retained from the general solution. For this purpose it is useful to know the following:

1 = P_{0}(cos θ)
cos θ = P_{1}(cos θ)
cos^{2} θ = \frac{2}{3}P_{2}(cos θ) + \frac{1}{3}P_{0}(cos θ)

2.3 Uniqueness Theorem for Poisson’s Equation

Consider Poisson’s equation

∇^{2}Φ = σ(x)  (2.37)

in a volume V with surface S, subject to so-called Dirichlet boundary conditions Φ(x) = f(x) on S, where f is a given function defined on the boundary.

From a physical point of view, we have a well-defined problem; say, find the steady-state temperature distribution throughout V, with heat sources given by σ(x), subject to a specified temperature distribution on the boundary. No further conditions are required in real life to ensure that there is only one solution. Mathematically, then, can we show that the problem above has a unique solution?

Suppose that there are actually two (or more) solutions Φ_{1}(x) and Φ_{2}(x). Let Ψ = Φ_{1} − Φ_{2}. Then

∇^{2}Ψ = ∇^{2}Φ_{1} − ∇^{2}Φ_{2} = σ − σ = 0 in V  (2.38)

subject to

Ψ = f − f = 0 on S.  (2.39)

One solution of this problem for Ψ is clearly Ψ = 0; is it unique? Consider

∇. (Ψ∇Ψ) = ∇Ψ . ∇Ψ + Ψ∇. (∇Ψ)
= |∇Ψ|^{2} + Ψ∇^{2}Ψ
= |∇Ψ|^{2}.

Hence

∫∫∫_{V} |∇Ψ|^{2} dV = ∫∫∫_{V} ∇. (Ψ∇Ψ) dV

= ∫∫_{S} Ψ∇Ψ . n dS

= 0

because Ψ = 0 on S. But |∇Ψ|^{2} ≥ 0 everywhere; its integral can only be zero if |∇Ψ| is zero everywhere, i.e., ∇Ψ ≡ 0, which implies that Ψ is constant throughout V. But Ψ = 0 on S, so Ψ ≡ 0 throughout V. Thus Φ_{1} = Φ_{2}, which demonstrates that our problem has a unique solution, as expected.

A similar theorem holds when instead of Dirichlet boundary conditions we have Neumann boundary conditions: that is to say instead of Φ being specified (by the function f) on the boundary S, ∂Φ/∂n is specified on S, where we use the notation

\frac{∂Φ}{∂n} ≡ n . ∇Φ.  (2.40)

NST IB Mathematical Methods II 37

2.4 Minimum and Maximum Properties of Laplace’s Equation

Suppose that Φ satisfies ∇^{2}Φ = 0 in a volume V with surface S. Then both the minimum and maximum values of Φ occur somewhere on S (and possibly also somewhere inside V).

Why is this? Suppose Φ has a local maximum somewhere in the interior of V. At that point we must have ∂Φ/∂x = ∂Φ/∂y = ∂Φ/∂z = 0 (stationary point); and as it is a maximum, ∂^{2}Φ/∂x^{2} < 0, ∂^{2}Φ/∂y^{2} < 0, ∂^{2}Φ/∂z^{2} < 0. But this cannot happen since 0 = ∇^{2}Φ = ∂^{2}Φ/∂x^{2} + ∂^{2}Φ/∂y^{2} + ∂^{2}Φ/∂z^{2}.

The same applies to minima.

This is not a formal proof since it is actually possible for a maximum to have

\frac{∂^{2}Φ}{∂x^{2}} = \frac{∂^{2}Φ}{∂y^{2}} = \frac{∂^{2}Φ}{∂z^{2}} = 0,  (2.41)

a case we haven’t considered: compare with the possibility in 1D that a maximum could have d^{2}y/dx^{2} = 0. However, the theorem can still be shown to hold.

Example: in the worked example of the steady-state temperature distribution in a cylinder, we can deduce that |T | ≤ T_{0} in r < a.

2.5 Green’s Function

The Delta Function in 3D

In 1D, δ(x − x_{0}) is a function which is zero everywhere except at x = x_{0}, and is infinite there in such a way that

∫_{a}^{b} δ(x − x_{0}) dx = 1  (2.42)

whenever x_{0} ∈ (a, b). As a consequence, ∫_{a}^{b} f(x)δ(x − x_{0}) dx = f(x_{0}). We extend the definition to 3D via

δ(x − x_{0}) = δ(x − x_{0})δ(y − y_{0})δ(z − z_{0})  (2.43)

where x_{0} = (x_{0}, y_{0}, z_{0}). Then

∫∫∫_{V} f(x)δ(x − x_{0}) dV = f(x_{0})  (2.44)

whenever x_{0} ∈ V (and the integral is 0 otherwise).

Green’s Function

Suppose that we wish to solve Poisson’s equation in a volume V with surface S on which Dirichlet boundary conditions are imposed. The Green’s function G(x; x_{0}) associated with this problem is a function of two variables: x, the position vector, and x_{0}, a fixed location. It is defined as the solution to

∇^{2}G(x; x_{0}) = δ(x − x_{0}) in V,
G = 0  on S.

(2.45)

(Physically, we can think of G as the “potential” from a point source at x_{0} with the boundary held at zero potential.)

It is possible to prove that G is symmetric, i.e., G(x; x_{0}) = G(x_{0}; x). This can be useful as a check that G has been correctly calculated. Physically, this corresponds to the remarkable fact that the potential at x due to a source at x_{0} is the same as the potential at x_{0} due to a source at x, regardless of the shape of S.

When V is all space (i.e., the limit of a sphere whose radius tends to ∞), Green’s function is known as the fundamental solution.

For a problem with Neumann boundary conditions, G(x; x_{0}) is defined to satisfy ∂G/∂n = 1/A on S, where A = ∫∫_{S} dS is the surface area of S, rather than G = 0 there. In many cases S is infinitely large, in which case the boundary condition reduces to ∂G/∂n = 0.

The Fundamental Solution in 3D

Consider first x_{0} = 0. Then ∇^{2}G = δ(x) and G → 0 as |x| → ∞. The problem is spherically symmetric about the origin, so we expect G to be a function of r alone. Try G = g(r). By the definition of δ(x), if V_{R} is the sphere of radius R with surface S_{R},

1 = ∫∫∫_{V_{R}} δ(x) dV  = ∫∫∫_{V_{R}} ∇. (∇G) dV

= ∫∫_{S_{R}} ∇G. n dS  = ∫∫_{S_{R}} g′(r) dS

(n is just the unit radial vector)

= g′(R) ∫∫_{S_{R}} dS  = 4πR^{2}g′(R)

⟹ g′(R) = \frac{1}{4πR^{2}}  for all R
⟹ g′(r) = \frac{1}{4πr^{2}}
⟹ g(r) = −\frac{1}{4πr} + A,  (2.36)

where A is a constant. As r → ∞, G → 0, so A = 0. Hence the solution is −1/4π|x|.

Shifting the origin to a non-zero x_{0}, we see that in general the fundamental solution in 3D is

G(x; x_{0}) = −\frac{1}{4π|x − x_{0}|}.  (2.46)

Example: an electron located at x_{0} is an electrostatic point source, so the charge distribution in space is ρ(x) = −e δ(x − x_{0}). Hence the electrostatic potential obeys

∇^{2}Φ = (e/ϵ_{0}) δ(x − x_{0})  (2.47)

using a result from §2.1. The solution Φ is therefore just a factor e/ϵ_{0} times the fundamental solution, i.e., −e/4πϵ_{0}|x − x_{0}|. This is the standard formula for the potential due to an electron.

The Fundamental Solution in 2D

Again, we solve ∇^{2}G = δ(x), where the delta-function is now in 2D. We will see that a solution with G → 0 as |x| → ∞ is impossible; instead we will find a solution such that |∇G| → 0.

As before, G = g(r) (where r is now the plane polar radius). Applying the Divergence Theorem in 2D to a circle of radius R,

1 = ∫∫_{r≤R} δ(x) dV = ∫∫_{r≤R} ∇. (∇G) dV

= ∮_{r=R} ∇G. n dl = ∮_{r=R} g′(r) dl
= 2πRg′(R)
⟹ g′(r) = \frac{1}{2πr}
⟹ g(r) = \frac{1}{2π} ln r + constant.

(Note that g′(r) → 0 as r → ∞, but g(r) → ∞, whatever the constant.)

Shifting the origin, we see that the fundamental solution in 2D is

G(x; x_{0}) = \frac{1}{2π} ln |x − x_{0}| + constant.

(2.48)

Example: consider an infinitely long charged wire in three dimensions lying along the z-axis, with a charge density of μ per unit length. What is the electric potential Φ around the wire?

We assume the wire to be mathematically perfect, i.e., of infinitesimal width. Then the electric charge distribution, in 3D, is ρ = μδ(x)δ(y). (Check that this gives the correct result for the amount of charge in a unit length of the wire.) But it is clear that this problem is fundamentally two-dimensional, with ρ = μδ(x) where x = (x, y); and the potential satisfies ∇^{2}Φ = −μδ(x)/ϵ_{0}. Hence the potential is (up to an arbitrary additional constant) just given by an appropriate multiple of the two-dimensional fundamental solution, namely

Φ = −\frac{μ}{2πϵ_{0}} ln |x| = −\frac{μ}{2πϵ_{0}} ln \sqrt{x^{2} + y^{2}} = −\frac{μ}{2πϵ_{0}} ln r  (2.49)

where r is the perpendicular distance to the wire (i.e., the “r” of cylindrical polar coordinates rather than of spherical polars).

2.6 The Method of Images

We can use the fundamental solution to find Green’s function in some simple geometries, using the “method of images”. We shall find a function which satisfies the equation and the boundary conditions; by uniqueness, this must be the Green’s function.

Example: A 3D half-space x > 0

Suppose that the domain D is the half-space of ℝ^{3} with x > 0. The Green’s function obeys

∇^{2}G = δ(x − x_{0}) ∀ x ∈ D,
G = 0  on x = 0,
G → 0  as |x| → ∞, x ∈ D.

Consider the solution in all space for the point source at x = x_{0} together with another (imaginary) source of strength −1 at the “image point” x = x_{1} as shown:

Φ = −\frac{1}{4π|x − x_{0}|} − \frac{−1}{4π|x − x_{1}|}  (2.50)

and

∇^{2}Φ = δ(x − x_{0}) − δ(x − x_{1})  (2.51)

by superposition of the two fundamental solutions. This certainly satisfies the requirement ∇^{2}Φ = δ(x − x_{0}) for all x ∈ D, because δ(x − x_{1}) ≡ 0 ∀ x ∈ D. It also satisfies Φ → 0 as |x| → ∞; and on x = 0, |x−x_{0}| = |x−x_{1}| so that Φ = 0. Hence by uniqueness,

G(x; x_{0}) = Φ = −\frac{1}{4π} (\frac{1}{|x − x_{0}|} − \frac{1}{|x − x_{1}|}) .  (2.52)

Example: A 2D quarter-plane x > 0, y > 0

In this case, we need to find G such that

∇^{2}G = δ(x − x_{0}) ∀ x ∈ D  (2.53)

with G = 0 on both x = 0 and y = 0. We find that we need 3 image sources as shown: x_{1} and x_{2} with strength −1, and x_{3} with strength +1. Then

G = \frac{1}{2π} ln |x − x_{0}| − \frac{1}{2π} ln |x − x_{1}| − \frac{1}{2π} ln |x − x_{2}| + \frac{1}{2π} ln |x − x_{3}| + constant
= \frac{1}{2π} ln \frac{|x − x_{0}| |x − x_{3}|}{|x − x_{1}| |x − x_{2}|} + constant.

Clearly ∇^{2}G = δ(x − x_{0}) in D (all the other delta-functions are zero there); on x = 0, |x−x_{0}| = |x−x_{1}| and |x−x_{2}| = |x−x_{3}|, so choosing the constant to be zero ensures that G = 0; similarly on y = 0. By uniqueness, then, this is the required Green’s function.

Further extensions to this idea are possible; e.g., planes inclined at 60^{∘} to each other, or a pair of parallel planes.

Example: Heat flow from a source in a 3D half-space with a wall at constant temperature

Suppose that the ambient temperature is T_{0} and that a wall at x = 0 is held at that temperature, with a heat source of strength Q at x_{0}. Then

T = T_{0} − \frac{Q}{k}G(x; x_{0}),  (2.54)

where G is the Green’s function for the 3D half-space x > 0. (Why? Because we need to solve ∇^{2}T = −\frac{Q}{k}δ(x − x_{0}) here.)

What is the total heat flux across the wall S? It is

∫∫_{S} (−k∇T) . n dS = k ∫_{−∞}^{∞} ∫_{−∞}^{∞} \frac{∂T}{∂x} dy dz = −Q ∫_{−∞}^{∞} ∫_{−∞}^{∞} \frac{∂}{∂x}G(x; x_{0})|_{x=0} dy dz (2.55)

which we can evaluate with some effort (see the worked example in the next section for an example of this sort of evaluation).

Alternatively, we can use the Divergence Theorem on the surface consisting of the wall plus the hemisphere at ∞. Since ∇T tends to zero on the hemisphere,

∫∫_{S} (−k∇T) . n dS = − ∫∫∫_{V} ∇. (k∇T) dV

= −k ∫∫∫_{V} ∇^{2}T dV

= −k ∫∫∫_{V} (−\frac{Q}{k}δ(x − x_{0})) dV

= Q,

so the total heat radiated across the wall is Q.

Example: A point charge near an earthed boundary plate

Here

Φ = −\frac{e}{ϵ_{0}} G(x; x_{0})  (2.56)

where G is the Green’s function for the 3D half-space x > 0.

Now the surface charge density induced on the plate is μ = ϵ_{0}E_{x} (standard result from electrostatics, where E_{x} is the x-component of E). The normal force (per unit area) on the plate, towards the charge, is

\frac{1}{2}μE_{x} = \frac{1}{2}ϵ_{0}E_{x}^{2} = \frac{1}{2}ϵ_{0} (−\frac{∂Φ}{∂x})^{2} = \frac{e^{2}}{2ϵ_{0}} (\frac{∂G}{∂x})^{2} ,  (2.57)

and we calculate ∂G/∂x as in the worked example in the next section. We can integrate this over the whole plate (with considerable effort) to obtain the total force:

\frac{e^{2}}{2ϵ_{0}} ∫_{−∞}^{∞} ∫_{−∞}^{∞} \frac{x_{0}^{2}}{4π^{2}(x_{0}^{2} + (y − y_{0})^{2} + (z − z_{0})^{2})^{3}} dy dz = ⋯ = \frac{e^{2}}{16πϵ_{0}x_{0}^{2}} .  (2.58)

The force on the charge from the plate is equal and opposite, i.e., e^{2}/4πϵ_{0}(2x_{0})^{2} towards the wall. Note that we could also have found this directly by considering the force on the charge due to the image charge, ignoring the plate!

Worked Example 7: Heat Source near an Insulated Wall

Hold a heat source of strength Q at x_{0} = (x_{0}, y_{0}, z_{0}) near an insulated plane wall, i.e., one through which no heat can pass, at x = 0. We must then have no component of heat flux through the wall; i.e., n.(−k∇T) = 0 on the wall. Therefore we must solve

∇^{2}T = −\frac{Q}{k}δ(x − x_{0}) in x > 0

subject to

\frac{∂T}{∂n} = 0 on x = 0.

This is a problem with Neumann (rather than Dirichlet) boundary conditions.

We use the method of images. Introduce an image source of strength +Q at x_{1} = (−x_{0}, y_{0}, z_{0}). (Note that for Dirichlet boundary conditions we would have used −Q for the strength of the image.) Because ∇T is radial from each source, the total ∇T (from the two sources combined) must have zero component perpendicular to the wall. Hence we have ∂T/∂n = 0 as required. Therefore (by uniqueness) the solution is

T = \frac{Q}{4πk} \{\frac{1}{|x − x_{0}|} + \frac{1}{|x − x_{1}|}\}.

NST IB Mathematical Methods II 45

Example: Images in a sphere What is the Green’s function for the domain r < a in 3D? We need

∇^{2}G = δ(x − x_{0}) in r < a,
G = 0  on r = a.

The image point turns out to be at the inverse point

x_{1} = \frac{a^{2}}{|x_{0}|^{2}} x_{0}  (2.59)

(so that a/|x_{1}| = |x_{0}|/a) with strength −a/|x_{0}|, so Green’s function is

G(x; x_{0}) = \frac{1}{4π} (−\frac{1}{|x − x_{0}|} + \frac{a/|x_{0}|}{|x − x_{1}|}) .  (2.60)

(Check this by first showing that |x − x_{1}|^{2} = (x − x_{1}) . (x − x_{1}) = (a^{2}/|x_{0}|^{2})|x − x_{0}|^{2} when |x| = a.)

Note that the same result holds if we consider the domain r > a instead.

Example: Images in a circle

This is the 2D equivalent of the above. The image point is at x_{1} = (a^{2}/|x_{0}|^{2})x_{0} as before, but now the strength of the image can be chosen to be −1, so the Green’s function is

G(x; x_{0}) = \frac{1}{2π} ln |x − x_{0}| − \frac{1}{2π} ln |x − x_{1}| + constant

= \frac{1}{2π} ln \frac{|x − x_{0}|}{|x − x_{1}|} + constant. (2.52)

Choosing the constant correctly, we can ensure that G = 0 on the circle r = a.

NST IB Mathematical Methods II 46

2.7 The Integral Solution of Poisson’s Equation

The most important application of Green’s function is that it can be used to find the solution of Poisson’s equation with an arbitrary source distribution.

Green’s Identity

For any smooth functions Φ and Ψ, Green’s Identity is ∫∫∫_{V} (Φ∇^{2}Ψ − Ψ∇^{2}Φ) dV = ∫∫_{S} (Φ∇Ψ − Ψ∇Φ) . n dS

(2.61)

where V is a volume with surface S. Prove this by applying the Divergence Theorem to the vector field F = Φ∇Ψ − Ψ∇Φ, and using ∇. (Φ∇Ψ) = ∇Φ . ∇Ψ + Φ∇^{2}Ψ.

The RHS is also written

∫∫_{S} (Φ\frac{∂Ψ}{∂n} − Ψ\frac{∂Φ}{∂n}) dS.  (2.62)

NST IB Mathematical Methods II 47

The Integral Solution

Consider the general problem of Poisson’s equation with Dirichlet boundary conditions:

∇^{2}Φ = σ in V,
Φ = f  on S.

Apply Green’s Identity, taking Ψ to be the Green’s function G(x; x_{0}) for the problem:

∫∫∫_{V} (Φ∇^{2}G − G∇^{2}Φ) dV = ∫∫_{S} (Φ∇G − G∇Φ) . n dS

δ(x − x_{0}) σ  f  0
⟹  ∫∫∫_{V} Φδ(x − x_{0}) dV = ∫∫∫_{V} Gσ dV + ∫∫_{S} f \frac{∂G}{∂n} dS
⟹  Φ(x_{0}) = ∫∫∫_{V} σ(x)G(x; x_{0}) dV + ∫∫_{S} f(x)\frac{∂G}{∂n} dS.  (2.55)

This is the Integral Solution of Poisson’s equation.

Notes:

(i) We can also use the integral solution to solve Laplace’s equation with Dirichlet boundary conditions, by taking σ(x) = 0.

(ii) A similar result (but with technical differences) can be derived for Neumann boundary conditions, in which case G is defined differently (see §2.5).

(iii) We sometimes wish to take V to be “all space”, by taking the limit of a sphere whose radius tends to ∞. In this case we simply use the fundamental solution for G; but (strictly speaking) we need to ensure that the surface integral tends to zero (by requiring, for example, that on the surface of the sphere, Φ(x) → 0 sufficiently quickly as the radius increases). Then

Φ(x_{0}) = ∫∫∫_{ℝ^{3}} σ(x)G(x; x_{0}) dV.  (2.63)

This latter result is easy to understand in many physical situations. For instance, consider an arbitrary electrostatic charge distribution ρ(x). Then

∇^{2}Φ = −ρ/ϵ_{0} in ℝ^{3},
Φ → 0  as |x| → ∞.

(We assume here that the charge distribution decays rapidly far from the origin.) Using the integral solution of Poisson’s equation, with V = ℝ^{3}, and setting G to be the fundamental solution in 3D,

Φ(x_{0}) = ∫∫∫_{ℝ^{3}} \frac{ρ(x)}{4πϵ_{0}|x − x_{0}|} dV.  (2.64)

We can interpret this physically as the superposition of many infinitesimal charge elements ρ(x) dV. Each of these is effectively a point charge, and the potential at x_{0} from such a point charge (using the standard formula for the electrostatic potential due to a point charge) is just ρ(x) dV/4πϵ_{0}|x − x_{0}|. Summing over all such infinitesimal elements gives the above result.

Worked Example 8: Electrostatics - Using the Integral Solution of Poisson’s Equation

Consider a wire of length 2L carrying a charge density μ per unit length, lying along the z-axis from z = −L to +L. What is the electric potential Φ?

The charge distribution is ρ(x) = μδ(x)δ(y) for −L ≤ z ≤ L (and zero for |z| > L). We shall use the integral solution of Poisson’s equation in the whole of space to obtain the potential at a point (x_{0}, y_{0}, z_{0}). We need Green’s function, which is simply the fundamental solution here.

Φ(x_{0})  = ∫∫∫_{ℝ^{3}} \frac{ρ(x)}{4πϵ_{0}|x−x_{0}|} dV
= ∫_{−L}^{L} \frac{μ}{4πϵ_{0}|(0,0,z)−x_{0}|} dz
= \frac{μ}{4πϵ_{0}} ∫_{−L}^{L} _{\frac{dz}{\sqrt{x_{0}^{2}+y_{0}^{2}+(z−z_{0})^{2}}}}
= \frac{μ}{4πϵ_{0}} [sinh^{−1} _{\frac{z−z_{0}}{\sqrt{x_{0}^{2}+y_{0}^{2}}}} ]_{−L}^{L}
= \frac{μ}{4πϵ_{0}} \{sinh^{−1} _{\frac{L−z_{0}}{\sqrt{x_{0}^{2}+y_{0}^{2}}}} + sinh^{−1} _{\frac{L+z_{0}}{\sqrt{x_{0}^{2}+y_{0}^{2}}}} \}.

This is true for arbitrary locations x_{0}, so replacing x_{0} by x we obtain

Φ(x, y, z) = \frac{μ}{4πϵ_{0}} \{sinh^{−1} \frac{L − z}{\sqrt{x^{2} + y^{2}}} + sinh^{−1} \frac{L + z}{\sqrt{x^{2} + y^{2}}} \}.

In particular, the potential at a point in the (x, y)-plane is given by

Φ(x, y, 0) = \frac{μ}{2πϵ_{0}} sinh^{−1}(L/\sqrt{x^{2} + y^{2}}).

Note, for completeness, that for very large L, i.e., in the limit as L → ∞, it is possible to check (using sinh^{−1} x ∼ ln x as x → ∞) that

Φ → −\frac{μ}{2πϵ_{0}} ln \sqrt{x^{2} + y^{2}} + constant,

which verifies an earlier result we obtained for the two-dimensional field around an infinitely long wire.

Worked Example 9: Solution of Laplace’s Equation in a 3D Half-Space

We wish to solve ∇^{2}Φ = 0 in the half-space x > 0 of ℝ^{3}, with Φ = f(y, z) on the boundary x = 0.

We use the integral solution of Poisson’s equation (with σ ≡ 0) in the half-space, with S being the plane x = 0 (strictly speaking, together with the hemisphere at ∞):

Φ(x_{0}) = ∫∫∫_{V} σ(x)G(x; x_{0}) dV + ∫∫_{S} f(x)\frac{∂G}{∂n} dS
= − ∫_{−∞}^{∞} ∫_{−∞}^{∞} f(y, z)\frac{∂}{∂x}G(x; x_{0}) dy dz

(because \frac{∂}{∂n} = −\frac{∂}{∂x} on S). To calculate this we need to evaluate

\frac{∂G}{∂x}|_{x=0} = \frac{∂}{∂x}\{−\frac{1}{4π|x − x_{0}|} + \frac{1}{4π|x − x_{1}|}\}|_{x=0}
= \frac{1}{4π} \frac{∂}{∂x}\{−\frac{1}{\sqrt{(x − x_{0})^{2} + (y − y_{0})^{2} + (z − z_{0})^{2}}}

+ \frac{1}{\sqrt{(x + x_{0})^{2} + (y − y_{0})^{2} + (z − z_{0})^{2}}} \}|_{x=0}

= \frac{1}{4π} \{\frac{x − x_{0}}{\{(x − x_{0})^{2} + (y − y_{0})^{2} + (z − z_{0})^{2}\}^{3/2}}

− \frac{x + x_{0}}{\{(x + x_{0})^{2} + (y − y_{0})^{2} + (z − z_{0})^{2}\}^{3/2}} \}|_{x=0}

= −\frac{x_{0}}{2π\{x_{0}^{2} + (y − y_{0})^{2} + (z − z_{0})^{2}\}^{3/2}} .

Therefore

Φ(x_{0}) = \frac{x_{0}}{2π} ∫_{−∞}^{∞} ∫_{−∞}^{∞} \frac{f(y, z)}{\{x_{0}^{2} + (y − y_{0})^{2} + (z − z_{0})^{2}\}^{3/2}} dy dz

or alternatively (swapping x and x_{0}),

Φ(x, y, z) = \frac{x}{2π} ∫_{−∞}^{∞} ∫_{−∞}^{∞} \frac{f(y_{0}, z_{0})}{\{x^{2} + (y − y_{0})^{2} + (z − z_{0})^{2}\}^{3/2}} dy_{0} dz_{0}.

This is the solution for:

(i) Steady-state temperature distribution with a wall heated to a specified temperature distribution; (ii) Steady-state concentration of solute with a wall kept at given concentration; (iii) Electrostatic potential with a conducting wall held at given potential.

Cartesian Tensors

3.1 [Revision] Suffix Notation and the Summation Convention

We will for the moment consider vectors in 3D, though the notation we shall introduce applies (mostly) just as well to n dimensions. For a general vector

x = (x_{1}, x_{2}, x_{3})  (3.1)

we shall henceforth refer to x_{i}, the i^{th} component of x. The index i may take any of the values 1, 2 or 3, and we refer to “the vector x_{i}” to mean “the vector whose components are (x_{1}, x_{2}, x_{3})”. Similarly we write [x]_{i} = x_{i}, and similarly [x + y]_{i} = x_{i} + y_{i}. Note that the expression y_{i} = x_{i} implies that y = x; the statement in suffix notation is implicitly true for all three possible values of i (separately).

We will use the summation convention whereby if a particular suffix (e.g., i) appears twice in a single term of an expression then it is implicitly summed. For example, in summation notation we simply write

x . y = x_{i}y_{i}.

Rules of summation convention

Summation convention does not allow any one suffix to appear more than twice within a single term; so x_{i}y_{i}z_{i} is meaningless. Care should be taken to avoid this. (For example, the vector relation y = (a . b)x must be written y_{i} = a_{j}b_{j}x_{i}, rather than y_{i} = a_{i}b_{i}x_{i}.)

In any given term, then, there are two possible types of suffix: one that appears precisely once, e.g., i in a_{j}b_{j}x_{i}, which is known as a free suffix; and one that appears precisely twice, e.g., j in a_{j}b_{j}x_{i}, which is known as a dummy suffix. It is an important precept of summation convention that the free suffixes must match precisely in every term (though dummy suffixes can be anything you like so long as they do not clash with the free suffixes). So in the equation

a_{j}b_{j}z_{k} = x_{k} + a_{i}a_{i}y_{k}b_{j}b_{j}  (3.2)

every term has a free suffix k, and all other suffixes are dummy ones. In vector notation, this equation reads

(a . b)z = x + |a|^{2}|b|^{2}y.

There need not be any free suffixes at all, as in the equation a_{i}z_{i} = (x_{i} + y_{i})a_{i} (which reads a . z = (x + y) . a in vector notation). Replacing two free suffixes (e.g. i, j in c_{ij}) by a single dummy suffix (c_{ii}) is known as contraction.

Note that (1) the order of variables written in suffix notation is unimportant; the final term of equation (3.2) could equally well have been written b_{j}y_{k}a_{i}b_{j}a_{i}; and (2) the role of the dummy suffix is analogous to that of the dummy variable in an integration.

More examples:

(i) y = Ax is written y_{i} = [Ax]_{i} = a_{ij}x_{j}.

(ii) C = AB (where A and B are 3 × 3 matrices) is written c_{ij} = [AB]_{ij} = a_{ik}b_{kj}.

(iii) A matrix C has trace Tr C = c_{ii}, so the trace of AB becomes Tr(AB) = a_{ik}b_{ki}.

[Not all expressions written in suffix notation can be recast in vector or matrix notation. For example, a_{ijk} = x_{i}y_{j}z_{k} is a valid equation in suffix notation (each term has three free suffixes, i, j and k), but there is no vector equivalent.]

The Kronecker delta and the alternating tensor:

The Kronecker delta is defined by

δ_{ij} = \begin{cases}1 & i = j \\ 0 & i ≠ j \end{cases}

and the alternating tensor (also called Levi-Civita or permutation tensor) is defined by

ϵ_{ijk} = \begin{cases}1 & if (i, j, k) is a cyclic permutation of (1, 2, 3) \\ −1 & if (i, j, k) is an anti-cyclic permutation of (1, 2, 3) \\ 0 & if any of i, j, k are equal \end{cases}

If I is the n × n identity matrix then [I]_{ij} = δ_{ij}, and x_{i} = δ_{ij}x_{j} for any vector x_{i} (because x = Ix; or because the Kronecker delta just “selects” entries: e.g., δ_{ik}a_{jk} is equal to a_{ji}).

Using the alternating tensor the expression z = x × y can be written z_{i} = [x × y]_{i} = ϵ_{ijk}x_{j}y_{k}. Similarly the determinant of a 3 × 3 matrix A = (a_{ij}) is given by ϵ_{ijk}a_{1i}a_{2j}a_{3k}. This can be written in several other ways; for example,

det A = ϵ_{ijk}a_{1i}a_{2j}a_{3k} = ϵ_{jik}a_{1j}a_{2i}a_{3k}
[swapping i and j] = −ϵ_{ijk}a_{2i}a_{1j}a_{3k}.

This proves that swapping two rows of a matrix changes the sign of the determinant.

The relation ϵ_{ijk}ϵ_{klm} = δ_{il}δ_{jm} − δ_{im}δ_{jl} holds, and simplifies the proof of many vector identities, such as the vector triple product a × (b × c) = (a . c)b − (a . b)c:

[a × (b × c)]_{i} = ϵ_{ijk}a_{j}[b × c]_{k}
= ϵ_{ijk}a_{j}ϵ_{klm}b_{l}c_{m}
= (δ_{il}δ_{jm} − δ_{im}δ_{jl})a_{j}b_{l}c_{m}
= a_{j}b_{i}c_{j} − a_{j}b_{j}c_{i}
= (a . c)b_{i} − (a . b)c_{i}
= [(a . c)b − (a . b)c]_{i} .

3.2 What is a Vector?

A vector is more than just 3 real numbers. It is also a physical entity: if we know its 3 components with respect to one set of Cartesian axes then we know its components with respect to any other set of Cartesian axes. (The vector stays the same even if its components do not.)

For example, suppose that {e_{1}, e_{2}, e_{3}} is a right-handed orthogonal set of unit vectors, and that a vector v has components v_{i} relative to axes along those vectors. That is to say,

v = v_{1}e_{1} + v_{2}e_{2} + v_{3}e_{3} = v_{j}e_{j}.  (3.3)

What are the components of v with respect to axes which have been rotated to align with a different set of unit vectors {e′_{1}, e′_{2}, e′_{3}}? Let

v = v′_{1}e′_{1} + v′_{2}e′_{2} + v′_{3}e′_{3} = v′_{j}e′_{j}.  (3.4)

Now e′_{i} . e′_{j} = δ_{ij}, so

v . e′_{i} = v′_{j}e′_{j} . e′_{i} = v′_{j}δ_{ij} = v′_{i}  (3.5)

but also

v . e′_{i} = v_{j}e_{j} . e′_{i} = v_{j}l_{ij}  (3.6)

where we define the matrix L = (l_{ij}) by

l_{ij} = e′_{i} . e_{j}.  (3.7)

Then

v′_{i} = l_{ij}v_{j}  (3.8)

(or, in matrix notation, v′ = Lv where v′ is the column vector with components v′_{i}). L is called the rotation matrix.

This looks like, but is not quite the same as, rotating the vector v round to a different vector v′ using a transformation matrix L. In the present case, v and v′ are the same vector, just measured with respect to different axes. The transformation matrix corresponding to the rotation {e_{1}, e_{2}, e_{3}} ↦ {e′_{1}, e′_{2}, e′_{3}} is not L (in fact it is L^{−1}).

Now consider the reverse of this argument. Exactly the same discussion would lead to

v_{i} = l̂_{ij}v′_{j}  (3.9)

where

l̂_{ij} = e_{i} . e′_{j}  (3.10)

(we swap primed and unprimed quantities throughout the argument). We note that l̂_{ij} = l_{ji} from their definitions; hence

L̂ = L^{T}  (3.11)

and so

v = L̂v′ = L^{T} v′.  (3.12)

We can deduce that

v = L^{T} Lv,  (3.13)

and furthermore, this is true for all vectors v. We conclude that

L^{T} L = I,  (3.14)

i.e.,

L^{T} = L^{−1}.  (3.15)

(Hence LL^{T} = I also.) L is therefore an orthogonal matrix. In suffix notation, the equation L^{T} L = I reads

l_{ki}l_{kj} = δ_{ij},  (3.16)

and LL^{T} = I reads

l_{ik}l_{jk} = δ_{ij};  (3.17)

both of these identities will be useful.

Another way of seeing that LL^{T} = I (or, equivalently, L^{T} L = I) is to consider the components of L directly:

L = (\begin{matrix}e′_{1} . e_{1} & e′_{1} . e_{2} & e′_{1} . e_{3} \\ e′_{2} . e_{1} & e′_{2} . e_{2} & e′_{2} . e_{3} \\ e′_{3} . e_{1} & e′_{3} . e_{2} & e′_{3} . e_{3} \end{matrix})

= ( \begin{matrix}e′_{1} & T \\ e′_{2} & T \\ e′_{3} & T \end{matrix}  )  [measured with respect to frame 1].

Alternatively, the i^{th} column consists of the components of e_{i} with respect to the second frame.

3.3 Tensors

Tensors are a generalisation of vectors. We think informally of a tensor as something which, like a vector, can be measured component-wise in any Cartesian frame; and which also has a physical significance independent of the frame, like a vector.

Physical Motivation

Recall the conductivity law, J = σE, where E is the applied electric field and J is the resulting electric current. This is suitable for simple isotropic media, where the conductivity is the same in all directions. But a matrix formulation may be more suitable in anisotropic media; for example,

J = (\begin{matrix}5 & 0 & 0 \\ 0 & 4 & 0 \\ 0 & 0 & 0 \end{matrix}) E  (3.18)

might represent a medium in which the conductivity is high in the x-direction but in which no current at all can flow in the z-direction. (For instance, a crystalline lattice structure where vertical layers are electrically insulated.)

More generally, in suffix notation we have

J_{i} = σ_{ij}E_{j}  (3.19)

where σ is the conductivity tensor.

What happens if we measure J and E with respect to a different set of axes? We would expect the matrix σ to change too: let its new components be σ′_{ij}. Then

J′_{i} = σ′_{ij}E′_{j}.  (3.20)

But J and E are vectors, so

J′_{i} = l_{ij}J_{j}  (3.21)

and

E_{i} = l_{ji}E′_{j}  (3.22)

from the results regarding the transformation of vectors in §3.2. Hence

σ′_{ij}E′_{j} = J′_{i}
= l_{ip}J_{p}
= l_{ip}σ_{pq}E_{q}
= l_{ip}σ_{pq}l_{jq}E′_{j}
⟹ (σ′_{ij} − l_{ip}l_{jq}σ_{pq})E′_{j} = 0.

This is true for all vectors E′, and hence the bracket must be identically zero; hence σ′_{ij} = l_{ip}l_{jq}σ_{pq}. This tells us how σ transforms.

Compare this argument with the corresponding argument for the case Ax = 0 where A is a matrix; if it is true for all x then A must be zero, though this is not the case if it is only true for some x’s.

σ is a second rank tensor, because it has two suffixes (σ_{ij}).

Definition: In general, a tensor of rank n is a mathematical object with n suffixes, T_{ijk...}, which obeys the transformation law

T′_{ijk...} = l_{ip}l_{jq}l_{kr} … T_{pqr...}  (3.23)

where L is the rotation matrix between frames.

Note: for second rank tensors such as σ, the transformation law

T′_{ij} = l_{ip}l_{jq}T_{pq}  (3.24)

can be rewritten in matrix notation as T′ = LT L^{T} – check this yourself!

Examples of Tensors

(i) Any vector v (e.g., velocity) is a tensor of rank 1, because v′_{i} = l_{ip}v_{p}.

(ii) Temperature T is a tensor of rank 0 – known as a scalar – because it is the same in all frames (T′ = T).

(iii) The inertia tensor. Consider a mass m which is part of a rigid body, at a location x within the body. If the body is rotating with angular velocity ω then the mass’s velocity is v = ω × x, and its angular momentum is therefore

mx × v = mx × (ω × x) = m(|x|^{2}ω − (ω . x)x).  (3.25)

Changing from a single mass m to a continuous mass distribution with density ρ(x), so that an infinitesimal mass element is ρ(x) dV, we see that the total angular momentum of a rigid body V is given by

h = ∫∫∫_{V} ρ(x)(|x|^{2}ω − (ω . x)x) dV,  (3.26)

NST IB Mathematical Methods II 58

or, in suffix notation,

h_{i} = ∫∫∫_{V} ρ(x)(x_{k}x_{k}ω_{i} − ω_{j}x_{j}x_{i}) dV

= ∫∫∫_{V} ρ(x)(x_{k}x_{k}δ_{ij} − x_{j}x_{i})ω_{j} dV

= I_{ij}ω_{j}

where

I_{ij} = ∫∫∫_{V} ρ(x)(x_{k}x_{k}δ_{ij} − x_{i}x_{j}) dV  (3.27)

is the inertia tensor of the rigid body. Note that the tensor I does not depend on ω, only on properties of the body itself; so it may be calculated once and for all for any given body. To see that it is indeed a tensor, note that both h and ω are vectors, and apply arguments previously used for the conductivity tensor.

(iv) Susceptibility χ. If M is the magnetization (magnetic moment per unit volume) and B is the applied magnetic field, then for a simple medium we have M = χ^{(m)}B where χ^{(m)} is the magnetic susceptibility. This generalises to M_{i} = χ_{ij}^{(m)}B_{j} where χ_{ij}^{(m)} is the magnetic susceptibility tensor. Similarly for polarization density P in a dielectric: P_{i} = χ_{ij}^{(e)}E_{j} where E is the electric field and χ_{ij}^{(e)} is the electric susceptibility tensor.

(v) The Kronecker delta itself. We have defined δ_{ij} without reference to frame; i.e., its components are by definition the same in all frames (δ′_{ij} ≡ δ_{ij}). Surprisingly, then, we can show that it is a tensor:

l_{ip}l_{jq}δ_{pq} = l_{ip}l_{jp} = δ_{ij} = δ′_{ij}  (3.28)

(from §3.2), which is exactly the right transformation law. We can also show that ϵ_{ijk} is a tensor of rank 3.

Both δ_{ij} and ϵ_{ijk} are isotropic tensors: that is, their components are the same in all frames.

(vi) Stress and strain tensors. In an elastic body, stresses (forces) result from displacements of small volume elements within the body. Let this displacement at a location x be u; then the strain tensor is defined to be

e_{ij} = \frac{1}{2} (\frac{∂u_{i}}{∂x_{j}} + \frac{∂u_{j}}{∂x_{i}} ) .  (3.29)

The stress tensor p_{ij} is defined as the j^{th} component of the forces within the body acting on an imaginary plane perpendicular to the i^{th} axis. Hooke’s law for simple isotropic media says that stress ∝ strain. We can now generalise this to the tensor formulation

p_{ij} = k_{ijkl}e_{kl}  (3.30)

where k_{ijkl} is a fourth rank tensor, which expresses the linear (but possibly anisotropic) relationship between p and e.

3.4 Properties of Tensors

Linear Combination of Tensors

If A_{ij} and B_{ij} are second rank tensors, and α, β are scalars, then T_{ij} = αA_{ij} + βB_{ij} is a tensor.

Proof:

T′_{ij} = αA′_{ij} + βB′_{ij}
= αl_{ip}l_{jq}A_{pq} + βl_{ip}l_{jq}B_{pq}
= l_{ip}l_{jq}(αA_{pq} + βB_{pq})
= l_{ip}l_{jq}T_{pq}

as required.

This result clearly extends to tensors of rank n.

Contraction (also known as the Inner Product)

If T_{ij} is a tensor then T_{ii} is a scalar. Proof:

T′_{ii} = l_{ip}l_{iq}T_{pq} = δ_{pq}T_{pq} = T_{pp} = T_{ii},  (3.31)

so T_{ii} has the same value in all frames as required.

We can extend this result: if T_{ijk...lmn...} is a tensor of rank N then S_{jk...mn...} = T_{ijk...imn...} is a tensor of rank N − 2. Proof:

S′_{jk...mn...} = T′_{ijk...imn...}
= l_{ip}l_{jq}l_{kr} … l_{iα}l_{mβ}l_{nγ} … T_{pqr...αβγ...}
= (l_{ip}l_{iα})l_{jq}l_{kr} … l_{mβ}l_{nγ} … T_{pqr...αβγ...}
= δ_{pα}l_{jq}l_{kr} … l_{mβ}l_{nγ} … T_{pqr...αβγ...}
= l_{jq}l_{kr} … l_{mβ}l_{nγ} … S_{qr...βγ...}.

Outer Product

If a and b are vectors then the outer product T_{ij} defined by T_{ij} = a_{i}b_{j} is a tensor of rank two. Proof:

T′_{ij} = a′_{i}b′_{j} = l_{ip}a_{p}l_{jq}b_{q} = l_{ip}l_{jq}a_{p}b_{q} = l_{ip}l_{jq}T_{pq}  (3.32)

as required.

Similarly (left as an exercise for the reader) we can show that if A_{ijk...} is a tensor of rank M and B_{lmn...} is a tensor of rank N, then T_{ijk...lmn...} = A_{ijk...}B_{lmn...} is a tensor of rank M + N.

Example: if a and b are vectors then a . b is a scalar. Proof: T_{ij} = a_{i}b_{j}, being an outer product of two vectors, is a tensor of rank two. Then T_{ii} = a_{i}b_{i}, being a contraction of a tensor, is a scalar, as required. Note that |a|^{2} = a . a and |b|^{2} are also scalars; hence a.b/|a| |b| = cos θ is a scalar, so that the angle between vectors is unaffected by a change of frame.

Extra example of outer product of (5,1,2) with (1,2,3)

3.5 Symmetric and Anti-Symmetric Tensors

A tensor T_{ijk...} is said to be symmetric in a pair of indices (say i, j) if

T_{ijk...} = T_{jik...}  (3.33)

or anti-symmetric in i, j if

T_{ijk...} = −T_{jik...}.  (3.34)

For example, δ_{ij} is symmetric; ϵ_{ijk} is anti-symmetric in any pair of indices. Another example: outer product a_{i}a_{j} of vector with itself is obviously symmetric.

Note: if A_{ij} is a symmetric second rank tensor then the matrix corresponding to A is symmetric, i.e. A = A^{T}. Similarly for an anti-symmetric tensor.

Suppose that S_{ij} is a symmetric tensor and A_{ij} an anti-symmetric tensor. Then S_{ij}A_{ij} = 0. Proof:

S_{ij}A_{ij} = −S_{ij}A_{ji} = −S_{ji}A_{ji}
= −S_{ij}A_{ij} (swapping dummy i and j)
⟹  2S_{ij}A_{ij} = 0,

as required. Try to work out also how to see this “by inspection”, by considering appropriate pairs of components.

Example: for any vector a, a × a = 0 because

[a × a]_{i} = ϵ_{ijk}a_{j}a_{k}  (3.35)

and ϵ_{ijk} is anti-symmetric in j, k whilst a_{j}a_{k} is symmetric.

The properties of symmetry and anti-symmetry are invariant under a change of frame: that is, they are truly tensor properties. For example, suppose that T_{ij} is symmetric. Then

T′_{ij} = l_{ip}l_{jq}T_{pq}
= l_{ip}l_{jq}T_{qp}
= l_{jq}l_{ip}T_{qp} = T′_{ji},

so that T′_{ij} is also symmetric.

NST IB Mathematical Methods II 62

(Alternative, and simpler, proof for second rank tensors:

T′ = LT L^{T} ⟹ T′^{T} = (LT L^{T})^{T} = LT^{T} L^{T} = LT L^{T} = T′ (3.36)

using T^{T} = T.)

Symmetry and anti-symmetry occur frequently in practical applications. For example, the strain tensor e_{ij} = \frac{1}{2}(∂u_{i}/∂x_{j} + ∂u_{j}/∂x_{i}) is clearly symmetric. In most situations the stress tensor is also symmetric; but in some circumstances (for instance in crystallography or geodynamics) it is forced to be anti-symmetric while the strain remains symmetric. Inertia tensors are always symmetric; conductivity and susceptibility tensors usually are.

Decomposition into Symmetric and Anti-Symmetric Parts

Any second rank tensor T_{ij} can be uniquely expressed as the sum of a symmetric and an anti-symmetric tensor; for

T_{ij} = S_{ij} + A_{ij}  (3.37)

where

S_{ij} = \frac{1}{2}(T_{ij} + T_{ji}), A_{ij} = \frac{1}{2}(T_{ij} − T_{ji}) (3.38)

are symmetric and anti-symmetric respectively. Exercise: prove that S and A are tensors.

Furthermore, any anti-symmetric tensor A_{ij} can be expressed in terms of a vector ω (sometimes known as the dual vector) such that

A_{ij} = ϵ_{ijk}ω_{k}.  (3.39)

Proof: define ω by

ω_{k} = \frac{1}{2}ϵ_{klm}A_{lm}.  (3.40)

Then

ϵ_{ijk}ω_{k} = \frac{1}{2}ϵ_{ijk}ϵ_{klm}A_{lm}
= \frac{1}{2}(δ_{il}δ_{jm} − δ_{im}δ_{jl})A_{lm}
= \frac{1}{2}(A_{ij} − A_{ji}) = A_{ij}

as required. ω is a tensor as it is a contraction of two tensors.

This definition of ω actually corresponds to setting

A = (\begin{matrix}0 & ω_{3} & −ω_{2} \\ −ω_{3} & 0 & ω_{1} \\ ω_{2} & −ω_{1} & 0 \end{matrix}) .  (3.41)

Worked Example: Decomposition of Second Rank Tensors

Consider an elastic body subjected to a simple shear, so that the displacement u at a location x = (x, y, z) is given by

u = (γy, 0, 0)

for some constant γ. Consider the differential of the displacement, ∂u_{i}/∂x_{j}, which is given by the matrix

(\begin{matrix}0 & γ & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}) .  (3.42)

We can decompose this tensor into symmetric and anti-symmetric parts,

\frac{∂u_{i}}{∂x_{j}} = (\begin{matrix}0 & \frac{1}{2}γ & 0 \\ \frac{1}{2}γ & 0 & 0 \\ 0 & 0 & 0 \end{matrix}) + (\begin{matrix}0 & \frac{1}{2}γ & 0 \\ −\frac{1}{2}γ & 0 & 0 \\ 0 & 0 & 0 \end{matrix}),  (3.43)

in which the symmetric part is just the strain tensor e_{ij}. The anti-symmetric part can also be written in the form ϵ_{ijk}ω_{k} where ω = (0, 0, \frac{1}{2}γ).

This decomposition corresponds to writing

u = (γy, 0, 0) = (\frac{1}{2}γy, \frac{1}{2}γx, 0) + (\frac{1}{2}γy, −\frac{1}{2}γx, 0).

The first term is a stretch at 45^{∘} to the (x, y)-axes, while the second is a rotation. In fact, any vector field u which has zero divergence can be decomposed using this method into a suitable stretch and a solid-body rotation.

Example: suppose that two symmetric second rank tensors R_{ij} and S_{ij} are linearly related. Then there must be a relationship between them of the form R_{ij} = c_{ijkl}S_{kl}. It is clear that c_{ijkl} must be symmetric in i, j (for otherwise, R_{ij} would not be). It is not necessarily the case that it must also be symmetric in k, l, but without loss of generality we may assume that it is, by the following argument. Decompose c_{ijkl} into a part c_{ijkl}^{(s)} which is symmetric in k, l and a part c_{ijkl}^{(a)} which is anti-symmetric. Then

R_{ij} = c_{ijkl}^{(s)}S_{kl} + c_{ijkl}^{(a)}S_{kl} = c_{ijkl}^{(s)}S_{kl}  (3.44)

because the second term is the contraction of an anti-symmetric tensor with a symmetric one, which we showed was zero above. Hence we can ignore any anti-symmetric part of c_{ijkl}.

3.6 Diagonalization of Symmetric Second Rank Tensors

Suppose T_{ij} is a symmetric second rank tensor. We shall show that there exists a frame such that, if we transform T to that frame, it has components given by

T′ = (\begin{matrix}λ_{1} & 0 & 0 \\ 0 & λ_{2} & 0 \\ 0 & 0 & λ_{3} \end{matrix}) .  (3.45)

This process is known as diagonalization. The values λ_{1}, λ_{2} and λ_{3} are known as the principal values of T, and the Cartesian coordinate axes of the corresponding frame are known as the principal axes. We will see that in fact the principal values are just the eigenvalues of the matrix corresponding to T, and the principal axes are the eigenvectors.

Because T is symmetric, we know that there are 3 real eigenvalues and that we can find 3 corresponding eigenvectors which are orthogonal and of unit length. Let λ_{1}, λ_{2}, λ_{3} be the eigenvalues and e′_{1}, e′_{2}, e′_{3} be the eigenvectors (arranged as a right-handed set of orthonormal vectors). Change frame to one in which the coordinate axes are aligned with {e′_{1}, e′_{2}, e′_{3}}. What is T′?

Recall that L^{T} = ( e′_{1} | e′_{2} | e′_{3} ); i.e., the three columns of L^{T} are the vectors e′_{1}, e′_{2} and e′_{3} (measured relative to the first frame). Hence in matrix notation,

T L^{T} = T ( e′_{1} | e′_{2} | e′_{3} )
= ( λ_{1}e′_{1} | λ_{2}e′_{2} | λ_{3}e′_{3} ).

So

T′ = LT L^{T} = ( \begin{matrix}e′_{1} & T \\ e′_{2} & T \\ e′_{3} & T \end{matrix}  ) (λ_{1}e′_{1} λ_{2}e′_{2} λ_{3}e′_{3})

= (\begin{matrix}λ_{1} & 0 & 0 \\ 0 & λ_{2} & 0 \\ 0 & 0 & λ_{3} \end{matrix})

because, for example, the top LHS entry is given by e′_{1} . λ_{1}e′_{1}, and the top RHS entry is e′_{1} . λ_{3}e′_{3}.

There is another way of seeing that T′ = (\begin{matrix}λ_{1} & 0 & 0 \\ 0 & λ_{2} & 0 \\ 0 & 0 & λ_{3} \end{matrix}). The equation T e′_{1} = λ_{1}e′_{1} is true in any frame

(because T is a tensor, e′_{1} a vector and λ_{1} a scalar). In particular it is true in the frame with {e′_{1}, e′_{2}, e′_{3}} as coordinate axes. But, measured in this frame, e′_{1} is just (1, 0, 0)^{T}, and T has components T′; so

T′ (\begin{matrix}1 \\ 0 \\ 0 \end{matrix}) = (\begin{matrix}λ_{1} \\ 0 \\ 0 \end{matrix})  (3.46)

which shows that the first column of T′ is (λ_{1}, 0, 0)^{T}. Similarly for the other columns.

Note: the three principal values are invariants of T. That is, whatever frame we start from, when we diagonalize T we will obtain the same values of λ. The eigenvalues are properties of the tensor, not of the coordinate system.

3.7 Isotropic Tensors

An isotropic tensor is one whose components are the same in all frames, i.e.,

T′_{ijk...} = T_{ijk...}.  (3.47)

We can classify isotropic tensors up to rank four as follows:

Rank 0: All scalars are isotropic, since the tensor transformation law states that T′ = T for tensors of rank zero.

NST IB Mathematical Methods II 66

Rank 1: There are no non-zero isotropic vectors. Rank 2: The most general isotropic second rank tensor is λδ_{ij} where λ is any scalar,

as proved below. Rank 3: The most general isotropic third rank tensor is λϵ_{ijk}. Rank 4: The most general isotropic fourth rank tensor is

λδ_{ij}δ_{kl} + μδ_{ik}δ_{jl} + νδ_{il}δ_{jk}  (3.48)

where λ, μ, ν are scalars.

What is the physical significance of an isotropic tensor? Consider the conductivity tensor σ_{ij} in an isotropic medium. As the medium is the same in all directions, we expect that σ_{ij} will be isotropic too. Hence σ_{ij} = λδ_{ij} and

J_{i} = σ_{ij}E_{j} = λδ_{ij}E_{j} = λE_{i},  (3.49)

i.e., J = λE. So we recover the “simple version” of the conductivity law, as we might expect.

Isotropic Second Rank Tensors

Consider a general tensor T of rank two, with components T_{ij} with respect to some set of axes {e_{1}, e_{2}, e_{3}}. Suppose that T is isotropic. Its components should then be unaltered by a rotation of 90^{∘} about the 3-axis, i.e., with respect to new axes

e′_{1} = e_{2}, e′_{2} = −e_{1}, e′_{3} = e_{3}. (3.50)

The matrix of this rotation is

L = (\begin{matrix}0 & 1 & 0 \\ −1 & 0 & 0 \\ 0 & 0 & 1 \end{matrix}) .  (3.51)

Using the matrix formulation of the transformation law for tensors, we see that

(\begin{matrix}T′_{11} & T′_{12} & T′_{13} \\ T′_{21} & T′_{22} & T′_{23} \\ T′_{31} & T′_{32} & T′_{33} \end{matrix}) = (\begin{matrix}0 & 1 & 0 \\ −1 & 0 & 0 \\ 0 & 0 & 1 \end{matrix}) (\begin{matrix}T_{11} & T_{12} & T_{13} \\ T_{21} & T_{22} & T_{23} \\ T_{31} & T_{32} & T_{33} \end{matrix}) (\begin{matrix}0 & −1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{matrix})

= (\begin{matrix}T_{22} & −T_{21} & T_{23} \\ −T_{12} & T_{11} & −T_{13} \\ T_{32} & −T_{31} & T_{33} \end{matrix}) .

But, because T is isotropic, T′_{ij} = T_{ij}. Hence, comparing matrix entries, we have:

T_{11} = T_{22};
T_{13} = T_{23} = −T_{13} so that T_{13} = T_{23} = 0;
T_{31} = T_{32} = −T_{31} so that T_{31} = T_{32} = 0.

Similarly, considering a rotation of 90^{∘} about the 2-axis, we find that T_{11} = T_{33} and that T_{12} = T_{32} = 0, T_{21} = T_{23} = 0. Therefore all off-diagonal elements of T are zero, and all diagonal elements are equal, say λ. Thus

T = (\begin{matrix}λ & 0 & 0 \\ 0 & λ & 0 \\ 0 & 0 & λ \end{matrix}) ,  (3.52)

or in suffix notation, T_{ij} = λδ_{ij}.

In summary, we have shown that any isotropic second rank tensor must be equal to λδ_{ij} for some scalar λ.

Worked Example: Evaluation of an Isotropic Integral

We wish to calculate

T_{ij} =∫ ∫ ∫_{All space}x_{i}x_{j}e^{−r^{2}}dV

for each value of i and j.

There are no special directions involved either in the domain of integration or in the integrand; so T must be isotropic. Hence T_{ij} = λδ_{ij} for some λ. To calculate λ, consider NST IB Mathematical Methods II 68

T_{ii} = λδ_{ii} = 3λ. But we know that

T_{ii} = ∫∫∫_{ℝ^{3}} x_{i}x_{i}e^{−r^{2}}dV  (3.53)
= ∫∫∫_{ℝ^{3}} r^{2}e^{−r^{2}}dV  (3.54)
= ∫_{0}^{2π}∫_{0}^{π}∫_{0}^{∞} r^{2}e^{−r^{2}}r^{2} sin θ drdθdϕ  (3.55)
= 4π ∫_{0}^{∞} r^{4}e^{−r^{2}}dr  (3.56)
= 4π(\frac{3}{8}\sqrt{π} ).

Hence we conclude that

∫∫∫_{ℝ^{3}} x_{i}x_{j}e^{−r^{2}}dV = \frac{1}{2}π\sqrt{π} δ_{ij}.

Such calculations are often of use when a physical situation has symmetry which can be exploited; for example, consider calculating the inertia tensor of a sphere.

3.8 Tensor Differential Operators

A tensor field is a tensor which depends on the location x. For example:

(i) Temperature is a scalar field (a tensor field of rank zero), because T = T (x).

(ii) Any vector field F(x), such as a gravitational force field, is a tensor field of rank one. In particular, x is itself a vector field, because it is a vector function of position!

(iii) In a conducting material where the conductivity varies with location, we have σ_{ij} = σ_{ij}(x), a tensor field of rank two.

NST IB Mathematical Methods II 69

We are interested here in calculating the derivatives of tensor fields; we start with scalars and vectors.

Recall that grad, div, and curl can be written using suffix notation:

Grad:  [∇Φ]_{i} = \frac{∂Φ}{∂x_{i}}  (3.57)
Div:  ∇. F = \frac{∂F_{1}}{∂x_{1}} + \frac{∂F_{2}}{∂x_{2}} + \frac{∂F_{3}}{∂x_{3}} = \frac{∂F_{i}}{∂x_{i}}  (3.58)
Curl: [∇ × F]_{i} = ϵ_{ijk} \frac{∂F_{k}}{∂x_{j}}  (3.59)

There is another useful notation: if u, v are vectors then we define the vector

(u . ∇)v = (u_{1} \frac{∂}{∂x_{1}} + u_{2} \frac{∂}{∂x_{2}} + u_{3} \frac{∂}{∂x_{3}} ) v.  (3.60)

In suffix notation,

[(u . ∇)v]_{i} = u_{j} \frac{∂v_{i}}{∂x_{j}} .  (3.61)

Laplace’s equation ∇^{2}Φ = 0 becomes

\frac{∂^{2}Φ}{∂x_{i}∂x_{i}} = 0  (3.62)

in suffix notation. Similarly,

[∇^{2}F]_{i} = \frac{∂^{2}F_{i}}{∂x_{j}∂x_{j}}  (3.63)

(note that we only use Cartesian coordinates here).

Worked Example: Proving Vector Differential Identities

To prove that ∇. (Φu) = u . ∇Φ + Φ ∇. u where Φ is a scalar field and u is a vector field:

∇. (Φu) = \frac{∂}{∂x_{i}} (Φu_{i})  (3.64)
= \frac{∂Φ}{∂x_{i}} u_{i} + Φ\frac{∂u_{i}}{∂x_{i}}  (3.65)
= [∇Φ]_{i}u_{i} + Φ\frac{∂u_{i}}{∂x_{i}}  (3.66)
= u . ∇Φ + Φ ∇. u.

To prove that ∇ × (u × v) = (∇. v)u − (∇. u)v + (v . ∇)u − (u . ∇)v where u and v are vector fields:

[∇ × (u × v)]_{i} = ϵ_{ijk} \frac{∂}{∂x_{j}} ϵ_{klm}u_{l}v_{m}
= (δ_{il}δ_{jm} − δ_{im}δ_{jl})\frac{∂}{∂x_{j}} (u_{l}v_{m})
= \frac{∂}{∂x_{j}} (u_{i}v_{j}) − \frac{∂}{∂x_{j}} (u_{j}v_{i})
= u_{i} \frac{∂v_{j}}{∂x_{j}} + v_{j} \frac{∂u_{i}}{∂x_{j}} − v_{i} \frac{∂u_{j}}{∂x_{j}} − u_{j} \frac{∂v_{i}}{∂x_{j}}
= [(∇. v)u + (v . ∇)u − (∇. u)v − (u . ∇)v]_{i}.

We sometimes find it useful to use the differential operator ∂_{i} defined by

∂_{i} = \frac{∂}{∂x_{i}} .  (3.67)

Then

[∇Φ]_{i} = ∂_{i}Φ;  ∇. F = ∂_{i}F_{i};  [∇ × F]_{i} = ϵ_{ijk}∂_{j}F_{k}.  (3.68)

It turns out that ∂_{i} is in fact a tensor of rank one. We know that x_{j} = l_{ij}x′_{i} (from x = L^{T} x′) so that

\frac{∂x_{j}}{∂x′_{i}} = \frac{∂}{∂x′_{i}} (l_{kj}x′_{k}) = l_{kj} \frac{∂x′_{k}}{∂x′_{i}} = l_{kj}δ_{ik} = l_{ij}.  (3.69)

(This looks obvious but has to be proved very carefully!) Now let T be some quantity (perhaps a scalar or a tensor of some rank). Then

∂′_{i}T = \frac{∂T}{∂x′_{i}} = \frac{∂T}{∂x_{1}} \frac{∂x_{1}}{∂x′_{i}} + \frac{∂T}{∂x_{2}} \frac{∂x_{2}}{∂x′_{i}} + \frac{∂T}{∂x_{3}} \frac{∂x_{3}}{∂x′_{i}}

= \frac{∂T}{∂x_{j}} \frac{∂x_{j}}{∂x′_{i}}

= l_{ij} \frac{∂T}{∂x_{j}} = l_{ij}∂_{j}T.

This is true for any quantity T, so

∂′_{i} = l_{ij}∂_{j},  (3.70)

i.e., ∂_{i} transforms like a vector, and is hence a tensor of rank one.

NST IB Mathematical Methods II 71

This result allows us to prove that ∇Φ, ∇ . F and ∇ × F are scalars or vectors (as appropriate). For example, to show that if F is a vector field then ∇ × F is a vector field:

[∇ × F]′_{i} = ϵ′_{ijk}∂′_{j}F′_{k}

= l_{ip}l_{jq}l_{kr}ϵ_{pqr}l_{js}∂_{s}l_{kt}F_{t}

[ϵ, ∂ and F are tensors] = l_{ip}(l_{jq}l_{js})(l_{kr}l_{kt})ϵ_{pqr}∂_{s}F_{t} = l_{ip}δ_{qs}δ_{rt}ϵ_{pqr}∂_{s}F_{t} = l_{ip}ϵ_{pqr}∂_{q}F_{r} = l_{ip}[∇ × F]_{p},

as required.

Alternatively, we can just state that ∇×F is a contraction of the tensor outer product T_{ijklm} = ϵ_{ijk}∂_{l}F_{m} (because [∇ × F]_{i} = T_{ijkjk}).

As an example of a tensor field of rank three, consider the derivative of the conductivity tensor, ∂_{i}σ_{jk}. This cannot be written using ∇.

NST IB Mathematical Methods II 72

Complex Analysis

4.1 Complex Differentiation

Recall the definition of differentiation for a real function f(x):

f′(x) = lim_{δx→0} \frac{f(x + δx) − f(x)}{δx}.  (4.1)

In this definition, it is important that the limit is the same whichever direction we approach from. Consider |x| at x = 0 for example; if we approach from the right (δx → 0^{+}) then the limit is +1, whereas if we approach from the left (δx → 0^{−}) the limit is −1. Because these limits are different, we say that |x| is not differentiable at x = 0.

Now extend the definition to complex functions f(z):

f′(z) = lim_{δz→0} \frac{f(z + δz) − f(z)}{δz}.  (4.2)

Again, the limit must be the same whichever direction we approach from; but now there is an infinity of possible directions.

Definition: if f′(z) exists and is continuous in some region R of the complex plane, we say that f is analytic in R. If f(z) is analytic in some small region around a point z_{0}, then we say that f(z) is analytic at z_{0}. The term regular is also used instead of analytic.

Note: the property of analyticity is in fact a surprisingly strong one! For example, two consequences include:

(i) If a function is analytic then it is differentiable infinitely many times. (Cf. the existence of real functions which can be differentiated N times but no more, for any given N.)

(ii) If a function is analytic and bounded in the whole complex plane, then it is constant. (Liouville’s Theorem.)

The Cauchy–Riemann Equations

Let f(z) be a differentiable complex-valued function of complex variable z. Separate f and z into real and imaginary parts:

f(z) = u(x, y) + iv(x, y)  (4.3)

where z = x + iy and u, v are real functions. Suppose that f is differentiable at z. We can take δz in any direction; first take it to be real, δz = δx. Then

f′(z) = lim_{δx→0} \frac{f(z + δx) − f(z)}{δx}

= lim_{δx→0} \frac{u(x + δx, y) + iv(x + δx, y) − u(x, y) − iv(x, y)}{δx}

= lim_{δx→0} \frac{u(x + δx, y) − u(x, y)}{δx} + i lim_{δx→0} \frac{v(x + δx, y) − v(x, y)}{δx}

= \frac{∂u}{∂x} + i\frac{∂v}{∂x}.

Now take δz to be pure imaginary, δz = iδy. Then

f′(z) = lim_{δy→0} \frac{f(z + iδy) − f(z)}{iδy}

= lim_{δy→0} \frac{u(x, y + δy) + iv(x, y + δy) − u(x, y) − iv(x, y)}{iδy}

= −i\frac{∂u}{∂y} + \frac{∂v}{∂y} .

The two values for f′(z) are the same since f is differentiable, so

\frac{∂u}{∂x} + i\frac{∂v}{∂x} = \frac{∂v}{∂y} − i\frac{∂u}{∂y}

⟹  \frac{∂u}{∂x} = \frac{∂v}{∂y} ,  \frac{∂u}{∂y} = −\frac{∂v}{∂x}

– the Cauchy–Riemann equations. It is also possible to show that if the Cauchy–Riemann equations hold at a point z, then f is differentiable there (subject to certain technical conditions on the continuity of the partial derivatives).

If we know the real part u of an analytic function, the Cauchy–Riemann equations allow us to find the imaginary part v (up to a constant), and vice versa. For example, if u(x, y) = x^{2} − y^{2} then

\frac{∂v}{∂y} = \frac{∂u}{∂x} = 2x ⟹  v = 2xy + g(x)  (4.4)

for some function g(x); so

−2y = \frac{∂u}{∂y} = −\frac{∂v}{∂x} = −2y − g′(x) ⟹ g′(x) = 0 ⟹ g = const. = α, say. (4.5)

Hence

f(z) = x^{2} − y^{2} + 2ixy + iα = (x + iy)^{2} + iα = z^{2} + iα.  (4.6)

Examples of Analytic Functions

(i) f(z) = z is analytic in the whole of ℂ. Here u = x, v = y, and the Cauchy–Riemann equations are satisfied (1 = 1; 0 = 0).

(ii) f(z) = z^{n} (n a positive integer) is analytic in ℂ. Here we write z = r(cos θ + i sin θ) and by de Moivre’s theorem, z^{n} = r^{n}(cos nθ + i sin nθ). Hence u = r^{n} cos nθ and v = r^{n} sin nθ: we can check the Cauchy–Riemann equations (using r = \sqrt{x^{2} + y^{2}}, θ = tan^{−1}(y/x)). The derivative is nz^{n−1}, as expected.

(iii) f(z) = e^{z} = e^{x}e^{iy} = e^{x}(cos y + i sin y). So

\frac{∂u}{∂x} = e^{x} cos y = \frac{∂}{∂y} (e^{x} sin y) = \frac{∂v}{∂y};  \frac{∂u}{∂y} = −e^{x} sin y = −\frac{∂v}{∂x}.  (4.7)

The derivative is

f′(z) = \frac{∂u}{∂x} + i\frac{∂v}{∂x} = e^{x} cos y + ie^{x} sin y = e^{z},  (4.8)

again as expected.

(iv) f(z) = 1/z: check that this is analytic with derivative −1/z^{2} in any region R which does not include the origin.

(v) More generally, any rational function – i.e., f(z) = P (z)/Q(z) where P and Q are polynomials – is analytic except at points where Q(z) = 0. For instance, f(z) = (z − i)/(z + i) is analytic except at z = −i.

(vi) Many standard functions obey the usual rules for their derivatives; e.g.,

\frac{d}{dz} sin z = cos z,  \frac{d}{dz} sinh z = cosh z,
\frac{d}{dz} cos z = − sin z,  \frac{d}{dz} cosh z = sinh z,
\frac{d}{dz} log z = \frac{1}{z} (when log z is defined as later).

The product, quotient and chain rules hold in exactly the same way as for real functions.

Examples of Non-Analytic Functions

(i) f(z) = Re(z). Here u = x, v = 0, but 1 ≠ 0. Re(z) is nowhere analytic.

(ii) f(z) = |z|; here u = \sqrt{x^{2} + y^{2}}, v = 0.

(iii) f(z) = |z|^{2} = x^{2} + y^{2}. The Cauchy–Riemann equations are only satisfied at the origin, so f is only differentiable at z = 0. However, it is not analytic there because there is no small region containing the origin within which f is differentiable.

(iv) f(z) = z̄ = x − iy (complex conjugate, also denoted z^{∗}). Here u = x, v = −y, so ∂u/∂x = 1 ≠ −1 = ∂v/∂y.

Harmonic Functions

Suppose f(z) = u + iv is analytic. Then

\frac{∂^{2}u}{∂x^{2}} = \frac{∂}{∂x} (\frac{∂u}{∂x}) = \frac{∂}{∂x} (\frac{∂v}{∂y} )

= \frac{∂}{∂y} (\frac{∂v}{∂x}) = \frac{∂}{∂y} (−\frac{∂u}{∂y}) = −\frac{∂^{2}u}{∂y^{2}}.

Hence

\frac{∂^{2}u}{∂x^{2}} + \frac{∂^{2}u}{∂y^{2}} = 0,  (4.9)

i.e., u satisfies Laplace’s equation in two dimensions. Similarly, v does too. Such functions u and v are said to be harmonic.

4.2 Zeros of Complex Functions

The zeros of f(z) are the points z_{0} where f(z_{0}) = 0. A zero is of order n if

0 = f′(z_{0}) = f′′(z_{0}) = ⋯ = f^{(n−1)}(z_{0}), but f^{(n)}(z_{0}) ≠ 0.  (4.10)

A zero of order one (i.e., one where f′(z_{0}) ≠ 0) is called a simple zero. Examples:

(i) f(z) = z has a simple zero at z = 0.

(ii) f(z) = (z − i)^{2} has a zero of order two at z = i.

(iii) f(z) = z^{2}−1 = (z−1)(z+1) has two simple zeros at z = ±1.

(iv) f(z) = (z − w)^{N}g(z), where w is a complex constant, N a positive integer and g(z) an analytic function satisfying g(w) ≠ 0, has a zero of order N at z = w.

(v) Where are the zeros of f(z) = sinh z? We know there is a simple zero at z = 0. The others are where

0 = sinh z = \frac{e^{z} − e^{−z}}{2} ⟺ e^{z} = e^{−z} ⟺ e^{2z} = 1 ⟺ 2z = 2nπi,

(4.11)

where n is an integer. (Check that e^{x+iy} = 1 ⇔ x = 0 and y = 2nπ.) So the zeros are on the imaginary axis at …, −2πi, −πi, 0, πi, 2πi, 3πi, …, and they are all simple.

Another way of defining the order of a zero is by the first non-zero power of (z − z_{0}) in its Taylor series. For example, consider the zero of sinh^{3} z at z = πi. Now sinh z = − sinh(z − πi) = − sinh ζ where ζ = z − πi, and close to z = πi the Taylor series for sinh z is therefore

−(ζ + \frac{1}{3!}ζ^{3} + ⋯).  (4.12)

Hence the Taylor series for sinh^{3} z at z = πi is

−(ζ + \frac{1}{3!}ζ^{3} + ⋯)^{3} = −(z − πi)^{3} − \frac{1}{2}(z − πi)^{5} + · · · .  (4.13)

The zero is therefore of order 3.

4.3 Laurent Expansions

Suppose that f(z) is analytic at z_{0}. Then we can expand f in a Taylor Series about z_{0}:

f(z) = ∑_{n=0}^{∞} a_{n}(z − z_{0})^{n}  (4.14)

NST IB Mathematical Methods II 78

for suitable complex constants a_{n}.

Example: e^{z} has a Taylor Series about z = i given by

e^{z} = e^{i}e^{z−i} = e^{i} ∑_{n=0}^{∞} \frac{(z − i)^{n}}{n!} ,  (4.15)

so a_{n} = e^{i}/n!.

Now consider an f(z) which is not analytic at z_{0}, but for which (z−z_{0})f(z) is analytic. (E.g., f(z) = e^{z}/(z − z_{0}).) Then, for suitable b_{n},

(z − z_{0})f(z) = ∑_{n=0}^{∞} b_{n}(z − z_{0})^{n}
⟹ f(z) = \frac{b_{0}}{z − z_{0}} + b_{1} + b_{2}(z − z_{0}) + b_{3}(z − z_{0})^{2} + ⋯

= ∑_{n=−1}^{∞} a_{n}(z − z_{0})^{n}

where a_{n} = b_{n+1}. Generalising this, if (z − z_{0})^{N}f(z) is analytic at z_{0} then for suitable a_{n},

f(z) = ∑_{n=−N}^{∞} a_{n}(z − z_{0})^{n}.  (4.16)

But what if however large N is, (z−z_{0})^{N}f(z) is still not analytic at z_{0}? We might try taking N to be “infinite”, and in fact this does always work.

Formally, it is possible to show that if f(z) is analytic in an annulus a < |z − z_{0}| < b for some a, b (regardless of whether f is analytic at z_{0} itself) then f has a unique Laurent expansion

f(z) = ∑_{n=−∞}^{∞} a_{n}(z − z_{0})^{n}  (4.17)

in the annulus.

Examples:

(i) f(z) = e^{z}/z^{3} = ∑_{n=0}^{∞} z^{n−3}/n! about z_{0} = 0 has a_{n} = 0 for n < −3 and a_{n} = 1/(n + 3)! for n ≥ −3.

(ii) f(z) = e^{z}/(z^{2} − 1) about z_{0} = 1 (where it has a singularity). Here we write everything in terms of ζ = z − z_{0}, so

f(z) = \frac{e^{ζ}e^{z_{0}}}{ζ(ζ + 2)} = \frac{e^{z_{0}}}{2ζ} e^{ζ}(1 + \frac{1}{2}ζ)^{−1}

= \frac{e}{2ζ} (1 + ζ + \frac{1}{2!}ζ^{2} + ⋯)(1 − \frac{1}{2}ζ + ⋯)

= \frac{e}{2ζ} (1 + \frac{1}{2}ζ + ⋯)

= \frac{e}{2} (\frac{1}{z − z_{0}} + \frac{1}{2} + ⋯) .

Hence a_{−1} = e/2, a_{0} = e/4, etc.

(iii) f(z) = exp(1/z) about z_{0} = 0 has

e^{1/z} = 1 + \frac{1}{z} + \frac{1}{2! z^{2}} + \frac{1}{3! z^{3}} + ⋯,  (4.18)

so that a_{n} = 1/(−n)! for n ≤ 0.

(iv) This doesn’t seem to work for f(z) = z^{− /_{2}^{1}} – why? We shall see later that it is impossible to find an annulus around z_{0} = 0 in which z^{− /_{2}^{1}} is analytic.

If f(z) is in fact analytic at z = z_{0}, then its Laurent expansion about z_{0} is just its Taylor series.

4.4 Classification of Singularities

Suppose that f has a singularity at z = z_{0}, but is analytic within some circle |z − z_{0}| < r except at z_{0} itself. Such a singularity is called an isolated singularity. Choosing any annulus inside the circle (e.g., r/2 < |z − z_{0}| < r), we see that f has a Laurent expansion about z_{0}. We can use the coefficients a_{n} of the expansion to classify the singularity; there are three possible cases.

Essential Isolated Singularities

If there is no integer N such that a_{n} = 0 for all n < −N – i.e., if however far n goes towards −∞ there are always some non-zero a_{n}’s – then f is said to have an essential isolated singularity. Examples:

(i) exp(1/z) has an essential isolated singularity at z = 0, because all the a_{n}’s are non-zero for n ≤ 0 (we showed above that a_{n} = 1/(−n)!).

(ii) sin(1/z) also has an essential isolated singularity at z = 0, because

a_{n} = \begin{cases}(−1)^{(n+1)/2}/(−n)! & n negative and odd, \\ 0 & n positive or even. \end{cases}  (4.19)

However negative n is, there are some non-zero a_{n}’s for still more negative n.

Near an essential isolated singularity of a function f(z), it can be shown that f takes all possible complex values (bar at most one). For example, sin(1/z) takes all possible complex values near the origin; exp(1/z) takes all except zero.

Poles

If a_{n} = 0 for all n < −N (where N is some specific positive integer) but a_{−N} ≠ 0, then f is said to have a pole of order N. (If N = 1, then we call this a simple pole.) This is the most common, and the most important, of the three cases. Examples:

(i) 1/(z − i) has a simple pole at z = i.

(ii) (cos z)/z (which has Laurent expansion z^{−1} − \frac{1}{2!}z + \frac{1}{4!}z^{3} − ⋯) has a simple pole at z = 0.

(iii) 1/(z^{2} + 1) has two simple poles, at z = ±i.

(iv) z^{2}/{(z − 3)^{3}(z − i)^{2}} has a pole of order 2 at z = i and a pole of order 3 at z = 3. To show formally that, for instance, there is a pole of order 2 at z = i, notice first that z^{2}/(z − 3)^{3} is analytic there so has a Taylor series

b_{0} + b_{1}(z − i) + b_{2}(z − i)^{2} + · · · .  (4.20)

Hence

\frac{z^{2}}{(z − 3)^{3}(z − i)^{2}} = \frac{b_{0}}{(z − i)^{2}} + \frac{b_{1}}{(z − i)} + b_{2} + · · · .  (4.21)

(v) If g(z) has a zero of order m at z = z_{0} then 1/g(z) has a pole of order m there (and vice versa). Hence cot z has a simple pole at the origin, because tan z has a simple zero there. To prove this, note that g(z) = (z − z_{0})^{m}G(z) for some function G(z) satisfying G(z_{0}) ≠ 0. The function 1/G(z) is analytic at z = z_{0}, so it has a Taylor series c_{0} + c_{1}(z − z_{0}) + c_{2}(z − z_{0})^{2} + · · ·. Hence

1/g(z) = c_{0}(z − z_{0})^{−m} + c_{1}(z − z_{0})^{−m+1} + c_{2}(z − z_{0})^{−m+2} + ⋯  (4.22)

as required.

Removable Singularities

If a_{n} = 0 for all n < 0 (so that the Laurent expansion is just a_{0} + a_{1}(z − z_{0}) + ⋯), then f is said to have a removable singularity. By redefining f(z_{0}) = a_{0} we can remove the singularity completely. Examples:

(i) (Somewhat contrived.)

f(z) = \begin{cases}0 & z = 0 \\ 1 & z ≠ 0 \end{cases}  (4.23)

has a singularity at z = 0. Because the origin is not part of any annulus around itself, so that f(z) = 1 everywhere in the annulus, the Laurent expansion has a_{0} = 1 and all other a_{n} = 0, so f has a removable singularity; by redefining f(0) = 1 we obtain an analytic function.

(ii) f(z) = (sin z)/z is not defined at z = 0, but by defining f(0) = 1 we obtain an analytic function.

(iii) A rational function f(z) = P (z)/Q(z) (where P and Q are polynomials) has a singularity at any point z_{0} where Q(z_{0}) = 0; but if P (z_{0}) = 0 as well then the singularity is removable by redefining f(z_{0}) = P′(z_{0})/Q′(z_{0}), assuming that Q′(z_{0}) ≠ 0.

Note: the singularity of z^{−1/2} at z = 0 is not classified under this scheme (though it might look as if it has a pole), as it has no Laurent expansion. (See §4.7 for further explanation: in fact z = 0 is a branch point singularity.)

4.5 Residues

We shall see in Chapter 5 that it is important to be able to calculate the coefficient a_{−1} of the Laurent expansion of a function f(z) about a pole at z_{0}. This coefficient is called the residue of the pole, which we shall denote by res_{z=z_{0}} f(z).

At a simple pole, the residue is given by a_{−1} = lim_{z→z_{0}}{(z − z_{0})f(z)}, because:

lim_{z→z_{0}}\{(z − z_{0})f(z)\} = lim_{z→z_{0}} \{(z − z_{0}) (\frac{a_{−1}}{z − z_{0}} + a_{0} + a_{1}(z − z_{0}) + ⋯)\}

= lim_{z→z_{0}}{a_{−1} + a_{0}(z − z_{0}) + a_{1}(z − z_{0})^{2} + ⋯ }

= a_{−1}.

In general, at a pole of order N, the residue is given by the useful formula

a_{−1} = lim_{z→z_{0}} \{\frac{1}{(N − 1)!} \frac{d^{N−1}}{dz^{N−1}} ((z − z_{0})^{N}f(z))\}  (4.24)

which can easily be proved by first writing down the Laurent expansion of f(z) and then evaluating the right-hand side of the above formula.

Worked Example: Calculating Residues

Example: e^{z}/z^{3} By expanding e^{z} as a Taylor series, we see that f(z) = e^{z}/z^{3} has a Laurent expansion about z = 0 given by

z^{−3} + z^{−2} + \frac{1}{2}z^{−1} + \frac{1}{3!} + · · · .

Hence the residue is \frac{1}{2} (the coefficient of z^{−1}).

Alternatively, we note that f has a pole of order 3 at z = 0, so we can use the general formula for the residue at a pole:

res_{z=0} f(z) = lim_{z→0}\{\frac{1}{2!} \frac{d^{2}}{dz^{2}} (z^{3}f(z))\} = \frac{1}{2} lim_{z→0}\{\frac{d^{2}}{dz^{2}} e^{z}\} = \frac{1}{2}.

Example: e^{z}/(z^{2} − 1) We have already calculated the Laurent expansion of g(z) = e^{z}/(z^{2} − 1) at z = 1:

\frac{e^{z}}{z^{2} − 1} = \frac{e}{2}(\frac{1}{z − 1} + \frac{1}{2} + ⋯),

so the residue is e/2.

Alternatively, we use the formula for the residue at a simple pole:

res_{z=1} g(z) = lim_{z→1} \frac{(z − 1)e^{z}}{z^{2} − 1} = lim_{z→1} \frac{e^{z}}{z + 1} = e/2.

Example: 1/(z^{8} − w^{8}) For any complex constant w, h(z) = (z^{8} − w^{8})^{−1} has 8 simple poles, at z = we^{nπi/4} (n = 0, 1, …, 7). The residue at z = w, say, could be evaluated by factorizing z^{8} − w^{8} into its eight linear factors, but is most easily calculated using L’Hôpital’s Rule:

res_{z=w} h(z) = lim_{z→w} \frac{z − w}{z^{8} − w^{8}} = lim_{z→w} \frac{1}{8z^{7}} = 1/8w^{7}.

Example: 1/ sinh πz 1/ sinh πz has a simple pole at z = ni for all integers n (because the zeros of sinh z are at nπi and are simple). We could use the Taylor series for sinh πz, or the general residue formula: again using L’Hôpital’s Rule, the residue is

lim_{z→ni} \frac{z − ni}{sinh πz} = lim_{z→ni} \frac{1}{π cosh πz} = \frac{1}{π cosh nπi} = \frac{1}{π cos nπ} = (−1)^{n}/π.

Example: 1/ sinh^{3} z

We have seen that sinh^{3} z has a zero of order 3 at z = πi, with Taylor series

sinh^{3} z = −(z − πi)^{3} − \frac{1}{2}(z − πi)^{5} + · · · .

Therefore

1/ sinh^{3} z = −(z − πi)^{−3}(1 + \frac{1}{2}(z − πi)^{2} + ⋯)^{−1}
= −(z − πi)^{−3}(1 − \frac{1}{2}(z − πi)^{2} + ⋯)
= −(z − πi)^{−3} + \frac{1}{2}(z − πi)^{−1} + · · · .

The residue is therefore \frac{1}{2}.

4.6 The Point at Infinity

In the complex plane, we can reach the “point at infinity” by going off in any direction. Conceptually, we may use the Riemann Sphere, which is a sphere resting on the complex plane with its “South Pole” at z = 0.

For any point in ℂ, drawing a line through the “North Pole” of the sphere to the point, and noting where this line intersects the sphere, specifies an equivalent point on the sphere. Then the point at infinity is equivalent to the “North Pole” of the sphere itself.

We can use all the concepts of §4.4 on the point at infinity by using the transformation ζ = 1/z. Let g(ζ) = f(1/ζ) and find the Laurent expansion of g at ζ = 0. Any type of singularity that g has at ζ = 0 is also said to apply to f “at infinity”. Examples:

(i) f(z) = z^{n} has a pole of order n at ∞ (because g(ζ) = f(1/ζ) = ζ^{−n} which has a pole of order n at ζ = 0).

(ii) f(z) = e^{z} has an essential singularity at ∞.

The residue at infinity of f is, similarly, defined to be the residue of g(ζ) at ζ = 0; so, for example, f(z) = 2z has a simple pole at ∞ (because g(ζ) = 2/ζ) with residue 2.

4.7 Multi-Valued Functions

In this section we shall use as our archetypal example log z, the natural logarithm function. For a complex number z = re^{iθ}, we define log z = log r + iθ. There are thus infinitely many values, or “branches”, of log z, for θ may take an infinity of values. For example,

log i = \frac{πi}{2} or \frac{5πi}{2} or − \frac{3πi}{2} or…,  (4.25)

depending on which choice of θ we make.

Branch Points

Consider the three curves shown in the diagram. On C_{1}, we could choose θ to be always in the range (0, \frac{π}{2}), and then log z would be continuous and single-valued going round C_{1}. On C_{2}, we could choose θ ∈ (\frac{π}{2}, \frac{3π}{2}) and log z would again be continuous and single-valued. But for C_{3}, which encircles the origin, there is no such choice; whatever we do, log z cannot be made continuous around C_{3} (it must either “jump” somewhere or be multi-valued). A branch point of a function – here, the origin – is a point which it is impossible to encircle with a curve upon which the function is continuous and single-valued. The function is said to have a branch point singularity at that point.

Examples:

(i) log(z − a) has a branch point at z = a.

(ii) log(z^{2} − 1) = log(z + 1) + log(z − 1) has two branch points, at ±1.

(iii) z /_{2}^{1} = \sqrt{r} e^{iθ/2} has a branch point at the origin. (Useful exercise: verify this.) The same is true of z^{α} = r^{α}e^{iαθ} when α is any non-integer.

Branch Cuts

If we wish to make log z continuous and single-valued, therefore, we must stop any curve from encircling the origin. We do this by introducing a branch cut from −∞ on the real axis to the origin. No curve is allowed to cross this cut. We can then decide to fix on values of θ lying in the range −π < θ ≤ π only, and we have defined a branch of log z which is single-valued and continuous on any curve (which doesn’t cross the cut). This branch is analytic everywhere (with \frac{d}{dz} log z = 1/z) except on the negative real axis.

(If a curve did cross the cut, from above it to below it say, then θ would jump from +π to −π, and log z would be discontinuous and therefore not analytic. This is not allowed.)

This cut is the canonical (i.e., “standard”) branch cut for log z, and the resulting value of log z is called the principal value of the logarithm.

What are the values of log z just above and below the branch cut? Consider a point on the negative real axis, z = x, x < 0. Just above the cut, at z = x + i0^{+}, θ = +π, so log z = log |x| + iπ. Just below it, at z = x + i0^{−}, log z = log |x| − iπ.

Note that many different branch cuts are possible: any cut which stops curves wrapping round the branch point will do. In diagram (a), we could choose −3π/2 < θ ≤ π/2; the exact choice is more difficult to write down in case (b), but this is an equally valid cut.

Exactly the same choices of branch cut can be made for z^{α} (when α is not an integer). Note that this implies that neither log z nor z^{α} have Laurent expansions about the origin: for any annulus a < |z| < b would have to be crossed by a branch cut, so the function would not be analytic in the annulus.

Multiple Branch Cuts

When there is more than one branch point, we may need more than one branch cut. For f(z) = {z(z − 1)} /_{3}^{1} there are two branch points, at 0 and 1.

We need two branch cuts; a possibility is shown in the diagram. Then no curve can wrap round either 0 or 1. For any z, we write z = r_{0}e^{iθ_{0}} where −π < θ_{0} ≤ π and z − 1 = r_{1}e^{iθ_{1}} where 0 ≤ θ_{1} < 2π. Then we define

\{z(z − 1)\} /_{3}^{1} = \sqrt[3]{r_{0}r_{1}} e^{i(θ_{0}+θ_{1})/3}.  (4.26)

This is continuous so long as we don’t cross either cut:

The value of f(z) just above the cut on the positive real axis at z = x is \sqrt[3]{x(x − 1)} (since θ_{0} = θ_{1} = 0 there); just below it is \sqrt[3]{x(x − 1)} e^{2πi/3} (θ_{0} = 0, θ_{1} = 2π). For the cut on the negative real axis we have \sqrt[3]{|x|(|x| + 1)} e^{2πi/3} just above and \sqrt[3]{|x|(|x| + 1)} just below.

Worked Example: Branch Cuts for Multiple Branch Points

What branch cuts would we require for the function

f(z) = log \frac{z − 1}{z + 1} ?

It is clear that there are branch points at ±1, but we have a non-trivial choice of branch cuts. Define z − 1 = r_{1}e^{iθ_{1}} and z + 1 = r_{2}e^{iθ_{2}}, as shown in the following diagram.

The most straightforward choice is to take two branch cuts, one emanating from each branch point to infinity. In the case shown, we choose 0 ≤ θ_{1} < 2π and −π < θ_{2} ≤ π, and the consequent single-valued definition of f(z) is

f(z) = log(z − 1) − log(z + 1)
= (log r_{1} + iθ_{1}) − (log r_{2} + iθ_{2})
= log(r_{1}/r_{2}) + i(θ_{1} − θ_{2}).

The two cuts make it impossible for z to “wind around” either of the two branch points, so we have obtained a single-valued function which is analytic except along the branch cuts.

The second possible choice is to take only one branch cut, between −1 and 1, as shown. This time, we choose both 0 ≤ θ_{1} < 2π and 0 ≤ θ_{2} < 2π (note that this seems at odds with the location of the branch cut, but this is not a problem as we will explain). The definition of f(z) is as before, but with these different ranges for θ_{1} and θ_{2}.

If z were to cross the branch cut, from above to below say, then θ_{1} would be unchanged (at π) but θ_{2} would “jump” from 0 to 2π. This is, of course, not allowed, as we may not cross branch cuts. So z cannot wind round just one of the branch points.

But it is now possible for z to wind around both of the branch points together. Consider a curve C which does so. Starting from the point of C on the positive real axis (where θ_{1} = θ_{2} = 0) and moving anti-clockwise, both θ_{1} and θ_{2} increase. When we have made one complete revolution and returned to the positive real axis, having encircled both branch points exactly once, θ_{1} and θ_{2} both suddenly “jump” from 2π back to 0. But this jump does not result in a jump in the value of θ_{1} − θ_{2}; so f(z) is not affected, and is indeed single-valued as claimed.

Exactly the same choice of branch cuts occurs for the function

g(z) = (z^{2} − 1)^{1/2}.

With the appropriate definitions of θ_{1} and θ_{2}, as above, the single-valued choice is

g(z) = (z − 1)^{1/2}(z + 1)^{1/2} = \sqrt{r_{1}r_{2}} e^{i(θ_{1}+θ_{2})/2}.

This time the single branch cut works because, when both θ_{1} and θ_{2} jump by 2π, \frac{1}{2}(θ_{1}+θ_{2}) jumps by 2π also; and e^{2πi} = 1. The cut prevents either θ_{1} or θ_{2} jumping on its own.

This idea can be extended to higher numbers of branch points in the right circumstances.

Example: consider

Φ(x, y) = Im (\frac{2T_{0}}{π} log \frac{z + a}{z − a}) + 2T_{0}  (4.27)

where z = x + iy, a is a real constant and we use the branch cuts shown below to make log analytic.

Then

Φ = \frac{2T_{0}}{π} (θ_{0} − θ_{1}) + 2T_{0}  (4.28)

where −π < θ_{0} ≤ π and 0 ≤ θ_{1} < 2π. On the circle |z| = a, what is Φ? Above the real axis, we know that θ_{0} + (π − θ_{1}) = π/2 (property of circles), so Φ = +T_{0}. Below the real axis, (−θ_{0}) + (θ_{1} − π) = π/2 (same property), so Φ = −T_{0}. We also note that everywhere in the circle |z| < a, ∇^{2}Φ = 0 as Φ is the imaginary part of an analytic function.

Hence Φ is the steady-state temperature distribution in a cylinder heated on one side to +T_{0} and on the other to −T_{0}; we solved this problem using separation of variables in Chapter 2. To see the connection between the solutions, write θ_{0} = tan^{−1}(y/(x + a)) and θ_{1} = π + tan^{−1}(y/(x − a)); then construct the Fourier Series for

\frac{2T_{0}}{π} (tan^{−1} \frac{y}{x + a} − tan^{−1} \frac{y}{x − a}) .  (4.29)

Contour Integration and Transform Theory

5.1 Path Integrals

For an integral ∫_{a}^{b} f(x) dx on the real line, there is only one way of getting from a to b. For an integral ∫ f(z) dz between two complex points a and b we need to specify which path or contour C we will use. As an example, consider

I_{1} = ∫_{C_{1}} \frac{dz}{z}  and I_{2} = ∫_{C_{2}} \frac{dz}{z}  (5.1)

where in both cases we integrate from z = −1 to z = +1 round a unit semicircle: C_{1} above, C_{2} below the real axis. Substitute z = e^{iθ}, dz = ie^{iθ} dθ:

I_{1} = ∫_{π}^{0} \frac{ie^{iθ} dθ}{e^{iθ}} = −iπ  (5.2)

but

I_{2} = ∫_{π}^{2π} i dθ = +iπ.  (5.3)

The result of a contour integration may depend on the contour.

To formally define the integral, divide C into small intervals, separated at points z_{k} (k = 0, …, N) on C, where z_{0} = a and

z_{N} = b. Let δz_{k} = z_{k+1} − z_{k} and let Δ = max_{k=0,...,N−1} |δz_{k}|. Then we

define

∫_{C} f(z) dz = lim_{Δ→0} ∑_{n=0}^{N−1} f(z_{k}) δz_{k}  (5.4)

where, as Δ → 0, N → ∞. Note that if C lies along the real axis then this definition is exactly the normal definition of a real integral.

Elementary properties

If C_{1} is a contour from w_{1} to w_{2} in ℂ, and C_{2} a contour from w_{2} to w_{3}, and C is the combined contour from w_{1} to w_{3} following first C_{1} then C_{2}, we have that ∫_{C} f(z) dz = ∫_{C_{1}} f(z) dz + ∫_{C_{2}} f(z) dz. (Obvious from definition; compare with the equivalent result on the real line, ∫_{a}^{c} f(x) dx = ∫_{a}^{b} f(x) dx + ∫_{b}^{c} f(x) dx.)

If C^{+} is a contour from w_{1} to w_{2}, and C^{−} is exactly the same contour traversed backwards, then clearly ∫_{C^{+}} f(z) dz = − ∫_{C^{−}} f(z) dz. (Cf. ∫_{a}^{b} f(x) dx = − ∫_{b}^{a} f(x) dx.)

Integration by substitution and by parts work in ℂ also.

If C has length L, then

|∫_{C} f(z) dz| ≤ L max_{C} |f(z)|  (5.5)

because at each point on C, |f(z)| ≤ max_{C} |f(z)|.

Closed contours

If C is a closed curve, then it doesn’t matter where we start from on C: ∮_{C} f(z) dz means the same thing in any case. (Notation: ∮ denotes an integral round a closed curve.) Note that if we traverse C in a negative sense (clockwise) we get negative the result we would have obtained had we traversed it in a positive sense (anticlockwise).

5.2 Cauchy’s Theorem

A simply-connected domain is a region R of the complex plane without any holes; formally, it is a region in which any closed curve encircles only points which are also in R. By a simple closed curve we mean one which is continuous, of finite length and does not intersect itself.

Cauchy’s Theorem states simply that if f(z) is analytic in a simply-connected domain R, then for any simple closed curve C in R,

∮_{C} f(z) dz = 0.  (5.6)

The proof is simple and follows from the Cauchy–Riemann equations and the Divergence Theorem in 2D:

∮_{C} f(z) dz = ∮_{C}(u + iv)(dx + i dy)

= ∮_{C}(u dx − v dy) + i ∮_{C}(v dx + u dy)

= ∫∫_{S} (−\frac{∂v}{∂x} − \frac{∂u}{∂y}) dxdy + i ∫∫_{S} (\frac{∂u}{∂x} − \frac{∂v}{∂y} ) dxdy,

by applying the Divergence Theorem, where S is the region enclosed by C. But the Cauchy–Riemann equations show that both brackets vanish, since f is analytic throughout S. The result follows.

This result is, of course, not true if C encircles a singularity (we could not then use the Cauchy–Riemann equations throughout S).

Changing the Contour

Suppose that C_{1} and C_{2} are two contours from a to b and that there are no singularities of f on or between the contours. Let C be the contour consisting of C_{1} followed by the reverse of C_{2}. C is a simple closed contour, so

∮_{C} f(z) dz = 0  (5.7)

(no singularities are enclosed). Hence

∫_{C_{1}} f(z) dz − ∫_{C_{2}} f(z) dz = 0,  (5.8)

i.e.

∫_{C_{1}} f(z) dz = ∫_{C_{2}} f(z) dz.  (5.9)

So if we have one contour, we can move it around so long as we don’t cross any singularities as we move it.

If f has no singularities anywhere, then ∫_{a}^{b} f(z) dz does not depend at all on the path chosen.

The same idea of “moving the contour” applies to closed contours; if C_{1} and C_{2} are closed contours as shown, then

∮_{C_{1}} f(z) dz = ∮_{C_{2}} f(z) dz  (5.10)

so long as there are no singularities between C_{1} and C_{2}. We prove this by considering the closed contour C shown: clearly

0 = ∮_{C} f(z) dz = ∮_{C_{1}} f(z) dz − ∮_{C_{2}} f(z) dz  (5.11)

(the two integrals along the “joins” shown cancel).

5.3 The Integral of f′(z)

For a real function f(x), ∫_{a}^{b} f′(x) dx = f(b) − f(a). This result extends immediately to complex functions, so long as both f and f′ are analytic in some simply-connected region R and the integration contour C lies entirely in R. Then

∫_{a}^{b} f′(z) dz = f(b) − f(a)  (5.12)

for any complex points a, b in R.

Note that the specified conditions ensure that the integral on the LHS is independent of exactly which path in R is used from a to b, using the results of §5.2.

Examples:

(i) ∫_{0}^{i} z dz = \frac{1}{2}(i^{2} − 0^{2}) = −\frac{1}{2}. (f and f′ are analytic in the whole of ℂ, so the LHS is path-independent.)

(ii) ∫_{C} e^{z} dz, where C is the semicircular contour joining −1 to +1 along |z| = 1 above the real axis, is equal to e − e^{−1}.

(iii) ∫_{1+i}^{−1+i} z^{−1} dz via a straight contour. Note that z^{−1} is not analytic everywhere, so we do need to specify the contour; but we can define a simply-connected region R, given by Im z > \frac{1}{2} say, in which it is analytic, and C lies entirely in R. Let f(z) = log z with the standard branch cut, so that f(z) is also analytic in R; then

∫_{1+i}^{−1+i} z^{−1} dz = log(−1 + i) − log(1 + i)
= log \sqrt{2} + \frac{3}{4}πi − (log \sqrt{2} + \frac{1}{4}πi)
= \frac{1}{2}πi.

(iv) Now consider ∫_{1+i}^{−1+i} z^{−1} dz via the contour shown. Define R as in the diagram; we cannot now choose the standard branch cut for log z (since C would cross it), so we choose a cut along the positive imaginary axis, and define log re^{iθ} = log r + iθ where −\frac{3π}{2} < θ ≤ \frac{π}{2} . Then

∫_{C} z^{−1} dz = log(−1 + i) − log(1 + i)
= log \sqrt{2} + (−\frac{5}{4}π)i − (log \sqrt{2} + \frac{1}{4}πi)
= −\frac{3}{2}πi.

5.4 The Calculus of Residues

The Contour Integral of a Laurent Expansion

Consider a single term a_{n}(z−z_{0})^{n} of an expansion, integrated round a closed curve C which encircles z_{0} in a positive sense (i.e., anticlockwise) once. For n ≥ 0, we can use Cauchy’s Theorem to obtain immediately

∮_{C} a_{n}(z − z_{0})^{n} dz = 0.  (5.13)

For n < 0, first change the contour C to C_{ε}, a circle of radius ε about z_{0}, using the ideas of §5.2. On C_{ε}, z = z_{0} + εe^{iθ} and so

∮_{C} a_{n}(z − z_{0})^{n} dz = ∫_{0}^{2π} a_{n}ε^{n}e^{inθ} iεe^{iθ} dθ

= ia_{n}ε^{n+1} ∫_{0}^{2π} e^{i(n+1)θ} dθ

= \begin{cases}ia_{n}ε^{n+1} [\frac{e^{i(n+1)θ}}{i(n + 1)}]_{0}^{2π} & n ≠ −1 \\ ia_{n}ε^{n+1}(2π) & n = −1 \end{cases}

= \begin{cases}0 & n ≠ −1 \\ 2πia_{−1} & n = −1 \end{cases}

We deduce that for a function f(z) with a singularity at z_{0}, and a contour C encircling the singularity in a positive sense,

∮_{C} f(z) dz = ∑_{n=−∞}^{∞} ∮_{C} a_{n}(z − z_{0})^{n} dz = 2πia_{−1} = 2πi res_{z=z_{0}} f(z).  (5.14)

We can also obtain the result as follows, using the method of §5.3:

∮_{C} a_{n}(z − z_{0})^{n} dz = \begin{cases}\frac{a_{n}}{n + 1}[(z − z_{0})^{n+1}]_{C} & n ≠ −1 \\ a_{n}[log(z − z_{0})]_{C} & n = −1 \end{cases}

= \begin{cases}0 & n ≠ −1 (because (z − z_{0})^{n+1} is single-valued) \\ 2πia_{−1} & n = −1 (because θ changes by 2π) \end{cases}

The Residue Theorem

Suppose that f(z) is analytic in a simply-connected region R except for a finite number of poles at z_{1}, z_{2}, …, z_{n}; and that a simple closed curve C encircles the poles anticlockwise. Then

∮_{C} f(z) dz = 2πi ∑_{k=1}^{n} res_{z=z_{k}} f(z).  (5.15)

(We have just proved this in the case of a single pole.)

Proof: Consider the curve Ĉ shown. Ĉ encircles no poles, so

∮Ĉ f(z) dz = 0  (5.16)

by Cauchy’s Theorem. But we can also work out the integral round Ĉ by adding together several contributions: the large outer curve (which is the same as C), the small circles around each pole, and the contributions from the lines joining the outer curve to the inner circles. For each k, the contribution from the small circle round z_{k} is −2πi res_{z=z_{k}} f(z) because the small circle goes clockwise round z_{k}. Also, the contribution from the line joining the outer curve to the small circle cancels exactly with the contribution from the line going back. Hence

0 = ∮Ĉ f(z) dz = ∮_{C} f(z) dz + ∑_{k=1}^{n}(−2πi res_{z=z_{k}} f(z))  (5.17)

from which the result follows.

5.5 Cauchy’s Formula for f(z)

Suppose that f(z) is analytic in a region R and that z_{0} lies in R. Then Cauchy’s formula states that

f(z_{0}) = \frac{1}{2πi} ∮_{C} \frac{f(z)}{z − z_{0}} dz  (5.18)

where C is any closed contour in R encircling z_{0} once anticlockwise.

Proof: f(z)/(z − z_{0}) is analytic except for a simple pole at z_{0}, where it has residue f(z_{0}). Using the Residue Theorem,

∮_{C} \frac{f(z)}{z − z_{0}} dz = 2πif(z_{0})  (5.19)

as required.

Note: Cauchy’s formula says that if we know f on C then we know it at all points within C. We can see that this must be so by the uniqueness theorem of Chapter 2: u and v, the real and imaginary parts of f, are harmonic, so if they are specified on C (Dirichlet boundary conditions), then there is a unique solution for u and v inside C.

Exercise: show that if instead f is analytic except for a singularity at z_{0}, and has a Laurent expansion ∑_{m=−∞}^{∞} a_{m}(z − z_{0})^{m}, then the coefficients of the expansion are given by

a_{n} = \frac{1}{2πi} ∮_{C} \frac{f(z)}{(z − z_{0})^{n+1}} dz.  (5.20)

If we differentiate Cauchy’s formula with respect to z_{0} (differentiating under the ∮ sign on the RHS), we see that

f′(z_{0}) = \frac{1}{2πi} ∮_{C} \frac{f(z)}{(z − z_{0})^{2}} dz.  (5.21)

So f′(z_{0}) is known for all z_{0} inside C. Continuing this process,

f^{(n)}(z_{0}) = \frac{n!}{2πi} ∮_{C} \frac{f(z)}{(z − z_{0})^{n+1}} dz,  (5.22)

and f^{(n)}(z_{0}) is known. So at any point where f is analytic, i.e. differentiable once, all its derivatives exist; hence it is differentiable infinitely many times.

5.6 Applications of the Residue Calculus

Suppose we wish to evaluate

I = ∫_{0}^{∞} \frac{dx}{1 + x^{2}}  (5.23)

(which we can already do using trigonometric substitutions).

Consider

∮_{C} \frac{dz}{1 + z^{2}}  (5.24)

where C is the contour shown: from −R to R along the real axis (C_{0}) then returning to −R via a semicircle of radius R in the upper half-plane (C_{R}). Now (1 + z^{2})^{−1} = (z + i)^{−1}(z − i)^{−1}, so the only singularity enclosed by C is a simple pole at z = i, where the residue is lim_{z→i}(z + i)^{−1} = 1/2i. Hence

∫_{C_{0}} \frac{dz}{1 + z^{2}} + ∫_{C_{R}} \frac{dz}{1 + z^{2}} = ∮_{C} \frac{dz}{1 + z^{2}} = 2πi \frac{1}{2i} = π.  (5.25)

Now

∫_{C_{0}} \frac{dz}{1 + z^{2}} = ∫_{−R}^{R} \frac{dx}{1 + x^{2}} → 2I as R → ∞.  (5.26)

Consider ∫_{C_{R}} dz/(1 + z^{2}): the integrand (1 + z^{2})^{−1} is of order R^{−2} on the semicircle, but the length of the contour is πR. Hence

|∫_{C_{R}} \frac{dz}{1 + z^{2}} | ≤ πR × O(R^{−2}) = O(R^{−1}) → 0 as R → ∞.  (5.27)

Combining all these results and taking the limit as R → ∞,

2I + 0 = π,  (5.28)

i.e. I = π/2.

This example is not in itself impressive. But the power of the method is clear when we see how easily it adapts to other such integrals (for which it would not be easy, or would be impossible, to use substitutions). Examples:

(i) We wish to calculate

I = ∫_{0}^{∞} \frac{dx}{(x^{2} + a^{2})^{2}}  (5.29)

where a > 0 is a real constant. We consider ∮_{C} dz/(z^{2} + a^{2})^{2}; most of the above analysis is unchanged. The poles now occur at z = ±ia, and they both have order 2; only the pole at +ia is enclosed by C. The residue there is

lim_{z→ia} \frac{d}{dz} \frac{1}{(z + ia)^{2}} = lim_{z→ia} \frac{−2}{(z + ia)^{3}} = \frac{−2}{−8ia^{3}} = −\frac{1}{4}ia^{−3}.  (5.30)

The integral round the semicircle still vanishes as R → ∞, since now

|∫_{C_{R}} \frac{dz}{(z^{2} + a^{2})^{2}} | ≤ πR × O(R^{−4}) = O(R^{−3}).  (5.31)

Therefore

2I = 2πi(−\frac{1}{4}ia^{−3}) = π/2a^{3},  (5.32)

i.e., I = π/4a^{3}.

(ii) For I = ∫_{0}^{∞} dx/(1 + x^{4}), the (simple) poles are at e^{πi/4}, e^{3πi/4}, e^{−πi/4} and e^{−3πi/4}.

Only the first two poles are enclosed. The residue at e^{πi/4} is

lim_{z→e^{πi/4}} \frac{z − e^{πi/4}}{1 + z^{4}} = lim_{z→e^{πi/4}} \frac{1}{4z^{3}} = \frac{1}{4}e^{−3πi/4} = −\frac{1}{4}e^{πi/4}  (5.33)

using L’Hôpital’s Rule, and at e^{3πi/4} it is (similarly) \frac{1}{4}e^{−πi/4}. Hence

2I = 2πi(−\frac{1}{4}e^{πi/4} + \frac{1}{4}e^{−πi/4}) = 2πi(−\frac{1}{4})(2i sin \frac{π}{4}) = π sin \frac{π}{4},  (5.34)

i.e., I = π/2\sqrt{2}.

(iii) For I = ∫_{0}^{∞} x^{2}dx/(1 + x^{4}), the poles are as in (ii) but with residues +\frac{1}{4}e^{−πi/4} and −\frac{1}{4}e^{πi/4} respectively (check for yourself). So the value of the integral is unchanged.

(iv) For I = ∫_{0}^{∞} dx/(1 + x^{4}) again, an alternative to the method used in example (ii) above (and similarly in example (iii) above) is to use a contour which is just a quarter-circle, as shown.

Let C consist of the real axis from 0 to R (C_{0}); the arc of circle from R to iR (C_{1}); and the imaginary axis from iR to 0 (C_{2}). Now ∫_{C_{0}} dz/(1 + z^{4}) → I as R → ∞; and, along C_{2}, we substitute z = iy to obtain

∫_{C_{2}} \frac{dz}{1 + z^{4}} = ∫_{R}^{0} \frac{i dy}{1 + (iy)^{4}} = −i ∫_{0}^{R} \frac{dy}{1 + y^{4}} → −iI  as R → ∞.  (5.35)

The integral along C_{1} vanishes as R → ∞, using the same argument as for C_{R} above, but this time we only enclose one pole, which makes the calculation easier. Hence

I − iI = 2πi(−\frac{1}{4}e^{πi/4}) = −\frac{1}{2}πe^{3πi/4} ⟹ I = π/2\sqrt{2}

as before.

Worked Example: Contour Integration - Integrals of Trigonometric Functions

We wish to evaluate

I = ∫_{0}^{2π} \frac{dθ}{a + cos θ}

where a > 1 (so that the integrand is always finite). Substitute z = e^{iθ}, so that dz = iz dθ and cos θ = \frac{1}{2}(z + z^{−1}). As θ increases from 0 to 2π, z moves round the circle C of radius 1 in the complex plane. Hence

I = ∮_{C} \frac{(iz)^{−1} dz}{a + \frac{1}{2}(z + z^{−1})} = −2i ∮_{C} \frac{dz}{z^{2} + 2az + 1}.

The integrand has poles at

z_{±} = −a ± \sqrt{a^{2} − 1},

both on the real axis. Note that z_{+} is inside the unit circle (check that a−1 < \sqrt{a^{2} − 1} < a, so −1 < z_{+} < 0) whereas z_{−} is outside it. The integrand is equal to

\frac{1}{(z − z_{+})(z − z_{−})}

so the residue at z = z_{+} is 1/(z_{+} − z_{−}) = 1/2\sqrt{a^{2} − 1}. Hence

I = −2i (\frac{2πi}{2\sqrt{a^{2} − 1}}) = \frac{2π}{\sqrt{a^{2} − 1}}.

Worked Example: Contour Integration - Integration Round a Branch Cut

We wish to evaluate

I = ∫_{0}^{∞} \frac{x^{α}}{1 + \sqrt{2} x + x^{2}} dx

where −1 < α < 1 so that the integral converges. We will need a branch cut for z^{α}; we take this along the positive real axis and define

z^{α} = r^{α}e^{iαθ}

where z = re^{iθ} and 0 ≤ θ < 2π.

Consider

∮_{C} \frac{z^{α}}{1 + \sqrt{2} z + z^{2}} dz

where the keyhole contour C consists of a large circle C_{R} of radius R, a small circle C_{ϵ} of radius ϵ (to avoid the singularity of z^{α} at z = 0) and two lines just above and below the branch cut, as shown.

The contribution from C_{R} is O(R^{α−2}) × 2πR = O(R^{α−1}) → 0 as R → ∞.

The contribution from C_{ϵ} is (substituting z = ϵe^{iθ} on C_{ϵ})

∫_{2π}^{0} \frac{ϵ^{α}e^{iαθ}}{1 + \sqrt{2} ϵe^{iθ} + ϵ^{2}e^{2iθ}} iϵe^{iθ} dθ = O(ϵ^{α+1}) → 0

as ϵ → 0.

The contribution from just above the branch cut is

∫_{ϵ}^{R} \frac{x^{α}}{1 + \sqrt{2} x + x^{2}} dx → I

as ϵ → 0 and R → ∞. The contribution from just below the branch cut is

∫_{R}^{ϵ} \frac{x^{α}e^{2απi}}{1 + \sqrt{2} x + x^{2}} dx → −e^{2απi}I

as ϵ → 0 and R → ∞.

Hence

∮_{C} \frac{z^{α}}{1 + \sqrt{2} z + z^{2}} dz → (1 − e^{2απi})I

as ϵ → 0 and R → ∞.

But the integrand is equal to

\frac{z^{α}}{(z − e^{3πi/4})(z − e^{5πi/4})}

(by finding the roots of the quadratic), so the poles inside C are at e^{3πi/4} with residue e^{3απi/4}/(\sqrt{2} i) and at e^{5πi/4} with residue e^{5απi/4}/(−\sqrt{2} i). Hence, taking the limits ϵ → 0 and R → ∞,

(1 − e^{2απi})I = 2πi (\frac{e^{3απi/4}}{\sqrt{2} i} + \frac{e^{5απi/4}}{−\sqrt{2} i} ) ,

i.e.,

e^{απi}(e^{−απi} − e^{απi})I = \sqrt{2} πe^{απi}(e^{−απi/4} − e^{απi/4}).

We conclude that

I = \sqrt{2} π \frac{sin(απ/4)}{sin(απ)}.

Jordan’s Lemma

For many applications (in particular, ones involving Fourier transforms) we need to show that

∫_{C_{R}} f(z)e^{iλz} dz → 0  (5.36)

as R → ∞, where λ > 0 is some real constant and f is an analytic function (except possibly for a finite number of poles). Jordan’s Lemma states that this is true so long as f(z) → 0 as |z| → ∞. For λ < 0, the same conclusion holds for the semicircular contour C′_{R} in the lower half-plane.

Note that this result is obvious if f(z) = O(|z|^{−2}) as |z| → ∞ – i.e., if f(z) = O(R^{−2}) on C_{R} – by the following argument. First note that e^{iλz} = e^{iλ(x+iy)} = e^{−λy}e^{iλx}, and y ≥ 0 on C_{R}, so |e^{iλz}| = e^{−λy} ≤ 1 on C_{R}. Hence

|∫_{C_{R}} f(z)e^{iλz} dz| ≤ πR max_{C_{R}} |f(z)|
= πR × O(R^{−2}) → 0 as R → ∞.

Jordan’s Lemma simply extends the result from functions satisfying f(z) = O(|z|^{−2}) to any function satisfying f(z) → 0 as |z| → ∞. Examples:

∫_{C_{R}} \frac{e^{2iz}}{z} dz → 0 as R → ∞;  ∫_{C′_{R}} \frac{e^{−iz}}{z^{2}} dz → 0 as R → ∞.  (5.37)

The proof of Jordan’s Lemma stems from the fact that for 0 ≤ θ ≤ π/2, sin θ ≥ 2θ/π. Now

|∫_{C_{R}} f(z)e^{iλz} dz| ≤ max_{C_{R}} |f(z)| ∫_{0}^{π} |e^{iλz}| |Re^{iθ}| dθ
= R max |f(z)| ∫_{0}^{π} e^{−λR sin θ} dθ

[using y = R sin θ]

= 2R max |f(z)| ∫ /_{2}^{π} ; 0 e^{−λR sin θ} dθ

≤ 2R max |f(z)| ∫ /_{2}^{π} ; 0 e^{−2λRθ/π} dθ

= \frac{π}{λ} (1 − e^{−λR}) max |f(z)|

→ 0 as R → ∞.

A similar proof holds on C′_{R} for λ < 0.

Worked Example: Contour Integration - Inverse Fourier Transforms

Consider the real function

f(x) = \begin{cases}0 & x < 0 \\ e^{−ax} & x > 0 \end{cases}

where a > 0 is a real constant. The Fourier Transform of f(x) is

f(k) = ∫_{−∞}^{∞} f(x)e^{−ikx} dx
= ∫_{0}^{∞} e^{−ax−ikx} dx
= −\frac{1}{a + ik} [e^{−ax−ikx}]_{0}^{∞}
= \frac{1}{a + ik} .  (5.38)

We shall verify the Inverse Fourier Transform by evaluating

\frac{1}{2π} ∫_{−∞}^{∞} f(k)e^{ikx} dk.

In the complex k-plane, let C_{0} be the contour from −R to R on the real axis, C_{R} be the semicircle of radius R in the upper half plane and C′_{R} be the semicircle of radius R in the lower half plane. Let C be C_{0} followed by C_{R} (this is known as closing in the upper half plane), and let C′ be C_{0} followed by C′_{R} (closing in the lower half plane).

Now f(k) has only one pole, at k = ia, which is simple, so

∮_{C} f(k)e^{ikx} dk = 2πi res_{k=ia} \frac{e^{ikx}}{i(k − ia)} = 2πe^{−ax},

whereas

∮_{C′} f(k)e^{ikx} dk = 0.

(Note that C′ is traversed in a negative sense, so if there had been any poles within C′ we would have had to introduce a minus sign.)

Now, if x > 0, we can apply Jordan’s Lemma (with λ = x) to C_{R} to show that ∫_{C_{R}} f(k)e^{ikx} dk → 0 as R → ∞, since f(k) = O(1/k) as |k| → ∞. Hence for x > 0,

\frac{1}{2π} ∫_{−∞}^{∞} f(k)e^{ikx} dk = \frac{1}{2π} lim_{R→∞} ∫_{C_{0}} f(k)e^{ikx} dk
= \frac{1}{2π} lim_{R→∞} (∮_{C} f(k)e^{ikx} dk −∫_{C_{R}} f(k)e^{ikx} dk)
= e^{−ax}.  (5.39)

For x < 0 we close in the lower half plane instead, and the same analysis applies to C′:

\frac{1}{2π} ∫_{−∞}^{∞} f(k)e^{ikx} dk = \frac{1}{2π} lim_{R→∞} (∮_{C′} f(k)e^{ikx} dk −∫_{C′_{R}} f(k)e^{ikx} dk)
= 0.  (5.40)

Hence, combining the above results, we obtain

\frac{1}{2π} ∫_{−∞}^{∞} f(k)e^{ikx} dk = \begin{cases}0 & x < 0 \\ e^{−ax} & x > 0 \end{cases}

as expected.

Note that by taking real and imaginary parts of this equality we can deduce the values of particular real integrals. The imaginary part gives

\frac{1}{2π} ∫_{−∞}^{∞} \frac{a sin kx − k cos kx}{a^{2} + k^{2}} dk = 0,

which is obvious anyway as the integrand is an odd function of k. But the real part gives

\frac{1}{2π} ∫_{−∞}^{∞} \frac{a cos kx + k sin kx}{a^{2} + k^{2}} dk = \begin{cases}0 & x < 0 \\ e^{−ax} & x > 0 \end{cases}

and in particular

∫_{−∞}^{∞} \frac{a cos θ + θ sin θ}{a^{2} + θ^{2}} dθ = 2πe^{−a}.

Worked Example: Contour Integration - Singular Point on the Real Axis

We wish to evaluate

∫_{−∞}^{∞} \frac{sin x}{x} dx.

This integrand is well-behaved at the origin, so the integral is non-singular. But the obvious approach via contour integration using

\frac{1}{2i} ∫_{−∞}^{∞} \frac{e^{iz} − e^{−iz}}{z} dz

runs into trouble because we cannot apply Jordan’s Lemma to the integrand. To get round this we might try to split it into two separate integrals, to each of which Jordan’s Lemma does apply, but then we find that our contour passes through a pole of the integrand.

Instead, we write

∫_{−∞}^{∞} \frac{sin x}{x} dx = lim_{{ϵ→0 \atop R→∞}}(∫_{−R}^{−ϵ} \frac{sin x}{x} dx + ∫_{ϵ}^{R} \frac{sin x}{x} dx)

= Im lim_{{ϵ→0 \atop R→∞}}(∫_{−R}^{−ϵ} \frac{e^{iz}}{z} dz + ∫_{ϵ}^{R} \frac{e^{iz}}{z} dz).  (5.41)

Let C be the contour from −R to −ϵ, then round a semicircle C_{ϵ} of radius ϵ, then from ϵ to R, and returning via a semicircle C_{R} of radius R. Then C encloses no poles of e^{iz}/z, so

∫_{−R}^{−ϵ} \frac{e^{iz}}{z} dz + ∫_{ϵ}^{R} \frac{e^{iz}}{z} dz = − ∫_{C_{ϵ}} \frac{e^{iz}}{z} dz − ∫_{C_{R}} \frac{e^{iz}}{z} dz.

Jordan’s Lemma tells us that the integral round C_{R} vanishes as R → ∞. On C_{ϵ}, z = ϵe^{iθ} and e^{iz} = 1 + O(ϵ); so

∫_{C_{ϵ}} \frac{e^{iz}}{z} dz = ∫_{π}^{0} \frac{1 + O(ϵ)}{ϵe^{iθ}} iϵe^{iθ} dθ = −iπ + O(ϵ).

Hence, taking the limit as ϵ → 0 and R → ∞,

∫_{−∞}^{∞} \frac{sin x}{x} dx = Im(iπ) = π.

A similar method works for

∫_{−∞}^{∞} \frac{sin^{2} x}{x^{2}} dx;

write sin^{2} x = \frac{1}{2} Re(1 − e^{2ix}), and then

∫_{−∞}^{∞} \frac{sin^{2} x}{x^{2}} dx = \frac{1}{2} Re lim_{{ϵ→0 \atop R→∞}}(∫_{−R}^{−ϵ} \frac{1 − e^{2iz}}{z^{2}} dz + ∫_{ϵ}^{R} \frac{1 − e^{2iz}}{z^{2}} dz)

= \frac{1}{2} Re lim_{{ϵ→0 \atop R→∞}}(− ∫_{C_{ϵ}} \frac{1 − e^{2iz}}{z^{2}} dz − ∫_{C_{R}} \frac{1 − e^{2iz}}{z^{2}} dz). (5.42)

The integral round C_{R} can be shown to vanish as R → ∞ by standard techniques (Jordan’s Lemma is not, however, applicable), and the integral round C_{ϵ} can be evaluated as before (expanding e^{2iz} to slightly higher order in ϵ than before), giving

∫_{C_{ϵ}} \frac{1 − e^{2iz}}{z^{2}} dz = −2π + O(ϵ).

Hence

∫_{−∞}^{∞} \frac{sin^{2} x}{x^{2}} dx = π

as well!

An alternative approach for both examples is to note that, for instance, (sin z)/z has a removable singularity at the origin. Having removed the singularity, we have an analytic integrand, and therefore the original contour along the real axis can be moved to one which does not pass through the origin. It is now possible to write sin z = (e^{iz} −e^{−iz})/2i, split the integrand in two, and apply Jordan’s Lemma to each part separately.

Worked Example: Contour Integration - Using the Residue at Infinity

We wish to evaluate

I = ∫_{−1}^{1} \sqrt{1 − x^{2}} dx.

Consider (z^{2} − 1) /_{2}^{1} with a branch cut from −1 to 1 (as explained in a worked example in Chapter 4). For z = x on the real axis, just above the cut we have

(z^{2} − 1) /_{2}^{1} = i\sqrt{1 − x^{2}}

(θ_{1} = π and θ_{2} = 0 in the notation of the branch cut example), whereas just below

(z^{2} − 1) /_{2}^{1} = −i\sqrt{1 − x^{2}}

(θ_{1} = π, θ_{2} = 2π).

Let C be the closed contour running along the x-axis from +1 to −1 above the cut and back again to +1 below it; at each end C runs round a circle of radius ϵ to avoid the branch point. For small ϵ the integrals round the small circles are negligible (check); so

∮_{C}(z^{2} − 1) /_{2}^{1} dz = ∫_{1}^{−1} i\sqrt{1 − x^{2}} dx + ∫_{−1}^{1}(−i\sqrt{1 − x^{2}}) dx
= −2iI.  (5.43)

We can deform the contour of integration from C to C_{R}, the circle of radius R, as there are no singularities between C and C_{R}. Hence

I = \frac{1}{2}i ∮_{C_{R}}(z^{2} − 1) /_{2}^{1} dz.  (1)

Let ζ = 1/z, and let C_{1/R} be the circle of radius 1/R traversed clockwise, so that C_{1/R} is the image of C_{R} under the transformation z ↦ ζ. Then dz = −ζ^{−2} dζ, so making the substitution,

I = \frac{1}{2}i ∮_{C_{1/R}}(ζ^{−2} − 1) /_{2}^{1}(−ζ^{−2}) dζ.

Now

−ζ^{−2}(ζ^{−2} − 1) /_{2}^{1} = −ζ^{−3}(1 − ζ^{2}) /_{2}^{1} = −ζ^{−3}(1 − \frac{1}{2}ζ^{2} + ⋯),

so the integrand has a singularity at ζ = 0 with residue \frac{1}{2}. Therefore (introducing a minus sign because the contour C_{1/R} is traversed in a negative sense),

I = −\frac{1}{2}i × 2πi × \frac{1}{2} = \frac{π}{2} .

Note that we could also have argued as follows: (z^{2} − 1) /_{2}^{1} has residue −\frac{1}{2} at infinity (using the definition of the residue at infinity and the above working); since C_{R} “encircles infinity”, from (1) we immediately obtain I = \frac{1}{2}i × 2πi × (−\frac{1}{2}) = π/2 as before. Comparison of the two techniques reveals the reason for the unintuitive definition of the residue at infinity.

5.7 Laplace Transforms

The Fourier transform is a powerful technique for solving differential equations and for investigating many physical problems, but not all functions have a Fourier transform: the integral defining the transform does not converge unless the function tends to zero at infinity.

To get around this restriction, we can use another kind of transform known as the Laplace transform. The price we pay is a different restriction: it is only defined for functions which are zero for t < 0 (by convention). From now on, we shall make this assumption, so that if we refer to the function f(t) = e^{t} for instance,

we really mean the function f(t) = \begin{cases}0 & t < 0, \\ e^{t} & t ≥ 0. \end{cases}

The Laplace transform of a function f(t) is defined by

f̄(p) = ∫_{0}^{∞} f(t)e^{−pt} dt  (5.44)

where p may be complex. The notation ℒ [f] or ℒ [f(t)] is also used for f̄(p); and the symbol s is often used instead of p. Many functions – for instance, t and e^{t} – which do not have Fourier transforms do have Laplace transforms; however, there are still exceptions (e.g., e^{t^{2}}). Laplace transforms are particularly useful in initial value problems, where we are given the state of a system at t = 0 and desire to find its state for t > 0.

Examples:

(i) ℒ [1] = ∫_{0}^{∞} e^{−pt} dt = \frac{1}{p}.

(ii) ℒ [t] = ∫_{0}^{∞} te^{−pt} dt = [−\frac{1}{p}te^{−pt}]_{0}^{∞} + \frac{1}{p} ∫_{0}^{∞} e^{−pt} dt = \frac{1}{p^{2}} .

(iii) ℒ [e^{λt}] = ∫_{0}^{∞} e^{(λ−p)t} dt = \frac{1}{p − λ}.

(iv) ℒ [sin t] = ℒ [\frac{1}{2i}(e^{it} − e^{−it})] = \frac{1}{2i} (\frac{1}{p − i} − \frac{1}{p + i}) = \frac{1}{p^{2} + 1}.

Note that, strictly speaking, in example (iii), the integral only converges for Re p > Re λ (otherwise the integrand, e^{(λ−p)t}, diverges as t → ∞). However, once we have calculated the integral for Re p > Re λ we can consider f̄(p) to exist everywhere in the complex p-plane (except for singularities such as at p = λ in this example). This process of extending a complex function which is initially only defined in some part of the complex plane to the whole of the plane is known as analytic continuation.

It is useful to have a “library” of Laplace transforms to hand; some common ones are listed below.

f(t)  f̄(p)  f(t)  f̄(p)
1  \frac{1}{p}  t^{n}  \frac{n!}{p^{n+1}}

e^{λt}  \frac{1}{p − λ}  t^{n}e^{λt}  \frac{n!}{(p − λ)^{n+1}}

sin ωt  \frac{ω}{p^{2} + ω^{2}}  cos ωt  \frac{p}{p^{2} + ω^{2}}

sinh λt  \frac{λ}{p^{2} − λ^{2}}  cosh λt  \frac{p}{p^{2} − λ^{2}}

e^{λt} sin ωt  \frac{ω}{(p − λ)^{2} + ω^{2}}  e^{λt} cos ωt  \frac{p − λ}{(p − λ)^{2} + ω^{2}}

δ(t)  1  δ(t − t_{0})  e^{−pt_{0}}

Elementary Properties of the Laplace Transform

(i) Linearity: ℒ [αf(t) + βg(t)] = αf̄(p) + βḡ(p).

(ii) Change of scale: using the substitution t′ = λt,

ℒ [f(λt)] = ∫_{0}^{∞} f(λt)e^{−pt} dt = \frac{1}{λ} ∫_{0}^{∞} f(t′)e^{−(p/λ)t′} dt′ = \frac{1}{λ}f̄(\frac{p}{λ}) .  (5.45)

(iii) Shifting theorem: ℒ [e^{λt}f(t)] = f̄(p − λ). (Easy to check.)

(iv) Derivative of a Laplace transform:

ℒ [tf(t)] = −\frac{d}{dp}f̄(p).  (5.46)

Proof:

f̄(p) = ∫_{0}^{∞} f(t)e^{−pt} dt ⟹ \frac{d}{dp}f̄(p) = − ∫_{0}^{∞} tf(t)e^{−pt} dt.

By repeating this trick n times, we see that the Laplace transform of t^{n}f(t) is (−1)^{n}f̄^{(n)}(p). Examples:

ℒ [t sin t] = −\frac{d}{dp} \frac{1}{p^{2} + 1} = \frac{2p}{(p^{2} + 1)^{2}} ;  ℒ [t^{n}] = (−1)^{n} \frac{d^{n}}{dp^{n}} \frac{1}{p} = \frac{n!}{p^{n+1}} .  (5.47)

(v) Laplace transform of a derivative:

ℒ [\frac{df}{dt} ] = pf̄(p) − f(0).  (5.48)

Proof:

∫_{0}^{∞} \frac{df}{dt} e^{−pt} dt = [f(t)e^{−pt}]_{0}^{∞} + p ∫_{0}^{∞} f(t)e^{−pt} dt = pf̄(p) − f(0).

We can deduce that

ℒ [\frac{d^{2}f}{dt^{2}} ] = pℒ [\frac{df}{dt} ] − ḟ(0) = p^{2}f̄(p) − pf(0) − ḟ(0)

and so on.

(vi) Asymptotic limits: pf̄(p) → f(0) as p → ∞, and pf̄(p) → lim_{t→∞} f(t) as p → 0. Proofs: from (v) above,

pf̄(p) = f(0) + ∫_{0}^{∞} \frac{df}{dt} e^{−pt} dt,

so as p → ∞ (and therefore e^{−pt} → 0 for all t > 0), pf̄(p) → f(0). Similarly, as p → 0, e^{−pt} → 1 so that

pf̄(p) → f(0) + ∫_{0}^{∞} \frac{df}{dt} dt = f(0) + [f(t)]_{0}^{∞} = lim_{t→∞} f(t).  (5.49)

Solving Differential Equations using Laplace Transforms

The Laplace transform is particularly suited to the solution of initial value problems. Example: solve

ÿ + 5ẏ + 6y = 0  (5.50)

for y(t) subject to y(0) = 1, ẏ(0) = −4. Taking Laplace transforms, and using the results for the Laplace transform of a derivative, we see that

(p^{2}ȳ(p) − p + 4) + 5(pȳ(p) − 1) + 6ȳ(p) = 0,  (5.51)

which we may solve for ȳ(p):

ȳ(p) = \frac{p + 1}{p^{2} + 5p + 6} = \frac{p + 1}{(p + 2)(p + 3)} = \frac{2}{p + 3} − \frac{1}{p + 2}  (5.52)

using partial fractions. We now need to invert ȳ(p) to find y(t); in general we must use the inversion formula described below, but in many cases (such as this one) it is possible to “spot” the answer using the “library” of transforms given above (and taking advantage of the fact that inverse Laplace transforms are unique). Here, we know that ℒ [e^{λt}] = 1/(p − λ); hence

y(t) = 2e^{−3t} − e^{−2t}.  (5.53)

The Convolution Theorem for Laplace Transforms

The convolution of two functions f(t) and g(t) is

(f ∗ g)(t) = ∫_{−∞}^{∞} f(t − t′)g(t′) dt′.  (5.54)

We are dealing here with functions which vanish for t < 0, so this reduces to

(f ∗ g)(t) = ∫_{0}^{t} f(t − t′)g(t′) dt′  (5.55)

since g(t′) = 0 for t′ < 0 and f(t−t′) = 0 for t′ > t. The convolution theorem for Laplace transforms then states that

ℒ [f ∗ g] = f̄(p)ḡ(p).  (5.56)

Proof:

ℒ [f ∗ g] = ∫_{0}^{∞} {∫_{0}^{t} f(t − t′)g(t′) dt′} e^{−pt} dt
= ∫_{0}^{∞} {∫_{0}^{t} f(t − t′)g(t′)e^{−pt} dt′} dt.

From the diagram, we see that we can change the order of integration in the (t, t′)- plane, giving

ℒ [f ∗ g] = ∫_{0}^{∞} {∫_{t′}^{∞} f(t − t′)g(t′)e^{−pt} dt} dt′
= ∫_{0}^{∞} {∫_{0}^{∞} f(t′′)e^{−pt′′}e^{−pt′} dt′′} g(t′) dt′
[substituting t′′ = t − t′]
= ∫_{0}^{∞}{f̄(p)e^{−pt′}}g(t′) dt′
= f̄(p) ∫_{0}^{∞} g(t′)e^{−pt′} dt′
= f̄(p)ḡ(p)

as required.

The Inverse Laplace Transform

Inverting Laplace transforms is more difficult than inverting Fourier transforms because it is always necessary to perform a contour integration. Given f̄(p), we can calculate f(t) using the Bromwich inversion formula

f(t) = \frac{1}{2πi} ∫_{γ−i∞}^{γ+i∞} f̄(p)e^{pt} dp.  (5.57)

Here γ is a real constant, and the Bromwich inversion contour Γ runs from γ − i∞ to γ + i∞ along a straight line. Γ must lie to the right of all the singularities of f̄(p).

Note that it is possible to derive the Bromwich inversion formula from the inverse Fourier transform by substituting p = ik and noting that f̄(p) = f(−ip) where f(k) is the Fourier transform of f(t). The only difference is in the detail of the inversion contour.

Suppose that f̄(p) has only poles, and no other singularities; all these poles lie to the left of Γ. When t < 0, consider the integral round the contour C shown consisting of C_{0} followed by C′_{R}. C encloses no poles, so

∮_{C} f̄(p)e^{pt} dp = 0.  (5.58)

Now on C′_{R}, Re p ≥ γ, so Re(pt) ≤ γt (since t < 0) and hence |e^{pt}| ≤ e^{γt}. Therefore if f̄(p) = O(|p|^{−2}) as |p| → ∞ – i.e., if f̄(p) = O(R^{−2}) on C′_{R} – then

|∫_{C′_{R}} f̄(p)e^{pt} dp| ≤ πRe^{γt} × O(R^{−2}) → 0  (5.59)

as R → ∞. In fact the same is true even if we only have f̄(p) → 0 as |p| → ∞, by a slight modification of Jordan’s Lemma. So in either case,

∫_{Γ} f̄(p)e^{pt} dp = lim_{R→∞} ∫_{C_{0}} f̄(p)e^{pt} dp

= lim_{R→∞}(∮_{C} f̄(p)e^{pt} dp − ∫_{C′_{R}} f̄(p)e^{pt} dp)

= 0 − 0 = 0,

and therefore for t < 0 the inversion formula gives

f(t) = \frac{1}{2πi} ∫_{Γ} f̄(p)e^{pt} dp = 0  (5.60)

(as it must do, since f(t) = 0 for t < 0 by our initial assumption).

When t > 0, we close the contour to the left instead, and once again we can show that

∫_{C_{R}} f̄(p)e^{pt} dp → 0  (5.61)

as R → ∞, so long as f̄(p) → 0 as |p| → ∞. Hence in the limit R → ∞ we obtain

∫_{Γ} f̄(p)e^{pt} dp = 2πi ∑_{k=1}^{n} res_{p=p_{k}}(f̄(p)e^{pt})  (5.62)

by the Residue Theorem, where p_{1}, …, p_{n} are the poles of f̄(p). We deduce that

f(t) = ∑_{k=1}^{n} res_{p=p_{k}}(f̄(p)e^{pt})  (5.63)

for t > 0, so long as f̄(p) → 0 as |p| → ∞.

Examples:

(i) f̄(p) = 1/(p − 1). This has a pole at p = 1, so we must use γ > 1. We have f̄(p) → 0 as |p| → ∞, so Jordan’s Lemma applies as above. For t < 0, therefore, f(t) = 0, and for t > 0,

f(t) = res_{p=1} (\frac{e^{pt}}{p − 1}) = e^{t}.  (5.64)

This agrees with our earlier result for the Laplace transform of e^{λt} when λ = 1.

(ii) f̄(p) = p^{−n}. Here we need γ > 0, because there is a pole of order n at p = 0. For t < 0, f(t) = 0 as usual. For t > 0,

f(t) = res_{p=0} (\frac{e^{pt}}{p^{n}} ) = lim_{p→0} \{\frac{1}{(n − 1)!} \frac{d^{n−1}}{dp^{n−1}} e^{pt}\}

= lim_{p→0} \{\frac{1}{(n − 1)!}(t^{n−1}e^{pt})\}

= \frac{t^{n−1}}{(n − 1)!}.

(iii) What if f̄(p) ↛ 0 as |p| → ∞? Consider the example

f̄(p) = \frac{e^{−p}}{p};  (5.65)

here, as p → −∞ on the real axis, f̄(p) → ∞. We need to calculate

f(t) = \frac{1}{2πi} ∫_{Γ} \frac{e^{−p}}{p} e^{pt} dp,  (5.66)

but Jordan’s Lemma does not immediately apply. Note, however, that e^{−p}e^{pt} = e^{p(t−1)} = e^{pt′} where t′ = t − 1; so

f(t) = \frac{1}{2πi} ∫_{Γ} \frac{e^{pt′}}{p} dp.  (5.67)

Now we can use Jordan’s Lemma: when t′ < 0, close to the right, and when t′ > 0, close to the left, picking up the residue from the pole at p = 0. Hence

f(t) = \begin{cases}0 & t′ < 0, \\ 1 & t′ > 0 \end{cases}

= \begin{cases}0 & t < 1, \\ 1 & t > 1. \end{cases}

Worked Example: Solving Differential Equations using the Laplace Transform and its Inverse

We shall solve

ẍ + x = 2 sin t

for x(t), with initial conditions x(0) = 0, ẋ(0) = 2. Taking the Laplace transform with respect to time,

(p^{2}x̄(p) − px(0) − ẋ(0)) + x̄(p) = \frac{2}{p^{2} + 1}.

Using the initial conditions, we obtain

p^{2}x̄ − 2 + x̄ = \frac{2}{p^{2} + 1}

from which we deduce that

x̄ = \frac{2p^{2} + 4}{(p^{2} + 1)^{2}} .

To invert this we write down the Bromwich inversion formula

x(t) = \frac{1}{2πi} ∫_{γ−i∞}^{γ+i∞} \frac{2p^{2} + 4}{(p^{2} + 1)^{2}} e^{pt} dp.

The integrand has poles of order two at p = ±i, so we must have γ > 0 in order that the integration contour lies to the right of the singularities.

What are the residues at the poles? At p = i, the residue is

lim_{p→i} \frac{d}{dp} (\frac{2p^{2} + 4}{(p + i)^{2}} e^{pt}) = lim_{p→i} (\frac{(p + i)(4p + (2p^{2} + 4)t) − 2(2p^{2} + 4)(p + i)}{(p + i)^{3}} e^{pt})
= −\frac{1}{2}(t + 3i)e^{it}.  (5.68)

Similarly, at p = −i the residue is −\frac{1}{2}(t − 3i)e^{−it}.

As |p| → ∞, x̄(p) = O(|p|^{−2}) → 0; hence for t > 0 we close the integration contour to the left, picking up the residues from the poles to obtain

x(t) = −\frac{1}{2}(t + 3i)e^{it} − \frac{1}{2}(t − 3i)e^{−it}
= −\frac{1}{2}(2t cos t + 3i(2i sin t))
= 3 sin t − t cos t.  (5.69)

What function f(t) has Laplace transform f̄(p) = p^{− /_{2}^{1}}? We need to find

f(t) = \frac{1}{2πi} ∫_{γ−i∞}^{γ+i∞} p^{− /_{2}^{1}}e^{pt} dp.  (5.70)

For t < 0 we can close the contour to the right as usual and obtain f(t) = 0. For t > 0, however, the branch cut gets in the way.

Use a contour as shown, with a small circle of radius ε round the origin and two large quarter-circles of radius R. Substituting p = εe^{iθ} on the small circle gives a contribution of

∫_{π}^{−π} ε^{− /_{2}^{1}}e^{−iθ/2}e^{εe^{iθ}t} iεe^{iθ} dθ = O(ε /_{2}^{1}) → 0 as ε → 0.  (5.71)

Similarly, the integrals round the two large quarter-circles vanish as R → ∞, using the method used to prove Jordan’s Lemma. Hence the required integral is equal to the sum of the integrals on either side of the branch cut: i.e., for t > 0,

f(t) = \frac{1}{2πi} \{− ∫_{∞}^{0} r^{− /_{2}^{1}}e^{−iπ/2}e^{−rt}(−dr) − ∫_{0}^{∞} r^{− /_{2}^{1}}e^{iπ/2}e^{−rt}(−dr)\}
[substituting p = re^{iπ} and p = re^{−iπ} respectively]

= \frac{1}{2πi} \{2i ∫_{0}^{∞} r^{− /_{2}^{1}}e^{−rt} dr\}

= \frac{2}{π} ∫_{0}^{∞} e^{−s^{2}t} ds

[substituting r = s^{2}]

= \frac{1}{π} \sqrt{\frac{π}{t}}

= \frac{1}{\sqrt{πt}}.

So ℒ [t^{− /_{2}^{1}}] = \sqrt{π} p^{− /_{2}^{1}}. This is a generalisation of the result that ℒ [t^{n}] = n!/p^{n+1} to ℒ [t^{α}] = Γ(α + 1)/p^{α+1} where the Gamma function is defined by

Γ(α) = ∫_{0}^{∞} x^{α−1}e^{−x} dx  (5.72)

and can easily be shown to be equal to (α − 1)! when α is a positive integer.