Dynamics Lecture Notes

Spring Semester 2021

Prof Romeel Davé
Syllabus
• Part 1 (5 lectures)
L1 - Newton’s Laws. Projectiles with resistive drag.
L2 - Circular motion. Lorentz force.
L3 - Introduction to Ordinary Differential Equations.
L4 - Linear Momentum. Variable mass problems. Rocket equation.
L5 - Work & Power. Kinetic & Potential Energy. Conservative Forces.
• Part 2 (6 lectures)
L6 - Small oscillations. Simple Harmonic Motion. Pendulum.
L7 - Energy in SHM. Damped Oscillations.
L8 - Second order ODEs.
L9 - Forced SHM. Resonance.
L10 - Coupled oscillations. Matrix methods.
L11 - Normal modes. Double pendulum.
• Part 3 (6 lectures)
L12 - Central forces. Angular momentum. Orbit equation.
L13 - Orbital energies. Elliptical orbits. Kepler’s Laws.
L14 - Transfer orbits. Two-body orbits. Precession.
L15 - Hyperbolic orbits. Unbound comet.
L16 - Scattering. Cross-sections.
L17 - Hard-body scattering. Rutherford scattering.

• Part 4 (3 lectures)
L18 - Rotating frames. Coriolis & centrifugal forces.
L19 - Earth as a rotating frame. Foucault’s pendulum.
L20 - Rigid bodies. Inertia tensor. Rolling Motion.

• Part 5 (2 lectures)
L21 - Euler-Lagrange Equation. Lagrangian dynamics.
L22 - Examples using Lagrangian dynamics

1
Recommended Texts
There are no required textbooks for this class. All examinable material is
contained in these notes, plus the workshop assignments and homeworks.

RDG: “Classical Mechanics”,

R. Douglas Gregory, Cambridge University Press (2006).

JRT: “Classical Mechanics”, John R. Taylor, UCB (2005).

FE: “Introduction to Classical Mechanics”,

A. P. French & M. G. Ebison, Springer (1987).

FC: “Analytical Mechanics”,

G. R. Fowles & G. L. Cassiday, 7th Edition, Brookes/Cole (2005).

The first half of:

McComb: “Dynamics and Relativity”,
W. D. McComb, Oxford University Press (1999).

and for Simple Harmonic Motion:

French: “Vibrations and Waves”, A. P. French, CRC Press (1971).

Recommended Texts on Mathematical Methods

RHB: “Mathematical Methods for Physics and Engineering”,

K. F. Riley, M. P. Hobson, S. J. Bence, Cambridge University Press (1997).

Boas: “Mathematical Methods in the Physical Sciences”,

Mary L. Boas, Wiley (1966).

BD: “Elementary Differential Equations”,

W. E. Boyce & R. C. Diprima, Wiley (1965).

2
1 Introduction to Dynamics
Dynamics, also called classical mechanics, is a branch of physics that is con-
cerned with the study of forces and their effect on the motion of objects.
This is arguably the most essential physics class you will take, in terms of
real-world applicability, and the concept introduced here underlie essentially
all future physics courses.
The concept of force is fundamental to dynamics. Much of this course will
be about how objects move under different types of forces. The equations
that connect force to motion are called Newton’s Laws, first published in
Principia in 1687. The idea of solving a dynamics problem is typically to
apply Newton’s Laws (or relations derived from them) to a physical situation,
and obtain a mathematical description for how the object(s) in the problem
with move over time. In other words, solving a dynamics problem means
determining the position ⃗x(t) and/or velocity ⃗v(t) for each object in the
system.
Since most of this course will involve studying applications of Newton’s
Laws to a wide variety of physical systems to determine mathematically
how they respond under common types of forces, let us begin by reviewing
Newton’s Laws.

1.1 Newton’s Laws

1.1.1 Newton’s First Law: Definition of Inertial Frame
N1: When viewed in an inertial reference frame,
an object continues to move with a constant velocity ⃗v
unless acted on by a net external force.

• This law implicitly defines what is meant by an inertial frame. There

are infinite choices of inertial frames, each with its own ⃗v. The choice
⃗v = 0 is the rest frame of the object. Given two inertial frames S and
S ′ with constant relative velocity between them ⃗u, then ⃗v′ = ⃗v + ⃗u.
• Galilean invariance: The laws of physics are the same in all inertial
frames of reference.
• Uniform motion: ⃗v is constant, and the object travels in a straight line
at a constant speed.

3
1.1.2 Newton’s Second Law: Definition of Force
N̄2: ⃗ is defined as the rate of change of momentum ⃗p ≡ m⃗v:
F̄orce F
⃗ =
F d⃗
p
.
dt

• In the case where the inertial mass m is constant, we obtain the familiar
formula:
⃗ = m d⃗v = m⃗a,
F
dt
while the more general form is

⃗ = d⃗p = m⃗a + ⃗v dm
F
dt dt

• N2 applies within any chosen inertial reference frame. Within the frame
the acceleration of an object is attributed to the presence of a net force.

• Force and acceleration are vectors. So in three dimensions, N2 repre-

sents three distinct equations, in any chosen (orthogonal) coordinate
system.

• The force F
⃗ is the vector sum of all forces acting on the body. When
more than one force is present, it is known as a superposition of forces.

• Equilibrium is defined as several external forces acting on an object

such that the net force is zero.

• The force and acceleration do not have to be in the direction of the

original motion ⃗v. As we will see, a parallel force leads to linear accel-
eration (rectilinear motion), while a perpendicular force leads to cen-
tripetal acceleration (circular motion).

• Inertial mass is a scalar quantity that measures the reluctance of an

object to be accelerated by a force. Einstein’s (weak) equivalence prin-
ciple states that this is identical to the gravitational mass, i.e. the mass
that attracts via gravity.

• Force is measured (in mks) in units of Newtons: N = kg m s−2 , where

1 kg is defined by a reference mass in Paris.

4
1.1.3 Newton’s Third Law: Conservation of Momentum
N3: When one object exerts force upon another (F ⃗ 21 ), the other object
⃗ 12 that is equal in magnitude
object exerts a force on the first F
and opposite in direction:
⃗ 12 = −F
F ⃗ 21

• From N2, we have ddtp⃗1 p⃗2

= − ddt , which implies that the total momentum
in a (closed) system is constant: dtd (p⃗1 + p⃗2 ) = 0. N3 is thus equivalent
to the Law of Conservation of Momentum.

• For constant masses m1 and m2 , we have m1⃗a1 = −m2⃗a2 . In this case,

the acceleration of each object is inversely proportional to its mass.

• N3 is symmetric and simultaneous between the two objects. Neither

object initiates the mutual interaction. The common phrase “action
causes reaction” is misleading in this sense.

• The mutual force can be attractive or repulsive, corresponding to ac-

celerations of the objects towards each other or away from each other.

• For central forces the mutual force acts along the line between the
centres-of-mass (or charge) of the objects. Examples include gravity
and the Coulomb electrostatic force.

• For normal forces the mutual force acts perpendicular to the surfaces at
the point of contact. The actual force is being provided by electrostatic
bonds in the surface material.

• For resistive forces such as friction and drag, the mutual force acts
along the line of relative motion between the objects, and opposite to
the direction of motion.

• N3 can be extended to the interactions of three (or more) objects. For

each object, the net force is the vector sum of the mutual interactions
with each of the other objects. This result assumes that the mutual
interactions are independent of each other.

5
1.2 Modelling Resistive Drag
We will now discuss some example applications of Newton’s Laws. Our first
example will be motion under a drag force.

1.2.1 Dependence on Velocity

An object moving through a medium experiences a drag force, always in
the direction opposing the motion. In the case of a plane slab face moving
through a medium, the magnitude of the force is proportional to v 2 (or ⃗v ·⃗v).
This is because each atom of the medium imparts a momentum ∝ v, while
the number of atoms intercepted by the face in a given time is also ∝ v.
Thus the momentum change in the medium is ∝ v 2 . By N3, the momentum
change in the object (and thus by N2, the force acting on it) must also be
∝ v2.
In the case of laminar (surface) flow, the magnitude of the force is typically
proportional to v (or |⃗v|). Here, the medium is not being displaced, but only
provides a surface resistance as the object passes. The faster the object
moves, the linearly more surface interactions it has, and hence the drag
depends linearly on v.
The purpose of making cars and planes more aerodynamically-shaped is
to be closer to the laminar regime of ∝ v rather than ∝ v 2 , so the drag stays
modest as the velocity increases.

1.2.2 Linear (Laminar) Drag Model

Let us consider laminar drag first, as it is mathematically simpler. Here,
⃗ = −k⃗v. Note that the coefficient k has dimensions of kg/s. Since the drag
F
only operates in the direction of motion, we can work in one linear dimension,
namely the direction of motion x. We assume the object’s mass is constant.
In the 1-D case, N2 gives
dvx
m = −kvx
dt
We can rearrange terms to put all the vx terms on one side, and t on the
other, and take the integral of both sides:
Z
dvx k Z
=− dt
vx m

6
which has the solution
kt
ln(vx ) = − +C
m
To evaluate the constant of integration C, we require some boundary condi-
tion. Let us assume a given initial velocity of vx,0 at t = 0:

vx = vx,0 e−kt/m

Integrating with respect to time with x = x0 at t = 0 gives the position:

vx,0 m
x(t) = x0 + [1 − e−kt/m ]
k
Hence the object starts out at vx (0) (as we assumed), and then slows down.
The object only stops as t → ∞, after travelling a distance vx (0)m/k from
x(0). Note that even though it takes infinite time to stop, the distance to
stop is finite!
Adding gravity: Now instead of the horizontal case, let us consider the
vertical case with laminar drag. Here, an object is let go from rest, and
it starts to free-fall, but there is a laminar drag that slows the fall. Here,
there is a superposition of two forces, mg downwards from gravity, and kvz
upwards from drag.
To write down the equation of motion, one has to be very careful of the
signs! It is typical to work in the convention that upwards corresponds to
+z, and thus +vz and +az . Hence in this case, vz and az will be negative,
since the object is falling and being accelerated downwards.
Gravity pulls downwards, so we assign it a negative sign, −mg. We will
use the convention that g is a scalar with a value +9.8 m/s2 on the surface
of the Earth. The drag force remains −kvz . Although this has a negative
sign, the force will actually be positive, since vz will be negative as the object
falls; hence −kvz will act positively (upwards) as our intuition says it should.
Hence we must solve the following first order differential equation:
dvz
m = −mg − kvz
dt
As before we can treat the differentials dvz and dt as variables, and rewrite
as
dvz
= −dt
g + (k/m)vz

7
This is known as separation of variables; we’ve separated all the vz terms
onto the LHS, and all the t terms on the RHS. This allows us to integrate
both sides:
Z
dvz Z
m
= − dt =⇒ ln [g + (k/m)vz ] = −t + C
g + (k/m)vz k
If we specify an initial condition vz = vz,0 at t = 0, we can determine the
constant of integration C:
m
C= ln [g + (k/m)vz,0 ]
k
Hence
m
−t = (ln [g + (k/m)vz ] − ln [g + (k/m)vz,0 ])
k
g + (k/m)vz
e−kt/m =
g + (k/m)vz,0
mg  −kt/m 
vz = vz,0 e−kt/m + e −1
k
This is then the equation of motion.
Let us examine some limits. As t → ∞, e−kt/m → 0, so the object
approaches a terminal velocity:
mg
vterm = − .
k
At terminal velocity, the balancing of gravity and drag forces leads to an
equilibrium, and thus by N1, uniform motion.
One can also compute the acceleration by differentiating with respect to
(w.r.t.) time:
dvz
az = = − [g + (k/m)vz,0 ] e−kt/m
dt
As the velocity approaches terminal, the acceleration (which is always nega-
tive) approaches zero exponentially.
Two-dimensional motion: Now imagine a projectile shot at an angle
θ with a linear drag force. Since forces are vectors, we can separate the
problem into a horizontal (x) component and a vertical (z) component.
N2 in this case is written
⃗ = −k⃗v − mgẑ.
F

8
This represents two equations, one each for x (horizontal) and z (vertical)
directions.
In fact, these are exactly the two problems we have solved above: The x
direction is given by the first case, where there is no gravity. The z direction
is given by the gravity case. In this case, these two equations are completely
independent, since the linear drag force acts in each direction independently.
Note that this is only true because the linear drag force depends only on the
velocity in that particular direction.
Two equations now require two initial conditions. For instance, vx,0 =
v0 cos θ and vz,0 = v0 sin θ, specified by two variables v0 (the speed at launch)
and θ. From this, one obtains the full solution.

1.2.3 Quadratic Drag Model

Now let us consider the case of a quadratic drag force, where the force is given
by −k ′ |⃗v|2 , with a quadratic drag coefficient k ′ (having dimensions kg/m).
Horizontal: As before, horizontally the only force is drag:
dvx Z
dvx k′ Z
m = −k ′ vx2 =⇒ = − dt
dt vx2 m
1 k′t
=−
− +C
vx m
Putting in an initial condition v = vx,0 at t = 0 and rearranging:

1 1 k′t vx,0
− =− =⇒ vx =
vx,0 vx m 1 + (k ′ vx,0 /m)t
Integrating one more time gives the displacement:
m
x = x0 + ′
ln[1 + (k ′ vx,0 /m)t]
k
This result differs from the linear drag model. The velocity goes to zero like
1/t rather than exponentially, and the displacement increases logarithmically.

Vertical: For vertical motion under gravity there are two forces, mg down-
wards and k ′ vz2 upwards:
dvz
m = −mg + k ′ vz2
dt
9
Note that the sign of the drag term is opposite to the linear case, since vz2 is
always positive and the drag force acts in the positive (upwards) direction.
It is not as easy to solve this, but we can still obtain the terminal ve-
locity just by setting dvdtz = 0, since we intuitively know terminal velocity
corresponds to zero acceleration. This gives
mg
r
vterm = −
k′
The full motion requires a more difficult integral which yields a hyperbolic
tangent function for vz . Assuming it starts at rest (i.e. vz (t = 0) = 0), then

vz = vterm tanh(gt/vterm )

2-D projectile: When using a quadratic drag force to describe the two-
dimensional motion of a projectile it is necessary to write down the vector
form of the drag force correctly:
⃗ drag = −k ′ |⃗v|2 v̂ = −k ′ |⃗v|⃗v
F
q
The magnitude |⃗v| = vx2 + vz2 , depends on both vx and vz . This leads to
coupled differential equations for the two components vx and vz :
dvx q
m = −k ′ vx2 + vz2 vx
dt
dvz q
m = −mg − k ′ vx2 + vz2 vz
dt
These cannot be solved analytically, and numerical methods must be used.

Note how quickly even a simple scenario already leads to an analytically

intractable problem! This is a common circumstance, and highlights the im-
portance of numerical methods in dynamics. In this course we will not delve
into this, but rather focus on the suite of problems that have analytic solu-
tions. This is already sufficiently rich and full of important results to keep
us occupied for the term.

10
1.3 Circular Motion
In the drag force case, the motion was rectilinear, i.e. the direction of the
force was aligned with the motion. We now consider several examples forces
that are perpendicular to the velocity. Such a force is called centripetal,
because as we will see it results in circular motion.
In contrast to being given a force and solving for the motion as in the
previous drag case, let us approach this dynamics problem in the opposite
way: Let us assume an object is in circular motion at a constant radius r
and a constant angular velocity ω, namely:
dr dθ
= 0; ω= = constant =⇒ θ = ωt
dt dt
and then ask: What force is required to generate this motion?

Coordinate choice is important – any coordinate system is valid, but the

right choice can greatly simplify the solution. Here, it should be evident that
circular motion is best described in polar coordinates r̂ and θ̂:

r̂ = cos θî + sin θĵ

ˆ
θ̂ = − sin θî + cos θĵ,
where î and ĵ are the usual basis vectors in x and y.

Figure 1: Polar coordinates, with unit vectors.

11
 Since θ = ωt is time-dependent, the polar basis vectors r̂ and θ̂ are also
time-dependent (unlike î and ĵ). They rotate round with the motion, with r̂
perpendicular to and θ̂ parallel to the tangential velocity ⃗v. The derivatives
of the basis vectors with respect to θ are:
dr̂ dr̂ dr̂ dθ
= − sin θî + cos θĵ = θ̂ =⇒ = = ω θ̂
dθ dt dθ dt
dθ̂ dθ̂ dθ̂ dθ
= − cos θî − sin θĵ = −r̂ =⇒ = = −ωr̂
dθ dt dθ dt
This allows us to simplify the derivatives of polar basis vectors. These rela-
tions will be used throughout the term.

Using these, we want to determine the force corresponding to constant cir-

cular motion. To get a force, we need to determine m d⃗
v
dt
, for which we need
d⃗r
⃗v = dt in polar coordinates.
We can write the radius in vector form as ⃗r = rr̂, where r ≡ |⃗r| is a
constant in our case. To get the tangential velocity ⃗v, we take the time
derivative of ⃗r:
d⃗r dr̂
⃗v = =r =⇒ ⃗v = rω θ̂.
dt dt
The velocity is thus tangential (θ̂), and has a constant magnitude rω.
This can be generalised to the 3-D case of cylindrical coordinates, with an
axis ẑ perpendicular to the plane of circular motion. Here, the vector (cross)
product is useful, because it picks out the plane in which the circular motion
is happening. Hence the tangential velocity straightforwardly generalises to
⃗ × ⃗r
⃗v = ω
The direction of ω
⃗ is given by the corkscrew (or right-hand) rule. It is up
(+ẑ) for counter-clockwise motion, down (−ẑ) for clockwise motion. You
can convince yourself that the cross product leads to a velocity in the corre-
sponding direction.

Given ⃗v, we can obtain an acceleration by differentiating w.r.t. time:

d⃗v dθ̂
⃗a = = rω ,
dt dt
hence

12
 ⃗a = −rω 2 r̂
The acceleration is centripetal (towards the centre), i.e. in the −r̂ direction,
with a constant magnitude rω 2 .
This can also be generalised in terms of vector products, where we note
⃗ is a constant and hence only ⃗r differentiates into ⃗v:
that ω
d⃗
v
⃗a = dt
⃗ × ⃗v = ω
=ω ⃗ × (⃗ω × ⃗r)

Note that the order of cross products specified by the parantheses is impor-
tant! It is not the same as (⃗ω × ω
⃗ ) × ⃗r = 0.
Thus we arrive at the force required for circular motion: a constant
centripetal (inwards) force of F = mrω 2 produces circular motion. Using
ω = v/r, we can also write this as F = mv 2 /r.

Centrifugal force: The above equations are valid within any inertial frame
of reference, where N2 is valid. But in many cases, we would like to solve the
problem in the object’s frame (i.e. as if we were sitting on the object) as it
moves around in a circle. Trouble is, this is no longer an inertial frame! Can
we still use N2 somehow?
The answer is yes, by introducing a “fictitious force” known as centrifugal
force. A fictitious force is an omnipresent force felt by objects in a non-inertial
frame. As you know, if you place a ball down on a rotating carousel, it will
not stand still, it will seem to feel a force “coming from nowhere” that pushes
it outwards. This is centrifugal force.
It should be evident that the centrifugal force must be equal and opposite
to the centripetal force required to keep the frame rotating. In the carousel
case, even though the ball feels no force, the frame itself feels a force, so
within that frame the ball will feel the opposite force.
Thus to do dynamics calculations in a non-inertial frame, we simply have
to add the appropriate fictitious force into N2, and then we can happily do
calculations as if we are in an inertial frame. In the case of a frame rotating
with constant circular motion, we must add a centrifugal force mω 2 r (or
mv 2 /r) that acts outwardly from the axis of rotation.

1.3.1 Example: Car going round a banked bend

Let us use N2 including a centrifugal force to solve a dynamics problem. A
car makes a turn at a constant speed, which can be approximated as circular

13
motion with a radius r through the turn. As is often the case on highways,
the road has been banked at an angle θ through the turn. Intuitively, you
probably know that if the car goes too slow, it will slip down the bank towards
the turn, while if it go too fast it will experience a force outwards from the
turn. We want to know: How fast must the car go in order not to slip off the
road?

Figure 2: Force diagram for a car on a banked curve. In the frame of the
car, there is additionally a centrifugal force that acts horizontally outwards
with an amplitude mv 2 /r.

Frictionless road: Here we assume the road is frictionless (e.g. covered

in ice). In dynamics problems, it is always beneficial to begin by drawing a
force diagram, as shown in the Figure 2. This shows the weight of the car mg
acting downwards, and a normal force N acting perpendicular to the banked
slope. If we now make this bank as part of a turn that the car does with
2
a velocity v, there is also a centripetal force F = mvr acting horizontally
inwards. To not slide along the bank, we need all forces to balance.

14
 Since our problem is 2-D, this involves two separate equations we must
balance: vertical and horizontal. For the vertical case, gravity downwards
must balance the component of the normal force upwards, which is N cos θ:
mg
mg = N cos θ =⇒ N =
cos θ
This gives us N , which is otherwise not specified.
Next, we do the horizontal case. Here, the centrifugal force must be
balanced by the horizontal component of the normal force:

mv 2
= N sin θ = mg tan θ
r
where we have substituted N from the vertical equation. Solving for v, we
get q
v = gr tan θ
To be concrete, let us take r = 300 m, θ = 30◦ , and g = 9.8 m/s2 . Then

v 2 = (300 m)(9.8 m/s2 )(0.5774) ≈ 1700 m2 /s2 =⇒ v = 41 m/s

This is 150 km/h or about 90 mph, which is quite fast since a 30◦ bank is
quite steep (comparable to what is found on NASCAR racetracks). Note
that the car’s mass m is irrelevant.

With friction: In reality, tyres provide a resistance against sideways slip-

page. We can include this as a frictional force between the tyres and the
road, with a coefficient of static friction µ. This adds a frictional force of up
to µN , opposing the direction of slippage (perpendicular to N). ⃗ It is static
friction because even though the car is rolling forwards, it is not moving
sideways (we hope), which is the direction we are concerned with here. With
friction, now what range of speeds can the car go around the bend without
sliding?
√
There are two cases here: One where the velocity is above gr tan θ, in
which case the frictional force points inwards to prevent slippage outwards,
and the other case where the velocity is below the threshold so the frictional
force points outwards. This can be represented by placing a − or + in front
of µ, where we assume outwards is positive as before.

15
 The vertical equation, now including a frictional force perpendicular to
⃗
N, becomes
mg
mg = N cos θ ∓ µN sin θ =⇒ N =
cos θ ∓ µ sin θ
In the horizontal plane the centripetal force is balanced by

mv 2
" #
sin θ ± µ cos θ
= N sin θ ± µN cos θ = mg
r cos θ ∓ µ sin θ

Solving for v and dividing through by cos θ we have

" #
2 tan θ ± µ
v = gr
1 ∓ µ tan θ

If we take µ = 0.4 and other quantities as before, then the − and + cases
give values of 21 m/s and 61 m/s, respectively. Thus with this frictional force
present there is a range of safe speeds between 21 and 61 m/s where the car
won’t slide in either direction. Without the frictional force (i.e. using µ = 0
in the above equation) there is only one safe speed of 41 m/s.

1.3.2 Example: Motion in Electric and Magnetic Fields

Another classic problem resulting in circular motion arises in electromagnet-
⃗ and magnetic field
ics. An electric charge q inside a constant electric field E
⃗ experiences a Lorentz force given by
B
⃗ = qE
F ⃗ + q⃗v × B
⃗

⃗ = 0 the Lorentz force is a cross product of the velocity, so it

In the case of E
will always be perpendicular to ⃗v. This means it is a centripetal force, and
⃗ direction. For a constant B,
will cause circular motion about the B ⃗ this will
be constant circular motion at a velocity v and frequency ω. What will be v
and ω for this charge?
We work in Cartesian coordinates, where B ⃗ ≡ Bz ẑ, and the velocity v
around a circle in the x − y plane. The Lorentz force qvBz provides the
centripetal force mrω 2 = mvω (using v = ωr):

qBz
qvBz = mvω =⇒ ωc =
m
16
Figure 3: Motion of a charge in parallel magnetic and electric fields along
the y axis. The helix shold be getting stretched out as the particle moves,
which is not depicted well in this figure.

This is known as the cyclotron frequency. The velocity is then ωc r, where r

corresponds to the radius of the ring in which the charge moves.
This is the basic idea of a particle accelerator, like a cyclotron. A very
large ring is lined with strong magnets, which results in a high ωc and hence
a high v that can come very close to the speed of light. Oppositely charged
particles move in beams in opposite directions, until they are brought to-
gether in a high-energy collision that produces many new particles.

Now let us add an electric field. In general, it can be in any direction relative
to the magnetic field. But let us consider the two cases where E ⃗ and B ⃗ are
parallel and perpendicular, respectively.

⃗ ∥ B:
E ⃗ Here, E ⃗ = E k̂, so the electric field accelerates the charge along
the z-axis. Since B ⃗ only acts in the x − y plane, nothing is changed in re-
gards to the magnetic field, which still causes cyclic motion at a frequency
⃗
ωc . In the k̂ direction, meanwhile, there is now an acceleration q E/m. Hence
the motion is helical, with the helix getting stretched out vertically as the
charge is accelerated.

⃗ ⊥ B:
E ⃗ Here the situation is more complex. Let us take B ⃗ = Bz k̂ (so
⃗ = Ey ĵ. In this case, for the x direction we have no
Bx = By = 0), and E

17
electric force, but we have a magnetic force:
dvx
Fx = m = qvy Bz
dt
In the y direction we have both a contribution from Ey and Bz :
dvy
Fy = m = q(Ey − vx Bz )
dt
Note that these equations are coupled, in the sense that the vy differential
depends on vx and vice versa.

Solving this requires a bit of cleverness. One can show that the solution
to these simultaneous equations with the initial condition vy (t = 0) = 0 is
Ey
vx = −v cos ωt +
Bz
vy = v sin ωt
where v is the constant tangential velocity. This corresponds to uniform
linear motion in the x direction given by the ratio of electric to magnetic
fields, superimposed on circular motion due to the magnetic field.
For the special case of v = Ey /Bz , vx varies between 0 at t = 0 and 2v at
ωt = π, and vy varies between +v at ωt = π/2 and −v at ωt = 3π/2. This
motion is known as a cycloid. A nice animation of various cases is shown at
[Link]

1.4 Summary: Solving for the Equation of Motion

Here is a good plan of attack for solving dynamics problems to obtain the
equations of motion ⃗x(t) and/or ⃗v(t):

1. Sketch the system, and draw a force diagram showing all forces acting
on the body, and their directions. Remember to choose a convenient
coordinate system to make your life easier.

2. Write down N2 for each direction, by separating the forces into compo-
nents in your chosen coordinates. Problems in N dimensions will result
in N equations. Each equation will represent an ordinary differential
equation (ODE).

18
 3. Solve the system of ODEs. This typically involves manipulating the
ODE into an integral form, and then solving the integrals.

4. Use supplied boundary conditions to find the constants of integration.

5. From the solution to each ODE, obtain the equation of motion ⃗x(t)
and ⃗v(t), or whatever quantities are desired.

6. Do a sanity check that the object behaves as you would expect physi-
cally; it is often useful to examine limiting cases for this.

Much of this course will involve discussing both specific and generalized
approaches to solving such dynamics problems in a variety of circumstances.
In the next section, we discuss Step 3 of this process, namely systematic
approaches to solving ODEs.

19
2 Ordinary Differential Equations
Ordinary Differential Equations (ODEs) are used to describe a vast range of
scientific phenomena. As we have already seen, they are particularly impor-
tant in dynamics, where motion is typically described by functions of time
and their derivatives. ODEs are the subject of four chapters (12 to 15) in
Riley, Hobson & Bence (RHB). The first few chapters of Boyce & Diprima
(BD) are also recommended.
As a reminder, solutions to integrals can be found in many places, so you
don’t normally need to solve them yourself. Some online tools are:
[Link] (Mathematica)
[Link] (Maple)
For exam purposes you will be expected to handle standard integrals, but
for more complicated cases the result will be given to you.

2.1 Introduction to ODEs

An ODE is a functional relation between a dependent variable y and its
derivatives with respect to a set of independent variables xi .
dn y
!
dy
f y(xi ), , ... n , xi = 0.
dxi dxi
In dynamics, y usually becomes a vector quantity representing the three spa-
tial dimensions (e.g. ⃗r), leading to three equations in the three components,
which may be separate or coupled depending on the forces involved. In most
cases (but not all), the derivatives will be taken with respect to one variable,
time xi = t.
d⃗r dn⃗r
!
f ⃗r, , ... n , t = 0
dt dt
For the moment we simplify the discussion to one dimension y, i.e. we assume
that a general three-dimensional dynamics problem is separable, and can be
reduced to one or more one-dimensional equations. However, this is not
always possible, as we saw in the Lorentz force problem; in this case, the
approach is often to manipulate the coupled ODEs into a separable form.
The solution of an ODE is the most general function y(t) that satisfies
the ODE. The most general aspect is important, because it encompasses
all possible solutions to the ODE, often added in superposition. Only after
applying external information can some parts of the solution be excluded.

20
 Usually, to solve an ODE, the idea will be to convert it into one or more
integrals, since these can be solved straightforwardly (if not always easily).
Each order of the ODE will result in one integration to reduce the order of
the differential, which thus requires one constant of integration. As such, an
nth order ODE will require n constants of integration. If there are m dimen-
sions with m separate ODEs, then each will require n separate constants of
integration to fully specify, for a total of mn constants of integration.
The constants of integration are determined by the application of suitable
boundary conditions, which specify the state of the system at some time t.
Boundary conditions at t = 0 are often called initial conditions.
A solution y(t) can always be checked by taking the derivatives and sub-
stituting back into the original ODE. This can be a useful way to figure out
what you did wrong if your resulting y(t) doesn’t make physical sense.
In this section we discuss some analytical methods for solving first order
ODEs. Usually there are several different methods for arriving at the same
solution, all of which are equally valid, and it is your choice which is the
easiest.

2.1.1 Classification of ODEs

Some definitions and terminology to set the stage:
• The order n of an ODE is the highest derivative dn y/dtn that appears.
First order ODE stop at dy/dt, while second order ODE go up to
d2 y/dt2 . Examples of first and second order equations are radioactive
decay and Newton’s second law, respectively:
dN (t)
= −λN (t) radioactive decay
dt
d2⃗r ⃗
!
d⃗r
m 2 = F ⃗r, , t N2
dt dt
• A homogeneous ODE is one which does not have terms which depend
only on the independent variable t:
dy dy
+P (t)y = 0 homogeneous; +P (t)y = Q(t) inhomogeneous
dt dt
This is a 1st order example, but this classification applies to any order
ODE. We will use this P and Q notation for first order ODEs in the
following sections.

21
 • An ODE is linear if neither y nor the derivatives dn y/dtn is raised to a
power, or appears in a product with another derivative. Linear ODEs
are generally the easiest to solve. Non-linear ODEs can contain terms
such as (dy/dt)2 , y 2 or y(dy/dt), such as the logistic equation which
leads to the well-known S-curve describing population growth:
dN (t) N
 
= rN 1 − logistic equation
dt K

• Partial differential equations (PDEs) have derivatives with respect

to more than one variable. An example involving ∂y/∂x and ∂y/∂t is
the wave equation:
∂ 2y 1 ∂ 2y
= 2 2 wave equation
∂x2 v ∂t
where y ≡ y(x, t). Note the symbol for a partial derivative ∂, which
indicates that the derivative is only to be taken w.r.t. that variable,
assuming all other variables are constant. We will not discuss PDEs
(much) in this course.

2.2 Separation of variables

The easiest ODE to solve is a separable, 1st order ODE. This type of ODE
has the form:
dy
= f (t)g(y)
dt
Such an equation can be directly solved by integrating both sides:
Z
dy Z
= f (t)dt + C
g(y)
where a single constant of integration C is determined by a boundary condi-
tion. Of course the integrals might be complicated, but integrals are generally
solvable so there is a path to a solution. A goal of solving ODEs is often to
reduce more complex ODEs into a separable (integrable) form.
Radioactive decay: Consider a population of N atoms with decay rate
λ. The decay rate specifies the fraction of remaining atoms that will decay
in time dt. Hence the ODE is
dN (t) Z
dN Z
= −λN (t) =⇒ = − λdt + C =⇒ ln(N ) = −λt + C
dt N
22
The integration constant is given by the initial population N0 at t = 0, giving
C = ln N0 . The particular solution is then N (t) = N0 e−λt .

Generalised drag: We return to the drag problem in the previous chapter,

but now have a generalised drag force including both a laminar and quadratic
term (but no gravity). The velocity v of an object is thus described by:

dv
= av − bv 2
dt
where a and b are non-zero coefficients with dimension s−1 and m−1 respec-
tively. This is solvable by separating variables:
dv
= dt
v(a − bv)

With two factors including v in the denominator, it is useful to use the trick of
partial fractions, which separates out the factors into their own denominators:
!
1 1 1 b
= + .
v(a − bv) a v (a − bv)

Now, the integral is more straightforward, comprising of two integrals of the

form dx/x. This yields, for the LHS integral:
Z
dv 1 Z dv 1 Z bdv 1 1
= + = ln(v) − ln(a − bv)
v(a − bv) a v a (a − bv) a a

Now we equate this to the RHS integral:

v a − bv a
ln = at − ln(K) =⇒ = Ke−at =⇒ v = .
a − bv v b + Ke−at
As t → ∞, this results in a constant terminal velocity of a/b.

EXTRA: Let us apply some initial conditions and relate it to the results
from our previous drag force calculations. We assume v(t = 0) = bca (for
reasons that will become clear), where c is a dimensionless constant. In that
case !
a 1
K = b(c − 1) =⇒ v =
b 1 + (c − 1)e−at

23
For the special case c = 1 we have uniform motion at v = a/b. For c < 1
(c > 1) the velocity reduces (increases) asymptotically towards v = a/b. For
any initial velocity the terminal velocity is always the same. This can also
be seen directly by setting dv
dt
= 0 in the original ODE.
To relate this to the drag examples in the previous chapter, we can set
a = −k/m and b = k ′ /m. In this case c should be chosen to be negative, and
the initial positive velocity goes asymptotically to zero due to the positive
exponential term.

2.3 Solution by substitution (i.e. change of variables)

If the ODE is not separable, a common trick is to try to change variables to
obtain a form that is separable (or can otherwise be easily solved). One can
change either the dependent variable (the one in the top differential) or the
independent one. This is where solving ODEs veers into an art, as one has
to cleverly choose the right variable to substitute.
One can change the independent variable by using the chain rule for
differentiation. For example, suppose the acceleration is specified as dv/dt =
f (x) rather than f (t), where x(t) is the position variable and v = dx/dt.
Then one can use the chain rule to get:
dv dv dx dv
f (x) = = =v
dt dx dt dx
which leads to the separable form:

vdv = f (x)dx

which you can integrate to solve for v(x) and thus v(t).

More generally, the RHS of the linear ODE can be a function of both y
and t, as in
dy
= f (y, t)
dt
Here we change the dependent variable from y(t) to some other cleverly-
chosen variable u(y, t). Again, we can apply the chain rule:

du du dy du
= = f (y, t)
dt dy dt dy

24
The trick is to choose u such that the RHS comes into a form that is more
tractable, typically a separable ODE.
Solving this will give u(t), which then must be substituted using y =
u(y, t) to get y(t). Do not forget to change back to the original dependent
variable y!

Let’s consider some particular cases.

• The form f (y/t): Suppose a first order homogeneous ODE can be written
in the form:
dy y
 
=f .
dt t
A good choice for a change of variable is u = y/t (or y = ut). Note that u is
in general a function of t, it is not a constant!
For example, consider
dy y 2 + 2yt
=
dt t2
This can be written in terms of ratios of y/t:
 2
dy y y
 
= +2
dt t t
The change of variable is
dy du
y = ut =⇒ =u+t
dt dt
which results in the equation
du du
u+t = u2 + 2u =⇒ t = u2 + u
dt dt
This is now a separable equation! We can thus write
Z
du Z
dt
=
u(1 + u) t

and integrate the LHS by partial fractions:

ln(u) − ln(1 + u) = ln(t) + ln(C)

25
where we choose to write the constant of integration as a natural log. Ex-
pressing the result in terms of the original variable y:

u y Ct2
= = Ct =⇒ y =
1+u y+t (1 − Ct)

• The form f (at + by + c): Here we have:

dy
= f (at + by + c)
dt
where a, b and c are constants, and f is any function of the given linear
combination of the variables. This can be solved by a change of variable
u = at + by + c.
du dy
= a + b = a + bf (u)
dt dt
This is now separable and can be integrated:
Z
du Z
= dt
a + bf (u)

• From the above examples, you can see that if you have a form dy/dt =
f (g(y, t)), it is often useful to choose u = g(y, t).

2.4 Integrating factors

For the case of a linear first order ODE, there is a fully general way to solve
the problem, using what is called an integrating factor. Often it is easier to
just think of a clever substitution, but if one is at a loss for this, it is always
possible to fall back on an integrating factor.
To see how this works, let us take our linear 1st order ODE and multiply
it by some (to be determined) function µ(t):
" #
dy
µ(t) + P (t)y = µ(t)Q(t)
dt

µ(t) is known as an integrating factor, and must be chosen so that it satisfies:

d dy dµ(t)
µ(t)Q(t) = [µ(t)y] = µ(t) + y .
dt dt dt
26
How can we choose µ(t)? Don’t worry, we don’t have to guess. Let us go
back to our ODE, but now substitute the RHS with our µ(t) equation:
" #
dy dy dµ(t)
µ(t) + P (t)y = µ(t) + y
dt dt dt

Note that the first term is the same on the LHS and RHS, and thus they
cancel! Further dividing by y, we obtain:

dµ(t) Z
dµ(t) Z Z
P (t)µ(t) = =⇒ = P (t)dt =⇒ ln[µ(t)] = P (t)dt
dt µ(t)

So an integral over P (t) determines the required integrating factor µ(t):

R
µ(t) = exp [ P (t)dt]

Note that we have discarded the constant of integration here, which would be
a multiplicative factor to the exponential. This is with no loss of generality,
because recall that we are multiplying the entire ODE by µ, so any constant
factor must cancel out and not impact the solution.

Now that we know µ(t), we return to the equation where we defined µ,

and integrate it: Z
µ(t)y = µ(t)Q(t)dt + C
Here we must include a constant of integration C, since we are now solving
for y! The solution is thus:

1 R
y= µ(t)
[ µ(t)Q(t)dt + C]

Notice that the constant of integration appears inside the brackets – do not
forget this! This is important because C is multiplied by µ−1 (t), and so
cannot be added in to y at the very end.
In short, by using this magic integrating factor µ(t), we are able to con-
vert any linear 1st order ODE into two integrals, one to obtain µ(t), and
another to obtain y by integrating µ(t)Q(t). This method will always give
you a solution with no guesswork, so long as you have a linear first order ODE.

27
Example: Consider the linear 1st order ODE
dy
+ 2yt = 4t
dt
Identifying P = 2t and Q = 4t, the integrating factor is:
Z
2
ln[µ(t)] = 2tdt = t2 =⇒ µ(t) = et
Z 
2 2
y = e−t et 4tdt + C
2
The integral here gives 2et , so the final result is:
2
y = 2 + Ce−t

2.5 Non-linear ODEs: Bernoulli’s equation

For non-linear ODEs, things get trickier, and the optimal substitution is not
always evident. However, there is one form that appears commonly and is
amenable to a clever substitution, known as Bernoulli’s equation:
dy
+ P (t)y = Q(t)y n .
dt
If n = 0, we have a linear ODE that we can solve by substitution or an
integrating factor, while if n = 1 we have a separable ODE. But here we
consider the general case where n is any real number – not necessarily an
integer – besides 0 or 1.
This can be linearised (i.e. turned into a linear ODE) via the substitution
u = y 1−n . For this, the change of variable gives
du dy dy y n du
= (1 − n)y −n =⇒ =
dt dt dt 1 − n dt
Note that y = y n u. So we can substitute back into the original ODE:
y n du
+ P (y)y n u = Q(t)y n
1 − n dt

28
The magic of Bernoulli’s substitution is that the y n term can now be divided
out of the entire equation. This leaves
1 du
+ P (t)u = Q(t) n ̸= 0 or 1.
1 − n dt
which is now a linear ODE!

Example with n = 2: In this case, the change of variable is u = 1/y.

dy
+ P (t)y = Q(t)y 2
dt
1 du dy 1
u= =−
y dt dt y 2
and the ODE is linear in the new variable u:
du
− P (t)u = −Q(t)
dt
√
Example with n = 0.5: Here, we set u = y:
dy √
+ P (t)y = Q(t) y
dt
√ du dy 1
u= y = √
dt dt 2 y
du
2 + P (t)u = Q(t)
dt
In each case, we have reduced the nonlinear ODE to a linear one that can be
solved by the usual means.

2.6 The Logistic Equation

As a example, consider the Logistic equation for population growth:
dN (t) N r 2
 
= rN 1 − = rN − N
dt K K
where r and K are constants. In population growth, r represents the repro-
duction rate, while K represents some limiting factor (e.g. food) that must

29
be shared among the population of N . The Logistic equation determines
how the population N evolves with time under these conditions.
This equation can be solved via separation of variables or Bernoulli’s
equation. Let us see how each works.

• Separation of variables: We can write the ODE as

Z
dN Z
= r dt
N − N 2 /K

Via partial fractions, the left integrand is

K 1 1
 
= + .
N (K − N ) N K −N

Integration leads to the solution

K −N K
ln N − ln(K − N ) = − ln = rt + C =⇒ − 1 = −Ce−rt
N N
Thus the solution is
1 1 K
= (1 − Ce−rt ), or N (t) = .
N (t) K 1 − Ce−rt

We will set some initial conditions to determine C in a minute.

• Bernoulli’s equation: We recognize that the Logistic equation is a Bernoulli

equation with n = 2. We identify P (t) = −r and Q(t) = −r/K, which are
constants. Hence the required substitution is u = 1/N , which yields:

du r du
− − ru = − =⇒ K = r(1 − Ku)
dt K dt
This is now a first order linear ODE in u, and is separable:
Z
Kdu Z
= rdt
(1 − Ku)

− ln(1 − Ku) = rt − ln(C) 1 − Ku = Ce−rt

1 1
u= = (1 − Ce−rt ).
N K
30
This is the same solution as before.

• Integrating factor: The u ODE above can also be solved using an inte-
grating factor. With P (t) = r and Q(t) = r/K, the integrating factor is:
Z 
µ(t) = exp P (t)dt = ert

and the solution for u is:

1 r rt C 1 1 1
Z   
u= = e−rt e dt − = re−rt ert − Ce−rt = (1 − Ce−rt ),
N K K K r K
again as before.

To evaluate C, let us assume an initial population of N0 , i.e. u(t = 0) = 1/N0 .

Then
1 1 K
= (1 − C) =⇒ C = 1 −
N0 K N0
So
K
N (t) =
1 − (1 − K/N0 )e−rt
As t → ∞, N → K which is the saturation population (see Figure 4). The
parameter r determines the rate of the population growth at small t.

Figure 4: Solution to the logistic equation modeling a population N (t). N (t)

grows exponentially at first, then levels off as limiting factors kick in.

31
2.6.1 EXTRA: Existence and uniqueness of ODE solution
The existence theorem states that:

If the functions P and Q in a 1st order linear inhomogeneous ODE

are continuous over an interval t1 < t < t2 containing the point t = t0 ,
then there exists a unique solution y = f (t) which satisfies the equation,
and the initial condition y = y0 at t = t0 .

A proof of this statement can be found on P.13-18 of BD.

This statement is Good News. In classical dynamics, for a physical system
modelled by an ODE, we would not normally expect there to be discontinu-
ities as a function of t, so there should be a physically meaningful solution of
the ODE. The uniqueness of the solution comes from the specification of the
initial conditions. When you have found a solution that satisfies the ODE
and the initial conditions, then you are done. You do not need to look for
alternative solutions!

32
3 Momentum Conservation and Variable Mass
We have seen how N2 can be used to write down a differential equation for
the velocity and/or position of a system, which can then be used to obtain
the equation of motion. An alternative approach is to use N3, or equivalently,
the idea of Conservation of Momentum. A circumstance where this can be
particularly simplifying is when the mass of the system is non-constant. We
look at some examples of this in this section.

3.1 Variable Mass Problems

Up to now we have assumed that the inertial (or gravitational) mass of an
object is a constant, and that it is independent of time. If it is now, then N2
will have terms for both the acceleration and the rate of change of mass (in
vectorial form):
⃗ = d⃗p = M(t) d⃗v + ⃗v dM(t)
F
dt dt dt
The second term in this expression is non-zero if the total mass of a system
M(t) varies with time.
⃗ = 0, we have
For an isolated system with no external force applied, i.e. F
d⃗p
dt
= 0, hence momentum is conserved. However, the CoM of the system will
still accelerate or decelerate if the mass changes:

d⃗v ⃗v dM(t)
=−
dt M(t) dt

To solve a variable mass problem, you thus need to know the function M(t),
describing how the mass varies as a function of time.
The rest of this section will consider examples of variable mass problems.
First we will consider a case where F ⃗ = 0, and then we will include gravity
as an external force.

3.2 Rainfall into a moving train

Example 1: Rain on a train. A train moves along a flat, frictionless track.
Initially the train has mass M0 , with a constant velocity v0 . At time t = 0,
a heavy rain shower begins, and rainwater collects in an open-sided wagon
at a rate of µ. Assume the rain falls vertically, so it provides no horizontal

33
momentum to the system. What is the equation of motion for the velocity?

First, we will solve this using N2. The mass at any time t is

M(t) = M0 + µt.

Since there is no external force, N2 gives

dv v
=− µ
dt M0 + µt
where we have used dM/dt = µ, and v(t) is in the direction of the track.
This is separable, and after integration yields the equation

ln v = − ln(M0 + µt) + ln C

The initial condition v(t = 0) = v0 leads to C = M0 v0 . Hence we get

M0 v0
v(t) =
M0 + µt
So the train slows down with time, and as t → ∞, v → 1/t. This is effectively
because the additional rain must be accelerated up to the velocity of the train,
thus it acts to slow the train down.
Now let us solve this using N3. Equating the initial momentum M0 v0
with the momentum at some time t, we get

M0 v0 = (M0 + µt)v

which immediately yields the same result! This illustrates how conservation
of momentum can be a powerful tool for more easily solving problems in
certain situations.
Example 2: Rain on a train down the drain. Now instead of accumulat-
ing the rain, there is a drainage spout on the side of the wagon such that,
after some accumulation, the water drips out. Assume the drip-rate is also
µ, so there is no net accumulation of rain in the wagon (after the initial
accumulation), and the mass of the wagon is fixed at M0 . What is v(t)?
In this problem, we can still use conservation of momentum, but now we
have to account for the momentum of the material that went out the drain
so that we retain a closed system.

34
 Let us apply N3 over some small time interval t → t + dt, over which the
change in velocity is dv. The drained mass in dt is µdt, and thus the drained
momentum is µvdt. The momentum at t is M0 v, and the momentum at t+dt
is that of the train plus the drained mass, which must be equal:
M0 v = M0 (v + dv) + µvdt =⇒ M0 dv = −µvdt
This is a separable ODE:
dv µ µt
= − dt =⇒ ln(v) = − + ln(v0 )
v M0 M0
which gives
v(t) = v0 e−µt/M0
Hence the velocity of a train with rain down the drain on a plain drops
exponentially with time.
Remarkably, at large times, the exponential slowing is faster than the
no-drip case, which has v(t → ∞) ∝ 1/t. Perhaps counter-intuitively, it is
more efficient just to keep the water on the train than to let it drip out! The
reason is that, in the accumulating case, the train’s momentum increases over
time, so the momentum of the additional rain has fractionally less impact in
slowing down the train.

3.3 Raindrops in free fall

Raindrops condense out of clouds. As they fall under gravity, they pick up
additional mass from the surrounding water vapour. In detail this is a quite
complicated process, but we will assume a physicist’s raindrop: It is purely
spherical, uniform density, and has no drag forces (so there is no terminal
velocity from drag). We want to determine the velocity of the falling rain.
To solve this, we must specify m(t), the rate at which the raindrop ac-
cretes water. We look at two different models for m(t).
Example 1: First, we assume that the mass accretion rate dm/dt is pro-
portional to the volume (or equivalently, for a uniform density, mass) of the
drop:
dm
= αm(t)
dt
By integrating this we obtain the function m(t)
m(t) = m0 eαt

35
where m0 is the initial (seed) mass of the drop.
Here, we must include the effects of gravity. Because it’s difficult to
determine how much momentum gravity provides, it is simpler in this case
to go back to N2:
dvz dm(t)
−m(t)g = m(t) + vz
dt dt
Substituting dm/dt = αm(t), all the mass terms m(t) happily drop out, and
we are left with a standard linear ODE:
dvz
= −g − αvz
dt
The mass accretion term acts exactly like a drag force, reducing the accel-
eration, with α replacing the linear drag coefficient k/m. We have already
seen the solution of this equation:
g 
vz = − 1 − e−αt az = −ge−αt
α
As t → ∞ the terminal velocity is g/α, and the acceleration goes to zero.
Example 2: A more realistic model is a constant mass accretion rate:
dm
=k
dt
m(t) = m0 + kt
dvz
F = −m(t)g = m(t) + kvz
dt
The mass term m(t) no longer drops out and must be included explicitly:
dvz k
+ vz = −g
dt m0 + kt
This is an inhomogeneous 1st order ODE. No clever substitution is obvious, so
let us solve this using an integrating factor. We identify P (t) = k/(m0 + kt),
Q(t) = −g, thus the required integrating factor is:
Z 
µ(t) = exp P (t)dt = exp(ln(m0 + kt)) = m0 + kt

and the solution with the initial condition v = 0 at t = 0 is:

gt 1
 
vz (t) = − m0 + kt
m0 + kt 2
At small t, vz ≈ −gt as if no drag, but as t → ∞, it asymptotes to − 12 gt.

36
3.4 Rockets
A classic variable-mass scenario is a rocket burning fuel, which leads to the
famous rocket equation. Here we will derive this.
A rocket consists of a fuselage of mass M0 , and a fuel mass of m(t) which
is initially m0 . The total mass is thus M (t) = M0 + m(t). We will work in
the frame of the ground, so it is inertial and we can use Newton’s Laws. At
launch the rocket is stationary so the initial momentum is zero,
During the burn, fuel is ejected with constant velocity u relative to the
rocket, and thus lost at a rate dm/dt. Since the fuselage mass M0 is constant
M also goes down by this amount, so dM/ dt = dm/dt.
We want to determine the rocket’s velocity v(t). In the inertial frame
of the ground, the exit velocity of the fuel is v − u the rate of change of
momentum carried off by the fuel is (v − u)(dm/dt) = (v − u)(dM/dt). This
is thus the external force on the rocket. We assume this force is so large that
gravity is negligible (we’ll add gravity later). Thus N2 becomes
dv dM (t) dM
M (t) +v = (v − u) .
dt dt dt
Cancelling out the v dM
dt
term yields a separable ODE:
dv dM (t) dM
M (t) = −u =⇒ dv = −u
dt dt M
Regardless of the exact form of M (t) (but assuming u is constant), we can
integrate this to give
M (t = 0)
v(t) = −u ln(M ) + C =⇒ v(t) = u ln
M (t)
This is the so-called rocket equation, from Tsiolkovsky (1903).
Let us compute the velocity boost given to the rocket by time T after all
the fuel has been burned. At this time, M (T ) = M0 (since the fuel is spent):
m0
 
v(T ) = u ln 1 +
M0
The boost to the rocket v(T ) depends only on the speed of the fuel ejection, u,
and the fuel-to-fuselage mass ratio m0 /M0 . Remarkably, it does not depend
on the functional form of M (t) or v(t). This is why rockets are designed to
eject fuel at as high a velocity as possible.

37
 Perhaps counter-intuitively, it is possible to have the final velocity exceed
u: this requires (M0 + m0 )/M0 > e. This means it is possible to achieve
Earth’s escape velocity of 11.2 km/s, even though the exhaust is slower, if
one has substantially more fuel mass than rocket mass. Figure 5 shows the
final rocket velocity as a function of the fuel-to-rocket mass ratio, for various
values of the ejection velocity.

Figure 5: Solutions to the rocket equation, showing final velocity as a function

of the mass ratio for various values of the ejection velocity u (called ve in the
legend).

The rocket equation can be integrated to give the rocket’s height z(t):
" #
Z Z
(M0 + m0 )
z(t) = v(t)dt = u ln dt
M (t)
Note that the post-burn height z(T ) depends on the functional form of M (t),
even though the post-burn velocity doesn’t.
We didn’t consider gravity so far, but it is easy to do so, as long as we
assume the acceleration of gravity is a constant. This just adds a term to
the RHS in the rocket equation derivation:
dv dM (t)
M (t) =− u − M (t)g
dt dt
38
The first term on the RHS, known as the thrust, has to have a minimum
value of M (0)g to overcome gravity, otherwise the rocket stays on the ground
(remember, u is negative). Assuming sufficient thrust, the solution to this
ODE is much like the previous case, except with an extra term for gravity:
" #
dM (M0 + m0 )
dv = −u − gdt =⇒ v(t) = u ln − gt
M M (t)

More realistically, g will vary with time as the rocket moves away from the
Earth. This requires numerical integration to solve.
The final velocity at burnout time T now depends on T , owing to the last
term, not just u and the mass ratio. This then depends on the rate at which
the fuel is burnt. The faster it is used, the smaller T is, and the greater is
the final velocity. Thus it is advantageous to burn fuel quickly and eject it
as fast as possible. This is the engineering challenge of designing rockets.

3.5 EXTRA: Hot air balloon

As another example of a variable mass problem, we consider a hot air balloon.
A hot air balloon stays at a constant altitude if the force of gravity, −M g,
is balanced by the uplift of the air R – this is Archimedes’ Principle. We
assume that R is constant, and neglect drag forces.
The height of the balloon can be changed by throwing sand overboard.
The balloon contains a sandbag of initial mass m0 , and the mass of the
balloon without the sandbag is M0 . Initially the balloon is in equilibrium
with no velocity in the vertical direction, so R = (M0 + m0 )g and vz (0) = 0.
At time t = 0 the sand is trickled out of the sandbag at a constant rate
α. The relative velocity between the sand and the balloon is assumed to be
zero at the point of release u = 0. The amount of sand in the balloon as a
function of time is:
t m0
 
m = m0 − αt m = m0 1− α= = k(M0 + m0 )
T T
where T is the total time to trickle all the sand out of the sandbag, and we
have defined a new constant k for convenience.
The mass of the balloon + remaining sand as a function of time is:

M (t) = M0 + m0 − αt =⇒ M (t) = (M0 + m0 )(1 − kt)

39
The net upward external force on the balloon is given by:

Fext = R − M (t)g = [M0 + m0 − M (t)]g = (M0 + m0 )gkt

With u = 0 and dm = dM the momentum change of the balloon is:

dM (t) dvz dm dvz

Fext = vz + M (t) − vz = M (t)
dt dt dt dt
which simplifies to:
dvz gkt g
= =⇒ vz = −gt − ln(1 − kt)
dt (1 − kt) k

where we have used the initial condition vz = 0 at t = 0.

For small kt << 1 this gives:

k 2 t2 gkt2
" #
g
vz = −gt − −kt − =
k 2 2

and the initial acceleration can be seen to be consistent with the expression
for Fext :
az (kt << 1) = gkt

3.6 EXTRA: Motion of Centre of Mass

We now take a break from dynamics problems to apply conservation of mo-
mentum to a more abstract situation. This will allow us to introduce the
idea of a centre of mass.
We take for granted that we can push on an object at a single point (e.g.
a cue hitting a billiard ball), and the entire object will move in unison in
response. But why is this the case? How does the object know to distribute
the force received at one spot over the entire body?
To see how this works, we begin by generalizing N3 for an ensemble of
objects, to show that total momentum is conserved. We further show that
if there is a net external force acting on an object, then the entire object
obeys N2, and the internal forces in between the molecules of the object can
be neglected (i.e. they formally sum to zero). This proves that we can treat
solid bodies as single objects when it comes to bulk motion.

40
 ⃗ is defined as:
The Centre of Mass (CoM) of an N-body system R

⃗ = Pi mi⃗ri = 1
P
X
R mi⃗ri
i mi M i
P
where M= i mi is the total mass of the system. Note that with the excep-
tion of the total mass, all the other sums over the N-bodies in this section are
vector sums, i.e. three separate sums over the three components in space.
⃗ w.r.t. time gives the total velocity of the system:
Differentiating R

⃗ = Pi mi⃗vi = 1
P
X d⃗ri
V mi
i mi M i dt

This can be interpreted as the motion of the CoM:

⃗ ⃗
⃗ = dR = P
V
dt M
⃗ is
where the total momentum of the system P
⃗ =
X X
P ⃗pi = mi⃗vi
i i

We can write a version of N2 for each of the N-bodies, by dividing the force
⃗ i , and mutual forces between the bodies F
up into external components F ⃗ ij :

⃗i + ⃗ ij = mi d⃗vi
X
F F
j̸=i dt

When we sum this over all the N-bodies it becomes:

X
⃗i =
X d⃗vi
F mi
i i dt

where the contributions from the mutual forces have cancelled out in the sum
as a consequence of N3. This leads to the generalized result for N2 for an
N-body system:
⃗ ⃗ 2⃗
⃗ ext = ⃗ i = dP = M dV = M d R
X
F F
i dt dt dt2

41
For reference we give the simplified forms of these results for a two-body
system:
M = m1 + m2 ⃗ = 1 (m1⃗r1 + m2⃗r2 )
R
M
!
⃗ = 1 d⃗r1 d⃗r2 1
V m1 + m2 = (⃗p1 + ⃗p2 )
M dt dt M
⃗
⃗ = ⃗p1 + ⃗p2
P ⃗ ext = F
F ⃗1 +F⃗ 2 = dP
dt
We also give the integral forms of these results. These are relevant when
considering the motion of bodies made up of infinitesimal volume elements
dV , each of mass dmi = ρdV , where the density ρ may be uniform or a
function of ⃗r. Z
⃗ 1 Z
M = ρdV R= ρ⃗rdV
M
⃗ 1 Z d⃗r
V= ρ dV
M dt
⃗ = d⃗p = ρ d⃗r dV
Z Z
P
dt
⃗
⃗ = dP
Z
⃗ ext = dF
F
dt
Finally we restate the law of conservation of linear momentum for an
N-body system.

• If F
⃗ ext = 0 the system is said to be isolated. All three of its components
must be zero, and hence all three components of the total momentum
⃗ are separately conserved.
P

• If F
⃗ ext is non-zero it points in a particular direction. The component
of the total linear momentum P ⃗ in this direction changes according to
the generalized form of N2. The two components of P ⃗ perpendicular
⃗ ext are conserved.
to F

42
4 Conservation of Energy
We have seen that conservation of momentum can be a powerful tool to solve
dynamics problems. Likewise, conservation of energy can be similarly power-
ful, under the appropriate circumstances. In this section, we will define work,
kinetic and potential energy, scalar potentials, and how these all combine to
give energy conservation. Then we will show examples of applying energy
conservation to solve dynamical systems.

4.1 Impulse, Work and Power

From N2, we see that the overall change in momentum of a body due to a
net force acting on it can be obtained by integrating the force w.r.t. time:

∆⃗p = ⃗p2 − ⃗p1 =

R t2
⃗
F(t)dt
t1

The RHS of this equation is known as an impulse. It is a vector quantity

because it is the integral of a force vector over a finite time interval. If the
force varies in direction the integral must be performed separately for the
three components.
Work is defined as force times distance. For an arbitrary path from point
A to B, the work done is

W =
RB
F(⃗ ⃗
⃗ r).dl
A

This is a scalar quantity, because it is an integral over a scalar product. It

has units of energy: Joules (J) = kg m2 s−2 .
The work done (W ) can be positive or negative, depending on whether
the force is parallel or anti-parallel to the motion. If the force is perpendicular
to the motion, as in a centripetal force causing circular motion, then no work
is done. The work in general depends on the path taken between A and B.
The path length dl ⃗ corresponds to the motion owing to the velocity dv ⃗
over a time dt, such that dl⃗ = ⃗vdt. We can thus write the work as an integral
over time: Z tB
W = ⃗ v)dt
(F.⃗
tA

The scalar quantity:

⃗ v = dW/dt
P = F.⃗

43
is the rate at which work is being done, which is known as power. Again
it can be either positive or negative depending on the direction of the force
relative to the motion. It has units of Watts (W) = kg m2 s−3 , name after
James Watt, the Scottish inventor of the steam engine.

4.2 Kinetic Energy

We can use N2 to replace F ⃗ with md⃗v/dt, to get an integral over d⃗v rather
than dt: Z tB Z vB
d⃗v ⃗
W = m .⃗vdt = m⃗[Link]
tA dt vA

1
W = m(|⃗vB |2 − |⃗vA |2 )
2
If we define the kinetic energy as:

T = 12 m|⃗v|2 ,

then the work done by a force can be interpreted as a change in the kinetic
energy, T , of an object:

W = ∆T = TB − TA .

Below we consider some simple examples.

Free fall: For the case of free fall under gravity we have:
1
v = −gt =⇒ h(t) = h(0) − gt2
2
1 1
T = m|⃗v|2 = mg 2 t2 = mg[h(0) − h(t)]
2 2
The RHS is the integral of the force over the distance travelled, i.e. it is
the work done. You might also recognize this as the change in gravitational
potential energy mg∆h!

Circular motion: The centripetal force does no work, because F ⃗ is perpen-

⃗
dicular to dl and ⃗v. Hence there the kinetic energy is constant, with a value:
1 1
T = m|⃗v|2 = mr2 ω 2
2 2

44
4.3 Scalar Potential
The concept of a potential is extremely useful in physics. In many circum-
stances, a force can be represented as the spatial derivative of a scalar po-
tential, V (⃗r), which is only a function of position in space:

⃗ r) = − ∂V î − ∂V ĵ − ∂V k̂
F(⃗
∂x ∂y ∂z

It is not always true that V can be found for any F. ⃗ However, if it can, we
will see that this implies the force has special conservative properties, and V
becomes a valuable tool for solving dynamics problems.
Employing the ∇ notation from vector calculus, the force is the gradient
of the potential:

⃗ r) = −∇V (⃗r)
F(⃗

From a comparison with the integral definition of work done it can be seen
that the potential has units of energy. It can be either positive or negative,
and any constant offset V = V0 is removed when the derivative is taken.
The choice of V0 is an example of a gauge, which is a mathematical choice
that does not affect the physical behaviour of the system. For inverse square
forces, it is conventional to set V (∞) = 0. But any additive V0 can be used
that is convenient for the particular problem.
Note the minus sign in the definition of V ! This is purely convention, to
keep with our intuitive notion of potential energy and gravity – the gradient
represents the steepest direction up a hill, while the force of gravity pulls on
in the steepest direction down a hill. The potential thus defines a field in
space, analogous to rolling hills, where at any given location one feels a force
that depends on the local gradient of that potential (i.e. the steepness of the
hill where you’re standing).
An example of a scalar potential is a gravitational field:
GM m ⃗ r) = − GM m r̂
V (⃗r) = − F(⃗
|⃗r| |⃗r|2
where the negative sign indicates an attractive force.
Note that a −1/r potential is attractive, and a +1/r potential is repul-
sive. The latter would be like the electrostatic potential between like-sign
charges. Gravitational and electrostatic forces are examples of central forces

45
that have a 1/r2 dependence, and no dependence on time or the other spa-
tial coordinates. All central forces can be represented as potentials, which is
partly why potentials are useful in so many circumstances.

4.4 Example: A Perturbed Harmonic Potential

4
V(x)
3 F(x)

2
1
V(x),F(x)

0
1
2
3
41.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0
x

Figure 6: Potential V (x) (blue curve) and force F (x) (green) as described
in text. The unstable equilibrium points at x ± 1 are locations where V (x)
is maximized and F (x) = 0, while the stable equilibrium point is at x = 0,
where the force restores small displacements in either direction towards the
origin.

As an example of how to work with a potential, consider a 1-D V (x):

x2
!
2
V (x) = kx 1− .
2

The force is obtained by differentiation w.r.t x (remembering the − sign):

dV
F (x) = − = −2kx + 2kx3
dx
It is often illustrative to sketch potentials, as shown for this case in Figure 6
above, because one can think of them much like hills. For this potential, the
main features of a sketch are:

46
 • There are three equilibrium points where F = 0 at x = ±1 and x = 0.
These correspond to extrema in the potential curve, i.e. peaks and
valleys. Peaks are unstable equilibria, while valleys are stable.

• For small |x| ≪ 1 the x3 term can be neglected and F ∝ −x. This
is a restoring force that moves objects back to the origin. You might
recognise this as a simple harmonic oscillator (SHO). From this it fol-
lows that x = 0 is a point of stable equilibrium, with a local minimum
at V (0) = 0

• For x < −1 the force is negative, while for x > +1 it is positive. From
this it follows that the points x = ±1 are points of unstable equilibrium
with maxima at V (±1) = k/2.

• As x → ±∞ the higher power terms in x dominate. The potential goes

to −∞ at both limits, and the force goes to ±∞.

• The turning points where F (x) changes sign are obtained from dF/dx =
0:
1 4k
x = ±√ =⇒ F = ∓ √
3 3 3

4.5 Conservative Forces

The work done by a force along a path from A to B can be expressed as a
potential difference:
Z B Z B
W = F. ⃗ =−
⃗ dl ⃗
∇[Link]
A A

⃗ is simpler than it looks. First, note that dl

The term ∇[Link] ⃗ = dxî+dyĵ+dzk̂.
Taking the dot product with ∇, we get that the integrand is ∂V∂x
dx + ∂V
∂y
dy +
∂V
∂z
dz. The partial derivatives sum up to simply give dV , leaving an integral
over dV from A to B. Thus

W = VA − VB = −∆V

This holds if V (⃗r) is a scalar function of position with a unique value at any
point in space. It follows that for any path between the points A and B, the

47
Figure 7: Explanation of a conservative force, as being one where the work
done around any closed path is zero.

RHS gives the same change in the potential energy. If we consider a closed
loop on which A and B are any two points:
I Z B Z A
W = F. ⃗ =−
⃗ dl ⃗ −
∇[Link] ⃗ =0
∇[Link]
A B

This is zero because the change in V between points A and B is independent

of the path taken. This leads to the definition of a conservative force which
must satisfy two conditions:
I. The force F (⃗r) must only be a function of position in space.
It cannot depend on time t or on the velocity ⃗v of an object.
II. The integral of the force around a closed loop must be zero.
I
⃗ =0
⃗ dl
F.

It turns out, both of these conditions are satisfied if we can write:

⃗ r) = −∇V (⃗r)
F(⃗
This can be seen by noting that, in vector calculus, the curl of a gradient is
always zero. One can think of the gradient as the direction of steepest ascent,
while the curl picks out the direction around the hill with no elevation change;
these two are orthogonal, and hence applying both operators in succession
always yields zero. Hence an alternative definition of a conservative force is:

48
 ⃗ =0
∇×F
where we have used the curl operator ∇×, which is given by
! ! !
⃗ = ∂Fy ∂Fz ∂Fz ∂Fx ∂Fx ∂Fy
∇×F − î + − ĵ + − k̂
∂z ∂y ∂x ∂z ∂y ∂x

To prove that this corresponds to the same definition as (II) above, we

note another result from vector calculus: Stokes’s Law. This relates a line
integral over a closed loop and a surface integral over the area enclosed by
⃗
the loop, for any vector field F:
I Z
F. ⃗ =
⃗ dl ⃗ dS
(∇ × F). ⃗
S

If one can write F⃗ = −∇V , then the RHS integrand (∇ × ∇)V = 0, which
immediately proves that the LHS (corresponding to condition II above) is
also identically zero, therefore the force is conservative.

4.6 Energy Conservation

In the previous sections, we have found that W = ∆T and W = −∆V .
Equating the two, we arrive at the law of conservation of energy:

∆T + ∆V = 0 E = T + V = constant

The sum of the kinetic energy and the potential energy

remains constant in the presence of a conservative force.
This statement can be generalised to several forces:
⃗ i acting on an object are conservative, each of which
If all the forces F
can be expressed in terms of a corresponding potential energy Ui (⃗r),
then the total energy is a constant.

X
E=T+ Ui = constant
i

The potential is a scalar quantity, and the sum of the individual potentials Ui
gives the overall potential V (⃗r). The scalar addition is a convenient feature
of potentials that simplifies the mathematics.

49
 ⃗ r) can then be obtained from the gradient of V :
The net external force F(⃗

⃗ = ⃗ i = −∇V
X X
V = Ui F F
i i

Conservation of energy only works for conservative forces!

Finally, energy conservation can be generalised to a system of many bod-

ies:
⃗ i acting on an N body system,
If both the external forces F
and the mutual forces F⃗ ij acting between them, are conservative,
then the total energy of the N body system is a constant.

X X XX
E= Ti + Ui + Uij = constant.
i i i j>i

The first term represents the total kinetic energy of all objects, the second
term represents the external potential that the objects fell, and the third
term represents the interactions between the objects that gives rise to a
(conservative) force. Note that the condition j > i is required to avoid
double counting the potential energy associated with the interaction forces
between the bodies.
If the total energy E is a constant, then any derivative of E must be
zero. A specific example is dE/dt = 0, but the derivative can be taken w.r.t
any independent variable. This can be used to check if you are dealing with
conservative forces.

Energy conservation can be a powerful tool for solving dynamics problems.

An example is the determination of the Earth’s escape velocity. One could
solve the dynamics problem of an object with an initial velocity, under a
height-varying gravitational force, which is long and difficult. Alternatively,
with conservation of energy it is almost trivial.
Assume an object is launched with velocity v0 . The initial kinetic energy
is T (0) = 21 mv02 , while the initial potential energy is V (0) = −GME m/RE
(in the gauge where V (∞) = 0), where ME and RE are the mass and radius
of the Earth. To just barely escape, we want the kinetic energy to be 0 at
infinity, where by definition V = 0 and thus the final total energy is zero.

50
Equating the initial and final energy gives:
1 2 GME m
mv − =0
2 0 RE
2GME q
v02 = =⇒ v0 = 2gRE ≈ 11.2 km/s
RE
Note that the escape velocity v0 at the Earth’s surface can be in any direc-
tion (if it doesn’t hit the Earth, obviously). It does not have to be vertically
upwards. This is because T is a scalar, and thus is does not matter which
direction v0 is pointed.

Conservation laws are often useful if you know the system state at an ini-
tial time, and are only interested in the system state at some final time,
without wanting to know the intermediate details (i.e. the full position or
velocity evolution). Working with energies can then be simpler because they
are scalars that sum, unlike dealing with vector forces.

4.7 EXTRA: Non-conservative Forces

We give a few examples of non-conservative forces.
• If a force is time-dependent then the integral of the force round a closed
loop is no longer zero.
I I  
⃗
F(t). ⃗ =
dl ⃗ v dt ̸= 0
F(t).⃗

It may still be possible to write the Force as a gradient of a time-

dependent potential:
⃗ r, t) = −∇V (⃗r, t)
F(⃗
but when you go round the loop from A → B and then back from
B → A, the potential at A will change because of the time dependence.
• If a force is velocity-dependent then the work done depends on the
speed with which the loop is travelled around. The example we have
seen of this type of force is the drag force due to air resistance. This
⃗ v is always negative,
always acts against the direction of motion, so F.⃗
⃗
and the integral of F(t). ⃗ round a closed loop is no longer zero. Work
dl
has to be done to overcome the drag force at all points round the loop,
and this work can not then be used to change the Kinetic Energy.

51
 • Another familiar example of a non-conservative force is dynamical fric-
⃗ v is negative.
tion. This always acts against the motion, and again F.⃗
For these types of forces energy conservation can be restored if we in-
troduce a transfer of mechanical energy into thermal energy. As the
object slows down the friction heats up the surfaces.
As an example consider a block of mass m sliding down an incline at an angle
θ against a frictional force µN . Taking components of the forces normal and
parallel to the slope:
dv
N = mg cos θ m = mg(sin θ − µ cos θ)
dt
v = gt(sin θ − µ cos θ)
and the distance travelled down the slope is:
gt2
d= (sin θ − µ cos θ)
2
The kinetic energy is:
mg 2 t2
T = (sin θ − µ cos θ)2
2
and the potential energy changes by:
mg 2 t2
V = mgd sin θ = (sin2 θ − µ sin θ cos θ)
2
The kinetic and potential energy changes are only equal if µ = 0. The
difference between them:
µmg 2 t2
∆E = T − V = − (sin θ − µ cos θ) cos θ = −µmgd cos θ
2
This can be interpreted as the work done against the frictional force:
Z
W = F. ⃗ = −µN d = −µmgd cos θ
⃗ dl

Note that, had we also included the energy associated with motion of
the atoms owing to the frictional force (i.e. the sound or heat energy gener-
ated), then the system would again be conservative. It is essentially because
we are transferring energy out of our defined system that we obtain non-
conservation.

END OF PART I

52
5 Simple Harmonic Motion
Oscillatory motion occurs in a wide range of physical circumstances. Gener-
ically, it results from any circumstance where displacement from an equilib-
rium value results in a restoring force. In this part of the course we will
mainly discuss small or linear oscillations, which can be described by simple
harmonic motion (SHM).
We will develop the theory of free, damped and forced oscillations of
mechanical systems. In the last part of this section we will discuss coupled
oscillations of systems with more than one degree of freedom, and use matrix
algebra to solve the coupled equations. This will lead to the concept of normal
modes.

5.1 Equilibrium and Small Oscillations

Consider an object at rest, in equilibrium, with no net force acting on it.
Restricting ourselves to conservative forces, we can define the force in terms
of a scalar potential:
⃗ = −∇V (⃗r) = 0
F
then this condition is satisfied if V is either a maximum (unstable), or a
minimum (stable). If the object is displaced from the equilibrium point by
some small distance x (in 1-D), then the potential can be expressed as a
Taylor expansion:
dV (0) x2 d2 V (0)
V (x) = V (0) + x + + ...
dx 2 dx2
Since the normalisation of V is arbitrary, we can freely choose V (0) = 0.
Also, at x = 0 our 1-D system satisfies F (x = 0) = −dV /dx = 0, since x = 0
is an equilibrium point. So the first non-vanishing term in the expansion is:
1
V (x) ≈ kx2
2
where k is now identified as the second derivative of the potential at the
equilibrium point x = 0. If k is positive, then the potential forms a parabolic
“cup”, which means any displacement will result in a restoring force back
towards x = 0. This is seen by noting that the force as a function of the
(small) displacement x is
dV
Fx = − = −kx
dx
53
Such a force gives rise to simple harmonic motion (SHM), and the quadratic
form of V (x) is known as a harmonic potential.
The equation F = −kx may be familiar as Hooke’s law for a spring, where
k is known as the spring constant. Another familiar system that performs
SHM is the simple pendulum, where we will later show that k = mg/L.
Note that the linear force relationship is only true in the limit of small
x, where the higher order terms in the Taylor expansion of V (x) can be
neglected. The higher order terms can be regarded as perturbations of SHM.

5.2 Description of SHM

SHM is described by

d2 x d2 x
Fx = m 2
= −kx =⇒ 2
+ ω2x = 0
dt dt
where we have defined the natural frequency of the system defined as
s
k
ω≡
m
From a dimensional analysis one can see that this will have the units of
inverse time, which is why we call it a frequency.
To solve this ODE, we need a function whose second time derivative is
itself (modulo constant factors) – only in this case will we be able to divide
out the time dependence from the system to satisfy the LHS being zero at
all times. Examples of such a function are an exponential, sine, and cosine;
indeed, all are valid solutions.
Let us first try x = eλt (we will worry about constants of integration
later). Here, d2 x/dt2 = λ2 eλt . Hence we can cancel the time dependent eλt
terms to get
λ2 + ω 2 = 0 =⇒ λ = ±iω
There are thus two possible functional forms that satisfy this, e+iωt and e−iωt .
Since this is a linear ODE, the most general solution is a linear superposi-
tion of the two possible solutions. This yields the exponential form for the
SHM solution. Using Euler’s formula and trigonometric relations, one can
straightforwardly show that there are three mathematically equivalent forms:

54
 • Exponential form:
x(t) = C1 e+iωt + C2 e−iωt
where C1 and C2 are complex conjugate coefficients that have to satisfy
C1 = C2∗ in order to make x real.

• Sine/Cosine form:

x(t) = A1 sin ωt + A2 cos ωt

where A1 and A2 are real: A1 = i(C1 − C2 ) and A2 = (C1 + C2 ).

• Amplitude/Phase form:

x(t) = A sin(ωt + δ)
q
where the amplitude A = A21 + A22 and the phase tan δ = A2 /A1 .

Since this is a 2nd order ODE, as expected there are two constants of in-
tegration. The values of these constants (C1 , C2 ) or (A1 , A2 ) or (A, δ) are
set by the boundary conditions. Which form you choose to use is typically
determined by how the boundary conditions are specified.

Let’s consider the amplitude/phase form. Here, the velocity of the motion is
given by:
dx
vx = = Aω cos(ωt + δ)
dt
This is π/2 out of phase with the displacement x. The velocity is zero at the
maximum displacements x = ±A, and has a maximum value vx = ±Aω at
x = 0.
If we apply the initial condition that the object starts from its equilibrium
point x(0) = 0, then we immediately see that δ = 0, and the solution can be
written as:
x(t) = A sin ωt
If we instead apply the initial condition that the object starts at its maximum
amplitude with vx (0) = 0, then δ = π/2, and the solution can be written as:

x(t) = A cos ωt

55
 In the Sine/Cosine form, for any initial conditions x(0) and vx (0), the
coefficients can be identified as:
vx (0)
A2 = x(0) A1 =
ω
Next we consider some specific examples of SHM.

5.2.1 A Vertical Spring

A spring hangs from a ceiling with a mass m attached to it. If z0 is the
vertical extension of the spring at equilibrium, then gravity must balance
the restoring force from the spring to leave no net force:

mg = kz0

We extend the spring downwards by some small displacement z from its

equilibrium position. Note that positive is downwards here, so the force of
gravity is mg, giving an equation of motion:

d2 z
Fz = m 2 = mg − k(z0 + z)
dt
Using the static equilibrium equation this can be simplified to:

d2 z
Fz = m = −kz
dt2
We can immediatelyq identify the solution to this as a harmonic function with
a frequency ω = k/m, as above. Hence the mass performs SHM about the
equilibrium point with a frequency that depends only on the spring constant
k and the mass m. The effect of gravity g and the static extension of the
spring z0 do not appear in the equation of motion.
This illustrates a key point in SHM: the only thing you need to consider in
small oscillations is the change in the net force compared to the equilibrium
position; any static equilibrium force is irrelevant!

5.2.2 The Simple Pendulum

Imagine a mass (bob) hanging from a rod of length L tied to the ceiling,
and it is swinging back and forth. Polar coordinates are useful here, since

56
Figure 8: Set up of simple pendulum, with a rod of length L. If θ is small, the
solution θ(t) is a simple harmonic oscillator,
q identical to the spring solution
above with a characteristic frequency ω = g/L.

if we define the origin to the be the suspension point, the radius of the bob
is constant at L. The position and velocity are thus (recalling your polar
coordinate derivatives from §1.3):
d⃗r dr̂ dθ
⃗r = Lr̂ =⇒ ⃗v = = L = L θ̂
dt dt dt
We take the time derivative of this to get the acceleration, which has com-
ponents in both the r̂ and θ̂ directions:
!2
d2 θ dθ dθ̂ d2 θ dθ
⃗a = L 2 θ̂ + L = L 2 θ̂ − L r̂
dt dt dt dt dt
We can set m⃗a equal to the force to obtain the equation of motion. To do
this, we need to draw a force diagram for the bob at an angle θ. We will
⃗ along the rod is in the
choose to work in (r̂, θ̂) coordinates. The tension T
−r̂ direction, but gravity will have both a radial and a theta component. A
bit of geometry shows that we can write a θ̂ force and a r̂ force as follows:
Fθ = −mg sin θ Fr = mg cos θ − T
⃗ Thus we have two force equations, corresponding to the θ̂
where T ≡ |T|.
and r̂ components of the above acceleration:
!2
d2 θ dθ
mL 2 = −mg sin θ, mL = mg cos θ − T
dt dt

57
The second of these might be familiar. Recall the centripetal force for circular
motion is mrω 2 = mL( dθdt
)2 in our case. Thus the r̂ equation is simply saying
that the inward centripetal force must be balanced by the outward radial
component of gravity minus the tension, in order for it to move along a
circular path.
The θ̂ equation then gives us the angular motion of the bob. If we assume
small perturbations around equilibrium, then in the small angle approxima-
tion sin θ ≈ θ, the tangential force is proportional to θ. This leads to a SHM
equation, like we had before:
d2 θ d2 θ g
−mL = −mgθ =⇒ = θ
dt2 dt2 L
q
The natural frequency of the oscillation is read off as ω0 = g/L.
Using the amplitude/phase form and assuming an angular amplitude Aθ ,
the motion is described by:

θ = Aθ sin(ω0 t + δ)
dθ
v=L = Aθ Lω0 cos(ω0 t + δ).
dt
Note that, unlike in the uniform circular motion case we considered earlier,
here ω = dθ/dt is time dependent. Thus the tension in the rod also must be
time dependent in order to balance the time-variable centripetal force.
We can use the radial (r̂) force equation to compute the tension T (t). We
solve for T , employing the small-angle approximation cos θ ≈ 1 − θ2 /2:
θ2 mv 2
!
T = mg cos θ − Fr = mg 1 − + ,
2 L
where we have written the centripetal force in the v form. As one expects
intuitively, the tension is minimized when v = 0 at the turnaround points
where θ = ±Aθ , with a value T = mg(1 − A2θ /2). One can show using the
formula for v above that it is maximized when v is largest at θ = 0, with
T = mg(1 + A2θ ).

5.3 Energy in Simple Harmonic Motion

Now that we have computed the velocity, we can determine the energetics
of SHM. Going back to our original SHM described by V (x) = kx2 /2, and

58
Figure 9: Diagram showing how potential and kinetic energy vary as a func-
tion of the displacement a. At the origin, the object is moving its fastest,
so K is maximised, while at the ends the energy is entirely in V . The total
energy is constant throughout.

x(t) = A sin(ω0 t + δ), the kinetic energy K is given by:

1 1 1
K = mv 2 = mA2 ω02 cos2 (ω0 t + δ) = kA2 cos2 (ω0 t + δ)
2 2 2
where we have used the amplitude/phase form and k = mω02 . Meanwhile the
potential energy is
1 1
V = kx2 = kA2 sin2 (ω0 t + δ)
2 2
At x = 0 the kinetic energy is a maximum kA2 /2 and the potential energy
is zero. At x = ±A the kinetic energy is zero and the potential energy is
kA2 /2. The two forms of energy are in antiphase with each other, as shown
in Figure 9 above. When averaged over a full oscillation there are equal
amounts of kinetic and potential energy, because < cos2 >=< sin2 >= 1/2.
The total energy is a constant, independent of time, since cos2 + sin2 = 1:
1 1
E = K + V = kA2 = mA2 ω02
2 2
Intuitively, this makes sense, as the total energy increases with the mass,
frequency, and amplitude of the oscillations.

59
5.3.1 EXTRA: Energetics of the simple pendulum
We can do an analogous exercise for the simple pendulum. For the simple
pendulum, the kinetic energy of the bob is:
!2
1 dθ 1
K = mL2 = mL2 A2θ ω02 cos2 (ω0 t + δ)
2 dt 2
For the potential energy, we need to know the height of the bob, and then we
case use mgh to describe the potential energy gained. Some simple geometry
shows that h = L − L cos θ above the bottom (θ = 0) point, hence
θ2 1
V = mg(L − L cos θ) ≈ mgL = mL2 A2θ ω02 sin2 (ω0 t + δ)
2 2
where the last step uses g = ω02 L. Hence the expression for the total energy
of small oscillations of a simple pendulum is:
1
E = K + V = mω02 L2 A2θ
2
which is the same expression as the SHM case with A = LAθ .

5.3.2 From energy conservation to an equation of motion

An interesting application of energy conservation is that one can actually use
it to derive the equation of motion. This might be useful in a situation where
one only knows the energetics of a system, but wants to determine the full
equation of motion.
For this exercise, let us pretend we don’t know the equation of motion,
but we know how to calculate the kinetic and potential energies. For the
pendulum, we can write !2
1 2 dθ
K = mL
2 dt
while a bit of geometry shows that the height of the bob is L − L cos θ, so
V = mg(L − L cos θ)
Energy conservation tells us that d(T + V )/dt should be zero:
 !2 
dE d 1 dθ
= mL2 + mg(L − L cos θ) = 0
dt dt 2 dt

60
 dE d2 θ dθ dθ
= 2 mL2 + mgL sin θ = 0
dt dt dt dt
d2 θ
L 2 + g sin θ = 0
dt
which is exactly the equation of motion (EoM) for the pendulum!
We can even use energetics to determine the period of oscillation. You
probably recall that the period of oscillation is given by 2π/ω0 , but let us
derive this.
Using energy conservation the kinetic energy can be written in terms of
the potential V (x):
1
K = mvx2 = E − V (x)
2
s
dx 2(E − V (x))
vx = =
dt m
This can be integrated over dx to give the time taken by the oscillations:
dx mZ dx
Z r
t= = q
vx 2 E − V (x)

The period of the oscillations is double the time for the mass to go from
x = −A to x = +A:
√ Z +A
dx
τ = 2m q
−A E − V (x)
So far this is completely general, for any oscillator. Now we specify a SHM,
for which V = kx2 /2 and E = kA2 /2:

m Z +A dx
r
τ =2 √
k −A A 2 − x2
This can be integrated using the substitution x = A sin θ, dx = A cos θdθ:

m Z π/2 m 2π
r r
τ =2 dθ = 2π =
k −π/2 k ω0
which is the expected result for SHM. This shows that we can use energy
conservation to obtain interesting information about the system, without
having to solve the equation of motion.

61
5.4 Damped Oscillations
We now consider a harmonic oscillator where there is a damping force, which
like the drag force depends on velocity. From the drag force discussion in
§1.2, we saw that drag forces are typically linearly or quadratically dependent
on the velocity. In the small-angle approximation, the quadratic terms are
expected to be small, hence we will only consider linear drag. This also makes
the solutions analytically tractable.
We introduce a drag force proportional to the velocity

d2 x dx
F =m 2
= −kx − b (1)
dt dt
where b is a constant coefficient setting the strength of damping.
Let us try an exponential solution, since any derivative of an exponential
is still an exponential. Trying

λt dx λt d2 x
x=e =⇒ = λe =⇒ 2
= λ2 eλt
dt dt
For convenience we define
b k
γ≡ ω02 ≡
2m m
which then results in a quadratic equation in λ (since the eλt cancels out of
each term):
λ2 + 2γλ + ω02 = 0
The solution is thus x(t) = eλt where λ can take two values:
q
λ = −γ ± γ 2 − ω02

As a result, there are two independent solutions, corresponding to the + and

− versions of λ. Since the ODE is linear, the fully general solution is a linear
combination of these two:
√ 2 2 √ 2 2
x(t) = C1 e(−γ+ γ −ω0 )t + C2 e(−γ− γ −ω0 )t
q
Defining ω ≡ ω02 − γ 2 , we get the solution
h i
x(t) = e−γt C1 eiωt + C2 e−iωt

62
This can also be written in amplitude/phase form:

x(t) = A0 e−γt sin(ωt + ϕ)

This is just like a SHO, except with a e−γt term which represents a decaying
amplitude owing to the damping constant γ.
This oscillator will have different behaviours depending on ω:

• Underdamped if ω is real (ω0 > γ).

• Overdamped if ω is imaginary (ω0 < γ).

• Critically damped if ω = 0 (ω0 = γ).

We consider each of these in turn.

5.4.1 Underdamped SHM

Figure 10: Evolution of underdamped oscillator. The amplitude of the oscil-

lations drops exponentially with time.
q The oscillation frequency is lowered
from the natural frequency by ω = ω02 − γ 2 .

If ω > 0, then we have an underdamped SHO. There are still oscillations

but its amplitude decreases with time.

63
 The oscillator frequency is ω, superimposed with an exponential decay
in amplitude from an initial value A0 , with a decay constant γ, as shown in
Figure 10. Note that the oscillation frequency ω is different than the
q natural
frequency ω0 ; it is lowered by the introduction of damping to ω = ω02 − γ 2 .
In the limit of very small damping γ ≪ ω0 there are many oscillations
with a frequency ω ≈ ω0 , and the amplitude decays away slowly. One can
show that the total energy of the oscillations in this case is:
1
E = mA20 ω02 e−2γt
2
Even with such minimal damping, energy is not conserved in a damped SHO.
It decays away with a decay constant 2γ, because work has to be done to
overcome the damping term.

5.4.2 Overdamping
The opposite case is when γ > ω0 and ω is imaginary, which results in
overdamping. Here we define
q
λ = −γ ± δ, δ≡ γ 2 − ω02 ,

where 0 < δ < γ represents the strength of the damping. The general solution
for overdamping is thus:

x = (C1 eδt + C2 e−δt )e−γt = C1 e−(γ−δ)t + C2 e−(γ+δ)t

Since γ > δ, both terms result in a decaying amplitude. The C2 term decays
faster than the C1 term.
Note that there are no oscillations! Overdamped systems only show tran-
sient responses to the initial conditions, so the motion approaches zero after
sufficient time.

5.4.3 Critical damping

Critical damping occurs when γ = ω0 . In this case there are repeated roots
with λ = −γ. This is a problem, since we now only have a single solution,
where we should have two. In the next section, we will show that the general
solution in this case is:
x = (C1 + C2 t)e−γt

64
Figure 11: Evolution of damped oscillators with various damping strengths.
The blue line shows the original undamped SHO. The red underdamped line
shows oscillations with a decaying amplitude. The green critically damped
case is shown here starting from rest which leads to exponential decay. The
magenta overdamped case approaches zero via a combination of two decaying
exponentials.

As with overdamping, there are no oscillations, and the motion is damped

away with a decay constant γ. However, in critical damping, x(t → ∞) →
te−γt . Unlike with overdamping where the oscillations were damped to zero
exponentially, here the decay is not generally exponential.
The detailed behaviour depends on the initial conditions, which specify
C1 and C2 . For zero initial velocity, the oscillations do decay exponentially:
dx
x(0) = x0 , (0) = 0 =⇒ x = x0 e−γt .
dt
However, with a non-zero initial velocity an object can move away from the
initial position before decaying back:
dx
x0 = 0 (0) = v0 x = v0 te−γt .
dt
At small t, x grows linearly, but eventually the exponential term takes over
and the object returns towards equilibrium.
The various cases of a damped oscillator are illustrated in Figure 11.

Finally, we would like to add a driving term to our damped SHO. However,
the increased complexity of the resulting ODEs requires us taking another
maths interlude to learn how to solve 2nd order ODEs more generally.

65
6 Second order ODEs
As you’ve seen, a simple harmonic oscillator (SHO) results in a 2nd order
ODE. Fortunately, this case had a rather obvious trial solution. The more
general case does not. In this section we develop some techniques to solve
2nd order ODEs. This is a rather lengthy subject worthy of a full course
in mathematics, so we will only cover some very basic cases applicable to
dynamics. For one, we will restrict ourselves to linear 2nd order ODEs.

6.1 Simple harmonic motion as a 2nd order ODE

We begin with the homogeneous case of a linear 2nd order ODE:

d2 y dy
2
+ p(t) + q(t)y = 0
dt dt
A SHO (including damping) corresponds to the case where p(t) and q(t) are
constants, p and q.
As we have seen, the solution is to try the form y = Aeλt and substitute
it into the 2nd order ODE. This leads to a polynomial in λ which has two
solutions:
1
 q 
λ2 + λp + q = 0 =⇒ λ = −p ± p2 − 4q
2
For constant p and q this gives two constant values of λ. and the solution is
a linear sum of the two cases:

yc (t) = C1 eλ1 t + C2 eλ2 t

If λ1 = λ2 , i.e. p2 = 4q, then the two terms are degenerate. This was the
case for the critically damped SHM, and the solution is to modify one of the
terms by multiplying it by t:

yc (t) = C1 eλt + C2 teλt

We will prove this in §6.4.

In general, even for a time-dependent p(t) and q(t), the solution to the
homogenous case will be a superposition of two linearly independent solutions:

yc (t) = C1 y1 (t) + C2 y2 (t)

66
 Each of y1 (t) and y2 (t) are independent solutions to the original ODE,
and can be simply added because the ODE is linear. We can prove this by
differentiating the above and substituting into the 2nd order ODE:
d2 y 1 d2 y 2
! !
dy1 dy2
C1 2 + C2 2 + p C1 + C2 + q (C1 y1 + C2 y2 ) = 0
dt dt dt dt

d2 y 1 d2 y2
! !
dy1 dy2
C1 + p + qy1 + C 2 + p + qy2 = 0
dt2 dt dt2 dt
which is satisfied if y1 and y2 are separately solutions of the original ODE.
C1 and C2 are the fractions of each separate solution. These coefficients
are determined by two initial conditions, which are usually on y(t0 ) and
dy/dt(t0 ), or less commonly on y(t1 ) and y(t2 ), or on dy/dt(t1 ) and dy/dt(t2 ).
EXTRA: One can check that the functions y1 and y2 are independent
by applying the initial conditions:
dy1 (t0 ) dy2 (t0 ) dy(t0 )
C1 y1 (t0 ) + C2 y2 (t0 ) = y(t0 ) C1 + C2 =
dt dt dt
These are a pair of simultaneous equation, which can be solved uniquely for
the unknowns C1 and C2 provided that the determinant of the coefficients
doesn’t vanish:
y1 (t0 ) y2 (t0 ) dy2 (t0 ) dy1 (t0 )
W= = y1 (t0 ) − y2 (t0 ) ̸= 0
dy1 (t0 )/dt dy2 (t0 )/dt dt dt
This determinant is known as the Wronskian, W, of the functions y1 and y2 .
If it is non-zero for any chosen t0 then y1 and y2 are a linearly independent
set of functions, and their superposition gives a unique solution for the pa-
rameters C1 and C2 for a given set of initial conditions. Examples of such
linearly independent pairs of functions are:

y1 = sin ωt y2 = cos ωt W = −ω

y1 = eiωt y2 = e−iωt W = −2iω

6.2 D’Alambert’s Method: Partially known solutions

In the general case where p(t) and q(t) are time-dependent, there is no sys-
tematic approach to obtaining a full solution. However, if one can somehow

67
guess the first solution y1 (t), then it is always possible to determine the
second solution y2 (t) using D’Alambert’s method.
We assume that we know a solution y1 (t) that satisfies

d2 y1 dy1
2
+ p(t) + q(t)y1 = 0
dt dt
and we want a second linearly-independent solution y2 . We begin by defining
a function u(t) such that

dy2 du dy1
y2 (t) = u(t)y1 (t) =⇒ = y1 +u
dt dt dt
d2 y 2 d2 u dy1 du d2 y 1
=⇒ = y 1 + 2 + u
dt2 dt2 dt dt dt2
Since y2 is also a solution to the original ODE, we can substitute in the
above terms to obtain
d2 u dy1 du d2 y1 du dy1
y1 2 + 2 + u 2 + p(t)y1 + p(t)u + q(t)uy1 = 0
dt dt dt dt dt dt
Rearranging terms, we have

d2 u d2 y1
!
dy1 du du dy1
y1 2 + 2 + p(t)y1 +u 2
+ p(t) + q(t)y1 = 0
dt dt dt dt dt dt

The magic of D’Alambert is that the term inside the parantheses is zero,
because it is exactly the original ODE for y1 (t)!
The remaining terms then represent a linear 1st order ODE in w(t) ≡
du/dt: !
dw dy1
y1 + p(t)y1 + 2 w=0
dt dt
Since y1 is known, this ODE can always be solved via techniques we developed
in §2. The solution yields w, which can be integrated to get u(t) = wdt.
R

Multiplying this by y1 (t) then gives the desired second solution y2 (t). Re-
markably, finding y2 doesn’t depend on knowing q(t) at all!
Thus D’Alambert’s method provides a guaranteed method to obtain a
second solution to a linear 2nd order homogeneous ODE if a first solution is
already known.

68
 Example: Consider the ODE
d2 y dy
t2 2
− 2t + (2 − t2 )y = 0
dt dt
for which you can verify that y1 = te−t is a solution by differentiation and
substitution. Identifying p(t) = −2/t and setting y2 (t) = u(t)te−t , we arrive
at a first order ODE in du/dt:

d2 u du
te−t 2
− 2te−t =0
dt dt
The te−t cancels out, which leaves a simple ODE that has the solution:
du C
= Ce2t =⇒ u = e2t
dt 2
If C2 = C/2, then the second solution is thus:

y2 (t) = C2 e2t te−t = C2 tet

Hence the final solution is

y(t) = C1 te−t + C2 tet

One might have guessed this second solution given the form of the first, but
in any case this shows how to get it using D’Alambert’s method.

6.3 The inhomogenous case: Particular functions

We finally consider the inhomogeneous case when f (t) ̸= 0, which we will
need when we consider driven harmonic oscillators. For simplicity, we return
to the assumption that p and q are constants, but the discussion here would
apply for time-dependent functions as well.
Our general inhomogeneous linear 2nd order ODE is

d2 y dy
2
+ p + qy = f (t)
dt dt
Let us break up the solution into two parts: y = yc + yp , where yc is the
complementary solution, while yp is the particular solution – the reason for
these names will become evident.

69
 Since this is a linear ODE, we can consider yc to be the solution to the case
where f (t) = 0, i.e. the homogenous version of this ODE. This must be part
of the solution to this ODE, because f (t) = 0 is a perfectly valid choice. Since
p and q are assumed to be constant, we have already determined yc = eλt ,
where λ is given by the solutions to the quadratic equation in §6.1. If p
and q are not constant, one must guess a solution and/or use D’Alambert’s
method.
Now, we just need to find the particular solution yp , for which the RHS is
f (t). Naturally, this is the difficult bit. The function yp must satisfy the ODE
with f (t) included, and must be linearly independent of the complementary
function yc . There is no generally applicable method for finding yp . But
the nice part is, it turns out we now only need to find one solution; the
complementary solution has already given us two other solutions!
In many cases, the functional form of f (t) can give us a big hint. Here
are some common cases:

• If f (t) = aekt , try the form yp = Ae(kt−δ) .

Note that the trial solution includes a possible phase shift δ.
See below for a caveat about this being linearly independent of yc .

• If f (t) = a sin kt or b cos kt, try the form yp = A cos kt + B sin kt.
Note that in general you have to allow for both sin and cos terms in yp ,
even if f (t) only contains one of them. This is equivalent to the phase
shift allowed for in the previous example.

• If f (t) = a0 + a1 t... + an tn try the form yp = A0 + A1 t... + An tn , where

the order n of the polynomials for yp and f (t) is the same.
Note that in general you should include all of the Ai terms even if some
of the ai coefficients are zero.

• If f (t) is the sum or product of any of the above, try the sum or product
of the corresponding functions.

The method is to substitute the trial function into the ODE and solve for
the coefficient (or coefficients) A and B. For the A sin +B cos example these
are found by matching separately the sin and cos terms in the ODE. For the
polynomial example the different powers of t are matched.
Note that unlike for the complementary function, all of the coefficients
will be uniquely determined from the ODE. Only the coefficients C1 and

70
C2 in the complementary function yc have to be determined from the initial
conditions. The physical reason for this will become evident later, when we
discuss transient solutions.
Example: Consider
d2 yp dyp
+ 2γ + ω02 yp = F0 cos Ωt
dt2 dt
for which the particular integral has the usual sine/cosine form solution:

yp = A cos Ωt + B sin Ωt

Differentiating and substituting into the ODE gives:

dyp d2 y p
= −AΩ sin Ωt + BΩ cos Ωt = −AΩ2 cos Ωt − BΩ2 sin Ωt
dt dt2
−Ω2 (A cos Ωt−B sin Ωt)+2γΩ(−A sin Ωt+B cos Ωt)+ω02 (A cos Ωt+AB sin Ωt) = F0 cos Ωt
Equating the cosine terms:

−Ω2 A + 2γΩB + ω02 A = F0

and the sine terms:

−Ω2 B − 2γΩA + ω02 B = 0
These are simultaneous equations that can be solved for A and B. From the
second equation we get the ratio of the coefficients:
A (ω 2 − Ω2 )
= 0
B 2γΩ
If γ = 0, the coefficient B = 0 and there is no phase shift between the driving
term and the particular integral.
Remember that one must also include the complementary function to get
the entire solution!

6.4 EXTRA: Dealing with repeated roots

If any combination of y1 , y2 or yp is not linearly independent they are said
to be repeated roots. As a specific example consider:
d2 y dy
2
− 2 + y = et
dt dt
71
Taking a trial function y = eλt for the complementary function leads to:

λ2 − 2λ + 1 = 0

which has only one root λ = 1, giving the same functional form for both y1
and y2 . An independent second solution can be obtained using D’Alembert’s
method, with y2 = uy1 . After taking the simplest possible form for u this
gives:

d2 u
=0 u=t y2 = tet yc (t) = (C1 + C2 t)et
dt2
By differentiating this and substituting into the ODE it can be shown that
both y = et and y = tet satisfy the ODE, and that these two parts of yc are
now linearly independent.

Moving on to the particular integral, it can be seen that it cannot be of

the form et or tet , because these are the two parts of the complementary
function. Instead we try:
yp (t) = At2 et
By differentiating twice and substituting into the ODE it can be shown that
the coefficient A = 1/2. The full solution of the given ODE is thus:

t2 t
!
y(t) = C1 + C2 t + e
2

72
7 Forced Oscillations
Forced oscillations are produced when an external force F (t) is applied to a
harmonic oscillator. Although F (t) can have any time-dependent form, by
far the most important case in practice is an oscillatory force, which can be
written in one of these forms:
F = F0 eiΩt F = F0 sin Ωt F = F0 cos Ωt
Note that the frequency of the applied force Ω is unrelated to the natural fre-
quency ω0 . The situation where Ω = ω0 is a special case known as resonance,
which is discussed in detail below.
The equation of motion including both damping and driving is just like
our previous damped SHO ODE except with a driving term on the RHS:
d2 x dx
2
+ 2γ + ω02 x = A0 eiΩt
dt dt
q
where we introduce the constants γ = b/2m, ω0 = k/m and A0 = F0 /m.
As we saw in the previous chapter, the general solution of the forced
oscillator will have a complementary function xc and a particular integral xp :
x = xc (t) + xp (t) = C1 eλ1 t + C2 eλ2 t + xp (t)

7.1 Complementary function in forced oscillations

The complementary function is the solution of the ODE ignoring the driving
term F (t). We have solved this damped harmonic oscillator problem already:
h i q
xc = e−γt C1 eωt + C2 e−ωt , ω≡ ω02 − γ 2
When γ > 0, the complementary solution is decaying. It can be overdamped,
underdamped, or critically damped, but in each case, after sufficient time,
the amplitude of xc will be negligibly small.
Intuitively, this must mean that the system behaviour is now set by the
driving force F (t), which is oscillatory and thus does not decay. Physically,
this means that after a sufficiently long time, the driving term overwhelms
any natural oscillations of the system, and thus determines the evolution. In
many situations, the behavior of interest is the stable, long-term solution, not
the initial transient state. Therefore in the case of forced oscillations with
non-zero damping, it is often possible to ignore the complementary part, and
just focus on the particular solution xp driven by F (t).

73
7.2 Particular integral in forced oscillations
The particular integral should have a similar form to the driving force. We
assume an amplitude/phase exponential form:

xp = Aei(Ωt−δ)

where δ allows for a possible phase shift between the driving force and the
response of the system. Note that the forced oscillations of the system are
at the frequency of the driving force Ω, not at the natural frequency of the
system ω0 . The driving force frequency can have any value Ω > 0.

Differentiating and substituting back into the 2nd order ODE:

A(−Ω2 + 2iγΩ + ω02 )ei(Ωt−δ) = A0 eiΩt

The eiΩt cancels, and we can collect real and imaginary terms into

(ω02 − Ω2 )A + 2iγΩA = A0 eiδ = A0 cos δ + iA0 sin δ

Equating the real and imaginary terms separately, we get two equations:

(ω02 − Ω2 )A = A0 cos δ 2γΩA = A0 sin δ

We can thus solve this for the amplitude A and phase δ. δ is given by
2γΩ
tan δ =
(ω02 − Ω2 )
To get the amplitude, we take the square of each equation and add them.
This is a handy trick when one has a sin and cos term, since sin2 + cos2 = 1.
This gives
(ω02 − Ω2 )2 A2 + (2γΩA)2 = A20
A0
=⇒ A(Ω) = q
(ω02 − Ω2 )2 + (2γΩ)2
Note that this amplitude is only a function of the driving frequency Ω, and
does not decay away as a function of time. This is because the driving force
supplies energy to the system that compensates for any loss of energy due
to damping. After a long time, the system will follow xp (t), and will have
forgotten about the decaying solutions xc (t).

74
EXTRA: This can also be done using the sine/cosine or amplitude/phase
form. For the sine/cosine form, we already did this as the example of a 2nd
order inhomogenous ODE in the previous section. To recap, the ODE is
d2 x dx
2
+ 2γ + ω02 x = A0 sin Ωt
dt dt
The assumed form for the particular integral is:
xp = A1 sin Ωt + A2 cos Ωt
Differentiating and substituting into the ODE gives:
dxp d2 xp
= A1 Ω cos Ωt − A2 Ω sin Ωt = −A1 Ω2 sin Ωt − A2 Ω2 cos Ωt
dt dt2
−Ω2 (A1 sin Ωt+A2 cos Ωt)+2γΩ(A1 cos Ωt−A2 sin Ωt)+ω02 (A1 sin Ωt+A2 cos Ωt) = A0 sin Ωt
Equating the sine terms:
−Ω2 A1 − 2γΩA2 + ω02 A1 = A0
and the cosine terms:
−Ω2 A2 + 2γΩA1 + ω02 A2 = 0
These are simultaneous equations that can be solved for A1 and A2 :
(ω02 − Ω2 )A0
A1 =
[(ω02 − Ω2 )2 + (2γΩ)2 ]
2γΩA0
A2 =
[(ω02 − Ω2 )2 + (2γΩ)2 ]

For the amplitude and phase form, the solution takes the form:
xp = A sin(Ωt + δ)
with
q A0 A2 2γΩ
A= A21 + A22 = q tan δ = = 2
(ω02 − Ω2 )2 + (2γΩ)2 A1 (Ω − ω02 )

These are the same as in the exponential version, since we assumed an

amplitude-phase form for the exponential xp to match the driving force.

75
7.3 Properties of forced oscillations
The damped, driven oscillator is a common occurence in many areas of
physics. Resonance, in particular, is an important situation. Here we ex-
amine some properties of forced oscillators, using the solution we derived in
the previous section.

7.3.1 Limits
At very low driving frequencies Ω ≪ ω0 , the natural frequency decays away
quickly, and the system then responds only to the slow driving fluctuations
with frequency Ω. This is known as the static limit. The amplitude and
phase are
A0 2γΩ
A(Ω ≪ ω0 ) ≈ 2 , δ ≈ tan−1 2 ≈ 0
ω0 ω0
The solution (ignoring the decaying complementary solution) is thus

A0 iΩt
xp (t) = e
ω02
This behaves like simple harmonic oscillator with a spring constant, except
with a frequency set by the driving frequency rather than the natural fre-
quency of the system.

At very high driving frequencies Ω, we have a situation where the natural

frequency is very slow, and thus it cannot respond to the rapidly fluctuating
driving force. To see this, we take the limit Ω ≫ ω0 :
A0 A2 2γ
A(Ω ≫ ω0 ) = tan δ = =−
Ω2 A1 Ω
In the limit Ω → ∞ both A and tan δ go to zero because the system cannot
follow the high frequency of the driving force, and we get

A(∞) = 0 δ(∞) = π

In this limit, xp = 0, and the solution is given only by the decaying com-
plementary solution. After the decay, the system is stationary, as the slow
natural frequency never allows a response to the extremely fast variations in
the driving force (which average out to zero on short timescales).

76
Figure 12: The amplitude of oscillations peaks at resonance where the driv-
ing frequency Ω is matched to the natural frequency ω0 of the system. As
the damping constant γ is increased, the resonant peak is lowered, and the
resonance width increases.

7.3.2 Resonance
The case when Ω = ω0 is known as resonance. Here, the system is being
driven exactly at the natural frequency. You probably have the intuition that
pushing a swing exactly at the top of its arc generates a greater amplitude;
this is an example of resonance. Hence resonance can drive large amplitudes,
limited only by the damping. This has important applications in atomic
physics and engineering, among other areas.
When Ω = ω0 , we have
A0
A(ω0 ) =
2γω0
Intuitively, the amplitude is inversely proportional to the damping constant
γ: as γ → 0, A → ∞.
In engineering, it is often interesting to consider the boost in amplitude at
resonance, as compared to the amplitude far from resonance. This is known
as the quality factor
A(ω0 ) ω0
Q≡ =
A(Ω ≪ ω0 ) 2γ

77
One can show that Q also represents the energy stored in the system divided
by the energy lost per radian. The quality is higher when the damping factor
γ is lower, and γ must be below 0.5ω0 to have any boost at all.

In atomic physics, a key concept related to resonance is the resonance width.

As the driving frequency Ω approaches the natural frequency ω0 , the system
increases its amplitude. A measure of this behavior is to ask: at which values
of Ω does the system achieve half of its maximum power?
We will show in the next section that power is proportional to A2 . So
we want the value of ∆Ω = Ω+ − Ω− where Ω± are the driving frequencies
on either side of ω0 where A2 = A2 (ω0 )/2, i.e. where the power is at half its
peak.
Looking at the denominator of the expression for A, half the peak power
is achieved when
ω02 − Ω2± = ±2γΩ±
(ω0 − Ω± )(ω0 + Ω± ) = ±2γΩ±
For small deviations from resonance, ω0 ≈ Ω± , so the sum term is (ω0 +
Ω± ) ≈ 2Ω± . The difference term is half the resonance width, since Ω+ and
Ω− are symmetric about ω0 . So dividing out 2Ω± from both sides leaves
1
∆Ω = ±γ
2
This is the resonance half-width. It is exactly the damping constant. Hence
the resonance width scales inversely with γ: more damping not only lowers
the peak amplitude, but smears out the resonance as well.
In terms of Q, we see that Q = ω0 /∆Ω. Hence the resonance width,
relative to the natural frequency ω0 , is just the inverse of Q.

Interestingly, the maximum value of the amplitude as a function of Ω is not

exactly at ω0 . Intuitively, you might have realized that the limits Ω ≪ ω0 and
Ω ≫ ω0 gave rise to different amplitudes, so the resonance curve cannot be
exactly symmetric. To determine the driving frequency where the amplitude
is maximized, we can differentiate the amplitude expression with respect to
Ω to find where it has its maximum value. This is equivalent to minimizing
the term inside the radical, i.e.
d h 2 i
(ω0 − Ω2 )2 + (2γΩ)2 = 0
dΩ
78
 2(ω02 − Ω2 )(−2Ω) + 8γ 2 Ω = 0 =⇒ −ω02 + Ω2 + 2γ 2 = 0
s
q 1
=⇒ Ωmax = ω02 − 2γ 2 = ω0 1 −
2Q2
So the driving frequency where the amplitude is maximized is slightly smaller
than the natural frequency, in the presence of damping. The amplitude itself
is
A0 Q
A(Ωmax ) = q = A(Ω = 0) q
2 1 − 1/4Q2
2γ ω0 − γ 2
So the boost in amplitude relative to when the driving frequency is zero is
somewhat larger than Q. Indeed, this is expected, since the boost at exactly
Ω = ω0 is just Q, while Ωmax maximises the boost.

Finally, one can compute the phase shift at Ω = ω0 :

π
tan δ = ∞ δ=
2
Hence the response of the system lags behind the driving force by π/2. In
fact, once can show that there is a continuous change in the phase between
0 (when Ω ≪ ω0 ) and π (when Ω ≫ ω0 ) as the frequency is scanned across
the resonance. The rate of this rotation is proportional to the width of the
resonance.

7.4 Power in Forced Oscillations

Without damping, there is no net energy required. However, the driving force
still provides power over the course of the oscillation. Recall that power is the
force F (t) times the velocity (ignoring the decaying complementary term):

dxp
P =F = mA0 eiΩt (AiΩ)ei(Ωt−δ)
dt
where we will take the real component of this to get the power. Without
damping (γ = 0), from §7.2 we have δ = 0 and A = A0 /(ω02 − Ω2 ):

mA20 Ω i2Ωt mA20 Ω

P = ie = [i cos 2Ωt − sin 2Ωt]
ω02 − Ω2 ω02 − Ω2

79
Taking the real component gives

mA20 Ω
P = sin 2Ωt
(Ω2 − ω02 )

This shows that the power depends on the amplitude of the driving force
squared. The power alternates between positive and negative over the course
of an oscillation, but the time-average of sin 2Ωt is zero, so there is no net
power transfer to the system. This is expected because there is no damping.

With damping, it is difficult to compute P (t) in general, but one can compute
A0
it at resonance where Ω = ω0 . Here, A = 2γΩ , and δ = π/2. Hence

dxp mA20 Ω iΩt i(Ωt−π/2)

P =F = e e
dt 2γΩ

One can show that eiΩt ei(Ωt−π/2) = cos2 Ωt, and setting Ω = ω0 , we get:

mA20
P (ω0 ) = cos2 ω0 t
2γ

The average of cos2 over an oscillation is 21 , so now there is a net power input
from the driving force:

F02 QF02
< P (ω0 ) >= =
4mγ 2mω0
This is the power required over each oscillation in order to overcome the
damping.

7.5 Summary: Forced, Damped Oscillations

The general forced (or driven), damped harmonic oscillator equation is:

d2 y dy
2
+ 2γ + ω02 y = F (t)
dt dt
q
where the natural frequency ω0 ≡ k/m and the damping frequency γ ≡
b/2m. The behavior of such a system depends on the relative values of ω0 ,
γ, and Ω.

80
• The solution consists of complementary and particular functions.

y(t) = yc (t) + yp (t)

• In the case of constant γ and ω0 , the complementary solution (with

RHS=0) corresponds to a damped oscillator:

yc (t) = A0 e−γt sin(ωt + ϕ)

where ω 2 = ω02 − γ 2 . A0 and ϕ are obtained from initial conditions.

• The case for a non-constant γ(t) and ω0 (t) does not have a general
solution, but if one solution y1 (t) is known, then y2 (t) can be obtained
using D’Alambert’s method.

• The complementary solution is transient, i.e. after some time it decays

away owing to the e−γt term. After some time, the system responds
only to the driving force F (t). If one is only interested in the long-term
behaviour, yc (t) can be neglected.

• The particular solution often takes a form similar to F (t). For instance,
if F (t) = F0 sin Ωt, then yp = A sin(Ωt + δ) is a good guess.

• F (t) fully specifies yp (t)! The initial conditions are not used for yp be-
cause they have already been used for the transient solution yc . Phys-
ically, only the transient term cares about the initial conditions; after
that is damped away, the system responds solely to the driving force.

• For the above F (t), one finds

F0 1 2γΩ
A= tan δ =
− ω02 )
q
m (ω02 − Ω2 )2 + (2γΩ)2 (Ω2

• In the static limit when the natural frequency is much smaller than the
driving frequency (ω0 ≪ Ω), the natural oscillations decay quickly, and
the system oscillates purely at the driving frequency with an amplitude
modulated by ω0 :
F0 iΩt
xp = e .
mω02

81
• Resonance occurs when ω0 = Ω, i.e. the system is driven at its natural
frequency. Here δ = π/2, so the system is 90◦ out of phase with the
driving. Without damping, this would drive oscillations to infinite
amplitude. With damping, the maximum amplitude is
F0
Amax =
2mγω0

• The resonant width ∆Ω spans the values of Ω that yield half the max
power: A2 = ±0.5A2max , i.e. ω02 − Ω2± = ±2γΩ± . For small damping,
this gives ∆Ω = Ω+ − Ω− = 2γ = ωQ0 , where Quality Q ≡ ω0 /2γ.

• The power required to drive the oscillator (after yc decays) when aver-
aged over its period 2π/Ω scales with the driving force squared:

dyp mA20
< P >= F v = F =
dt 4γ

82
8 Coupled Oscillations
Many physical systems can be described a collection of oscillators. Often, the
oscillators are not isolated, but are connected to each other in some way. Thus
the behaviour of coupled oscillators, and their response to perturbations, is
an important problem to solve.
Note: We will be using matrix algebra in this section, as this is the
most compact and intuitive way to represent coupled oscillations. It is worth
brushing up on basic matrix notation and manipulation; we will not cover
that here.

8.1 Coupled horizontal masses

Figure 13: Coupled oscillators, connected by springs to the walls and to each
other.

We consider a horizontal system of two masses m1 and m2 linked by three

springs between two fixed walls. We assume that the masses are resting on
a frictionless surface so we can ignore gravity, and only consider motion in
the horizonal direction of the springs.

The first mass m1 is connected to a wall by a spring constant k1 , and the sec-
ond to the other wall by a spring constant k3 . The two masses are connected
to each other by a spring constant k2 . The displacements of the masses from
their static equilibrium positions are x1 and x2 . The signs of the displace-
ments matter, and it is simplest to assume they are both positive to the right.
As in the earlier discussion of simple harmonic motion of a single mass, for
small oscillations we are interested in the changes in the forces due to the

83
displacements from equilibrium, i.e. we must solve for x1 (t) and x2 (t).

We can solve this using N2, by computing the forces on m1 and m2 from
their connecting springs. For m1 , there is a force to the left from k1 , and a
force to the right from k2 . For k1 , one can see that a positive displacement
in x1 will stretch spring 1, and thus result in a negative force:

F1L = −k1 x1

For spring 2, this depends on the stretching of the middle spring, which is
given by (x2 − x1 ). If this is positive (spring 2 is stretched), then this will
result in a positive force on m1 :

F1R = k2 (x2 − x1 )

Note that m1 doesn’t care about k3 ; it only responds to the forces directly
on its own centre of mass. The eventual solution for m1 will depend on k3
through the coupling, but it is important to recognize that, as a general rule,
the only forces a body cares about for N2 are the forces that are directly
acting on it.

Combining the two contributions gives the overall force in m1 :

d 2 x1
F1 = F1L + F1R = m1 = −(k1 + k2 )x1 + k2 x2
dt2
A similar exercise for m2 shows that the forces from springs 2 and 3 will be:

F2L = −k2 (x2 − x1 ) F2R = −k3 x2

Note that the F1R = −F2L must be true according to N3. So m2 feels a force

d 2 x2
F2 = F2L + F2R = m2 = −(k2 + k3 )x2 + k2 x1
dt2
The equations for F1 and F2 are coupled equations in x1 and x2 . It is con-
venient to represent these using using matrices, by writing them in the com-
bined form:
d2 x
!
x1
M 2 = −Kx x=
dt x2

84
where the 2 × 2 mass matrix M and the 2 × 2 force matrix K are given by:
! !
m1 0 k1 + k2 −k2
M= K=
0 m2 −k2 k2 + k3

This equation is like Hooke’s Law, except in matrix form. Note that M is
diagonal and both M and K are symmetric – you can convince yourself that
N3 implies that K will always be symmetric. We must solve this for the
vector x(t).

To solve this, we use our now-familiar technique: We assume there is a

harmonic form for the solution, with some amplitude A and frequency ω,
and then we obtain equations for these that will satisfy the original ODE.
We will assume no damping, so this is relatively straightforward. There will
be two modes, corresponding to two natural frequencies of the system; we’ll
call them A and B. But we will see that the coupling will result in these
modes being linear combinations of x1 (t) and x2 (t).

Let us assume the form x = Ceiωt for both modes, except with different
values of C and ω for each mode. Taking d2 x/dt2 results in a factor of −ω 2
as usual. Thus ω must satisfy

(−ω 2 M + K)x = 0

For this to be zero for any general x, the determinant of this matrix must
vanish:
|K − ω 2 M| = 0
It is conventional to define the eigenvalues of the system λ ≡ ω 2 . Writing
the determinant out explicitly, we have

k1 + k2 − λm1 −k2
=0
−k2 k2 + k3 − λm2

(k1 + k2 − λm1 )(k2 + k3 − λm2 ) − k22 = 0

m1 m2 λ2 − [(k1 + k2 )m2 + (k2 + k3 )m1 ]λ + (k1 + k2 )(k2 + k3 ) − k22 = 0
This is a quadratic equation on λ, which gives two eigenvalues λA and λB as
solutions, corresponding to the two natural frequencies of this system.

85
8.2 Equal masses and spring constants
We now consider some examples to develop intuition. We begin with the
simplest case, where we have all masses and springs the same: m1 = m2 = m
and k1 = k2 = k3 = k:
2k − λm −k
=0
−k 2k − λm

m2 λ2 − 4kmλ + 3k 2 = 0 =⇒ (mλ − 3k)(mλ − k) = 0

There are two solutions, with the eigenvalues:
k 3k
λA = ωA2 = λB = ωB2 =
m m
These represent two distinct normal modes of the coupled oscillations, with
two different frequencies ωA and ωB .

To solve for x(t), we now need the amplitudes. It turns out that there
will be only one independent equation to solve for two amplitudes, hence we
can only solve for the relative amplitudes of each mode. We call the result-
ing solution for x(t) as the eigenvectors of the system, with each eigenvector
associated with the particular eigenvalue. As we will see, the eigenvectors
tell us the way in which the two masses move relative to each other.

To determine the eigenvectors, let us go back to our original matrix equation

!
0
(−λM + K)x =
0
This matrix equation thus corresponds to two equations, one corresponding
to the top row and one for the bottom. The first equation comes from
multiplying the top row of the matrix −λM + K by x. This gives
x2 (−λm + 2k)
(−λm + 2k)x1 − kx2 = 0 =⇒ =
x1 k
In this case, you can show that the equation for the bottom row is identical,
so there is no separate information from this. We can obtain the specific
eigenvectors for each solution by inserting the eigenvalues λA and λB :
k 3k
λA = =⇒ x2 = x1 and λB = =⇒ x2 = −x1
m m
86
Hence the eigenvectors associated with these eigenvalues are:
! !
1 1 1 1
xA = √ xB = √
2 1 2 −1
√
where we have divided by 2 because it is conventional for the eigenvectors
to be normalised like unit vectors, so their amplitude is unity.

Thus we have the solution for the normal modes or eigenmodes of our sys-
tem. The eigenmodes correspond to linear combinations of x1 (t) and x2 (t)
that correspond to the natural behaviour of the coupled oscillator. In this
case, xA (t) = √12 (x1 + x2 ), and xB (t) = √12 (x1 − x2 ) Conversely, x1 (t) and
x2 (t) can be expressed as a linear combination of eigenmodes.
The two eigenmodes can be described physically as follows:
• Mode A has a the two displacements x1 and x2 being exactly the same:
x1 = x2 = AeiωA t
Thus the masses are moving unison, back and forth. The central spring
never stretches so it plays no role, while the outer springs are alternately
compressed and extended in antiphase with each other. The frequency
of this corresponds to the natural frequency of a single mass, ωA =
q
k/m, because both masses are moving effectively as a single mass.
• Mode B has the two displacements x1 and x2 being equal and opposite:
x1 = −x2 = BeiωB t
q
at a higher frequency ωB = 3k/m. Here, the central spring is alter-
nately compressed and extended, while the outer springs are alternately
extended and compressed in phase with each other, but in exact an-
tiphase with the central spring.
One requires initial conditions for the two masses to determine A and B.

8.3 Equal masses and different spring constants

Check for yourself that making the spring constants k1 = k3 but with a
different k2 leads to the solutions:
k1 (k1 + 2k2 )
λA = ωA2 = λB = ωB2 =
m m
87
and the eigenvectors are still:
! !
1 1 1 1
xA = √ xB = √
2 1 2 −1

For Mode A, where the two masses moved in unison, the middle spring played
no role, so changing k2 has no effect. The only change here is in the frequency
of Mode B.

It is instructive to check that the limits are behaving as expected physi-

cally. If k1 → 0, then the outer springs provide no force, so the system
reduces to two masses connected by a single spring. ωA → 0, so there is no
mode for motion in unison (since there are no walls to keep the masses in
place). The disappearance of Mode A can be thought of as the introduction
of a freedom of translation of the whole system in the x direction. Mode B
remains as the only
q vibrational mode with the two masses in antiphase, and
the frequency is 2k2 /m which is lower than the Mode B frquency previously.

In the limit k2 ≪ k1 , ωB ≈ ωA , the weak coupling of the central spring

can be ignored and each mass is separately connected to its own wall. Mode
B disappears, and only Mode A remains with each mass performing SHM
independently. The phase and amplitudes of the motions of the two masses
can be anything, since they are effectively not connected.

8.4 Coupled pendulums

We now add gravity to this problem, by coupling two pendulums. Here, we
assume that pendulum bobs with masses m1 = m2 = m are connected by
a spring with spring constant k. The equations of motion are the same as
we had for θ(t) for a single pendulum, except that we must now add a force
from the spring on each pendulum. If the displacements are very small, then
x = sin θ ≈ θ, so we can write the coupled ODEs for x1 and x2 as
d2 x 1 d 2 x2
m = −mω02 x1 + k(x2 − x1 ) m = −mω02 x2 − k(x2 − x1 )
dt2 dt2
q
where the natural frequency ω0 = g/L accounts for gravity. Note that for
pendulum 1, stretching the spring with a positive (x2 −x1 ) result in a positive
force, and by N3 the corresponding force on pendulum 2 must be negative.

88
Figure 14: Coupled pendulums, with a spring constant k connecting them.

As usual we take a trial solution of the form:

x1 = C1 eiωt x2 = C2 eiωt

From which we must determine C1 and C2 , and solve for possible values of
ω. Substituting this into the equations of motion gives:

−mω 2 C1 + mω02 C1 + k(C1 − C2 ) = 0

−mω 2 C2 + mω02 C2 − k(C1 − C2 ) = 0

These are two simultaneous equations in the amplitudes C1 and C2 , but also
with another two unknowns ω1 and ω2 which are the two possible values that
ω can take. Note that previously when we had no coupling spring, C1 and
C2 would have cancelled out – now, because of the coupling, the amplitudes
of these two oscillations are interdependent.

It is easiest to see the solutions by inspection. If we find two independent

solutions, we are done. The first solution is easily seen as:

ω1 = ω0 , C1 = C2

In this mode, the two pendulums move in phase with each other, so the
spring has constant length. Thus the spring constant k is irrelevant, and the

89
system oscillates in unison at the natural frequency.

The second solution can be obtained with a bit of algebra, but also can
be seen slightly less trivially by inspection:
s
2k
ω2 = ω02 + C1 = −C2
m
In this mode, the two pendulums move in antiphase with each other, with the
spring alternately compressed and expanded. As with the coupled springs,
the frequency of the antiphase mode is higher than the natural frequency.

EXTRA: We now consider unequal masses m1 ̸= m2 and equal spring con-

stants k1 = k2 = k3 = k. The matrices K and M have the symmetric
forms: ! !
m1 0 2k −k
M= K=
0 m2 −k 2k
and the determinant used for extracting the eigenvalues is:
!
2k − λm1 −k
|K − λM| =
−k 2k − λm2

(2k − λm1 )(2k − λm2 ) − k 2 = 0 3k 2 − 2k(m1 + m2 )λ + m1 m2 λ2 = 0

k
 q 
2
λ=ω = m1 + m2 ± m21 + m22 − m1 m2
m1 m2
We leave the eigenvectors relating x1 and x2 in the form:
m2 /m1
x1 = [−m1 λ + 2k] x2
k
m1 /m2
x2 = [−m2 λ + 2k] x1
k
Note that these ratios scale with the masses. This is significantly more com-
plicated than the case of equal masses, and becomes even more so with
unequal spring constants!

In the limit m1 ≪ m2 the motion of m2 can be neglected, and only the

higher mode survives from the +ve square root. A quick check shows that
the results in this section are still correct for m1 = m2 .

90
Figure 15: A simple model for a CO2 molecule. The springs both have spring
constant k.

8.5 Carbon Dioxide Molecule

A carbon dioxide molecule can be well modeled as a C atom linked symmet-
rically to two O atoms by spring constants k. In this case we have three
masses whose motion we must solve for, x1 , x2 , and x3 . But the basic pro-
cedure remains the same.
First, we must determine the forces. The force on the leftmost O atom
is k(x2 − x1 ), since an extension of the first spring pulls in the +x direction.
The force on the C atom x2 is −k(x2 − x1 ) + k(x3 − x2 ), and the force on the
right O atom is −k(x3 − x2 ). This will result in 3 × 3 matrices for M and
K. The M matrix is trivial; it’s just the three masses along the diagonal.
The K matrix can be constructed by looking at the force terms above.
The force for x1 has a −k dependence for x1 , and a +k for the depen-
dence on x2 . This means the upper left element of K will be +k (since
Md2 x/dt2 = −Kx), while the upper middle element will be −k. Similarly
the C atom’s force for x2 has −2k, so the diagonal term will be +2k, while
the off-diagonal terms for x1 and x3 are both +k, meaning the K matrix will
have −k for those terms. Similarly, the third O atom has +k for the diagonal
K term associated with x3 , and −k for its off-diagonal term for x2 .

Hence we have 3 × 3 matrices for the mass and spring constant:

   
1 0 0 1 −1 0
M = mO  0 µ 0  K = k  −1 2 −1 
   

0 0 1 0 −1 1

where µ = mC /mO . To obtain the eigenvalues, as usual we set the determi-

91
nant |K − λM| = 0:
1−λ −1 0
−1 2 − µλ −1 = 0
0 −1 1−λ
Note that we have divided out by k/mO for convenience, which means that
ω 2 = mkO λ in our case.

The determinant of a 3 × 3 matrix is in general not trivial, but in this case

it evaluates fairly simply:
0 = (1−λ)2 (2−µλ)−2(1−λ) = (1−λ)[(1−λ)(2−µλ)−2] = λ(1−λ)(−2−µ+µλ)
The three eigenvalues are thus λA = 0, λB = 1 and λC = 1 + 2/µ.

The first mode with λA = 0 is a trivial translation of the Centre of Mass

(CoM) with no oscillations and no stored energy. It is an important feature
of a free-standing system of coupled masses that the motion of the CoM re-
moves one of the degrees of freedom for the oscillations.

The frequencies of the two non-trivial modes B and C are:

s s
k k k
ωB = ωC = +2
mO mO mC
Getting the eigenvectors requires setting up three equations from the equa-
tion of motion [λM + K]⃗x = 0, and then solving these for the relations
between x1 , x2 , and x3 by plugging in the particular eigenvalues. Doing this,
we find that mode B (λB = 1) is represented by the eigenvector:
 
1
1 
xB = √  0 

2 −1

which has the two O atoms oscillating in antiphase with the central C atom
at rest.

Mode C is represented by the eigenvector:

 
1
1
xC = q  −2mO /mC 
 
2 2
2 + (4mO /mC ) 1

92
This has the two O atoms oscillating in unison and the central C atom
oscillating in antiphase and with a greater amplitude.
Note that neither Mode B or C displaces the CoM, because the transla-
tional eigenvalue is accounted for by Mode A. In Mode C, the motion of the
Carbon atom is just what is needed to counterbalance the motion of the two
Oxygen atoms.

93
9 Normal Modes
9.1 Superposition of normal modes: Beat frequencies
We’ve seen that coupled oscillators are described by normal modes that estab-
lish the characteristic frequencies of oscillation. Each mode is an independent
solution to the coupled equations of motion. The modes are linear combina-
tions of the motions of the individual masses. One can think of these normal
modes as setting the natural frequencies of the system. However, the masses
themselves do not (usually) oscillate at these frequencies!

It is conversely possible to describe the motion of the masses as a linear

combination, or superposition, of the normal modes. Let us return to the
case of two coupled masses. If xA and xB are the two eigenvectors of the
system, then the motion of the masses can be described by
x = AxA e−i(ωA t+δA ) + BxB e−i(ωB t+δB )
where A and B are the (unknown) amplitudes of the modes, and the δ’s are
the phases of the two modes at t = 0. This represents two separate equations,
written more explicitly as
x1 = Ae−i(ωA t+δA ) + Be−i(ωB t+δB ) and x2 = Ae−i(ωA t+δA ) − Be−i(ωB t+δB )
Both the amplitudes and the phases can be determined by the initial condi-
tions on x1 , x2 and their derivatives. With two masses, each described by a
2nd order ODE, we require a total of 4 boundary conditions.

So what does the motion of the masses look like? As a concrete example,
consider the case where mass m1 has initially been displaced to its maximum
amplitude A0 , while mass m2 is at its equilibrium position, and both masses
are initially at rest. The initial conditions at t = 0 are thus x1 = A0 , x2 = 0,
v1 = 0, v2 = 0. From these four constraints, we must determine the two
amplitudes and phases.

We can determine A and B by applying the initial conditions to the sum

x1 + x2 and difference x1 − x2 (this is a useful trick in many coupled equa-
tions). At t = 0, both quantities have a value of A0 . Thus
A0 iδA A0 iδB
A0 = 2Ae−iδA =⇒ A = e , A0 = 2Be−iδA =⇒ B = e
2 2
94
Figure 16: Combination of two modes with frequencies ω2 > ω1 . The total
amplitude is the sum of the two, which alternates from constructive interfer-
ence to destructive interference, forming a beat pattern. In the beat pattern,
the fast oscillations are at the mean frequency ω̄ = (ω1 + ω2 )/2. while the
enveloping beat frequency is ∆ω = ω2 − ω1 .

Plugging back in nicely cancels out the δ’s, leaving

A0 −iωA t A0 −iωA t
x1 = (e + e−iωB t ) x2 = (e − e−iωB t )
2 2
We can rewrite these using the mean frequency ω̄ = (ωA + ωB )/2 and the
beat frequency ∆ω = (ωA − ωB )/2, so ωA = ω̄ + ∆ω and ωA = ω̄ − ∆ω:

A0 iω̄t −i∆ωt A0 iω̄t −i∆ωt

x1 = e (e + e+i∆ωt ), x2 = e (e − e+i∆ωt )
2 2
Using e−ix = cos x + i sin x = cos x − ei3π/2 sin x, we can simplify this to

x1 = A0 eiω̄t cos ∆ωt, x2 = A0 eiω̄t+3π/2 sin ∆ωt

95
This shows that for each mode there is an oscillatory term at the average
frequency ω̄, which is modulated in amplitude by a lower-frequency envelope
given by the beat frequency ∆ω. The modes are π/2 out of phase.
Physically, the solution represents two characteristic frequencies, a high
frequency ω̄ and a low frequency ∆ω, and the motion is the superposition
of the two amplitudes at any given time. At one time, the two modes add
constructively, giving an overall maximum, while at a later time they add
destructively. The resulting pattern in shown in the Figure 16 above.
The low frequency ∆ω is known as the beat frequency. This beat pattern
is a generic outcome of superposing two normal modes. Systems with more
modes will have a more complicated pattern, but the general solution will
always be the sum of the amplitudes of the various modes at any given time.
In AM radio, the beat frequency is set to the station frequency, requiring
the radio to tune to the beat frequency in order to hear the station. Within
that, the signal is broadcast via the higher mean frequency. In contrast, in
FM radio, the beat frequency is applied to the frequency of the signal rather
than the amplitude.

9.2 Energy of Coupled Oscillations

Each mode conserves energy independently. To see this, consider the case
of two equal masses connected by unequal springs that we discussed in §8.3.
There we showed that the eigenvalues

k1 (k1 + 2k2 )
λA = ωA2 = λB = ωB2 =
m m
and the eigenvectors are just as in the equal k case:
! !
1 1 1 1
xA = √ xB = √
2 1 2 −1

We need to consider both the kinetic energy of the masses ( 21 mA2 ω 2 ) and
the potential energy stored in the springs ( 21 kx2 ), where A is the amplitude
of the displacement of each mass (and thus the maximum compression or
extension of each spring).
q
In Mode A with frequency ωA = k1 /m, the two masses move in unison.
Hence spring 2 does not extend, and all the ∆V is in the two outer springs

96
with constants k1 . This is balanced by a change in kinetic energy ∆T between
zero at xA = xB = A and a maximum at xA = xB = 0:
1 1
∆V = 2 × k1 A2 = k1 A2 ∆T = 2 × mA2 ωA2 = mA2 ωA2 = k1 A2
2 2
Hence energy is conserved between these two points.
For Mode B, the k1 springs are moving in opposition, while the middle k2
spring extends and contracts. At maximum displacement we have the max-
imum potential energy, corresponding to the k1 springs being (oppositely)
displaced by A, while the k2 spring is extended or contracted by 2A. Mean-
while the kinetic energy is difference between zero at xA = −xB = A and a
maximum at xA = xB = 0, is given by:
1 1 1
∆V = 2× k1 A2 + k2 (2A)2 = (k1 +2k2 )A2 ∆T = 2× mA2 ωB2 = mA2 ωB2
2 2 2
Since
(k1 + 2k2 )
ωB2 = =⇒ ∆T = (k1 + 2k2 )A2 = ∆V,
m
we see that energy is conserved in this mode as well.
This illustrates for this example that energy is conserved separately in
each of the normal modes. The higher frequency mode contains more energy.
This turns out to be a general result.

Conservation of energy can be proven more generally using matrix algebra.

For the normal coordinates xi for each mode, the kinetic energy is given by:
!T ! !
1 dxi dxi m1 0
T = M M=
2 dt dt 0 m2

where M is the symmetric time-independent mass matrix.

The potential energy is given by:
!
1 k1 + k2 −k2
V = xTi Kxi K=
2 −k2 k1 + k2

where K is the symmetric time-independent force matrix.

To prove conservation of energy we must differentiate the total energy E =

97
T + V with respect to time, and show that it vanishes. Because each mode
xi is independent, we can take the derivative of each and sum over i’s:
!T !T !T
d2 xi d2 xi
! ! !
dE 1 dxi 1 dxi 1 dxi 1 dxi
= M 2
+ M + xTi K + Kxi
dt 2 dt dt 2 dt2 dt 2 dt 2 dt
!T "  !T 
2 2
! # !
dE 1 dxi d xi 1 d xi dxi
= M + Kxi +  M+ xTi K
dt 2 dt dt2 2 dt2 dt
Note that the quantity in the first square bracket is just our N2 equation,
which vanishes! The second bracket is simply its transpose, which also van-
ishes: !T
d2 xi d2 xi
M 2 + Kxi = 0 M + xTi K = 0
dt dt2
This proves that energy is conserved separately for each normal mode i.
Note that it is not conserved for each mass independently, because the time
derivative assumed that there was no cross-terms between xi ’s.

9.3 Double Pendulum via Energy Conservation

Figure 17: Double pendulum. Angles are exagerrated for clarity.

As a final example, we consider a rather complicated looking problem

that turns out to have a relatively simple path to a solution using energy
conservation.

98
A double pendulum is formed from a simple pendulum of length L1 with
a bob of mass m1 , to which is attached a second simple pendulum of length
L2 with a bob of mass m2 . The first pendulum can oscillate about the
suspension point with angular coordinate θ1 , and the second pendulum can
oscillate about the mass m1 with angular coordinate θ2 . For simplicity we
just consider the case of equal lengths L1 = L2 = L here.

The complicated approach to a solution would be to use N2 for each case.

Polar coordinates are not applicable for m2 , so this would require a coupled
set of equations in x and y. However, we can use conservation of energy to
get the equations of motion and find the normal modes, without ever hav-
ing to explicitly construct the N2 force equations. The advantage of using
energies is that they are scalars, and hence additive without worrying about
directions.

Energy conservations requires that d(T + V )/dt = 0. The potential energy

V is the sum of terms from m1 and m2 :

V = m1 g(L − L cos θ1 ) + m2 g(L − L cos θ1 + L − L cos θ2 )

We work in the small angle approximation, as usual. Thus, cos θ ≈ 1 −

θ2 /2 =⇒ L − L cos θ ≈ Lθ2 /2. Thus this equation simplifies to
gL h i
V = (m1 + m2 )θ12 + m2 θ22
2
The kinetic energy T is a bit trickier. For the upper mass it is simple:
!2
1 dθ1
T1 = m1 L2
2 dt

For the lower mass, it is better to work in Cartesian coordinates x and y,

since the motion of m2 is no longer at fixed radius relative to the origin:
" ! !#
dx dθ1 dθ2
x = L sin θ1 + L sin θ2 =⇒ = L cos θ1 + cos θ2
dt dt dt
" ! !#
dy dθ1 dθ2
y = L cos θ1 + L cos θ2 =⇒ = −L sin θ1 + sin θ2
dt dt dt

99
Now  !2 !2 
1 dx dy
T2 = m2  + 
2 dt dt
 !2 !2 ! !
1 dθ1 dθ2 dθ1 dθ2 
T2 = m2 L2  + + 2(cos θ1 cos θ2 + sin θ1 sin θ2 )
2 dt dt dt dt
We can simplify this using the small angle approximation. In this case,
cos θ1 cos θ2 ≈ 1, and sin θ1 sin θ2 ≈ 0 to first order. Since the derivative
terms are already second order, the sin terms would be a 4th order correction,
which we can neglect. Hence

cos θ1 cos θ2 + sin θ1 sin θ2 ≈ 1 + O(θ2 ) ≈ 1

Thus the total kinetic energy is:

 !2 ! ! !2 
1 dθ1 dθ1 dθ2 dθ2
T = T1 + T2 = L2 (m1 + m2 ) + 2m2 + m2 
2 dt dt dt dt

We can write the kinetic and potential energies written in matrix form, using
the transpose (i.e. rows and columns interchanged):
!T ! !
1 dθi dθi m1 + m2 m2
T = L2 MT MT =
2 dt dt m2 m2
! !
1 m1 + m2 0 θ1
V = gLθiT MV θi MV = where θi =
2 0 m2 θ2
By energy conservation, dT /dt + dV /dt = 0. We take the derivatives of T
and V to get:
!T !T !T
L2 d2 θ i L2 d2 θi
! ! !
dE dθi dθi gL dθi gL dθi
= MT + MT + θiT MV + MV θi
dt 2 dt dt2 2 dt2 dt 2 dt 2 dt
!T "  !T 
L2 d2 θi L2 d2 θi
! # !
dE dθi g g dθi
= MT + MV θi +  MT + θiT MV 
dt 2 dt dt2 L 2 dt2 L dt
For this to always vanish, the terms inside the two brackets must be zero.
These give the same relation, because MT and MV are symmetric so the

100
terms are just transposes of each other (recall (AB)T = B T AT ). Taking the
first term, we obtain the desired equation of motion:

d2 θi g
MT 2
+ MV θi = 0
dt L
This corresponds to two coupled SHO-like equations, except with different
mass coefficients!

The solution works as before with the coupled springs, except that our deter-
minant will look a bit different. To compute a specific case, consider equal
masses m1 = m2 = m. As usual we take a trial solution of θ = Aeiωt , so the
second derivative pulls out a factor of −ωi2 = −λi . The eigenvalues λi thus
must satisfy (noting that all m’s cancel out):

g 2(g/L − λ) −λ
MV − λMT = =0
L −λ g/L − λ

g 2 g g2
 
2 − λ − λ2 = 0 =⇒ λ2 − 2 λ + 2 2 = 0
L L L
The two solutions of this for the normal modes are:
g √
λi = ωi2 = (2 ± 2)
L
The eigenvectors are obtained as usual by substituting these values back into
the original equation of motion. For the top row in the matrix, we have
g
 
2 − λ θ1 − λθ2 = 0
L
√
Plugging in the lower-frequency mode λA = Lg (2 − 2), we get
√ √
2(1 − (2 − 2))θ1 − (2 − 2)θ2 = 0
√ √ √ √
2(2 − 2)θ1 = (2 − 2)θ2 =⇒ θ2 = 2θ1
√
Similarly one can show that for λB , θ2 = − 2θ1 . Normalizing the eigenvec-
tors to unity, we have
! !
1 √1 1 1
√
θA = √ θB = √
3 2 3 − 2

101
Physically, this√shows that the lower pendulum has an amplitude of oscil-
lation that is 2 times the upper pendulum, for both modes. In the first
mode, they are in phase, while in the second, they are in antiphase. For large
enough amplitude, the system goes unstable and becomes chaotic.

EXTRA: For different masses m1 ̸= m2 the results get modified by a fac-

tor:
m2
µ2 =
m1 + m2
and the equation derived from the determinant reads:
2
g

−λ − λ2 µ2 = 0
L
The two solutions of this are:
g
λi = ωi2 =
L(1 ± µ)

The corresponding eigenvectors are:

s ! s !
1 µ 1 µ
θA = θB =
1 + µ2 1 1 + µ2 −1

One can see that this reduces to the equal-mass case when µ2 = 12 .

When m2 ≪ m1 , µ → 0 the two modes become degenerate and represent a

simple pendulum with mass m1 and length L.

When m2 ≫ m1 , µ → 1 and the higher mode ω2 → ∞. In this case the

lower mode represents a simple pendulum with mass m2 and length 2L.

END OF PART II

102
10 Central Forces
Many familiar forces in physics are central forces. A force can be represented
by a central field if it has the form:
⃗ r) = F (r)r̂
F(⃗
where r = |⃗r|, i.e. the force depends only on the distance from the origin, and
acts in the radial direction. It is spherically symmetric, and can be either
attractive or repulsive. The function F (r) can have any dependence on r, but
the two most familiar cases, gravity and electromagnetism, have F (r) ∝ 1/r2 .

A central force as defined above is necessarily a conservative force:

⃗ =0
∇×F ⃗ = −∇V
F
where the potential V (r) is well-defined at all points in space.

Previously we have considered central forces in the context of circular motion.

Now, we generalize this to non-circular motion, which allows us to describe
orbital dynamics. We begin by defining a number of useful quantities for this
situation. As with circular motion, it is most convenient to work in polar
coordinates.

10.1 Radial acceleration

We want to determine the equations of motion in the radial and tangential
directions. We did this previously for circular motion, assuming r and ω ≡
dθ/dt are constant. Now we repeat the exercise relaxing those assumptions.
In the plane of motion, recall that the velocity vector can be expressed
in terms of the cylindrical polar coordinates r and θ as:
d(rr̂) dr dθ
⃗v = = r̂ + r θ̂
dt dt dt
where we recall the derivatives of the unit polar vectors discussed earlier:
dr̂ dθ dθ̂ dθ
= θ̂ = − r̂.
dt dt dt dt
The second term of v represents a tangential velocity, which in the case of
no radial motion (dr/dt = 0) would be a constant velocity v = rω.

103
Differentiating the expression for the velocity:

d2 r dr dr̂ dr dθ d2 θ dθ dθ̂
⃗a = 2
r̂ + + θ̂ + r 2
θ̂ + r
dt dt dt dt dt dt dt dt
Collecting terms in the r and θ components of the acceleration:
 !2 
d2 r d2 θ
! !!
dθ dr dθ
⃗a =  2 − r  r̂ + r
2
+2 θ̂
dt dt dt dt dt

You can check that for constant radius d2 r/dt2 = dr/dt = 0, this reduces to
a centripetal acceleration and a tangential acceleration:
!2
dθ d2 θ
⃗a = −r r̂ + r θ̂
dt dt2

If the tangential acceleration is zero, then ω = dθ/dt is a constant and we

have circular motion. From now onwards we will assume that both ω and r
are variables which can be functions of both θ and t.

We now put the above equation for acceleration into N2, assuming a radial
force. The radial (r̂) acceleration equation is:

F (r) d2 r
ar = = 2 − rω 2
m dt
2
Note that there the radial acceleration is not just ddt2r , but there is an ad-
ditional term −rω 2 . This should be familiar: it is exactly the centripetal
acceleration discussed earlier. This comes in because we are working in a
non-inertial (rotating) frame.

10.2 Conservation of Angular Momentum

Meanwhile, the tangential (θ̂) acceleration is zero for a central force, so the
equation is
dω dr
aθ = 0 = r + 2ω
dt dt

104
The equation for aθ can be rewritten as the derivative of a product:
" #
1 2 dω dr 1d  2 
r + 2rω = r ω =0
r dt dt r dt

To be true for any r, the term inside the parantheses r2 ω must be conserved.
This is more familiar when multiplied by the mass m of the object:

mr2 ω = L = constant

You should recognise this as the angular momentum L of the mass m.

Not only can a central force not provide any θ̂ acceleration, it also can-
not provide an acceleration in the direction k̂ perpendicular to the initial
motion. This means that an object moving in a central force is confined to
move in the orbital plane of its initial motion.

Figure 18: Angular momentum of a planar orbit.

Angular momentum can be expressed as a vector quantity, as follows:

⃗ = m⃗r × (⃗ω × ⃗r) = m⃗r × ⃗v = ⃗r × ⃗p
L

The last definition is the most general, and can even be applied to non-
constant masses1 .
1
Some textbooks (including RDG) use r2 ω rather than mr2 ω for angular momentum.
In that case the vector product form is ⃗r × ⃗v. This is non-standard, however.

105
The orbital plane is defined by the directions of ⃗r and ⃗v. The angular mo-
mentum is then a vector pointing along the axis k̂ normal to the plane of the
radius and velocity. The angular velocity vector ω ⃗
⃗ = dθ/dt is also in that
direction.

We can now prove conservation of angular momentum more generally. In

⃗ are
the presence of a central force both the magnitude and direction of L
conserved. This can be shown as follows:
⃗
dL d⃗r d⃗p
= × ⃗p + ⃗r × ⃗ =0
= ⃗v × m⃗v + ⃗r × F
dt dt dt
⃗ × ⃗r = 0 since F
since ⃗v × ⃗v = 0, and F ⃗ = F (r)r̂ .

10.3 Moment of inertia and torque

In the linear case, momentum is velocity times mass. In the angular case
here, we can define an analogous quantity to mass known as the moment of
inertia I, that satisfies:
⃗ = I⃗ω
L
For the scalar case, I = mr2 . In general, the moment of inertia is a tensor
(matrix), which multiplies the vector ω ⃗ When we deal
⃗ to yield the vector L.
with solid bodies later, this will become important, but for a point mass the
scalar case is usually all we need.

If there is a non-central force in the problem, such as gravity from a sec-

ond object, then there can be a component Fθ . This introduces a change in
the angular momentum, which is thus no longer conserved. We define the
torque
⃗
dL d⃗r d⃗p
⃗τ ≡ = × ⃗p + ⃗r × ⃗
= ⃗r × F
dt dt dt
The torque is thus the moment of the force acting about the origin, such as
to change the angular momentum vector. For a constant moment of inertia:
⃗
dL d⃗ω
⃗τ = =I
dt dt
⃗ = m d⃗v for linear motion.
which is the angular equivalent to the equation F dt

106
10.4 Rotational Energy
We now consider the energy associated with angular motion. The kinetic
energy can be obtained from the velocity in cylindrical polar coordinates,
derived above. In vectorial form, T = 12 m⃗v · ⃗v, so
!2
1 1 dr 1
T = m⃗v · ⃗v = m + mr2 ω 2 = Tr + Tθ
2 2 dt 2

where we have split T into components associated with radial (Tr ) and an-
gular (Tθ ) motion. Tθ can be expressed in terms of the angular momentum:

1 1 1 L2
Tθ = mr2 ω 2 = Lω =
2 2 2 mr2
Because L is conserved, Tθ only has a radial dependence. Note that for a
central force, L is conserved, and Tθ has a 1/r2 dependence for any function
F (r) (not just for a 1/r2 force).

Because F⃗ and thus V (r) are radial, it is intuitive to think of the kinetic
energy Tr as being the energy associated with rolling down the (radial) po-
tential. Mathematically, we can thus introduce an effective potential:

1 L2
Veff (r) ≡ V (r) + Tθ = V (r) +
2 mr2
This effective potential depends only on radius. Thus Tr takes care of the
angular energy, while Veff (r) accounts for the radial energy. Conservation of
energy thus balances Veff (r) against Tr . This can be intuitively useful. Again,
this is valid for a central force of any form.

107
11 Orbital Motion
Let us choose a common form for the central force, an attractive 1/r2 force:
F (r) = −k/r2 , where k is a positive constant. In that case, one can imme-
diately see that the potential will be
Z ∞
k
V (r) = − F (r)dr = −
r r
A useful convention is to set the constant of integration via the boundary
condition V (∞) = 0, as we have done. Note that there is no potential dif-
ference in the θ or ϕ directions; the potential depends only on the radius.

The above form describes gravity for k = GM m between two masses M

and m, or the electrostatic force for k = Zze2 /(4πϵ0 ), where Ze and ze are
the two charges and ϵ0 = 8.854 × 10−12 (SI units) is the electric constant. For
now, we assume that the central object is fixed, and only the orbiting object
moves. This is a good approximation in many circumstances, for example
the Solar System. The goal is to determine a general solution for the motion
of the orbiting object.

11.1 Orbit Equation

Earlier, we derived equations for circular motion, and showed that the force
required to keep an object in a circular orbit was an attractive 1/r2 force.
Now we consider the general case of a non-circular orbit, and derive the key
equation that describes the family of orbits in a 1/r2 force field.

We begin with the general form for the acceleration in polar coordinates.
From N2, the force must provide both a centripetal acceleration and a change
of radial coordinate. We must return to the full equation for ⃗a in §10.1.1,
but we can no longer assume d2 r/dt2 = 0 as in circular motion. We equate
the radial component to ar = F/m to get
d2 r k
2
= − 2 + rω 2
dt mr
Putting this into the radial acceleration equation, and using ω = L/mr2 :
d2 r k L2
= − +
dt2 mr2 m2 r3
108
Note the trick we have used to remove the (unknown) ω from the equation,
in favor of the constant L and radius r; we will use this trick repeatedly.

We must now solve this differential equation for r(t). Since we have r in
the denominator, let us try changing variables to u ≡ 1/r:

1 dr 1 du
r= =− 2
u dt u dt
A handy trick is to change the independent variable to θ, which introduces
a factor of ω:
du dθ du du
= =ω
dt dt dθ dθ
which allows us to convert a time derivative into a θ derivative. Now we can
use our trick to get rid of ω, via ω = L/(mr2 ) = u2 L/m. This is particularly
convenient because the u’s cancel out, so we are left with an equation in
terms of constants L and m:
dr ω du L du
=− 2 =−
dt u dθ m dθ
To get the second time derivative, we again convert to a θ derivative and
replace ω by L:

d2 r L d2 u dθ Lω d2 u L2 u2 d2 u
= − = − = −
dt2 m dθ2 dt m dθ2 m2 dθ2
We can now substitute this back into our original ODE for r, and express it
in terms of u:
L2 u2 d2 u ku2 L2 u3
− 2 = − +
m dθ2 m m2
Multiplying through by −m2 /(L2 u2 ) leads to the simple form:

d2 u km
+ u =
dθ2 L2
This is known as the orbit equation. It is a second order ODE in the inverse
radius u = 1/r, as a function of the angle θ(t).

The orbit equation is just a forced (undamped) harmonic oscillator, with

109
the RHS corresponding to a constant driving term. To solve this, we try a
solution:
d2 u
u = A cos θ + B =⇒ = −A cos θ
dθ2
A more general solution would include a sin term, but since we are uncon-
cerned about the exact phase but rather just the periodic behavior, we follow
custom and only include a cos term. Plugging this back in to the ODE, we
see that
km
B= 2
L
Let us write the orbit equation back in terms of the radius r:
1 1 R
u= = (1 + e cos θ) =⇒ r =
r R 1 + e cos θ
where we have defined
1 L2 A
R≡ = e≡
B km B

Figure 19: The orbit equation gives conic sections for various choices of e.

This equation describes conic sections. The constant e is known as the ec-
centricity, which describes the variation in radius as a function of θ. The
constant R is a radius, corresponding to a circular orbit with e = 0, for
which L2 = kmR.

Sketch for yourself the behavior of r(θ), for various cases of e. One can
classify the type of orbits based on the eccentricity:

110
 • Circular: e = 0 when A = 0

• Elliptical: 0 < e < 1 when A < B

• Parabolic: e = 1 when A = B results in r → ∞ when θ = π

• Hyperbolic: e > 1 when A > B results in r → ∞ when θ < π

The closest distance of approach to the origin of the central force for any
orbit is at cos θ = +1, where r = R/(1 + e). This is known as the pericentre
Elliptical orbits have an apocentre at cos θ = −1 and r = R/(1−e). Parabolic
and hyperbolic orbits have no apocentre, as they are unbound, i.e. the motion
goes towards infinite radius. We discuss elliptical and hyperbolic orbits in
more detail in the next few sections.

11.2 Energy of Orbits

The total energy of an orbit in a 1/r2 field is given by E = Tr + V (r) + Tθ
(from §10.3), which in our case is
!2
1 dr k 1 L2
E= m − +
2 dt r 2 mr2
dr L du
Expressing this in terms of u = 1/r, and recalling that dt
= −m dθ
:
!2
L2 du L2 u2
E= − ku +
2m dθ 2m

The first term is the kinetic energy due to the angular velocity. which in-
creases the total energy above the minimum value for a circular orbit. The
last two terms represent the effective potential, written in terms of u:

L2 u2
Veff (u) = −ku +
2m
The effective potential (Figure 20) allows us to visualise how an object’s
energy changes with radius, enabling analogous intuition as for the normal
potential in the non-rotating case.

111
 Figure 20: Effective potential Veff (r).

• Going towards the origin (r → 0, or u → ∞), the angular momentum

term dominates and the effective potential Veff increases rapidly as r →
0 (like 1/r2 ). This forms a potential barrier at small radii, owing to
the centrifugal force. This should familiar since it is difficult to move
towards the centre of a spinning merry-go-roung. Only when L = 0 is
it possible for a mass m to reach the origin of the force.

• Far out (r → ∞, or u → 0) the potential term dominates, the angular

momentum term is small, and the effective potential gradually tends
to zero like −k/r at large radii.

• As a result, there is a effective potential minimum when dVeff /du = 0, or

L2 = km/u. Here, Veff = −ku/2. This is exactly the energy of a circular
orbit (§11.1). Hence a circular orbit is the lowest possible potential
energy orbit, and thus most stable. If an object is able to dissipate
energy without losing (much) angular momentum, it will settle towards
a circular orbit.

112
Related to the last point, one can (as we’ve done before) obtain the equation
of motion from requiring conservation of energy (dE/dt = 0):

L2 d2 u L2
" ! ! ! !#
dE dE du du du
=ω =ω 2 −k + 2u
dt dθ 2m dθ2 dθ dθ 2m dθ

We rearrange and set dE/dt = 0:

L2 ω du d2 u
" ! #
km
0= 2
− 2 +u
m dθ dθ L

The term inside the brackets must vanish, which straightforwardly results in
the orbit equation from before:

d2 u km
2
+u= 2
dθ L

113
12 Elliptical Orbits
Let us now consider the case of an elliptical orbit, where 0 < e < 1. This is
the typical general case for a bound orbit, and is applicable to a wide variety
of physical situations in an attractive potential. For this, we can write the
orbit in terms of R and e. But in many situations, it is useful to consider the
geometry of the ellipse, which (if you recall from Geometry class) is described
by the semi-major and semi-minor axes a and b. Therefore it is useful to be
able to relate R and e to a and b.

12.1 Geometry of the Orbit

Figure 21: Parameters of an ellipse. The ellipse has two foci, and the at-
tracting object sits at one of these (F ). The major and minor axes are the
longest and shortest distances across the ellipse, respectively.

An ellipse is defined geometrically as the locus of points where the sum of

the distances to the foci F and F ′ is equal to the major axis (so r + r′ = 2a in
Figure 21 above). Defining x and y with respect to the centre of the ellipse,
the equation for an ellipse is:
x2 y 2 b2
+ 2 =1 ϵ2 ≡ 1 −
a2 b a2
where 0 < ϵ < 1 is the ellipticity, and a and b are the semi-major and semi-
minor axes. Comparing to our polar coordinate system defined with the
attracting object at the origin, the object moving in an ellipse ranges from
pericentre at R(1 − e) to apocentre at R(1 + e), with respect to the focus.

114
This means the attracting object sits at one focus of the ellipse (F or F ′ in
Figure 21), not at the centre!

We first show that the ellipticity ϵ is identical to the eccentricity e from the
orbit equation. Consider the point y = b (i.e. along the minor axis). Call the
distance from the focus to this point rb . Since this distance is symmetric to
both F and F ′ , the definition of an ellipse means that rb +rb = 2a =⇒ rb = a.
We also know the distance from the centre to the focus is ea (see Fig-
ure 21). Together with the semi-minor axis, this makes a right triangle with
rb2 = b2 + e2 a2 . Therefore

b2 rb2 − e2 a2 rb2
ϵ2 ≡ 1 − = 1 − = 1 − + e2 = e2 .
a2 a2 a2
Hence ϵ = e, and from now on we drop ϵ and just use e for both.

Now we can relate R to a and b (and e). From geometry, the distance from
the focus to the pericentre rp = a−ae = a(1−e) (see Figure 21). Meanwhile,
from the orbit equation, the distance from the focus to pericentre occurs at
θ = 0, which is rp = R/(1 + e). We can equate these two:

R R
rp = a(1 − e) = =⇒ a =
1+e 1 − e2
Using e2 = 1 − b2 /a2 , we have b2 = a2 (1 − e2 ), which gives

R2 R2 R
b2 = 2 2
(1 − e 2
) = 2
=⇒ b = √
(1 − e ) 1−e 1 − e2
We can also express the angular momentum L in terms of geometrical pa-
2
rameters. One can see from the above that R = ba . Recall that R = L2 /km
from the orbit equation, so solving for L2 gives

kmb2
L2 = kmR =
a
Note that setting b = a = R (i.e. e = 0) recovers L = kmR for circular orbits.

Finally, let us determine the orbital energy in terms of geometrical param-

eters. Since we’ve proven that total energy is conserved for any orbit, we

115
can choose to compute E at any point in the orbit. It’s convenient to choose
the pericentre rp = R/(1 + e). (It’s equally straightforward to choose the
apocentre ra ; this is left as an exercise to the reader.)

The energy is E = T + V . At pericentre, for our inverse-square force,

V (rp ) = − rkp . The kinetic term is as usual 21 mvp2 , where vp is the veloc-
ity at pericentre. We can determine vp using the angular momentum. At
pericentre, ⃗r and ⃗v are perpendicular, so L = mvp rp . But we know that
L2 = kmR. Hence
kR kR
L2 = kmR = m2 rp2 vp2 =⇒ vp2 = =⇒ T (rp ) =
mrp2 2rp2

To get rid of R, we use R = rp (1 + e), so the pericentre kinetic energy is

krp (1 + e) k (1 + e)
T (rp ) = =
2rp2 2 rp

Hence !
k k (1 + e) k e−1
E=− + =
rp 2 rp 2 rp
Using rp = a(1 − e), this reduces to the very simple relation

k
E=−
2a
Hence the total energy of an elliptical orbit is negative (bound) and only
dependent on the semi-major axis! A circular orbit or an orbit with any
ellipticity with the same semi-major axis will have the same total energy.
You can convince yourself that using the apocentre E = − rka + 21 mva2 gives
the same answer.

12.2 Kepler’s Laws

The laws for planetary motion were deduced by Kepler in the early 1600s,
based on celestial observations by Tycho Brahe. They eventually led to
Newton’s formulation of the law of gravity later in the century. Here we
will show how Kepler’s three laws emerge from orbital motion under gravity.
Kepler’s Laws state:

116
 K1: The orbits of the planets in the solar system are ellipses
with the Sun located at one focus of the ellipse.

We’ve already seen that this is a general property of bounded orbits in an

attractive 1/r2 central force field.

K2: The radius vector ⃗r from the Sun to a given planet

sweeps out equal areas in equal times.

The area of a thin triangle along the orbit between θ and θ + dθ is given by:
1 dA 1 L
dA = r2 dθ =⇒ = r2 ω = = constant
2 dt 2 2m
Hence K2 is a direct consequence of the conservation of angular momentum.

K3: The square of the planet’s orbital period τ around the Sun
is proportional to the cube of the planet’s semi-major axis a.

The period of the orbit can be obtained by integrating around the orbit using
the angular momentum, using dA = 12 r2 dθ:

dθ m Z 2π 2 2m Z 2mA
L = mr2 =⇒ τ = r dθ = dA =
dt L 0 L L
where A is the area swept out over the orbit. The area of an ellipse is
A = πab, hence
2πmab
τ=
L
From the previous section the angular momentum of an elliptical orbit is
L2 = kmb2 /a, so

4π 2 m2 a2 b2 4π 2 ma3 4π 2 3
τ2 = = = a
kmb2 /a k GM

for the case of gravity where k = GM m. Hence τ 2 ∝ a3 . The constant

of proportionality depends only on the mass of the central object M ; it is
independent of the mass of the planet m or ellipticity e, which explains why
every planet orbiting the Sun obeys exactly the same relation. However, the
Moon orbiting the Earth would have a different constant of proportionality.

117
12.3 Hohmann Transfer Orbits
As an example of an orbital mechanics problem, let us consider orbit transfer.
We wish to send a probe to another planet, say Jupiter. Assuming we can
get the probe out of Earth’s gravitational field, how do we send it to Jupiter?

Figure 22: Hohmann transfer orbit from Earth to Mars. A boost is given at
perihelion to send the Mars probe into an elliptical orbit whose aphelion will
coincide with Mars’s arrival. Another boost is then needed to get the probe
into orbit with Mars, enabling landing.

It turns out the most energy-efficient way to do this is to put the probe onto
an elliptical orbit, such that the closest point to the Sun is at the Earth, and
the farthest point is at Jupiter’s orbit. This is known as a Hohmann (1925)
transfer orbit. Assuming the probe is already outside Earth’s gravity, we
must therefore do two velocity boosts:
1. Boost the rocket tangentially from Earth’s orbit into the required el-
liptical orbit;
2. When it reaches Jupiter, boost the spacecraft’s velocity tangential to
Jupiter’s orbit, such that it has Jupiter’s velocity (otherwise Jupiter
will smash into it at high speed!).

118
What velocity boosts do we need to provide at each time, and how long will
the journey take2 ?

We first need to determine the required elliptical orbit for the probe. We
want perihelion at Earth’s orbit, and aphelion at Jupiter’s. We assume
for simplicity that the planetary orbits are circles at radii RE = 1 AU
and RJ = 5.2 AU (1 AU≈ 1.5 × 1011 m is an “astronomical unit”). Since
L2 = kmR,
q and L = mvR for a circular orbit at radius R, for k = GM m we
get v = GM/R. Hence the velocities of Earth and Jupiter are given by
s s
GMSun GMSun
vE = = 29.8 km/s; vJ = = 13.1 km/s
RE RJ

using G = 6.67 × 10−11 m3 kg−1 s−2 and MSun = 2 × 1030 kg.

Our probe’s elliptical orbit must be such that perihelion is at Earth’s radius,
and aphelion is at Jupiter’s radius. Recall that for an ellipse the pericentre
distance is a(1 − e), while the apocentre distance is a(1 + e). Therefore our
probe’s orbit must satisfy

RE = a(1 − e) = 1 AU RJ = a(1 + e) = 5.2 AU

We can solve these two equations for a and e to get the geometric parameters
of the probe’s orbit:
a = 3.1 AU e = 0.677

To determine the boost we need, we must find the velocity that the probe
should be at in order to be on the required orbit. Let us first determine
this velocity v1 at perihelion (RE ). Recall that v(RE )2 = GM (1 + e)/RE .
2
Meanwhile
√ for the Earth, e = 0, so vE = GM/RE . Hence we can write
v1 = vE 1 + e = 38.6 km/s. Therefore, the boost in velocity needed to get
from vE up to v1 is vboost,1 = v1 − vE = 8.8 km/s. This is the first boost we
need.

We do the analogous computation at RJ , now at aphelion√ra = RJ . The ve-

locity that the probe has at aphelion is, similarly, v2 = vJ 1 − e = 7.4 km/s
2
The launch must be done at a very special time, when the Earth and Jupiter are in
exactly the right position so that launching the probe from Earth will exactly meet Jupiter
when it gets to the probe’s aphelion.

119
(not the minus sign compared to perihelion). Hence the second boost needs
to be vboost,2 = vJ − v2 = 5.7 km/s. These are the two velocity boosts we
need from our rocket, the first at Earth’s orbit and the second at Jupiter’s, to
put our probe into the same orbit as Jupiter and be able to land it smoothly.

How long will it take to make the journey? Since the orbit goes from peri-
helion to aphelion, it is half the orbital period of the elliptical orbit. Using
a = 3.1 AU and Kepler’s 3rd Law,

4π 2 a3 τ
τ2 = =⇒ t = = 2.7 years
GMSun 2
This Hohmann transfer orbit can be shown to be the least energetic way to
transfer orbits (see RDG P.187 for proof). NASA uses this type of orbit to
send probes to other planets, after an initial boost to escape Earth’s gravity.
In detail, several TCMs (trajectory correction maneuvers) are used along
with way, since planetary orbits are not perfect circles.

12.4 Two-Body Orbits

So far we have considered orbiting bodies whose mass is negligible. We now
extend this to consider the mutual orbit of two bodies, where the masses
m1 ̸= m2 ̸= 0. The only force in the system is the mutual attraction of a
central force (e.g. gravity).

The Centre-of-Mass (CoM) of a two-body system is:

⃗ = 1
R (m1⃗r1 + m2⃗r2 )
(m1 + m2 )

In the absence of external forces the CoM is at rest, and we set this to be
⃗ = 0). The separation of the two bodies is:
the origin (R

⃗r = ⃗r1 − ⃗r2

By setting ⃗r1 = ⃗r + ⃗r2 one can straightforwardly show

m2 m1
⃗r1 = ⃗r ⃗r2 = − ⃗r
m1 + m2 m1 + m2

120
There is a central force F (r) acting between the two bodies, and from N3 we
know that
⃗ 1 = F (r)r̂
F ⃗ 2 = −F (r)r̂
F
The motion of each body is thus specified by:
d2⃗r1 d2⃗r2
m1 = F (r)r̂ m2 = −F (r)r̂
dt2 dt2
We can rewrite these both in terms of ⃗r, as follows:
d2⃗r1 d2⃗r2 m1 m2 d2⃗r m1 m2 d2⃗r m1 m2 d2⃗r
2F (r)r̂ = m1 − m 2 = + = 2
dt2 dt2 m1 + m2 dt2 m1 + m2 dt2 m1 + m2 dt2
Hence the relative motion of the two bodies is:
d2⃗r 1 m1 m2
2
= F (r)r̂ µ≡
dt µ m1 + m2
µ is known as the reduced mass of the two-body system, called so because
this mass is always less than the more massive of the two bodies.

The above equation shows that the two-body system behaves exactly like
the one-body system we have already studied extensively, except with the
(constant) reduced mass µ replacing the mass m, and both objects orbiting
a common centre of mass. For equal masses, µ = m/2, while for m1 ≪ m2 ,
µ ≈ m1 and the motion becomes a one-body problem with m2 at rest.

In terms of energetics, we can define the potential

2
⃗ d ⃗r
F(r) = µ 2 = −∇V (r)
dt
while the kinetic energy of the two-body system can be written in terms of
the reduced mass and the change in the separation vector:
 !2 !2  !2
1 d⃗r1 d⃗r2 1 d⃗r
T = m1 + m2 = µ
2 dt dt 2 dt

There can be both radial and tangential contributions to the kinetic energy,
corresponding to changes in the magnitude and direction of ⃗r. Again, these
are exactly as for the single mass case, with m → µ.

121
12.5 Extra: L and E of Two-body Orbits
We can show more explicitly that two-body orbits behave like one-body or-
bits with a reduced mass, and the separation vector replacing the radius, by
computing the angular momentum and energy of the system.

For a circular orbit of the two bodies about their CoM the centripetal forces
are:
F1 = m1 r1 ω12 F2 = m2 r2 ω22
These must be equal and opposite forces. Substituting for r1 and r2 :
F1 = −F2 = µrω 2 ω = ω1 = ω2
The angular velocity ω has to be the same for both masses. Putting in a
gravitational force between the two masses:
G(m1 + m2 ) r3
ω2 = τ 2 = 4π 2
r3 G(m1 + m2 )
The total angular momentum of the circular orbit is:
L = (m1 r12 + m2 r22 )ω = µr2 ω
The total energy of a circular orbit is given by the effective potential:
k L2 k
U =− + 2
=−
r 2µr 2r
The CoM is the centre of the orbit, the reduced mass replaces the single
masses, and the separation of the two masses acts as the radius of the orbit.

This treatment can be extended to the discussion of elliptical two-body orbits

where there is an additional radial kinetic energy term:
!2
1 dr
Tr = µ
2 dt
An application of these results is the search for exoplanets. As they orbit
around a distant star, they create a small “wobble” of the star, with an
amplitude proportional to the mass of the exoplanet and the distance between
the star and the exoplanet. The wobble can be measured as a Doppler shift
using extremely high-resolution spectroscopy. This is how exoplanets were
first discovered.

122
12.6 Precession
In the Solar System, the Sun is the dominant mass, but the other planets
(particularly Jupiter) are not negligible. The impact of the planets on each
others’ orbits causes a gradual precession of the perihelion, i.e. the orbits do
not exactly close on themselves after one revolution. To treat this properly,
one would have to solve the system including a non-central force from other
planets, which is quite challenging, and usually done numerically.

A common trick in physics when considering small perturbations is to as-

sume a more mathematically tractable form that slightly departs from the
original form, and hope that for small perturbations this provides a sufficient
effective theory. This is how we treat precession3 .

In our model, the effect of the other planets is treated as a screening correc-
tion which modifies the gravitational force due to the Sun, analogous to the
one used in the Bohr model to account for a reduction of the Coulomb force
due to the inner electron shells. The modification of the 1/r2 force law can
be written as a small change α in the power of r:
k
F (r) = −
r(2+α)
A common alternative choice is to write this as an exponential:
ke−αr
F (r) = −
r2
We leave it to the reader to show that the Taylor expansions of these two
forms are equivalent to first order in α. We will use the former version here.

We must now obtain the orbit equation for our new force, F (u) = −ku2+α ,
where u ≡ 1/r as before. This can be derived in a manner nearly identically
to before. However, in the last step before, we cancelled the u2 terms, which
happened because we had assumed F (u) = −ku2 . Now, we end up with an
extra term of uα , as in

d2 u mF (u) d2 u kmuα uα
2
+u=− 2 2 =⇒ + u = =
dθ Lu dθ2 L2 R
3
For a more general treatment of this problem see RDG P167-169.

123
with R = L2 /km as before. For α = 0, we know the solution is u = (1 +
e cos θ)/R. But now, let us assume a more general solution:
1 d2 u 1 d2 ϵ
u= [1 + ϵ(θ)] =⇒ =
R dθ2 R dθ2
For α = 0, ϵ(θ) = e cos θ. But now we must solve for ϵ(θ). Substituting the
2
u and ddθu2 terms, we have

1 d2 ϵ 1 R−α (1 + ϵ)α d2 ϵ 1
 
2
+ (1+ϵ) = =⇒ 2
+1+ϵ ≈ R−α 1 + αϵ + α2 ϵ2 + ...
R dθ R R dθ 2
We keep only first-order terms in α, and we assume that R−α ≈ 1. This gives
d2 ϵ
= −(1 − α)ϵ
dθ2
This should be familiar as a SHO for ϵ(θ). As usual, the solution is
√
ϵ(θ) = A cos(ωθ + δ) ω = 1−α

We need boundary conditions to determine the amplitude A and phase δ.

We can use the condition that if we set α = 0, we should recover our usual
1/r2 force solution: ϵ = e cos θ. For this to be true, we must have A = e and
δ = 0. Hence the solution for our perturbed orbit is:
1
u= (1 + e cos ωθ)
R
with (to first order in α):
α
ω ≈1−
2
The frequency of the perturbed orbit is not quite unity, but is slightly lower
(if α > 0). This means that when θ = 2π, the planet will not have come back
to the same location when θ = 0, but rather will need to go a bit farther in
θ to get back to perihelion where cos(1 − α2 )θ = 1.

Thus the perihelion at θ ≈ 2π will precess by an angle

∆θ = απ radians/orbit

The precise value of α must be determined based on the parameters of

the perturbing planet(s). Note that this is different from precession of the

124
Figure 23: Precession of a planetary orbit. The orbital phase lags a bit each
time around, requiring the planet to go further before reaching perihelion.

equinoxes, which owes to a precession of the Earth’s spin.

EXTRA: An interesting application of this is the precession of the peri-

helion of Mercury. This is one of the key early tests of General Relativity.
For Mercury it is observed that:

∆θM = 574”/century

In the 19th century it was found that, assuming Newtonian gravity, the effects
from the other planets accounted for 531”/century, leaving a mystery as to
what accounted for the other 43”/century. General Relativity naturally does
this, by introducing a relativistic term in the orbit equation that is inversely
proportional to c2 , where c is the speed of light:

d2 u GM m2 3GM 2
+ u = + 2 u
dθ2 L2 c
The relativistic correction is ∝ 1/r2 , so this can be interpreted as a perturbing
1/r4 force. Since this drops rapidly with radius, the precessions for more
distant planets are much harder to detect; hence GR also explained why
the other planets do not show significant deviations from Newtonian in their
precession.

125
13 Hyperbolic orbits
13.1 Orbital geometry

Figure 24: Diagram of hyperbola. An attractive force at focus F2 yields the

left orbit; for now we only consider this part of the diagram.

Hyperbolic orbits describe an unbound orbit in a central force, with e > 1.

The equation of a hyperbola is given in Cartesian coordinates by:
x2 y 2 b2
− 2 =1 e2 = 1 +
a2 b a2
where x and y are taken relative to the symmetry point C of the hyperbola.
Note the changes of signs compared to the ellipse!

As shown in Figure 24, a hyperbola has two branches around two foci F1
and F2 . By convention we place the attracting mass at F2 . Similar to an
ellipse, the distance from the centre to F2 (or F1 ) is ae. An orbiting object
only uses one branch depending on the sign of the force: An attracting force
results in the left hyperbola, while a repulsive force yields the right hyper-
bola. For now let us consider an attractive force (left branch).

The closest approach is when y = 0, and thus x = ±a. From this, the

126
distance of closest approach between the orbiting and attracting mass is
ae − a = a(e − 1). This is like the analogous equation for an ellipse, only
with e − 1 instead of 1 − e.

The orbit equation is independent of the nature of the orbit. Hence, just
as in the elliptical case, rmin = R/(1 + e). This allows us to relate the
semi-major and semi-minor axes to R and e via
R R b2
a= 2
b= √ 2 =⇒ R =
e −1 e −1 a
which is exactly as in the ellipse except with (1 − e2 ) → (e2 − 1).

What is the geometric intepretation of b? To see this, consider the asymptote

(blue dashed line in above diagram) to F2 . Define the angle of the asymptote
to the major axis as α. Then, the asymptote is defined by y = −x tan α.
Now, as x → ∞, the asymptote joins the hyperbola. So if we rewrite the
hyperbola equation as
1 tan2 α 1
2
− 2
= 2
a b x
and take x → ∞ (so the RHS vanishes), we get
b b b b
tan α = =⇒ sin α = √ 2 = √ =
a a + b2 a2 e 2 ae
Now, the sin represents the length of the opposite side divided by the hy-
potenuse. Going back to Figure 24, we see that b then represents the vertical
distance between the asymptote and F2 . In other words, b is the closest dis-
tance the object would have had if it passed by F2 without any deflection. b
is known as the impact parameter.

Another important quantity is the deflection angle θ. The object initially

approaches along an asymptote, but leaves along a different asymptote. The
difference in angle between these asymptotes is the deflection angle. Simple
geometry shows that the deflection angle is θ = π − 2α.

13.2 Angular momentum and energy

We now compute the angular momentum and energy of a hyperbolic orbit.
L2 of a hyperbolic orbit does not change from the elliptical case, since it

127
comes directly from the orbit equation, independent of e:
kmb2
L2 = kmR =
a
Since L is conserved, we can evaluate it at any point on the orbit. Let us
consider r = ∞, where the perpendicular distance from the asymptote is the
impact parameter b. Here we have

⃗ = m⃗r × ⃗v =⇒ L = mbv∞ =⇒ v 2 = k
L ∞
ma
We can determine the velocity at the minimum distance rmin = a(e − 1) via
conservation of angular momentum:
bv∞
L = mbv∞ = ma(e − 1)v(rmin ) =⇒ v(rmin ) =
a(e − 1)
Now we can use b2 /a2 = e2 − 1 = (e + 1)(e − 1) and v∞
2
= k/ma , to get

2 b2 v∞
2
k (e + 1)
v (rmin ) = 2 2
=
a (e − 1) ma (e − 1)
Note again that these differ from the ellipse by the change from 1 ± e to e ± 1.

Finally, we compute the (conserved) energy of the orbit. We can choose

to evaluate E at either rmin or ∞. At rmin , the velocity is maximized, so
dr/dt = 0, which means Tr = 0. Hence E = V (r) + Tθ :
k 1 k k e+1 k k
E=− + mv 2 (rmin ) = − + = (−2 + e + 1) = +
rmin 2 a(e − 1) 2a e − 1 2a(e − 1) 2a
The total energy of a hyperbolic orbit is positive and inversely proportional
to the semi-major axis, as in the elliptical case.

We can equivalently derive this by noting that as r → ∞, the energy be-

comes purely kinetic since the potential goes to zero. Hence:
1 2 m k k
E = mv∞ = =
2 2 ma 2a
This shows that the E and L of an orbit are directly relatable to a and b.
Indeed, it is customary to specify orbits in central forces by (E, L).

128
Extra: The repulsive force branch of a hyperbolic orbit has a point of closest
approach to F2 at:
k (e − 1)
rmin = a(e + 1) v 2 (rmin ) =
ma (e + 1)
Remembering that the potential is now repulsive and hence positive, the total
energy is given by:
k L2 k
E=+ + 2
=+
rmin 2mrmin 2a
The total energy is the same for both branches of the hyperbola!

The asymptotes are still at the half angle α, and b is again the impact pa-
rameter (the perpendicular distance between the asymptotes and the source).
Hence these relationships are also the same for both branches:

2 k
v∞ = L = mbv∞
ma

13.3 Example: Unbounded Comet

A comet of mass m is approaching the Sun from a great distance. At this
time it has a constant speed V and is moving in a straight line whose per-
pendicular distance from the Sun is b. It completes one unbounded orbit
around the Sun and returns to a great distance. For the special case when
V 2 = 4GM/3b, show that the distance of closest approach to the Sun is
rmin = b/2, find the velocity at this point, and determine the angle through
which the comet is deflected by the Sun.

Since the comet has a positive total (kinetic) energy as r → ∞ the path
must be a hyperbola with e > 1. It is described by:
1 1 b2
= (1 + e cos θ) R=
r R a
The total energy is (specifying k = GM m for gravity):
k GM m 1 GM
E=+ = = mV 2 =⇒ V 2 =
2a 2a 2 a
129
Comparing with the specified initial condition V 2 = 4GM/3b, as see that

3b b2 25 5
a= =⇒ e2 = 1 + 2 = =⇒ e =
4 a 9 3
Thus the distance of closest approach is:
b
rmin = a(e − 1) =
2
The velocity v of the comet at this point can be found using angular momen-
tum conservation:

L = mvrmin = mV b =⇒ v = 2V

Finally the angle through which the comet is deflected by the Sun can be
determined by looking at the half angle between the asymptotes:
b 4
tan α = = =⇒ α = 53◦
a 3
The angle of deflection is:

θ = π − 2α = 74◦

130
14 Elastic Scattering
Particle collisions, and the resultant scattering, are a fundamental aspect for
many problems in physics, particularly in nuclear and particle physics. In this
process, an object is fired at an another object, the two interact or collide,
and one or both objects are scattered. The interaction can be modeled as a
repulsive central force, allowing us to use the machinery we have developed
for hyperbolic orbits in the previous chapter, except now we will use the
repulsive branch of the hyperbola. Scattering processes with no energy loss
are known as elastic collisions; we will only consider elastic collisions here.

14.1 Impact Parameter and Scattering Angle

In particle physics, the incoming particles are part of a beam. A beam parti-
cle starts at r = ∞ with a velocity v∞ . It travels towards the target parallel
to a beam axis that passes through the centre of the target, at a transverse
distance given by the impact parameter, b.

The beam particle is scattered by the target through a scattering angle θ.

Assuming a central force, the scattering is symmetric about the beam axis,
so the polar angle ϕ can always be integrated over 2π. If θ = 0, there is no
scattering, while θ = π occurs from a head-on collision. If the force depends
inversely on distance, then the smaller the impact parameter, the more scat-
tering there will be. Hence b and θ are anti-correlated.

The initial conditions are usually specified as the total energy E and an-
gular momentum L of the beam particle. These are related to b and the
initial velocity v∞ via
1 2
E = mv∞ L = mbv∞
2
These two quantities are conserved during motion in a central force. For a
hyperbolic orbit, E = k/2a, so these two quantities can be used to solve for
the geometric parameters of the orbit a and b.

14.2 The Differential Cross-Section

The total cross-section of a target is the incident area over which a beam par-
ticle will scatter (i.e. where θ > 0). For instance, for a hard target of radius

131
Figure 25: Scattering experiment setup. The incident beam at an impact
parameter b is deflected by an angle θ. A change of db encompassing an
incident area dσ results in a change of −dθ in the scattering angle, which
results in an outgoing differential area dΩ indicated by the shaded region.

R (and an infinitesimally small beam particle), this is simply the face area
of the target, e.g. πR2 . However, for a subatomic particle, the cross-section
can have a more complicated form. A key goal of atomic and nuclear physics
is to determine cross-sections, which can help constrain internal structure.
This can be done by sending beams of particles towards the target at a range
of impact parameters, and measuring the resulting deflection angles.

Consider a small range of impact parameters from b → b+db, as in Figure 25.

The differential area of the beam is dσ = 2πbdb; this would be the scattering
cross-section if the entire ring was scattered. The particle then scatters into
an angle θ → θ + dθ. Far from the target, simple geometry as shown in
Figure 25 shows that the corresponding area is dA = 2π(r sin θ)(rdθ). The
solid angle covered by this area is dΩ = dA/r2 = 2π sin θdθ. Hence beams
sent into the target within the area dσ end up coming out spread out over a
solid angle dΩ.

We define a differential cross-section σd (θ), which as shown in Figure 26 is a

132
Figure 26: The differential cross section σd (θ) quantifies how the area dσ in
the incident plane maps onto an outgoing solid angle dΩ.

function that relates dσ to dΩ:

dσ 2πbdb b db
σd (θ) ≡ − =− =−
dΩ 2π sin θdθ sin θ dθ
We include a negative sign because b and θ are anti-correlated, so db/dθ < 0;
the negative sign thus makes σd (θ) positive. Note that since b (and db) have
units of distance, σd (θ) has the units of area.

The total (or integrated) cross-section of the target σT is the integral of

dσ over the full solid angle dΩ = dϕdθ, where the scattering angle goes from
π → 0 since it runs opposite to the polar angle:
Z 0 Z 2π Z π
Z
dσ
σT = dσ = dϕdθ = 2π σd (θ) sin θdθ
π 0 dΩ 0

where the integral over dϕ gives 2π since there is no ϕ dependence (i.e. the
scattering is symmetric around the beam axis), and the negative sign goes
away by flipping the limits from π → 0 to 0 → π. σT , like σd (θ), has dimen-
sions of area.

The total cross section is often infinity for infinite-range forces like grav-
ity or electromagnetism, since any b value will give at least some deflection.
So it is not so useful. But σd (θ) can still be instructive because the beam
is not infinite, but typically only extends up to some bmax . In this case, one
can use σd (θ) to determine what the expected cross-section is as a function

133
of the (measurable) scattering angle, by integrating from b = 0 → bmax , or
alternatively from θ = π → θmin . This tells you the relative area of the
target that will scatter into a given angle θ. Hence one can calculate σd (θ)
based on some assumed physical scenario and a given beam size, and test
that prediction by measuring how often a particular value of θ is seen in the
detector. This is how σd (θ) is used to test a particular physical scenario via
a scattering experiment.

14.3 Hard-body scattering with a heavy target

Let’s start by considering the simple where the beam and target particles are
rigid spheres, and the target is infinitely heavy so undergoes no recoil. This is
like a billiard ball hitting another billiard ball that is affixed to the table. We
assume the billiard balls have radii R1 for the “beam” billiard ball, and R2
for the (fixed) target ball. It is hopefully obvious that the total cross-section
in this case should be π(R1 + R2 )2 , since that is the area over which the two
balls will hit each other. But to build intuition, let’s determine this using
differential cross-sections by finding σd (θ) and then integrating to get σT .

We define α as the angle to the normal at the point of contact between

the two spheres, relative to the beam axis (i.e. the direction of the incoming
ball). In the case of an infinitely heavy target and elastic scattering, the
angle of incidence must be the angle of reflection, so m1 reflects at the same
angle α. We thus get the scattering angle to be (as before) θ = π − 2α.

The central force here acts instantaneously while the balls are in contact.
It is still a repulsive radial force – the force is in the direction from the cen-
tre of the target to the centre of beam ball. The angle between radial force
and the beam axis is α. At the exact time of contact, the ball centres are
a distance R1 + R2 apart, at an angle α relative to the beam. We can thus
construct a right triangle with angle α, having R1 + R2 as the hypotenuse
and the impact parameter b as the opposite side. This gives us a relation
between b and α = (π − θ)/2:
!
π θ θ
b = (R1 + R2 ) sin α = sin − = (R1 + R2 ) cos
2 2 2
This is the relationship between b and θ that we need to get σd (θ). Note that
since θ1 ∈ [0, π], so b ∈ [R1 + R2 , 0]; for larger impact parameters there is no

134
contact and the scattering probability is zero. Taking the derivative gives

db (R1 + R2 ) θ
=− sin
dθ 2 2
Recalling that sin θ = 2 sin 2θ cos 2θ , the relationship between differential cross-
section and b and θ is:
b db (R1 + R2 )2 θ θ (R1 + R2 )2
σd (θ) = − = sin cos =
sin θ dθ 4 sin 2θ cos 2θ 2 2 4

In this case, the differential cross-section is independent of θ. The total cross

section is thus:
Z π
(R1 + R2 )2
σT = 2πσd sin θdθ = 2π (2) = π(R1 + R2 )2
0 4
This is the expected result.

14.4 Two body scattering: LAB and CM frames

Figure 27: Two-body scattering: Mass m1 comes in and strikes mass m2 ,

causing m1 to scatter at an angle θ1 while mass m2 recoils at angle θ2 .

135
 A more realistic situation is that the second ball will recoil, as in billiards.
How do we solve this? Recall we discussed §12.4 how two-body orbits work
just like the one-body case if we work in the centre-of-mass (CM) frame, so
let’s try that here. We will refer to the frame of the billiards table as the
laboratory frame (LAB).

A useful way to approach the two-body problem is thus to solve the one-
body case in the CM frame, and then transform the solution into the LAB
frame. We will find that considering momentum conservation in the two
frames is helpful for this, so let’s relate the momenta of the balls in the CM
and LAB frames. Let’s define the CM frame velocity as ⃗u in the LAB.

It is conventional to denote CM frame quantities with an asterisk (∗ ). Recall

from §12.4 that the positions of the beam (1) and target (2) particles in the
CM frame are given by:
µ µ
⃗r1 ∗ = ⃗r ⃗r2 ∗ = − ⃗r
m1 m2
where µ = m1 m2 /(m1 + m2 ) is the reduced mass, and ⃗r = ⃗r1 − ⃗r2 is the
relative position of the two bodies which is the same in both the LAB and
CM frames. Differentiating w.r.t. time, we get analogous relations for ⃗v1∗
and ⃗v2∗ :
µ µ
⃗v1 ∗ = ⃗v ⃗v2 ∗ = − ⃗v
m1 m2
with ⃗v = ⃗v1 −⃗v2 in place of ⃗r. One can immediately see that m1⃗v1∗ = −m2⃗v2∗ .
In other words, the momenta of the two balls are equal and opposite in the
CM frame, and the total momentum in the frame is zero. Let’s call the beam
momentum in the CM frame as ⃗p∗ = m1⃗v1∗ , which means the target momen-
tum in the CM frame is −⃗p∗ = m2⃗v2∗ .

Because the CM frame is boosted w.r.t. LAB, the momentum of an ob-

ject changes when measured in LAB vs. CM: ⃗p = ⃗p∗ + m⃗u. So the beam
ball momentum is:
⃗p1 = m1⃗v1 = m1⃗v1∗ + m1 ⃗u.
The target ball meanwhile is stationary in the LAB frame, so

0 = m2⃗v2∗ + m2 ⃗u =⇒ ⃗u = −⃗v2∗ .

136
This makes sense; the CM frame is boosted by +⃗u thus the (LAB-stationary)
target must acquire a CM velocity −⃗u. Combining these, we have
m1 ∗
⃗p1 = ⃗p∗ − m1⃗v2∗ = ⃗p∗ + ⃗p .
m2
This relates the LAB beam momentum to the CM beam momentum.

14.5 Hard-body scattering of unequal masses

Armed with our LAB/CM formalism, let us consider the billiard ball problem
again, this time allowing the second ball to recoil. The “beam” ball has mass
m1 , and an initial velocity ⃗v1,i in the LAB frame. It strikes a ball will mass
m2 . that is initially stationary. m1 scatters out at an angle θ1 with a final
velocity ⃗v1,f , while the second ball travels off at an angle θ2 with a final
velocity ⃗v2,f . Note that this is a central force, since we can represent the
scattering force as an infinite potential barrier:

V (r) = 0 for r > (R1 + R2 ) V (r) = ∞ for r ≤ (R1 + R2 )

Again, it is clear that the total cross section should be π(R1 + R2 )2 ; we will
use our cross section formalism to confirm this. But the more interesting
question is: What is the scattering angle of the two balls?

Here the CM frame proves useful, because the two objects can be regarded as
approaching a stationary center of mass point, and scattering off that point
as if it’s an infinitely heavy mass. Thus in the CM frame, we can use much
of the one-body case’s maths that we’ve already solved. The tricky part is,
we want the scattering angles in the LAB frame, so we will have to revert
the angles back into the LAB frame.

In the CM frame, the two bodies approach each other with equal and oppo-
site momenta, and thus must scatter at an equal angle in order for the total
momentum to be conserved (and zero). We define this CM frame scattering
angle as θ∗ for the beam (and thus −θ∗ for the target). The same geometry
applies as in the heavy target case, with the balls striking at an angle α (this
is the same in either frame), except now in the CM frame:

π θ∗
θ∗ = π − 2α =⇒ α = −
2 2
137
Figure 28: Billiard ball problem in the CM frame. The momenta must be
equal and opposite, but the velocities can be different since m1 ̸= m2 . The
angle of scattering θ∗ in the CM frame is the same for both, but in the LAB
frame the scattering angles are in general different.

The cross section works as before, except in terms of CM angles:

θ∗ db (R1 + R2 ) θ∗
b = (R1 + R2 ) sin α = (R1 + R2 ) cos =⇒ =− sin
2 dθ 2 2
b db (R1 + R2 )2
σd (θ∗ ) = − =
sin θ∗ dθ∗ 4
The differential cross-section is independent of the scattering angle as be-
fore, so it doesn’t depend on any CM frame quantities. Hence the total
cross-section, obtained by integrating over all angles, still gives the geomet-
rical overlap π(R1 + R2 )2 . So far, this is all as in the one-body case.

Next we want to find the LAB scattering angles θ1 and θ2 , given θ∗ . A

good way to approach this problem is to think about the momenta in the
two frames. For the beam ball, initially p1x = mv1,i and p1y = 0, while the
target has ⃗p2 = 0. After the collision, they must conserve momentum in the
x and y directions, in each frame. Hence the tangent of the scattering angle
must be equal to py /px for each ball in each frame. So we need to find py
and px for the scattered balls.

In our setup, ⃗u = ux̂ (the CM vs. LAB velocity) only has a component
in the x direction. Hence the y-momentum is unchanged by transformation
from LAB to CM: p1y = p∗1y = m1 v1∗ sin θ∗ . In contrast, the x-momentum is

138
boosted: p1x = p∗1x +m1 u = m1 v1∗ cos θ∗ +m1 u. But we know that, in the CM
frame, the momenta are equal and opposite: m1 v1∗ = −m2 v2∗ . Also, we know
that the the CM velocity u = −v2∗ , since m2 was initially stationary. Hence
m1 u = −m1 v2∗ = m21 v1∗ /m2 . Thus the scattering angle in the LAB frame is
p1y m1 v1∗ sin θ∗ sin θ∗
tan θ1 = = =
p1x m1 v1∗ cos θ∗ + m21 v1∗ /m2 cos θ∗ + m1 /m2
This relates the LAB scattering angle θ1 to the CM scattering angle θ∗ , which
we can compute from α or b as in the 1-body case.

What about the recoil angle θ2 of the target particle? This proceeds much
like above, except that for m2 the scattering angle is −θ∗ , and the momentum
boost between the frames is m2 u = −m2 v2∗ :
p2y −m2 v2∗ sin θ∗ sin θ∗ θ∗
tan θ2 = = = = cot
p2x m2 v2∗ cos θ∗ − m2 v2∗ 1 − cos θ∗ 2
Interestingly, the recoil angle of the stationary ball is independent of mass.

For m1 = m2 , the relationship between the LAB and CM scattering angles

is (using some trig identities):
∗ ∗
sin θ∗ 2 cos θ2 sin θ2 θ∗ θ∗
tan θ1 = = ∗ = tan =⇒ θ 1 =
1 + cos θ∗ 1 + 2 cos2 θ2 − 1 2 2
θ∗ π θ∗ π θ∗
!
tan θ2 = cot = tan − =⇒ θ2 = −
2 2 2 2 2
Hence for m1 = m2 , the opening angle θ1 + θ2 = π/2. In other words for
same-mass spheres, the balls always scatter at 90◦ with respect to each other.
Those who have played billiards know this is approximately the case, since
the cue ball is almost (although not exactly) the same mass as the billiard
balls. Note that the LAB scattering angle θ is thus restricted to the range 0
to π/2 for identical spheres (i.e. it is impossible to back-scatter of an object
of the same mass), but the CM angle θ∗ spans from 0 to π.

As an exercise to the reader, one can show that the LAB opening angle
between the scattered objects is:
(m1 + m2 ) θ∗
tan(θ1 + θ2 ) = cot
(m1 − m2 ) 2

139
14.6 Rutherford Scattering
In the early 1900’s, the internal structure of the atom was unknown. Thom-
son’s model argued for a “’plum pudding” arrangement, where the positive
protons and negative electrons were intermingled. Bohr favoured an orbital
configuration, with protons concentrated in a nucleus orbited by the lighter
electrons. Thomson’s model was widely accepted. But Rutherford’s famous
scattering experiment demonstrated that the positively charged portion of
the atom is concentrated in a massive nucleus, as argued by Bohr. He did
so by shooting α-particles (i.e. helium nuclei) at gold foil and measuring
the distribution of scattering angles, and thus inferring that the scattering
cross-section was much smaller than the atom itself.

We will compute the cross-section using our orbit formalism. We assume

that the gold nucleus is sufficiently massive compared to the α particle that
it can be taken to be at rest. Since this is a repulsive 1/r2 (electrostatic) force
with an unbound orbit, the trajectory will be a (repulsive branch) hyperbola.

Figure 29: The hyperbolic orbit of a charged particle approaching a massive

target nucleus. The scattering angle is θ = π − 2α, where α is the angle of
the asymptote.

To determine the differential cross section σ(θ), we need b(θ). Recall that
the angle α is related to b and a simply by
b π θ
tan α = α= − =⇒ b = a cot(θ/2)
a 2 2
The initial conditions are set by the impact parameter b and the initial ve-
locity v∞ , which can also be expressed as an energy E = 21 mv∞
2
. However,

140
we can’t know the exact b for any give α particle, so we’d like to express σ(θ)
in terms E and the measurable scattering angle θ. We first replace a with E
as follows:
k k k θ
E= =⇒ a = =⇒ b = cot
2a 2E 2E 2
Here, k = qα qG /4πϵ0 , where qα = +2e is the charge of the α particles, qG =
+79e is the charge of the gold nucleus, and ϵ0 is the vacuum permittivity.
From this, we can differentiate to get:
db k θ
=− csc2
dθ 4E 2
which gives us our differential cross section:
!2
b db 1 k cot 2θ csc2 2θ
σd (θ) = − =
sin θ dθ 2 2E 2 sin 2θ cos 2θ

k2 0.766 × 10−28
σd (θ) = = m2
16E 2 sin4 θ/2 sin4 θ/2
where we have used E = 6.5 MeV as in Rutherford’s original experiment
(1912). The value in the numerator gives an estimate of the cross-section of
the nucleus; in this case, 0.766 barns, where 1 barn is defined as 10−28 m2 is
the unit of cross-section used by atomic physicists. This cross-section is much
smaller than that of an entire atom (∼ 108 barns), thus strongly favoring the
Bohr model.

Note that the full integral of σd (θ) gives σT = ∞, since it is strongly di-
vergent at θ = 0. This is not surprising since there is always some deflection
even at very large b. But the integral gives a finite result when integrated
from b = 0 → bmax (or θ = π → θmin ), where bmax is the extent of the beam.

Rutherford also saw that the scattering probability followed the dependence
on θ as expected from this cross-section formula. This is quite different than
what the Thomson model predicted. For instance, in the Thomson model,
the α particle would only ever interact with similar-mass particles, meaning
that it could never back-scatter. In the heavy nucleus case, back-scattering
was rare, but possible. Rutherford’s famous quote describing his result was:
“It was almost as incredible as if you fired a 15-inch shell at a piece of tissue
paper and it came back and hit you.” Later, quantum mechanics would show

141
that the situation is not so simplistic as the Bohr model, instead favoring
Schrödinger’s probabilistic model, but to this day the Bohr model remains a
useful effective theory for atomic orbitals.

14.7 Rutherford scattering from momentum change

As a final demonstration of using the orbit equation, let’s solve Rutherford
scattering in a different way, by considering the impulse (i.e. change of mo-
mentum) imparted by the central force during the hyperbolic orbit. Going
back to the geometry in our original repulsive-branch hyperbolic orbit as
shown in Figure 24), the particle comes in at an angle −α with |⃗p| = mv∞ ,
and leaves at an angle +α with the same |⃗p|. It’s straightforward to see that
the change of momentum must be entirely in the x direction, hence

∆px = 2mv∞ cos α.

This total change in momentum must be equal to the sum of the infinites-
imal changes in x momentum, integrated over the entire orbit. During a
⃗ = kr̂/r2 , the beam particle will change its
time interval dt, given a force F
momentum by:

⃗ = kdt r̂ = kdt (cos ϕx̂ + sin ϕŷ)

d⃗p = Fdt
r2 r2
where we use ϕ instead of θ as our polar coordinate in order to distinguish
it from the scattering angle θ. The integral of the x-component is
Z +∞
k
∆px = cos ϕ dt
−∞ r2

We need to change variables from t to ϕ. Recall that L = mr2 dϕ

dt
, while for
a hyperbolic orbit L = mv∞ b. Equating these, we get

dϕ r2
mr2 = mv∞ b =⇒ dt = dϕ
dt v∞ b
Fortuitously, the r2 in dt cancels the 1/r2 inside the integral, so we don’t
need to know r(ϕ). This leaves a simple form:
k Z +α 2k sin α
∆px = cos ϕdϕ =
bv∞ −α bv∞

142
The integration now runs from −α to +α, because these are the values of θ′
at t = ∓∞. We can now equate this to the ∆px we found from considering
the orbit:
2k sin α k
= 2mv∞ cos α =⇒ b = 2
tan α
bv∞ mv∞
2
Recall that mv∞ = 2E and α = π/2 − θ/2, therefore

k k
b= tan α = cot(θ/2)
2E 2E
Recall that this is the exact same b(θ) as we derived in the previous section.
From this, the derivation to get σ(θ) follows just as before. This illustrates
how using differential momentum change can be a useful way to solve an
orbit problem.

143
15 Rotating Frames
We usually assume that experiments on Earth are conducted in an inertial
frame. But, in actuality, it is not. The surface of the Earth is a rotating
frame, and hence we and everything on the surface of the Earth is constantly
being accelerated as the Earth rotates. Nonetheless, we would like to be
able to define quantities within the Earth’s rotating frame, and be able to
understand the dynamics of objects as measured in this frame (even when
the Earth’s rotation impacts the motion). In this chapter, we will cover how
to handle dynamics in a rotating frame. The net result is that the rotating
frame can be regarded as inertial, provided that one adds several “fictitious”
forces to the dynamics. This idea has wide applicability beyond the Earth,
such as the car turning on a banked road that we considered early in the
course. The main difference in comparison to the 2-D rotating case before,
where we introduced the idea of a (fictitious) centrifugal force, is that now
the situation is 3-dimensional.

15.1 Non-inertial frames

Let us begin by setting up some definitions for non-inertial frames. A non-
inertial frame S ′ is defined as the “rest frame” of an object m1 that is accel-
erating at A⃗ within an inertial frame S. If |A|⃗ → 0 we recover an inertial
frame with ⃗v = 0. For simplicity assume that A ⃗ is time-independent, i.e.
there is a constant acceleration between the frames.

Consider a mass m1 at rest in S ′ , and a mass m2 at rest in S. When viewed

in the inertial frame S, an object in S ′ will have an acceleration boosted by
⃗
A:
⃗ =⇒ ⃗a′ = ⃗a − A
⃗a = ⃗a′ + A ⃗

Integrating w.r.t. time gives

⃗ t + ⃗u(0),
⃗v ′ = ⃗v − A

where ⃗u(0) = ⃗v′ (0) − ⃗v(0) is the relative velocity between the frames S and
S ′ at t = 0. Note that the velocity between frames ⃗u is time-dependent:
⃗ t
⃗u(t) = ⃗v ′ (t) − ⃗v(t) = u(0) − A

144
An observer in S attributes thus the acceleration of m1 to a real force, asso-
ciated with whatever is causing frame S ′ to rotate:
⃗ = m1 A
F ⃗ (real)

In contrast, an observer in S ′ views m2 (at rest in S) as accelerating, with

⃗ Thus the S ′ observer, in order to explain the motion of
an acceleration −A.
m2 , must introduce a fictitious force:
⃗ ′ = −m2 A
F ⃗ (fictitious)

Note that F⃗ ̸= −F ⃗ ′ (unless m1 = m2 ); it is the real and apparent accelera-

tions, not forces, that are equal and opposite.

This discussion can be generalised to the case where both objects have ac-
celerations ⃗a1 and ⃗a2 in an inertial frame, i.e. both S and S ′ are non-inertial
frames. In this case the two observers in S ′ and S attribute forces acting on
m2 and m1 respectively, which are partly real and partly fictitious:
⃗ ′ = m2 (⃗a2 − ⃗a1 )
F ⃗ = m1 (⃗a1 − ⃗a2 )
F

⃗ = ⃗a1 − ⃗a2 .
which can also be seen by defining A

For gravity, the special case of A ⃗ = −⃗g is known as free fall. A person
inside a lift whose cable breaks experiences this, as the frame is now acceler-
ating downwards at g. The observer thus feels a fictitious force that exactly
counteracts gravity, and there is no net force; the observer is weightless.
A satellite stably orbiting Earth is in the same situation – an observer on
the satellite experiences a centrifugal acceleration that balances gravity, and
hence is weightless. This is a fictitious force seen because the observer is in
the non-inertial frame of the orbiting satellite.

15.2 Constant angular velocity

A useful case is a rotating frame with a constant angular speed ω, such as
the Earth. An object at “rest” within this rotating frame experiences a cen-
tripetal acceleration with a magnitude of rω 2 r̂ owing to gravity. But the
acceleration direction is not constant, as r̂ varies with time. Thus the rotat-
ing frame will experience fictitious forces. You might guess that there will

145
be a centrifugal force owing to the rotation. Here we will show that, in this
3-D case, there is another fictitious force called the coriolis force, which is
actually the dominant fictitious force in most circumstances.

Let’s call the rotating frame S ′ and the inertial frame S, as before. The
S frame is often called the “fixed stars” frame. Consider an arbitrary time-
⃗
dependent vector B(t). In S, we can represent it via its Cartesian compo-
nents, as usual:
⃗ = Bx î + By ĵ + Bz k̂
B
⃗ in the rotating S ′ frame:
We can similarly write B
⃗ = B ′ î′ + B ′ ĵ′ + B ′ k̂′
B x y z

where here we introduce a basis vector set (î′ , ĵ′ , k̂′ ) that is rotating with the
S ′ frame. These S ′ basis vectors change with time! For instance, k̂′ might
be defined as the direction vertical to the Earth’s surface, but relative to
the “fixed stars” frame, this basis vector rotates with the Earth, and so is
not constant. This is analogous to the circular motion case where the polar
coordinate basis vectors were time-dependent.

⃗
Now, dB/dt should be the same in the two frames. In the S frame, the
derivative is straightforward, since the basis vectors are constant:
 
⃗
dB dBx dBy dBz
  = î + ĵ + k̂.
dt dt dt dt
S

The subscript here denotes the frame of the basis vectors against which the
derivative is computed.

In the S ′ frame, the derivative will need to include terms that account for
the time variation of the S ′ basis vectors:
 
⃗
dB dBx′ ′ dBy′ ′ dBz′ ′ dî ′ dĵ ′ dk̂ ′
  = î + ĵ + k̂ + Bx′ + By′ + Bz′
dt dt dt dt dt dt dt
S

The first three terms on the RHS are just [dB ⃗ ′ /dt]S ′ , i.e. what an observer
⃗ But this observer would
in S ′ would naturally regard as the derivative of B.
see that the derivative doesn’t fully account for how B ⃗ changes. To account

146
for those last three terms, we need to determine what the derivatives of the
S ′ basis vectors are.

Consider dĵ′ , where ĵ′ is the east-west direction (mostly because it’s the
easiest one to illustrate). Over some time dt, ĵ′ changes by an angle of dθ,
in a direction towards the axis of rotation. Hence dĵ′ /dt = dθ/dt = ω in
magnitude, and the direction is given by ω̂ × ĵ′ . One can similarly show that
this is true for each of the basis vectors, so:
dî ′ dĵ ′ dk̂ ′
⃗ × î ′
=ω ⃗ × ĵ ′
=ω ⃗ × k̂ ′
=ω
dt dt dt
Thus
dĵ ′
By′ = By′ ω
⃗ × ĵ ′ = ω
⃗ × By ĵ ′
dt
and similarly for the î′ (north-south) and k̂′ (up-down) terms. Combining
⃗ ×B
these gives ω ⃗ for all the basis vector derivative terms.

⃗ in the
Putting this together, we can now relate the time derivative of B
′
S frame to quantities in the S frame:
   
⃗
dB ⃗
dB
  =  +ω ⃗
⃗ ×B
dt dt ′
S S

This tells us how to take a time derivative of a quantity in a rotating frame:

we take the derivative w.r.t. the S ′ basis vectors, then add an extra term
⃗ ×B
ω ⃗ to account for the frame’s rotation.

15.3 Coriolis and Centrifugal Forces

With the above generic formula, we can determine the fictitious forces intro-
⃗ = ⃗r. Hence
duced by the rotation of the S ′ frame. To do this, we first set B
⃗
dB/dt = ⃗vS . The formula then lets us relate the velocities in the S and S ′
frames:
⃗ × ⃗r
⃗vS = ⃗vS ′ + ω
To get the acceleration, we must take another time derivative. So we now
⃗ = ⃗vS , and hence dB/dt
set B ⃗ = ⃗aS :
" # " #
d⃗vS d⃗vS
⃗aS = = ⃗ × ⃗vS
+ω
dt S
dt S′

147
Replacing ⃗vS by ⃗vS ′ on the RHS:
" #
⃗ × ⃗r)
d(⃗vS ′ + ω
⃗aS = ⃗ × (⃗vS ′ + ω
+ω ⃗ × ⃗r)
dt S′

⃗aS = ⃗aS ′ + 2⃗ω × ⃗vS ′ + ω

⃗ × (⃗ω × ⃗r)
This shows that the rotating frame has introduced two new terms, each of
which can be interpreted as a fictitious force. Here we explain all the terms
in this equation in more detail:

• ⃗aS is the acceleration of the object in the inertial frame S. It corre-

⃗ = m⃗aS acting on the object.
sponds to a real force F

• ⃗aS ′ = d⃗vS ′ /dt is the acceleration of the object in the rotating frame
S ′ , which is a combination of its inertial-frame acceleration plus the
frame’s acceleration.

• 2⃗ω × ⃗vS ′ represents the fictitious coriolis acceleration in the S ′ frame.

It acts perpendicular to the rotation axis, and to the velocity of the
object in the S ′ frame. It is zero for objects at rest in S ′ or moving
parallel to the rotation axis. It has a magnitude ∝ ω.

• ω
⃗ ×(⃗ω ×⃗r) is recognisable as the fictitious centrifugal acceleration. This
force acts radially outwards so as to compensate for the real centripetal
acceleration. It depends on the radial position of the object, but not its
motion. One can see that the cross products pick out the component
of ⃗r that is perpendicular to the rotation axis, so this is equivalent
to treating the Earth as a giant (flattened) merry-go-round. It has a
magnitude ∝ ω 2 .

We can thus write a modified version of N2 which describes the motion in

the rotating S ′ frame, purely in rotating-frame quantities:
⃗ − 2m⃗ω × ⃗vS ′ − m⃗ω × (⃗ω × ⃗r)
m⃗aS ′ = F

Remember that with ⃗aS ′ on the LHS and ⃗aS on the RHS, the fictional forces
appear with negative signs. With ω⃗ = 0, we recover the usual N2. But to
use N2 in the rotating frame, we must add the coriolis and centrifugal forces.

148
15.4 The Earth as a Rotating Frame
The Earth is an example of a rotating frame with a constant angular velocity:
2π
ω
⃗E = k̂ = 7.3 × 10−5 k̂ rad/s
24 × 60 × 60

where we define k̂ as the unit vector towards the North pole from the Earth’s
centre. Since ωE ≪ 1, the coriolis force (first order in ω) usually dominates
over the centrifugal force (second order in ω). Note that a second-order ac-
curate solution can be obtained by making a perturbative correction to the
value of g, analogous to how we did precession; we will not consider this here.

The coriolis force can have a significant effect on objects moving relative to
the Earth’s surface (⃗vS ′ ̸= 0), particularly over long distances. Consider
ocean currents. As the Earth turns, the equator moves faster than the poles,
so currents moving south will find its angular velocity slower than the land,
and hence will be deflected towards the west (see image) in the Earth’s frame.
Conversely, currents moving north will find their angular velocity faster than
the land, and so will deflect towards the east in the Earth’s frame. This is the
effect which gives rise to the Gulf Stream circulation pattern in the Atlantic,
in which water coming northwards from the warm Carribean and is pushed
eastwards towards Scotland and Europe. This gives Scotland much milder
temperatures than comparable latitudes in e.g. Canada or Russia.

149
Example: Consider a skydiver jumping out of a plane. For simplicity, we
assume zero initial velocity in the Earth’s frame, and we ignore air resistance.
Because of the coriolis force, the skydiver in free fall towards the Earth does
not fall vertically, but rather ends up is displaced from the vertical position
at the jump. How much is that displacement?

We define x′ as North (along the Earth’s surface), y ′ as East, and z ′ as

the height above the Earth in the rotating frame. Since the Earth rotates
to the East, the displacement will be in the y ′ direction. We thus need to
consider the ĵ′ component of the Coriolis force:
d2 y ′ ′
m ĵ = −2m⃗ω × ⃗vS ′
dt2
For the RHS, the ĵ′ of the cross product is ωx′ vz′ − ωz′ vx′ . But we assume
that, in the S ′ frame, the skydiver falls only vertically (in the z ′ direction),
hence vx′ = 0. Meanwhile, a bit of geometry shows that ωx′ = ωE cos ϕ where
′
ϕ is the latitude. We write vz′ = dz
dt
, thus:

d2 y ′ ′ dz ′
m ĵ = −2mωE cos ϕ
dt2 dt
dz ′
To determine dt
, we note that the acceleration of gravity is −gẑ ′ , so

dz ′ d2 y ′
= −gt =⇒ = 2ωE gt cos ϕ
dt dt2
where we have applied the initial condition that the vertical velocity is zero
at t = 0. Integrating twice gives the displacement to the East as a function
of t:
1
y ′ = ωE gt3 cos ϕ
3
where for simplicity we have assumed y ′ = dy ′ /dt = 0 at t = 0. Hence the
skydiver moves in the y ′ (East) direction by an amount that depends on the
cube of the fall time. At the North Pole (ϕ = 90◦ ), the skydiver experiences
no Coriolis displacement, since they are falling along the axis of the Earth’s
rotation.

Extra: If we can’t assume ω

⃗ is small, then we must write down the equation

150
of motion for the velocity of a falling object under gravity, in the rotating
frame as follows:
d⃗v′
!
m + 2ω⃗E × ⃗v = −mgẑ′
′
dt
2 ′
In this fully general case, besides the ddty2 equation above, we additionally
have equations including the Coriolis force in the North/South (x′ ) and
up/down (z ′ ) direction, which must be solved as a coupled set of ODEs:

d2 x′ dy ′ d2 z ′ dy ′
m = 2mωE sin ϕ m = −mg + 2mωE cos ϕ
dt2 dt dt2 dt

15.5 Foucault’s Pendulum

In 1851, Léon Foucault had an ingeneous idea to prove that the Earth rotated
relative to the fixed stars, by using the coriolis force. Let us solve Foucault’s
pendulum to see how he did it. We set up the acceleration equation for a
simple pendulum as we did before, but now we include the coriolis force (but
neglecting centrifugal):

d⃗v′
!
m ⃗ = −mgẑ′ − Tr̂
+ 2⃗ωE × ⃗v′ = F
dt

Previously, we have solved this problem in an inertial frame. In the small

angle approximation, we were able to reduce this to a simple harmonic oscil-
2
lator equation ddt2θ + ω02 θ = 0, giving the usual oscillatory solution having a
frequency ω02 = g/L.

In the rotating frame, we now have an extra term corresponding to cori-

olis force. This provides an acceleration perpendicular to the plane of the
swinging. Recall that a perpendicular force results in circular motion. So in
addition to the usual back and forth oscillation, we expect that the plane of
oscillation of the pendulum will actually spin around, owind to the coriolis
force! Let us calculate this rotation rate.

As before, we define y ′ to be east, and z ′ to be normal to the ground (and

thus x′ is to the south), all in the rotating (Earth) frame. The cross product
⃗ × ⃗v′ is:
ω

⃗ × ⃗v′ = [ωy′ vz′ − ωz′ vy′ ]î′ + [ωz′ vx′ − ωx′ vz′ ]ĵ′ + [ωx′ vy′ − ωy′ vx′ ]k̂′
ω

151
For small displacements, the vertical motion is small (vz′ ≈ 0), so we can
ignore vz′ terms. This also means the k̂′ component is irrelevant, since the
coriolis term will not add any motion in that direction. So this reduces to:
⃗ × ⃗v′ = (−ωz′ vy′ )î′ + (ωz′ vx′ )ĵ′
ω
Now, ωz′ is the component of the Earth’s rotation in the vertical direction
at latitude ϕ, which simple geometry shows is ωz′ = ωE sin ϕ. Hence:
dy ′ ′ dx′ ′
⃗ E × ⃗v′ = −ωE sin ϕ
ω î + ωE sin ϕ ĵ
dt dt
Thus, our vector force equation including the Coriolis term gives us two
coupled equations for the î′ and ĵ′ directions:
d2 x′ dy ′ d2 y ′ dx′
− 2ωE sin ϕ + ω02 x′ = 0 + 2ω E sin ϕ + ω02 y ′ = 0
dt2 dt dt2 dt

We will guess a trial solution for x′ (t) and y ′ (t). We note that to satisfy the
above equations, a single derivative of x′ must have the same functional form
as y ′ or its double derivative, and vice versa. It is evident that that x′ as a
cosine and y ′ as a sine would satisfy this. This suggests a trial solution of
the form:
x′ = cos αt y ′ = sin αt
where we choose a simple phase so that the pendulum begins swinging in
the N-S direction (along î′ ); we’ll consider the amplitude later. Note that we
could have equivalently used η = eiαt = x′ + iy ′ , and considered the real and
imaginary parts. For our trial solution,
dx′ dy ′
= −α sin αt = α cos αt
dt dt
d 2 x′ 2 d2 y ′
= −α cos αt = −α2 sin αt
dt2 dt2
Substituting back in to the x′ equation above, we get
−α2 − 2αωE sin ϕ + ω02 = 0
(you can check that the y ′ equation gives exactly the same relation.) This is
a quadratic equation for our unknown α, which has the solution:
q
α = −ωE sin ϕ ± (ωE sin ϕ)2 + ω02

152
We define ω and Ω as follows, and get a simple form for α:
ω ≡ −ωE sin ϕ, Ω2 ≡ ω02 + ω 2 =⇒ α = ω ± Ω
Hence the full solution is
x′ = C1 cos(ω + Ω)t + C2 cos(ω − Ω)t y ′ = C1 sin(ω + Ω)t + C2 sin(ω − Ω)t
Using the trig identities cos(A + B) = cos A cos B − sin A sin B, and sin(A +
B) = sin A cos B − cos A sin B, we have
x′ = C1 [cos Ωt cos ωt − sin Ωt sin ωt] + C2 [cos Ωt cos ωt + sin Ωt sin ωt]
y ′ = C1 [sin Ωt cos ωt + cos Ωt sin ωt] + C2 [− sin Ωt cos ωt + cos Ωt sin ωt]
We pick a set of initial conditions such that C1 = C2 = 12 C, which simplifies
these equations to:
x′ = C cos ωt cos Ωt y ′ = C sin ωt cos Ωt
where ω and Ω are defined above. Note that if we set ωE = 0 (so ω = 0), we
recover the inertial frame SHO solution with Ω = ω0 , i.e. x′ = cos ω0 t, y ′ = 0.

Foucault’s pendulum is thus a superposition of two harmonic oscillators:

1. An oscillation with a frequency Ω corresponding to the pendulum bob
swinging back and forth at slightly above the natural frequency ω0 ;
2. A precession rotating the oscillation in the (x′ , y ′ ) plane (i.e. the surface
of the Earth) with a frequency ω that depends on the latitude.
The direction of the precession (clockwise or counter-clockwise) is determined
by whether you are in the northern (ϕ > 0, so ω < 0) or southern (ϕ < 0, so
ω > 0) hemisphere. In the north, at t = 0 we have x′ = C and y ′ = 0, and
as t increases y ′ goes negative since ω < 0; this means clockwise rotation.
In the south, ω > 0, hence we get counter-clockwise motion. At the equator
(ϕ = 0), ω = 0, so there is no precession. At Edinburgh (latitude ϕ = 56◦ N),
the rate of rotation is 298◦ per day. At the Paris Observatory where Foucault
originally set up his pendulum, the latitude is 49◦ , and a full rotation occurs
in about 32 hours. Such precession is easily noticeable, thus demonstrating
that the Earth rotates and measuring ωE .

See [Link] [Link] for

further description and cool animations.

153
16 Rigid Bodies
So far we have generally treated objects as point masses. But real objects
have a finite extent and a fixed shape that they maintain as a rigid body. For
the purpose of discussing their dynamics, rigid bodies are best described via
moments of inertia I about the principal symmetry axes of the body.

There are three additional degrees of freedom in the dynamics of rigid bodies
associated with rotational motion about the three principal symmetry axes.
These are orthogonal and pass through the Centre of Mass of the body.

16.1 Moments of Inertia

The moment of inertia of a rigid body about one of its principal axes A is
given by: Z
IA = |⃗rA |2 ρ(⃗r)dV
where the distance ⃗rA is taken perpendicular to A.

Figure 30: Moments of inertia of some common shapes.

Moments of inertia for some shapes are shown in Figure 30 above. Here is
how to derive the first two examples:
• A hoopRabout the symmetry axis through its centre: Here, rA = R, so
I = R2 ρdV = M R2 (Note: ρ(⃗r) does not have to be constant here).

154
 • A solid cylinder about its symmetry axis, with ρ =constant, height H,
radius R, and M = ρV = πR2 Hρ:
Z R
2 R4 1
I = Hρ rA 2πrA drA = (2πHρ) = M R2
0 4 2

The numerical factor in front of each of these can be absorbed into the
definition of a radius of gyration, kA , for each object:
1 Z
kA2 = |⃗rA |2 ρ(⃗r)dV IA = M kA2
M
This represents the distance from the axis where all the mass would need to
be put (like a hoop) to give the same moment of inertia.

16.1.1 Parallel and Perpendicular Axes

The parallel axis theorem states that:

IA′ = IA + M d2

where d = |⃗d| is the distance of an axis A′ from A, assuming A passes through

the centre of mass. This can be proved by changing ⃗rA → ⃗rA + ⃗d:
Z Z Z Z
IA′ = |⃗rA + ⃗d|2 ρ(⃗r)dV = |⃗rA |2 ρ(⃗r)dV + 2⃗d · ⃗rA ρ(⃗r)dV + |⃗d|2 ρ(⃗r)dV

The first term is IA . The second term vanishes since A passes through the
CoM (one can show this from the definition of the center of mass). The third
term is trivially M d2 . This proves the parallel axis theorem.

In Figure 30, the rightmost lower case (a rod about an axis at one end) can
be seen from taking the rightmost upper case (a rod about an axis through
the middle) and noting that the d = L/2. Hence I in the end-axis case is
1
12
M L2 + 14 M L2 = 31 M L2 .

The perpendicular axes theorem states that for a planar (2-D) body, the
moment of inertia along an axis perpendicular to the plane is the sum of the
moments about two perpendicular axes in the plane (call them A and A′ ).

I⊥ = Ix + Iy

155
This can be proved as follows. Let x be the distance to A and y be the
distance to A′ . Then the perpendicular distance to the third (z) axis is
|⃗rA |2 = x2 + y 2 , hence
Z Z Z
2 2 2
Itot = (x + y )ρdV = x ρdV + y 2 ρdV = IA + IA′

For instance, in Figure 30, the moment of inertia diametrically through a

hoop is 21 M L2 , so summing these in the two x and y directions gives M L2 ,
which is the moment of inertia through the perpendicular axis.

16.1.2 The Inertia Tensor

For complicated and non-uniform shapes the moments of inertia can be gen-
eralised in the form of a 3 × 3 inertia tensor:
   
Z
x2k  − xi xj  ρdV
X
Iij = δij 
k=1,3

Setting i = j recovers Ix , Iy and Iz as the diagonal elements of the tensor:

Z Z Z
I11 = Ix = (y 2 +z 2 )ρdV I22 = Iy = (x2 +z 2 )ρdV I33 = Iz = (x2 +y 2 )ρdV

The off-diagonal elements are given by:

Z Z Z
I12 = I21 = − (xy)ρdV I23 = I32 = − (yz)ρdV I31 = I13 = − (xz)ρdV

The inertia tensor is a symmetric matrix that can be diagonalised in the

standard way by finding the eigenvalues and eigenvectors, and constructing
a transformation matrix M:
I′ = M · I · MT
The matrix M represents a rotation in the 3D space between the original
cartesian coordinates and the coordinate system of the principal axes.

16.2 Rolling Motion

An object being a rigid body only comes into play when there is some rotation
about an axis, since in the case where the rigid body doesn’t rotate, it can
be treated as a point mass. Recall that the angular momentum about A is:
⃗ = IA ω
L ⃗A

156
where L and ω point along the direction of A.

Spherical and cylindrical objects can roll across surfaces. For an object of
radius R, in order to be rolling, the distance travelled as the object rotates by
∆θ must be given by ∆x = R∆θ. Dividing by ∆t, we get v = Rω. Expressed
more generally, the roll condition is:

⃗v = R(⃗ω × k̂)

since the cross-product with k̂ (normal to the surface) gives the direction of
motion.

Hence in order to roll, the object must separately satisfy two dynamics equa-
tions: one for the linear velocity ⃗v set by the linear force on the object, and
another for the angular velocity ω ⃗ set by the torque on the object. If the
force and torque are such that the roll condition is not satisfied, the object
will slide in addition to roll.

Example: Consider a pool ball at rest on a table. If you strike it with

the cue through its centre, then initially there is a force but no torque. It
is evident that the roll condition cannot be satisfied immediately because ⃗v
is nonzero but ω⃗ = 0, so the ball will slide across the table. If you want it
to roll immediately you need to hit it above the centre, such that your cue
provides both a force and a torque in exactly the right proportion to satisfy
the roll equation. On a frictionless table, what impact parameter b above
the centre of the ball do you need to strike the ball in order to get pure rolling?

The linear force is

dv dω
F =M = MR
dt dt
using the roll condition v = Rω.

Meanwhile, the torque is τ = ⃗r × F. ⃗ The cross product picks out the per-
pendicular distance at which the force is applied relative to the ball centre,
which is b. Hence τ = bF = bM R dω dt
. We can equate this to the torque
dω 2 2
τ = I dt . For I = 5 M R , this gives:

dω dω 2 dω 2
bM R =I = M R2 =⇒ b = R
dt dt 5 dt 5

157
Hence in order to ensure pure rolling, a horizontal force must be applied
at a distance of 0.4R above the centre of the ball. You can convince your-
self that for a cylinder, this would be b = 12 R, and for a hollow sphere b = 23 R.

Example: Let’s now consider a more realistic situation, where the pool table
has friction. Say you strike the ball in the middle. At first, it will slide, i.e.
purely translation motion. Friction will cause it to slow down, and also will
apply a torque to the ball since the friction is acting at the point of contant
with the table. Hence eventually, the roll condition will be satisfied, and ball
will start rolling. At what velocity, relative to the initial velocity V , does
this rolling happen?

For simplicity we assume that, until the roll condition is satisfied, the ball’s
motion is purely translational. (In reality, the ball will already slowly start
turning until it can fully roll.) After the ball starts moving, it experiences a
⃗ = −Ffrict î. The linear force equation is thus:
frictional force F
d⃗v
M = −Ffrict î
dt
As before, we substitute this into the torque equation τ = ⃗r × F. ⃗ The
frictional force is applied at a radius ⃗r = −Rk̂, where the ball meets the
table. Hence
2 d⃗ω ⃗ = (−Rk̂) × (−Ffrict⃗i)
τ = M R2 = ⃗r × F
5 dt
So we have a translation equation for ⃗v, and a rotational equation for ω
⃗ . To
solve these simultaneously, we can cross-product the translational equation
with −Rk̂, which you will see is exactly equal to the RHS of rotational
equation. Hence !
d⃗v 2 d⃗ω
−M R k̂ × = M R2
dt 5 dt
We now want to apply the roll condition. Let us first think about the di-
rections. The velocity is in the x-direction, ⃗v = vx î. The rotation is in the
⃗ = ωy ĵ. Since k̂ × î = ĵ, the result of both sides is thus in the ĵ
y-direction, ω
direction, which is good. Substituting for ⃗v and ω ⃗ (and cancelling M R), we
have:
dvx 2R dωy
=−
dt 5 dt
158
We can integrate both sides with time. On the LHS, we get vx − vx (t = 0) =
vx − V . On the RHS, the integral of the differential gives ωy − ωy (t = 0) = ωy
since there is no initial rotation. Thus:
2R
vx − V = − ωy
5
For the roll condition to be satisfied, we need vx = Rωy . Solving for vx gives

2 5
vx − V = − vx =⇒ vx,roll = V
5 7
So if you strike a ball in the middle, it loses 2/7 of its initial speed before it
starts to roll. Notice that this is independent of frictional force Ffrict , or the
radius R! Try this the next time you are playing billiards.

159
17 Lagrangian Dynamics
As you probably noticed over this course, much of the challenge in a solv-
ing dynamics problem is in correctly setting up the force equations. So
far we have most commonly used N2, and sometimes energy conservation
(dE/dt = 0) or total momentum conservation (d⃗ptot /dt = 0), to obtain
the equations of motion. Lagrangian dynamics is an alternative approach
to doing this that elegantly avoids some of the difficulties with these other
methods. The idea is to simplify obtaining the equations of motion by only
requiring one to write down the energies (potential and kinetic) of a system,
rather than the forces, and then obtaining the equations of motion from the
ansatz that the dynamical system satistfies the Principle of Least Action,
which we discuss below. It turns out, this reformulation of classical mechan-
ics has wide applicability beyond dynamics, including in quantum mechanics
and electromagnetism.

17.1 The Euler-Lagrange Equation

We begin with a derivation of the Euler-Lagrange equation, which is the
mathematics that underlies Lagrangian dynamics. Consider a function f (q, q̇),
where q(t) is some generalized coordinate, and q̇ = dq/dt is like a generalized
velocity (or momentum). The action S is defined an integral of f between
times t1 and t2 : Z t2
S= f (q, q̇)dt
t1
The Principle of Least Action states that for any dynamical system, the path
taken by f from t1 → t2 will minimize S. Historically, this was an outgrowth
of Fermat’s Principle of Least Time, which states that light takes the shortest
possible time between two points.

To formulate this mathematically, let us assume that the qtrue (t) is the path
that minimises S. Now consider a deviation from this path given by ϵη(t),
where ϵ is a small time-independent number. The deviated path is then given
by
dη
q(t) = qtrue (t) + ϵη(t) q̇(t) = q̇true (t) + ϵ
dt
Note that η(t1 ) = η(t2 ) = 0, since the true and deviated paths must both
begin and end at the same place.

160
The Principle of Least Action implies that S must satisfy the condition
dS/dϵ = 0. In detail, this only requires that S be extremised, not neces-
sarily minimised, and the mathematics only requires an extremum; however,
in most realistic situations, S will be minimised not maximised by this con-
dition. To take the derivative of S, we sum the partial derivatives of f w.r.t.
q and q̇:
dS Z t2  ∂f dq ∂f dq̇  Z t2 
∂f ∂f dη 
= + dt = η+ dt
dϵ t1 ∂q dϵ ∂ q̇ dϵ t1 ∂q ∂ q̇ dt
The second term inside the integral can be rewritten via the product rule:
! !
d ∂f ∂f dη d ∂f
η = + η
dt ∂ q̇ ∂ q̇ dt dt ∂ q̇
Thus we can write the second term as :
Z t2 Z t2 Z t2 h ∂f it2 Z t2 d  ∂f 
∂f dη d  ∂f  d  ∂f 
dt = η dt − ηdt = η − ηdt
t1 ∂ q̇ dt t1 dt ∂ q̇ t1 dt ∂ q̇ ∂ q̇ t1 t1 dt ∂ q̇

∂f
Now, η(t1 ) = η(t2 ) = 0 since at the endpoints, q(t) = qtrue (t). Hence the ∂ q̇
η
term above evaluates to zero. This leaves
dS Z t2  ∂f d ∂f 
= − ηdt
dϵ t1 ∂q dt ∂ q̇

The Principle of Least Action says that the above must vanish, i.e. dS
dϵ
= 0.
For this to be true for any arbitrary η(t), the integrand must be identically
zero. Hence !
∂f d ∂f
=
∂q dt ∂ q̇
This is known as the Euler-Lagrange equation. This is true for any function
f (q, q̇) that is only a function of a coordinate and its time derivative, and
not for instance an explicit function of time, that satisfies the Principle of
Least Action.

17.2 The Lagrangian

The Euler-Lagrange equation itself has nothing to do with dynamics; we
have simply used a variational principle to determine a condition on f that

161
minimizes the action, given that f is a function of some variable p and its
time derivative. Lagrange (1788) realised, however, that a particular choice
of f leads to dynamical equations of motion!

The function that turns out to be useful is what we now call the Lagrangian:
L=T −V
where T is the kinetic energy and V the potential energy of the system. Con-
trast this with total energy E = T + V . Note that the forces in the problem
must be describable by a potential V (e.g. conservative forces).

We can identify q as coordinate positions, and q̇ as the corresponding ve-

locities. The chosen coordinates can be Cartesian, polar, spherical, or even
some generalised coordinate that conveniently describes the system. The
Euler-Lagrange equation must be satisfied for each coordinate separately.
Thus in general, we can write the Euler-Lagrange equations for L for each
coordinate qi as: !
∂L d ∂L
=
∂qi dt ∂ q̇i
For a single body problem in N dimensions there are N equations, one for
each of the generalised coordinates {qi }. For a second body, there will be
another set of N coordinates describing that body, which will have another
N ODEs, and so on.

{qi } can be any set of coordinates that describe the system. It is often
convenient to choose {qi } as the dependent variables of interest in the sys-
tem, even if they are not orthogonal. Then one obtains equations of motion
for each of the dependent variables w.r.t. time t, which is what one wants.

17.3 Relationship to Newtonian mechanics

We now demonstrate that the Lagrangian formulation recovers the same
equations of motion as the more familiar Newtonian mechanics.

Linear motion in 3D: Here, the qi are the Cartesian (x, y, z), with q̇i =
(vx , vy , vz ), and the Lagrangian is given by
1
L = m(vx2 + vy2 + vz2 ) − V (x, y, z)
2
162
We now substitute this into the Euler-Lagrange equation. Remember, in a
partial derivative, we compute the derivative only of terms that explicitly
depend on the variable. So ∂L/∂x = −dV /dx = Fx , even though vx may
have some implicit dependence on x. Similarly, ∂L/∂vx = mvx . Putting
these into the Euler-Lagrange equation, we get
d
Fx = (mvx ) = max
dt
(for constant m). The y, z coordinates yield analogous equations. Hence in
this case, the Euler-Lagrange equations reduce to N2 in Cartesian coordi-
nates!

Rotational motion: We choose qi as the polar coordinates (r,θ), for which

the q̇i are the radial and tangential velocities vr and rω. Hence
m 2 2 
L= r ω + vr2 − V (r, θ)
2
Taking the partial derivatives, ∂L/∂r = mrω 2 −dV /dr, while ∂L/∂vr = mvr .
Hence the Euler-Lagrange equation for r is:
dV dvr dvr
mrω 2 − =m =⇒ Fr = m − mrω 2
dr dt dt
You should recognise this as the radial force in a rotating system including
the centrifugal acceleration term −rω 2 .

The second coordinate is qi = θ, and q̇i = ω. Here, ∂L/∂θ = −dV /dθ

and ∂L/∂ω = mr2 ω, giving:

dω dV
mr2 =− = rFθ
dt dθ
since Fθ = −dV /d(rθ), i.e. the derivative of the potential w.r.t. a tangential
line element at a radius r. This equation should be familiar as the torque
equation:
⃗ = I d⃗ω
⃗τ = ⃗r × F
dt
Hence by taking the coordinates (r, θ), we have derived both the centripetal
acceleration term of the radial force equation, and the torque equation for

163
the tangential component!

Simple pendulum: Let us now do an actual physical problem. In a simple

pendulum, there is only one coordinate, θ, with a corresponding velocity ω.
The pendulum bob’s velocity is Lω, while the height of the bob is L(1−cos θ).
Hence the Lagrangian for this system is:
1
L = mL2 ω 2 − mgL(1 − cos θ)
2
Taking the partial derivatives, we have:
∂L ∂L
= −mgL sin θ = mL2 ω
∂θ ∂ω
Putting these partial derivatives into the Euler-Lagrange equation gives the
familiar result:
dω d2 θ g
mL2 = −mgL sin θ =⇒ 2
≈− θ
dt dt L
in the small angle approximation.
q This should be familiar as SHM with a
natural frequency ω0 = g/L, as we obtained earlier in the course using N2.

These examples illustrate that, if one can write down the kinetic and poten-
tial energies for a system in terms of some general orthonormal coordinates,
it is straightforward to derive the equations of motion in those coordinates
from the Euler-Lagrange equation. Since the kinetic and potential energies
are scalars, these are often easier to compute that vector force balance in
each coordinate. It is rather remarkable that the simple choice of L = T − V
leads to this rich behavior; while not as physically intuitive as N2, it is often
mathematically quite a bit simpler, especially in more complex situations.

17.4 More complicated systems

We now consider some complex situations that can be challenging to solve
with traditional Newtonian dynamics, but are relatively straightforward via
Lagrangian dynamics.

164
17.4.1 A ball on a cylinder
We start with a a ball constrained to move on a cylindrical surface under
the influence of a central force F = −kr acting towards the symmetry centre
of the cylinder. There are two contributions to the kinetic energy coming
from motion parallel to the axis of the cylinder, z, and motion round the
axis in ϕ. The surface of the cylinder constraint prevents motion in the third
coordinate ρ.  !2 !2 
1 dϕ dz 
T = m R2 +
2 dt dt
The potential energy due to the central force is:
1 1
V = kr2 = k(R2 + z 2 )
2 2
Using the Euler-Lagrange equation with derivatives taken w.r.t. ϕ:

d2 ϕ dϕ
mR2 2
= 0 =⇒ L = mR2 = constant
dt dt
which is just angular momentum conservation in a central force.
Using the Euler-Lagrange equation with derivatives taken w.r.t. z:

d2 z
m = −kz
dt2
The ball can perform SHM in the z direction in combination with circular
motion in the ϕ direction.

17.4.2 Motion of ball on table and ball below table

We return to a problem from the orbits part of this course and solve it using
Lagrangian methods. A ball of mass m can move on a smooth horizontal
table. It is connected by an inextensible string of length ℓ to an identical
ball of mass m that hangs vertically below the table, with the string passing
through a frictionless hole in the table. We represent the position of the mass
on the table in polar coordinates (r, θ), and the hanging mass has a distance
z = ℓ − r below the table.

165
The kinetic energy thus has contributions from the polar coordinate mo-
tion of the ball on the table, and from the vertical motion of the ball below
the table:
 !2 !2   !2 
1 dr dz 1 dr
T = m r2 θ̇2 + +  = m r 2 θ̇ 2 + 2 
2 dt dt 2 dt

since z = ℓ − r hence dz/dt = −dr/dt. The gravitational potential energy

only comes from the lower ball:

V = −mgz = mg(r − ℓ)

We can use the derivatives of L = T − V w.r.t. r and ṙ:

∂L ∂L dr
= mθ̇2 r − mg = 2m
∂r ∂ ṙ dt
and using the Euler-Lagrange equation, we obtain the equation of motion:
d2 r
2 2 = θ̇2 r − g
dt
You may recall that this agrees with the result obtained using N2. Meanwhile,
the θ equation once again reduces to mr2 θ̇ =constant, i.e. conservation of
angular momentum.

17.4.3 A bead on a rotating circular wire

A bead of mass m is constrained to move along a (massless) wire which is
bent into a vertical loop of radius R. The circular wire rotates about its ver-
tical axis with a constant angular velocity Ω, and gravity acts on the bead.
Find an equation of motion for θ(t), where θ is the angle measured from the
bottom of the wire loop to the loop centre.

The kinetic energy is a combination of the motion of the bead along the
wire in θ, and the bead spinning around with the wire. A bit of geometry
shows that at an angle θ, the bead is at a radius of R sin θ, so its spin velocity
is RΩ sin θ. Hence the kinetic energy is
 !2 
1 dθ
T = mR2  + Ω2 sin2 θ
2 dt

166
The potential energy due to gravity acting on the bead is:

V = mgR(1 − cos θ)

We take the derivative of the Lagrangian w.r.t. θ:

∂L ∂L dθ
= mR2 Ω2 sin θcosθ − mgR sin θ = mR2
∂θ ∂ θ̇ dt
So via the Euler-Lagrange equation we have:

d2 θ
mR2 Ω2 sin θ cos θ − mgR sin θ = mR2
dt2
This directly gives the equation of motion of the bead:

d2 θ g
 
2
= sin θ Ω cos θ −
dt2 R
While this is complicated to solve fully, it can readily be seen that there
are stationary points with d2 θ/dt2 = 0 at θ = 0 and θ = π. For Ω2 > g/R
there are two additional stationary points satisfying cos θ = g/RΩ2 . At these
points, the force from gravity pulling the bead down the wire is equal to the
centrifugal force pushing it up the wire. One can find whether these are
points of stable or unstable equilibrium by considering small displacements
in θ around the points.

17.4.4 Block on a wedge

Another type of dynamical problem that is best solved by Lagrangian meth-
ods is one in which there is a superposition of two difference types of motion.
Here we look at a combination of a wedge of mass M sliding along a smooth
horizontal surface, and a small block of mass m sliding down the smooth
inclined upper surface of the wedge at an angle of α to the horizontal.

We must first define our set of Euler-Lagrange coordinates {qi }. Let us

define the horizontal direction as x, and the direction down the slope as ℓ.
Note that these are not orthogonal, but that’s fine; they are the two variables
of interest in the system, so it’s easiest to choose these as our “coordinates”.
Being able to directly use the quantities of interest as generalised coordi-
nates, without having to resolve into orthonormal components, is another

167
advantage of Lagrangian dynamics.

The kinetic energy comes from the wedge moving in x, plus the block moving
in ℓ. We must remember that ⃗ℓ is defined in the opposite sense to ⃗x, since
the block and wedge move in opposite directions.
!2  2
1 dx 1 d⃗ℓ d⃗x 
T = M + m − +
2 dt 2 dt dt

Note that although we are mathematically free to choose {qi } however we

want, when we write down the physical energies, we must still account for
the fact that they are not orthonormal. Geometrically, it is simple to show
that ℓ̂ · x̂ = − cos α. Hence we can expand out the second term to get:
!2  !2 ! !
1 dx 1 dℓ dℓ dx 
T = (M + m) + m + 2 cos α
2 dt 2 dt dt dt

The potential energy comes from m’s height on the wedge ℓ sin α, with a pos-
itive change in ℓ resulting in a lower potential (since it is moving downwards
along the wedge):
V = −mgℓ sin α
So now we have our Lagrangian. Applying Euler-Lagrange to the x coordi-
nate, we see that there is no explicit x dependence in L, so:
d ∂L d2 x d2 ℓ
= (M + m) 2 + m cos α 2 = 0
dt ∂ ẋ dt dt
∂L
Hence ∂ ẋ
is conserved. If we integrate w.r.t. t we get:

∂L dx dℓ
= (M + m) + m cos α = constant
∂ ẋ dt dt
This you might recognise as conservation of momentum in the x direction.
One can thus think of ∂∂L q˙i
as a generalised momentum for a given coordinate
qi , which is conserved if L has no explicit dependence on qi .

We now apply Euler-Lagrange to the ℓ coordinate:

d2 ℓ d2 x
m + m cos α = mg sin α
dt2 dt2
168
The two Euler-Lagrange equations represent coupled equations in x and ℓ. If
we multiply the second equation by cos α, we can subtract the two resulting
equations to remove the ℓ dependence, which (after some algebra) gives:
d2 x mg sin α cos α
2
=−
dt M + m(1 − cos2 α)
Similarly, one can multiply the second equation by (M + m)/(m cos α) and
subtract, which gives:
d2 ℓ (M + m)g sin α
2
=
dt M + m(1 − cos2 α)
The RHS of each equation is a constant. Hence these can be independently
solved to get the motion of the wedge dx/dt and the block dℓ/dt. With La-
grangian dynamics, this complicated problem turns out to be quite straight-
forward to solve.

As a check we look at some specific cases. For α = 90◦ , the wedge is vertical,
and we have d2 ℓ/dt2 = g and d2 x/dt2 = 0 which is just free fall, as expected.
If we take M → ∞, we have d2 ℓ/dt2 → g sin α and d2 x/dt2 → 0, which is
the well known result for a block sliding down a stationary slope. In the
limit M → 0 we have d2 x/dt2 = cos α d2 ℓ/dt2 which states that as the block
moves down, it pushes the wedge with an acceleration exactly matching its
downwards motion.

17.4.5 Double pendulum

For another example of the Lagrangian approach, we revisit the double pen-
dulum, this time more generally. A bob of mass m2 is suspended by a length
L2 from a mass m1 which is itself suspended by a length L1 from a fixed
point. For the generalised “coordinates” of the system, we choose θ1 and θ2 ,
since these are the quantities for which we are interested in obtaining the
equations of motion.

The potential energy of this system corresponds to the height above the
equilibrium position for each bob:

V = (m1 + m2 )gL1 (1 − cos θ1 ) + m2 gL2 (1 − cos θ2 )

169
With a small angle approximation cos θ = 1 − θ2 /2 this becomes:
1 1
V = (m1 + m2 )gL1 θ12 + m2 gL2 θ22
2 2
Meanwhile, the kinetic energy of the system corresponds to the velocities of
the bobs, where m2 ’s velocity is a combination of its own plus m1 ’s:
!2 " ! !#2
1 dθ1 1 dθ1 dθ2
T = m1 L21 + m2 L1 + L2
2 dt 2 dt dt
!2 ! ! !2
1 dθ1 dθ1 dθ2 1 dθ2
T = (m1 + m2 )L21 + m2 L1 L2 + m2 L22
2 dt dt dt 2 dt
Taking derivatives w.r.t. θ1 , and using Euler-Lagrange, we get:

d2 θ1 d2 θ2
(m1 + m2 )L21 + m2 L 1 L 2 = −(m1 + m2 )gL1 θ1
dt2 dt2
Now we do the same for θ2 :

d2 θ1 2
2 d θ2
m2 L1 L2 + m L
2 2 = −m2 gL2 θ2
dt2 dt2
These are our coupled ODEs for this system. We can write this in matrix
form as:
d2 θ
M 2 = −Kθ
dt
where the matrices are:
! ! !
(m1 + m2 )L21 m2 L1 L2 (m1 + m2 )L1 0 θ1
M= K=g θ=
m2 L1 L2 m2 L22 0 m2 L2 θ2

This exactly reduces to what we got earlier, for L1 = L2 and m1 = m2 . Using

a Lagrangian, even the more general case was considerably easier to derive.
From here, recall that one assumes a solution of the form eλt , which reduces
to solving the determinant of the matrix to obtain the eigenvalues λ. This
part proceeds exactly as we did before.

170
17.4.6 Extra: The top
A top spinning about its axis while precessing is one of the most complex
systems in classical dynamics. The rotation axis is â with angular velocity
ω, and it precesses about the vertical axis k̂ with angular velocity Ω.

ω
⃗ = ωâ + Ωk̂

The angle θ between â and k̂ is constant for steady precession, and the
rotation axes are related by

k̂ = cos θâ − sin θâ⊥

The angular momentum vector is given by:

⃗ = −I⊥ Ω sin θâ⊥ + I∥ (ω + Ω cos θ)â
L

where the I are the principal moments of inertia perpendicular and parallel
to â. The total angular momentum about the axis â is conserved because
there is no torque about this axis:

La = I∥ (ω + Ω cos θ) = I∥ Ωa = constant

The kinetic energy is given by:

1 1
T = I⊥ Ω2 sin2 θ + I∥ (ω 2 + 2ωΩ cos θ + Ω2 cos2 θ)
2 2
The potential energy is just:

V = M gh cos θ

where h is the distance along the axis â between the CoM of the top and the
point of contact with the table.
Taking derivatives of the Lagrangian w.r.t. θ:

I⊥ Ω2 sin θ cos θ − I∥ (ωΩ sin θ + Ω2 cos θ sin θ) + M gh sin θ = 0

Removing a common factor of sin θ and simplifying the second term using
La we are left with a quadratic in the frequency Ω:

I⊥ Ω2 cos θ − La Ω + M gh = 0

171
 q
La ± L2a − 4M ghI⊥ cos θ
Ω=
2I⊥ cos θ
There is a small solution with the −ve sign which can be approximated for
a high spinning rate of the top Ωa ≫ Ω:

La − La + 2M ghI⊥ cos θ/La M gh

Ω= =
2I⊥ cos θ La

This is the steady precession frequency of the axis â about the k̂ direction.
It can be interpreted as being due to the torque of the force M g acting round
the point of contact with the table.

The large solution with the +ve sign and the same high spin limit is:
La
Ω=
I⊥ cos θ
This is the projection of the La component of the angular momentum onto
the k̂ direction.

To get the more general solution with non-steady precession you need to
add another term to the kinetic energy:
!2
1 dθ
T = Tsteady + I⊥
2 dt

leading to the equation of motion in θ:

d2 θ
I⊥ = (I⊥ Ω2 cos θ − La Ω + M gh) sin θ
dt2
The variation in θ as the top spins is known as nutation. For small θ the
motion is still stable, but at larger θ and lower ω it becomes unstable.

17.5 Extra: Hamiltonian Dynamics

Another oft-used formulation of dynamics is known as Hamiltonian dynam-
ics, which is an extension of Lagrangian dynamics. In Hamiltonian dynamics,

172
we introduce the Hamiltonian H = T + V (in mechanics, H is just the total
energy), along with a generalized momentum

dL
pi ≡
dq˙i
The time evolution of the system is then given by two Hamilton’s equations:
dpi ∂H dqi ∂H
=− , =
dt ∂qi dt ∂pi
These equations can be written for each independent coordinate qi , with its
associated generalized momentum pi .
A simple application to dynamics illustrates the physical interpretation of
Hamiltonian mechanics. For T = p2 /2m and V = V (q), Hamilton’s equations
are
dpi dV dqi p
=− = Fi , = vi =
dt dqi dt m
The first equation is N2, while the second shows that in dynamics the gen-
eralized momentum is actually just the momentum.

173