Lectures On Differential Equations: Department of Mathematics University of California Davis, CA 95616

Lectures on Differential Equations1
Craig A. Tracy2
Department of Mathematics
University of California
Davis, CA 95616
January 2015
1
c Craig A. Tracy, 2000, 2014 Davis, CA 95616
2
email: tracy@math.ucdavis.edu
2
Contents
1 Introduction 1
1.1 What is a differential equation? . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Differential equation for the pendulum . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Introduction to computer software . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2 First Order Equations & Conservative Systems 15

2.1 Linear first order equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2 Conservative systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3 Level curves of the energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3 Second Order Linear Equations 43

3.1 Theory of second order equations . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.2 Reduction of order . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.3 Constant coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.4 Forced oscillations of the mass-spring system . . . . . . . . . . . . . . . . . . 54
3.5 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4 Difference Equations 61
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.2 Constant coefficient difference equations . . . . . . . . . . . . . . . . . . . . . 62
4.3 Inhomogeneous difference equations . . . . . . . . . . . . . . . . . . . . . . . 64
4.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
i
ii CONTENTS
5 Matrix Differential Equations 67

5.1 The matrix exponential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.2 Application of matrix exponential to DEs . . . . . . . . . . . . . . . . . . . . 70
5.3 Relation to earlier methods of solving constant coefficient DEs . . . . . . . . 74
5.4 Problem from Markov processes . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.5 Application of matrix DE to radioactive decays . . . . . . . . . . . . . . . . . 77
5.6 Inhomogenous matrix equations . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.7 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
6 Weighted String 89
6.1 Derivation of differential equations . . . . . . . . . . . . . . . . . . . . . . . . 90
6.2 Reduction to an eigenvalue problem . . . . . . . . . . . . . . . . . . . . . . . 92
6.3 Computation of the eigenvalues . . . . . . . . . . . . . . . . . . . . . . . . . . 93
6.4 The eigenvectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
6.5 Determination of constants . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
6.6 Continuum limit: The wave equation . . . . . . . . . . . . . . . . . . . . . . . 101
6.7 Inhomogeneous problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
6.8 Vibrating membrane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
6.9 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
7 Quantum Harmonic Oscillator 121

7.1 Schrödinger equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
7.2 Harmonic oscillator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
7.3 Some properties of the harmonic oscillator . . . . . . . . . . . . . . . . . . . . 132
7.4 The Heisenberg Uncertainty Principle . . . . . . . . . . . . . . . . . . . . . . 135
7.5 Comparison of three problems . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
7.6 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
8 Heat Equation 141

8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
8.2 Fourier transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
8.3 Solving the heat equation by use of the Fourier transform . . . . . . . . . . . 143
8.4 Heat equation on the semi-infinite rod . . . . . . . . . . . . . . . . . . . . . . 146
8.5 Heat equation on the circle . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
8.6 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
CONTENTS iii
9 Laplace Transform 151

9.1 Matrix version . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
9.2 Structure of (sIn − A)−1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
9.3 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
iv CONTENTS
Preface
Figure 1: Sir Isaac Newton, December 25, 1642–March 20, 1727 (Julian Calendar).
These notes are for a one-quarter course in differential equations. The approach is to tie
the study of differential equations to specific applications in physics with an emphasis on
oscillatory systems.
Mathematics is a part of physics. Physics is an experimental science, a part of

natural science. Mathematics is the part of physics where experiments are cheap.
In the middle of the twentieth century it was attempted to divide physics and
mathematics. The consequences turned out to be catastrophic. Whole gener-
ations of mathematicians grew up without knowing half of their science and,
of course, in total ignorance of any other sciences. They first began teaching
their ugly scholastic pseudo-mathematics to their students, then to schoolchil-
dren (forgetting Hardy’s warning that ugly mathematics has no permanent place
under the Sun).
Since scholastic mathematics that is cut off from physics is fit neither for teaching
nor for application in any other science, the result was the universal hate towards
mathematicians—both on the part of the poor schoolchildren (some of whom in
the meantime became ministers) and of the users.
V. I. Arnold, On Teaching Mathematics
CONTENTS v
Newton’s fundamental discovery, the one which he considered necessary to keep

secret and published only in the form of an anagram, consists of the following:
Data aequatione quotcunque fluentes quantitae involvente fluxions invenire et
vice versa. In contemporary mathematical language, this means: “It is useful to
solve differential equations”.
V. I. Arnold, Geometrical Methods in the Theory of Ordinary Differential
Equations.
I thank Eunghyun (Hyun) Lee for his help with these notes during the 2008–09 academic
year. Also thanks to Andrew Waldron for his comments on the notes.
Craig Tracy, Sonoma, California
vi CONTENTS
Notation
Symbol Definition of Symbol
R field of real numbers
Rn the n-dimensional vector space with each component a real number
C field of complex numbers
ẋ the derivative dx/dt, t is interpreted as time
ẍ the second derivative d2 x/dt2 , t is interpreted as time
:= equals by definition
Ψ = Ψ(x, t) wave function in quantum mechanics
ODE ordinary differential equation
PDE partial differential equation
KE kinetic energy
PE potential energy
det determinant
δij the Kronecker delta, equal to 1 if i = j and 0 otherwise
L
the Laplace transform operator
n
k The binomial coefficient n choose k.
Maple is a registered trademark of Maplesoft.
Mathematica is a registered trademark of Wolfram Research.
MatLab is a registered trademark of the MathWorks, Inc.
Chapter 1
Introduction
Figure 1.1: Galileo Galilei, 1564–1642. From The Galileo Project : “Galileo’s discovery was
that the period of swing of a pendulum is independent of its amplitude–the arc of the swing–
the isochronism of the pendulum. Now this discovery had important implications for the
measurement of time intervals. In 1602 he explained the isochronism of long pendulums in
a letter to a friend, and a year later another friend, Santorio Santorio, a physician in Venice,
began using a short pendulum, which he called “pulsilogium,” to measure the pulse of his
patients. The study of the pendulum, the first harmonic oscillator, date from this period.”
See the You Tube video http://youtu.be/MpzaCCbX-z4.
1
2 CHAPTER 1. INTRODUCTION
1.1 What is a differential equation?

From Birkhoff and Rota [3]
A differential equation is an equation between specified derivative on an unknown

function, its values, and known quantities and functions. Many physical laws are
most simply and naturally formulated as differential equations (or DEs, as we
will write for short). For this reason, DEs have been studied by the greatest
mathematicians and mathematical physicists since the time of Newton.
Ordinary differential equations are DEs whose unknowns are functions of a single
variable; they arise most commonly in the study of dynamical systems and elec-
trical networks. They are much easier to treat that partial differential equations,
whose unknown functions depend on two or more independent variables.
Ordinary DEs are classified according to their order. The order of a DE is
defined as the largest positive integer, n, for which an nth derivative occurs in
the equation. Thus, an equation of the form
φ(x, y, y ) = 0
is said to be of the first order.
From Wikipedia
A differential equation is a mathematical equation that relates some function of

one or more variables with its derivatives. Differential equations arise whenever
a deterministic relation involving some continuously varying quantities (mod-
eled by functions) and their rates of change in space and/or time (expressed as
derivatives) is known or postulated. Because such relations are extremely com-
mon, differential equations play a prominent role in many disciplines including
engineering, physics, economics, and biology.
Differential equations are mathematically studied from several different perspec-
tives, mostly concerned with their solutions the set of functions that satisfy
the equation. Only the simplest differential equations admit solutions given
by explicit formulas; however, some properties of solutions of a given differential
equation may be determined without finding their exact form. If a self-contained
formula for the solution is not available, the solution may be numerically approx-
imated using computers. The theory of dynamical systems puts emphasis on
qualitative analysis of systems described by differential equations, while many
numerical methods have been developed to determine solutions with a given
degree of accuracy.
Many fundamental laws of physics and chemistry can be formulated as differential
equations. In biology and economics, differential equations are used to model the
behavior of complex systems. The mathematical theory of differential equations
first developed together with the sciences where the equations had originated
and where the results found application. However, diverse problems, sometimes
originating in quite distinct scientific fields, may give rise to identical differential
equations. Whenever this happens, mathematical theory behind the equations
1.1. WHAT IS A DIFFERENTIAL EQUATION? 3
can be viewed as a unifying principle behind diverse phenomena. As an example,

consider propagation of light and sound in the atmosphere, and of waves on the
surface of a pond. All of them may be described by the same second-order partial
differential equation, the wave equation, which allows us to think of light and
sound as forms of waves, much like familiar waves in the water. Conduction
of heat, the theory of which was developed by Joseph Fourier, is governed by
another second-order partial differential equation, the heat equation. It turns
out that many diffusion processes, while seemingly different, are described by
the same equation; the Black–Scholes equation in finance is, for instance, related
to the heat equation.
1.1.1 Examples
1. A simple example of a differential equation (DE) is
dy
= λy
dx
where λ is a constant. The unknown is y and the independent variable is x. The
equation involves both the unknown y as well as the unknown dy/dx; and for this
reason is called a differential equation. We know from calculus that
y(x) = c eλx , c = constant,
satisfies this equation since

dy d λx
= c e = c λeλx = λy(x).
dx dx
The constant c is uniquely specified once we give the initial condition which in this
case would be to give the value of y(x) at a particular point x0 . For example, if we
impose the initial condition y(0) = 3, then the constant c is now determined, i.e. c = 3.
2. Consider the DE
dy
= y2
dx
subject to the initial condition y(0) = 1. This DE was solved in your calculus courses
using the method of separation of variables:
• First rewrite DE in differential form:
dy = y 2 dx
• Now separate variables (all x’s on one side and all y’s on the other side):
dy
= dx
y2
• Now integrate both sides

1
− =x+c
y
where c is a constant to be determined.
• Solve for y = y(x)

1
y(x) = −
x+c
• Now require y(0) = −1/c to equal the given initial condition:
1
− =1
c
Solving this gives c = −1 and hence the solution we want is
1
y(x) = −
x−1
3. An example of a second order ODE is
d2 x
F (x) = m (1.1)
dt2
where F = F (x) is a given function of x and m is a positive constant. Now the
unknown is x and the independent variable is t. The problem is to find functions
x = x(t) such that when substituted into the above equation it becomes an identity.
Here is an example; choose F (x) = −kx where k > 0 is a positive number. Then (1.1)
reads
d2 x
−kx = m 2
dt
We rewrite this ODE as
d2 x k
+ x = 0. (1.2)
dt2 m
You can check that
k
x(t) = sin t
m
satisfies (1.2). Can you find other functions that satisfy this same equation? One of the
problems in differential equations is to find all solutions x(t) to the given differential
equation. We shall soon prove that all solutions to (1.2) are of the form

k k
x(t) = c1 sin t + c2 cos t (1.3)
m m
where c1 and c2 are arbitrary constants. Using differential calculus1 one can verify
that (1.3) when substituted into (1.2) satisfies the differential equation (show this!).
It is another matter to show that all solutions to (1.2) are of the form (1.3). This is
a problem we will solve in this class.
1 Recall the differentiation formulas
d d
cos(ωt) = −ω sin(ωt)
sin(ωt) = ω cos(ωt),
dt dt

where ω is a constant. In the above the constant ω = k/m.
1.2. DIFFERENTIAL EQUATION FOR THE PENDULUM 5
1.2 Differential equation for the pendulum

Newton’s principle of determinacy
The initial state of a mechanical system (the totality of positions and velocities
of its points at some moment of time) uniquely determines all of its motion.
It is hard to doubt this fact, since we learn it very early. One can imagine a world
in which to determine the future of a system one must also know the acceleration
at the initial moment, but experience shows us that our world is not like this.
V. I. Arnold, Mathematical Methods of Classical Mechanics[1]
Many interesting ordinary differential equations (ODEs) arise from applications. One
reason for understanding these applications in a mathematics class is that you can combine
your physical intuition with your mathematical intuition in the same problem. Usually the
result is an improvement of both. One such application is the motion of pendulum, i.e. a
ball of mass m suspended from an ideal rigid rod that is fixed at one end. The problem
is to describe the motion of the mass point in a constant gravitational field. Since this is
a mathematics class we will not normally be interested in deriving the ODE from physical
principles; rather, we will simply write down various differential equations and claim that
they are “interesting.” However, to give you the flavor of such derivations (which you will see
repeatedly in your science and engineering courses), we will derive from Newton’s equations
the differential equation that describes the time evolution of the angle of deflection of the
pendulum.
Let
= length of the rod measured, say, in meters,

m = mass of the ball measured, say, in kilograms,
g = acceleration due to gravity = 9.8070 m/s2 .
The motion of the pendulum is confined to a plane (this is an assumption on how the rod
is attached to the pivot point), which we take to be the xy-plane (see Figure 1.2). We treat
the ball as a “mass point” and observe there are two forces acting on this ball: the force
due to gravity, mg, which acts vertically downward and the tension T in the rod (acting in
the direction indicated in figure). Newton’s equations for the motion of a point x in a plane
are vector equations2
F = ma
where F is the sum of the forces acting on the the point and a is the acceleration of the
point, i.e.
d2 x
a = 2 .
dt
Since acceleration is a second derivative with respect to time t of the position vector, x,
Newton’s equation is a second-order ODE for the position x. In x and y coordinates Newton’s
2 In your applied courses vectors are usually denoted with arrows above them. We adopt this notation
when discussing certain applications; but in later chapters we will drop the arrows and state where the
quantity lives, e.g. x ∈ R2 .
equations become two equations
d2 x d2 y
Fx = m , Fy = m ,
dt2 dt2
where Fx and Fy are the x and y components, respectively, of the force F . From the figure
(note definition of the angle θ) we see, upon resolving T into its x and y components, that
Fx = −T sin θ, Fy = T cos θ − mg.
(T is the magnitude of the vector T .)
θ
T
mass m
mg
Figure 1.2: Simple pendulum
Substituting these expressions for the forces into Newton’s equations, we obtain the
differential equations
d2 x
−T sin θ = m , (1.4)
dt2
d2 y
T cos θ − mg = m 2. (1.5)
dt
From the figure we see that
x = sin θ, y = − cos θ. (1.6)
(The origin of the xy-plane is chosen so that at x = y = 0, the pendulum is at the bottom.)
Differentiating3 (1.6) with respect to t, and then again, gives
ẋ = cos θ θ̇,
ẍ = cos θ θ̈ − sin θ (θ̇)2 , (1.7)
ẏ = sin θ θ̇,
ÿ = sin θ θ̈ + cos θ (θ̇)2 . (1.8)
Substitute (1.7) in (1.4) and (1.8) in (1.5) to obtain
−T sin θ = m cos θ θ̈ − m sin θ (θ̇)2 , (1.9)

T cos θ − mg = m sin θ θ̈ + m cos θ (θ̇)2 . (1.10)
3 We use the dot notation for time derivatives, e.g. ẋ = dx/dt, ẍ = d2 x/dt2 .
1.2. DIFFERENTIAL EQUATION FOR THE PENDULUM 7
Now multiply (1.9) by cos θ, (1.10) by sin θ, and add the two resulting equations to obtain
−mg sin θ = mθ̈,
or
g
θ̈ + sin θ = 0. (1.11)

Remarks
• The ODE (1.11) is called a second-order equation because the highest derivative ap-
pearing in the equation is a second derivative.
• The ODE is nonlinear because of the term sin θ (this is not a linear function of the
unknown quantity θ).
• A solution to this ODE is a function θ = θ(t) such that when it is substituted into the
ODE, the ODE is satisfied for all t.
• Observe that the mass m dropped out of the final equation. This says the motion will
be independent of the mass of the ball. If an experiment is performed, will we observe
this to be the case; namely, the motion is independent of the mass m? If not, perhaps
in our model we have left out some forces acting in the real world experiment. Can
you think of any?
• The derivation was constructed so that the tension, T , was eliminated from the equa-
tions. We could do this because we started with two unknowns, T and θ, and two
equations. We manipulated the equations so that in the end we had one equation for
the unknown θ = θ(t).
• We have not discussed how the pendulum is initially started. This is very important
and such conditions are called the initial conditions.
We will return to this ODE later in the course. At this point we note that if we were
interested in only small deflections from the origin (this means we would have to start out
near the origin), there is an obvious approximation to make. Recall from calculus the Taylor
expansion of sin θ
θ3 θ5
sin θ = θ − + + ··· .
3! 5!
For small θ this leads to the approximation sin θ ≈ θ . Using this small deflection approxi-
mation in (1.11) leads to the ODE
g
θ̈ + θ = 0. (1.12)

We will see that (1.12) is mathematically simpler than (1.11). The reason for this is that
(1.12) is a linear ODE. It is linear because the unknown quantity, θ, and its derivatives
appear only to the first or zeroth power. Compare (1.12) with (1.2).
1.3 Introduction to MatLab, Mathematica and Maple
In this class we may use the computer software packages MatLab, Mathematica or Maple
to do routine calculations. It is strongly recommended that you learn to use at least one of
these software packages. These software packages take the drudgery out of routine calcula-
tions in calculus and linear algebra. Engineers will find that MatLab is used extenstively
in their upper division classes. Both MatLab and Maple are superior for symbolic compu-
tations (though MatLab can call Maple from the MatLab interface).
1.3.1 MatLab
What is MatLab ? “MatLab is a powerful computing system for handling the calculations
involved in scientific and engineering problems.”4 MatLab can be used either interactively
or as a programming language. For most applications in Math 22B it suffices to use MatLab
interactively. Typing matlab at the command level is the command for most systems to
start MatLab . Once it loads you are presented with a prompt sign >>. For example if I
enter
>> 2+22
and then press the enter key it responds with
ans=24
Multiplication is denoted by * and division by / . Thus, for example, to compute
(139.8)(123.5 − 44.5)
125
we enter
>> 139.8*(123.5-44.5)/125
gives
ans=88.3536
MatLab also has a Symbolic Math Toolbox which is quite useful for routine calculus
computations. For example, suppose you forgot the Taylor expansion of sin x that was used
in the notes just before (1.12). To use the Symbolic Math Toolbox you have to tell MatLab
that x is a symbol (and not assigned a numerical value). Thus in MatLab
>> syms x
>> taylor(sin(x))
4 Brian D. Hahn, Essential MatLab for Scientists and Engineers.
1.3. INTRODUCTION TO COMPUTER SOFTWARE 9
gives
ans = x -1/6*x^3+1/120*x^5
Now why did taylor expand about the point x = 0 and keep only through x5 ? By
default the Taylor series about 0 up to terms of order 5 is produced. To learn more about
taylor enter
>> help taylor
from which we learn if we had wanted terms up to order 10 we would have entered
>> taylor(sin(x),10)
If we want the Taylor expansion of sin x about the point x = π up to order 8 we enter
>> taylor(sin(x),8,pi)
A good reference for MatLab is MatLab Guide by Desmond Higham and Nicholas Higham.
1.3.2 Mathematica
There are alternatives to the software package MatLab. Two widely used packages are
Mathematica and Maple. Here we restrict the discussion to Mathematica . Here are
some typical commands in Mathematica .
1. To define, say, the function f (x) = x2 e−2x one writes in Mathematica
f[x_]:=x^2*Exp[-2*x]
∞2. One can now use f in other Mathematica

2 −2x
commands. For example, suppose we want
0
f (x) dx where as above f (x) = x e . The Mathematica command is
Integrate[f[x],{x,0,Infinity}]
Mathematica returns the answer 1/4.

3. In Mathematica to find the Taylor series of sin x about the point x = 0 to fifth order
you would type
Series[Sin[x],{x,0,5}]
4. Suppose we want to create the 10 × 10 matrix

1
M= .
i + j + 1 1≤i,j≤10
In Mathematica the command is

M=Table[1/(i+j+1),{i,1,10},{j,1,10}];
(The semicolon tells Mathematica not to write out the result.) Suppose we then want the
determinant of M . The command is
Det[M]
Mathematica returns the answer
1/273739709893086064093902013446617579389091964235284480000000000
If we want this number in scientific notation, we would use the command N[· ] (where the
number would be put in place of ·). The answer Mathematica returns is 3.65311 × 10−63 .
The (numerical) eigenvalues of M are obtained by the command
N[Eigenvalues[M]]
Mathematica returns the list of 10 distinct eigenvalues. (Which we won’t reproduce here.)
The reason for the N[·] is that Mathematica cannot find an exact form for the eigenvalues,
so we simply ask for it to find approximate numerical values. To find the (numerical)
eigenvectors of M , the command is
N[Eigenvectors[M]]
5. Mathematica has nice graphics capabilities. Suppose we wish to graph the function
f (x) = 3e−x/10 sin(x) in the interval 0 ≤ x ≤ 50. The command is
Plot[3*Exp[-x/10]*Sin[x],{x,0,50},PlotRange->All,
AxesLabel->{x},PlotLabel->3*Exp[-x/10]*Sin[x]]
The result is the graph shown in Figure 1.3.

1.4. EXERCISES 11
3 x10 sinx
x
10 20 30 40 50
1
2
Figure 1.3:
1.4 Exercises
#1. MatLab and/or Mathematica Exercises
1. Use MatLab or Mathematica to get an estimate (in scientific notation) of 9999 . Now
use
>> help format
to learn how to get more decimal places. (All MatLab computations are done to a
relative precision of about 16 decimal places. MatLab defaults to printing out the
first 5 digits.) Thus entering
>> format long e
on a command line and then re-entering the above computation will give the 16 digit
answer.
In Mathematica to get 16 digits accuracy the command is
N[99^(99),16]
Ans.: 3.697296376497268 × 10197 .


2. Use MatLab to compute sin(π/7). (Note that MatLab has the special symbol pi;
that is pi ≈ π = 3.14159 . . . to 16 digits accuracy.)
In Mathematica the command is
N[Sqrt[Sin[Pi/7]],16]
3. Use MatLab or Mathematica to find the determinant, eigenvalues and eigenvectors

of the 4 × 4 matrix  
√1 −1 2 0
 2 1 
A= √0 −2 
 0 1 2 −1 
1 2 2 0
Hint: In MatLab you enter the matrix A by
>> A=[1 -1 2 0; sqrt(2) 1 0 -2;0 1 sqrt(2) -1; 1 2 2 0]
To find the determinant
>> det(A)
and to find the eigenvalues
>> eig(A)
If you also want the eigenvectors you enter
>> [V,D]=eig(A)
In this case the columns of V are the eigenvectors of A and the diagonal elements of
D are the corresponding eigenvalues. Try this now to find the eigenvectors. For the
determinant you should get the result 16.9706. One may also calculate the determi-
nant symbolically. First we tell MatLab that A is to be treated as a symbol (we are
assuming you have already entered A as above):
>> A=sym(A)
and then re-enter the command for the determinant
det(A)
and this time MatLab returns
ans =
12*2^(1/2)
√
that is, 12 2 which is approximately equal to 16.9706.
4. Use MatLab or Mathematica to plot sin θ and compare this with the approximation
sin θ ≈ θ. For 0 ≤ θ ≤ π/2, plot both on the same graph.
1.4. EXERCISES 13
#2. Inverted pendulum
This exercise derives the small angle approximation to (1.11) when the pendulum is nearly
inverted, i.e. θ ≈ π. Introduce
φ=θ−π
and derive a small φ-angle approximation to (1.11). How does the result differ from (1.12)?
Chapter 2
First Order Equations &

Conservative Systems
2.1 Linear first order equations
2.1.1 Introduction
The simplest differential equation is one you already know from calculus; namely,
dy
= f (x). (2.1)
dx
To find a solution to this equation means one finds a function y = y(x) such that its
derivative, dy/dx, is equal to f (x). The fundamental theorem of calculus tells us that all
solutions to this equation are of the form
x
y(x) = y0 + f (s) ds. (2.2)
x0
Remarks:
• y(x0 ) = y0 and y0 is arbitrary. That is, there is a one-parameter family of solutions;

y = y(x; y0 ) to (2.1). The solution is unique once we specify the initial condition
y(x0 ) = y0 . This is the solution to the initial value problem. That is, we have found
a function that satisfies both the ODE and the initial value condition.
• Every calculus student knows that differentiation is easier than integration. Observe
that solving a differential equation is like integration—you must find a function such
that when it and its derivatives are substituted into the equation the equation is
identically satisfied. Thus we sometimes say we “integrate” a differential equation. In
the above case it is exactly integration as you understand it from calculus. This also
suggests that solving differential equations can be expected to be difficult.
• For the integral to exist in (2.2) we must place some restrictions on the function f
appearing in (2.1); here it is enough to assume f is continuous on the interval [a, b].
It was implicitly assumed in (2.1) that x was given on some interval—say [a, b].
15
16 CHAPTER 2. FIRST ORDER EQUATIONS & CONSERVATIVE SYSTEMS
A simple generalization of (2.1) is to replace the right-hand side by a function that

depends upon both x and y
dy
= f (x, y).
dx
Some examples are f (x, y) = xy 2 , f (x, y) = y, and the case (2.1). The simplest choice in
terms of the y dependence is for f (x, y) to depend linearly on y. Thus we are led to study
dy
= g(x) − p(x)y,
dx
where g(x) and p(x) are functions of x. We leave them unspecified. (We have put the minus
sign into our equation to conform with the standard notation.) The conventional way to
write this equation is
dy
+ p(x)y = g(x). (2.3)
dx
It’s possible to give an algorithm to solve this ODE for more or less general choices of p(x)
and g(x). We say more or less since one has to put some restrictions on p and g—that they
are continuous will suffice. It should be stressed at the outset that this ability to find an
explicit algorithm to solve an ODE is the exception—most ODEs encountered will not be
so easily solved.
But before we give the general solution to (2.3), let’s examine the special case p(x) = −1
and g(x) = 0 with initial condition y(0) = 1. In this case the ODE becomes
dy
=y (2.4)
dx
and the solution we know from calculus
y(x) = ex .
In calculus one typically defines ex as the limit

x n
ex := lim 1 +
n→∞ n
or less frequently as the solution y = y(x) to the equation
y
dt
x= .
1 t
In calculus courses one then proves from either of these starting points that the derivative
of ex equals itself. One could also take the point of view that y(x) = ex is defined to be the
(unique) solution to (2.4) satisfying the initial condition y(0) = 1. Taking this last point of
view, can you explain why the Taylor expansion of ex ,
∞
x xn
e = ,
n=0
n!
follows almost immediately?

2.1. LINEAR FIRST ORDER EQUATIONS 17
2.1.2 Method of integrating factors

If (2.3) were of the form (2.1), then we could immediately write down a solution in terms
of integrals. For (2.3) to be of the form (2.1) means the left-hand side is expressed as the
derivative of our unknown quantity. We have some freedom in making this happen—for
instance, we can multiply (2.3) by a function, call it µ(x), and ask whether the resulting
equation can be put in form (2.1). Namely, is
dy d
+ µ(x)p(x)y =
µ(x) (µ(x)y) ? (2.5)
dx dx
Taking derivatives we ask can µ be chosen so that
dy dy dµ
µ(x) + µ(x)p(x)y = µ(x) + y
dx dx dx
holds? This immediately simplifies to1
dµ
µ(x)p(x) = ,
dx
or
d
log µ(x) = p(x).
dx
Integrating this last equation gives

log µ(x) = p(s) ds + c.
Taking the exponential of both sides (one can check later that there is no loss in generality
if we set c = 0) gives2 x
µ(x) = exp p(s) ds . (2.6)
Defining µ(x) by (2.6), the differential equation (2.5) is transformed to

d
(µ(x)y) = µ(x)g(x).
dx
This last equation is precisely of the form (2.1), so we can immediately conclude
x
µ(x)y(x) = µ(s)g(s) ds + c,
and solving this for y gives our final formula
x
1 c
y(x) = µ(s)g(s) ds + , (2.7)
µ(x) µ(x)
where µ(x), called the integrating factor, is defined by (2.6). The constant c will be deter-
mined from the initial condition y(x0 ) = y0 .
1 Notice y and its first derivative drop out. This is a good thing since we wouldn’t want to express µ in
terms of the unknown quantity y.

2 By the symbol x f (s) ds we mean the indefinite integral of f in the variable x.
An example
Suppose we are given the DE

dy 1
+ y = x2 , x > 0
dx x
with initial condition
y(1) = 2.
This is of form (2.3) with p(x) = 1/x and g(x) = x2 . We apply formula (2.7):
• First calculate the integrating factor µ(x):

1
µ(x) = exp p(x) dx = exp dx = exp(log x) = x.
x
• Now substitute into (2.3)

1 c 1 x4 c x3 c
y(x) = x · x2 dx + = · + = + .
x x x 4 x 4 x
• Impose the initial condition y(1) = 2:

1 7
+ c = 2, solve for c, c = .
4 4
• Solution to DE is
x3 7
y(x) = + .
4 4x
2.1.3 Application to mortgage payments
Suppose an amount P , called the principal, is borrowed at an interest I (100I%) for a period
of N years. One is to make monthly payments in the amount D/12 (D equals the amount
paid in one year). The problem is to find D in terms of P , I and N . Let
y(t) = amount owed at time t (measured in years).
We have the initial condition
y(0) = P (at time 0 the amount owed is P ).
We are given the additional information that the loan is to be paid off at the end of N years,
y(N ) = 0.
We want to derive an ODE satisfied by y. Let ∆t denote a small interval of time and ∆y
the change in the amount owed during the time interval ∆t. This change is determined by
• ∆y is increased by compounding at interest I; that is, ∆y is increased by the amount

Iy(t)∆t.
2.1. LINEAR FIRST ORDER EQUATIONS 19
• ∆y is decreased by the amount paid back in the time interval ∆t. If D denotes this
constant rate of payback, then D∆t is the amount paid back in the time interval ∆t.
Thus we have
∆y = Iy∆t − D∆t,
or
∆y
= Iy − D.
∆t
Letting ∆t → 0 we obtain the sought after ODE,
dy
= Iy − D. (2.8)
dt
This ODE is of form (2.3) with p = −I and g = −D. One immediately observes that this
ODE is not exactly what we assumed above, i.e. D is not known to us. Let us go ahead and
solve this equation for any constant D by the method of integrating factors. So we choose
µ according to (2.6),
t
µ(t) := exp p(s) ds
t
= exp − I ds
= exp(−It).
Applying (2.7) gives

t
1 c
y(t) = µ(s)g(s) ds +
µ(t) µ(t)
t
= eIt e−Is (−D) ds + ceIt

1 −It
= −De It
− e + ceIt
I
D
= + ceIt .
I
The constant c is fixed by requiring
y(0) = P,
that is
D
+ c = P.
I
Solving this for c gives c = P − D/I. Substituting this expression for c back into our solution
y(t) gives
D D
y(t) = − − P eIt .
I I
First observe that y(t) grows if D/I < P . (This might be a good definition of loan sharking!)
We have not yet determined D. To do so we use the condition that the loan is to be paid
off at the end of N years, y(N ) = 0. Substituting t = N into our solution y(t) and using
this condition gives
D D
0= − − P eN I .
I I
Solving for D,
eN I
D = PI , (2.9)
eN I − 1
gives the sought after relation between D, P , I and N . For example, if P = $100, 000,
I = 0.06 (6% interest) and the loan is for N = 30 years, then D = $7, 188.20 so the
monthly payment is D/12 = $599.02. Some years ago the mortgage rate was 12%. A quick
calculation shows that the monthly payment on the same loan at this interest would have
been $1028.09.
We remark that this model is a continuous model—the rate of payback is at the continuous
rate D. In fact, normally one pays back only monthly. Banks, therefore, might want to take
this into account in their calculations. I’ve found from personal experience that the above
model predicts the bank’s calculations to within a few dollars.
Suppose we increase our monthly payments by, say, $50. (We assume no prepayment
penalty.) This $50 goes then to paying off the principal. The problem then is how long does
it take to pay off the loan? It is an exercise to show that the number of years is (D is the
total payment in one year)
1 PI
− log 1 − . (2.10)
I D
Another questions asks on a loan of N years at interest I how long does it take to pay off
one-half of the principal? That is, we are asking for the time T when
P
y(T ) = .
2
It is an exercise to show that

1 1 NI
T = log (e + 1) . (2.11)
I 2
For example, a 30 year loan at 9% is half paid off in the 23rd year. Notice that T does not
depend upon the principal P .
2.2 Conservative systems
2.2.1 Energy conservation
Consider the motion of a particle of mass m in one dimension, i.e. the motion is along a
line. We suppose that the force acting at a point x, F (x), is conservative. This means there
exists a function V (x), called the potential energy, such that
dV
F (x) = − .
dx
2.2. CONSERVATIVE SYSTEMS 21
Amount owed
500 000
400 000
300 000
200 000
100 000
time in years
5 10 15 20 25 30
Loan of $500,000 for 30 years

Monthly payment
6000
5000
4000
3000
2000
1000
Interest rate
0.05 0.10 0.15
Figure 2.1: The top figure is the graph of the amount owned, y(t), as a function of time t
for a 30-year loan of $500,000 at interest rates 3%, 6%, 9% and 12%. The horizontal line in
the top figure is the line y = $250, 000; and hence, its intersection with the y(t)-curves gives
the time when the loan is half paid off. The lower the interest rate the lower the y(t)-curve.
The bottom figure gives the monthly payment on a 30-year loan of $500,000 as a function
of the interest rate I.
(Tradition has it we put in a minus sign.) In one dimension this requires that F is only a
function of x and not ẋ (= dx/dt) which physically means there is no friction. In higher
spatial dimensions the requirement that F is conservative is more stringent. The concept of
conservation of energy is that
E = Kinetic energy + Potential energy
does not change with time as the particle’s position and velocity evolves according to New-
ton’s equations. We now prove this fundamental fact. We recall from elementary physics
that the kinetic energy (KE) is given by
1
KE = mv 2 , v = velocity = ẋ.
2
Thus the energy is
2
1 dx
E = E(x, ẋ) = m + V (x).
2 dt
To show that E = E(x, ẋ) does not change with t when x = x(t) satisfies Newton’s equations,
we differentiate E with respect to t and show the result is zero:
dE dx d2 x dV dx
= m + (by the chain rule)
dt dt dt2 dx dt
dx d2 x dV (x)
= m 2 +
dt dt dx
2

dx d x
= m 2 − F (x) .
dt dt
Now not any function x = x(t) describes the motion of the particle—x(t) must satisfy
d2 x
F =m ,
dt2
and we now get the desired result
dE
= 0.
dt
This implies that E is constant on solutions to Newton’s equations.
We now use energy conservation and what we know about separation of variables to solve
the problem of the motion of a point particle in a potential V (x). Now
2
1 dx
E= m + V (x) (2.12)
2 dt
is a nonlinear first order differential equation. (We know it is nonlinear since the first
derivative is squared.) We rewrite the above equation as
2
dx 2
= (E − V (x)) ,
dt m
or
dx 2
=± (E − V (x)) .
dt m
(In what follows we take the + sign, but in specific applications one must keep in mind the
possibility that the − sign is the correct choice of the square root.) This last equation is of
the form in which we can separate variables. We do this to obtain
dx
= dt.
2
m (E − V (x))
This can be integrated to

1
± dx = t − t0 .
2
(E − V (x)) (2.13)
m
2.2.2 Hooke’s Law
Figure 2.2: Robert Hooke, 1635–1703.
Consider a particle of mass m subject to the force
F = −kx, k > 0, (Hooke’s Law). (2.14)
The minus sign (with k > 0) means the force is a restoring force—as in a spring. Indeed,
to a good approximation the force a spring exerts on a particle is given by Hooke’s Law. In
Figure 2.3: The mass-spring system: k is the spring constant in Hook’s Law, m is the mass
of the object and c represents a frictional force between the mass and floor. We neglect this
frictional force. (Later we’ll consider the effect of friction on the mass-spring system.)
this case x = x(t) measures the displacement from the equilibrium position at time t; and the
constant k is called the spring constant. Larger values of k correspond to a stiffer spring.
Newton’s equations are in this case
d2 x
m + kx = 0. (2.15)
dt2
This is a second order linear differential equation, the subject of the next chapter. However,
we can use the energy conservation principle to derive an associated nonlinear first order
equation as we discussed above. To do this, we first determine the potential corresponding
to Hooke’s force law.
One easily checks that the potential equals
1
V (x) = k x2 .
2
(This potential is called the harmonic potential.) Let’s substitute this particular V into
(2.13):
1
dx = t − t0 . (2.16)
2E/m − kx2 /m
Recall the indefinite integral

dx x
√ = arcsin + c.
a − x2
2 |a|
Using this in (2.16) we obtain

1 1 dx
dx =
2E/m − kx2 /m k/m 2E/k − x2

1 x
= arcsin + c.
k/m 2E/k
Thus (2.16) becomes3

x k
arcsin = t + c.
2E/k m
Taking the sine of both sides of this equation gives

x k
= sin t+c ,
2E/k m
or
2E k
x(t) = sin t+c . (2.17)
k m
Observe that there are two constants appearing in (2.17), E and c. Suppose one initial
condition is
x(0) = x0 .
Evaluating (2.17) at t = 0 gives
2E
x0 = sin(c). (2.18)
k
Now use the sine addition formula,
sin(θ1 + θ2 ) = sin θ1 cos θ2 + sin θ2 cos θ1 ,
in (2.17):

2E k k
x(t) = sin t cos c + cos t sin c
k m m

2E k k
= sin t cos c + x0 cos t (2.19)
k m m
where we use (2.18) to get the last equality.

Now substitute t = 0 into the energy conservation equation,
1 1 1
E= mv02 + V (x0 ) = mv02 + k x20 .
2 2 2
(v0 equals the velocity of the particle at time t = 0.) Substituting (2.18) in the right hand
side of this equation gives
1 1 2E
E = mv02 + k sin2 c
2 2 k
or
1
E(1 − sin2 c) = mv02 .
2
Recalling the trig identity sin2 θ + cos2 θ = 1, this last equation can be written as
1
E cos2 c = mv02 .
2
3 We use the same symbol c for yet another unknown constant.
Solve this for v0 to obtain the identity

2E
v0 = cos c.
m
We now use this in (2.19)

m k k
x(t) = v0 sin t + x0 cos t .
k m m
To summarize, we have eliminated the two constants E and c in favor of the constants x0
and v0 . As it must be, x(0) = x0 and ẋ(0) = v0 . The last equation is more easily interpreted
if we define
k
ω0 = . (2.20)
m
Observe that ω0 has the units of 1/time, i.e. frequency. Thus our final expression for the
position x = x(t) of a particle of mass m subject to Hooke’s Law is
v0
x(t) = x0 cos(ω0 t) + sin(ω0 t). (2.21)
ω0
Observe that this solution depends upon two arbitrary constants, x0 and v0 .4 In (2.7),
the general solution depended only upon one constant. It is a general fact that the number
of independent constants appearing in the general solution of a nth order5 ODE is n.
Period of mass-spring system satisfying Hooke’s Law
The sine and cosine are periodic functions of period 2π, i.e.
sin(θ + 2π) = sin θ, cos(θ + 2π) = cos θ.
This implies that our solution x = x(t) is periodic in time,
x(t + T ) = x(t),
where the period T is

2π m
T = = 2π . (2.23)
ω0 k
Observe that the period T , for the mass-spring system following Hooke’s law, depends on
the mass m and the spring constant k but not on the initial conditions .
4 ω is a constant too, but it is a parameter appearing in the differential equation that is fixed by the
0
mass m and the spring constant k. Observe that we can rewrite (2.15) as
ẍ + ω02 x = 0. (2.22)
Dimensionally this equation is pleasing: ẍ has the dimensions of d/t2 (d is distance and t is time) and so
does ω02 x since ω0 is a frequency. It is instructive to substitute (2.21) into (2.22) and verify directly that it
is a solution. Please do so!
5 The order of a scalar differential equation is equal to the order of the highest derivative appearing in
the equation. Thus (2.3) is first order whereas (2.15) is second order.
2.2.3 Period of the nonlinear pendulum
In this section we use the method of separation of variables to derive an exact formula for the
period of the pendulum. Recall that the ODE describing the time evolution of the angle of
deflection, θ, is (1.11). This ODE is a second order equation and so the method of separation
of variables does not apply to this equation. However, we will use energy conservation in a
manner similar to the previous section on Hooke’s Law.
To get some idea of what we should expect, first recall the approximation we derived for
small deflection angles, (1.12). Comparing this differential equation with (2.15), we see that
under the identification x → θ and m k
→ g , the two equations are identical. Thus using the
period derived in the last section, (2.23), we get as an approximation to the period of the
pendulum
2π
T0 = = 2π . (2.24)
ω0 g
An important feature of T0 is that it does not depend upon the amplitude of the oscillation.6
That is, suppose we have the initial conditions7
θ(0) = θ0 , θ̇(0) = 0, (2.25)
then T0 does not depend upon θ0 . We now proceed to derive our formula for the period, T ,
of the pendulum.
We claim that the energy of the pendulum is given by
1
E = E(θ, θ̇) = m2 θ̇2 + mg(1 − cos θ). (2.26)
2
Proof of (2.26)
We begin with
E = Kinetic energy + Potential energy

1
= mv 2 + mgy. (2.27)
2
(This last equality uses the fact that the potential at height h in a constant gravitational
force field is mgh. In the pendulum problem with our choice of coordinates h = y.) The x
and y coordinates of the pendulum ball are, in terms of the angle of deflection θ, given by
x = sin θ, y = (1 − cos θ).
Differentiating with respect to t gives
ẋ = cos θ θ̇, ẏ = sin θ θ̇,

6 Of course, its validity is only for small oscillations.
7 For simplicity we assume the initial angular velocity is zero, θ̇(0) = 0. This is the usual initial condition
for a pendulum.
from which it follows that the velocity is given by
v2 = ẋ2 + ẏ 2
= 2 θ̇2 .
Substituting these in (2.27) gives (2.26).

The energy conservation theorem states that for solutions θ(t) of (1.11), E(θ(t), θ̇(t)) is
independent of t. Thus we can evaluate E at t = 0 using the initial conditions (2.25) and
know that for subsequent t the value of E remains unchanged,
1
E = m2 θ̇(0)2 + mg (1 − cos θ(0))
2
= mg(1 − cos θ0 ).
Using this (2.26) becomes

1
mg(1 − cos θ0 ) = m2 θ̇2 + mg(1 − cos θ),
2
which can be rewritten as
1
m2 θ̇2 = mg(cos θ − cos θ0 ).
2
Solving for θ̇,
2g
θ̇ = (cos θ − cos θ0 ) ,

followed by separating variables gives
dθ
= dt. (2.28)
2g
(cos θ − cos θ0 )
We now integrate (2.28). The next step is a bit tricky—to choose the limits of integration
in such a way that the integral on the right hand side of (2.28) is related to the period T .
By the definition of the period, T is the time elapsed from t = 0 when θ = θ0 to the time T
when θ first returns to the point θ0 . By symmetry, T /2 is the time it takes the pendulum
to go from θ0 to −θ0 . Thus if we integrate the left hand side of (2.28) from −θ0 to θ0 the
time elapsed is T /2. That is,
θ0
1 dθ
T = .
2 −θ0 2g

Since the integrand is an even function of θ,
θ0
dθ
T =4 . (2.29)
2g
0
This is the sought after formula for the period of the pendulum. For small θ0 we expect
that T , as given by (2.29), should be approximately equal to T0 (see (2.24)). It is instructive
to see this precisely.
We now assume |θ0 | 1 so that the approximation
1 2 1
cos θ ≈ 1 − θ + θ4
2! 4!
is accurate for |θ| < θ0 . Using this approximation we see that
1 2 1
cos θ − cos θ0 ≈ (θ − θ2 ) − (θ04 − θ4 )
2! 0 4!
1 2 2 1 2 2
= (θ − θ ) 1 − (θ0 + θ ) .
2 0 12
From Taylor’s formula8 we get the approximation, valid for |x| 1,
1 1
√ ≈ 1 + x.
1−x 2
Thus

1 1 1
≈
2g
(cos θ − cos θ0 ) g θ02 − θ2 1 − 1 (θ2 + θ2 )
12 0

1 1 2 2
≈ 1 + (θ + θ ) .
g θ02 − θ2 24 0
Now substitute this approximate expression for the integrand appearing in (2.29) to find

T θ0 1 1 2 2
= 1 + (θ + θ ) + higher order corrections.
4 g 0 θ02 − θ2 24 0
Make the change of variables θ = θ0 x, then

θ0 1
dθ dx π
= √ = ,
θ 2 − θ2 1 − x2 2
0 0 0
θ0 1
θ2 dθ x2 dx π

2
= θ02 √ = θ02 .
θ0 − θ 2 1−x 2 4
0 0
Using these definite integrals we obtain

T π 1 2π 2π
= + (θ0 + θ0 )
4 g 2 24 2 4

π θ2
= 1 + 0 + higher order terms.
g 2 16
8 You should be able to do this without resorting to MatLab . But if you wanted higher order terms
MatLab would be helpful. Recall to do this we would enter

>> syms x
>> taylor(1/sqrt(1-x))
Recalling (2.24), we conclude

θ02
T = T0 1 + + ··· (2.30)
16
where the · · · represent the higher order correction terms coming from higher order terms
in the expansion of the cosines. These higher order terms will involve higher powers of θ0 .
It now follows from this last expression that
lim T = T0 .
θ0 →0
Observe that the first correction term to the linear result, T0 , depends upon the initial
amplitude of oscillation θ0 .
In Figure 2.4 shows the graph of the ratio T (θ0 )/T0 as a function of the initial displace-
ment angle θ0 .
Figure 2.4: Graph of the

the exact period T (θ0 ) of the pendulum divided by the linear
approximation T0 = 2π g as a function of the initial deflection angle θ0 . It can be proved
that as θ0 → π, the period T (θ0 ) diverges to +∞. Even so, the linear approximation is quite
good for moderate values of θ0 . For example at 45◦ (θ0 = π/4) the ratio is 1.03997. At 20◦
(θ0 = π/9) the ratio is 1.00767. The approximation (2.30) predicts for θ0 = π/9 the ratio
1.007632.
Remark: To use MatLab to evaluate symbolically these definite integrals you enter (note
the use of ’)
>> int(’1/sqrt(1-x^2)’,0,1)
and similarly for the second integral
>> int(’x^2/sqrt(1-x^2)’,0,1)
Numerical example
Suppose we have a pendulum of length = 1 meter. The linear theory says that the period
of the oscillation for such a pendulum is

1
T0 = 2π = 2π = 2.0071 sec.
g 9.8
If the amplitude of oscillation of the of the pendulum is θ0 ≈ 0.2 (this corresponds to roughly
a 20 cm deflection for the one meter pendulum), then (2.30) gives

1 2
T = T0 1 + (.2) = 2.0121076 sec.
16
One might think that these are so close that the correction is not needed. This might well
be true if we were interested in only a few oscillations. What would be the difference in one
week (1 week=604,800 sec)?
One might well ask how good an approximation is (2.30) to the exact result (2.29)? To
answer this we have to evaluate numerically the integral appearing in (2.29). Evaluating
(2.29) numerically (using say Mathematica’s NIntegrate) is a bit tricky because the end-
point θ0 is singular—an integrable singularity but it causes numerical integration routines
some difficulty. Here’s how you get around this problem. One isolates where the problem
occurs—near θ0 —and takes care of this analytically. For ε > 0 and ε 1 we decompose
the integral into two integrals: one over the interval (0, θ0 − ε) and the other one over the
interval (θ0 − ε, θ0 ). It’s the integral over this second interval that we estimate analytically.
Expanding the cosine function about the point θ0 , Taylor’s formula gives
cos θ0
cos θ = cos θ0 − sin θ0 (θ − θ0 ) − (θ − θ0 )2 + · · · .
2
Thus
1
cos θ − cos θ0 = sin θ0 (θ − θ0 ) 1 − cot θ0 (θ − θ0 ) + · · · .
2
So
1 1 1
√ = + ···
cos θ − cos θ0 sin θ0 (θ − θ0 ) 1 − 1 cot θ0 (θ0 − θ)
2

1 1
= 1 + cot θ0 (θ0 − θ) + · · ·
sin θ0 (θ0 − θ) 4
Thus
θ0 θ0
dθ dθ 1
√ = 1 + cot θ0 (θ − θ0 ) dθ + · · ·
θ0 −ε cos θ − cos θ0 θ0 −ε sin θ0 (θ0 − θ) 4
ε ε
1 1
= √ u−1/2 du + cot θ0 u1/2 du + · · · (u := θ0 − θ)
sin θ0 0 4 0

1 1
= √ 2ε1/2 + cot θ0 ε3/2 + · · · .
sin θ0 6
Choosing ε = 10−2 , the error we make in using the above expression is of order ε5/2 = 10−5 .
Substituting θ0 = 0.2 and ε = 10−2 into the above expression, we get the approximation
θ0
dθ
√ ≈ 0.4506
θ0 −ε cos θ − cos θ0
where we estimate the error lies in fifth decimal place. Now the numerical integration routine
in MatLab quickly evaluates this integral:
θ0 −ε
dθ
√ ≈ 1.7764
0 cos θ − cos θ0
for θ0 = 0.2 and ε = 10−2 . Specifically, one enters
>> quad(’1./sqrt(cos(x)-cos(0.2))’,0,0.2-1/100)
Hence for θ0 = 0.2 we have

θ0
dθ
√ ≈ 0.4506 + 1.77664 = 2.2270
0 cos θ − cos θ0
This implies
T ≈ 2.0121.
Thus the first order approximation (2.30) is accurate to some four decimal places when
θ0 ≤ 0.2. (The reason for such good accuracy is that the correction term to (2.30) is of order
θ04 .)
Remark: If you use MatLab to do the integral from 0 to θ0 directly, i.e.
>> quad(’1./sqrt(cos(x)-cos(0.2))’,0,0.2)
what happens? This is an excellent example of what may go wrong if one uses software
packages without thinking first ! Use help quad to find out more about numerical integration
in MatLab .
The attentive reader may have wondered how we produced the graph in Figure 2.4. It
turns out that the integral (2.29) can be expressed in terms of a special function called
“elliptic integral of the first kind”. The software Mathematica has this special function
and hence graphing it is easy to do: Just enter
Integrate[1/Sqrt[Cos[x]-Cos[x0]],{x,0,x0},Assumptions->{0<x0<Pi}]
to get the integral in terms of this special function. You can now ask Mathematica to plot
the result.
2.3. LEVEL CURVES OF THE ENERGY 33
2.3 Level curves of the energy
For the mass-spring system (Hooke’s Law) the energy is
1 1
E= mv 2 + kx2 (2.31)
2 2
which we can rewrite as
x 2 v 2
+ =1
a b

where a = 2E/k and b = 2E/m. We recognize this last equation as the equation of an
ellipse. Assuming k and m are fixed, we see that for various values of the energy E we get
different ellipses in the (x, v)-plane. Thus the values of x = x(t) and v = v(t) are fixed to
lie on various ellipses. The ellipse is fixed once we specify the energy E of the mass-spring
system.
For the pendulum the energy is
1
E= m2 ω 2 + mg(1 − cos θ) (2.32)
2
where ω = dθ/dt. What do the contour curves of (2.32) look like? That is we want the
curves in the (θ, ω)-plane that obey (2.32).
1
To make things simpler, we set 2 m2 = 1 and mg = 1 so that (2.32) becomes
E = ω 2 + (1 − cos θ) (2.33)
We now use Mathematica to plot the contour lines of (2.33) in the (θ, ω)-plane (see Fig-
ure 2.5). For small E the contour lines look roughly like ellipses but as E gets larger the
ellipses become more deformed. At E = 2 there is a curve that separates the deformed
elliptical curves from curves that are completely different (those contour lines corresponding
to E > 2). In terms of the pendulum what do you think happens when E > 2?
Figure 2.5: Contour lines for (2.33) for various values of the energy E.
2.4. EXERCISES 35
2.4 Exercises for Chapter 2

#1. Radioactive decay
Figure 2.6:
Carbon 14 is an unstable (radioactive) isotope of stable Carbon 12. If Q(t) represents

the amount of C14 at time t, then Q is known to satisfy the ODE
dQ
= −λQ
dt
where λ is a constant. If T1/2 denotes the half-life of C14 show that
log 2
T1/2 = .
λ
Recall that the half-life T1/2 is the time T1/2 such that Q(T1/2 ) = Q(0)/2. It is known for
C14 that T1/2 ≈ 5730 years. In Carbon 14 dating9 it becomes difficult to measure the levels
of C14 in a substance when it is of order 0.1% of that found in currently living material.
How many years must have passed for a sample of C14 to have decayed to 0.1% of its original
value? The technique of Carbon 14 dating is not so useful after this number of years.
9 From Wikipedia: The Earth’s atmosphere contains various isotopes of carbon, roughly in constant pro-
portions. These include the main stable isotope C12 and an unstable isotope C14. Through photosynthesis,
plants absorb both forms from carbon dioxide in the atmosphere. When an organism dies, it contains the
standard ratio of C14 to C12, but as the C14 decays with no possibility of replenishment, the proportion of
carbon 14 decreases at a known constant rate. The time taken for it to reduce by half is known as the half-life
of C14. The measurement of the remaining proportion of C14 in organic matter thus gives an estimate of
its age (a raw radiocarbon age). However, over time there are small fluctuations in the ratio of C14 to C12
in the atmosphere, fluctuations that have been noted in natural records of the past, such as sequences of
tree rings and cave deposits. These records allow fine-tuning, or “calibration”, of the raw radiocarbon age,
to give a more accurate estimate of the calendar date of the material. One of the most frequent uses of
radiocarbon dating is to estimate the age of organic remains from archaeological sites. The concentration
of C14 in the atmosphere might be expected to reduce over thousands of years. However, C14 is constantly
being produced in the lower stratosphere and upper troposphere by cosmic rays, which generate neutrons
that in turn create C14 when they strike nitrogen–14 atoms. Once produced, the C14 quickly combines with
the oxygen in the atmosphere to form carbon dioxide. Carbon dioxide produced in this way diffuses in the
atmosphere, is dissolved in the ocean, and is taken up by plants via photosynthesis. Animals eat the plants,
and ultimately the radiocarbon is distributed throughout the biosphere.
#2: Mortgage payment problem
In the problem dealing with mortgage rates, prove (2.10) and (2.11). Using either a hand
calculator or some computer software, create a table of monthly payments on a loan of
$200,000 for 30 years for interest rates from 1% to 15% in increments of 1%.
#3: Discontinuous forcing term
Solve
y + 2y = g(t), y(0) = 0,
where
1, 0 ≤ t ≤ 1
g(t) =
0, t>1
We make the additional assumption that the solution y = y(t) should be a continuous
function of t. Hint: First solve the differential equation on the interval [0, 1] and then on
the interval [1, ∞). You are given the initial value at t = 0 and after you solve the equation
on [0, 1] you will then know y(1).10 Plot the solution y = y(t) for 0 ≤ t ≤ 4. (You can use
any computer algebra program or just graph the y(t) by hand.)
#4. Application to population dynamics
In biological applications the population P of certain organisms at time t is sometimes

assumed to obey the equation
dP P
= aP 1 − (2.34)
dt E
where a and E are positive constants. This model is sometimes called the logistic growth
model.
1. Find the equilibrium solutions. (That is solutions that don’t change with t.)
2. From (2.34) determine the regions of P where P is increasing (decreasing) as a function

of t. Again using (2.34) find an expression for d2 P/dt2 in terms of P and the constants
a and E. From this expression find the regions of P where P is convex (d2 P/dt2 > 0)
and the regions where P is concave (d2 P/dt2 < 0).
3. Using the method of separation of variables solve (2.34) for P = P (t) assuming that
at t = 0, P = P0 > 0. Find
lim P (t)
t→∞
Hint: To do the integration first use the identity

1 1 1
= +
P (1 − P/E) P E−P
4. Sketch P as a function of t for 0 < P0 < E and for E < P0 < ∞.

10 This is problem #32, pg. 74 (7th edition) of the Boyce & DiPrima [4].
2.4. EXERCISES 37
#5: Mass-spring system with friction
We reconsider the mass-spring system but now assume there is a frictional force present and
this frictional force is proportional to the velocity of the particle. Thus the force acting on
the particle comes from two terms: one due to the force exerted by the spring and the other
due to the frictional force. Thus Newton’s equations become
−kx − β ẋ = mẍ (2.35)

where as before x = x(t) is the displacement from the equilibrium position at time t. β and
k are positive constants. Introduce the energy function
1 1
E = E(x, ẋ) = mẋ2 + kx2 , (2.36)
2 2
and show that if x = x(t) satisfies (2.35), then
dE
< 0.
dt
What is the physical meaning of this last inequality?
#6: Nonlinear mass-spring system
Consider a mass-spring system where x = x(t) denotes the displacement of the mass m
from its equilibrium position at time t. The linear spring (Hooke’s Law) assumes the force
exerted by the spring on the mass is given by (2.14). Suppose instead that the force F is
given by
F = F (x) = −kx − ε x3 (2.37)
where ε is a small positive number.11 The second term represents a nonlinear correction
to Hooke’s Law. Why is it reasonable to assume that the first correction term to Hooke’s
Law is of order x3 and not x2 ? (Hint: Why is it reasonable to assume F (x) is an odd
function of x?) Using the solution for the period of the pendulum as a guide, find an exact
integral expression for the period T of this nonlinear mass-spring system assuming the initial
conditions
dx
x(0) = x0 , (0) = 0.
dt
Define
εx2
z = 0.
2k
Show that z is dimensionless and that your expression for the period T can be written as
1
4 1
T = √ du (2.38)
ω0 0 1 − u + z − zu4
2

where ω0 = k/m. We now assume that z 1. (This is the precise meaning of the
parameter ε being small.) Taylor expand the function
1
√
1 − u + z − zu4
2
11 One could also consider ε < 0. The case ε > 0 is a called a hard spring and ε < 0 a soft spring.
in the variable z to first order. You should find

1 1 1 + u2
√ = √ − √ z + O(z 2 ).
1 − u2 + z − zu4 1 − u2 2 1 − u2
Now use this approximate expression in the integrand of (2.38), evaluate the definite integrals
that arise, and show that the period T has the Taylor expansion

2π 3 2
T = 1 − z + O(z ) .
ω0 4
#7: Motion in a central field
A (three-dimensional) force F is called a central force 12 if the direction of F lies along the
the direction of the position vector r. This problem asks you to show that the motion of a
particle in a central force, satisfying
2
d r
F = m 2 , (2.39)
dt
lies in a plane.
1. Show that
:= r × p with p := mv
M (2.40)
is constant in t for r = r(t) satisfying (2.39). (Here v = dr/dt is the velocity vector
and p is the momentum vector. In words, the momentum vector is mass times the
velocity vector.) The × in (2.40) is the vector cross product. Recall (and you may
assume this result) from vector calculus that
d da db
(a × b) = × b + a × .
dt dt dt
is called the angular momentum vector.
The vector M
2. From the fact that M is a constant vector, show that the vector r(t) lies in a plane
. Also you may find helpful the vector identity
perpendicular to M . Hint: Look at r · M
a · (b × c) = b · (c × a) = c · (a × b).
#8: Motion in a central field (cont)
From the preceding problem we learned that the position vector r(t) for a particle moving
in a central force lies in a plane. In this plane, let (r, θ) be the polar coordinates of the point
r, i.e.
x(t) = r(t) cos θ(t), y(t) = r(t) sin θ(t) (2.41)
12 For an in depth treatment of motion in a central field, see [1], Chapter 2, §8.
2.4. EXERCISES 39
1. In components, Newton’s equations can be written (why?)

x y
Fx = f (r) = mẍ, Fy = f (r) = mÿ (2.42)
r r
where f (r) is the magnitude of the force F . By twice differentiating (2.41) with
respect to t, derive formulas for ẍ and ÿ in terms of r, θ and their derivatives. Use
these formulas in (2.42) to show that Newton’s equations in polar coordinates (and
for a central force) become
1
f (r) cos θ = r̈ cos θ − 2ṙθ̇ sin θ − rθ̇2 cos θ − rθ̈ sin θ, (2.43)
m
1
f (r) sin θ = r̈ sin θ + 2ṙθ̇ cos θ − rθ̇2 sin θ + rθ̈ cos θ. (2.44)
m
Multiply (2.43) by cos θ, (2.44) by sin θ, and add the resulting two equations to show
that
1
r̈ − rθ̇2 = f (r). (2.45)
m
Now multiply (2.43) by sin θ, (2.44) by cos θ, and substract the resulting two equations
to show that
2ṙ θ̇ + rθ̈ = 0. (2.46)
Observe that the left hand side of (2.46) is equal to
1 d 2
(r θ̇).
r dt
Using this observation we then conclude (why?)
r2 θ̇ = H (2.47)
for some constant H. Use (2.47) to solve for θ̇, eliminate θ̇ in (2.45) to conclude that
the polar coordinate function r = r(t) satisfies
1 H2
r̈ = f (r) + 3 . (2.48)
m r
2. Equation (2.48) is of the form that a second derivative of the unknown r is equal to
some function of r. We can thus apply our general energy method to this equation.
Let Φ be a function of r satisfying
1 dΦ
f (r) = − ,
m dr
and find an effective potential V = V (r) such that (2.48) can be written as
dV
r̈ = − (2.49)
dr
H2
(Ans: V (r) = Φ(r) + 2r 2 ). Remark: The most famous choice for f (r) is the inverse
square law
mM G0
f (r) = −
r2
which describes the gravitational attraction of two particles of masses m and M . (G0
is the universal gravitational constant.) In your physics courses, this case will be
analyzed in great detail. The starting point is what we have done here.
3. With the choice

mM G0
f (r) = −
r2
the equation (2.48) gives a DE that is satisfied by r as a function of t:
G H2
r̈ = − 2
+ 3 (2.50)
r r
where G = M G0 . We now use (2.50) to obtain a DE that is satisfied by r as a function
of θ. This is the quantity of interest if one wants the orbit of the planet. Assume that
H = 0, r = 0, and set r = r(θ). First, show that by chain rule
r̈ = r θ̇2 + r θ̈. (2.51)
(Here, implies the differentiation with respect to θ, and as usual, the dot refers to
differentiation with respect to time.) Then use (2.47) and (2.51) to obtain
H2 2 2H
2
r̈ = r − (r ) (2.52)
r4 r5
Now, obtain a second order DE of r as a function of θ from (2.50) and (2.52). Finally,
by letting u(θ) = 1/r(θ), obtain a simple linear constant coefficient DE
G
u + u = (2.53)
H2
which is known as Binet’s equation.13
#9: Euler’s equations for a rigid body with no torque
In mechanics one studies the motion of a rigid body14 around a stationary point in the
absence of outside forces. Euler’s equations are differential equations for the angular velocity
vector Ω = (Ω1 , Ω2 , Ω3 ). If Ii denotes the moment of inertia of the body with respect to the
ith principal axis, then Euler’s equations are
dΩ1
I1 = (I2 − I3 )Ω2 Ω3
dt
dΩ2
I2 = (I3 − I1 )Ω3 Ω1
dt
dΩ3
I3 = (I1 − I2 )Ω1 Ω2
dt
Prove that
M = I12 Ω21 + I22 Ω22 + I32 Ω23
and
1 1 1
I1 Ω21 + I2 Ω22 + I3 Ω23
E=
2 2 2
are both first integrals of the motion. (That is, if the Ωj evolve according to Euler’s equa-
tions, then M and E are independent of t.)
13 For further discussion of Binet’s equation see [8].
14 For an in-depth discussion of rigid body motion see Chapter 6 of [1].
2.4. EXERCISES 41
#10. Exponential function
In calculus one defines the exponential function et by

t n
et := lim (1 + ) , t ∈ R.
n→∞ n
Suppose one took the point of view of differential equations and defined et to be the (unique)
solution to the ODE
dE
=E (2.54)
dt
that satisfies the initial condition E(0) = 1.15 Prove that the addition formula
et+s = et es
follows from the ODE definition. [Hint: Define
φ(t) := E(t + s) − E(t)E(s)
where E(t) is the above unique solution to the ODE satisfying E(0) = 1. Show that φ
satisfies the ODE
dφ
= φ(t)
dt
From this conclude that necessarily φ(t) = 0 for all t.]
Using the above ODE definition of E(t) show that
t
E(s) ds = E(t) − 1.
0
Let E0 (t) = 1 and define En (t), n ≥ 1 by

t
En+1 (t) = 1 + En (s) ds, n = 0, 1, 2, . . . . (2.55)
0
Show that
t2 tn
+ ···+ .
En (t) = 1 + t +
2! n!
By the ratio test this sequence of partial sums converges as n → ∞. Assuming one can take
the limit n → ∞ inside the integral (2.55),16 conclude that
∞ n
t
et = E(t) =
n=0
n!
15 That is, we are taking the point of view that we define et to be the solution E(t). Here is a proof that
given a solution to (2.54) satisfying the initial condition E(0) = 1, that such a solution is unique. Suppose
we have found two such solutions: E1 (t) and E2 (t). Let y(t) = E1 (t)/E2 (t), then
dy 1 dE1 E1 dE2 E1 E1
= − 2 = − 2 E2 = 0
dt E2 dt E2 dt E2 E2
Thus y(t) = constant.
But we know that y(0) = E1 (0)/E2 (0) = 1. Thus y(t) = 1, or E1 (t) = E2 (t).
16 The series n
n≥0 s /n! converges uniformly on the closed interval [0, t]. From this fact it follows that
one is allowed to interchange the sum and integration. These convergence topics are normally discussed in
an advanced calculus course.
#11. Addition formula for the tangent function
Suppose we wish to find a real-valued, differentiable function F (x) that satisfies the func-
tional equation
F (x) + F (y)
F (x + y) = (2.56)
1 − F (x)F (y)
1. Show that such an F necessarily satisfies F (0) = 0. Hint: Use (2.56) to get an
expression for F (0 + 0) and then use fact that we seek F to be real-valued.
2. Set α = F (0). Show that F must satisfy the differential equation
dF
= α(1 + F (x)2 ) (2.57)
dx
Hint: Differentiate (2.56) with respect to y and then set y = 0.
3. Use the method of separation of variables to solve (2.57) and show that
F (x) = tan(αx).
Chapter 3
Second Order Linear Equations
Figure 3.1: eix = cos +i sin x, Leonhard Euler, Introductio in Analysin Infinitorum, 1748
43
44 CHAPTER 3. SECOND ORDER LINEAR EQUATIONS
3.1 Theory of second order equations
3.1.1 Vector space of solutions
First order linear differential equations are of the form
dy
+ p(x)y = f (x). (3.1)
dx
Second order linear differential equations are linear differential equations whose highest
derivative is second order:
d2 y dy
2
+ p(x) + q(x)y = f (x). (3.2)
dx dx
If f (x) = 0,
d2 y dy
2
+ p(x) + q(x)y = 0, (3.3)
dx dx
the equation is called homogeneous. For the discussion here, we assume p and q are contin-
uous functions on a closed interval [a, b]. There are many important examples where this
condition fails and the points at which either p or q fail to be continuous are called singular
points. An introduction to singular points in ordinary differential equations can be found
in Boyce & DiPrima [4]. Here are some important examples where the continuity condition
fails.
Legendre’s equation
2x n(n + 1)
p(x) = − , q(x) = .
1 − x2 1 − x2
At the points x = ±1 both p and q fail to be continuous.
Bessel’s equation
1 ν2
, q(x) = 1 − 2 .
p(x) =
x x
At the point x = 0 both p and q fail to be continuous.
We saw that a solution to (3.1) was uniquely specified once we gave one initial condition,
y(x0 ) = y0 .
In the case of second order equations we must give two initial conditions to specify uniquely
a solution:
y(x0 ) = y0 and y (x0 ) = y1 . (3.4)
This is a basic theorem of the subject. It says that if p and q are continuous on some
interval (a, b) and a < x0 < b, then there exists an unique solution to (3.3) satisfying the
initial conditions (3.4).1 We will not prove this theorem in this class. As an example of the
1 See Theorem 3.2.1 in the [4], pg. 131 or chapter 6 of [3]. These theorems dealing with the existence and
uniqueness of the initial value problem are covered in an advanced course in differential equations.
3.1. THEORY OF SECOND ORDER EQUATIONS 45
appearance to two constants in the general solution, recall that the solution of the harmonic
oscillator
ẍ + ω02 x = 0
contained x0 and v0 .
Let V denote the set of all solutions to (3.3). The most important feature of V is that
it is a two-dimensional vector space. That it is a vector space follows from the linearity of
(3.3). (If y1 and y2 are solutions to (3.3), then so is c1 y1 + c2 y2 for all constants c1 and c2 .)
To prove that the dimension of V is two, we first introduce two special solutions. Let Y1
and Y2 be the unique solutions to (3.3) that satisfy the initial conditions
Y1 (0) = 1, Y1 (0) = 0, and Y2 (0) = 0, Y2 (0) = 1,
respectively.
We claim that {Y1 , Y2 } forms a basis for V. To see this let y(x) be any solution to (3.3).2
Let c1 := y(0), c2 := y (0) and
∆(x) := y(x) − c1 Y1 (x) − c2 Y2 (x).
Since y, Y1 and Y2 are solutions to (3.3), so too is ∆. (V is a vector space.) Observe
∆(0) = 0 and ∆ (0) = 0. (3.5)
Now the function y0 (x) :≡ 0 satisfies (3.3) and the initial conditions (3.5). Since solutions
are unique, it follows that ∆(x) ≡ y0 ≡ 0. That is,
y = c1 Y1 + c2 Y2 .
To summarize, we’ve shown every solution to (3.3) is a linear combination of Y1 and Y2 .

That Y1 and Y2 are linearly independent follows from their initial values: Suppose
c1 Y1 (x) + c2 Y2 (x) = 0.
Evaluate this at x = 0, use the initial conditions to see that c1 = 0. Take the derivative of
this equation, evaluate the resulting equation at x = 0 to see that c2 = 0. Thus, Y1 and Y2
are linearly independent. We conclude, therefore, that {Y1 , Y2 } is a basis and dimV = 2.
3.1.2 Wronskians
Given two solutions y1 and y2 of (3.3) it is useful to find a simple condition that tests
whether they form a basis of V. Let ϕ be the solution of (3.3) satisfying ϕ(x0 ) = ϕ0 and
ϕ (x0 ) = ϕ1 . We ask are there constants c1 and c2 such that
ϕ(x) = c1 y1 (x) + c2 y2 (x)
for all x? A necessary and sufficient condition that such constants exist at x = x0 is that
the equations
ϕ0 = c1 y1 (x0 ) + c2 y2 (x0 ),
ϕ1 = c1 y (x0 ) + c2 y2 (x0 ),
2 We assume for convenience that x = 0 lies in the interval (a, b).
have a unique solution {c1 , c2 }. From linear algebra we know this holds if and only if the
determinant
y1 (x0 ) y2 (x0 )

y1 (x0 ) y2 (x0 ) = 0.
We define the Wronskian of two solutions y1 and y2 of (3.3) to be

y (x) y2 (x)
W (y1 , y2 ; x) := 1 = y1 (x)y (x) − y (x)y2 (x). (3.6)
y1 (x) y2 (x). 2 1
From what we have said so far one would have to check that W (y1 , y2 ; x) = 0 for all x to
conclude {y1 , y2 } forms a basis.
We now derive a formula for the Wronskian that will make the check necessary at only
one point. Since y1 and y2 are solutions of (3.3), we have
y1 + p(x)y1 + q(x)y1 = 0, (3.7)

y2 + p(x)y2 + q(x)y2 = 0. (3.8)
Now multiply (3.7) by y2 and multiply (3.8) by y1 . Subtract the resulting two equations to
obtain
y1 y2 − y1 y2 + p(x) (y1 y2 − y1 y2 ) = 0. (3.9)
Recall the definition (3.6) and observe that
dW
= y1 y2 − y1 y2 .
dx
Hence (3.9) is the equation
dW
+ p(x)W (x) = 0, (3.10)
dx
whose solution is x
W (y1 , y2 ; x) = c exp − p(s) dx . (3.11)
Since the exponential is never zero we see from (3.11) that either W (y1 , y2 ; x) ≡ 0 or
W (y1 , y2 ; x) is never zero.
To summarize, to determine if {y1 , y2 } forms a basis for V, one needs to check at only
one point whether the Wronskian is zero or not.
Applications of Wronskians
1. Claim: Suppose {y1 , y2 } form a basis of V, then they cannot have a common point of
inflection in a < x < b unless p(x) and q(x) simultaneously vanish there. To prove
this, suppose x0 is a common point of inflection of y1 and y2 . That is,
y1 (x0 ) = 0 and y2 (x0 ) = 0.
Evaluating the differential equation (3.3) satisfied by both y1 and y2 at x = x0 gives
p(x0 )y1 (x0 ) + q(x0 )y1 (x0 ) = 0,

p(x0 )y2 (x0 ) + q(x0 )y2 (x0 ) = 0.
3.1. THEORY OF SECOND ORDER EQUATIONS 47
Assuming that p(x0 ) and q(x0 ) are not both zero at x0 , the above equations are a set
of homogeneous equations for p(x0 ) and q(x0 ). The only way these equations can have
a nontrivial solution is for the determinant

y1 (x0 ) y1 (x0 )
= 0.
y2 (x0 ) y2 (x0 )
That is, W (y1 , y2 ; x0 ) = 0. But this contradicts that {y1 , y2 } forms a basis. Thus
there can exist no such common inflection point.
2. Claim: Suppose {y1 , y2 } form a basis of V and that y1 has consecutive zeros at x = x1
and x = x2 . Then y2 has one and only one zero between x1 and x2 . To prove this we
first evaluate the Wronskian at x = x1 ,
W (y1 , y2 ; x1 ) = y1 (x1 )y2 (x1 ) − y1 (x1 )y2 (x1 ) = −y1 (x1 )y2 (x1 )
since y1 (x1 ) = 0. Evaluating the Wronskian at x = x2 gives
W (y1 , y2 ; x2 ) = −y1 (x2 )y2 (x2 ).
Now W (y1 , y2 ; x1 ) is either positive or negative. (It can’t be zero.) Let’s assume it
is positive. (The case when the Wronskian is negative is handled similarly. We leave
this case to the reader.) Since the Wronskian is always of the same sign, W (y1 , y2 ; x2 )
is also positive. Since x1 and x2 are consecutive zeros, the signs of y1 (x1 ) and y1 (x2 )
are opposite of each other. But this implies (from knowing that the two Wronskian
expressions are both positive), that y2 (x1 ) and y2 (x2 ) have opposite signs. Thus there
exists at least one zero of y2 at x = x3 , x1 < x3 < x2 . If there exist two or more such
zeros, then between any two of these zeros apply the above argument (with the roles
of y1 and y2 reversed) to conclude that y1 has a zero between x1 and x2 . But x1 and
x2 were assumed to be consecutive zeros. Thus y2 has one and only one zero between
x1 and x2 .
In the case of the harmonic oscillator, y1 (x) = cos ω0 x and y2 (x) = sin ω0 x, and the
fact that the zeros of the sine function interlace those of the cosine function is well
known.
Here is a second example: Consider the Airy differential equation
d2 y
− xy = 0 (3.12)
dx2
Two linearly independent solutions to the Airy DE are plotted in Figure 3.2. We denote
these particular linearly independent solutions by y1 (x) := Ai(x) and y2 (x) := Bi(x).
The function Ai(x) is the solution approaching zero as x → +∞ in Figure 3.2. Note
the interlacing of the zeros.
Figure 3.2: The Airy functions Ai(x) and Bi(x) are plotted. Note that the between any two
zeros of one solutions lies a zero of the other solution.
3.2 Reduction of order

Suppose y1 is a solution of (3.3). Let
y(x) = v(x)y1 (x).
Then
y = v y1 + vy1 and y = v y1 + 2v y1 + vy1 .
Substitute these expressions for y and its first and second derivatives into (3.3) and make
use of the fact that y1 is a solution of (3.3). One obtains the following differential equation
for v:
y
v + p + 2 1 v = 0,
y1
or upon setting u = v ,
y1
u + p+2 u = 0.
y1
This last equation is a first order linear equation. Its solution is

y c
u(x) = c exp − p + 2 1 dx = 2 exp − p(x) dx .
y1 y1 (x)
This implies
v(x) = u(x) dx,
3.3. CONSTANT COEFFICIENTS 49
so that
y(x) = cy1 (x) u(x) dx.
The point is, we have shown that if one solution to (3.3) is known, then a second solution
can be found—expressed as an integral.
3.3 Constant coefficients

We assume that p(x) and q(x) are constants independent of x. We write (3.3) in this case
as3
ay + by + cy = 0. (3.13)
We “guess” a solution of the form
y(x) = eλx .
Substituting this into (3.13) gives
aλ2 eλx + bλeλx + ceλx = 0.
Since eλx is never zero, the only way the above equation can be satisfied is if
aλ2 + bλ + c = 0. (3.14)
Let λ± denote the roots of this quadratic equation, i.e.

√
−b ± b2 − 4ac
λ± = .
2a
We consider three cases.
1. Assume b2 − 4ac > 0 so that the roots λ± are both real numbers. Then exp(λ+ x) and
exp(λ− x) are two linearly independent solutions to (3.14). That they are solutions
follows from their construction. They are linearly independent since
W (eλ+ x , eλ− x ; x) = (λ− − λ+ )eλ+ x eλ− x = 0
Thus in this case, every solution of (3.13) is of the form
c1 exp(λ+ x) + c2 exp(λ− x)
for some constants c1 and c2 .
2. Assume b2 − 4ac = 0. In this case λ+ = λ− . Let λ denote their common value. Thus
we have one solution y1 (x) = eλx . We could use the method of reduction of order
to show that a second linearly independent solution is y2 (x) = xeλx . However, we
choose to present a more intuitive way of seeing this is a second linearly independent
solution. (One can always make it rigorous at the end by verifying that that it is
3 This corresponds to p(x) = b/a and q(x) = c/a. For applications it is convenient to introduce the
constant a.
indeed a solution.) Suppose we are in the distinct root case but that the two roots are
very close in value: λ+ = λ + ε and λ− = λ. Choosing c1 = −c2 = 1/ε, we know that
1 (λ+ε)x 1 λx
c1 y 1 + c2 y 2 = e − e
ε ε
λx e
εx
−1
= e
ε
is also a solution. Letting ε → 0 one easily checks that
eεx − 1
→ x,
ε
so that the above solution tends to
xeλx ,
our second solution. That {eλx , xeλx } is a basis is a simple Wronskian calculation.
3. We assume b2 − 4ac < 0. In this case the roots λ± are complex. At this point we
review the the exponential of a complex number.
Complex exponentials
Let z = x + iy (x, y real numbers, i2 = −1) be a complex number. Recall that x is

called the real part of z, z, and y is called the imaginary part of z, z. Just as we
picture real numbers as points lying in a line, called the real line R; we picture complex
numbers as points lying in the plane, called the complex plane C. The coordinates
of z in the complex plane are (x, y). The absolute value of z, denoted |z|, is equal to

x2 + y 2 . The complex conjugate of z, denoted z, is equal to x − iy. Note the useful
relation
2
z z = |z| .
In calculus, or certainly an advanced calculus class, one considers (simple) functions
of a complex variable. For example the function
f (z) = z 2
takes a complex number, z, and returns it square, again a complex number. (Can you
show that f = x2 − y 2 and f = 2xy?). Using complex addition and multiplication,
one can define polynomials of a complex variable
an z n + an−1 z n−1 + · · · + a1 z + a0 .
The next (big) step is to study power series

∞

an z n .
n=0
With power series come issues of convergence. We defer these to your advanced calculus
class.
With this as a background we are (almost) ready to define the exponential of a complex
number z. First, we recall that the exponential of a real number x has the power series
expansion
∞
xn
ex = exp(x) = (0! := 1).
n=0
n!
In calculus classes, one normally defines the exponential in a different way4 and then
proves ex has this Taylor expansion. However, one could define the exponential func-
tion by the above formula and then prove the various properties of ex follow from this
definition. This is the approach we take for defining the exponential of a complex
number except now we use a power series in a complex variable:5
∞
zn
ez = exp(z) := , z∈C (3.15)
n=0
n!
We now derive some properties of exp(z) based upon this definition.

• Let θ ∈ R, then
∞
(iθ)n
exp(iθ) =
n=0
n!
∞ ∞
(iθ)2n (iθ)2n+1
= +
n=0
(2n)! n=0
(2n + 1)!
∞
∞
θ2n θ2n+1
= (−1)n +i (−1)n
n=0
(2n)! n=0
(2n + 1)!
= cos θ + i sin θ.
This last formula is called Euler’s Formula. See Figure 3.3. Two immediate
consequences of Euler’s formula (and the facts cos(−θ) = cos θ and sin(θ) =
− sin θ) are
exp(−iθ) = cos θ − i sin θ

exp(iθ) = exp(−iθ)
Hence
2
|exp(iθ)| = exp(iθ) exp(−iθ) = cos2 θ + sin2 θ = 1
That is, the values of exp(iθ) lie on the unit circle. The coordinates of the point
eiθ are (cos θ, sin θ).
• We claim the addition formula for the exponential function, well-known for real
values, also holds for complex values
4A common definition is ex = limn→∞ (1 + x/n)n .
5 It can be proved that this infinite series converges for all complex values z.
Figure 3.3: Euler’s formula.
exp(z + w) = exp(z) exp(w), z, w ∈ C. (3.16)
We are to show
∞
1
exp(z + w) = (z + w)n
n=0
n!
∞ n
1 n k n−k
= z w (binomial theorem)
n=0
n! k
k=0
is equal to
∞
∞
1 k 1 m
exp(z) exp(w) = z w
k! m=0 m!
k=0
∞
1 k m
= z w
k!m!
k,m=0
∞ n
1
= z k wn−k n := k + m
n=0 k=0
k!(n − k)!
∞
n
1 n!
= z k wn−k .
n=0
n! k!(n − k)!
k=0
Since
n n!
= ,
k k!(n − k)!
we see the two expressions are equal as claimed.

• We can now use these two properties to understand better exp(z). Let z = x+ iy,
then
exp(z) = exp(x + iy) = exp(x) exp(iy) = ex (cos y + i sin y) .
Observe the right hand side consists of functions from calculus. Thus with a
calculator you could find the exponential of any complex number using this for-
mula.6
A form of the complex exponential we frequently use is if λ = σ + iµ and x ∈ R,
then
exp(λx) = exp ((σ + iµ)x)) = eσx (cos(µx) + i sin(µx)) .
Returning to (3.13) in case b2 − 4ac < 0 and assuming a, b and c are all real, we see
that the roots λ± are of the form7
λ+ = σ + iµ and λ− = σ − iµ.
Thus eλ+ x and eλ− x are linear combinations of
eσx cos(µx) and eσx sin(µx).
That they are linear independent follows from a Wronskian calculuation. To summa-
rize, we have shown that every solution of (3.13) in the case b2 − 4ac < 0 is of the
form
c1 eσx cos(µx) + c2 eσx sin(µx)
for some constants c1 and c2 .
Remarks: The MatLab function exp handles complex numbers. For example,
>> exp(i*pi)
ans =
-1.0000 + 0.0000i
The imaginary unit i is i in MatLab . You can also use sqrt(-1) in place of i. This is
sometimes useful when i is being used for other purposes. There are also the functions
abs, angle, conj, imag real
For example,
>> w=1+2*i
w =
1.0000 + 2.0000i
6 Of course, this assumes x

7σ
√ your calculator doesn’t overflow or underflow in computing e .
2
= −b/2a and µ = 4ac − b /2a.
>> abs(w)
ans =
2.2361
>> conj(w)
ans =
1.0000 - 2.0000i
>> real(w)
ans =
>> imag(w)
ans =
>> angle(w)
ans =
1.1071
3.4 Forced oscillations of the mass-spring system

The forced mass-spring system is described by the differential equation
d2 x dx
m +γ + k x = F (t) (3.17)
dt2 dt
where x = x(t) is the displacement from equilibrium at time t, m is the mass, k is the
constant in Hooke’s Law, γ > 0 is the coefficient of friction, and F (t) is the forcing term. In
these notes we examine the solution when the forcing term is periodic with period 2π/ω. (ω
is the frequency of the forcing term.) The simplest choice for a periodic function is either
sine or cosine. Here we examine the choice
F (t) = F0 cos ωt
where F0 is the amplitude of the forcing term. All solutions to (3.17) are of the form
x(t) = xp (t) + c1 x1 (t) + c2 x2 (t) (3.18)

3.4. FORCED OSCILLATIONS OF THE MASS-SPRING SYSTEM 55
where xp is a particular solution of (3.17) and {x1 , x2 } is a basis for the solution space of
the homogeneous equation.
The homogeneous solutions have been discussed earlier. We know that both x1 and x2
will contain a factor
e−(γ/2m)t
times factors involving sine and cosine. Since for all a > 0, e−at → 0 as t → ∞, the
homogeneous part of (3.18) will tend to zero. That is, for all initial conditions we have for
large t to good approximation
x(t) ≈ xp (t).
Thus we concentrate on finding a particular solution xp .
With the right-hand side of (3.17) having a cosine term, it is natural to guess that the
particular solution will also involve cos ωt. If one guesses
A cos ωt
one quickly sees that due to the presence of the frictional term, this cannot be a correct
since sine terms also appear. Thus we guess
xp (t) = A cos ωt + B sin ωt (3.19)
We calculate the first and second dervatives of (3.19) and substitute the results together
with (3.19) into (3.17). One obtains the equation

−Aω 2 m + Bωγ + kA cos ωt + −Bω 2 m − Aωγ + kB sin ωt = F0 cos ωt
This equation must hold for all t and this can happen only if

−Aω 2 m + Bωγ + kA = F0 and −Bω 2 m − Aωγ + kB = 0
These last two equations are a pair of linear equations for the unknown coefficients A and B.
We now solve these linear equations. First we rewrite these equations to make subsequent
steps clearer:

k − ω 2 m A + ωγ B = F0 ,

−ωγ A + k − ω 2 m B = 0.
Using Cramer’s Rule we find (check this!)
k − mω 2
A = F0
(k − mω 2 )2 + γ 2 ω 2
γω
B = F0
(k − mω 2 )2 + γ 2 ω 2
We can make these results notationally simpler if we recall that the natural frequency of a
(frictionless) oscillator is
k
ω02 =
m
and define
∆(ω) = m2 (ω 2 − ω02 )2 + γ 2 ω 2 (3.20)
so that
m(ω02 − ω 2 ) γω
A= 2
F0 and B = F0
∆(ω) ∆(ω)2
Using these expressions for A and B we can substitute into (3.19) to find our particular
solution xp . The form (3.19) is not the best form in which to understand the properties of
the solution. (It is convenient for performing the above calculations.) For example, it is not
obvious from (3.19) what is the amplitude of oscillation. To answer this and other questions
we introduce polar coordinates for A and B:
A = R cos δ and B = R sin δ.
Then
xp (t) = A cos ωt + B sin ωt

= R cos δ cos ωt + R sin δ sin ωt
= R cos(ωt − δ)
where in the last step we used the cosine addition formula. Observe that R is the amplitude
of oscillation. The quantity δ is called the phase angle. It measures how much the oscillation
lags (if δ > 0) the forcing term. (For example, at t = 0 the amplitude of the forcing term is
a maximum, but the maximum oscillation is delayed until time t = δ/ω.)
Clearly,
A2 + B 2 = R2 cos2 δ + R2 sin2 δ = R2
and
B
tan δ =
A
Substituting the expressions for A and B into the above equations give
m2 (ω02 − ω 2 ) 2 γ 2 ω 2 2
R2 = F0 + F
∆4 ∆4 0
∆2 2
= F
∆4 0
F02
=
∆2
Thus
F0
R= (3.21)
∆
where we recall ∆ is defined in (3.20). Taking the ratio of A and B we see that
γω
tan δ =
m(ω02 − ω 2 )
3.4.1 Resonance
We now examine the behavior of the amplitude of oscillation, R = R(ω), as a function of

the frequency ω of the driving term.
3.4. FORCED OSCILLATIONS OF THE MASS-SPRING SYSTEM 57
1/Delta
omega
Figure 3.4: 1/∆(ω) as a function of ω.
Low frequencies: When ω → 0, ∆(ω) → mω02 = k. Thus for low frequencies the amplitude
of oscillation approaches F0 /k. This result could have been anticipated since when
ω → 0, the forcing term tends to F0 , a constant. A particular solution in this case is
itself a constant and a quick calculation shows this constant is eqaul to F0 /k.
High frequencies: When ω → ∞, ∆(ω) ∼ mω 2 and hence the amplitude of oscillation
R → 0. Intuitively, if you shake the mass-spring system too quickly, it does not have
time to respond before being subjected to a force in the opposite direction; thus, the
overall effect is no motion. Observe that greater the mass (inertia) the smaller R is
for large frequencies.
Maximum Oscillation: The amplitude R is a maximum (as a function of ω) when ∆(ω)
is a minimum. ∆ is a minimum when ∆2 is a minimum. Thus to find the frequency
corresponding to maximum amplitude of oscillation we must minimize
2
m2 ω 2 − ω02 + γ 2 ω 2 .
To find the minimum we take the derivative of this expression with respect to ω and
set it equal to zero:
2m2 (ω 2 − ω02 )(2ω) + 2γ 2 ω = 0.
Factoring the left hand side gives

ω γ 2 + 2m2 (ω 2 − ω02 ) = 0.
Since we are assuming ω = 0, the only way this equation can equal zero is for the
expression in the square brackets to equal zero. Setting this to zero and solving for ω 2
gives the frequency at which the amplitude is a maximum. We call this ωmax :

2 γ2 γ2
ωmax = ω02 − = ω 2
0 1 − .
2m2 2km
Taking the square root gives

γ2
ωmax = ω0 1− .
2km
Assuming γ 1 (the case of very small friction), we can expand the square root to
get the approximate result

γ2
ωmax = ω0 1 − + O(γ 4 ) .
4km
That is, when ω is very close to the natural frequency ω0 we will have maximum
oscillation. This phenomenon is called resonance. A graph of 1/∆ as a function of ω
is shown in Fig. 3.4.
3.5 Exercises
#1. Euler’s formula
Using Euler’s formula prove the trig identity
cos(4θ) = cos4 θ − 6 cos2 θ sin2 θ + sin4 θ.
Again using Euler’s formula find a formula for cos(2nθ) where n = 1, 2, . . .. In this way one
can also get identities for cos(2n + 1)θ as well as sin nθ.
#2. Roots of unity
Show that the n (distinct) solutions to the polynomial equation

xn − 1 = 0
are e2πik/n for k = 1, 2, . . . , n. For n = 6 draw a picture illustrating where these roots lie in
the complex plane.
#3. Constant coefficient ODEs
In each case find the unique solution y = y(x) that satisfies the ODE with stated initial
conditions:
1. y − 3y + 2y = 0, y(0) = 1, y (0) = 0.
2. y + 9y = 0, y(0) = 1, y (0) = −1.
3. y − 4y + 4y = 0, y(0) = 2, y (0) = 0.
3.5. EXERCISES 59
#4. Higher order equations
The third order homogeneous differential equation with constant coefficients is

a3 y + a2 y + a1 y + a0 y = 0 (3.22)
where ai are constants. Assume a solution of the form
y(x) = eλx
and derive an equation that λ must satisfy in order that y is a solution. (You should get a
cubic polynomial.) What is the form of the general solution to (3.22)?
#5. Euler’s equation
A differential equation of the form

t2 y + aty + by = 0, t > 0 (3.23)
where a, b are real constants, is called Euler’s equation.8 This equation can be transformed
into an equation with constant coefficients by letting x = ln t. Solve
t2 y + 4ty + 2y = 0 (3.24)
#6 Forced undamped system
Consider a forced undamped system described by

ÿ + y = 3 cos(ωt)
with initial conditions y(0) = 1 and ẏ(0) = 1. Find the solution for ω = 1.
#7. Driven damped oscillator
Let
ÿ + 3ẏ + 2y = 0
be the equation of a damped oscillator. If a forcing term is F (t) = 10 cos t and the oscillator
is initially at rest at the origin, what is the solution of the equation for this driven damped
oscillator? What is the phase angle?
#8. Damped oscillator
A particle is moving according to

ÿ + 10ẏ + 16y = 0
with the initial condition y(0) = 1 and ẏ(0) = 4. Is this oscillatory ? What is the maximum
value of y?
8 There is perhaps no other mathematician whose name is associated to so many functions, identities,
equations, numbers, . . . as Euler.

#9. Wronskian
Consider (3.3) with p(x) and q(x) continuous on the interval [a, b]. Prove that if two solutions
y1 and y2 have a maximum or minimum at the same point in [a, b], they cannot form a basis
of V.
#10. Euler’s equation (revisited) from physics
In Exercise 2.3.9 we obtained a set of three first-order differential equations for Ω1 , Ω2 and
Ω3 , which are called the Euler equations when there is no torque. Let us assume that
I1 = I2 = I3 . (The body with these moments of inertia is called a free symmetric top.) In
this case we have
I1 Ω̇1 = (I2 − I3 )Ω2 Ω3 (3.25)

I2 Ω̇2 = (I3 − I1 )Ω3 Ω1 (3.26)
I3 Ω̇3 = 0 (3.27)
Notice that Ω3 is a constant from (3.27). Show that Ω1 and Ω2 have the form of
Ω1 (t) = A sin(ωt + θ0 );
Ω2 (t) = A cos(ωt + θ0 )
where A and θ0 are some constants. Here Ω1 , Ω2 and Ω3 are three components of the angular
Show that it follows that the magnitude (length) of Ω
velocity vector Ω. is a constant. Find
an explicit expression for ω in terms of Ii and the constant Ω3 .
Chapter 4
Difference Equations
Figure 4.1: Leonardo Fibonacci, c. 1170–c. 1250.
Science is what we understand well enough to explain to a computer. Art is

everything else we do.
D.E. Knuth in the preface of A=B by H. Wilf & D. Zeilberger
61
62 CHAPTER 4. DIFFERENCE EQUATIONS
4.1 Introduction
We have learned that the general inhomogeneous second order linear differential equation is
of the form
d2 y dy
a(x) 2 + b(x) + c(x)y = f (x).
dx dx
The independent variable, x, takes values in R. (We say x is a continuous variable.) Many
applications lead to problems where the independent variable is discrete; that is, it takes
values in the integers. Instead of y(x) we now have yn , n an integer. The discrete version
of the above equation, called an inhomogeneous second order linear difference equation, is
an yn+2 + bn yn+1 + cn yn = fn (4.1)
where we assume the sequences {an }, {bn }, {cn } and {fn } are known. For example,
3
(n2 + 5)yn+2 + 2yn+1 + yn = en , n = 0, 1, 2, 3, . . .
n+1
is such a difference equation. Usually we are given y0 and y1 (the initial values), and the
problem is to solve the difference equation for yn .
In this chapter we consider the special case of constant coefficient difference equations:
a yn+2 + b yn+1 + c yn = fn
where a, b, and c are constants independent of n. If fn = 0 we say the difference equation

is homogeneous. An example of a homogeneous second order constant coefficient difference
equation is
1
6yn+2 + yn+1 + 2yn = 0.
3
4.2 Constant coefficient difference equations
4.2.1 Solution of constant coefficient difference equations
In this section we give an algorithm to solve all second order homogeneous constant coeffi-
cient difference equations
a yn+2 + b yn+1 + c yn = 0. (4.2)
The method is the discrete version of the method we used to solve contant coefficient differ-
ential equations. We first guess a solution of the form
yn = λn , λ = 0.
(For differential equations we guessed y(x) = eλx .) We now substitute this into (4.2) and
require the result equal zero,
0 = aλn+2 + bλn+1 + cλn

= λn aλ2 + bλ + c .
4.2. CONSTANT COEFFICIENT DIFFERENCE EQUATIONS 63
This last equation is satisfied if and only if
aλ2 + bλ + c = 0. (4.3)
Let λ1 and λ2 denote the roots of this quadratic equation. (For the moment we consider
only the case when the roots are distinct.) Then
λn1 and λn2
are both solutions to (4.2). Just as in our study of second order ODEs, the linear combination
c1 λn1 + c2 λn2
is also a solution and every solution of (4.2) is of this form. The constants c1 and c2 are
determined once we are given the initial values y0 and y1 :
y0 = c1 + c2 ,
y1 = c1 λ1 + c2 λ2 ,
are two equation that can be solved for c1 and c2 .
4.2.2 Fibonacci numbers
Consider the sequence of numbers
0 1 1 2 3 5 8 13 21 34 55 89 144 233 · · ·
that is, each number is the sum of the preceding two numbers starting with
0 1
as initial values. These integers are called Fibonacci numbers and the nth Fibonacci number
is denoted by Fn . The numbers grow very fast, for example,
F100 = 354 224 848 179 261 915 075.
From their definition, Fn satisfies the difference equation
Fn+1 = Fn + Fn−1 for n ≥ 1
with
F0 = 0, F1 = 1.
The quadratic equation we must solve is
λ2 = λ + 1,
whose roots are √

1± 5
λ1,2 = .
2
Setting
Fn = c1 λn1 + c2 λn2 ,
we see that at n = 0 and 1 we require
0 = c1 + c2 ,
1 = c1 λ1 + c2 λ2 .
Solving these we find

1 1
c1 = √ , c2 = − √ ,
5 5
and hence √ n √ n
1 1+ 5 1 1− 5
Fn = √ −√ .
5 2 5 2
Since λ1 > 1 and |λ2 | < 1, λn1 grows with increasing n whereas λn2 → 0 as n → ∞. Thus for
large n
1
Fn ∼ √ λn1 ,
5
and
Fn−1 1
lim = := ω.
n→∞ Fn λ1
The number √
5−1
ω= = 0.61803398 . . . . . .
2
is called the golden mean.1
4.3 Inhomogeneous difference equations
In a completely analogous way to the ODE case, one proves that every solution to the
inhomogeneous linear difference equation (4.1) is of the form
(yn )homo + (yn )part
where (yn )homo is a solution to the homogeneous equation (4.1) with fn = 0 and (yn )part is
a particular solution to the inhomogeneous equation (4.1).
1 More often the number √
1+ 5
φ = 1/ω = = 1.6180339887 . . .
2
is called the golden mean or golden ratio. Two quantities a and b are said to be in the golden ratio φ if
a+b a
= = φ.
a b
In words, a + b is to a as a is to b. A golden rectangle is a rectangle whose ratio of the longer side to the
shorter side is the golden mean φ. Since the ratio of Fibonacci numbers Fn /Fn−1 converges to φ as n → ∞,
rectangles whose long side is of length Fn and whose short side is of length Fn−1 are approximate golden
rectangles. For example, here are some approximate golden rectangles: 1 × 1, 1 × 2, 2 × 3, 3 × 5, 5 × 8,
8 × 13, and so on. The higher we go in this sequence the closer we come to the golden rectangle.
4.4. EXERCISES 65
4.4 Exercises
#1. Degenerate roots
Consider the constant coefficient difference equation (4.2) but now assume the two roots
λ1,2 are equal. Show that
nλn1
is a second linearly independent solution to (4.2).
√
#2. Rational approximations to 2
Solve the difference equation

xn+1 = 2xn + xn−1 , n ≥ 1
with initial conditions x0 = 0 and x1 = 1 that corresponds to the sequence 0, 1, 2, 5, 12,
29,. . . . Show that
xn+1 − xn √
lim = 2.
n→∞ xn
The rational numbers
xn+1 − xn
xn
provide us with very good rational approximations to the square root of two.2
#3. Catalan numbers
Many times nonlinear recurrence relations arise. For example, Catalan numbers Tn satisfy
the nonlinear recurrence relation

n−1
Tn = Tk Tn−1−k , n = 1, 2, . . .
k=0
2 The
√
square
√ root of two is irrational. Here is a proof that 2 is irrational. Suppose not; that is, we
suppose that 2 is a rational number. All rational
√ numbers can be written as a ratio of two integers. Thus
we assume there are integers p and q such that 2 = p/q. By canceling out common factors in p and q we
can assume that the fraction is reduced. By definition of the square root,
p2
2=
q2
which implies p2 = 2q 2 . Thus p2 is an even integer since it is two times q 2 . If p2 is even then it follows that
p is even. (The square of an odd integer is an odd integer.) Since p is even it can be written as p = 2n for
some integer n. Substituting this into p2 = 2q 2 and canceling the common factor of two, we obtain
q 2 = 2p2
But this last equation means q2 ;
and hence q, is an even integer. Thus both p and q have a common factor
of two. But we√assumed that all common factors were cancelled. Thus we arrive at a contradiction to the
assertion that 2 is rational.
From Wikipedia: Pythagoreans discovered that the diagonal of a square is incommensurable with its side,
or in modern language, that the square root of two is irrational. Little is known with certainty about the
time or circumstances of this discovery, but the name of Hippasus of Metapontum is often mentioned. For a
while, the Pythagoreans treated as an official secret the discovery that the square root of two is irrational,
and, according to legend, Hippasus was murdered for divulging it. The square root of two is occasionally
called “Pythagoras’ number” or “Pythagoras’ Constant.”
where T0 := 1. Define
∞

T (z) = Tn z n .
n=0
Show that √
1− 1 − 4z
T (z) = .
2z
From this prove that
1 2n
Tn =
n+1 n

where nk is the binomial coefficient. Catalan numbers arise in a variety of combinatorial
problems. Here is one example:
Suppose 2n points are placed in fixed positions, evenly distributed on the cir-
cumference of a circle. Then there are Tn ways to join n pairs of the points so
that the resulting chords do not intersect.
Figure 4.2: The Catalan numbers T1 = 1, T2 = 2, T3 = 5 and T4 = 14 are illustrated in the

above counting problem of ways to join n pairs so that the chords do not intersect. Here
n = 1, 2, 3, 4. A more pictorial presentation is suppose there are 8 people at a table, then
ask how many ways can they shake hands with one other person and so that no arms are
crossed? The answer is T4 = 14.
One can easily make a table of values of Tn using, say, the Mathematica command (this
gives T1 through T10 ).
Table[{n, Binomial[2*n, n]/(n + 1)}, {n, 1, 10}]
Chapter 5
Matrix Differential Equations
Figure 5.1: Vladimir Arnold, 1937–2010.
Linear systems are almost the only large class of differential equations for which
there exists a definitive theory. This theory is essentially a branch of linear
algebra, and allows us to solve all autonomous linear equations.
V. A. Arnold, Ordinary Differential Equations [2]
67
68 CHAPTER 5. MATRIX DIFFERENTIAL EQUATIONS
5.1 The matrix exponential
Let A be a n×n matrix with constant entries. In this chapter we study the matrix differential
equation
dx
= Ax where x ∈ Rn . (5.1)
dt
We will present an algorithm that reduces solving (5.1) to problems in linear algebra.
The exponential of the matrix tA, t ∈ R, is defined by the infinite series1
t2 2 t3 3
etA = exp(tA) := I + tA + A + A + ··· . (5.2)
2! 3!
Remark: In an advanced course you will prove that this infinite series of matrices converges
to a n × n matrix.
It is important to note that for matrices A and B, in general,
exp(tA) exp(tB) = exp(tA + tB).
If A and B commute (AB = BA) then it is the case that
exp(tA) exp(tB) = exp(tA + tB).
This last fact can be proved by examining the series expansion of both sides—on the left
hand side one has to multiply two infinite series. You will find that by making use of
AB = BA the result follows precisely as in the case of complex exponentials.
Here are some examples:
1.
A = D = diagonal matrix = diag (λ1 , λ2 , . . . , λn ) .
Observe that

Dk = diag λk1 , λk2 , . . . , λkn .
Thus
∞ k
t
Dk = diag etλ1 , etλ2 , . . . , etλn .
k!
k=0
1 We put the scalar factor t directly into the definition of the matrix exponential since it is in this form
we will use the matrix exponential.

5.1. THE MATRIX EXPONENTIAL 69
2. Suppose that A is a diagonalizable matrix; that is, there exist matrices U and D with
U invertible and D diagonal such that
A = U DU −1 .
Observe
A2 = (U DU −1 )(U DU −1 ) = U D2 U −1 ,
and more generally,
Ak = U Dk U −1 .
Thus
∞ k
t
exp(tA) = Ak
k!
k=0
∞ k
t
= U Dk U −1
k!
k=0
∞ k
t k −1
= U D U
k!
k=0
= U exp(tD)U −1 . (5.3)
In the next to the last equality, we used the fact that U and U −1 do not depend
upon the summation index k and can therefore be brought outside of the sum. The
last equality makes use of the previous example where we computed the exponential
of a diagonal matrix. This example shows that if one can find such U and D, then
the computation of the exp(tA) is reduced to matrix multiplications. This last result,
(5.3), is quite suitable for using MatLab or Mathematica.
3. Let

0 −1
A= .
1 0
Matrix multiplication shows
A2 = −I,
and thus
k
A2k = A2 = (−I) = (−1)k I,
k
A2k+1 = A2k A = (−1)k A.

Hence
∞ k
t
exp(tA) = Ak (5.4)
k!
k=0
∞ ∞
t2k 2k t2k+1
= A + A2k+1
(2k)! (2k + 1)!
k=0 k=0
∞ 2k ∞
t2k+1
t
= (−1)k I + (−1)k A
(2k)! (2k + 1)!
k=0 k=0
= cos t I + sin t A

cos t 0 0 − sin t
= +
0 cos t sin t 0

cos t − sin t
= . (5.5)
sin t cos t
Remark: You can also compute

0 −1
exp t
1 0
by the method of Example #2. Try it!
5.2 Application of etA to differential equations

5.2.1 Derivative of etA with respect to t
The following is the basic property of the exponential that we apply to differential equations.
As before, A denotes a n × n matrix with constant coefficients.
d
exp(tA) = A exp(tA) = exp(tA)A. (5.6)
dt
Here is the proof: Differentiate

t 2 2 t3 3
etA = I + tA + A + A + ···
2! 3!
term-by-term2 with the result
d tA t2
e = A + tA2 + A3 + · · ·
dt 2!
t2
= A I + tA + A2 + · · ·
2!
= A etA
= etA A.
2 In a complex analysis course you will prove that convergent complex power series can be differentiated
term-by-term and the resulting series has the same radius of convergence. Note there really is something to
prove here since there is an interchange of two limits.
5.2. APPLICATION OF MATRIX EXPONENTIAL TO DES 71
The last equality follows by factoring A out on the right instead of the left.
5.2.2 Solution to matrix ODE with constant coefficients
We now use (5.6) to prove

Theorem: Let
dx
= Ax (5.7)
dt
where x ∈ Rn and A is a n × n matrix with constant coefficients. Then every solution of
(5.7) is of the form
x(t) = exp(tA)x0 (5.8)
for some constant vector x0 ∈ Rn .
Proof : (i) First we show that x(t) = etA x0 is a solution:

dx d tA d tA
= e x0 = e x0
dt dt dt
= AetA x0
= Ax(t).
(ii) We now show that every solution of (5.7) is of the form (5.8). Let y(t) be any solution
to (5.7). Let
∆(t) := e−tA y(t).
If we can show that ∆(t) is independent of t—that it is a constant vector which we call x0 ,
then we are done since multiplying both sides by etA shows
etA x0 = etA ∆(t) = etA e−tA y(t) = y(t).
(We used the fact that tA and −tA commute so that the addition formula for the matrix
exponential is valid.) To show that ∆(t) is independent of t we show its derivative with
respect to t is zero:
d∆ d −tA
= e y(t)
dt dt

d −tA dy
= e y(t) + e−tA (product rule)
dt dt
−tA
= −e A y(t) + e−tA (Ay(t)) (y(t) satisfies ODE)
= 0.
The next theorem relates the solution x(t) of (5.7) to the eigenvalues and eigenvectors of
the matrix A (in the case A is diagonalizable).
Theorem: Let A be a diagonalizable matrix. Any solution to (5.7) can be written as
x(t) = c1 etλ1 ψ1 + c2 etλ2 ψ2 + · · · + cn etλn ψn (5.9)

where λ1 , . . . , λn are the eigenvalues of A with associated eigenvectors ψ1 , . . . , ψn , and

c1 , . . . , cn are constants.
Proof : All solutions of (5.7) are of the form (5.8). Since A is diagonalizable, the eigen-
vectors of A can be used to form a basis: {ψ1 , . . . , ψn }. Since this is a basis there exist
constants c1 , . . . , cn such that
x0 = c1 ψ1 + c2 ψ2 + · · · + cn ψn .
(x0 is the constant vector appearing in (5.8).)

For any eigenvector ψ of A with eigenvalue λ we have
etA ψ = etλ ψ.
(This can be proved by applying the infinite series (5.2) to the eigenvector ψ and noting
Ak ψ = λk ψ for all positive integers k.) Thus
etA x0 = c1 etA ψ1 + · · · cn etA ψn

= c1 etλ1 ψ1 + · · · + cn etλn ψn .
Here are two immediate corollaries of this theorem:
1. If A is diagonalizable and has only real eigenvalues, then any solution x(t) of (5.1) will
have no oscillations.
2. If A is diagonalizable and the real part of every eigenvalue is negative, then
x(t) → 0 (zero vector), as t → +∞
To see this recall that if λ = σ + iµ (σ and µ both real), then
eλt = eσt eiµt .
If σ < 0, eσt → 0 as t → +∞. Now apply preceding theorem.
Example: Here is an example using MatLab or Mathematica to solve a system of equa-

tions. Consider  
1 −1 2
dx
= Ax where A =  2 2 0 .
dt
3 1 2
The solution x(t) to this DE that satisfies the initial condition x(0) = x0 is
x(t) = etA x0
so we must compute the 3 × 3 matrix etA . Here is a step-by-step approach:
Step 1: Find the eigenvalues of A

1−λ −1 2

p(λ) = det (A − λI) = 2 2−λ 0 = −λ3 + 5λ2 − 4λ = −λ(λ − 1)(λ − 4)

3 1 2−λ
5.2. APPLICATION OF MATRIX EXPONENTIAL TO DES 73
Thus the eigenvalues of A are
λ1 = 4; λ2 = 1, λ3 = 0.
Since the matrix is 3 × 3 and there are three distinct eigenvalues, we know that A is
diagonalizable.
Step 2: Find the eigenvectors of A. That is we must solve the linear equations
Aψ1 = λ1 ψ1 , Aψ2 = λ2 ψ2 , Aψ3 = λ3 ψ3
where eigenvalues are given in Step 1. Note that the solutions will be determined only
up to an overall constant, e.g. if Aψ = λψ, then cψ is also an eigenvector where c is a
constant. Solving the above equations one finds
     
1 −1 −1
ψ1 =  1  , ψ2 =  2  , ψ3 =  1  .
2 1 1
Step 3: Form matrices U and D and compute U −1 . D is the diagonal matrix

   
λ1 0 0 4 0 0
D =  0 λ2 0 = 0 1 0 
0 0 λ3 0 0 0
and U is the matrix formed from the eigenvectors:

   
↑ ↑ ↑ 1 −1 −1
U =  ψ1 ψ2 ψ3  =  1 2 1 .
↓ ↓ ↓ 2 1 1
We now compute the inverse of U :

 
1/3 0 1/3
U −1 =  1/3 1 −2/3 
−1 −1 1
Thus we have
A = U DU −1
Step 4: Compute etA from formula

   4t   
1 −1 −1 e 0 0 1/3 0 1/3
etA = U etD U −1 =  1 2 1  ·  0 et 0  ·  1/3 1 −2/3 
2 1 1 0 0 1 −1 −1 1
 t 4t t 4t

1 − e3 + e3 1 − et −1 + 2e3 + e3
 e4t 4t 
=  −1 + 3 + 3 −1 + 2et 1 − 4e3 + e3 
2et t
 
et 2e4t 2et 2e4t
−1 + 3 + 3 −1 + e t
1− 3 + 3
Step 5: If the initial condition is vector

 
c1
x0 =  c2 
c3
then the solution to the differential equation dx/dt = Ax with initial condition x(0) =
x0 is
 t 4t t 4t
  
1 − e3 + e3 1 − et −1 + 2e3 + e3 c1
 
x(t) = etA x0 =  −1 + 2e3 + e3 −1 + 2et 1 − 4e3 + e3  ·  c2 
t 4t t 4t
t 4t t 4t
−1 + e3 + 2e3 −1 + et 1 − 2e3 + 2e3 c3
 t 4t
t 4t

1 − e3 + e3 c1 + (1 − et ) c2 + −1 + 2e3 + e3 c3
 
 t 4t t 4t 
=  −1 + 2e3 + e3 c1 + (−1 + 2et ) c2 + 1 − 4e3 + e3 c3 
 t 4t
t 4t

−1 + e3 + 2e3 c1 + (−1 + et ) c2 + 1 − 2e3 + 2e3 c3
Step 6: An alternative to Step 5 is to simply say the general solution is of the form
x(t) = c1 e4t ψ1 + c2 et ψ2 + c3 ψ3
where c1 , c2 and c3 are arbitrary constants (not the constants cj appearing in Step 5!)
and ψj are the eigenvectors found above.
5.3 Relation to earlier methods of solving constant co-

efficient DEs
Earlier we showed how to solve
aÿ + bẏ + cy = 0
where a, b and c are constants. Indeed, we proved that the general solution is of the form
y(t) = c1 etλ1 + c2 etλ2
where λ1 and λ2 are the roots to
aλ2 + bλ + c = 0.
(We consider here only the case of distinct roots.)
Let’s analyze this familiar result using matrix methods. The x ∈ R2 is

x1 y
x(t) = =
x2 dy/dt
Therefore,

dx dy/dt
=
dt d2 y/dt2

x2
=
− b x2 − ac x1
a
0 1 x1
= .
− ac − ab x2
5.4. PROBLEM FROM MARKOV PROCESSES 75
This last equality defines the 2 × 2 matrix A. The characteristic polynomial of A is

−λ
p(λ) = det (A − λI) = c
1 = λ2 + b λ + c .
−a −a − λ
b
a a
Thus the eigenvalues of A are the same quantities λ1 and λ2 appearing above. Since

etλ1 0
x(t) = etA x0 = S S −1 x0 ,
0 etλ2
x1 (t) is a linear combination of eλ1 t and eλ2 t .
5.4 Problem from Markov processes
In the theory of continuous time Markov processes (see, e.g., Chapter 2 in [7]), one has a set
of states and one asks for the transition probability from state i to state j in time t. Here is
an example3
Suppose we have three states that we label 1, 2 and 3. We are given that the rates of
transition from i −→ j are known numbers qij . In this example suppose
q12 = 1, q13 = 1, q21 = 1, q23 = 0, q31 = 2, q32 = 1.
The theory of Markov processes tells us that we form a Q-matrix whose off-diagonal elements
are qij , i = j, and whose row sums equal zero. In this example
 
−2 1 1
Q =  1 −1 0  (5.10)
2 1 −3
Denote by pij (t) the transition probability from state i to state j in time t. Form the matrix
P (t) = (pij (t)) ,
then the theory of Markov processes tells us that P (t) satisfies the matrix DE
dP
= QP (5.11)
dt
where P (0) = I. We know that the solution to this DE is
P (t) = exp(tQ). (5.12)
So we must compute exp(tQ) where Q is given by (5.10).

3 This example is taken from [7], Chapter 2.
5.4.1 Computing exp(tQ)
Step 1. We first compute the eigenvalues of Q. To do this we compute the characteristic

polynomial

−2 − λ 1 1

p(λ) = det(Q−λI) = 1 −1 − λ 0 = −8λ−6λ2 −λ3 = −λ(λ+2)(λ+4)

2 1 −3 − λ
Thus the eigenvalues are 0, −2 and −4. Since there are three distinct eigenvalues and
the matrix is 3 × 3, we know that Q is a diagonalizable matrix, see (5.3).4
Step 2. We next compute the eigenvectors corresponding to the eigenvalues. That is, we
must solve the linear equations
QΨk = λk Ψk , k = 1, 2, 3
where λ1 = 0, λ2 = −2, λ3 = −4. Linear algebra computations show5

     
1 1 −3
Ψ1 =  1  , Ψ2 =  −1  , Ψ3 =  1  . (5.13)
1 1 5
Step 3. Form matrices U and U −1 . We know that U is formed from the eigenvectors Ψk :
   
↑ ↑ ↑ 1 1 −3
U =  Ψ1 Ψ2 Ψ3  =  1 −1 1 
↓ ↓ ↓ 1 1 5
We now compute U −1 . A linear algebra computation gives6

 3 1 1 
8 2 8
 
U −1 =  1
4 − 21 1
4 
− 18 0 1
8
Step 4. We now have enough information to compute exp(tQ). We let D denote the
diagonal matrix
   
λ1 0 0 0 0 0
D =  0 λ2 0  =  0 −2 0 
0 0 λ3 0 0 −4
4 If one first defines the matrix Q, then the command Det[Q-λ*IdentityMatrix[3]] in Mathematica finds
the characteristic polynomial.

5 Using a computer software package makes this step easier. For example, the command Eigenvectors[Q]
in Mathematica gives the eigenvectors of the matrix Q.

6 After defining the matrix U , the command Inverse[U] in Mathematica computes the inverse matrix
U −1 .
5.5. APPLICATION OF MATRIX DE TO RADIOACTIVE DECAYS 77
From (5.3) we know that
P (t) = exp(tQ) = U exp(tD)U −1

   3 1 1 
1 1 −3 1 0 0 8 2 8
 
=  1 −1 1   0 e−2t 0   41 − 21 1
4 
1 1 5 0 0 e−4t − 18 0 1
8
 
3 3e−4t e−2t 1 e−2t 1 3e−4t e−2t
8 + 8 + 4 2 − 2 8 − 8 + 4
 
 −4t −2t
e−2t e−4t e−2t 
=  83 − e 8 − e 4 1
+ 1
+ −  (5.14)
 2 2 8 8 4 
3 5e−4t e−2t 1 e−2t 1 5e−4t e−2t
8 − 8 + 4 2 − 2 8 + 8 + 4
Thus, for example, the transition probability from state 1 to 1 in time t is given by the
(1, 1) matrix element of P (t):
3 1 −2t 3 −4t
p11 (t) = + e + e .
8 4 8
and the transition probability from state 3 to 2 in time t is given by the (3, 2) matrix element
of P (t):
1 1
p32 (t) = − e−2t .
2 2
5.4.2 Some observations

1. Observe that the row sums of P (t) are all equal to one. For example in the first row

3 3 −4t 1 −2t 1 1 1 3 1
p11 (t)+p12 (t)+p13 (t) = + e + e + − e−2t + − e−4t + e−2t = 1
8 8 4 2 2 8 8 4
This is simply the statement that the probability of starting in an initial state i and
being in any state j is one.
2. Observe that  
3 1 1
8 2 8
 3 1 1 
lim P (t) =  8 2 8 
t→∞
3 1 1
8 2 8
This means that after a long time no matter what is the initial state i, the transition
probability to any state j is independent of the initial state.
5.5 Application of matrix DE to radioactive decays

Consider a radioactive decay of type
A −→ B −→ C
where C is assumed to be stable. We denote the decay rate of A by λA and the decay
rate of B by λB . (Since C is assumed to be stable, the decay rate of C equals zero.) An
example of such a radioactive decay (by β-decay) is Iodine–133 decays to Xenon–133 which
then decays to Cesium–133.
Introduce pα (t) equal to the amount of material α = A, B, C present at time t. We start
with an initial amount of material A which for simplicity we take to be 1. Thus the initial
conditions are
pA (0) = 1, pB (0) = 0, pC (0) = 0.
The differential equations satisfied by the pα , α = A, B, C, are
dpA
= −λA pA ,
dt
dpB
= λA pA − λB pB ,
dt
dpC
= λB pB .
dt
We can write this in matrix form
dP
= AP
dt
where  
pA
P =  pB 
pC
and  
−λA 0 0
A =  λA −λB 0 .
0 λB 0
We know that
P (t) = exp(tA)P (0)
so we need to compute exp(tA).
Step 1. Since A is a lower triangular matrix, the eigenvalues of A are −λA , −λB and 0.
Step 2. The corresponding eigenvectors of A are
 λ −λ 
A B    

λB
 0 0
ΨA =  − λA
λB
 , ΨB =  −1  , ΨC =  0  .

1 1
1
(Mathematica does these calculations easily.)
Step 3. The columns of U are the vectors ΨA , ΨB , ΨC .
Step 4.
exp(tA) = U exp(tD)U −1
   −λ t  −1
(λA − λB )/λB 0 0 e A 0 0 (λA − λB )/λB 0 0
=  −λA /λB −1 0  0 e−λB t 0   −λA /λB −1 0 
1 1 1 0 0 1 1 1 1
 −λA t

e 0 0
=  −λB t
− e−λ At
e−λB t 0 
−1
λA (e −λB t
)
−λA t

(λA − λB ) λA (1 − e − λB 1 − e ) 1 − e−λB t 1
5.6. INHOMOGENOUS MATRIX EQUATIONS 79
Step 5. To find the solution P (t) we must apply the matrix exp(tA) to the initial vector
 
1
P (0) =  0  .
0
This results in
pA (t) = e−λA t ,
λA −λA t
pB (t) = e − e−λB t ,
λB − λA
1
pC (t) = 1− λA e−λB t − λB e−λA t .
λA − λB
It is an exercise in calculus to show that the maximum of pB (t) occurs at

log(λA /λB )
tmax = .
λA − λB
We now apply this to the example of Iodine–133 decay. The half-life of A:=Iodine-133
is 20.81 hours and the half-life of B:=Xenon–133 is 5.243 days. From this we can calculate
the rates after we decide whether our unit of time is hours or days—let’s take days. Then7
λA = 0.7994, λB = 0.1322.
Thus, for example, the maximum amount of Xenon–133 occurs in approximately 2.70 days.
In the figure we show the populations of the three nuclides as a function of time.
5.6 Inhomogenous matrix equations

Consider the inhomogenous equation
dx
= Ax + f (t), x(0) = x0 (5.15)
dt
where x is a vector of dimension n, A a n × n constant coefficient matrix, and f (t) is a given
vector which in general depends, say continuously, upon the independent variable t. We use
the method of variation of parameters to find a particular solution to (5.15). Let8
x(t) = etA y(t)
Then
dx dy
= AetA y(t) + etA
dt dt
To satisfy the differential equation this must equal
AetA y(t) + f (t)

7 Recall the formula T1/2 = log
λ
2
.
8 If
f (t) = 0 then we know that y(t) would be a constant vector. For nonzero f (t) we are allowing for the
possibility that y can depend upon t; hence the name variation of parameters.
1.0
0.8
0.6 pA t
pB t
0.4
pC t
0.2
t
5 10 15 20 25
Figure 5.2: The population of A=Iodine–133, B=Xenon–133 and C=Cesium–133 starting

with an initial amount of Iodine–133. The unit of time is in days.
5.6. INHOMOGENOUS MATRIX EQUATIONS 81
and hence we must have

dy
etA = f (t)
dt
Solving this for dy/dt:
dy
= e−tA f (t)
dt
The right hand side of the above equation is expressed in terms of known quantities. Inte-
grating gives
t
y(t) = e−sA f (s) ds
0
and hence the particular solution

t
xpart (t) = etA e−sA f (s) ds
0
Thus the solution satisfying the initial condition is
t
x(t) = etA e−sA f (s) ds + etA x0 (5.16)
0
Observe that the solution of (5.15) has been reduced in (5.16) to matrix calculations and
integration.
5.6.1 Nonautonomous linear equations
Consider the inhomogeneous equation
dx
= A(t)x + f (t) (5.17)
dt
where now we assume the n × n matrix A(t) depends upon the independent variable t and
f (t), as before, is a given column vector of size n which in general depends upon t. We
assume that the coefficients of A(t) are continuous on some closed interval [a, b] and for
simplicity we take a < 0 < b. Everything we say below is for t in the interval a ≤ t ≤ b.
The homogeneous version of (5.17) is
dx
= A(t)x. (5.18)
dt
Since A is now assumed to depend upon t, the solution to (5.18) satisfying the initial
condition x(0) = x0 is not exp(tA)x0 . So what plays the role of exp(tA) when A is a
function of t?
Let xj (t), a column vector of size n, denote the solution to (5.18) satisfying the initial
condition xj (0) = ej , j = 1, 2, . . . , n, where {ek }k=1,2,...n denotes the standard basis of
Rn . It can be shown that the functions {xj (t)}j=1,...,n form a basis for the vector space of
solutions to (5.18). Form the n × n matrix X (t) whose columns are the xj (t):
 
↑ ↑ ··· ↑
X (t) =  x1 (t) x2 (t) · · · xn (t)  .
↓ ↓ ··· ↓
X (t) satisfies
X (0) = In ,
where In is the n × n identity matrix. (This follows immediately from the fact that xj (0) =
ej .) Furthermore, it can be shown that X (t) is invertible for all t. Note that X (t) satisfies
the matrix differential equation
dX (t)
= A(t)X (t),
dt
a fact that follows immediately from the fact that for each j, dxj /dt = A(t)xj . The n × n
matrix X (t) is called the fundamental matrix. In the case when A has constant coefficients,
the fundamental matrix is exp(tA).
We now use the matrix X (t) to construct a particular solution, xpart (t), to the inhomo-
geneous equation (5.17) by following the ideas for the constant coefficient case. We seek
solutions to (5.17) of the form
xpart (t) = X (t)y(t).
We differentiate xpart (t) with respect to t:
dxpart d
= (X (t)y(t))
dt dt
dX dy
= y(t) + X (t)
dt dt
dy
= A(t)X (t)y(t) + X (t) .
dt
We want this to equal A(t)xpart (t) + f (t) so we get the equation
dy
A(t)X (t) y(t) + X (t) = A(t)xpart (t) + f (t) = A(t)X (t)y + f (t)
dt
or
dy
X (t)
= f (t).
dt
Multiplying both sides of the above equation by [X (t)]−1 gives
dy
= [X (t)]−1 f (t).
dt
Now integrate the above from 0 to t
t
y(t) = [X (s)]−1 f (s) ds,
0
and recalling xpart (t) = X (t)y gives

t
xpart (t) = X (t)y(t) = X (t) [X (s)]−1 f (s) ds.
0
5.7. EXERCISES 83
Since the general solution to (5.17) is a particular solution plus a solution to the homogeneous
equation, we have
t
x(t) = X (t) [X (s)]−1 f (s) ds + X (t)x0 (5.19)
0
solves (5.17) with initial condition x(0) = x0 . This formula should be compared with (5.16).
To summarize,
To solve the inhomogeneous equation (5.17), one first finds the basis of solutions
{xj (t)}j=1,...,n , xj (0) = ej , to the homogeneous equation (5.18) and constructs
the n × n fundamental matrix X (t). Then the solution to the inhomogeneous
equation satisfying the initial condition x(0) = x0 is given by (5.19). Thus the
real difficulty is in solving the homogeneous problem, i.e. finding the fundamental
matrix X (t).
5.7 Exercises
#1. Harmonic oscillator via matrix exponentials
Write the oscillator equation

ẍ + ω02 x = 0
as a first order system (5.1). (Explicitly find the matrix A.) Compute exp(tA) and show
that x(t) = exp(tA)x0 gives the now familiar solution. Note that we computed exp(tA) in
(5.5) for the case ω0 = 1.
#2. Exponential of nilpotent matrices
1. Using the series expansion for the matrix exponential, compute exp(tN ) where

0 1
N= .
0 0
Answer the same question for
 
0 1 1
N = 0 0 1 .
0 0 0
How do these answers differ from exp(tx) where x is any real number?
2. A n × n matrix N is called nilpotent9 if there exists a positive integer k such that
Nk = 0
where the 0 is the n × n zero matrix. If N is nilpotent let k be the smallest integer
such that N k = 0. Explain why exp(tN ) is a matrix whose entries are polynomials in
t of degree at most k − 1.
9 In an advanced course in linear algebra, it will be proved that every matrix A can be written uniquely
#3. Computing etA
Let  
1 −1 4
A= 3 2 −1  (5.20)
2 1 −1
1. Find the eigenvalues and eigenvectors of A. (You can use any software package you
like and merely quote the results.)
2. Use these to compute etA .
#4.
Consider the system of linear differential equations

dx
= Ax
dt
where A is the 4 × 4 matrix
 
−5/2 1 1/2 −1/2
 3/4 −5/2 0 3/4 
A=

 (5.21)
1 2 −3 1 
0 2 −1/2 −2
Prove that all solutions x(t) to this DE tend to zero as t → ∞. Hint: You need not compute
etA . You can prove this statement simply by computing the eigenvalues of A. (Why?)
#5.

dx
= Ax
dt
 
0 0 −3/2 2
 −3/4 1/2 0 −3/4 
A=
 −1 −2
 (5.22)
1 −1 
1/2 −3 3/2 −3/2
Find a subspace V of R4 such that if x(0) ∈ V , then x(t) → 0 as t → ∞. Hint: The subspace
V is described in terms of (some of) the eigenvectors of A.
as D + N where D is a diagonalizable matrix, N is a nilpotent matrix, and DN = N D. Furthermore, an
algorithm will be given to find the matrices D and N from the matrix A. Once this is done then one can
compute exp(tA) as follows
exp(tA) = exp(tD + tN ) = exp(tD) exp(tN ).
We showed above how to reduce the computation of exp(tD), D a diagonalizable matrix, to linear algebra.
This problem shows that exp(tN ) reduces to finitely many matrix multiplications. Thus the computation of
both exp(tD) and exp(tN ) are reduced to linear algebra and hence so is exp(tA). Observe that it is crucial
that we know DN = N D.
5.7. EXERCISES 85
#6.
dx
= Ax
dt
1 α
A= (5.23)
−α 3
For what values of α will the solutions exhibit oscillatory behavior?
#7. Radioactive decay & first introduction to Laplace transforms
Birth processes have been used since the time of Rutherford to model radioactive decay.
(Radioactive decay occurs when an unstable isotope transforms to a more stable isotope,
generally by emitting a subatomic particle.) In many cases a radioactive nuclide A decays
into a nuclide B which is also radioactive; and hence, B decays into a nuclide C, etc.
The nuclides B, C, etc. are called the progeny (formerly called daughters). This continues
until the decay chain reaches a stable nuclide. For example, uranium-238 decays through
α-emission to thorium-234 which in turn decays to protactinium-234 through β-emission.
This chain continues until the stable nuclide lead-206 is reached.
1. Let the decay states be E1 → E2 → · · · → EN where EN is the final stable state. We

can relabel these states to be simply 1, 2, . . . , N . (That is, we write Ej as simply j.)
Let N (t) denote the state of the nuclide at time t. N (t) is a random process (called
a Markov process) due to the fact that radioactive decay is inherently random. Thus
we introduce
pj (t) = P(N (t) = j|N (0) = 1)

= probability that nuclide is in state j at time t
given it starts in state 1 at time t = 0.
These probabilities pj (t) satisfy differential equations called the Kolmogorov forward equations:
dpj
= λj−1 pj−1 (t) − λj pj (t), j = 1, 2, . . . , N. (5.24)
dt
The constants λj are called the decay rates. A decay rate λ is related to the half-life,
T1/2 , of the nuclide by the well-known formula
log 2
T1/2 = , log 2 = 0.693147 . . . (5.25)
λ
We assume λi = λj for i, j = 1, . . . , N − 1. We set λ0 = 0 and λN = 0. (λN is set
equal to zero since the final state N is a stable nuclide and does not decay.)
In applications to radioactive decay, if N1 is the number of initial nuclides (the number
of nuclides in state E1 ), then N1 pj (t) is the number of nuclides in state Ej at time t.
2. Introduce the Laplace transform10

∞
p̂j (s) = e−ts pj (t) dt
0
and show that the Laplace transform of (5.24) is
sp̂j (s) − δj,1 = λj−1 p̂j−1 (s) − λj p̂j (s), j = 1, . . . , N. (5.26)
Solve these equations for p̂j (s) and show that
λ1 λ2 λj−1 1
p̂j (s) = ···
s + λ1 s + λ2 s + λj−1 s + λj
3. Using the above expression for p̂j (s) partial fraction the result:

j
cj,k
p̂j (s) =
s + λk
k=1
See if you can find expressions for cj,k . You might want to take some special cases to
see if you can make a guess for the cj,k . (The Mathematica command Apart will
prove useful.)
4. From the partial fraction decomposition of p̂j (s) explain why you can almost immedi-
ately conclude

j
pj (t) = cj,k e−λk t , j = 1, 2, . . . , N. (5.27)
k=1
5. For the special case of N = 4: E1 → E2 → E3 → E4 find explicitly the proba-

bilities pj (t). (You can use Mathematica if you wish. Note there is a command
InverseLaplaceTransform.)
6. Show that p2 (t) has a maximum at t = tm
log(λ1 /λ2 )
tm = > 0.
λ1 − λ2
In terms of the radioactive decay interpretation, this is the time when the first progeny
has a maximum population.
7. Using Mathematica (recall the command Series) show that as t → 0
p1 (t) = 1 − λ1 t + O(t2 )
p2 (t) = λ1 t + O(t2 )
1
p3 (t) = λ1 λ2 t2 + O(t3 )
2
1
p4 (t) = λ1 λ2 λ3 t3 + O(t4 )
3!
10 See Chapter 8 of these Notes and Boyce & Diprima, Chapter 6 [4].
5.7. EXERCISES 87
8. Radon 222 gas is a chemically inert radioactive gas that is part of the Uranium 238
decay chain. Radon and its radioactive progeny are known carcinogens. Here is part
of the decay chain11
· · · −→ Rn 222 −→ Po 218 −→ Pb 214 −→ Bi 214 −→ · · ·
The half-life of each nuclide is known (recall (5.25)):
Rn 222: T1/2 = 3.8235 days

Po 218: T1/2 = 3.10 minutes
Pb 214: T1/2 = 26.8 minutes
Bi 214: T1/2 = 19.9 minutes
Let NRn denote the initial amount of Rn 220 and assume the other nuclides are not
present at time t = 0. Solve the Kolmogorov forward equations for this particular
birth process. (Note that here the probabilities do not sum to one since the Bi 214
also decays.) This is not so messy if you use Mathematica. Find the times when
each of the progeny have maximum population. (Highest probability) You might want
to use Mathematica’s FindRoot.
11 Po=polonium, Pb=lead, Bi=bismuth.

Figure 5.3: From the EPA website: Radon is a radioactive gas. It comes from the natural
decay of uranium that is found in nearly all soils. It typically moves up through the ground
to the air above and into your home through cracks and other holes in the foundation. Your
home traps radon inside, where it can build up. Any home may have a radon problem.
This means new and old homes, well-sealed and drafty homes, and homes with or without
basements. Radon from soil gas is the main cause of radon problems. Sometimes radon
enters the home through well water. In a small number of homes, the building materials can
give off radon, too. However, building materials rarely cause radon problems by themselves.
Chapter 6
Weighted String
Figure 6.1: Hermann von Helmholtz, 1821–1894.
Because linear equations are easy to solve and study, the theory of linear oscil-
lations is the most highly developed area of mechanics. In many nonlinear prob-
lems, linearization produces a satisfactory approximate solution. Even when this
is not the case, the study of the linear part of a problem is often a first step, to
be followed by the study of the relation between motions in a nonlinear system
and its linear model.
V. I. Arnold, Mathematical Methods of Classical Mechanics [1]
89
90 CHAPTER 6. WEIGHTED STRING
6.1 Derivation of differential equations

The weighted string is a system in which the mass is concentrated in a set of equally spaced
mass points, N in number with spacing d, imagined to be held together by massless springs
of equal tension T . We further assume that the construction is such that the mass points
move only in the vertical direction (y direction) and there is a constraining force to keep
the mass points from moving in the horizontal direction (x direction). We call it a “string”
since these mass points give a discrete string—the tension in the string is represented by the
springs. The figure below illustrates the weighted string for N = 5.
0 1d 2d 3d 4d 5d 6d
The string is “tied down” at the endpoints 0 and (N + 1)d. The horizontal coordinates
of the mass points will be at x = d, 2d, . . . , N d. We let uj denote the vertical displacement
of the j th mass point and Fj the transverse force on the j th particle. To summarize the
variables introduced so far:
m = mass of particle,
N = total number of particles,
T = tension of spring,
d = horizontal distance between two particles,
uj = vertical displacement of j th particle, j = 1, 2, . . . N,
Fj = transverse force on j th particle, j = 1, 2, . . . N.
To impose the boundary conditions that the ends of the string are rigidly fixed at x = 0 and
x = (N + 1)d, we take
u0 = 0 and uN +1 = 0.
Newton’s equations for these mass points are
d2 uj
Fj = m , j = 1, 2, . . . , N.
dt2
This is a system of N second order differential equations. We now find an expression for the
transverse force Fj in terms of the vertical displacements.
In the diagram below, the forces acting on the j th particle are shown.
T
β
α
T
(j − 1)d jd (j + 1)d
6.1. DERIVATION OF DIFFERENTIAL EQUATIONS 91
From the diagram,

Fj = T sin β − T sin α.
We make the assumption that the angles α and β are small. (The string is not stretched
too much!) In this small angle approximation we have
sin α ≈ tan α and sin β ≈ tan β.
Therefore, in this small angle approximation
Fj ≈ T tan β − T tan α

uj+1 − uj uj − uj−1
= T −T .
d d
Thus,
d2 uj T
m = (uj+1 − 2uj + uj−1 ) , j = 1, 2, . . . , N. (6.1)
dt2 d
Note that these equations are valid for j = 1 and j = N when we interpret u0 = 0 and
uN +1 = 0. For example, for j = 1 the force F1 is determined from the diagram:
T
T
0 d 2d
(u2 − u1 ) u1
F1 = T −T
d d
T
= (u2 − 2u1 + u0 ) , u0 = 0.
d
Equation (6.1) is a system of N second order linear differential equations. Thus the
dimension of the vector space of solutions is 2N ; that is, it takes 2N real numbers to specify
the initial conditions (N initial positions and N initial velocities). Define the N × N matrix
 
2 −1 0 0 ··· 0 0 0
 −1 2 −1 0 · · · 0  0 0
 
 0 −1 2 −1 · · · 0  0 0
 
VN = .. .. .. .. ..  .. .. (6.2)
 . . . . ··· .  . .
 
 0 0 0 0 ··· −1 2 −1 
0 0 0 0 ··· 0 −1 2
and the column vector u  

u1
 u2 
 
u=
 · .
 (6.3)
 · 
uN
Then (6.1) can be written in the compact matrix form
d2 u T (6.4)
+ VN u = 0.
dt2 md
Note: We could also have written (6.1) as a first order matrix equation of the form
dx
= Ax (6.5)
dt
where A would be a 2N × 2N matrix. However, for this application it is simpler to develop
a special theory for (6.4) rather than to apply the general theory of (6.5) since the matrix
manipulations with VN will be a bit clearer than they would be with A.
6.2 Reduction to an eigenvalue problem
Equation (6.4) is the matrix version of the harmonic oscillator equation
d2 x k k
+ x = 0, ω02 = . (6.6)
dt2 m m
Indeed, we will show that (6.4) is precisely N harmonic oscillators (6.6)—once one chooses
the correct coordinates. We know that solutions to (6.6) are linear combinations of
cos ω0 t and sin ω0 t.
Thus we “guess” that solutions to (6.4) are linear combinations of the form
cos ωt f and sin ωt f
where ω is to be determined and f is a column vector of length N . (Such a “guess” can be

theoretically deduced from the theory of the matrix exponential when (6.4) is rewritten in
the form (6.5).)
Thus setting
u = eiωt f ,
we see that (6.4) becomes the matrix equation
md 2
VN f = ω f.
T
6.3. COMPUTATION OF THE EIGENVALUES 93
That is, we must find the eigenvalues and eigenvectors of the matrix VN . Since VN is a real
symmetric matrix, it is diagonalizable with real eigenvalues. To each eigenvalue λn , i.e.
VN fn = λn fn , n = 1, 2, . . . , N,
there will correspond a positive frequency

T
ωn2 = λn , n = 1, 2, . . . , N,
md
and a solution of (6.4) of the form
un = (an cos(ωn t) + bn sin(ωn t)) fn
where an and bn are constants. This can now be easily verified by substituting this above
expression into the differential equation. To see we have enough constants of integration we
observe that we have two constants, an and bn , for each (vector) solution un . And we have
N vector solutions un —thus 2N constants in all. We now turn to an explicit evaluation of
the frequencies ωn2 —such frequencies are called normal modes.
6.3 Computation of the eigenvalues of VN

We introduce the characteristic polynomial of the matrix VN :
DN (λ) = DN = det (VN − λI) .
Expanding the determinant DN in the last column, we see that it is a sum of two terms—
each a determinant of matrices of size (N − 1) × (N − 1). One of these determinants equals
(2 − λ)DN −1 and the other equals DN −2 as is seen after expanding again, this time by the
last row. In this way one deduces
DN = (2 − λ)DN −1 − DN −2 , N = 2, 3, 4, . . .
with
D0 = 1 and D1 = 2 − λ.
We now proceed to solve this constant coefficient difference equation (in N ). From earlier
work we know that the general solution is of the form
c 1 µN N
1 + c 2 µ2
where µ1 and µ2 are the roots of
µ2 − (2 − λ)µ + 1 = 0.
Solving this quadratic equation gives

λ 1
µ1,2 = 1 − ± (2 − λ)2 − 4 .
2 2
It will prove convenient to introduce an auxiliary variable θ through
2 − λ = 2 cos θ,
A simple computation now shows

µ1,2 = e±iθ .
Thus
DN = c1 eiN θ + c2 e−iN θ .
To determine c1 and c2 we require that
D0 = 1 and D1 = 2 − λ.
That is,
c1 + c2 = 1,
c1 e + c2 e−iθ
iθ
= 2 − λ = 2 cos θ.
Solving for c1 and c2 ,
eiθ
c1 = ,
eiθ − e−iθ
e−iθ
c2 = − iθ .
e − e−iθ
Therefore,
1
i(N +1)θ −i(N +1)θ
DN = e − e
eiθ − e−iθ
sin ((N + 1)θ)
= .
sin θ
The eigenvalues of VN are solutions to
DN (λ) = det(VN − λI) = 0.
Thus we require
sin ((N + 1)θ) = 0,
which happens when
nπ
θ = θn := , n = 1, 2, . . . , N.
N +1
Thus the eigenvalues of VN are
λn = 2 − 2 cos θn = 4 sin2 (θn /2), n = 1, 2, . . . , N. (6.7)
The eigenfrequencies are
T 2T
ωn2 = λn = (1 − cos θn )
md md
2T nπ 4T 2 nπ
= 1 − cos = sin . (6.8)
md N +1 md 2(N + 1)
Remark: We know there are at most N distinct eigenvalues of VN . The index n does
not start at zero because this would imply θ = 0, but θ = 0—due to the presence of
sin θ in the denominator of DN —is not a zero of the determinant and hence does not
correspond to an eigenvalue of VN . We conclude there are N distinct eigenvalues of VN .
These eigenfrequencies are also called normal modes or characteristic oscillations.
6.3. COMPUTATION OF THE EIGENVALUES 95
Eigenvalues for N=50 particles

4
n
10 20 30 40 50
Figure 6.2: Eigenvalues λn , (6.7), for N = 50 particles.
6.4 The eigenvectors
6.4.1 Constructing the eigenvectors fn
We now find the eigenvector fn corresponding to eigenvalue λn . That is, we want a column
vector fn that satisfies
VN fn = 2(1 − cos θn )fn , n = 1, 2, . . . , N.
Setting,  
fn1
 fn2 
 
fn = 
 · ,

 · 
fnN
the above equation in component form is
−fn,j−1 + 2fn,j − fn,j+1 = 2(1 − cos θn )fn,j
with
fn,0 = fn,N +1 = 0.
This is a constant coefficient difference equation in the j index . Assume, therefore, a solution
of the form
fn,j = eijϕ .
The recursion relation becomes with this guess
−2 cos ϕ + 2 = 2(1 − cos θn ),
i.e.
ϕ = ±θn .
The fn,j will be linear combinations of e±ijθn ,
fn,j = c1 sin(jθn ) + c2 cos(jθn ).
We require fn,0 = fn,N +1 = 0 which implies c2 = 0.

To summarize,
md 2
VN fn = ω fn , n = 1, 2, . . . , N,
T n
2T nπ
ωn2 = (1 − cos θn ), θn = ,
md N +1
 
sin(θn ) (6.9)
 sin(2θn ) 
 
fn =  ·  n = 1, 2, . . . , N.
 
 · 
sin(N θn )
6.5. DETERMINATION OF CONSTANTS 97
The general solution to (6.4) is

N
u(t) = (an cos(ωn t) + bn sin(ωn t)) fn ,
n=1
or in component form,

N
uj (t) = (an cos(ωn t) + bn sin(ωn t)) sin(jθn ). (6.10)
n=1
Thus every oscillation of the weighted string is a sum of characteristic oscillations.
6.4.2 Orthogonality of eigenvectors
The set of eigenvectors {fn }N

n=1 forms a basis for R
N
since the matrix VN is symmetric.
(Another reason they form a basis is that the eigenvalues of VN are distinct.) We claim the
eigenvectors have the additional (nice) property that they are orthogonal, i.e.
fn · fm = 0, n = m,
where · denotes the vector dot product. The orthogonality is a direct result of the fact that
VN is a symmetric matrix. Another way to prove this is to use (6.9) to compute

N
fn · fm = sin(jθn ) sin(jθm ). (6.11)
j=1
To see that this is zero for n = m, we leave as an exercise to prove the trigonometric identity

N
njπ mjπ 1
sin sin = (N + 1)δn,m
j=1
N +1 N +1 2
where δn,m is the Kronecker delta function. (One way to prove this identity is first to use
the formula sin θ = (eiθ − e−iθ )/2i to rewrite the above sum as a sum of exponentials. The
resulting sums will be finite geometric series.) From this identity we also get that the length
of each vector, fn , is

N +1
fn = .
2
6.5 Determination of constants an and bn
Given the initial vectors u(0) and u̇(0), we now show how to determine the constants an
and bn . At t = 0,
N
u(0) = an fn .
n=1
Dotting the vector fp into both sides of this equation and using the orthogonality of the
eigenvectors, we see that

2
N
pjπ
ap = sin uj (0), p = 1, 2, . . . , N. (6.12)
N + 1 j=1 N +1
Differentiating u(t) with respect to t and then setting t = 0, we have

u̇(0) = ωn bn fn .
n=1
Likewise dotting fp into both sides of this equation results in

1
N
2 pjπ
bp = sin u̇j (0), p = 1, 2, . . . , N. (6.13)
N + 1 ωp j=1 N +1
If we assume the weighted string starts in an initial state where all the initial velocities are
zero,
u̇j (0) = 0,
then the solution u(t) has components

N
uj (t) = an cos(ωn t) sin(jθn ) (6.14)
n=1
where the constants an are given by (6.12) in terms of the initial displacements uj (0). The
special solutions obtained by setting all the an except for one to zero, are called the normal
modes of oscillation for the weighted string. They are most interesting to graph as a function
both in space (the j index) and in time (the t variable). In figures we show a “snapshot” of
various normal mode solutions at various times t.
6.5. DETERMINATION OF CONSTANTS 99
1 1
0.5 0.5
5 10 15 20 25 j 5 10 15 20 25 j
-0.5 -0.5
-1 -1
1 1
0.5 0.5
5 10 15 20 25 j 5 10 15 20 25 j
-0.5 -0.5
-1 -1
1 1
0.5 0.5
5 10 15 20 25 j 5 10 15 20 25 j
-0.5 -0.5
-1 -1
Figure 6.3: Vertical displacements uj for the two lowest (n = 1 and n = 2) normal modes are
plotted as function of the horizontal position index j. Each column gives the same normal
mode but at different times t. System is for N = 25 particles.
1 1
0.5 0.5
20 40 60 80100 j 20 40 60 80100 j
-0.5 -0.5
-1 -1
1 1
0.5 0.5
20 40 60 80100 j 20 40 60 80100 j
-0.5 -0.5
-1 -1
1 1
0.5 0.5
20 40 60 80100 j 20 40 60 80100 j
-0.5 -0.5
-1 -1
Figure 6.4: Vertical displacements uj for the two normal modes n = 5 and n = 10 are
plotted as function of the horizontal position index j. Each column gives the same normal
mode but at different times t. System is for N = 100 particles.
6.6. CONTINUUM LIMIT: THE WAVE EQUATION 101
6.6 Continuum limit: The wave equation

As the pictures illustrate, when the number of particles N becomes large and the distance d
between the particles becomes small, there appear limiting curves that describe the oscilla-
tions of the entire system. These limiting curves describe the oscillations of the string. Let
us pursue this in more detail. We assume
N → ∞ and d → 0 such that N d → L
where L is the length of the string (under no tension). We assume that the mass of the
string is given by µL where µ is the mass per unit length. Thus we assume
mN → µL
The positions of the particles, jd, j = 1, 2, . . . , N , are then assumed to approach a continuous
position variable x:
jd → x
We now examine the continuum limit of the system of ordinary differential equations
d2 uj T
= (uj−1 − 2uj + uj+1 )
dt2 md
To do this we assume there exists a function u(x, t) such that
uj (t) = u(jd, t)
Then, since d is small,
∂u 1 ∂2u
uj−1 = u(jd − d, t) = u(x, t) − d (x, t) + d2 2 (x, t) + O(d3 )
∂x 2 ∂x
and similarly
∂u 1 ∂2u
uj+1 = u(jd + d, t) = u(x, t) + d (x, t) + d2 2 (x, t) + O(d3 )
∂x 2 ∂x
and hence
∂ 2u
uj−1 − 2uj + uj+1 = d2 (x, t) + O(d3 )
∂x2
Substituting this into our differential equations we obtain
∂2u T ∂2u
=
∂t2 µ ∂x2
Note that since m = µL/N ,

T d2 Td T dN T
= = =
md m µL µ
Also observe that T /µ has the dimensions of (velocity)2 . Thus let’s call
T
v2 =
µ
so that we have
∂2u 1 ∂2u
− = 0. (6.15)
∂x2 v 2 ∂t2
This is the one-dimensional wave equation. It is an example of a partial differential equation.

Given our analysis of the weighted string, we can anticipate that if we studied solutions of
the single partial differential equation (6.15), then u = u(x, t) would describe the oscillations
of a string. Note that we would have the two boundary conditions
u(0, t) = u(L, t) = 0
which corresponds to the statement that the string is tied down at x = 0 and at x = L for all
times t. In addition, we specify at t = 0 the initial displacement of the string: u(x, 0) = f (x)
where f is a given function as well as the initial velocity ∂u∂t (x, 0). The problem then is to
find the solution to (6.15) satisfying these conditions. In the next section we show how the
methods we’ve developed so far permit us to find such a solution.
6.6.1 Solution to the wave equation
We first look for solutions of the form (called separation of variables)
u(x, t) = X(x) T (t)
where X is only a function of x and T is only a function of t. Since
∂2u d2 X ∂2u d2 T
= T (t) and = X(x) ,
∂x2 dx2 ∂t2 dt2
we have, upon substituting these expressions into (6.15) and dividing by X T the condition
1 d2 X 1 1 d2 T
= .
X dx2 v 2 T dt2
The left-hand side of the above equation is a function only of x and the right-hand side of
the same equation is a function only of t. The only way this can be true is for both sides to
equal the same constant. (We will see below that this constant has to be negative to satisfy
the boundary conditions. Anticipating this fact we write the constant as −k 2 .) That is to
say, we have
1 d2 X 1 1 d2 T
2
= −k 2 = 2
X dx v T dt2
This gives us two ordinary differential equations:
d2 X 2 d2 T
+ k X = 0, + k 2 v 2 T = 0.
dx2 dt2
The solution to the first equation is
X(x) = c1 cos(kx) + c2 sin(kx).

6.6. CONTINUUM LIMIT: THE WAVE EQUATION 103
We want u(0, t) = 0 which implies c1 = 0. We also require u(L, t) = 0. If we set c2 = 0 then

X is identically zero and we have the trivial solution. Thus we must require
sin(kL) = 0.
This is satisfied if
kL = nπ, n = 1, 2, 3, . . . .
(Note that n = −1, −2, . . . give the same solution up to a sign and n = 0 corresponds to
X identically zero.) The solution to the T equation is also a linear combination of sines
and cosines. Thus for each value of n we have found a solution satisfying the conditions
u(0, t) = u(L, t) = 0 of the form
nπ nπv nπv !
un (x, t) = sin( x) an cos( t) + bn sin( t)
L L L
where an and bn are constants. Since the wave equation is linear, linear supposition of
solutions results in a solution. Thus
nπv !
∞
nπ nπv
u(x, t) = sin( x) an cos( t) + bn sin( t)
n=1
L L L
is a solution satisfying u(0, t) = u(L, t) = 0. We now require that u(x, 0) = f (x). That is
we want
∞
nπ
u(x, 0) = an sin( x) = f (x)
n=1
L
We now use the fact that the for m, n = 1, 2, . . .
L
mπ nπ L
sin( x) sin( x) dx = δm,n
0 L L 2
to find L
2 nπ
an = f (x) sin( x) dx. (6.16)
L 0 L
This determines the constants an . If we further assume (for simplicity) that
∂u
(x, 0) = 0
∂t
(initial velocity is zero), then a very similar calculation gives bn = 0. Thus we have shown
∞
nπ nπv
u(x, t) = an sin( x) cos( t) (6.17)
n=1
L L
where an are given by (6.16).

It is instructive to compare this solution of the wave equation to the solution (6.14) of
the weighted string. We take the N → ∞ limit directly in (6.14) and use the same scaling
as we have in the above derivation of the wave equation. In this limit we can replace
L µL xN
d −→ , m −→ , j −→
N N L
Thus
4T nπ 4T n2 π 2 T n2 π 2
ωn2 = sin2 ( )∼ ∼
md 2(N + 1) md 4(N + 1)2 µ L2
so that
nπ
ωn −→ v .
L

(Recall the definition v = T /µ.) Similarly,
njπ N nπ nπ
jθn = = x −→ x.
N +1 N +1 L L
Putting these limiting expressions into (6.14) and taking the N → ∞ limit we see that (6.14)
becomes (6.17). The only point that needs further checking is to show the an as given by
(6.12) approaches the an as given by (6.16). This requires the natural assumption that the
initial conditions uj (0) can be written in the form uj (0) = f (jd) for some smooth function
f . This is the f of u(x, 0) = f (x). A calculation then shows that (6.12) is the Riemann sum
approximation to (6.16) and approaches (6.16) as N → ∞.
The take home message is that the oscillations described by the solution to the wave
equation can be equivalently viewed as an infinite system of harmonic oscillators.
6.7 Inhomogeneous problem

The inhomogeneous version of (6.4) is
d2 u T
+ VN u = F(t) (6.18)
dt2 md
where F(t) is a given driving term. The j th component of F(t) is the external force acting
on the particle at site j. An interesting case of (6.18) is
F(t) = cos ωt f
where f is a constant vector. The general solution to (6.18) is the sum of a particular
solution and a solution to the homogeneous equation. For the particular solution we assume
a solution of the form
up (t) = cos ωt g.
Substituting this into the differential equation we find that g satisfies

md 2 md
VN − ω I g= f.
T T
For ω 2 = ωn2 , n = 1, 2, . . . , N , the matrix

md 2
VN − ω I
T
is invertible and hence −1
md md 2
g= VN − ω I f.
T T
6.8. VIBRATING MEMBRANE 105
Writing (possible since the eigenvectors form a basis)

N
f= αn fn ,
n=1
we conclude that

N
αn
g= f
2 − ω2 n
ω
n=1 n
for ω 2 = ωn2 , n = 1, 2, . . . , N . The solution with initial values
u(0) = 0, u̇(0) = 0 (6.19)
is therefore of the form

N
αn
u(t) = cos ωt fn + (an cos(ωn t) + bn sin(ωn t)) fn .
ω2
n=1 n
−ω 2
n=1
Imposing the initial conditions (6.19) we obtain the two equations

αn
+ a n fn = 0, (6.20)
n=1
ωn2 − ω 2

ωn bn fn = 0. (6.21)
n=1
From the fact that {fn }N

n=1 is a basis we conclude
αn
an = − , bn = 0 for n = 1, 2, . . . , N.
ωn2 − ω2
Thus the solution is

N
αn
u(t) = 2 − ω2
(cos(ωt) − cos(ωn t)) fn (6.22)
ω
n=1 n
N
2αn 1 1
= sin (ωn + ω)t sin (ωn − ω)t fn . (6.23)
ω2 − ω2
n=1 n
2 2
We observe that there is a beat whenever the driving frequency ω is close to a normal mode
of oscillation ωn . Compare this discussion with that of Boyce & DiPrima [4].
6.8 Vibrating membrane
6.8.1 Helmholtz equation
In the previous section we discussed the vibrating string. Recall that we have a string of
unstretched length L that is tied down at ends 0 and L. If u = u(x; t) denotes the vertical
displacement of the string at position x, 0 ≤ x ≤ L, at time t, then we showed that for small
displacements u satisfies the one-dimensional wave equation
∂2u 1 ∂2u
− =0
∂x2 v 2 ∂t2
where v 2 = T /µ, T equals the tension in string and µ is the density of the string. We solved
this equation subject to the boundary conditions u(0, t) = u(L, t) = 0 for all t and with
initial conditions u(x, 0) = f (x) and ∂u
∂t (x, 0) = g(x) where f and g are given.
Now we imagine a uniform, flexible membrane, of mass ρ per unit area, stretched under
a uniform tension T per unit length over a region Ω in the plane whose boundary ∂Ω is a
smooth curve (with a possible exception of a finite number of corners).
We now let U = U (x, y; t) denote the vertical displacement of the membrane at position
(x, y) ∈ Ω at time t from its equilibrium position. We again assume that the membrane is
tied down at the boundary; that is1
U (x, y; t) = 0 for (x, y) ∈ ∂Ω.
The motion of U = U (x, y; t) is governed by the two-dimensional wave equation:
∂2U ∂2U 1 ∂2U

+ − = 0 for (x, y) ∈ Ω (6.24)
∂x2 ∂y 2 v 2 ∂t2
∂2 U ∂2U
where v 2 = T /ρ. One recognizes ∂x2 + ∂y 2 as the two-dimensional Laplacian. So if we
introduce
∂2 ∂2
∆= +
∂x2 ∂y 2
the wave equation takes the form
1 ∂2U
∆U − = 0.
v 2 ∂t2
We proceed as before and look for solutions of (6.24) in which the variables separate
U (x, y; t) = u(x, y)T (t).
Substituting this into (6.24), and then dividing by uT gives
1 1 1 d2 T
∆u = 2 .
u v T dt2
The right-hand side depends only upon t where as the left-hand side depends only upon x, y.
Thus for the two sides to be equal they must equal the same constant. Call this constant
−k 2 . Thus we have the two equations
1 In one dimension Ω = (0, L) and the boundary of Ω consists of the two points 0 and L.
d2 T
+ ω2T = 0 where ω = kv,
dt2
∆u + k 2 u = 0. (6.25)
The differential equation for T has our well-known solutions
eiωt and e−iωt .
The second equation (6.25), called the Helmholtz equation, is a partial differential equation
for u = u(x, y). We wish to solve this subject to the boundary condition
u(x, y) = 0 for (x, y) ∈ ∂Ω.
6.8.2 Rectangular membrane

Consider the rectangular domain
Ω = {(x, y) : 0 < x < a, 0 < y < b} (6.26)
For this rectangular domain the Helmholtz equation can be solved by the method of sepa-
ration of variables. If one assumes a solution of the form (variables x and y separate)
u(x, y) = X(x)Y (y)
then the problem is reduced to two one-dimensional problems. It is an exercise to show that
the allowed frequencies are
1/2
m 2 n 2
ωm,n = πv + , m, n = 1, 2, 3, . . . (6.27)
a b
6.8.3 Circular membrane: The drum

We now consider the circular domain

Ω = (x, y) : x2 + y 2 < a2
so the boundary of Ω is the circle ∂Ω : x2 + y 2 = a2 . Even though the variables separate in
the Cartesian coordinates x and y, this is of no use since the boundary is circular and we
would not be able to apply the BC u = 0 on the circular boundary ∂Ω. Since the domain is
circular it is natural to introduce polar coordinates
x = r cos θ, y = r sin θ.
It is an exercise in the chain rule to show that in polar coordinates the 2D Laplacian is
∂ 2 u 1 ∂u 1 ∂2u
∆u =
2
+ + 2 ;
∂r r ∂r r ∂θ2
and hence, the Helmholtz equation in polar coodinates is
∂ 2 u 1 ∂u 1 ∂2u
+ + + k 2 u = 0. (6.28)
∂r2 r ∂r r2 ∂θ2
We write u = u(r, θ).2
2 Possible point of confusion: We wrote u = u(x, y) so really our new function of r and θ is u(r cos θ, r sin θ).
Technically we should give this function of r and θ a new name but that would be rather pedantic.
Separation of variables
We now show that the variables separate. So we look for a solution of the form
u(r, θ) = R(r)Θ(θ).
Substituting this into (6.28), multiplying by r2 /RΘ we have
r2 d2 R r dR 1 d2 Θ
2
+ + k 2 r2 = −
R dr R dr Θ dθ2
By the now familiar argument we see that each of the above sides must equal a constant,
call it m2 , to obtain the two differential equations
d2 Θ
+ m2 Θ = 0 (6.29)
dθ2
d2 R 1 dR m2
2
+ + (k 2 − 2 )R = 0 (6.30)
dr r dr r
Two linearly independent solutions to (6.29) are
eimθ and e−imθ
The point with polar coordinates (r, θ) is the same point as the one with polar coordinates
(r, θ + 2π). Thus our solution u(r, θ) and u(r, θ + 2π) must be the same solution. This
requires
eimθ+im2π = eimθ
or e2πim = 1. That is, m must be an integer. If m = 0 the general solution to (6.29) is
c1 + c2 θ. But the θ → θ + 2π argument requires we take c2 = 0. Thus the general solution
to (6.29) is
am cos(mθ) + bm sin(mθ), m = 0, 1, 2, . . .
We now return to (6.30), called the Bessel equation, which is a second order linear differential
equation. General theory tells us there are two linearly independent solutions. Tradition has
it we single out two solutions. One solution, called Jm (kr), is finite as r → 0 and the other
solution, called Ym (kr) goes to infinity as r → 0. Both of these functions are called Bessel
functions. It can be shown that the Bessel function Jm (z) is given by the series expansion
z m
∞
1 z 2j
Jm (z) = (−1)j (6.31)
2 j=0
j!(m + j)! 2
A plot of the Bessel function J0 (x) for 0 ≤ x ≤ 40 is given in Figure 6.5. In Mathematica,
Bessel functions Jm (z) are called by the command BesselJ[m,z]. Since u(r, θ) is well-
defined at r = 0 (center of the drum), this requires we only use the Jm solutions. Thus we
have shown that
Jm (kr) (am cos(mθ) + bm sin(mθ)) , m = 0, 1, 2, . . .
are solutions to (6.28). We now require that these solutions vanish on ∂Ω. That is, when
r = a and for all θ we require the above solution to vanish. This will happen if
Jm (ka) = 0.
Bessel Function J0
1.0
0.8
0.6
0.4
0.2
10 20 30 40
0.2
0.4
Figure 6.5: The Bessel function J0 (x). First zero occurs at approximately 2.4048, the second
zero at 5.5201, the third zero at 8.6537, . . . .
That is we have to be at a zero of the Bessel function Jm . It is known that Jm has an infinite
number of real zeros, call them jm,n , n = 1, 2, . . .. Thus the frequencies that the drum can
oscillate at are
v
ωm,n = jm,n , m = 0, 1, 2, . . . ; n = 1, 2, . . .
a
where jm,n is the nth zero of the Bessel function Jm (z). These zeros can be found in
Mathematica using the command BesselJZero[m,n].
6.8.4 Comments on separation of variables
For general domains Ω one cannot solve the Helmholtz equation (6.25) by the method
of separation of variables. In general if one makes the transformations x = f (ξ, η) and
y = g(ξ, η) then one would want the curves of constant ξ (or constant η) to describe the
boundary ∂Ω and for Helmholtz’s equation to separate variables in the new variables ξ and
η. In general there are no such coordinates. For an elliptical membrane the Helmholtz
equation does separate in what are called elliptic coordinates
c c
x= cosh µ cos θ, y = sinh µ sin θ
2 2
where c ∈ R+ , 0 < µ < ∞ and 0 ≤ θ ≤ 2π. The curves µ = constant and θ = constant are
confocal ellipses and hyperbolas, respectively. Qualitative new phenomena arise for elliptical
(and more generally convex) membranes: the existence of whispering gallery modes and
bouncing ball modes. In the whispering gallery mode the eigenfunction is essentially nonzero
only in a thin strip adjacent to the boundary of Ω. Thus a person who speaks near the wall
of a convex room can be heard across the room near the wall, but not in the interior of the
room. For further information see [5] and references therein.
6.9 Exercises
#1. Weighted string on a circle
We consider the same weighted string problem but now assume the masses lie on a circle;
this means that the first mass is coupled to the last mass by a string. The effect of this
is that (6.1) remains the same if we now interpret u0 = uN and uN +1 = u1 . Explain why
this is the case. What is the matrix VN in this case? Show that the differential equations
can still be written in the matrix form (6.4) where now the VN is your new VN . Does the
reduction to an eigenvalue problem, as in §6.2, remain the same? Explain.
#2. Diagonalization of VN from problem #1
Let VN be the N × N matrix found in the previous problem. Show that the eigenvalue
problem
VN f = λf
becomes in component form
−fj−1 + 2fj − fj+1 = λfj , j = 1, 2, . . . , N (6.32)
where f0 = fN and fN +1 = f1 . Let ω denote an N th root of unity; that is, any of the values
e2πin/N , n = 0, 1, . . . , N − 1. For each such choice of ω, define

N
fˆω = fj ω j (6.33)
j=1
Multiply (6.32) by ω j and sum the resulting equation over j = 1, 2, . . . , N . Show that the
result is
2(1 − cos φ)fˆω = λfˆω
where ω = eiφ . From this we conclude that the eigenvalues are

2πn
λn = 2 1 − cos , n = 0, 1, . . . , N − 1
N
Explain why this is so. This should be compared with (6.7). Find an eigenvector fn corre-
sponding to eigenvalue λn . (Hint: Follow the method in §6.4.1.)
6.9. EXERCISES 111
#3. Coupled pendulums
Figure 6.6: Coupled pendulums. Here we assume the damping force is zero.
Consider the system of two mathematical pendulums of lengths 1 and 2 and masses m1
and m2 , respectively, in a gravitional field mg which move in two parallel vertical planes
perpendicular to a common flexible support such as a string from which they are suspended.
Denote by θ1 (θ2 ) the angle of deflection of pendulum #1 (#2). The kinetic energy of this
system is
1 1
KE = m1 21 θ̇12 + m2 22 θ̇22 ,
2 2
and the potential energy is
PE = m1 g 1 (1 − cos θ1 ) + m2 g 2 (1 − cos θ2 ) + Vint
where Vint is the interaction potential energy.3 If there is no twist of the support, then there
is no interaction of the two pendulums. We also expect the amount of twist to depend upon
the difference of the angles θ1 and θ2 . It is reasonable to assume Vint to be an even function
of θ1 − θ2 . Thus

Vint (0) = 0, Vint (0) = 0.
3 These expressions should be compared with (2.26).
For small deflection angles (the only case we consider) the simplest assumption is then to
take
1
Vint (θ1 − θ2 ) = κ(θ1 − θ2 )2
2
where κ is a positive constant. Since we are assuming the angles are small, the potential
energy is then given, to a good approximation, by
1 1 1
PE = m1 g1 θ12 + m2 g2 θ22 + κ(θ1 − θ2 )2 .
2 2 2
Under these assumptions it can be shown that Newton’s equations are

m1 21 θ̈1 = −(m1 g1 + κ)θ1 + κθ2 ,
m2 22 θ̈2 = κθ1 − (m2 g2 + κ)θ2 .
Observe that for κ = 0 the ODEs reduce to two uncoupled equations for the linearized
mathematical pendulum. To simplify matters somewhat, we introduce
g g κ κ
ω12 = , ω2 = k1 = , k2 = .
1 2 m1 21 m2 22
Then it is not difficult to show (you need not do this) that the above differential equations
become
θ̈1 = −(ω12 + k1 )θ1 + k1 θ2

(6.34)
θ̈2 = k2 θ1 − (ω22 + k2 )θ2 .
We could change this into a system of first order DEs (the matrix A would be 4 × 4).
However, since equations of this form come up frequently in the theory of small oscillations,
we proceed to develop a “mini theory” for these equations. Define

θ1
Θ= .
θ2
Show that the equations (6.34) can be written as
Θ̈ = AΘ (6.35)
where A is a 2 × 2 matrix. Find the matrix A. Assume a solution of (6.35) to be of the form

a1
Θ(t) = eiωt . (6.36)
a2
Using (6.36) in (6.35) show that (6.35) reduces to
AΘ = −ω 2 Θ. (6.37)
This is an eigenvalue problem. Show that ω 2 must equal
2 1 2
ω± = (ω + ω22 + k1 + k2 )
2 1

1
± (ω12 − ω22 )2 + 2(ω12 − ω22 )(k1 − k2 ) + (k1 + k2 )2 . (6.38)
2
6.9. EXERCISES 113
2
Show that an eigenvector for ω+ is

1
f1 = 2 , (6.39)
−k2 (ω+ − ω22 − k2 )−1
2
and an eigenvector corresponding to ω− is
2

−k1 (ω− − ω12 − k1 )−1
f2 = . (6.40)
1
Now show that the general solution to (6.34) is

θ1 (t)
= (c1 cos(ω+ t) + c2 sin(ω+ t)) f1 + (c3 cos(ω− t) + c4 sin(ω− t)) f2 (6.41)
θ2 (t)
where ci are real constants. One can determine these constants in terms of the initial data
θ1 (0), θ̇1 (0), θ2 (0), θ̇2 (0).
To get some feeling for these rather complicated expressions, we consider the special case
θ1 (0) = θ0 , θ̇1 (0) = 0, θ2 (0) = 0, θ̇2 (0) = 0 (6.42)
with
m1 = m2 = m, 1 = 2 = . (6.43)
These last conditions imply
ω1 = ω2 := ω0 .
Explain in words what these initial conditions, (6.42), correspond to in the physical set
up.
If we define
κ
k= ,
m 2
show that in the special case (6.42) and (6.43) that

ω+ = ω02 + 2k and ω− = ω0 . (6.44)
In this same case solve for the coefficients c1 , c2 , c3 and c4 and show that
1 1
c1 = θ0 , c2 = 0, c3 = θ0 , c4 = 0,
2 2
and hence (6.41) becomes

1 1
θ1 (t) = θ0 cos (ω+ + ω− )t cos (ω+ − ω− )t ,
2 2

1 1
θ2 (t) = θ0 sin (ω+ + ω− )t sin (ω+ − ω− )t .
2 2
Suppose further that

k
1. (6.45)
ω02
What does this correspond to physically? Under assumption (6.45), show that approximately

k
θ1 (t) ≈ θ0 cos (ω0 t) cos t ,
2ω0

k
θ2 (t) ≈ θ0 sin (ω0 t) sin t . (6.46)
2ω0
Discuss the implications of (6.46) in terms of the periods

2π 2π
T0 = and T1 = .
ω0 k/2ω0
Show that in this approximation
T1 T0 .
Draw plots of θ1 (t) and θ2 (t) using the approximate expressions (6.46).
#4. The Toda chain and Lax pairs
Consider N particles on a circle (periodic boundary conditions) whose positions xn (t) at

time t satisfy the Toda equations
d2 xn
= exp (−(xn − xn−1 )) − exp (−(xn+1 − xn )) , n = 1, 2, . . . , N, (6.47)
dt2
where xN +1 = x1 and x0 = xN . These equations are nonlinear and admit certain solutions,
called solitons, which are stable pulses. This system of equations has been extensively
studied. Here we give only a brief introduction to some of these results.4
To make the problem easier we now set N = 5 but everything that follows can be
generalized to any positive integer N .
Define
1 1 dxn
exp (−(xn+1 − xn )/2) and bn =
an = , n = 1, . . . , 5. (6.48)
2 2 dt
Show that if xn satisfies the Toda equations (6.47), then an and bn satisfy the differential
equations
dan dbn
= an (bn − bn+1 ) and = 2 a2n−1 − a2n . (6.49)
dt dt
Define two 5 × 5 matrices L and B, they are called a Lax pair, by
   
b1 a1 0 0 a5 0 −a1 0 0 a5
 a1 b 2 a2 0 0   a1 0 −a 0 0 
   2 
L= 0 a2 b 3 a3 0   and B = 
 0 a 2 0 −a 3 0 .
 (6.50)
 0 0 a3 b 4 a4   0 0 a3 0 −a4 
a5 0 0 a4 b 5 −a5 0 0 a4 0
4 See, for example, Theory of Nonlinear Lattices by Morikazu Toda, Springer-Verlag, 1981.
6.9. EXERCISES 115
Show (6.49) can be written as the matrix equation

dL
= BL − LB (6.51)
dt
Define the matrix U = U (t) to be the solution to the differential equation
dU
= BU
dt
that satisfies the initial condition U (0) = I.
Show that U (t) is a unitary matrix; that is, U ∗ (t)U (t) = I for all t where U ∗ is the
adjoint matrix.5 Hint: Observe that B ∗ = −B. Use this to first show that
dU ∗
= −U ∗ B
dt
d ∗
and then show dt U (t)U (t) = 0.
Now prove that
d
(U ∗ L(t)U (t)) = 0
dt
and hence that
U ∗ (t)L(t)U (t) = L(0)
That is, L(0) and L(t) are unitarily equivalent. From this conclude
The eigenvalues of L(t) are independent of t
Thus the eigenvalues of the Lax matrix L are first integrals of motion of the Toda chain. For
general N this means that we have found N integrals of the motion. This is a remarkable
result since normally one can only find a limited number of integrals of the motion (energy,
angular momentum, etc.).
#5. Wave equation
In the section “Solution to the Wave Equation” it was claimed that a similar argument
shows that the coefficients bn are equal to zero. (See discussion between (6.16) and (6.17).)
Prove that bn = 0.
#6. Weighted string with friction
We now assume that the particles in the weighted string problem are subject to a force due
to the presence of friction. (Imagine the particles are moving in a medium which offers
resistance to the motion of the particles.) Assuming the frictional force is proportional to
the velocity, the system of differential equations describing the motion is
d2 uj T duj
m = (uj+1 − 2uj + uj−1 ) − γ , j = 1, 2, . . . , N (6.52)
dt2 d dt
where γ is positive and, as before, u0 = uN +1 = 0.
5 Recall that if X is any matrix then X ∗ is the matrix obtained by taking the complex conjugate of each
element in X and then taking the transpose.

1. Rewrite the system (6.52) in matrix form such that when γ = 0 the equation becomes
identical to the matrix equation (6.4).
2. Assume a solution of the form
u(t) = eiωt f (6.53)
where f is a column vector independent of t and ω is to be determined. For what
values of ω is (6.53) a solution to the matrix equation derived in part (1)?
Note: This will not require a complete reworking of the eigenvalues since you may use
the information we already have proved about VN to find the eigenvalues in this new
problem. You should not have to solve anything more complicated than a quadratic
equation.
3. Explain the significance of the fact that the ω’s you obtain are complex numbers.
4. For a large system N 1 explain why you expect some of the allowed ω’s to be purely
imaginary. Explain the significance of this result, i.e. what is the implication for the
motion?
#7. Rectangular membrane
Figure 6.7: Four lowest normal modes of a rectangular drum.
In this section we obtain the solution of (6.25) in the case of a rectangular domain (6.26).
1. By assuming that the solution can be written as u(x, y) = X(x)Y (y), obtain a 2nd
order DE for X with independent variable x and similarly a DE for Y with independent
variable y.
6.9. EXERCISES 117
2. We assume the membrane is tied down at the boundary of the domain Ω. (This implies
boundary conditions on the solutions we seek.)
3. Show that the eigenvalues and the corresponding eigenfunctions of the differential
equations with boundary conditions in parts (1) and (2) are
m2 π 2 mπx
µm = ; X m (x) = Am sin , m = 1, 2, · · · (6.54)
a2 a
n2 π 2 nπy
νn = ; Yn (y) = B n sin , n = 1, 2, · · · (6.55)
b2 b
4. Show that the eigenfrequencies (normal modes) of the rectangular membrane are given
by (6.27). (By dimensional analysis conclude where the factor v, which was set equal
to one here, must appear.)
5. Find the general solution to (6.25) for this rectangular domain.
#8. Alternating mass-spring: Acoustic and optical phonons
Consider 2N particles on a circle interacting via a spring connnecting adjacent particles.

We assume the particles on the odd sites have mass m1 and the particles on the even sites
have mass m2 . If uj denotes the displacement from equilibrium of particle j, the differential
equations describing the motion are
d2 uj
mj + k (−uj−1 + 2uj − uj+1 ) = 0 for j = 1, 2, 3, . . . , 2N, (6.56)
dt2
where because the particles are on a circle
u2N +1 = u1 and u0 = u2N .
Here k is the spring constant for the spring connecting any two particles. We are interested
in finding the frequencies at which the system can oscillate.
1. Assume a solution of the form

uj (t) = eiωt vj , vj independent of t,
and show that (6.56) becomes
−mj ω 2 vj + k (−vj−1 + 2vj − vj+1 ) = 0 for j = 1, 2, 3, . . . , 2N, (6.57)
2. For j = 1, 2, . . . , N define the vectors

v2j−1
Vj = .
v2j
Show that (6.57) can be written equivalently as
"
m1 0 0 1 2 −1 0 0
− ω 2 Vj +k − Vj−1 + Vj − Vj+1 =0
0 m2 0 0 −1 2 1 0
(6.58)
for j = 1, 2, . . . , N .
3. Let η denote any N th root of unity, i.e. η N = 1 so η is of the form η = eiφ = e2πij/N
for some integer j = 0, 1, . . . , N − 1. Define

N
V̂η = Vj η j
j=1
Show that V̂η satisfies the equation

"
m1 0 0 1 2 −1 0 0
− ω2 + k − η+ − η −1 V̂η = 0.
0 m2 0 0 −1 2 1 0
(6.59)
4. What is the condition for nontrivial solutions V̂η to exist for (6.59)? Hint: Equation
(6.59) is of the form AV̂η = 0 where the matrix A is the 2 × 2 matrix inside the curly
brackets of (6.59). Using the condition you just found, show that the normal modes
of vibration are given by
1/2
2 1 1 1 1 2 2πj
ω±,j =k + ±k + 2+ cos (6.60)
m1 m2 m21 m2 m1 m2 N
where j = 0, 1, 2, . . . N − 1.
5. Show that the frequencies derived in (6.60) lie on two curves, called dispersion curves.
These two curves should be compared with the one dispersion curve for the equal mass
problem. Plot the two dispersion curves.6 The curve that is zero at j = 0 is called
the acoustic mode and the other is called the optical mode.7 This is a model of a
one-dimensional lattice vibrations of a diatomic system.
#9. Energy of the vibrating string
The vibrating string has the total energy E(t) at time t

L 1
1
E(t) = µ u2t (x, t) + T u2x (x, t) dx
0 2 2
Explain why the first term is the kinetic energy and the second term is the potential energy
of the vibrating string. Recall the solution u(x, t) of the vibrating string problem, i.e. (6.17).
Above we use the notation
∂u(x, t) ∂u(x, t)
ut (x, t) := and ux (x, t) := .
∂t ∂x
You may assume as given the following integrals:
L mπ nπ 1
sin x sin x dx = L δm,n (6.61)
0 L L 2
6 Inplotting you might want to fix some values of m1 , m2 and k.
7 The acoustic modes correspond to sound waves in the lattice. The optical modes, which are nonzero at
j = 0, are called “optical” because in ionic crystals they are excited by light. The quantized version of these
excitations are called acoustic phonons and optical phonons.
6.9. EXERCISES 119
and L mπ nπ 1
cos x cos x dx = L δm,n (6.62)
0 L L 2
Use (6.61) and (6.62) to show
∞
π2 T 2 2
E(t) = n an . (6.63)
4L n=1
Note that the result is independent of t, i.e. the energy of the vibrating string is conserved.
Give a physical interpretation of this expression for E in terms of harmonic oscillators.
Chapter 7
Quantum Harmonic Oscillator
Figure 7.1: Erwin Schrödinger, 1887–1961 and Paul Dirac, 1902–1984.
A simple and interesting example of a dynamical system in quantum mechanics

is the harmonic oscillator. This example is of importance for general theory,
because it forms a corner-stone in the theory of radiation.
P. A. M. Dirac, The Principles of Quantum Mechanics
121
122 CHAPTER 7. QUANTUM HARMONIC OSCILLATOR
7.1 Schrödinger equation
In classical mechanics the state of a system consisting of N particles is specified by the

position x and momentum p = mv of each particle. The time evolution of this state is de-
termined by solving Newton’s Second Law (or equivalently, say Hamilton’s equations). Thus,
for example, a one particle system moving in three-dimensions (three degrees of freedom)
determines a curve in 6-dimensional space: namely, (x(t), p(t)). For the familiar harmonic
oscillator (mass-spring system) there is only one-degree of freedom (the movement of the
mass is in one dimension only) and the position and momentum are given by the now familiar
formulas1
p0
x(t) = x0 cos(ω0 t) + sin(ω0 t), (7.1)
mω0
p(t) = p0 cos(ω0 t) − mω0 x0 sin(ω0 t). (7.2)
In quantum mechanics the notion of the state of the system is more abstract. The state
is specified by a vector Ψ in some abstract vector space H. This vector space has an inner
product (·, ·).2 Thus every state Ψ ∈ H satisfies
1/2
Ψ := (Ψ, Ψ) < ∞. (7.3)
The importance of (7.3) is that in the Born interpretation |(Ψ, Φ)|2 is interpreted as a
probability; and hence, must be finite (and less than or equal to one).3 In what is called the
Schrödinger representation, one can describe the state Ψ as a function Ψ(x) where x is the
position (of say the particle). Then |Ψ(x)|2 is the probability density of finding the particle
in some small neighborhood of the point x. Integrating this over all possible positions must
then give one.
The evolution of the state Ψ with time is determined by solving the Schrödinger equation:
∂Ψ
i = HΨ. (7.4)
∂t
Here is the Planck’s constant4 (divided by 2π) and H is the quantum mechanical Hamil-
tonian, a linear self-adjoint operator on the space H.5
1 Actually, the second may look a little different from the earlier formulas. This is due to the fact that
we are using momentum p instead of velocity v to describe the second coordinate of (x, p). Here p0 is the
initial momentum and is related to the initial velocity by p0 = mv0 .
2 Such vector spaces are called Hilbert spaces.
3 This assumes that states Ψ are normalized so that their “length” is one, i.e. Ψ = 1.
4 In the cgs system, = 1.05457 × 10−27 erg-sec. A quantity that has the units of energy×time is called
an action. In modern particle physics a unit system is adopted such that in these units = 1. Max Planck
received the Nobel prize in 1919 for “his discovery of energy quanta”.
5 An operator H is self-adjoint if (Hψ, ψ) = (ψ, Hψ) for all ψ ∈ H. It is the generalization to Hilbert
spaces of the notion of a Hermitian matrix. There are some additional subtle questions regarding the domain
of the operator H. In these notes we ignore such questions and assume H is well-defined on all states Ψ ∈ H.
7.2. HARMONIC OSCILLATOR 123
7.2 Harmonic oscillator

7.2.1 Harmonic oscillator equation
We illustrate the notions of quantum mechanics and its relationship to differential equations
in the context of the harmonic oscillator. The harmonic oscillator is one of the most impor-
tant simple examples in quantum mechanics. In this case the vector space H is the space of
square-integrable functions. This space consists of all (complex valued) functions ψ(x) such
that ∞
2
|ψ(x)| dx < ∞.
−∞
This space is denoted by L2 (R) and it comes equipped with an inner product
∞
(ψ, ϕ) = ψ(x)ϕ̄(x) dx
−∞
where ϕ̄ is the complex conjugate of ϕ. (Note that in most physics books the complex
conjugation is on the first slot.) The first observation, and an important one at that, is
that the state space is infinite dimensional. For example, it can be proved that the infinite
sequence of functions
2
xj e−x , j = 0, 1, 2 . . .
are linearly independent elements of L2 (R). Thus in quantum mechanics one quickly goes
beyond linear algebra which is traditionally restricted to finite-dimensional vector spaces.
The operator H which describes the harmonic oscillator can be defined once we give
the quantization procedure—a heuristic that allows us to go from a classical Hamiltonian to
a quantum Hamiltonian. As mentioned above, classically the state is given by the vector
(x, p) ∈ R2 . In quantum mechanics the position and momentum are replaced by operators
x̂ and p̂. For the vector space of states H = L2 (R), the position operator acts on L2 (R) by
multiplication,
(x̂ψ)(x) = xψ(x), ψ ∈ L2 (R)
and the momentum operator p̂ acts by differentiation followed by multiplication by the
constant −i,
∂ψ
(p̂ψ)(x) = −i (x), ψ ∈ L2 (R).
∂x
Since x̂ is multiplication by x we usually don’t distinguish between x and x̂. From this we
observe that in quantum mechanics the position operator and the momentum operator do
not commute. To see this, let ψ ∈ L2 (R), then
∂ψ ∂
(x̂p̂ − p̂x̂) ψ(x) = −ix + i (xψ(x))
∂x ∂x
∂ψ ∂ψ
= −ix + ix + i ψ(x)
∂x ∂x
= i ψ(x).
Introducing the commutator; namely, for any two operators A and B we define [A, B] =
AB − BA, the above can be written more compactly as6
[x̂, p̂] = i id (7.5)
6 Just as in linear algebra, if A and B are two linear operators and it holds for all vectors ψ that Aψ = Bψ,
then we can conclude that as operators A = B.

where by id we mean the identity operator. Equation (7.5) is at the heart of the famous
Heisenberg Uncertainty Relation.
With these rules we can now define the quantum harmonic oscillator Hamiltonian given
the classical Hamiltonian (energy). Classically,7
E = KE + PE
1 2 1
= p + mω02 x2 .
2m 2
Replacing p → p̂ and x by multiplication by x we have
2 d2 1
H=− 2
+ mω02 x2
2m dx 2
so that Schrödinger’s equation becomes
∂Ψ 2 d2 Ψ 1
i = HΨ = − + mω02 x2 Ψ. (7.6)
∂t 2m dx2 2
We first look for solutions in which the variables x and t separate
Ψ(x, t) = A(t)ψ(x).
Substituting this into (7.6) and dividing the result by A(t)ψ(x) we find
1 dA 1
i = Hψ.
A dt ψ
Since the left hand side is only a function of t and the right hand side is only a function of
x both sides must equal a common constant. Calling this constant E (observe this constant
has the units of energy), we find
dA iE
= − A,
dt
Hψ = Eψ.
The first equation has solution

A(t) = e−iEt/
so that
Ψ(x, t) = e−iEt/ ψ(x). (7.7)
We now examine
Hψ = Eψ (7.8)
in detail. The first observation is that (7.8) is an eigenvalue problem in L2 (R). Thus the
eigenvalues of the operator H are interpreted as energies. It is convenient to introduce
dimensionless variables to simplify notationally the differential equation. Let

mω0 2E
ξ=x , ε= .
ω0
7 Recall 1 1
the potential energy for the harmonic oscillator is V (x) = 2
kx2 = 2
mω02 x2 .
Performing this change of variables, the Schrödinger equation Hψ = Eψ becomes

d2 ψ
− + (ξ 2 − ε)ψ = 0. (7.9)
dξ 2
We want solutions to (7.9) that are square integrable. It is convenient to also perform a
change in the dependent variable8
2
ψ(ξ) = e−ξ /2
v(ξ).
Then a straightforward calculation shows that v must satisfy the equation
d2 v dv
2
− 2ξ + (ε − 1)v = 0. (7.10)
dξ dξ
Observe that (7.10) is not a constant coefficient differential equation, so that the methods
we have developed do not apply to this equation.
7.2.2 Hermite polynomials

To find solutions of (7.10) we look for solutions that are of the form9
∞

v(ξ) = a0 + a1 ξ + a2 ξ 2 + a3 ξ 3 + · · · = ak ξ k . (7.11)
k=0
The idea is to substitute this into (7.10) and to find conditions that the coefficients ak must
satisfy. Since
∞

dv
= a1 + 2a2 ξ + 3a3 ξ 2 + · · · = kak ξ k−1
dξ
k=1
and
∞ ∞
d2 v
= 2a 2 + 6a 3 ξ + · · · = k(k − 1)a k ξ k−2
= (k + 1)(k + 2)ak+2 ξ k ,
dξ 2
k=2 k=0
we have
∞
d2 v dv
2
− 2ξ + (ε − 1)v = 2a2 + (ε − 1)a0 + {(k + 2)(k + 1)ak+2 + (ε − 1 − 2k)ak } ξ k .
dξ dξ
k=1
For a power series to be identically zero, each of the coefficients must be zero. Hence we
obtain10
(k + 2)(k + 1)ak+2 + (ε − 1 − 2k) ak = 0, k = 0, 1, 2, . . . (7.12)
Thus once a0 is specified, the coefficients a2 , a4 , a6 , . . . are determined from the above recur-
rence relation. Similarly, once a1 is specified the coefficients a3 , a5 , a7 , . . . are determined.
The recurrence relation (7.12) can be rewritten as
ak+2 2k − ε + 1
= , k = 0, 1, 2, . . . (7.13)
ak (k + 2)(k + 1)
8 This change of variables makes the recursion relation derived below simpler.
9 This is called the power series method.
10 Note that the k = 0 condition is 2a + (ε − 1)a = 0.
2 0
Our first observation from (7.13) is that

ak+2
lim =0
k→∞ ak
and so by the ratio test for power series, the radius of convergence of (7.11) is infinite. (This
is good since we want our functions ψ to be defined for all ξ.)
Now comes a crucial point. We have shown for any choices of a0 and a1 and for any
2
choice of the parameter (dimensionless energy) ε, that the function ψ(ξ) = e−ξ /2 v(ξ) solves
the differential equation (7.9) where v is given by (7.11) and the coefficients ak satisfy
(7.13). However, a basic requirement of the quantum mechanical formalism is that ψ(ξ)
is an element of the state space L2 (R); namely, it is square integrable. Thus the question
2
is whether e−ξ /2 v(ξ) is square integrable. We will show that we have square integrable
functions for only certain values of the energy ε; namely, we will find the quantization of
energy.
The ratio of the series coefficients, ak+1 /ak , in the function

∞
∞
αk
eαz = ak z k = zk
k!
k=0 k=0
is α/(k + 1) ∼ α/k as k → ∞. For the series (recall given a0 we can determine a2 , a4 , . . .)

∞
∞

v(ξ) = a2k ξ 2k = bk z k , bk = a2k , z = ξ 2 ,
k=0 k=0
the ratio of coefficients, bk+1 /bk , is (we use (7.13) to get the second equality)
bk+1 a2k+2 4k − ε + 1 1
= = ∼ , k → ∞.
bk a2k (2k + 2)(2k + 1) k
This suggests in comparing the series for v with the series for eαz , and it can be proved,11
that
2
v(ξ) ∼ eξ , ξ → ∞.
#∞ 2k+1
Similar remarks hold for the series k=0 a2k+1 ξ . This means our solution ψ(ξ) =
2 2
−ξ /2
v(ξ)e is not square integrable since it grows as eξ /2 . Hence ψ is not a valid state
in quantum mechanics. There is a way out of this: If the coefficients ak would vanish
identically from some point on, then the solution v(ξ) will be a polynomial and thus ψ will
be square integrable. From the recurrence relation (7.13) we see that this will happen if the
numerator vanishes for some value of k. That is, if
ε = 2n + 1
for some nonnegative integer n, then an+2 = an+4 = · · · = 0. It is traditional to choose a

normalization (which amounts to choices of a0 and a1 ) so that the coefficient of the highest
11 This asymptotic analysis can be made rigorous using the theory of irregular singular points.
power is 2n . With this normalization the polynomials are called Hermite polynomials and
are denoted by Hn (ξ). The first few polynomials are12
H0 (ξ) = 1,
H1 (ξ) = 2ξ,
H2 (ξ) = 4ξ 2 − 2,
H3 (ξ) = 8ξ 3 − 12ξ,
H4 (ξ) = 16ξ 4 − 48ξ 2 + 12,
H5 (ξ) = 32ξ 5 − 160ξ 3 + 120ξ,
H6 (ξ) = 64ξ 6 − 480ξ 4 + 720ξ 2 − 120.
Thus we have found solutions13

2
ψn (ξ) = Nn Hn (ξ)e−ξ /2
(7.14)
to (7.9); namely,
ω0
Hψn = (2n + 1)ψn , n = 0, 1, 2, . . .
2
We have solved an eigenvalue problem in the infinite dimensional space L2 (R). It is conve-
nient to choose the overall normalization constant Nn such that
ψn = 1, n = 0, 1, 2, . . .
That is, Nn is chosen so that

∞
2
Nn2 Hn (ξ)2 e−ξ dξ = 1. (7.15)
−∞
It can be shown that √ −1/2

Nn = π n! 2n .
7.2.3 Quantization of energy
The quantized energy levels are
1
En = ω0 εn = ω0 (n + 1/2), n = 0, 1, 2, . . . .
2
That is to say, the energy of the quantum oscillator cannot have arbitrary real values (as in
the case of the classical oscillator), but must be one of the discrete set of numbers
1 3 5
ω0 , ω0 , ω0 , . . .
2 2 2
12 One can compute a Hermite polynomial in Mathematica by the command HermiteH[n,x] where n is
a nonnegative integer.
13 Here N is an overall normalization constant which we choose below.
n
n=0
n=1
n=2
n=3
n=4
Figure 7.2: Harmonic Oscillator Wave Functions ψn (x) for n = 0, 1, 2, 3, 4.

1
The lowest energy, 2 ω0 , is called the ground state energy and has associated wave function
1 2
ψ0 (ξ) = e−ξ /2
.
π 1/4
Thus the ground state energy of the quantum harmonic oscillator is nonzero. In the classical
harmonic oscillator, we can have p = x = 0 which corresponds to E = 0.
7.2.4 Some properties of Hermite polynomials
Solution of recurrence relation
To obtain a more explicit formula for the Hermite polynomials we must solve the recurrence
relation (7.13). The polynomials Hn (x) are normalized so that the coefficient of the highest
power is 2n . This will determine a0 (when n is even) and a1 (when n is odd). We treat here
the case of n even and leave the case of n odd to the reader. First
an a2 a4 an
= ···
a0 a0 a2 an−2
The right hand side of this expression is determined from (7.13) and equals
2(n) 2(n − 2) 2(n − 4) 2(2)

···
1·2 3·4 5·6 (n − 1)n
This can be rewritten as
2n/2 n(n − 2)(n − 4) · · · 4 · 2
n!
This is the ratio an /a0 . Requiring that an = 2n gives
a0 = 2n/2 (n − 1)(n − 3)(n − 5) · · · 5 · 3 · 1
We now determine am —the coefficient of xm — (when n is even we can take m even too).
Proceeding in a similar manner we write
am a2 a4 am
= ···
a0 a0 a2 am−2
and again note the right hand side is determined from the recurrence relation (7.13); namely,
2(n) 2(n − 2) 2(n − 4) 2(n − m + 2)

(−1)m/2 ···
1·2 3·4 5·6 (m − 1) · m
Multiplying this by the value of a0 we get that am equals
2(n+m)/2
(−1)m/2 [n(n − 2)(n − 4) · · · (n − m + 2)] [(n − 1)(n − 3) · · · 5 · 3 · 1]
m!
The product of the two quantities in square brackets can be rewritten as
n!
(n − m − 1)(n − m − 3)(n − m − 5) · · · 5 · 3 · 1
(n − m)!
Now let m → n − m (so am is the coefficient of xn−m ) to find that am equals

n
(−1)m/2 2n−m/2 1 · 3 · 5 · · · (m − 1)
m
n
where m is the binomial coefficient. Since m is even and runs over 0, 2, 4, . . . n, we can let
m → 2m to get the final formula14
[n/2]
(−1)m
Hn (x) = n! (2x)n−2m . (7.16)
m=0
m!(n − 2m)!
This same formula holds for n odd if we interpret [n/2] = (n − 1)/2 when n is odd. From
(7.16) we can immediately derive the differentiation formula
dHn
= 2n Hn−1 (x). (7.17)
dx
Orthogonality properties
The harmonic oscillator Hamiltonian H is self-adjoint with distinct eigenvalues. Just as we

proved for matrices, it follows that that eigenfunctions ψn are orthogonal. The normalization
constant Nn is chosen so that they are orthonormal. That is if ψn are defined by (7.14),
then ∞
2
2
(ψn , ψm ) = Nn Hn (x) Hm (x) e−x dx = δm,n (7.18)
−∞
where Nn are defined in (7.15) and δm,n is the Kronecker delta function.15 The functions
ψn are called the harmonic oscillator wave functions.
From the orthogonality relations we can derive what is called the three-term recursion
relation; namely, we claim that
Hn+1 (x) − 2xHn (x) + 2nHn−1 (x) = 0. (7.19)
Since the highest power of Hn has coefficient 2n , we see that

Hn+1 (x) − 2xHn (x)
must be a polynomial of degree less than or equal to n. Using (7.16) we can see that the
highest power is the same as the highest power of 2n Hn−1 (x). Thus the left hand side of
(7.19) is a polynomial of degree less than or equal to n − 2. It can be written as the linear
combination
c0 H0 (x) + c1 H1 (x) + · · · + cn−2 Hn−2 (x).
2
We now multiply both sides of this resulting equation by Hk (x) e−x , 0 ≤ k ≤ n − 2, and
integrate over all of R. Using the orthogonality relation one concludes that ck = 0.16
14 We used the fact that
(2m − 1)!!/(2m)! = 1/(2m m!)
where (2m − 1)!! = (2m − 1)(2m − 3) · · · 5 · 3 · 1.
15 δ
m,n equals 1 if m = n and 0 otherwise.
16 Perhaps the only point that needs clarification is why

2
2xHn (x)Hk (x)e−x dx
R
For applications to the harmonic oscillator, it is convenient to find what (7.16) and (7.19)
imply for the oscillator wave functions ψn . It is an exercise to show that17

n n+1
x ψn (x) = ψn−1 (x) + ψn+1 (x) , (7.20)
2 2

dψn (x) n n+1
= ψn−1 (x) − ψn+1 (x) . (7.21)
dx 2 2
7.2.5 Completeness of the harmonic oscillator wave functions {ψn }n≥0
In finite-dimensional vector spaces, we understand the notion of a basis. In particular,

we’ve seen the importance of an orthonormal basis. In Hilbert spaces these concepts are
more subtle and a full treatment will not be given here. Here is what can be proved. For
any vector Ψ ∈ L2 (R) we can find coefficients an such that
∞

Ψ= an ψn . (7.22)
n=0
Since this is an infinite sum we must say in what sense this sum converges. If we define the
partial sums
n
Ψn = ak ψk ,
k=0
then we say Ψn → Ψ as n → ∞ if
lim Ψ − Ψn = 0.
n→∞
(Observe that Ψn , being a sum of a finite number of terms is well-defined.) Recall that the
norm · in L2 (R) is

2 2
Ψ − Ψn = |Ψ(x) − Ψn (x)| dx.
R
It is in this sense the series converges. Since ψn form an orthonormal sequence, the coeffi-
cients an are given simply by
an = (ψn , Ψ) .
Observe that since ψn form an orthonormal basis, the vector Ψ in (7.22) satisfies
∞

2 2
Ψ = |an | < ∞.
n=0
is zero for 0 ≤ k ≤ n − 2. Since 2xHk (x) is a polynomial of degree k + 1 ≤ n − 1, it too can be expanded
in terms of Hermite polynomials of degree less than or equal to n − 1; but these are all orthogonal to Hn .
Hence the expansion coefficients must be zero.
17 These formulas are valid for n = 0 if we define ψ
−1 (x) = 0.
7.3 Some properties of the harmonic oscillator

In quantum mechanics if O is an observable (mathematically, a self-adjoint operator on the
Hilbert space H), then the average (or expected) value of O in the state Ψ is
O = (OΨ, Ψ) .
For the quantum oscillator, the average position in the eigenstate ψn is

x = (xψn , ψn ) = xψn (x)2 dx = 0.
R
(The integral is zero since ψn2 is an even function of x so that xψn (x)2 is an odd function.)
The average of the square of the position in the eigenstate ψn is
$ 2% 2
x = x ψn , ψn .
This inner product (integral) can be evaluated by first using (7.20) twice to write x2 ψn as
a linear combination of the ψk ’s:

2 n n+1
x ψn = x ψn−1 + ψn+1
2 2

n n−1 n n+1 n+1 n+1
= ψn−2 + ψn + ψn + ψn+2
2 2 2 2 2 2
1 1 1
= n(n − 1) ψn−2 + (n + ) ψn + n + 1)(n + 2) ψn+2 .
2 2 2
The inner product can now be calculated using the orthonormality of the wave functions to
find18
$ 2% 1
x =n+ .
2
A very similar calculation with p̂ = −i dx

d
(but this time using (7.21)) gives19
p̂ = 0,
$ 2% 1
p̂ = (n + ) .
2
If we define

2
∆x = x2 − x ,

2
∆p = p̂2 − p̂ ,
18 This is the dimensionless result. Putting back in the dimensions, the average is

(n + 1/2)
mω0
.
19 Again these are the dimensionless results. Putting back the units the second average is
2 1
p̂ = (n + ) mω0 .
2
7.3. SOME PROPERTIES OF THE HARMONIC OSCILLATOR 133
then (in physical units) we have in state ψn

1
∆x ∆p = n + ≥ .
2 2
This is the Heisenberg Uncertainty Principle for the harmonic oscillator. The inequality
part of the statement can be shown to be valid under very general conditions.
7.3.1 Averages x̂(t) and p̂(t)
Let Ψ be any state of the system

∞

Ψ(x) = an ψn (x)
n=0
such that
n n+1
A := an ān−1 , B := an ān+1
2 2
n≥0 n≥0
are convergent sums. (Here, as throughout, the ψn are the harmonic oscillator wave func-
tions.) The time evolution of Ψ is then given by

Ψ(x, t) = an e−iEn t/ ψn (x) (7.23)
n≥0
which follows from the above discussion of separation of variables.

In the state Ψ(x, t) we are interested in computing the average values of x̂ and p̂. For
notationally convenience let’s define
xavg (t) = x̂ = (x̂Ψ(x, t), Ψ(x, t))
and
pavg (t) = p̂ = (p̂Ψ(x, t), Ψ(x, t)) .
Let
x0 := xavg (0) = (x̂Ψ(x, 0), Ψ(x, 0))
and
p0 := pavg (0) = (p̂Ψ(x, 0), Ψ(x, 0)) .
We first calculate x0 and p0 .

x0 = an ām (xψn , ψm )
m,n≥0
& '
n n+1
= an ām (ψn−1 , ψm ) + (ψn+1 , ψm )
2 2
m,n≥0
& '
n n+1
= an ām δn−1,m + δn+1,m
2 2
m,n≥0
= A+B
where we use the orthonormality of the functions ψn . Similarly,
p0 = −iA + iB.
We now calculate xavg (t) Now the state is (7.23). Proceeding as in the t = 0 case we see

xavg (t) = an ām e−i(En −Em )t/ (xψn , ψm ) .
m,n≥0
The calculation of the inner products (xψn , ψm ) was done in the t = 0 case. Noting that
En − En−1 = ω0 and En − En+1 = −ω0 ,
we see that
xavg (t) = e−iω0 t A + eiω0 t B. (7.24)
Similarly, we find
pavg (t) = −ie−iω0 t A + ieiω0 t B. (7.25)
Writing these averages in terms of sines and cosines and using the above expressions for x0
and p0 , we see that the average position and momentum in the state Ψ(x, t) evolve according
to20
p0
xavg (t) = x0 cos(ω0 t) + sin(ω0 t) (7.26)
mω0
pavg (t) = p0 cos(ω0 t) − mω0 x0 sin(ω0 t) (7.27)
One should now compare the time evolution of the quantum averages (7.26) and (7.27)
with the time evolution of the classical position and momentum (7.1) and (7.2). They are
identical. It is a special property of the quantum harmonic oscillator that the quantum
averages exactly follow the classical trajectories. More generally, one expects this to occur
only for states whose wave function remains localized in a region of space.
20 We restore the physical units in these last equations

7.4. THE HEISENBERG UNCERTAINTY PRINCIPLE 135
7.4 The Heisenberg Uncertainty Principle

The more precisely the position is determined, the less precisely the momentum
is known in this instant, and vice versa. Werner Heisenberg, 1927.
In §7.3 we proved the Heisenberg Uncertainty Principle for the special case of the harmonic
oscillator. Here we show this is a general feature of quantum mechanics.21 First we recall
some basic facts about complex vector spaces.
1. If Ψ and Φ are any two states in our Hilbert space of states H, we have an inner
product defined (Ψ, Φ) that satisfies the properties
(a) (Ψ, Φ) = (Φ, Ψ) where z denotes the complex conjugate of z.
(b) (c1 Ψ1 + c2 Ψ2 , Φ) = c1 (Ψ1 , Φ) + c2 (Ψ2 , Φ) for all states Ψ1 , Ψ2 and all complex
numbers c1 , c2 .
(c) The length or norm of the state Ψ is defined to be Ψ2 = (Ψ, Ψ) ≥ 0 with
Ψ = 0 if and only if Ψ = 0, the zero vector in H.
2. An operator A is called Hermitian (or self-adjoint) if
(AΨ, Φ) = (Ψ, AΦ)
for all states Ψ, Φ. In quantum mechanics observables are assumed to be Hermitian.
Note this makes the expected value of the observable A in state Ψ a real number
A := (AΨ, Ψ)
= (Ψ, AΨ)
= (AΨ, Ψ)
= A.
Sometimes one writes AΨ to denote the state in which the expected value is com-
puted.
3. Just as in linear algebra, we have the Cauchy-Schwarz inequality
2
|(Ψ, Φ)| ≤ Ψ2 Φ2 (7.28)
for all states Ψ, Φ ∈ H.
We now assume we have observables A and B that satisfy the commutation relation
AB − BA = i id (7.29)
where id is the identity operator and i is the imaginary number, i2 = −1. We showed earlier
that in units where = 1 the position and momentum operators satisfy such a commutation
relation. For a given state Ψ and observable A we define22
( )
2 2
∆A = (A − A) = A2 − A ≥ 0.
21 This is an optional section and is not covered in the class lectures.

22 In classical probability ∆A is called the standard deviation of A. The quantity ∆A is a measure of the
deviation of A from its expected value A .
We now prove that if observables A and B satisfy (7.29) then

1
∆A · ∆B ≥ . (7.30)
2
Proof: Let Ψ denote any normalized state vector, i.e. Ψ = 1. Apply (7.29) to Ψ to obtain
ABΨ − BAΨ = iΨ
Now take the inner product of each side with the state Ψ to obtain
(ABΨ, Ψ) − (BAΨ, Ψ) = i (Ψ, Ψ)
which simplifies to23
(BΨ, AΨ) − (AΨ, BΨ) = i (7.31)
Let t denote any real number, then by the Cauchy-Schwarz inequality (7.28)
2
|(Ψ, AΨ + itBΨ)| ≤ AΨ + itBΨ2 (7.32)
since Ψ = 1. Let’s simplify the left-hand side of (7.32)
(Ψ, AΨ + itBΨ) = (Ψ, AΨ) − it(Ψ, BΨ) = A − itB.
The absolute value squared of this is
A2 + t2 B2 .
We now examine the right-hand side of (7.32)
AΨ + itBΨ2 = (AΨ + itBΨ, AΨ + itBΨ)
= AΨ2 + it {(BΨ, AΨ) − (AΨ, BΨ)} + t2 BΨ2
= AΨ2 − t + t2 BΨ2 by use of (7.31). (7.33)
Thus the inequality (7.32) becomes
A2 + t2 B2 ≤ AΨ2 − t + t2 BΨ2 .
Using the fact that AΨ2 = (AΨ, AΨ) = (A2 Ψ, Ψ) = A2 (and similarly for B) and the
definition of ∆A (and similarly for ∆B), the above inequality can be rewritten as
2 2
t2 (∆B) − t + (∆A) ≥ 0.
This holds for all real t. The above is a quadratic polynomial in t that is always nonnegative.
This means that the discriminant of the quadratic polynomial must be nonpositive, i.e.
b2 − 4ac ≤ 0. That is,
2 2
1 − 4 (∆A) (∆B) ≤ 0
which implies that
1
∆A · ∆B ≥
2
which is what we want to prove.
When A is the position operator and B is the momentum operator we get the Heisenberg
Uncertainty Principle which states

∆x · ∆p ≥
2
where we have returned to physical units. The appearance of Planck’s constant in the
right hand side shows that sets the scale for quantum mechanical phenomena.
23 Note (ABΨ, Ψ) = (BΨ, AΨ) and (BAΨ, Ψ) = (AΨ, BΨ) since A and B are Hermitian. Also note on
the right hand side we used the fact that (Ψ, Ψ) = Ψ2 = 1.
7.5. COMPARISON OF THREE PROBLEMS 137
7.5 Comparison of three problems
Weighted String Vibrating String Quantum Harmonic Oscillator

Vector space RN Functions on [0, L] that vanish H = Square integrable
V at the endpoints 0 and L functions on R
L ∞
Inner product (f, g ) = f · g = N
i=1 fi gi (f, g) = 0
f (x)g(x) dx (f, g) = −∞
f (x)g(x) dx
1/2 1/2 1/2
∞
|f (x)|2 dx |f (x)|2 dx
L
Norm f = (f, f) f = 0
f = −∞
d2 2 2
Operator VN , N × N tridiagonal matrix L= dx2
H = − dx
d
2 + x

Symmetry VN f, g = f, VN g (Lf, g) = (f, Lg) (Hf, g) = (f, Hg)
Eigenvalues & VN fn = λnfn , n = 1, 2, . .. , N Lun = −kn2 un , n = 1, 2, 3, . . . Hψn = n ψn , n = 0, 1, 2, . . .

Eigenvectors λn = 2 1 − nπ
cos( N+1 ) kn = nπ n = 2n + 1
L 2
fn = (fn,j )1≤j≤N un (x) = L2 sin( nπ x) ψn (x) = Nn Hn (x)e−x /2
L
2 nπ
fn,j = N+1 sin( N+1 j) Hn (x) = Hermite polynomial
√ −1/2
Nn = [ π n! 2n ]
Orthonormal (fm , fn ) = δm,n (um , un ) = δm,n (ψm , ψn ) = δm,n

N ∞ ∞
Completeness u ∈ RN , u =
j=1 aj fj u ∈ V, u(x) = j=1 aj uj (x) Ψ ∈ H, Ψ(x) = j=0 aj ψj (x)

Coefficients aj = fj , u aj = (u, uj ) aj = (Ψ, ψj )
7.6 Exercises
#1.
Using (7.17) and (7.19), prove (7.20) and (7.21).
#2. Averages x̂4 and p̂4 :
For the state ψn , compute the averages x̂4 and p̂4 .
#3.
Prove (7.25). (The proof is similar to the proof of (7.24).)

#4.
Define the operators

1 d
a = √ x+
2 dx

1 d
a∗ = √ x−
2 dx
That is, if ψ = ψ(x), then

1 dψ
(aψ)(x) = √ xψ(x) +
2 dx
and similarly for a∗ . Using (7.20) and (7.21) show that for the harmonic oscillator wave
functions ψn
√
aψn = n ψn−1 , n ≥ 1, aψ0 = 0,
∗
√
a ψn = n + 1 ψn+1 , n = 0, 1, . . . ,
a∗ aψn = n ψn , n = 0, 1, . . . ,
(aa∗ − a∗ a) ψn = ψn , n = 0, 1, . . . .
Explain why this last equation implies the operator equation
[a, a∗ ] = id.
In quantum mechanics the operator a is called an annihilation operator and the operator
a∗ is called a creation operator. On the basis of this exercise, why do you think they have
these names?
#5. Hermite polynomials
We obtained the Hermite polynomials from the recurrence relation (7.13). Alternatively,
we have a generating formula for the Hermite polynomials. Starting with this (which many
books take as the definition of the Hermite polynomials), we may obtain the Schrödinger
equation.
1. Verify that the first three Hermite polynomials H0 (ξ), H1 (ξ) and H2 (ξ) are given using
the generating formula n
2 d −ξ2
Hn (ξ) = (−1)n eξ e (7.34)
dξ n
2. The generating function (7.34) can be used to give an alternative generating function
for Hermite polynomials. Show that
∞
2 zn
e−z +2zξ
= Hn (ξ) . (7.35)
n=0
n!
2 2
Hint: Let F (z) = e−z and consider the Taylor expansion of eξ F (z − ξ) about the
point z = 0.
7.6. EXERCISES 139
3. Derive (7.17) from (7.35). Now derive (7.19) using (7.34) and the newly derived (7.17).
4. Use (7.17) and (7.19) to show that the Hermite polynomials are solutions of the Her-
mite equation
d2 d
Hn (ξ) − 2ξ Hn (ξ) + 2nHn (ξ) = 0. (7.36)
dξ 2 dξ
5. We know that the Hermite polynomials satisfy
∞
2
2
Nn Hn (ξ)Hm (ξ)e−ξ dξ = δnm . (7.37)
−∞
2
Here by setting ψn (ξ) = Nn Hn (ξ)e−ξ /2 we see that ψn (ξ) are orthonormal in L2 (R).
Use (7.36) to obtain the differential equation that ψn (ξ) satisfy. You should obtain
2
(7.9) with ε = 2n + 1. This implies that ψn (ξ) = Nn Hn (ξ)e−ξ /2 is the eigenfunction
corresponding to the eigenvalue of the Hamiltonian operator.
Chapter 8
Heat Equation
Figure 8.1: Joseph Fourier (1768–1830) a French mathematician and physicist best known
for initiating the study of Fourier series with applications to the theory of oscillatory systems
and heat transfer.
141
142 CHAPTER 8. HEAT EQUATION
8.1 Introduction
The heat equation is the partial differential equation
∂u
− ∆u = 0 (8.1)
∂t
where ∆ is the Laplacian. In three-dimensions if (x, y, z) are Cartesian coordinates, the heat
equation reads
∂u ∂ 2 u ∂ 2 u ∂ 2 u
− 2 − 2 − 2 = 0,
∂t ∂x ∂y ∂z
whereas in one-dimension the heat equation reads
∂u ∂ 2 u
− 2 = 0. (8.2)
∂t ∂x
The heat equation has the following physical interpretation: the temperature u(x, t) at
position x and time t satisfies (8.1). This assumes that the medium is homogeneous and
the thermal diffusivity has been set equal to one. The heat equation is fundamental in
probability theory due to its connection with Brownian motion. In this chapter we show
how to solve the (8.2) subject to the initial condition u(x, 0) = f (x). The function f (x) can
be interpreted as the initial distribution of temperature.
8.2 Fourier transform

Suppose f = f (x), x ∈ R, is smooth (say continuous) and decays sufficiently fast at ±∞ so
that the integral
∞
f*(ξ) := f (x)e−2πixξ dx, ξ ∈ R, (8.3)
−∞
exists. The function f*(ξ) is called the Fourier transform of f . For example, if f (x) = e−a|x| ,
a > 0, then
2a
f*(ξ) = 2 .
a + 4π 2 ξ 2
The most important property of the Fourier transform is the Fourier inversion formula
which says given f*(ξ) we can find f (x) from the formula
∞
f (x) = f*(ξ)e2πixξ dξ, x ∈ R. (8.4)
−∞
Thus, for example, we have

∞
2a
e2πixξ dξ = e−a|x| .
−∞ a2 + 4π 2 ξ 2
8.3. SOLVING THE HEAT EQUATION BY USE OF THE FOURIER TRANSFORM143
8.3 Solving the heat equation by use of the Fourier

transform
We consider the 1D heat equation (8.2) for x ∈ R subject to the initial condition u(x, 0) =
f (x). Write u(x, t) in terms of its Fourier transform u *(ξ, t):
∞
u(x, t) = *(ξ, t)e2πixξ dξ.
u
−∞
Then ∞
∂u ∂ 2 u 2πixξ ∂*
u(ξ, t) 2 2
− 2 = e + 4π ξ u*(ξ, t) dξ
∂t ∂x −∞ ∂t
We want this to equal zero so this will certainly be the case if
∂*
u(ξ, t)
+ 4π 2 ξ 2 u
*(ξ, t) = 0.
∂t
We solve this last equation:
2 2
*(ξ, t) = e−4π
u ξ t
*(ξ, 0).
u
Thus we have ∞
2 2
u(x, t) = e2πixξ−4π ξ t
*(ξ, 0) dξ.
u (8.5)
−∞
Specializing to t = 0 and using the initial condition gives

∞
f (x) = u(x, 0) = e2πixξ u
*(ξ, 0) dξ.
−∞
By the Fourier inversion formula

∞
*(ξ, 0) =
u e−2πixξ f (x) dx = f*(ξ).
−∞
We now substitute this expression for u*(ξ, 0) into (8.5) to obtain

∞
e2πixξ−4π ξ t f*(ξ) dξ
2 2
u(x, t) =
−∞
∞ ∞
2πixξ−4π 2 ξ 2 t −2πix ξ
= e f (x )e dx dξ
−∞ −∞
∞ ∞
2 2
= e2πi(x−x )ξ−4π ξ t dξ f (x ) dx (8.6)
−∞ −∞
We now evaluate the integral

∞
2 2
I := e2πi(x−x )ξ−4π ξ t
dξ.
−∞
To do this we first evaluate the integral

∞
2
I1 = e−αξ dξ, α > 0.
−∞
We look at I12 and transform the resulting double integral to polar coordinates:
∞ ∞
2 2
2
I1 = e−αξ1 e−αξ2 dξ1 dξ2
−∞ −∞
∞ 2π
2
= e−αr r drdθ
0 0
∞
2
= 2π re−αr dr
0
π
=
α
Thus
π
. I1 =
α
We are now ready to evaluate I. Let y = x − x so the integral we wish to evaluate is
∞
2 2
e2πiyξ−4π ξ t dξ
−∞
We make the change of variables ξ = ξ + a and choose a so that the linear term in ξ
vanishes:

2πiyξ−4π 2 ξ 2 t = 2πiy(ξ +a)−4π 2 t(ξ 2 +2ξ a+a2 ) = −4π 2 tξ 2 + 2πiy − 8π 2 ta ξ +2πiya−4π 2 ta2
Thus we choose
2πiy iy
a= = .
8π 2 t 4πt
With this value of a we then have
y2
2πiyξ − 4π 2 ξ 2 t = −4π 2 tξ 2 −
4t
Thus
∞ ∞
2 2 2 2 2
e2πiyξ−4π ξ t
dξ = e−4π tξ −y /(4t) dξ
−∞ −∞
∞
−y 2 /(4t) 2 2 1 2
= e e−4π tξ dξ = √ e−y /(4t)
−∞ 4πt
We summarize what we have found in the following theorem:
Theorem. Let f (x) be continuous and decaying sufficiently fast at infinity so that its
Fourier transform exists. Then the solution to the one-dimensional heat equation (8.2)
satisfying the initial condition u(x, 0) = f (x) is
∞
u(x, t) = K(x, y; t)f (y) dy (8.7)
−∞
where
1 2
K(x, y; t) = √ e−(x−y) /(4t) . (8.8)
4πt
8.3. SOLVING THE HEAT EQUATION BY USE OF THE FOURIER TRANSFORM145
Heat kernel K0,y;t for fixed values of t
0.8
0.6
t110
t12
0.4 t1
t2
0.2
y
5 5
Figure 8.2: The heat kernel K(x, y; t), (8.8), for x = 0 as a function of y for various fixed
values of t. Note that for small t the functions is concentrated near x = 0 but for large t, K
is spread out over a wide range of values.
8.3.1 Comments on (8.7) and (8.8)
1. For all x ∈ R and all t > 0 we have

∞
K(x, y; t) dy = 1.
−∞
2. Let’s take the initial condition to be (see Figure 8.3)



 0 if x < −1,

1+x if −1 ≤ x < 0,
f (x) = (8.9)
 1−x
 if 0 ≤ x < 1,

0 if x ≥ 1.
In Figure 8.4 we plot the solution u(x, t) as a function of x for various values of t.
1.0
0.8
0.6
0.4
0.2
x
2 1 1 2
Figure 8.3: The initial temperature distribution f (x) as defined in (8.9). Note that the
distribution is not smooth at zero.
8.4 Heat equation on the semi-infinite rod
We now consider the heat equation but with the spatial variable x restricted to positive
axis—the semi-infinite rod. Thus we wish to solve (8.2) with an initial heat distribution
f (x), x > 0; and we further assume, that the end of the rod, x = 0, is fixed to zero
temperature for all t > 0.
For the infinite rod we found a solution of the form

∞
u(x, t) = K(x, x ; t)f (x ) dx .
−∞
For the semi-infinite rod we look for a solution of the form

∞
u(x, t) = Ksemi (x, x ; t)f (x ) dx .
0
We want Ksemi (x, x ; t) to satisfy the heat equation as a function of (x, t) but we also want
Ksemi(0, x ; t) = 0 for all t.
This last condition will guarantee that u(0, t) = 0 for all t > 0. Observe that
Ksemi(x, x ; t) = K(x, x ; t) − K(x, −x ; t)

1 −(x−x )2 /4t 1 −(x+x )2 /4t
= √ e − e
4πt 4πt
satisfies these conditions. What we have done is to place a second “heat pulse” at −x . This
technique is called the method of images.
8.4. HEAT EQUATION ON THE SEMI-INFINITE ROD 147
1.0
0.8
t0
0.6
t110
t12
0.4
t1
0.2
x
4 2 2 4
Figure 8.4: The solution u(x, t) with initial condition (8.9). Note that the solution is
“smoothed” by the action of the heat operator.
0.5
x
4 2 2 4
0.5
Figure 8.5: The kernel Ksemi (x, x ; t) as a function of x for the values x = 3 and t = 1/10.
0.4
0.3
0.2
0.1
x
2 4 6 8
Figure 8.6: Plot of u(x, t) as a function of x for three values of t = 1/10, 1/2, 1 with initial
heat distribution f (x) = e−x . As t increases the heat distribution becomes flatter and less
peaked.
8.5 Heat equation on the circle

Suppose now that the spatial variable x is restricted to a circle; that is, we want a solution
to (8.2) that is periodic in x
u(x + 1, t) = u(x, t)
(we’ve taken the period to equal 1 for simplicity). We are given an initial distribution f (x),
0 ≤ x ≤ 1. Can we use the method of images to find a
Kcircle(x, x ; t),
so that the solution we seek is of the form
1
u(x, t) = Kcircle (x, x ; t)f (x ) dx
0
and satisfies u(x, 0) = f (x) and u(x + 1, t) = u(x)?

We take the kernel K(x, x ; t) and place image pulses at all the integers:
∞

Kcircle (x, x ; t) = K(x, x + n; t)
n=−∞
∞
1 2
= √ e−(x−x −n) /(4πt)
4πt n=−∞
8.6. EXERCISES 149
8.6 Exercises
#1.
Consider the partial differential equation for u = u(x, t):
∂u ∂2u ∂u
= + (1 − a) , t > 0, −∞ < x < ∞, (8.10)
∂t ∂x2 ∂x
with initial condition u(x, 0) = f (x). When a = 1 this reduces to the heat equation.
Following the method (Fourier transforms) used to solve the heat equation in §8.3, find the
analogue of equations (8.7) and (8.8).
Chapter 9
Laplace Transform
Figure 9.1: Pierre-Simon Laplace (1749–1827) was a French mathematician known for his
contributions to celestial mechanics, probability theory and analysis. Both the Laplace
equation and the Laplace transform are named after Pierre-Simon Laplace.
151
152 CHAPTER 9. LAPLACE TRANSFORM
9.1 Matrix version of the method of Laplace transforms

for solving constant coefficient DEs
The Laplace transform of a function f (t) is
∞
F (s) = L(f )(s) = e−ts f (t) dt (9.1)
0
for (s) sufficiently large and positive. For the Laplace transform to make sense the function
f cannot grow faster that an exponential near infinity. Thus, for example, the Laplace
2
transform of ex is not defined.
We extend (9.1) to vector-valued functions f (t),
 
f1 (t)
 f2 (t) 
 
f (t) =  .  (9.2)
 .. 
fn (t)
by  ∞ −ts 
0∞ e−ts f1 (t) dt
 
 0 e f2 (t) dt 
F (s) = L(f )(s) =  .. . (9.3)
 . 
∞ −ts
0
e fn (t) dt
Integration by parts shows that
df
L( )(s) = sL(f )(s) − f (0). (9.4)
dt
We now explain how matrix Laplace transforms are used to solve the matrix ODE
dx
= Ax + f (t) (9.5)
dt
where A is a constant coefficient n × n matrix, f (t) is a vector-valued function of the
independent variable t (“forcing term”) with initial condition
x(0) = x0 . (9.6)
First, we take the Laplace transform of both sides of (9.5). From (9.4) we see that the
Laplace transform of the LHS of (9.5) is
dx
L( ) = sL(x) − x0 .
dt
9.1. MATRIX VERSION 153
The Laplace transform of the RHS of (9.5) is
L(Ax + f ) = L(Ax) + L(f )

= AL(x) + F (s)
where we set F (s) = L(f )(s) and we used the fact that A is independent of t to conclude1
L(Ax) = AL(x). (9.7)
Thus the Laplace transform of (9.5) is
sL(x) − x0 = AL(x) + F,
or
(sIn − A)L(x) = x0 + F (s) (9.8)
where In is the n×n identity matrix. Equation (9.8) is a linear system of algebraic equations
for L(x). We now proceed to solve (9.8). This can be done once we know that (sIn − A) is
invertible. Recall that a matrix is invertible if and only if the determinant of the matrix is
nonzero. The determinant of the matrix in question is
p(s) := det(sIn − A), (9.9)
which is the characteristic polynomial of the matrix A. We know that the zeros of p(s) are
the eigenvalues of A. If s is larger than the absolute value of the largest eigenvalue of A; in
symbols,
s > max|λi |, (9.10)
then p(s) cannot vanish and hence (sIn − A)−1 exists. We assume s satisfies this condition.
Then multiplying both sides of (9.8) by (sIn − A)−1 results in
L(x)(s) = (sIn − A)−1 x0 + (sIn − A)−1 F (s). (9.11)
Equation (9.11) is the basic result in the application of Laplace transforms to the solution
of constant coefficient differential equations with an inhomogeneous forcing term. Equation
(9.11) will be a quick way to solve initial value problems once we learn efficient methods
to (i) compute (sIn − A)−1 , (ii) compute the Laplace transform of various forcing terms
F (s) = L(f )(s), and (iii) find the inverse Laplace transform. Step (i) is easier if one
uses software packages such as MatLab . Steps (ii) and(iii) are made easier by the use of
extensive Laplace transform tables or symbolic integration packages such as Mathematica.
It should be noted that many of the DE techniques one learns in engineering courses can
be described as efficient methods to do these three steps for examples that are of interest to
engineers.
We now give two examples that apply (9.11).
1 You are asked to prove (9.7) in an exercise.
9.1.1 Example 1
Consider the scalar ODE

d2 y dy
2
+b + cy = f (t) (9.12)
dt dt
where b and c are constants. We first rewrite this as a system

x1 y(t)
x(t) = =
x2 y (t)
so that
dx 0 1 0
= x+ .
dt −c −b f (t)
Then
s −1
sI2 − A = ,
c s+b
and
−1 1 s+b 1
(sI2 − A) = 2 .
s + bs + c −c s
Observe that the characteristic polynomial
p(s) = det(sI2 − A) = s2 + bs + c
appears in the denominator of the matrix elements of (sI2 − A)−1 . (This factor in Laplace
transforms should be familiar from the scalar treatment—here we see it is the characteristic
polynomial of A.) By (9.11)

1 (s + b)y(0) + y (0) F (s) 1
L(x)(s) = 2 + 2
s + bs + c −cy(0) + sy (0) s + bs + c s
where F (s) = L(f )(s). This implies that the Laplace transform of y(t) is given by
(s + b)y(0) + y (0) F (s)

L(y)(s) = + 2 . (9.13)
s2 + bs + c s + bs + c
This derivation of (9.13) may be compared with the derivation of equation (16) on page 302
of Boyce and DiPrima [4] (in our example a = 1).
9.1.2 Example 2
We consider the system (9.5) for the special case of n = 3 with f (t) = 0 and A given by
 
1 0 −1
A = 1 2 1. (9.14)
1 −1 −1
The characteristic polynomial of (9.14) is
p(s) = s3 − 2s2 + s − 2 = (s2 + 1)(s − 2) (9.15)

9.1. MATRIX VERSION 155
and so the matrix A has eigenvalues ±i and 2. A rather long linear algebra computation
shows that  2 
s −s−1 1 −s + 2
1 
(sI3 − A)−1 = s+2 s2 s−2  . (9.16)
p(s)
s−3 −s + 1 s2 − 3s + 2
If one writes a partial fraction decomposition of each of the matrix elements appearing in
(9.16) and collects together terms with like denominators, then (9.16) can be written as
 
1/5 1/5 0
1 
(sI3 − A)−1 = 4/5 4/5 0
s−2
−1/5 −1/5 0
 
(3 + 4s)/5 −(2 + s)/5 −1
1 
+ 2 −(3 + 4s)/5 (2 + s)/5 1 . (9.17)
s +1
(7 + s)/5 (−3 + s)/5 −1 + s
We now apply (9.17) to solve (9.5) with the above A and f = 0 for the case of initial
conditions  
1
x0 = −2 . (9.18)
1
We find
     
−1/5 6/5 2/5
1 −4/5 + s −6/5 + 1 −2/5 . (9.19)
L(x)(s) = (sI3 − A)−1 x0 =
s−2 s2 + 1 s2 + 1
1/5 4/5 8/5
To find x(t) from (9.19) we use Table 6.2.1 on page 300 of Boyce and DiPrima [4]; in
particular, entries 2, 5, and 6. Thus
     
−1/5 6/5 2/5
x(t) = e2t −4/5 + cos t −6/5 + sin t −2/5 .
1/5 4/5 8/5
One can also use Mathematica to compute the inverse Laplace transforms. To do so use
the command InverseLaplaceTransform. For example if one inputs
InverseLaplaceTransform[1/(s-2),s,t] then the output is e2t .
We give now a second derivation of (9.19) using the eigenvectors of A. As noted above,
the eigenvalues of A are λ1 = 2, λ2 = i, and λ3 = −i. If we denote by φj an eigenvector
associated to eigenvalue λj (j = 1, 2, 3), then a routine linear algebra computation gives the
following possible choices for the φj :
     
−1 (1 + i)/2 (1 − i)/2
φ1 = −4 , φ2 = −(1 + i)/2 , φ3 = (−1 + i)/2 .
1 1 1
Now for any eigenvector φ corresponding to eigenvalue λ of a matrix A we have
(sIn − A)−1 φ = (s − λ)−1 φ.
To use this observation we first write
x0 = c1 φ1 + c2 φ2 + c3 φ3 .
A computation shows that
c1 = 1/5, c2 = 2/5 − 4i/5, and c3 = 2/5 + 4i/5.
Thus
1 2 − 4i 2 + 4i
(sI3 − A)−1 x0 = (s − 2)−1 φ1 + (s − i)−1 φ2 + (s + i)−1 φ3 .
5 5 5
Combining the last two terms gives (9.19).
9.2 Structure of (sIn − A)−1 when A is diagonalizable

In this section we assume that the matrix A is diagonalizable; that is, we assume a set of
linearly independent eigenvectors of A form a basis. Recall the following two theorems from
linear algebra: (1) If the n × n matrix A has n distinct eigenvalues, then A is diagonalizable;
and (2) If the matrix A is symmetric (hermitian if the entries are complex), then A is
diagonalizable.
Since A is assumed to be diagonalizable, there exists a nonsingular matrix P such that
A = P DP −1
where D is  
λ1 0 ... 0
0 λ2 ... 0 
 
D=. .. .. .. 
 .. . . . 
0 0 . . . λn
and each eigenvalue λi of A appears as many times as the (algebraic) multiplicity of λi .
Thus
sIn − A = sIn − P DP −1
= P (sIn − D)P −1 ,
so that
−1
(sIn − A)−1 = P (sIn − D)P −1 )
= P (sIn − D)−1 P −1 .
Since P and P −1 are independent of s, the s dependence of (sIn − A)−1 resides in the
diagonal matrix (sIn − D)−1 . This tells us that the partial fraction decomposition of the
matrix (sIn − A)−1 is of the form

n
1
(sIn − A)−1 = Pj
j=1
s − λj
where
Pj = P Ej P −1
9.3. EXERCISES 157
and Ej is the diagonal matrix with all zeros on the main diagonal except for 1 at the (j, j)th
entry. This follows from the fact that

n
1
(sIn − D)−1 = Ej
j=1
s − λj
Note that Pj have the property that

Pj2 = Pj .
Such matrices are called projection operators.
In general, it follows from Cramer’s method of computing the inverse of a matrix, that the
general structure of (sIn − A)−1 will be 1/p(s) times a matrix whose entries are polynomials
of at most degree n − 1 in s. When an eigenvalue, say λ1 , is degenerate and of (algebraic)
multiplicity m1 , then the characteristic polynomial will have a factor (s − λ1 )m1 . We have
seen that if the matrix is diagonalizable, upon a partial fraction decomposition only a single
power of (s − λ1 ) will appear in the denominator of the partial fraction decompostion.
Finally, we conclude by mentioning that when the matrix A is not diagonalizable, then this
is reflected in the partial fraction decomposition of (sIn − A)−1 in that some powers of
(s − λj ) occur to a higher degree than 1.
9.3 Exercises
#1.
Use the Laplace transform to find the solution of the initial value problem
     
1 −1 0 0 0
dx 
= 0 −1 1  x +  12  , x(0) =  0  .
dt
0 1 −1 0 0
#2.
Let A be a n × n matrix whose entries are real numbers and x ∈ Rn . Prove that
L(Ax) = AL(x)
where L denotes the Laplace transform.
#3.
Let Ej denote the diagonal n × n matrix with all zeros on the main diagonal except for 1
at the (j, j) entry.
• Prove that Ej2 = Ej .
• Show that if P is any invertible n × n matrix, then Pj2 = Pj where Pj := P Ej P −1 .

#4.
It is a fact that you will learn in an advanced linear algebra course, that if a 2 × 2 matrix
A is not diagonalizable, then there exists a nonsingular matrix P such that
A = P B P −1
where
λ 1
B=
0 λ
for some constant λ.
• Show that λ must be an eigenvalue of A with algebraic multiplicity 2.

• Find an eigenvector of A (in terms of the matrix P ), and show that A has no other
eigenvectors (except, of course, scalar multiples of the vector you have already found).
• Show that
1 1 1
(sI2 − A)−1 = P E1 P −1 + P E2 P −1 + P N P −1
s−λ s−λ (s − λ)2
where
0 1
N= .
0 0
• Relate what is said here to the remarks in the footnote in Exercise 5.5.2.
Bibliography
[1] V. I. Arnold, Mathematical Methods of Classical Mechanics, 2nd ed., Springer-

Verlag, N.Y. 1989.
[2] V. I. Arnold, Ordinary Differential Equations, Second printing of the 1992 edition,
Springer-Verlag, Berlin, 2006.
[3] G. Birkhoff and G-C Rota, Ordinary Differential Equations 4th ed., John Wiley &
Sons, 1989.
[4] W. E. Boyce and R. C. DiPrima, Elementary Differential Equations and Boundary

Value Problems, 7th–9th eds., John Wiley & Sons.
[5] G. Chen, P. J. Morris and J. Zhou, Visualization of special eigenmode shapes of a
vibrating elliptical membrane, SIAM Review 36 (1994), 453–469.
[6] M. Levi, Classical Mechanics with Calculus of Variations and Optimal Control: An
Intuitive Introduction, American Mathematical Society, 2014.
[7] J. R. Norris, Markov Chains, Cambridge University Press, 1997.
[8] B. Temple and C. A. Tracy, From Newton to Einstein, American Math. Monthly
99 (1992), 507–512.
159
Index
Acceleration due to gravity, 5 Golden ratio φ, 64

Acoustic phonons, 117
Airy differential equation, 47 Half-life, 35
Airy functions, 47 Harmonic oscillator, 26, 47, 83, 104
Angular momentum, 38 Forced oscillations, 54
Annihilation and creation operators, 138 Quantum, 123
Arnold, Vladimir, v, 67, 89 Harmonic oscillator wave functions, 130
Heat equation, 141
Bessel functions, 108 Semi-infinite rod, 146
Binet’s equation, 40 Heisenberg Uncertainty Principle, 135
Heisenberg Uncertainty Relation, 124
Catalan numbers, 65 Helmholtz equation, 105, 116
Classical harmonic oscillator, 26, 47 Hermite polynomials, 125, 127, 129, 130, 138
Forced oscillations, 54 Generating functions, 138
Friction, 37 Hilbert space, 122
Commutator, 123 Hooke’s Law, 23
Complex exponential function, 50 Forced oscillations, 54
Addition formula, 51 Resonance, 56
Euler’s formula, 51 Frequency, 26
Friction, 37
Diagonalizable matrix, 71 Nonlinear correction, 37
Difference equations, 62 Potential energy, 24
Constant coefficients, 62 Spring constant, 23
Catalan numbers, 65 Hooke, Robert, 23
Fibonacci numbers, 63
Differential equation, definition, 15 Integrating factor, 17
Dirac, Paul, 121
Kinetic energy, 22, 118, 124
Energy conservation, 20, 27 Kolmogorov forward equations, 85
Energy level curves, 33
Euler’s equation, 59 Laplace transform, 152
Euler’s formula, 43, 58 Application to radioactive decay, 85
Euler, Leonhard, 43 Laplacian, 106
Lax pairs, 114
Fibonacci numbers, 63 Linear first order ODE, 15
Fibonacci, Leonardo, 61 integrating factor, 17
Fourier inversion formula, 142 Loan payments, 18
Fourier transform, 142 Logistic growth model, 36
Fourier, Joseph, 141
Fundamental matrix, 81 Markov processes and Q-matrix, 75
160
INDEX 161
Mathematica, 9, 66, 76, 85, 123, 155 Toda chain, 114

MatLab, 8, 30, 53
Symbolic Math Toolbox, 8 Variation of parameters, 79
Matrix exponential, 68 Vibrating membrane, 105
Applied to matrix ODE, 70 Circular domain, 107
Diagonalizable matrix, 71 Elliptical domain, 109
Nilpotent matrix, 83 Rectangular domain, 107, 116
Matrix ODE, 70 Vibrating string, 101
Alternating masses-spring problem, 117 Energy, 118
Inhomogeneous, 79 von Helmholtz, Hermann, 89
Markov processes, 75
Radioactive decay, 77, 85 Wave equation, 101
Method of images, 146 Helmholtz equation, 105
Mortgage payment problem, 18 Relation to harmonic oscillators, 104
Separation of variables, 102
Newton’s equations, 5, 22, 90, 112, 122 Vibrating membrane, 105
Central field, 38 Vibrating string, 101
Inverse square law, 39 Weighted string, 90
Newton’s principle of determinacy, 5 Circle, 110
Nonautonomous linear equations, 81 Friction, 115
Normal modes, 94 Normal modes, 94
Wronskian, 45, 60
Optical phonons, 117
Pendulum, 5
Coupled, 111
Period, 27
Planck’s constant, 122
Population dynamics, 36
Potential energy, 20, 118, 124
Quantum harmonic oscillator, 123

Annihilation and creation operators, 138
Hermite polynomials, 127
Quantization of energy, 127
Quantum mechanics, 122
State space, 122
Radioactive decay, 35, 77, 85

Reduction of order, 48
Rigid body motion, 40
Roots of unity, 58
Schrödinger equation, 122

Schrödinger, Erwin, 121
Second order linear ODE, 44
Bessel functions, 108
Constant coefficients, 49
Vector space of solutions, 44
Wronskian, 45

Lectures On Differential Equations: Department of Mathematics University of California Davis, CA 95616

Uploaded by

Copyright:

Available Formats

Lectures On Differential Equations: Department of Mathematics University of California Davis, CA 95616

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Lectures On Differential Equations: Department of Mathematics University of California Davis, CA 95616

Uploaded by

Copyright:

Available Formats

Lectures on Diﬀerential Equations1

2 First Order Equations & Conservative Systems 15

3 Second Order Linear Equations 43

5 Matrix Diﬀerential Equations 67

7 Quantum Harmonic Oscillator 121

8 Heat Equation 141

9 Laplace Transform 151

Mathematics is a part of physics. Physics is an experimental science, a part of

Newton’s fundamental discovery, the one which he considered necessary to keep

See the You Tube video http://youtu.be/MpzaCCbX-z4.

1.1 What is a diﬀerential equation?

A diﬀerential equation is an equation between speciﬁed derivative on an unknown

is said to be of the ﬁrst order.

A diﬀerential equation is a mathematical equation that relates some function of

can be viewed as a unifying principle behind diverse phenomena. As an example,

y(x) = c eλx , c = constant,

satisﬁes this equation since

• Now integrate both sides

• Solve for y = y(x)

3. An example of a second order ODE is

1.2 Diﬀerential equation for the pendulum

V. I. Arnold, Mathematical Methods of Classical Mechanics[1]

= length of the rod measured, say, in meters,

equations become two equations

Fx = −T sin θ, Fy = T cos θ − mg.

(T is the magnitude of the vector T .)

Substitute (1.7) in (1.4) and (1.8) in (1.5) to obtain

−T sin θ = m cos θ θ̈ − m sin θ (θ̇)2 , (1.9)

1.3 Introduction to MatLab, Mathematica and Maple

and then press the enter key it responds with

Multiplication is denoted by * and division by / . Thus, for example, to compute

>> help taylor

 ∞2. One can now use f in other Mathematica

Mathematica returns the answer 1/4.

4. Suppose we want to create the 10 × 10 matrix

In Mathematica the command is

Mathematica returns the answer

The result is the graph shown in Figure 1.3.

#1. MatLab and/or Mathematica Exercises

>> help format

>> format long e

Ans.: 3.697296376497268 × 10197 .

3. Use MatLab or Mathematica to ﬁnd the determinant, eigenvalues and eigenvectors

>> A=[1 -1 2 0; sqrt(2) 1 0 -2;0 1 sqrt(2) -1; 1 2 2 0]

To ﬁnd the determinant

and to ﬁnd the eigenvalues

If you also want the eigenvectors you enter

and then re-enter the command for the determinant

and this time MatLab returns

#2. Inverted pendulum

First Order Equations &

• y(x0 ) = y0 and y0 is arbitrary. That is, there is a one-parameter family of solutions;

A simple generalization of (2.1) is to replace the right-hand side by a function that

In calculus one typically deﬁnes ex as the limit

follows almost immediately?

2.1.2 Method of integrating factors

Deﬁning µ(x) by (2.6), the diﬀerential equation (2.5) is transformed to

and solving this for y gives our ﬁnal formula

(T is the magnitude of the vector T .)

∞2. One can now use f in other Mathematica

Thus (2.16) becomes3

Recalling (2.24), we conclude

a · (b × c) = b · (c × a) = c · (a × b).