Introduction To Differential Equations: Jeffrey R. Chasnov
Introduction To Differential Equations: Jeffrey R. Chasnov
Introduction To Differential Equations: Jeffrey R. Chasnov
Equations
Jeffrey R. Chasnov
Copyright ○
c 2009–2014 by Jeffrey Robert Chasnov
This work is licensed under the Creative Commons Attribution 3.0 Hong Kong License.
To view a copy of this license, visit http://creativecommons.org/licenses/by/3.0/hk/
or send a letter to Creative Commons, 171 Second Street, Suite 300, San Francisco,
California, 94105, USA.
Preface
What follows are my lecture notes for a first course in differential equations,
taught at the Hong Kong University of Science and Technology. Included in
these notes are links to short tutorial videos posted on YouTube.
Much of the material of Chapters 2-6 and 8 has been adapted from the
widely used textbook “Elementary differential equations and boundary value
problems” by Boyce & DiPrima (John Wiley & Sons, Inc., Seventh Edition,
○2001).
c Many of the examples presented in these notes may be found in this
book. The material of Chapter 7 is adapted from the textbook “Nonlinear
dynamics and chaos” by Steven H. Strogatz (Perseus Publishing, ○1994).
c
All web surfers are welcome to download these notes, watch the YouTube
videos, and to use the notes and videos freely for teaching and learning. An
associated free review book with links to YouTube videos is also available from
the ebook publisher bookboon.com. I welcome any comments, suggestions or
corrections sent by email to [email protected]. Links to my website, these
lecture notes, my YouTube page, and the free ebook from bookboon.com are
given below.
Homepage:
http://www.math.ust.hk/~machas
YouTube:
https://www.youtube.com/user/jchasnov
Lecture notes:
http://www.math.ust.hk/~machas/differential-equations.pdf
Bookboon:
http://bookboon.com/en/differential-equations-with-youtube-examples-ebook
iii
Contents
1 Introduction to odes 11
1.1 The simplest type of differential equation . . . . . . . . . . . . . 11
2 First-order odes 13
2.1 The Euler method . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2 Separable equations . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3 Linear equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.4 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.4.1 Compound interest . . . . . . . . . . . . . . . . . . . . . . 20
2.4.2 Chemical reactions . . . . . . . . . . . . . . . . . . . . . . 21
2.4.3 Terminal velocity . . . . . . . . . . . . . . . . . . . . . . . 23
2.4.4 Escape velocity . . . . . . . . . . . . . . . . . . . . . . . . 24
2.4.5 RC circuit . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4.6 The logistic equation . . . . . . . . . . . . . . . . . . . . . 27
v
vi CONTENTS
5 Series solutions 61
5.1 Ordinary points . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.2 Regular singular points: Cauchy-Euler equations . . . . . . . . . 65
5.2.1 Real, distinct roots . . . . . . . . . . . . . . . . . . . . . . 67
5.2.2 Complex conjugate roots . . . . . . . . . . . . . . . . . . 67
5.2.3 Repeated roots . . . . . . . . . . . . . . . . . . . . . . . . 67
6 Systems of equations 69
6.1 Determinants and the eigenvalue problem . . . . . . . . . . . . . 69
6.2 Coupled first-order equations . . . . . . . . . . . . . . . . . . . . 71
6.2.1 Two distinct real eigenvalues . . . . . . . . . . . . . . . . 71
6.2.2 Complex conjugate eigenvalues . . . . . . . . . . . . . . . 75
6.2.3 Repeated eigenvalues with one eigenvector . . . . . . . . . 77
6.3 Normal modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
A short mathematical
review
1
2 CHAPTER 0. A SHORT MATHEMATICAL REVIEW
The exponential function exp (𝑥) = 𝑒𝑥 and natural logarithm ln 𝑥 are inverse
functions satisfying
𝑒ln 𝑥 = 𝑥, ln 𝑒𝑥 = 𝑥.
The usual rules of exponents apply:
ln (𝑥𝑦) = ln 𝑥 + ln 𝑦, ln (𝑥/𝑦) = ln 𝑥 − ln 𝑦, ln 𝑥𝑝 = 𝑝 ln 𝑥.
𝑓 (𝑥 + ℎ) − 𝑓 (𝑥)
𝑓 ′ (𝑥) = lim . (1)
ℎ→0 ℎ
(𝑓 + 𝑔)′ = 𝑓 ′ + 𝑔 ′ .
(𝑓 − 𝑔)′ = 𝑓 ′ − 𝑔 ′ .
(𝑓 𝑔)′ = 𝑓 ′ 𝑔 + 𝑓 𝑔 ′ ,
and should be memorized as “the derivative of the first times the second plus
the first times the derivative of the second.”
and should be memorized as “the derivative of the top times the bottom minus
the top times the derivative of the bottom over the bottom squared.”
0.5. DIFFERENTIATING ELEMENTARY FUNCTIONS 3
and should be memorized as “the derivative of the outside times the derivative
of the inside.”
𝑑 𝑝
𝑥 = 𝑝𝑥𝑝−1 .
𝑑𝑥
We thus say that “the derivative of sine is cosine,” and “the derivative of cosine
is minus sine.” Notice that the second derivatives satisfy
where 𝑁 = (𝑏 − 𝑎)/ℎ is the number of terms in the sum. The symbols on the
left-hand-side of (2) are read as “the integral from 𝑎 to 𝑏 of f of x dee x.” The
4 CHAPTER 0. A SHORT MATHEMATICAL REVIEW
Riemann Sum definition is extended to all values of 𝑎 and 𝑏 and for all values
of 𝑓 (𝑥) (positive and negative). Accordingly,
∫︁ 𝑎 ∫︁ 𝑏 ∫︁ 𝑏 ∫︁ 𝑏
𝑓 (𝑥)𝑑𝑥 = − 𝑓 (𝑥)𝑑𝑥 and (−𝑓 (𝑥))𝑑𝑥 = − 𝑓 (𝑥)𝑑𝑥.
𝑏 𝑎 𝑎 𝑎
which states (when 𝑓 (𝑥) > 0) that the total area equals the sum of its parts.
𝑥
∫︀ 𝑥+ℎ ∫︀ 𝑥
𝑓 (𝑠)𝑑𝑠 − 𝑓 (𝑠)𝑑𝑠
∫︁
𝑑 𝑎 𝑎
𝑓 (𝑠)𝑑𝑠 = lim
𝑑𝑥 𝑎 ℎ→0 ℎ
∫︀ 𝑥+ℎ
𝑓 (𝑠)𝑑𝑠𝑥
= lim
ℎ→0 ℎ
ℎ𝑓 (𝑥)
= lim
ℎ→0 ℎ
= 𝑓 (𝑥).
This result is called the fundamental theorem of calculus, and provides a con-
nection between differentiation and integration.
The fundamental theorem teaches us how to integrate functions. Let 𝐹 (𝑥)
be a function such that 𝐹 ′ (𝑥) = 𝑓 (𝑥). We say that 𝐹 (𝑥) is an antiderivative of
𝑓 (𝑥). Then from the fundamental theorem and the fact that the derivative of a
constant equals zero,
∫︁ 𝑥
𝐹 (𝑥) = 𝑓 (𝑠)𝑑𝑠 + 𝑐.
𝑎
∫︀ 𝑏
Now, 𝐹 (𝑎) = 𝑐 and 𝐹 (𝑏) = 𝑎 𝑓 (𝑠)𝑑𝑠 + 𝐹 (𝑎). Therefore, the fundamental
theorem shows us how to integrate a function 𝑓 (𝑥) provided we can find its
antiderivative:
∫︁ 𝑏
𝑓 (𝑠)𝑑𝑠 = 𝐹 (𝑏) − 𝐹 (𝑎). (3)
𝑎
The last expression has an interesting structure. All the values of 𝐹 (𝑥) eval-
uated at the points lying between the endpoints 𝑎 and 𝑏 cancel each other in
consecutive terms. Only the value −𝐹 (𝑎) survives when 𝑛 = 1, and the value
+𝐹 (𝑏) when 𝑛 = 𝑁 , yielding again (3).
𝑥3 7𝑥2
∫︁
(𝑥2 + 7𝑥 + 2)𝑑𝑥 = + + 2𝑥 + 𝑐,
3 2
and ∫︁
(5 cos 𝑥 + sin 𝑥)𝑑𝑥 = 5 sin 𝑥 − cos 𝑥 + 𝑐.
0.10 Substitution
More complicated functions can be integrated using the chain rule. Since
𝑑 (︀
𝑓 𝑔(𝑥) = 𝑓 ′ 𝑔(𝑥) · 𝑔 ′ (𝑥),
)︀ (︀ )︀
𝑑𝑥
we have ∫︁
𝑓 ′ 𝑔(𝑥) · 𝑔 ′ (𝑥)𝑑𝑥 = 𝑓 𝑔(𝑥) + 𝑐.
(︀ )︀ (︀ )︀
= 𝑓 (𝑦) + 𝑐
(︀ )︀
= 𝑓 𝑔(𝑥) + 𝑐.
𝑧 = 𝑥 + 𝑖𝑦,
where 𝑥 and 𝑦 are real numbers. We call 𝑥 the real part of 𝑧 and 𝑦 the imaginary
part and write
𝑥 = Re 𝑧, 𝑦 = Im 𝑧.
Two complex numbers are equal if and only if their real and imaginary parts
are equal.
The complex conjugate of 𝑧 = 𝑥 + 𝑖𝑦, denoted as 𝑧¯, is defined as
𝑧¯ = 𝑥 − 𝑖𝑦.
𝑧 𝑧¯ = (𝑥 + 𝑖𝑦)(𝑥 − 𝑖𝑦)
= 𝑥2 − 𝑖2 𝑦 2
= 𝑥2 + 𝑦 2 ;
|𝑧| = (𝑧 𝑧¯)1/2
√︀
= 𝑥2 + 𝑦 2 .
We can add, subtract, multiply and divide complex numbers to get new
complex numbers. With 𝑧 = 𝑥 + 𝑖𝑦 and 𝑤 = 𝑠 + 𝑖𝑡, and 𝑥, 𝑦, 𝑠, 𝑡 real numbers,
we have
𝑧𝑤 = (𝑥 + 𝑖𝑦)(𝑠 + 𝑖𝑡)
= (𝑥𝑠 − 𝑦𝑡) + 𝑖(𝑥𝑡 + 𝑦𝑠);
𝑧 𝑧𝑤¯
=
𝑤 𝑤𝑤¯
(𝑥 + 𝑖𝑦)(𝑠 − 𝑖𝑡)
=
𝑠2 + 𝑡2
(𝑥𝑠 + 𝑦𝑡) (𝑦𝑠 − 𝑥𝑡)
= 2 +𝑖 2 .
𝑠 + 𝑡2 𝑠 + 𝑡2
0.13. COMPLEX NUMBERS 9
Furthermore,
√︀
|𝑧𝑤| = (𝑥𝑠 − 𝑦𝑡)2 + (𝑥𝑡 + 𝑦𝑠)2
√︀
= (𝑥2 + 𝑦 2 )(𝑠2 + 𝑡2 )
= |𝑧||𝑤|;
and
Similarly
⃒𝑧⃒ |𝑧| 𝑧 𝑧¯
⃒ ⃒= , ( )= .
⃒ ⃒
𝑤 |𝑤| 𝑤 𝑤¯
Also, 𝑧 + 𝑤 = 𝑧 + 𝑤. However, |𝑧 + 𝑤| ≤ |𝑧| + |𝑤|, a theorem known as the
triangle inequality.
It is especially interesting and useful to consider the exponential function of
an imaginary argument. Using the Taylor series expansion of an exponential
function, we have
(𝑖𝜃)2 (𝑖𝜃)3 (𝑖𝜃)4 (𝑖𝜃)5
𝑒𝑖𝜃 = 1 + (𝑖𝜃) + + + + ...
2! )︂3! (︂ 4! 3 5!5
2 4
(︂ )︂
𝜃 𝜃 𝜃 𝜃
= 1− + − ... + 𝑖 𝜃 − + + ...
2! 4! 3! 5!
= cos 𝜃 + 𝑖 sin 𝜃.
Therefore, we have
cos 𝜃 = Re 𝑒𝑖𝜃 , sin 𝜃 = Im 𝑒𝑖𝜃 .
Since cos 𝜋 = −1 and sin 𝜋 = 0, we derive the celebrated Euler’s identity
𝑒𝑖𝜋 + 1 = 0,
that links five fundamental numbers, 0, 1, 𝑖, 𝑒 and 𝜋, using three basic mathe-
matical operations, addition, multiplication and exponentiation, only once.
Using the even property cos (−𝜃) = cos 𝜃 and the odd property sin (−𝜃) =
− sin 𝜃, we also have
𝑒−𝑖𝜃 = cos 𝜃 − 𝑖 sin 𝜃;
and the identities for 𝑒𝑖𝜃 and 𝑒−𝑖𝜃 results in the frequently used expressions,
𝑒𝑖𝜃 + 𝑒−𝑖𝜃 𝑒𝑖𝜃 − 𝑒−𝑖𝜃
cos 𝜃 = , sin 𝜃 = .
2 2𝑖
The complex number 𝑧 can be represented in the complex plane with Re 𝑧
as the 𝑥-axis and Im 𝑧 as the 𝑦-axis. This leads to the polar representation of
𝑧 = 𝑥 + 𝑖𝑦:
𝑧 = 𝑟𝑒𝑖𝜃 ,
where 𝑟 = |𝑧| and tan 𝜃 = 𝑦/𝑥. We define arg 𝑧 = 𝜃. Note that 𝜃 is not unique,
though it is conventional to choose the value such that −𝜋 < 𝜃 ≤ 𝜋, and 𝜃 = 0
when 𝑟 = 0.
10 CHAPTER 0. A SHORT MATHEMATICAL REVIEW
Useful trigonometric relations can be derived using 𝑒𝑖𝜃 and properties of the
exponential function. The addition law can be derived from
We have
yielding
cos(𝑥 + 𝑦) = cos 𝑥 cos 𝑦 − sin 𝑥 sin 𝑦, sin(𝑥 + 𝑦) = sin 𝑥 cos 𝑦 + cos 𝑥 sin 𝑦.
Therefore,
cos 2𝜃 = cos2 𝜃 − sin2 𝜃, sin 2𝜃 = 2 cos 𝜃 sin 𝜃.
With a little more manipulation using cos 𝜃 + sin2 𝜃 = 1, we can derive
2
1 + cos 2𝜃 1 − cos 2𝜃
cos2 𝜃 = , sin2 𝜃 = ,
2 2
which are useful formulas for determining
∫︁ ∫︁
1 1
cos2 𝜃 𝑑𝜃 = (2𝜃 + sin 2𝜃) + 𝑐, sin2 𝜃 𝑑𝜃 = (2𝜃 − sin 2𝜃) + 𝑐,
4 4
from which follows
∫︁ 2𝜋 ∫︁ 2𝜋
2
sin 𝜃 𝑑𝜃 = cos2 𝜃 𝑑𝜃 = 𝜋.
0 0
Chapter 1
Introduction to odes
𝑑𝑛 𝑥
= 𝐺(𝑡),
𝑑𝑡𝑛
where the derivative of 𝑥 = 𝑥(𝑡) can be of any order, and the right-hand-side
may depend only on the independent variable 𝑡. As an example, consider a mass
falling under the influence of constant gravity, such as approximately found on
the Earth’s surface. Newton’s law, 𝐹 = 𝑚𝑎, results in the equation
𝑑2 𝑥
𝑚 = −𝑚𝑔,
𝑑𝑡2
where 𝑥 is the height of the object above the ground, 𝑚 is the mass of the
object, and 𝑔 = 9.8 meter/sec2 is the constant gravitational acceleration. As
Galileo suggested, the mass cancels from the equation, and
𝑑2 𝑥
= −𝑔.
𝑑𝑡2
Here, the right-hand-side of the ode is a constant. The first integration, obtained
by antidifferentiation, yields
𝑑𝑥
= 𝐴 − 𝑔𝑡,
𝑑𝑡
11
12 CHAPTER 1. INTRODUCTION TO ODES
with 𝐴 the first constant of integration; and the second integration yields
1
𝑥 = 𝐵 + 𝐴𝑡 − 𝑔𝑡2 ,
2
with 𝐵 the second constant of integration. The two constants of integration 𝐴
and 𝐵 can then be determined from the initial conditions. If we know that the
initial height of the mass is 𝑥0 , and the initial velocity is 𝑣0 , then the initial
conditions are
𝑑𝑥
𝑥(0) = 𝑥0 , (0) = 𝑣0 .
𝑑𝑡
Substitution of these initial conditions into the equations for 𝑑𝑥/𝑑𝑡 and 𝑥 allows
us to solve for 𝐴 and 𝐵. The unique solution that satisfies both the ode and
the initial conditions is given by
1
𝑥(𝑡) = 𝑥0 + 𝑣0 𝑡 − 𝑔𝑡2 . (1.1)
2
For example, suppose we drop a ball off the top of a 50 meter building. How
long will it take the ball to hit the ground? This question requires solution of
(1.1) for the time 𝑇 it takes for 𝑥(𝑇 ) = 0, given 𝑥0 = 50 meter and 𝑣0 = 0.
Solving for 𝑇 ,
√︂
2𝑥0
𝑇 =
𝑔
√︂
2 · 50
= sec
9.8
≈ 3.2sec.
Chapter 2
First-order differential
equations
Reference: Boyce and DiPrima, Chapter 2
The general first-order differential equation for the function 𝑦 = 𝑦(𝑥) is written
as
𝑑𝑦
= 𝑓 (𝑥, 𝑦), (2.1)
𝑑𝑥
where 𝑓 (𝑥, 𝑦) can be any function of the independent variable 𝑥 and the depen-
dent variable 𝑦. We first show how to determine a numerical solution of this
equation, and then learn techniques for solving analytically some special forms
of (2.1), namely, separable and linear first-order equations.
13
14 CHAPTER 2. FIRST-ORDER ODES
𝑑𝑦
𝑔(𝑦) = 𝑓 (𝑥), 𝑦(𝑥0 ) = 𝑦0 , (2.2)
𝑑𝑥
where the function 𝑔(𝑦) is independent of 𝑥 and 𝑓 (𝑥) is independent of 𝑦. Inte-
gration from 𝑥0 to 𝑥 results in
∫︁ 𝑥 ∫︁ 𝑥
′
𝑔(𝑦(𝑥))𝑦 (𝑥)𝑑𝑥 = 𝑓 (𝑥)𝑑𝑥.
𝑥0 𝑥0
and since 𝑢 is a dummy variable of integration, we can write this in the equivalent
form ∫︁ 𝑦 ∫︁ 𝑥
𝑔(𝑦)𝑑𝑦 = 𝑓 (𝑥)𝑑𝑥. (2.3)
𝑦0 𝑥0
A simpler procedure that also yields (2.3) is to treat 𝑑𝑦/𝑑𝑥 in (2.2) like a fraction.
Multiplying (2.2) by 𝑑𝑥 results in
𝑔(𝑦)𝑑𝑦 = 𝑓 (𝑥)𝑑𝑥,
which is a separated equation with all the dependent variables on the left-side,
and all the independent variables on the right-side. Equation (2.3) then results
directly upon integration.
𝑑𝑦
Example: Solve 𝑑𝑥 + 21 𝑦 = 32 , with 𝑦(0) = 2.
We first manipulate the differential equation to the form
𝑑𝑦 1
= (3 − 𝑦), (2.4)
𝑑𝑥 2
and then treat 𝑑𝑦/𝑑𝑥 as if it was a fraction to separate variables:
𝑑𝑦 1
= 𝑑𝑥.
3−𝑦 2
We integrate the right-side from the initial condition 𝑥 = 0 to 𝑥 and the left-side
from the initial condition 𝑦(0) = 2 to 𝑦. Accordingly,
∫︁ 𝑦 ∫︁ 𝑥
𝑑𝑦 1
= 𝑑𝑥. (2.5)
2 3−𝑦 2 0
2.2. SEPARABLE EQUATIONS 15
3
y
0
0 1 2 3 4 5 6 7
x
𝑑𝑦
Figure 2.1: Solution of the following ode: 𝑑𝑥 + 21 𝑦 = 32 .
The integrals in (2.5) need to be done. Note that 𝑦(𝑥) < 3 for finite 𝑥 or the
integral on the left-side diverges. Therefore, 3 − 𝑦 > 0 and integration yields
]︀𝑦 1 ]︀𝑥
− ln (3 − 𝑦) 2 = 𝑥 0 ,
2
1
ln (3 − 𝑦) = − 𝑥,
2
− 21 𝑥
3−𝑦 =𝑒 ,
1
𝑦 = 3 − 𝑒− 2 𝑥 .
Since this is our first nontrivial analytical solution, it is prudent to check our
result. We do this by differentiating our solution:
𝑑𝑦 1 1
= 𝑒− 2 𝑥
𝑑𝑥 2
1
= (3 − 𝑦);
2
𝑑𝑦
Example: Solve 𝑑𝑥 + 12 𝑦 = 32 , with 𝑦(0) = 4.
This is the identical differential equation as before, but with different initial
conditions. We will jump directly to the integration step:
∫︁ 𝑦
1 𝑥
∫︁
𝑑𝑦
= 𝑑𝑥.
4 3−𝑦 2 0
16 CHAPTER 2. FIRST-ORDER ODES
The solution curves for a range of initial conditions is presented in Fig. 2.1.
All solutions have a horizontal asymptote at 𝑦 = 3 at which 𝑑𝑦/𝑑𝑥 = 0. For
𝑦(0) = 𝑦0 , the general solution can be shown to be 𝑦(𝑥) = 3+(𝑦0 −3) exp(−𝑥/2).
𝑑𝑦
Example: Solve 𝑑𝑥 = 23+2𝑦
cos 2𝑥
, with 𝑦(0) = −1. (i) For what values of
𝑥 > 0 does the solution exist? (ii) For what value of 𝑥 > 0 is 𝑦(𝑥)
maximum?
Notice that the solution of the ode may not exist when 𝑦 = −3/2, since 𝑑𝑦/𝑑𝑥 →
∞. We separate variables and integrate from initial conditions:
(3 + 2𝑦)𝑑𝑦 = 2 cos 2𝑥 𝑑𝑥
∫︁ 𝑦 ∫︁ 𝑥
(3 + 2𝑦)𝑑𝑦 = 2 cos 2𝑥 𝑑𝑥
−1 0
]︀𝑦 ]︀𝑥
3𝑦 + 𝑦 2 −1 = sin 2𝑥 0
𝑦 2 + 3𝑦 + 2 − sin 2𝑥 = 0
1 √
𝑦± = [−3 ± 1 + 4 sin 2𝑥].
2
Solving the quadratic equation for 𝑦 has introduced a spurious solution that
does not satisfy the initial conditions. We test:
{︂
1 -1;
𝑦± (0) = [−3 ± 1] =
2 -2.
Only the + root satisfies the initial condition, so that the unique solution to the
ode and initial condition is
1 √
𝑦 = [−3 + 1 + 4 sin 2𝑥]. (2.6)
2
To determine (i) the values of 𝑥 > 0 for which the solution exists, we require
1 + 4 sin 2𝑥 ≥ 0,
or
1
sin 2𝑥 ≥ − . (2.7)
4
Notice that at 𝑥 = 0, we have sin 2𝑥 = 0; at 𝑥 = 𝜋/4, we have sin 2𝑥 = 1;
at 𝑥 = 𝜋/2, we have sin 2𝑥 = 0; and at 𝑥 = 3𝜋/4, we have sin 2𝑥 = −1 We
therefore need to determine the value of 𝑥 such that sin 2𝑥 = −1/4, with 𝑥 in
the range 𝜋/2 < 𝑥 < 3𝜋/4. The solution to the ode will then exist for all 𝑥
between zero and this value.
2.3. LINEAR EQUATIONS 17
−0.2
−0.4
−0.6
−0.8
y
−1
−1.2
−1.4
−1.6
0 0.5 1 1.5
x
Figure 2.2: Solution of the following ode: (3 + 2𝑦)𝑦 ′ = 2 cos 2𝑥, 𝑦(0) = −1.
To solve sin 2𝑥 = −1/4 for 𝑥 in the interval 𝜋/2 < 𝑥 < 3𝜋/4, one needs to
recall the definition of arcsin, or sin−1 , as found on a typical scientific calculator.
The inverse of the function
is denoted by arcsin. The first solution with 𝑥 > 0 of the equation sin 2𝑥 = −1/4
places 2𝑥 in the interval (𝜋, 3𝜋/2), so to invert this equation using the arcsine
we need to apply the identity sin (𝜋 − 𝑥) = sin 𝑥, and rewrite sin 2𝑥 = −1/4 as
sin (𝜋 − 2𝑥) = −1/4. The solution of this equation may then be found by taking
the arcsine, and is
𝜋 − 2𝑥 = arcsin (−1/4),
or (︂ )︂
1 1
𝑥= 𝜋 + arcsin .
2 4
Therefore the solution exists for 0 ≤ 𝑥 ≤ (𝜋 + arcsin (1/4)) /2 = 1.6971 . . . ,
where we have used a calculator value (computing in radians) to find arcsin(0.25) =
0.2527 . . . . At the value (𝑥, 𝑦) = (1.6971 . . . , −3/2), the solution curve ends and
𝑑𝑦/𝑑𝑥 becomes infinite.
To determine (ii) the value of 𝑥 at which 𝑦 = 𝑦(𝑥) is maximum, we examine
(2.6) directly. The value of 𝑦 will be maximum when sin 2𝑥 takes its maximum
value over the interval where the solution exists. This will be when 2𝑥 = 𝜋/2,
or 𝑥 = 𝜋/4 = 0.7854 . . . .
The graph of 𝑦 = 𝑦(𝑥) is shown in Fig. 2.2.
The first-order linear differential equation (linear in 𝑦 and its derivative) can be
written in the form
𝑑𝑦
+ 𝑝(𝑥)𝑦 = 𝑔(𝑥), (2.8)
𝑑𝑥
with the initial condition 𝑦(𝑥0 ) = 𝑦0 . Linear first-order equations can be inte-
grated using an integrating factor 𝜇(𝑥). We multiply (2.8) by 𝜇(𝑥),
[︂ ]︂
𝑑𝑦
𝜇(𝑥) + 𝑝(𝑥)𝑦 = 𝜇(𝑥)𝑔(𝑥), (2.9)
𝑑𝑥
and try to determine 𝜇(𝑥) so that
[︂ ]︂
𝑑𝑦 𝑑
𝜇(𝑥) + 𝑝(𝑥)𝑦 = [𝜇(𝑥)𝑦]. (2.10)
𝑑𝑥 𝑑𝑥
Equation (2.9) then becomes
𝑑
[𝜇(𝑥)𝑦] = 𝜇(𝑥)𝑔(𝑥). (2.11)
𝑑𝑥
Equation (2.11) is easily integrated using 𝜇(𝑥0 ) = 𝜇0 and 𝑦(𝑥0 ) = 𝑦0 :
∫︁ 𝑥
𝜇(𝑥)𝑦 − 𝜇0 𝑦0 = 𝜇(𝑥)𝑔(𝑥)𝑑𝑥,
𝑥0
or (︂ ∫︁ 𝑥 )︂
1
𝑦= 𝜇0 𝑦0 + 𝜇(𝑥)𝑔(𝑥)𝑑𝑥 . (2.12)
𝜇(𝑥) 𝑥0
It remains to determine 𝜇(𝑥) from (2.10). Differentiating and expanding (2.10)
yields
𝑑𝑦 𝑑𝜇 𝑑𝑦
𝜇 + 𝑝𝜇𝑦 = 𝑦+𝜇 ;
𝑑𝑥 𝑑𝑥 𝑑𝑥
and upon simplifying,
𝑑𝜇
= 𝑝𝜇. (2.13)
𝑑𝑥
Equation (2.13) is separable and can be integrated:
∫︁ 𝜇 ∫︁ 𝑥
𝑑𝜇
= 𝑝(𝑥)𝑑𝑥,
𝜇0 𝜇 𝑥0
∫︁ 𝑥
𝜇
ln = 𝑝(𝑥)𝑑𝑥,
𝜇0 𝑥0
(︂∫︁ 𝑥 )︂
𝜇(𝑥) = 𝜇0 exp 𝑝(𝑥)𝑑𝑥 .
𝑥0
𝑑𝑦
Example: Solve 𝑑𝑥 + 2𝑦 = 𝑒−𝑥 , with 𝑦(0) = 3/4.
Note that this equation is not separable. With 𝑝(𝑥) = 2 and 𝑔(𝑥) = 𝑒−𝑥 , we
have
(︂∫︁ 𝑥 )︂
𝜇(𝑥) = exp 2𝑑𝑥
0
= 𝑒2𝑥 ,
and
(︂ ∫︁ 𝑥 )︂
3
𝑦 = 𝑒−2𝑥 + 𝑒2𝑥 𝑒−𝑥 𝑑𝑥
4 0
(︂ ∫︁ 𝑥 )︂
3
= 𝑒−2𝑥 + 𝑒𝑥 𝑑𝑥
4 0
(︂ )︂
3
= 𝑒−2𝑥 + (𝑒𝑥 − 1)
4
(︂ )︂
1
= 𝑒−2𝑥 𝑒𝑥 −
4
(︂ )︂
1
= 𝑒−𝑥 1 − 𝑒−𝑥 .
4
𝑑𝑦
Example: Solve 𝑑𝑥 − 2𝑥𝑦 = 𝑥, with 𝑦(0) = 0.
This equation is separable, and we solve it in two ways. First, using an inte-
grating factor with 𝑝(𝑥) = −2𝑥 and 𝑔(𝑥) = 𝑥:
(︂ ∫︁ 𝑥 )︂
𝜇(𝑥) = exp −2 𝑥𝑑𝑥
0
−𝑥2
=𝑒 ,
and ∫︁ 𝑥
2 2
𝑦 = 𝑒𝑥 𝑥𝑒−𝑥 𝑑𝑥.
0
1 𝑥2 (︁ 2
)︁
𝑦= 𝑒 1 − 𝑒−𝑥
2
1 (︁ 𝑥2 )︁
= 𝑒 −1 .
2
20 CHAPTER 2. FIRST-ORDER ODES
𝑑𝑦
− 2𝑥𝑦 = 𝑥,
𝑑𝑥
𝑑𝑦
= 𝑥(1 + 2𝑦),
∫︁ 𝑦𝑑𝑥 ∫︁ 𝑥
𝑑𝑦
= 𝑥𝑑𝑥,
0 1 + 2𝑦 0
1 1
ln (1 + 2𝑦) = 𝑥2 ,
2 2
2
1 + 2𝑦 = 𝑒𝑥 ,
1 (︁ 𝑥2 )︁
𝑦= 𝑒 −1 .
2
The results from the two different solution methods are the same, and the choice
of method is a personal preference.
2.4 Applications
2.4.1 Compound interest
view tutorial
The equation for the growth of an investment with continuous compounding
of interest is a first-order differential equation. Let 𝑆(𝑡) be the value of the
investment at time 𝑡, and let 𝑟 be the annual interest rate compounded after
every time interval Δ𝑡. We can also include deposits (or withdrawals). Let 𝑘 be
the annual deposit amount, and suppose that an installment is deposited after
every time interval Δ𝑡. The value of the investment at the time 𝑡 + Δ𝑡 is then
given by
𝑆(𝑡 + Δ𝑡) = 𝑆(𝑡) + (𝑟Δ𝑡)𝑆(𝑡) + 𝑘Δ𝑡, (2.14)
where at the end of the time interval Δ𝑡, 𝑟Δ𝑡𝑆(𝑡) is the amount of interest
credited and 𝑘Δ𝑡 is the amount of money deposited (𝑘 > 0) or withdrawn
(𝑘 < 0). As a numerical example, if the account held $10,000 at time 𝑡, and
𝑟 = 6% per year and 𝑘 = $12,000 per year, say, and the compounding and
deposit period is Δ𝑡 = 1 month = 1/12 year, then the interest awarded after
one month is 𝑟Δ𝑡𝑆 = (0.06/12) × $10,000 = $50, and the amount deposited is
𝑘Δ𝑡 = $1000.
Rearranging the terms of (2.14) to exhibit what will soon become a deriva-
tive, we have
𝑆(𝑡 + Δ𝑡) − 𝑆(𝑡)
= 𝑟𝑆(𝑡) + 𝑘.
Δ𝑡
The equation for continuous compounding of interest and continuous deposits
is obtained by taking the limit Δ𝑡 → 0. The resulting differential equation is
𝑑𝑆
= 𝑟𝑆 + 𝑘, (2.15)
𝑑𝑡
which can solved with the initial condition 𝑆(0) = 𝑆0 , where 𝑆0 is the initial
capital. We can solve either by separating variables or by using an integrating
2.4. APPLICATIONS 21
To have saved approximately one million US$ at retirement, the worker would
need to save about HK$50,000 per year over his/her working life. Note that the
amount saved over the worker’s life is approximately 40 × $50,000 = $2,000,000,
while the amount earned on the investment (at the assumed 6% real return) is
approximately $8,000,000 − $2,000,000 = $6,000,000. The amount earned from
the investment is about 3× the amount saved, even with the modest real return
of 6%. Sound investment planning is well worth the effort.
Similarly, the law of mass action enables us to write equations for the time-
derivatives of the reactant concentrations 𝐴 and 𝐵:
𝑑𝐴 𝑑𝐵
= −𝑘𝐴𝐵, = −𝑘𝐴𝐵. (2.18)
𝑑𝑡 𝑑𝑡
The ode given by (2.17) can be solved analytically using conservation laws. We
assume that 𝐴0 and 𝐵0 are the initial concentrations of the reactants, and that
no product is initially present. From (2.17) and (2.18),
𝑑
(𝐴 + 𝐶) = 0 =⇒ 𝐴 + 𝐶 = 𝐴0 ,
𝑑𝑡
𝑑
(𝐵 + 𝐶) = 0 =⇒ 𝐵 + 𝐶 = 𝐵0 .
𝑑𝑡
Using these conservation laws, (2.17) becomes
𝑑𝐶
= 𝑘(𝐴0 − 𝐶)(𝐵0 − 𝐶), 𝐶(0) = 0,
𝑑𝑡
which is a nonlinear equation that may be integrated by separating variables.
Separating and integrating, we obtain
∫︁ 𝐶 ∫︁ 𝑡
𝑑𝐶
=𝑘 𝑑𝑡
0 (𝐴0 − 𝐶)(𝐵0 − 𝐶) 0
= 𝑘𝑡. (2.19)
The remaining integral can be done using the method of partial fractions. We
write
1 𝑎 𝑏
= + . (2.20)
(𝐴0 − 𝐶)(𝐵0 − 𝐶) 𝐴0 − 𝐶 𝐵0 − 𝐶
The cover-up method is the simplest method to determine the unknown coeffi-
cients 𝑎 and 𝑏. To determine 𝑎, we multiply both sides of (2.20) by 𝐴0 − 𝐶 and
set 𝐶 = 𝐴0 to find
1
𝑎= .
𝐵 0 − 𝐴0
Similarly, to determine 𝑏, we multiply both sides of (2.20) by 𝐵0 − 𝐶 and set
𝐶 = 𝐵0 to find
1
𝑏= .
𝐴0 − 𝐵 0
Therefore,
(︂ )︂
1 1 1 1
= − ,
(𝐴0 − 𝐶)(𝐵0 − 𝐶) 𝐵 0 − 𝐴0 𝐴0 − 𝐶 𝐵0 − 𝐶
and the remaining integral of (2.19) becomes (using 𝐶 < 𝐴0 , 𝐵0 )
∫︁ 𝐶 (︃∫︁ )︃
𝐶 ∫︁ 𝐶
𝑑𝐶 1 𝑑𝐶 𝑑𝐶
= −
0 (𝐴0 − 𝐶)(𝐵0 − 𝐶) 𝐵0 − 𝐴0 0 𝐴0 − 𝐶 0 𝐵0 − 𝐶
(︂ (︂ )︂ (︂ )︂)︂
1 𝐴0 − 𝐶 𝐵0 − 𝐶
= − ln + ln
𝐵0 − 𝐴0 𝐴0 𝐵0
(︂ )︂
1 𝐴0 (𝐵0 − 𝐶)
= ln .
𝐵0 − 𝐴0 𝐵0 (𝐴0 − 𝐶)
2.4. APPLICATIONS 23
As one would expect, the reaction stops after one of the reactants is depleted;
and the final concentration of product is equal to the initial concentration of
the depleted reactant.
𝑑𝑣
𝑚 = −𝑚𝑔 − 𝑘𝑣. (2.21)
𝑑𝑡
The terminal velocity 𝑣∞ of the mass is defined as the asymptotic velocity after
air resistance balances the gravitational force. When the mass is at terminal
velocity, 𝑑𝑣/𝑑𝑡 = 0 so that
𝑚𝑔
𝑣∞ = − . (2.22)
𝑘
decays to zero.
As an example, a skydiver of mass 𝑚 = 100 kg with his parachute closed
may have a terminal velocity of 200 km/hr. With
2 2
𝑔 = (9.8 m/s )(10−3 km/m)(60 s/min)2 (60 min/hr)2 = 127, 008 km/hr ,
one obtains from (2.22), 𝑘 = 63, 504 kg/hr. One-half of the terminal velocity
for free-fall (100 km/hr) is therefore attained when (1 − 𝑒−𝑘𝑡/𝑚 ) = 1/2, or 𝑡 =
𝑚 ln 2/𝑘 ≈ 4 sec. Approximately 95% of the terminal velocity (190 km/hr ) is
attained after 17 sec.
𝐺𝑀
𝑔= ,
𝑅2
and 𝑔 ≈ 9.8 m/s2 . Newton’s law 𝐹 = 𝑚𝑎 for the mass 𝑚 is thus given by
𝑑2 𝑥 𝐺𝑀
2
=−
𝑑𝑡 (𝑅 + 𝑥)2
𝑔
=− , (2.23)
(1 + 𝑥/𝑅)2
i
(a)
(b)
+
C
ߝ
_
2.4.5 RC circuit
view tutorial
Consider a resister 𝑅 and a capacitor 𝐶 connected in series as shown in Fig. 2.3.
A battery providing an electromotive force, or emf ℰ, connects to this circuit
by a switch. Initially, there is no charge on the capacitor. When the switch
is thrown to a, the battery connects and the capacitor charges. When the
switch is thrown to b, the battery disconnects and the capacitor discharges,
with energy dissipated in the resister. Here, we determine the voltage drop
across the capacitor during charging and discharging.
The equations for the voltage drops across a capacitor and a resister are
given by
𝑉𝐶 = 𝑞/𝐶, 𝑉𝑅 = 𝑖𝑅, (2.24)
where 𝐶 is the capacitance and 𝑅 is the resistance. The charge 𝑞 and the current
𝑖 are related by
𝑑𝑞
𝑖= . (2.25)
𝑑𝑡
Kirchhoff’s voltage law states that the emf ℰ in any closed loop is equal to
the sum of the voltage drops in that loop. Applying Kirchhoff’s voltage law
when the switch is thrown to a results in
𝑉𝑅 + 𝑉𝐶 = ℰ. (2.26)
Using (2.24) and (2.25), the voltage drop across the resister can be written in
terms of the voltage drop across the capacitor as
𝑑𝑉𝐶
𝑉𝑅 = 𝑅𝐶 ,
𝑑𝑡
and (2.26) can be rewritten to yield the first-order linear differential equation
for 𝑉𝑐 given by
𝑑𝑉𝐶
+ 𝑉𝐶 /𝑅𝐶 = ℰ/𝑅𝐶, (2.27)
𝑑𝑡
with initial condition 𝑉𝐶 (0) = 0.
2.4. APPLICATIONS 27
𝜇(𝑡) = 𝑒𝑡/𝑅𝐶 ,
with solution (︁ )︁
𝑉𝐶 (𝑡) = ℰ 1 − 𝑒−𝑡/𝑅𝐶 .
The voltage starts at zero and rises exponentially to ℰ, with characteristic time
scale given by 𝑅𝐶.
When the switch is thrown to b, application of Kirchhoff’s voltage law results
in
𝑉𝑅 + 𝑉𝐶 = 0,
with corresponding differential equation
𝑑𝑉𝐶
+ 𝑉𝐶 /𝑅𝐶 = 0.
𝑑𝑡
Here, we assume that the capacitance is initially fully charged so that 𝑉𝐶 (0) = ℰ.
The solution, then, during the discharge phase is given by
𝑉𝐶 (𝑡) = ℰ𝑒−𝑡/𝑅𝐶 .
The voltage starts at ℰ and decays exponentially to zero, again with character-
istic time scale given by 𝑅𝐶.
𝑁 (𝑡) = 𝑁0 𝑒𝑟𝑡 ,
where 𝑁0 is the initial population size. However, when the population growth
is constrained by limited resources, a heuristic modification to the Malthusian
growth model results in the Verhulst equation,
(︂ )︂
𝑑𝑁 𝑁
= 𝑟𝑁 1 − , (2.28)
𝑑𝑡 𝐾
where 𝐾 is called the carrying capacity of the environment. Making (2.28)
dimensionless using 𝜏 = 𝑟𝑡 and 𝑥 = 𝑁/𝐾 leads to the logistic equation,
𝑑𝑥
= 𝑥(1 − 𝑥),
𝑑𝜏
28 CHAPTER 2. FIRST-ORDER ODES
where we may assume the initial condition 𝑥(0) = 𝑥0 > 0. Separating variables
and integrating ∫︁ 𝑥 ∫︁ 𝜏
𝑑𝑥
= 𝑑𝜏.
𝑥0 𝑥(1 − 𝑥) 0
The integral on the left-hand-side can be done using the method of partial
fractions:
1 𝑎 𝑏
= +
𝑥(1 − 𝑥) 𝑥 1−𝑥
𝑎 + (𝑏 − 𝑎)𝑥
= ;
𝑥(1 − 𝑥)
𝑥(1 − 𝑥0 )
= 𝑒𝜏 ,
𝑥0 (1 − 𝑥)
𝑥(1 − 𝑥0 ) = 𝑥0 𝑒𝜏 − 𝑥𝑥0 𝑒𝜏 ,
𝑥(1 − 𝑥0 + 𝑥0 𝑒𝜏 ) = 𝑥0 𝑒𝜏 ,
𝑥0
𝑥= . (2.29)
𝑥0 + (1 − 𝑥0 )𝑒−𝜏
lim 𝑁 (𝑡) = 𝐾.
𝑡→∞
The population, therefore, grows in size until it reaches the carrying capacity of
its environment.
Chapter 3
Second-order linear
differential equations with
constant coefficients
Reference: Boyce and DiPrima, Chapter 3
29
30 CHAPTER 3. SECOND-ORDER ODES, CONSTANT COEFFICIENTS
The values 𝑥1 and 𝑢1 at the time 𝑡1 = 𝑡0 + Δ𝑡 are then used as new initial values
to march the solution forward to time 𝑡2 = 𝑡1 + Δ𝑡. As long as 𝑓 (𝑡, 𝑥, 𝑢) is a
well-behaved function, the numerical solution converges to the unique solution
of the ode as Δ𝑡 → 0.
𝑥
¨ + 𝑝(𝑡)𝑥˙ + 𝑞(𝑡)𝑥 = 0; (3.4)
and suppose that 𝑥 = 𝑋1 (𝑡) and 𝑥 = 𝑋2 (𝑡) are solutions to (3.4). We consider
a linear combination of 𝑋1 and 𝑋2 by letting
= 𝑐1 × 0 + 𝑐2 × 0
= 0,
𝑐1 𝑋1 (𝑡0 ) + 𝑐2 𝑋2 (𝑡0 ) = 𝑥0 ,
𝑐1 𝑋˙ 1 (𝑡0 ) + 𝑐2 𝑋˙ 2 (𝑡0 ) = 𝑢0 , (3.7)
which is observed to be a system of two linear equations for the two unknowns
𝑐1 and 𝑐2 . Solution of (3.7) by standard methods results in
Our ansatz has thus converted a differential equation into an algebraic equation.
Equation (3.11) is called the characteristic equation of (3.9). Using the quadratic
formula, the two solutions of the characteristic equation (3.11) are given by
1 (︁ √︀ )︁
𝑟± = −𝑏 ± 𝑏2 − 4𝑎𝑐 .
2𝑎
There are three cases to consider: (1) if 𝑏2 − 4𝑎𝑐 > 0, then the two roots are
distinct and real; (2) if 𝑏2 −4𝑎𝑐 < 0, then the two roots are distinct and complex
conjugates of each other; (3) if 𝑏2 − 4𝑎𝑐 = 0, then the two roots are degenerate
and there is only one real root. We will consider these three cases in turn.
The unknown constants 𝑐1 and 𝑐2 can then be determined by the given initial
conditions 𝑥(𝑡0 ) = 𝑥0 and 𝑥(𝑡
˙ 0 ) = 𝑢0 . We now present two examples.
Example 1: Solve 𝑥
¨ + 5𝑥˙ + 6𝑥 = 0 with 𝑥(0) = 2, 𝑥(0)
˙ = 3, and find the
maximum value attained by 𝑥.
view tutorial
We take as our ansatz 𝑥 = 𝑒𝑟𝑡 and obtain the characteristic equation
𝑟2 + 5𝑟 + 6 = 0,
which factors to
(𝑟 + 3)(𝑟 + 2) = 0.
The general solution to the ode is thus
𝑥(𝑡)
˙ = −2𝑐1 𝑒−2𝑡 − 3𝑐2 𝑒−3𝑡 .
Use of the initial conditions then results in two equations for the two unknown
constant 𝑐1 and 𝑐2 :
𝑐1 + 𝑐2 = 2,
−2𝑐1 − 3𝑐2 = 3.
Adding three times the first equation to the second equation yields 𝑐1 = 9; and
the first equation then yields 𝑐2 = 2 − 𝑐1 = −7. Therefore, the unique solution
that satisfies both the ode and the initial conditions is
Note that although both exponential terms decay in time, their sum increases
initially since 𝑥(0)
˙ > 0. The maximum value of 𝑥 occurs at the time 𝑡𝑚 when
𝑥˙ = 0, or
𝑡𝑚 = ln (7/6) .
The maximum 𝑥𝑚 = 𝑥(𝑡𝑚 ) is then determined to be
𝑥𝑚 = 108/49.
and
𝑍1 − 𝑍2
𝑋2 (𝑡) =
2𝑖
𝑒 − 𝑒−𝑖𝜇𝑡
(︂ 𝑖𝜇𝑡 )︂
𝜆𝑡
=𝑒
2𝑖
= 𝑒𝜆𝑡 sin 𝜇𝑡.
Having found the two real solutions 𝑋1 (𝑡) and 𝑋2 (𝑡), we can then apply the
principle of superposition a second time to determine the general solution 𝑥(𝑡):
It is best to memorize this result. The real part of the roots of the characteristic
equation goes into the exponential function; the imaginary part goes into the
argument of cosine and sine.
Example 1: Solve 𝑥
¨ + 𝑥 = 0 with 𝑥(0) = 𝑥0 and 𝑥(0)
˙ = 𝑢0 .
view tutorial
The characteristic equation is
𝑟2 + 1 = 0,
with roots
𝑟± = ±𝑖.
The general solution of the ode is therefore
The derivative is
𝑥(𝑡)
˙ = −𝐴 sin 𝑡 + 𝐵 cos 𝑡.
𝑥(0) = 𝐴 = 𝑥0 , 𝑥(0)
˙ = 𝐵 = 𝑢0 ;
¨ − 𝑥 = 0 as 𝑥(𝑡) =
Recall that we wrote the analogous solution to the ode 𝑥
𝑥0 cosh 𝑡 + 𝑢0 sinh 𝑡.
Example 2: Solve 𝑥
¨ + 𝑥˙ + 𝑥 = 0 with 𝑥(0) = 1 and 𝑥(0)
˙ = 0.
The characteristic equation is
𝑟2 + 𝑟 + 1 = 0,
with roots √
1 3
𝑟± = − ± 𝑖 .
2 2
The general solution of the ode is therefore
(︃ √ √ )︃
− 12 𝑡 3 3
𝑥(𝑡) = 𝑒 𝐴 cos 𝑡 + 𝐵 sin 𝑡 .
2 2
The derivative is
(︃ √ √ )︃
1 1 3 3
𝑥(𝑡)
˙ = − 𝑒− 2 𝑡 𝐴 cos 𝑡 + 𝐵 sin 𝑡
2 2 2
√ (︃ √ √ )︃
3 −1𝑡 3 3
+ 𝑒 2 −𝐴 sin 𝑡 + 𝐵 cos 𝑡 .
2 2 2
𝐴 = 1,
√
−1 3
𝐴+ 𝐵 = 0;
2 2
or √
3
𝐴 = 1, 𝐵= .
3
Therefore,
(︃ √ √ √ )︃
− 21 𝑡 3 3 3
𝑥(𝑡) = 𝑒 cos 𝑡+ sin 𝑡 .
2 3 2
36 CHAPTER 3. SECOND-ORDER ODES, CONSTANT COEFFICIENTS
Example: Solve 𝑥
¨ + 2𝑥˙ + 𝑥 = 0 with 𝑥(0) = 1 and 𝑥(0)
˙ = 0.
The characteristic equation is
𝑟2 + 2𝑟 + 1 = (𝑟 + 1)2
= 0,
which has a repeated root given by 𝑟 = −1. Therefore, the general solution to
the ode is
𝑥(𝑡) = 𝑐1 𝑒−𝑡 + 𝑐2 𝑡𝑒−𝑡 ,
with derivative
𝑥(𝑡)
˙ = −𝑐1 𝑒−𝑡 + 𝑐2 𝑒−𝑡 − 𝑐2 𝑡𝑒−𝑡 .
Applying the initial conditions, we have
𝑐1 = 1,
−𝑐1 + 𝑐2 = 0,
so that 𝑐1 = 𝑐2 = 1. Therefore, the solution is
𝑥(𝑡) = (1 + 𝑡)𝑒−𝑡 .
so that (3.17) solves (3.15), and the two free constants in 𝑥ℎ can be used to
satisfy the initial conditions.
We will consider here only the constant coefficient case. We now illustrate
the solution method by an example.
𝑟2 − 3𝑟 − 4 = (𝑟 − 4)(𝑟 + 1)
= 0,
so that
𝑥ℎ (𝑡) = 𝑐1 𝑒4𝑡 + 𝑐2 𝑒−𝑡 .
Second, we find a particular solution of the inhomogeneous equation. The form
of the particular solution is chosen such that the exponential will cancel out of
both sides of the ode. The ansatz we choose is
4𝐴 − 6𝐴 − 4𝐴 = 3,
𝑥(𝑡)
˙ = 4𝑐1 𝑒4𝑡 − 𝑐2 𝑒−𝑡 − 𝑒2𝑡 .
or
3
𝑐1 + 𝑐2 = ,
2
4𝑐1 − 𝑐2 = 1.
This system of linear equations can be solved for 𝑐1 by adding the equations
to obtain 𝑐1 = 1/2, after which 𝑐2 = 1 can be determined from the first equa-
tion. Therefore, the solution for 𝑥(𝑡) that satisfies both the ode and the initial
3.5. INHOMOGENEOUS ODES 39
conditions is given by
1 4𝑡 1 2𝑡
𝑥(𝑡) = 𝑒 − 𝑒 + 𝑒−𝑡
2 2
1 4𝑡 (︀
= 𝑒 1 − 𝑒−2𝑡 + 2𝑒−5𝑡 ,
)︀
2
where we have grouped the terms in the solution to better display the asymptotic
behavior for large 𝑡.
We now find particular solutions for some relatively simple inhomogeneous
terms using this method of undetermined coefficients.
¨ − 3𝑥˙ − 4𝑥 = 2 sin 𝑡.
Example: Find a particular solution of 𝑥
view tutorial
We show two methods for finding a particular solution. The first more direct
method tries the ansatz
where the argument of cosine and sine must agree with the argument of sine in
the inhomogeneous term. The cosine term is required because the derivative of
sine is cosine. Upon substitution into the differential equation, we obtain
or regrouping terms,
This equation is valid for all 𝑡, and in particular for 𝑡 = 0 and 𝜋/2, for which
the sine and cosine functions vanish. For these two values of 𝑡, we find
5𝐴 + 3𝐵 = 0, 3𝐴 − 5𝐵 = 2;
and rewrite the differential equation in complex form. We can rewrite the equa-
tion in one of two ways. On the one hand, if we use sin 𝑡 = Re{−𝑖𝑒𝑖𝑡 }, then the
differential equation is written as
and by equating the real and imaginary parts, this equation becomes the two
real differential equations
¨ − 3𝑥˙ − 4𝑥 = 2 sin 𝑡,
𝑥 𝑦¨ − 3𝑦˙ − 4𝑦 = −2 cos 𝑡.
¨ − 3𝑥˙ − 4𝑥 = 2 cos 𝑡,
𝑥 𝑦¨ − 3𝑦˙ − 4𝑦 = 2 sin 𝑡.
𝑧(𝑡) = 𝐶𝑒𝑖𝑡 ,
−𝐶 − 3𝑖𝐶 − 4𝐶 = 2;
or solving for 𝐶:
−2
𝐶=
5 + 3𝑖
−2(5 − 3𝑖)
=
(5 + 3𝑖)(5 − 3𝑖)
−10 + 6𝑖
=
34
−5 + 3𝑖
= .
17
Therefore,
𝑥𝑝 = Im{𝑧𝑝 }
{︂ }︂
1
= Im (−5 + 3𝑖)(cos 𝑡 + 𝑖 sin 𝑡)
17
1
= (3 cos 𝑡 − 5 sin 𝑡).
17
¨ + 𝑥˙ − 2𝑥 = 𝑡2 .
Example: Find a particular solution of 𝑥
view tutorial
The correct ansatz here is the polynomial
𝑥(𝑡) = 𝐴𝑡2 + 𝐵𝑡 + 𝐶.
or
−2𝐴𝑡2 + 2(𝐴 − 𝐵)𝑡 + (2𝐴 + 𝐵 − 2𝐶)𝑡0 = 𝑡2 .
Equating powers of 𝑡,
−2𝐴 = 1, 2(𝐴 − 𝐵) = 0, 2𝐴 + 𝐵 − 2𝐶 = 0;
and solving,
1 1 3
𝐴=− , 𝐵=− , 𝐶=− .
2 2 4
The particular solution is therefore
1 1 3
𝑥𝑝 (𝑡) = − 𝑡2 − 𝑡 − .
2 2 4
𝑟 + 2 = 0, or 𝑟 = −2.
To find a particular solution, we try the ansatz 𝑥𝑝 (𝑡) = 𝐴𝑒−𝑡 , and upon substi-
tution
−𝐴 + 2𝐴 = 1, or 𝐴 = 1.
Therefore, the general solution to the ode is
3
𝑥(0) = = 𝑐 + 1,
4
so that 𝑐 = −1/4. Hence,
1
𝑥(𝑡) = 𝑒−𝑡 − 𝑒−2𝑡
(︂ 4 )︂
−𝑡 1
=𝑒 1 − 𝑒−𝑡 .
4
42 CHAPTER 3. SECOND-ORDER ODES, CONSTANT COEFFICIENTS
3.7 Resonance
view tutorial
Resonance occurs when the frequency of the inhomogeneous term matches the
frequency of the homogeneous solution. To illustrate resonance in its simplest
embodiment, we consider the second-order linear inhomogeneous ode
Our main goal is to determine what happens to the solution in the limit 𝜔 → 𝜔0 .
The homogeneous equation has characteristic equation
𝑟2 + 𝜔02 = 0,
−𝜔 2 𝐴 + 𝜔02 𝐴 = 𝑓,
or
𝑓
𝐴= .
𝜔02 − 𝜔 2
Therefore,
𝑓
𝑥𝑝 (𝑡) = cos 𝜔𝑡.
𝜔02 − 𝜔2
Our general solution is thus
𝑓
𝑥(𝑡) = 𝑐1 cos 𝜔0 𝑡 + 𝑐2 sin 𝜔0 𝑡 + cos 𝜔𝑡,
𝜔02 − 𝜔 2
with derivative
𝑓𝜔
𝑥(𝑡)
˙ = 𝜔0 (𝑐2 cos 𝜔0 𝑡 − 𝑐1 sin 𝜔0 𝑡) − sin 𝜔𝑡.
𝜔02 − 𝜔 2
so that
𝑓 𝑢0
𝑐1 = 𝑥0 − , 𝑐2 = .
𝜔02 − 𝜔 2 𝜔0
3.7. RESONANCE 43
Therefore, the solution to the ode that satisfies the initial conditions is
(︂ )︂
𝑓 𝑢0 𝑓
𝑥(𝑡) = 𝑥0 − 2 2
cos 𝜔0 𝑡 + sin 𝜔0 𝑡 + 2 cos 𝜔𝑡
𝜔0 − 𝜔 𝜔0 𝜔0 − 𝜔 2
𝑢0 𝑓 (cos 𝜔𝑡 − cos 𝜔0 𝑡)
= 𝑥0 cos 𝜔0 𝑡 + sin 𝜔0 𝑡 + ,
𝜔0 𝜔02 − 𝜔 2
¨ + 𝜔02 𝑥 = 𝑓 cos 𝜔0 𝑡,
𝑥
and find a particular solution directly. The particular solution is the real part
of the particular solution of
𝑧¨ + 𝜔02 𝑧 = 𝑓 𝑒𝑖𝜔0 𝑡 ,
𝑧𝑝 = 𝐴𝑡𝑒𝑖𝜔0 𝑡 .
We have
𝑧˙𝑝 = 𝐴𝑒𝑖𝜔0 𝑡 (1 + 𝑖𝜔0 𝑡) , 𝑧¨𝑝 = 𝐴𝑒𝑖𝜔0 𝑡 2𝑖𝜔0 − 𝜔02 𝑡 ;
(︀ )︀
= 2𝑖𝜔0 𝐴𝑒𝑖𝜔0 𝑡
= 𝑓 𝑒𝑖𝜔0 𝑡 .
Therefore,
𝑓
𝐴= ,
2𝑖𝜔0
44 CHAPTER 3. SECOND-ORDER ODES, CONSTANT COEFFICIENTS
and
𝑓 𝑡 𝑖𝜔0 𝑡
𝑥𝑝 = Re{ 𝑒 }
2𝑖𝜔0
𝑓 𝑡 sin 𝜔0 𝑡
= ,
2𝜔0
the same result as (3.23).
¨ − 3𝑥˙ − 4𝑥 = 5𝑒−𝑡 .
Example: Find a particular solution of 𝑥
view tutorial
If we naively try the ansatz
𝑥 = 𝐴𝑒−𝑡 ,
and substitute this into the inhomogeneous differential equation, we obtain
𝐴 + 3𝐴 − 4𝐴 = 5,
𝑥 = 𝐴𝑡𝑒−𝑡 ,
𝑥𝑝 = −𝑡𝑒−𝑡 .
𝑚¨
𝑥 + 𝛾 𝑥˙ + 𝑘𝑥 = 𝐹 cos 𝜔𝑡, (3.24)
where 𝑚 > 0 is the oscillator’s mass, 𝛾 > 0 is the damping coefficient, 𝑘 >
0 is the spring constant, and 𝐹 is the amplitude of the external force. The
homogeneous equation has characteristic equation
𝑚𝑟2 + 𝛾𝑟 + 𝑘 = 0,
3.8. DAMPED RESONANCE 45
so that
𝛾 1 √︀ 2
𝑟± = − ± 𝛾 − 4𝑚𝑘.
2𝑚 2𝑚
When 𝛾 2 − 4𝑚𝑘 < 0, the motion of the unforced oscillator is said to be under-
damped; when 𝛾 2 − 4𝑚𝑘 > 0, overdamped; and when 𝛾 2 − 4𝑚𝑘 = 0, critically
damped. For all three types of damping, the roots of the characteristic equa-
tion satisfy Re(𝑟± ) < 0. Therefore, both linearly independent homogeneous
solutions decay exponentially to zero, and the long-time asymptotic solution
of (3.24) reduces to the (non-decaying) particular solution. Since the initial
conditions are satisfied by the free constants multiplying the (decaying) homo-
geneous solutions, the long-time asymptotic solution is independent of the initial
conditions.
If we are only interested in the long-time asymptotic solution of (3.24), we
can proceed directly to the determination of a particular solution. As before,
we consider the complex ode
𝑧 + 𝛾 𝑧˙ + 𝑘𝑧 = 𝐹 𝑒𝑖𝜔𝑡 ,
𝑚¨
with 𝑥𝑝 = Re(𝑧𝑝 ). With the ansatz 𝑧𝑝 = 𝐴𝑒𝑖𝜔𝑡 , we have
−𝑚𝜔 2 𝐴 + 𝑖𝛾𝜔𝐴 + 𝑘𝐴 = 𝐹,
or
𝐹
𝐴= .
(𝑘 − 𝑚𝜔 2 ) + 𝑖𝛾𝜔
√︀
To simplify, we define 𝜔0 = 𝑘/𝑚, which corresponds to the natural frequency
of the undamped oscillator, and define Γ = 𝛾𝜔/𝑚 and 𝑓 = 𝐹/𝑚. Therefore,
𝑓
𝐴=
(𝜔02 − 𝜔 2 ) + 𝑖Γ
(︂ )︂ (3.25)
𝑓 (︀ 2
(𝜔0 − 𝜔 2 ) − 𝑖Γ .
)︀
= 2 2 2 2
(𝜔0 − 𝜔 ) + Γ
To determine 𝑥𝑝 , we utilize the polar form of a complex number. The complex
𝑖𝜑
number 𝑧 = 𝑥 + 𝑖𝑦 can
√︀ be written in polar form as 𝑧 = 𝑟𝑒 , where 𝑥 = 𝑟 cos 𝜑,
𝑦 = 𝑟 sin 𝜑, and 𝑟 = 𝑥2 + 𝑦 2 , tan 𝜑 = 𝑦/𝑥. We therefore write
(𝜔02 − 𝜔 2 ) − 𝑖Γ = 𝑟𝑒𝑖𝜑 ,
with
Γ
√︁
𝑟= (𝜔02 − 𝜔 2 )2 + Γ2 , tan 𝜑 = .
(𝜔 2 − 𝜔02 )
Using the polar form, 𝐴 in (3.25) becomes
(︃ )︃
𝑓
𝐴 = √︀ 2 𝑒𝑖𝜑 ,
(𝜔0 − 𝜔 2 )2 + Γ2
and 𝑥𝑝 = Re(𝐴𝑒𝑖𝜔𝑡 ) becomes
(︃ )︃
𝑓 (︁ )︁
𝑥𝑝 = √︀ Re 𝑒𝑖(𝜔𝑡+𝜑)
(𝜔02 − 𝜔 2 )2 + Γ2
(︃ )︃ (3.26)
𝑓
= √︀ cos (𝜔𝑡 + 𝜑).
(𝜔02 − 𝜔 2 )2 + Γ2
46 CHAPTER 3. SECOND-ORDER ODES, CONSTANT COEFFICIENTS
The particular solution given by (3.26), with 𝜔02 = 𝑘/𝑚, Γ = 𝛾𝜔/𝑚, 𝑓 = 𝐹/𝑚,
and tan 𝜑 = 𝛾𝜔/𝑚(𝜔 2 − 𝜔02 ), is the long-time asymptotic solution of the forced,
damped, harmonic oscillator equation (3.24).
We conclude with a couple of observations about (3.26). First, if the forcing
frequency 𝜔 is equal to the natural frequency 𝜔0 of the undamped oscillator, then
𝐴 = −𝑖𝐹/𝛾𝜔0 , and 𝑥𝑝 = (𝐹/𝛾𝜔0 ) sin 𝜔0 𝑡. The oscillator position is observed
to be 𝜋/2 out of phase with the external force, or in other words, the velocity
of the oscillator, not the position, is in phase with the force. Second, the value
of the forcing frequency 𝜔𝑚 that maximizes the amplitude of oscillation is the
value of 𝜔 that minimizes the denominator of (3.26). To determine 𝜔𝑚 we thus
minimize the function 𝑔(𝜔 2 ) with respect to 𝜔 2 , where
𝛾 2 𝜔2
𝑔(𝜔 2 ) = (𝜔02 − 𝜔 2 )2 + .
𝑚2
Taking the derivative of 𝑔 with respect to 𝜔 2 and setting this to zero to determine
𝜔𝑚 yields
2 𝛾2
2(𝜔𝑚 − 𝜔02 ) + 2 = 0,
𝑚
or
2 𝛾2
𝜔𝑚 = 𝜔02 − .
2𝑚2
We can interpret this result by saying that damping lowers the “resonance”
frequency of the undamped oscillator.
Chapter 4
The improper integral given by (4.1) diverges if 𝑓 (𝑡) grows faster than 𝑒𝑠𝑡 for
large 𝑡. Accordingly, some restriction on the range of 𝑠 may be required to
guarantee convergence of (4.1), and we will assume without further elaboration
that these restrictions are always satisfied.
The Laplace transform is a linear transformation. We have
∫︁ ∞
𝑒−𝑠𝑡 𝑐1 𝑓1 (𝑡) + 𝑐2 𝑓2 (𝑡) 𝑑𝑡
(︀ )︀
ℒ{𝑐1 𝑓1 (𝑡) + 𝑐2 𝑓2 (𝑡)} =
0
∫︁ ∞ ∫︁ ∞
= 𝑐1 𝑒−𝑠𝑡 𝑓1 (𝑡)𝑑𝑡 + 𝑐2 𝑒−𝑠𝑡 𝑓2 (𝑡)𝑑𝑡
0 0
= 𝑐1 ℒ{𝑓1 (𝑡)} + 𝑐2 ℒ{𝑓2 (𝑡)}.
47
48 CHAPTER 4. THE LAPLACE TRANSFORM
1. 𝑒𝑎𝑡 𝑓 (𝑡) 𝐹 (𝑠 − 𝑎)
1
2. 1
𝑠
1
3. 𝑒𝑎𝑡
𝑠−𝑎
𝑛!
4. 𝑡𝑛
𝑠𝑛+1
𝑛!
5. 𝑡𝑛 𝑒𝑎𝑡
(𝑠 − 𝑎)𝑛+1
𝑏
6. sin 𝑏𝑡
𝑠2 + 𝑏2
𝑠
7. cos 𝑏𝑡
𝑠2 + 𝑏2
𝑏
8. 𝑒𝑎𝑡 sin 𝑏𝑡
(𝑠 − 𝑎)2 + 𝑏2
𝑠−𝑎
9. 𝑒𝑎𝑡 cos 𝑏𝑡
(𝑠 − 𝑎)2 + 𝑏2
2𝑏𝑠
10. 𝑡 sin 𝑏𝑡
(𝑠2 + 𝑏2 )2
𝑠2 − 𝑏2
11. 𝑡 cos 𝑏𝑡
(𝑠2 + 𝑏2 )2
𝑒−𝑐𝑠
12. 𝑢𝑐 (𝑡)
𝑠
13. 𝑢𝑐 (𝑡)𝑓 (𝑡 − 𝑐) 𝑒−𝑐𝑠 𝐹 (𝑠)
15. 𝑥(𝑡)
˙ 𝑠𝑋(𝑠) − 𝑥(0)
16. 𝑥
¨(𝑡) 𝑠2 𝑋(𝑠) − 𝑠𝑥(0) − 𝑥(0)
˙
ℒ{𝑒𝑎𝑡 } = ℒ{𝑒𝑎𝑡 · 1}
1 (4.2)
= .
𝑠−𝑎
The transform of 𝑡𝑛 (line 4) can be found by successive integration by parts. A
more interesting method uses Taylor series expansions for 𝑒𝑎𝑡 and 1/(𝑠 − 𝑎). We
have
{︃ ∞ }︃
∑︁ (𝑎𝑡)𝑛
𝑎𝑡
ℒ{𝑒 } = ℒ
𝑛=0
𝑛!
∞
(4.3)
∑︁ 𝑎𝑛
= ℒ{𝑡𝑛 }.
𝑛=0
𝑛!
Using (4.2), and equating the coefficients of powers of 𝑎 in (4.3) and (4.4), results
in line 4:
𝑛!
ℒ{𝑡𝑛 } = 𝑛+1 .
𝑠
The Laplace transform of 𝑡𝑛 𝑒𝑎𝑡 (line 5) can be found from line 1 and line 4:
𝑛!
ℒ{𝑡𝑛 𝑒𝑎𝑡 } = .
(𝑠 − 𝑎)𝑛+1
The Laplace transform of sin 𝑏𝑡 (line 6) may be found from the Laplace transform
of 𝑒𝑎𝑡 (line 3) using 𝑎 = 𝑖𝑏:
ℒ{sin 𝑏𝑡} = Im ℒ{𝑒𝑖𝑏𝑡 }
{︀ }︀
{︂ }︂
1
= Im
𝑠 − 𝑖𝑏
{︂ }︂
𝑠 + 𝑖𝑏
= Im
𝑠2 + 𝑏2
𝑏
= 2 .
𝑠 + 𝑏2
Similarly, the Laplace transform of cos 𝑏𝑡 (line 7) is
ℒ{cos 𝑏𝑡} = Re ℒ{𝑒𝑖𝑏𝑡 }
{︀ }︀
𝑠
= 2 .
𝑠 + 𝑏2
The transform of 𝑒𝑎𝑡 sin 𝑏𝑡 (line 8) can be found from the transform of sin 𝑏𝑡
(line 6) and line 1:
𝑏
ℒ{𝑒𝑎𝑡 sin 𝑏𝑡} = ;
(𝑠 − 𝑎)2 + 𝑏2
and similarly for the transform of 𝑒𝑎𝑡 cos 𝑏𝑡:
𝑠−𝑎
ℒ{𝑒𝑎𝑡 cos 𝑏𝑡} = .
(𝑠 − 𝑎)2 + 𝑏2
The Laplace transform of 𝑡 sin 𝑏𝑡 (line 10) can be found from the Laplace trans-
form of 𝑡𝑒𝑎𝑡 (line 5 with 𝑛 = 1) using 𝑎 = 𝑖𝑏:
ℒ{𝑡 sin 𝑏𝑡} = Im ℒ{𝑡𝑒𝑖𝑏𝑡 }
{︀ }︀
{︂ }︂
1
= Im
(𝑠 − 𝑖𝑏)2
(𝑠 + 𝑖𝑏)2
{︂ }︂
= Im
(𝑠2 + 𝑏2 )2
2𝑏𝑠
= 2 .
(𝑠 + 𝑏2 )2
Similarly, the Laplace transform of 𝑡 cos 𝑏𝑡 (line 11) is
ℒ{𝑡 cos 𝑏𝑡} = Re ℒ{𝑡𝑒𝑖𝑏𝑡 }
{︀ }︀
(𝑠 + 𝑖𝑏)2
{︂ }︂
= Re
(𝑠2 + 𝑏2 )2
𝑠2 − 𝑏2
= 2 .
(𝑠 + 𝑏2 )2
4.2. SOLUTION OF INITIAL VALUE PROBLEMS 51
by parts and using the just derived result for the first derivative. We let
𝑢 = 𝑒−𝑠𝑡 𝑑𝑣 = 𝑥
¨𝑑𝑡
−𝑠𝑡
𝑑𝑢 = −𝑠𝑒 𝑑𝑡 𝑣 = 𝑥,
˙
so that
∫︁ ∞ ]︀∞
∫︁ ∞
−𝑠𝑡
𝑒 𝑥
¨𝑑𝑡 = ˙ −𝑠𝑡 0
𝑥𝑒 +𝑠 𝑒−𝑠𝑡 𝑥𝑑𝑡
˙
0 0
(︀ )︀
= −𝑥(0)
˙ + 𝑠 𝑠𝑋(𝑠) − 𝑥(0)
= 𝑠2 𝑋(𝑠) − 𝑠𝑥(0) − 𝑥(0),
˙
−𝑠𝑡
where similarly we assume lim𝑡→∞ 𝑥(𝑡)𝑒
˙ = 0.
¨ − 𝑥−2𝑥
Example: Solve 𝑥 ˙ = 0 with 𝑥(0) = 1 and 𝑥(0)
˙ = 0 by two different
methods.
view tutorial
The characteristic equation of the ode is determined from the ansatz 𝑥 = 𝑒𝑟𝑡
and is
𝑟2 − 𝑟 − 2 = (𝑟 − 2)(𝑟 + 1) = 0.
The general solution of the ode is therefore
1 2𝑡 2 −𝑡
𝑥(𝑡) = 𝑒 + 𝑒 . (4.5)
3 3
We now solve this example using the Laplace transform. Taking the Laplace
transform of both sides of the ode, using the linearity of the transform, and
applying our result for the transform of the first and second derivatives, we find
or
(𝑠 − 1)𝑥(0) + 𝑥(0)
˙
𝑋(𝑠) = .
𝑠2 − 𝑠 − 2
Note that the denominator of the right-hand-side is just the quadratic from the
characteristic equation of the homogeneous ode, and that this factor arises from
the derivatives of the exponential term in the Laplace transform integral.
Applying the initial conditions, we find
𝑠−1
𝑋(𝑠) = . (4.6)
(𝑠 − 2)(𝑠 + 1)
𝑠−1 𝑎 𝑏
= + . (4.7)
(𝑠 − 2)(𝑠 + 1) 𝑠−2 𝑠+1
The cover-up method can be used to solve for 𝑎 and 𝑏. We multiply both sides
of (4.7) by 𝑠 − 2 and put 𝑠 = 2 to isolate 𝑎:
]︂
𝑠−1
𝑎=
𝑠+1 𝑠=2
1
= .
3
4.2. SOLUTION OF INITIAL VALUE PROBLEMS 53
Example: Solve 𝑥
¨ + 𝑥 = sin 2𝑡 with 𝑥(0) = 2 and 𝑥(0)
˙ = 1 by Laplace
transform methods.
Taking the Laplace transform of both sides of the ode, we find
2𝑠 + 1 2
𝑋(𝑠) = + 2 .
𝑠 + 1 (𝑠 + 1)(𝑠2 + 4)
2
To determine the inverse Laplace transform from Table 4.1, we perform a partial
fraction expansion of the second term:
2 𝑎𝑠 + 𝑏 𝑐𝑠 + 𝑑
= 2 + . (4.8)
(𝑠2 2
+ 1)(𝑠 + 4) 𝑠 + 1 𝑠2 + 4
2𝑠 + 1 2/3 2/3
𝑋(𝑠) = 2
+ 2 − 2
𝑠 + 1 𝑠 + 1 (𝑠 + 4)
2𝑠 5/3 2/3
= 2 + 2 − 2 .
𝑠 + 1 𝑠 + 1 (𝑠 + 4)
From lines 6 and 7 of Table 4.1, we obtain the solution by taking inverse Laplace
transforms of the three terms separately, where 𝑏 = 1 in the first two terms, and
𝑏 = 2 in the third term:
5 1
𝑥(𝑡) = 2 cos 𝑡 + sin 𝑡 − sin 2𝑡.
3 3
54 CHAPTER 4. THE LAPLACE TRANSFORM
The step-up, step-down function—zero for 𝑡 < 𝑎, one for 𝑎 ≤ 𝑡 < 𝑏, and zero
for 𝑡 ≥ 𝑏—is defined as
⎧
⎨ 0, 𝑡 < 𝑎;
𝑢𝑎 (𝑡) − 𝑢𝑏 (𝑡) = 1, 𝑎 ≤ 𝑡 < 𝑏; (4.11)
0, 𝑡 ≥ 𝑏.
⎩
x x=f(t)
Figure 4.1: A linearly increasing function which turns into a sinusoidal function.
Using the Heaviside function 𝑢𝑡* , the function 𝑓 (𝑡) can be written in a single
line as (︀ )︀
𝑓 (𝑡) = 𝑓1 (𝑡) + 𝑓2 (𝑡) − 𝑓1 (𝑡) 𝑢𝑡* (𝑡).
This example can be generalized to piecewise functions defined on multiple
intervals.
As a concrete example, suppose the inhomogeneous term is represented by
a linearly increasing function, which then turns into a sinusoidal function for
56 CHAPTER 4. THE LAPLACE TRANSFORM
This specific form of 𝑓 (𝑡) enables a relatively easy Laplace transform. We can
write
𝑓 (𝑡) = 𝑎𝑡 + ℎ(𝑡 − 𝑡* )𝑢𝑡* (𝑡),
where we have defined
ℎ(𝑡) = 𝑏 sin 𝜔𝑡 − 𝑎𝑡.
Using line 13, the Laplace transform of 𝑓 (𝑡) is
1 𝑏𝜔 𝑎
ℒ{𝑡} = , ℒ{ℎ(𝑡)} = − 2.
𝑠2 𝑠2 + 𝜔 2 𝑠
The usual view of the shifted Dirac delta function 𝛿(𝑡 − 𝑐) is that it is zero
everywhere except at 𝑡 = 𝑐, where it is infinite, and the integral over the Dirac
delta function is one. The Dirac delta function is technically not a function,
but is what mathematicians call a distribution. Nevertheless, in most cases of
practical interest, it can be treated like a function, where physical results are
obtained following a final integration.
There are many ways to represent the Dirac delta function as a limit of a
well-defined function. For our purposes, the most useful representation makes
use of the step-up, step-down function of (4.11):
1
𝛿(𝑡 − 𝑐) = lim (𝑢𝑐−𝜖 (𝑡) − 𝑢𝑐+𝜖 (𝑡)).
𝜖→0 2𝜖
Before taking the limit, the well-defined step-up, step-down function is zero
except over a small interval of width 2𝜖 centered at 𝑡 = 𝑐, over which it takes
the large value 1/2𝜖. The integral of this function is one, independent of the
value of 𝜖.
4.4. DISCONTINUOUS OR IMPULSIVE TERMS 57
The Laplace transform of the Dirac delta function is easily found by inte-
gration using the definition of the delta function:
∫︁ ∞
ℒ{𝛿(𝑡 − 𝑐)} = 𝑒−𝑠𝑡 𝛿(𝑡 − 𝑐)𝑑𝑡
0
= 𝑒−𝑐𝑠 .
𝑒−5𝑠 𝑒−20𝑠
2 𝑠2 𝑋(𝑠) − 𝑠𝑥(0) − 𝑥(0)
(︀ )︀
˙ + 𝑠𝑋(𝑠) − 𝑥(0) + 2𝑋(𝑠) = − .
𝑠 𝑠
Using the initial values and solving for 𝑋(𝑠), we find
𝑒−5𝑠 − 𝑒−20𝑠
𝑋(𝑠) = .
𝑠(2𝑠2 + 𝑠 + 2)
To determine the solution for 𝑥(𝑡) we now need to find the inverse Laplace
transform. The exponential functions can be dealt with using line 13 of Table
4.1. We write
𝑋(𝑠) = (𝑒−5𝑠 − 𝑒−20𝑠 )𝐻(𝑠),
where
1
𝐻(𝑠) = .
𝑠(2𝑠2 + 𝑠 + 2)
Then using line 13, we have
1 𝑎 𝑏𝑠 + 𝑐
= + 2 ,
𝑠(2𝑠2 + 𝑠 + 2) 𝑠 2𝑠 + 𝑠 + 2
1/2 𝑠 + 12
𝐻(𝑠) = − 2
𝑠 2𝑠 + 𝑠 + 2
𝑠 + 21
(︂ )︂
1 1
= − .
2 𝑠 𝑠2 + 12 𝑠 + 1
Inspecting Table 4.1, the first term can be transformed using line 2, and the
second term can be transformed using lines 8 and 9, provided we complete the
square of the denominator and then massage the numerator. That is, first we
complete the square:
(︂ )︂2
1 1 15
𝑠2 + 𝑠 + 1 = 𝑠 + + ;
2 4 16
and next we write
√︁
𝑠 + 14 + √115 15
(︀ )︀
𝑠 + 21 16
2 1 = )︀2 15 .
𝑠 + 2𝑠 + 1 1
(︀
𝑠 + 4 + 16
Therefore, the function 𝐻(𝑠) can be written as
⎛ √︁ ⎞
(︀ 1
)︀ (︂ )︂ 15
1 ⎝1 𝑠+ 4 1 16
𝐻(𝑠) = − − √ ⎠.
2 𝑠 (𝑠 + 41 )2 + 15 16 15 (𝑠 + 1 2
4 ) + 15
16
The first term is transformed using line 2, the second term using line 9 and the
third term using line 8. We finally obtain
√ √
(︂ (︂ )︂)︂
1 1
ℎ(𝑡) = 1 − 𝑒−𝑡/4 cos ( 15𝑡/4) + √ sin ( 15𝑡/4) , (4.13)
2 15
which when combined with (4.12) yields the rather complicated solution for
𝑥(𝑡).
We briefly comment that it is also possible to solve this example without
using the Laplace transform. The key idea is that both 𝑥 and 𝑥˙ are continuous
functions of 𝑡. Clearly from the form of the inhomogeneous term and the initial
conditions, 𝑥(𝑡) = 0 for 0 ≤ 𝑡 ≤ 5. We then solve the ode between 5 ≤ 𝑡 ≤ 20
with the inhomogeneous term equal to unity and initial conditions 𝑥(5) = 𝑥(5)
˙ =
0. This requires first finding the general homogeneous solution, next finding a
particular solution, and then adding the homogeneous and particular solutions
and solving for the two unknown constants. To simplify the algebra, note that
the best ansatz to use to find the homogeneous solution is 𝑥(𝑡) = 𝑒𝑟(𝑡−5) , and
not 𝑥(𝑡) = 𝑒𝑟𝑡 . Finally, we solve the homogeneous ode for 𝑡 ≥ 20 using as
boundary conditions the previously determined values 𝑥(20) and 𝑥(20),
˙ where
we have made use of the continuity of 𝑥 and 𝑥. ˙ Here, the best ansatz to use
is 𝑥(𝑡) = 𝑒𝑟(𝑡−20) . The student may benefit by trying this as an exercise and
attempting to obtain a final solution that agrees with the form given by (4.12)
and (4.13).
4.4. DISCONTINUOUS OR IMPULSIVE TERMS 59
so that
(︂ )︂
1 −5𝑠 1
𝑋(𝑠) = 𝑒
2 𝑠2 + 21 𝑠 + 1
(︂ )︂
1 −5𝑠 1
= 𝑒
2 (𝑠 + 41 )2 + 15
⎛ √︁16 ⎞
√︂ 15
1 16 −5𝑠 ⎝ 16
= 𝑒 ⎠.
2 15 (𝑠 + 41 )2 + 15
16
The inverse Laplace transform may now be computed using lines 8 and 13 of
Table 4.1:
2 (︀√
𝑥(𝑡) = √ 𝑢5 (𝑡)𝑒−(𝑡−5)/4 sin 15(𝑡 − 5)/4 .
)︀
(4.14)
15
It is interesting to solve this example without using a Laplace transform.
Clearly, 𝑥(𝑡) = 0 up to the time of impulse at 𝑡 = 5. Furthermore, after the
impulse the ode is homogeneous and can be solved with standard methods. The
only difficulty is determining the initial conditions of the homogeneous ode at
𝑡 = 5+ .
When the inhomogeneous term is proportional to a delta-function, the solu-
tion 𝑥(𝑡) is continuous across the delta-function singularity, but the derivative
of the solution 𝑥(𝑡)
˙ is discontinuous. If we integrate the second-order ode across
the singularity at 𝑡 = 5 and consider 𝜖 → 0, only the second derivative term of
the left-hand-side survives, and
∫︁ 5+𝜖 ∫︁ 5+𝜖
2 𝑥
¨𝑑𝑡 = 𝛿(𝑡 − 5)𝑑𝑡
5−𝜖 5−𝜖
= 1.
Series solutions of
second-order linear
homogeneous differential
equations
Reference: Boyce and DiPrima, Chapter 5
where 𝑃 (𝑥), 𝑄(𝑥) and 𝑅(𝑥) are polynomials or convergent power series (around
𝑥 = 𝑥0 ), with no common polynomial factors (that could be divided out). The
value 𝑥 = 𝑥0 is called an ordinary point of (5.1) if 𝑃 (𝑥0 ) ̸= 0, and is called a
singular point if 𝑃 (𝑥0 ) = 0. Singular points will later be further classified as
regular singular points and irregular singular points. Our goal is to find two
independent solutions of (5.1), valid in a neighborhood about 𝑥 = 𝑥0 .
61
62 CHAPTER 5. SERIES SOLUTIONS
and
∞
∑︁
𝑦 ′′ (𝑥) = 𝑛(𝑛 − 1)𝑎𝑛 𝑥𝑛−2 .
𝑛=2
Substituting the power series for 𝑦 and its derivatives into the differential equa-
tion to be solved, we obtain
∞
∑︁ ∞
∑︁
𝑛−2
𝑛(𝑛 − 1)𝑎𝑛 𝑥 + 𝑎𝑛 𝑥𝑛 = 0. (5.2)
𝑛=2 𝑛=0
The power-series solution method requires combining the two sums on the left-
hand-side of (5.2) into a single power series in 𝑥. To shift the exponent of 𝑥𝑛−2
in the first sum upward by two to obtain 𝑥𝑛 , we need to shift the summation
index downward by two; that is,
∞
∑︁ ∞
∑︁
𝑛(𝑛 − 1)𝑎𝑛 𝑥𝑛−2 = (𝑛 + 2)(𝑛 + 1)𝑎𝑛+2 𝑥𝑛 .
𝑛=2 𝑛=0
For (5.3) to be satisfied, the coefficient of each power of 𝑥 must vanish separately.
(This can be proved by setting 𝑥 = 0 after successive differentiation.) We
therefore obtain the recurrence relation
𝑎𝑛
𝑎𝑛+2 = − , 𝑛 = 0, 1, 2, . . . .
(𝑛 + 2)(𝑛 + 1)
We observe that even and odd coefficients decouple. We thus obtain two inde-
pendent sequences starting with first term 𝑎0 or 𝑎1 . Developing these sequences,
we have for the sequence beginning with 𝑎0 :
𝑎0 ,
1
𝑎2 = − 𝑎0 ,
2
1 1
𝑎4 = − 𝑎2 = 𝑎0 ,
4·3 4·3·2
1 1
𝑎6 = − 𝑎4 = − 𝑎0 ;
6·5 6!
and the general coefficient in this sequence for 𝑛 = 0, 1, 2, . . . is
(−1)𝑛
𝑎2𝑛 = 𝑎0 .
(2𝑛)!
5.1. ORDINARY POINTS 63
𝑎1 ,
1
𝑎3 = − 𝑎1 ,
3·2
1 1
𝑎5 = − 𝑎3 = 𝑎1 ,
5·4 5·4·3·2
1 1
𝑎7 = − 𝑎5 = − 𝑎1 ;
7·6 7!
and the general coefficient in this sequence for 𝑛 = 0, 1, 2, . . . is
(−1)𝑛
𝑎2𝑛+1 = 𝑎1 .
(2𝑛 + 1)!
as expected.
In our next example, we will solve the Airy’s Equation. This differential
equation arises in the study of optics, fluid mechanics, and quantum mechanics.
We shift the first sum to 𝑥𝑛+1 by shifting the exponent up by three, i.e.,
∞
∑︁ ∞
∑︁
𝑛(𝑛 − 1)𝑎𝑛 𝑥𝑛−2 = (𝑛 + 3)(𝑛 + 2)𝑎𝑛+3 𝑥𝑛+1 .
𝑛=2 𝑛=−1
When combining the two sums in (5.4), we separate out the extra 𝑛 = −1 term
in the first sum given by 2𝑎2 . Therefore, (5.4) becomes
∞ (︂
∑︁ )︂
2𝑎2 + (𝑛 + 3)(𝑛 + 2)𝑎𝑛+3 − 𝑎𝑛 𝑥𝑛+1 = 0. (5.5)
𝑛=0
64 CHAPTER 5. SERIES SOLUTIONS
𝑎1 ,
1
𝑎4 = 𝑎1 ,
4·3
1
𝑎7 = 𝑎1 ,
7·6·4·3
1
𝑎10 = 𝑎1 .
10 · 9 · 7 · 6 · 4 · 3
The general solution for 𝑦 = 𝑦(𝑥), can therefore be written as
𝑥3 𝑥6 𝑥9 𝑥4 𝑥7 𝑥10
(︂ )︂ (︂ )︂
𝑦(𝑥) = 𝑎0 1 + + + + . . . + 𝑎1 𝑥 + + + + ...
6 180 12960 12 504 45360
= 𝑎0 𝑦0 (𝑥) + 𝑎1 𝑦1 (𝑥).
Suppose we would like to graph the solutions 𝑦 = 𝑦0 (𝑥) and 𝑦 = 𝑦1 (𝑥)
versus 𝑥 by solving the differential equation 𝑦 ′′ − 𝑥𝑦 = 0 numerically. What
initial conditions should we use? Clearly, 𝑦 = 𝑦0 (𝑥) solves the ode with initial
values 𝑦(0) = 1 and 𝑦 ′ (0) = 0, while 𝑦 = 𝑦1 (𝑥) solves the ode with initial values
𝑦(0) = 0 and 𝑦 ′ (0) = 1.
The numerical solutions, obtained using MATLAB, are shown in Fig. 5.1.
Note that the solutions oscillate for negative 𝑥 and grow exponentially for posi-
tive 𝑥. This can be understood by recalling that 𝑦 ′′ + 𝑦 = 0 has oscillatory sine
and cosine solutions and 𝑦 ′′ − 𝑦 = 0 has exponential hyperbolic sine and cosine
solutions.
5.2. REGULAR SINGULAR POINTS: CAUCHY-EULER EQUATIONS 65
Airy’s functions
2
1
y0
−1
−2
−10 −8 −6 −4 −2 0 2
1
y1
−1
−2
−10 −8 −6 −4 −2 0 2
x
if 𝑝(𝑥) and 𝑞(𝑥) have convergent Taylor series about 𝑥 = 𝑥0 , i.e., 𝑝(𝑥) and 𝑞(𝑥)
can be written as a power-series in (𝑥 − 𝑥0 ):
𝑝(𝑥) = 𝑝0 + 𝑝1 (𝑥 − 𝑥0 ) + 𝑝2 (𝑥 − 𝑥0 )2 + . . . ,
𝑞(𝑥) = 𝑞0 + 𝑞1 (𝑥 − 𝑥0 ) + 𝑞2 (𝑥 − 𝑥0 )2 + . . . ,
𝑥2 𝑦 ′′ + 𝛼𝑥𝑦 ′ + 𝛽𝑦 = 0, (5.8)
66 CHAPTER 5. SERIES SOLUTIONS
= 0.
As expected, the ode for 𝑌 = 𝑌 (𝜉) has constant coefficients, and with 𝑌 = 𝑒𝑟𝜉 ,
the characteristic equation for 𝑟 is given by (5.9). We now directly transfer
previous results obtained for the constant coefficient second-order linear homo-
geneous ode.
If (𝛼 − 1)2 − 4𝛽 = 0, there is one real root 𝑟 of (5.9). The general solution for
𝑌 is
𝑌 (𝜉) = 𝑒𝑟𝜉 (𝑐1 + 𝑐2 𝜉) ,
yielding
𝑦(𝑥) = |𝑥|𝑟 (𝑐1 + 𝑐2 ln |𝑥|) .
We now give examples illustrating these three cases.
0 = 2𝑟(𝑟 − 1) + 3𝑟 − 1
= 2𝑟2 + 𝑟 − 1
= (2𝑟 − 1)(𝑟 + 1).
Since the characteristic equation has two real roots, the general solution is given
by
1
𝑦(𝑥) = 𝑐1 𝑥 2 + 𝑐2 𝑥−1 .
68 CHAPTER 5. SERIES SOLUTIONS
We now encounter for the first time two-point boundary conditions, which can
be used to determine the coefficients 𝑐1 and 𝑐2 . Since y(0)=0, we must have
𝑐2 = 0. Applying the remaining condition 𝑦(1) = 1, we obtain the unique
solution √
𝑦(𝑥) = 𝑥.
Note that 𝑥 = 0 is called a singular point of the ode since the general solution
is singular at 𝑥 = 0 when 𝑐2 ̸= 0. Our boundary condition imposes that 𝑦(𝑥) is
finite at 𝑥 = 0 removing the singular solution. Nevertheless, 𝑦 ′ remains singular
at 𝑥 = 0. Indeed, this is why we imposed a two-point boundary condition rather
than specifying the value of 𝑦 ′ (0) (which is infinite).
0 = 𝑟(𝑟 − 1) + 𝑟 + 𝜋 2
= 𝑟2 + 𝜋2 ,
0 = 𝑟(𝑟 − 1) + 5𝑟 + 4
= 𝑟2 + 4𝑟 + 4
= (𝑟 + 2)2 ,
𝑐1 + 𝑐2 ln 𝑥
𝑦(𝑥) = .
𝑥2
The first boundary condition 𝑦(1) = 0 yields 𝑐1 = 0. The second boundary
condition 𝑦(𝑒) = 1 yields 𝑐2 = 𝑒2 . The solution is therefore
𝑒2 ln 𝑥
𝑦(𝑥) = .
𝑥2
Chapter 6
Systems of coupled linear equations can result, for example, from linear stability
analyses of nonlinear equations, and from normal mode analyses of coupled os-
cillators. Here, we will consider only the simplest case of a system of two coupled
first-order, linear, homogeneous equations with constant coefficients. These two
first-order equations are in fact equivalent to a single second-order equation, and
the methods of Chapter 3 could be used for solution. Nevertheless, viewing the
problem as a system of first-order equations introduces the important concept
of the phase space, and can easily be generalized to higher-order linear systems.
Ax = 0. (6.2)
When does there exist a nontrivial (not identically zero) solution for x? To
answer this question, we solve directly the system
𝑎𝑥1 + 𝑏𝑥2 = 0,
𝑐𝑥1 + 𝑑𝑥2 = 0.
Multiplying the first equation by 𝑑 and the second by 𝑏, and subtracting the
second equation from the first, results in
(𝑎𝑑 − 𝑏𝑐)𝑥1 = 0.
Similarly, multiplying the first equation by 𝑐 and the second by 𝑎, and subtract-
ing the first equation from the second, results in
(𝑎𝑑 − 𝑏𝑐)𝑥2 = 0.
69
70 CHAPTER 6. SYSTEMS OF EQUATIONS
Av = 𝜆v (6.3)
0 = det (A − 𝜆I)
⃒ ⃒
⃒ 𝑎−𝜆 𝑏 ⃒
= ⃒⃒ ⃒
𝑐 𝑑−𝜆 ⃒
= (𝑎 − 𝜆)(𝑑 − 𝜆) − 𝑏𝑐
= 𝜆2 − (𝑎 + 𝑑)𝜆 + (𝑎𝑑 − 𝑏𝑐).
where TrA is the trace, or sum of the diagonal elements, of the matrix A. If
𝜆 is an eigenvalue of A, then the corresponding eigenvector v may be found by
solving (︂ )︂ (︂ )︂
𝑎−𝜆 𝑏 𝑣1
= 0,
𝑐 𝑑−𝜆 𝑣2
where the equation of the second row will always be a multiple of the equation
of the first row. The eigenvector v has arbitrary normalization, and we may
always choose for convenience 𝑣1 = 1.
In the next section, we will see several examples of an eigenvector analysis.
6.2. COUPLED FIRST-ORDER EQUATIONS 71
ẋ = Ax. (6.7)
We will consider by example three cases separately: (i) eigenvalues of A are real
and there are two linearly independent eigenvectors; (ii) eigenvalues of A are
complex conjugates, and; (iii) A has only one linearly independent eigenvector.
These three cases are analogous to the cases considered previously when solving
the second-order, linear, constant-coefficient, homogeneous equation.
or using vector notation, written as (6.7). We take as our ansatz x(𝑡) = v𝑒𝜆𝑡 ,
where v and 𝜆 are independent of 𝑡. Upon substitution into (6.7), we obtain
𝜆v𝑒𝜆𝑡 = Av𝑒𝜆𝑡 ;
Av = 𝜆v. (6.8)
0 = det (A − 𝜆I)
= 𝜆2 − 2𝜆 − 3
= (𝜆 − 3)(𝜆 + 1).
1.5
0.5
0
2
x
−0.5
−1
−1.5
−2
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
x1
Figure 6.1: Phase space diagram for example with two real eigenvalues of oppo-
site sign.
¨1 − 2𝑥˙ 1 − 3𝑥1 = 0.
𝑥
The ansatz x = v𝑒𝜆𝑡 leads to the eigenvalue problem Av = 𝜆v, with A the
matrix above. The eigenvalues are determined from
0 = det (A − 𝜆I)
= 𝜆2 + 5𝜆 + 4
= (𝜆 + 4)(𝜆 + 1).
74 CHAPTER 6. SYSTEMS OF EQUATIONS
1.5
0.5
0
2
x
−0.5
−1
−1.5
−2
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
x1
Figure 6.2: Phase space diagram for example with two real eigenvalues of same
sign.
Taking the normalization 𝑣11 = 1 and 𝑣12 = 1, we obtain for the eigenvalues
and associated eigenvectors
(︂ )︂ (︂ )︂
𝜆1 = −4, v1 = √ 1 ; 𝜆2 = −1, v2 = √1 .
− 2/2 2
We show the phase space plot√in Fig. 6.2. If 𝑐2 = 0, the motion is along
the eigenvector v1 with 𝑥2 = −( 2/2)𝑥1 with eigenvalue √ 𝜆1 = −4 < 0. If
𝑐1 = 0, the motion is along the eigenvector v2 with 𝑥2 = 2𝑥1 with eigenvalue
𝜆2 = −1 < 0. When the eigenvalues are real and have the same sign, the origin
is called a node. A node may be attracting or repelling depending on whether
the eigenvalues are both negative (as is the case here) or positive. Observe that
the trajectories collapse onto the v2 eigenvector since 𝜆1 < 𝜆2 < 0 and decay is
more rapid along the v1 direction.
6.2. COUPLED FIRST-ORDER EQUATIONS 75
0 = det (A − 𝜆I)
5
= 𝜆2 + 𝜆 + .
4
Therefore, 𝜆 = −1/2±𝑖; and we observe that the eigenvalues occur as a complex
conjugate pair. We will denote the two eigenvalues as
1 ¯ = − 1 − 𝑖.
𝜆=− +𝑖 and 𝜆
2 2
¯ Therefore, the eigen-
Now, for A a real matrix, if Av = 𝜆v, then Av̄ = 𝜆v̄.
vectors also occur as a complex conjugate pair. The eigenvector v associated
with eigenvalue 𝜆 satisfies −𝑖𝑣1 + 𝑣2 = 0, and normalizing with 𝑣1 = 1, we have
(︂ )︂
1
v= .
𝑖
then two real functions can be constructed from the following linear combina-
tions of 𝑧 and 𝑧¯:
𝑧 + 𝑧¯ 𝑧 − 𝑧¯
= Re{𝑧(𝑡)} and = Im{𝑧(𝑡)}.
2 2𝑖
Thus the two real vector functions that can be constructed from our two complex
vector functions are
{︂(︂ )︂ }︂
1 1
Re{v𝑒𝜆𝑡 } = Re 𝑒(− 2 +𝑖)𝑡
𝑖
{︂(︂ )︂ }︂
1 1
= 𝑒− 2 𝑡 Re (cos 𝑡 + 𝑖 sin 𝑡)
𝑖
(︂ )︂
1 cos 𝑡
= 𝑒− 2 𝑡 ;
− sin 𝑡
76 CHAPTER 6. SYSTEMS OF EQUATIONS
1.5
0.5
0
2
x
−0.5
−1
−1.5
−2
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
x1
Figure 6.3: Phase space diagram for example with complex conjugate eigenval-
ues.
and
{︂(︂
)︂ }︂
𝜆𝑡 − 21 𝑡 1
Im{v𝑒 } = 𝑒 Im (cos 𝑡 + 𝑖 sin 𝑡)
𝑖
(︂ )︂
1 sin 𝑡
= 𝑒− 2 𝑡 .
cos 𝑡
Taking a linear superposition of these two real solutions yields the general so-
lution to the ode, given by
(︂ (︂ )︂ (︂ )︂)︂
1 cos 𝑡 sin 𝑡
x = 𝑒− 2 𝑡 𝐴 +𝐵 .
− sin 𝑡 cos 𝑡
The corresponding phase space diagram is shown in Fig. 6.3. We say the
origin is a spiral point. If the real part of the complex eigenvalue is negative,
as is the case here, then solutions spiral into the origin. If the real part of the
eigenvalue is positive, then solutions spiral out of the origin.
The direction of the spiral—here, it is clockwise—can be determined using
a concept from physics. If a particle of unit mass moves along a phase space
trajectory, then the angular momentum of the particle about the origin is equal
to the cross product of the position and velocity vectors: L = x × ẋ. With
both the position and velocity vectors lying in the two-dimensional phase space
plane, the angular momentum vector is perpendicular to this plane. With
x = (𝑥1 , 𝑥2 , 0), ẋ = (𝑥˙ 1 , 𝑥˙ 2 , 0),
then
L = (0, 0, 𝐿), with 𝐿 = 𝑥1 𝑥˙ 2 − 𝑥2 𝑥˙ 1 .
By the right-hand-rule, a clockwise rotation corresponds to 𝐿 < 0, and a coun-
terclockwise rotation to 𝐿 > 0.
6.2. COUPLED FIRST-ORDER EQUATIONS 77
1 1
𝑥1 𝑥˙ 2 − 𝑥2 𝑥˙ 1 = 𝑥1 (−𝑥1 − 𝑥2 ) − 𝑥2 (− 𝑥1 + 𝑥2 )
2 2
= −(𝑥21 + 𝑥22 )
< 0.
And since 𝐿 < 0, the spiral rotation is clockwise, as shown in Fig. 6.3.
0 = det (A − 𝜆I)
= 𝜆2 − 4𝜆 + 4
= (𝜆 − 2)2 .
and we need to find the missing second solution to be able to satisfy the initial
conditions. An ansatz of 𝑡 times the first solution is tempting, but will fail. Here,
we will cheat and find the missing second solution by solving the equivalent
second-order, homogeneous, constant-coefficient differential equation.
We already know that this second-order differential equation for 𝑥1 (𝑡) has a
characteristic equation with a degenerate eigenvalue given by 𝜆 = 2. Therefore,
the general solution for 𝑥1 is given by
1.5
0.5
0
2
x
−0.5
−1
−1.5
−2
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
x1
Figure 6.4: Phase space diagram for example with only one eigenvector.
so that
𝑥2 = 𝑥1 − 𝑥˙ 1
= (𝑐1 + 𝑡𝑐2 )𝑒2𝑡 − 2𝑐1 + (1 + 2𝑡)𝑐2 𝑒2𝑡
(︀ )︀
The second term of (6.11) is just 𝑡 times the first solution; however, this is
not sufficient. Indeed, the correct ansatz to find the second solution directly is
given by
x = (w + 𝑡v) 𝑒𝜆𝑡 , (6.12)
where 𝜆 and v is the eigenvalue and eigenvector of the first solution, and w
is an unknown vector to be determined. To illustrate this direct method, we
substitute (6.12) into ẋ = Ax, assuming Av = 𝜆v . Canceling the exponential,
we obtain
v + 𝜆 (w + 𝑡v) = Aw + 𝜆𝑡v.
Further canceling the common term 𝜆𝑡v and rewriting yields
(A − 𝜆I) w = v. (6.13)
6.3. NORMAL MODES 79
The first and second equation are the same, so that 𝑤2 = −(𝑤1 + 1). Therefore,
(︂ )︂
𝑤1
w=
−(𝑤1 + 1)
(︂ )︂ (︂ )︂
1 0
= 𝑤1 + .
−1 −1
Notice that the first term repeats the first found solution, i.e., a constant times
the eigenvector, and the second term is new. We therefore take 𝑤1 = 0 and
obtain (︂ )︂
0
w= ,
−1
as before.
The phase space diagram for this ode is shown in Fig. 6.4. The dark line is
the single eigenvector v of the matrix A. When there is only a single eigenvector,
the origin is called an improper node.
There is a definite counterclockwise rotation to the phase space trajectories,
and this can be confirmed from the calculation
𝐿 = 𝑥1 𝑥˙ 2 − 𝑥2 𝑥˙ 1
= 𝑥1 (𝑥1 + 3𝑥2 ) − 𝑥2 (𝑥1 − 𝑥2 )
= 𝑥21 + 2𝑥1 𝑥2 + 𝑥22
= (𝑥1 + 𝑥2 )2
> 0.
𝑥1 = −𝑘𝑥1 − 𝐾(𝑥1 − 𝑥2 ),
𝑚¨
𝑥2 = −𝑘𝑥2 − 𝐾(𝑥2 − 𝑥1 ).
𝑚¨
The equations for the coupled mass-spring system form a system of two second-
order linear homogeneous odes. In matrix form, 𝑚ẍ = Ax, or explicitly,
𝑑2
(︂ )︂ (︂ )︂ (︂ )︂
𝑥1 −(𝑘 + 𝐾) 𝐾 𝑥1
𝑚 2 = . (6.14)
𝑑𝑡 𝑥2 𝐾 −(𝑘 + 𝐾) 𝑥2
0 = det (A − 𝑚𝜆2 I)
⃒ ⃒
⃒ −(𝑘 + 𝐾) − 𝑚𝜆2 𝐾 ⃒
= ⃒⃒ ⃒
2 ⃒
𝐾 −(𝑘 + 𝐾) − 𝑚𝜆
= (𝑚𝜆2 + 𝑘 + 𝐾)2 − 𝐾 2 .
Since 𝜆21 , 𝜆22 < 0, both values of 𝜆 are imaginary, and 𝑥1 (𝑡) and 𝑥2 (𝑡) oscillate
with angular frequencies 𝜔1 = |𝜆1 | and 𝜔2 = |𝜆2 |, where
√︀ √︀
𝜔1 = 𝑘/𝑚, 𝜔2 = (𝑘 + 2𝐾)/𝑚.
−𝐾𝑣11 + 𝐾𝑣12 = 0,
so
√︀ that 𝑣11 = 𝑣12 . The normal mode associated with the frequency 𝜔1 =
𝑘/𝑚 thus follows a motion where 𝑥1 = 𝑥2 . Referring to Fig. 6.5, during this
motion the center spring length does not change, and the two masses oscillate
as if the center spring was absent (which is why the frequency of oscillation is
independent of 𝐾).
6.3. NORMAL MODES 81
𝐾𝑣21 + 𝐾𝑣22 = 0,
√︀ that 𝑣21 = −𝑣22 . The normal mode associated with the frequency 𝜔2 =
so
(𝑘 + 2𝐾)/𝑚 thus follows a motion where 𝑥1 = −𝑥2 . Again referring to
Fig. 6.5, during this motion the two equal masses symmetrically push or pull
against each side of the middle spring.
A general solution for x(𝑡) can be constructed from the eigenvalues and
eigenvectors. Our ansatz was x = v𝑒𝜆𝑡 , and for each of two eigenvectors v,
we have a pair of complex conjugate values for 𝜆. Accordingly, we first apply
the principle of superposition to obtain four real solutions, and √︀ then apply the
principle
√︀ again to obtain the general solution. With 𝜔1 = 𝑘/𝑚 and 𝜔2 =
(𝑘 + 2𝐾)/𝑚, the general solution is given by
(︂ )︂ (︂ )︂ (︂ )︂
𝑥1 1 1
= (𝐴 cos 𝜔1 𝑡 + 𝐵 sin 𝜔1 𝑡) + (𝐶 cos 𝜔2 𝑡 + 𝐷 sin 𝜔2 𝑡) ,
𝑥2 1 −1
where the now real constants 𝐴, 𝐵, 𝐶, and 𝐷 can be determined from the four
independent initial conditions, 𝑥1 (0), 𝑥2 (0), 𝑥˙ 1 (0), and 𝑥˙ 2 (0).
82 CHAPTER 6. SYSTEMS OF EQUATIONS
Chapter 7
Nonlinear differential
equations and bifurcation
theory
Reference: Strogatz, Sections 2.2, 2.4, 3.1, 3.2, 3.4, 6.3, 6.4, 8.2
𝑥˙ = 𝑓 (𝑥). (7.1)
𝜖˙ = 𝑓 (𝑥* + 𝜖)
= 𝑓 (𝑥* ) + 𝜖𝑓 ′ (𝑥* ) + . . .
= 𝜖𝑓 ′ (𝑥* ) + . . . .
83
84 CHAPTER 7. NONLINEAR DIFFERENTIAL EQUATIONS
The omitted terms in the Taylor series expansion are proportional to 𝜖2 , and
can be made negligible over a short time interval with respect to the kept term,
proportional to 𝜖, by taking 𝜖(0) sufficiently small. Therefore, at least over short
times, the differential equation to be considered, 𝜖˙ = 𝑓 ′ (𝑥* )𝜖, is linear and has
by now the familiar solution
′
𝜖(𝑡) = 𝜖(0)𝑒𝑓 (𝑥* )𝑡
.
The perturbation of the fixed point solution 𝑥(𝑡) = 𝑥* thus decays exponentially
if 𝑓 ′ (𝑥* ) < 0, and we say the fixed point is stable. If 𝑓 ′ (𝑥* ) > 0, the perturbation
grows exponentially and we say the fixed point is unstable. If 𝑓 ′ (𝑥* ) = 0, we
say the fixed point is marginally stable, and the next higher-order term in the
Taylor series expansion must be considered.
Example: Find all the fixed points of the logistic equation 𝑥˙ = 𝑥(1 − 𝑥)
and determine their stability.
where the function 𝐹 and all of its partial derivatives on the right-hand-side are
evaluated at the origin. Note that the Taylor series is constructed so that all
partial derivatives of the left-hand-side match those of the right-hand-side at
the origin.
We now consider the two-dimensional system given by
𝜖˙ = 𝑓 (𝑥* + 𝜖, 𝑦* + 𝛿)
𝜕𝑓 𝜕𝑓
=𝑓 +𝜖 +𝛿 + ...
𝜕𝑥 𝜕𝑦
𝜕𝑓 𝜕𝑓
=𝜖 +𝛿 + ....
𝜕𝑥 𝜕𝑦
𝛿˙ = 𝑔(𝑥* + 𝜖, 𝑦* + 𝛿)
𝜕𝑔 𝜕𝑔
=𝑔+𝜖 +𝛿 + ...
𝜕𝑥 𝜕𝑦
𝜕𝑔 𝜕𝑔
=𝜖 +𝛿 + ...,
𝜕𝑥 𝜕𝑦
where in the Taylor series 𝑓 , 𝑔 and all their partial derivatives are evaluated at
the fixed point (𝑥* , 𝑦* ). Neglecting higher-order terms in the Taylor series, we
thus have a system of odes for the perturbation, given in matrix form as
(︂ )︂ (︃ 𝜕𝑓 𝜕𝑓 )︃ (︂ )︂
𝑑 𝜖 𝜕𝑥 𝜕𝑦 𝜖
= 𝜕𝑔 𝜕𝑔 . (7.3)
𝑑𝑡 𝛿 𝜕𝑥 𝜕𝑦
𝛿
The two-by-two matrix in (7.3) is called the Jacobian matrix at the fixed point.
An eigenvalue analysis of the Jacobian matrix will typically yield two eigenvalues
𝜆1 and 𝜆2 . These eigenvalues may be real and distinct, complex conjugate pairs,
or repeated. The fixed point is stable (all perturbations decay exponentially)
if both eigenvalues have negative real parts. The fixed point is unstable (some
perturbations grow exponentially) if at least one of the eigenvalues has a positive
real part. Fixed points can be further classified as stable or unstable nodes,
unstable saddle points, stable or unstable spiral points, or stable or unstable
improper nodes.
Example: Find all the fixed points of the nonlinear system 𝑥˙ = 𝑥(3 −
𝑥 − 2𝑦), 𝑦˙ = 𝑦(2 − 𝑥 − 𝑦), and determine their stability.
view tutorial
The fixed points are determined by solving
There are four fixed points (𝑥* , 𝑦* ): (0, 0), (0, 2), (3, 0) and (1, 1). The Jacobian
matrix is given by
(︃ )︃ (︂
𝜕𝑓 𝜕𝑓 )︂
𝜕𝑥 𝜕𝑦 3 − 2𝑥 − 2𝑦 −2𝑥
𝜕𝑔 𝜕𝑔 = .
𝜕𝑥 𝜕𝑦
−𝑦 2 − 𝑥 − 2𝑦
Stability of the fixed points may be considered in turn. With J* the Jacobian
matrix evaluated at the fixed point, we have
(︂ )︂
3 0
(𝑥* , 𝑦* ) = (0, 0) : J* = .
0 2
86 CHAPTER 7. NONLINEAR DIFFERENTIAL EQUATIONS
1.5
y
0.5
0
0 0.5 1 1.5 2 2.5 3 3.5
x
The characteristic
√ equation of J* is given by (−1 − 𝜆)2 − 2 = 0, so that 𝜆 =
−1 ± 2. Since one eigenvalue is negative and the other positive the fixed point
(1, 1) is an unstable saddle point. From our analysis of the fixed points, one can
expect that all solutions will asymptote to one of the stable fixed points (0, 2)
or (3, 0), depending on the initial conditions.
It is of interest to sketch the phase space diagram for this nonlinear system.
The eigenvectors associated with the unstable saddle point (1, 1) determine the
directions of the flow into and away from this fixed √ point. The eigenvector
associated with the positive eigenvalue 𝜆1 = −1 + 2 can determined from the
first equation of (J* − 𝜆1 I)v1 = 0, or
√
− 2𝑣11 − 2𝑣12 = 0,
√
so that 𝑣12 = −( √2/2)𝑣11 . The eigenvector
√ associated with the negative eigen-
value 𝜆1 = −1 − 2 satisfies 𝑣22 = ( 2/2)𝑣21 . The eigenvectors give the slope
7.2. ONE-DIMENSIONAL BIFURCATIONS 87
of the lines with origin at the fixed point for incoming (negative eigenvalue) and
outgoing (positive
√ eigenvalue) trajectories. The outgoing trajectories have
√ neg-
ative slope − 2/2 and the incoming trajectories have positive slope 2/2. A
rough sketch of the phase space diagram can be made by hand (as demonstrated
in class). Here, a computer generated plot obtained from numerical solution of
the nonlinear coupled odes is presented in Fig. 7.1. The curve starting from
the origin and at infinity, and terminating at the unstable saddle point is called
the separatrix. This curve separates the phase space into two regions: initial
conditions for which the solution asymptotes to the fixed point (0, 2), and initial
conditions for which the solution asymptotes to the fixed point (3, 0).
𝑥˙ = 𝑓𝑟 (𝑥),
𝑥˙ = 𝑟 + 𝑥2 .
√
The fixed points are 𝑥* = ± −𝑟. Clearly, two real fixed points exist when
𝑟 < 0 and no real fixed points exist when 𝑟 > 0. The stability of the fixed
points when 𝑟 < 0 are determined by the derivative of 𝑓 (𝑥) = 𝑟 + 𝑥2 , given by
𝑓 ′ (𝑥) = 2𝑥. Therefore, the negative fixed point is stable and the positive fixed
point is unstable.
Graphically, we can illustrate this bifurcation in two ways. First, in Fig. 7.2(a),
we plot 𝑥˙ versus 𝑥 for the three parameter values corresponding to 𝑟 < 0, 𝑟 = 0
and 𝑟 > 0. The values at which 𝑥˙ = 0 correspond to the fixed points, and
arrows are drawn indicating how the solution 𝑥(𝑡) evolves (to the right if 𝑥˙ > 0
and to the left if 𝑥˙ < 0). The stable fixed point is indicated by a filled circle
and the unstable fixed point by an open circle. Note that when 𝑟 = 0, solutions
converge to the origin from the left, but diverge from the origin on the right.
88 CHAPTER 7. NONLINEAR DIFFERENTIAL EQUATIONS
dx/dt
(a)
→ ← → x → →
Second, in Fig. 7.2(b), we plot a bifurcation diagram illustrating the fixed point
𝑥* versus the bifurcation parameter 𝑟. The stable fixed point is denoted by a
solid line and the unstable fixed point by a dashed line. Note that the two fixed
points collide and annihilate at 𝑟 = 0, and there are no fixed points for 𝑟 > 0.
dx/dt
(a)
← → ← x ← ← ← → ←
𝑥˙ = 𝑟𝑥 − 𝑥3 .
Note that the linear term results in exponential growth when 𝑟 > 0 and the
nonlinear
√ term stabilizes this growth. The fixed points are 𝑥* = 0 and 𝑥* =
± 𝑟, the latter fixed points existing only when √ 𝑟 > 0. The derivative of 𝑓 is
𝑓 ′ (𝑥) = 𝑟 − 3𝑥2 so that 𝑓 ′ (0) = 𝑟 and 𝑓 ′ (± 𝑟) = −2𝑟. Therefore, the fixed
point 𝑥√ * = 0 is stable for 𝑟 < 0 and unstable for 𝑟 > 0 while the fixed points
𝑥 = ± 𝑟 exist and are stable for 𝑟 > 0. Notice that the fixed point 𝑥* =√0
becomes unstable as 𝑟 crosses zero and two new stable fixed points 𝑥* = ± 𝑟
are born. The supercritical pitchfork bifurcation is illustrated in Fig. 7.4.
(a) dx/dt
→ ← x → ← → ← → ←
x*
(b)
√
for 𝑟 < 0 and unstable for 𝑟 > 0 while the fixed points 𝑥 = ± −𝑟 exist and are
unstable for 𝑟 < 0. There are no stable fixed points when 𝑟 > 0.
The absence of stable fixed points for 𝑟 > 0 indicates that the neglect of
terms of higher-order in 𝑥 than 𝑥3 in the normal form may be unwarranted.
Keeping to the intrinsic symmetry of the equations (only odd powers of 𝑥) we
can add a stabilizing nonlinear term proportional to 𝑥5 . The extended normal
form (to order 𝑥5 ) is
𝑥˙ = 𝑟𝑥 + 𝑥3 − 𝑥5 ,
and is somewhat more difficult to analyze. The fixed points are solutions of
𝑥(𝑟 + 𝑥2 − 𝑥4 ) = 0.
The fixed point 𝑥* = 0 exists for all 𝑟, and four additional fixed points can be
found from the solutions of the quadratic equation in 𝑥2 :
√︂
1 √
𝑥* = ± (1 ± 1 + 4𝑟).
2
x*
points are
1
𝑟<− : 𝑥* = 0 (one fixed point);
4 √︂
1 1 √
− <𝑟<0: 𝑥* = 0, 𝑥* = ± (1 ± 1 + 4𝑟) (five fixed points);
4 2
√︂
1 √
𝑟>0: 𝑥* = 0, 𝑥* = ± (1 + 1 + 4𝑟) (three fixed points).
2
√
With 𝑥2* = 21 (1 ± 1 + 4𝑟), we have
√
(︂ )︂
1
𝑓 ′ (𝑥* ) = −2 2𝑟 + (1 ± 1 + 4𝑟)
2
(︀ √ )︀
= − (1 + 4𝑟) ± 1 + 4𝑟
√ (︀√ )︀
= − 1 + 4𝑟 1 + 4𝑟 ± 1 .
Clearly, the plus root is always stable since 𝑓 ′ (𝑥* ) < 0. The minus root exists
only for − 14 < 𝑟 < 0 and is unstable since 𝑓 ′ (𝑥* ) > 0. We summarize the
92 CHAPTER 7. NONLINEAR DIFFERENTIAL EQUATIONS
x
*
0
α 2
(1+α) /4
c
the maximum rate at which fish can be caught, and 𝐴 is a constant satisfying
𝐴 < 𝐾 that is used to model the idea that fish become harder to catch when
scarce.
We nondimensionalize (7.4) using 𝑥 = 𝑁/𝐾, 𝜏 = 𝑟𝑡, 𝑐 = 𝐶/𝑟𝐾, 𝛼 = 𝐴/𝐾:
𝑑𝑥 𝑐𝑥
= 𝑥(1 − 𝑥) − . (7.5)
𝑑𝜏 𝛼+𝑥
Note that 0 ≤ 𝑥 ≤ 1, 𝑐 > 0 and 0 < 𝛼 < 1. We wish to qualitatively describe
the equilibrium solutions of (7.5) and the bifurcations that may occur as the
nondimensional catch rate 𝑐 increases from zero. Practically, a government
would like to issue each year as large a catch quota as possible without adversely
affecting the number of fish that may be caught in subsequent years.
The fixed points of (7.5) are 𝑥* = 0, valid for all 𝑐, and the solutions to
𝑥2 − (1 − 𝛼)𝑥 + (𝑐 − 𝛼) = 0, or
1 [︁ √︀ ]︁
𝑥* = (1 − 𝛼) ± (1 + 𝛼)2 − 4𝑐 . (7.6)
2
The fixed points given by (7.6) are real only if 𝑐 < 41 (1 + 𝛼)2 . Furthermore, the
minus root is greater than zero only if 𝑐 > 𝛼. We therefore need to consider
three intervals over which the following fixed points exist:
1 [︁ √︀ ]︁
0≤𝑐≤𝛼: 𝑥* = 0, 𝑥* = (1 − 𝛼) + (1 + 𝛼)2 − 4𝑐 ;
2
1 1 [︁ √︀ ]︁
𝛼<𝑐< (1 + 𝛼)2 : 𝑥* = 0, 𝑥* = (1 − 𝛼) ± (1 + 𝛼)2 − 4𝑐 ;
4 2
1
𝑐 > (1 + 𝛼)2 : 𝑥* = 0.
4
94 CHAPTER 7. NONLINEAR DIFFERENTIAL EQUATIONS
The stability of the fixed points can be determined with rigor analytically or
graphically. Here, we simply apply biological intuition together with knowledge
of the types of one dimensional bifurcations. An intuitive argument is made
simpler if we consider 𝑐 decreasing from large values. When 𝑐 is large, that
is 𝑐 > 41 (1 + 𝛼)2 , too many fish are being caught and our intuition suggests
that the fish population goes extinct. Therefore, in this interval, the single
fixed point 𝑥* = 0 must be stable. As 𝑐 decreases, a bifurcation occurs at
𝑐 = 41 (1 + 𝛼)2 introducing two additional fixed points at 𝑥* = (1 − 𝛼)/2. The
type of one dimensional bifurcation in which two fixed points are created as a
square root becomes real is a saddlenode bifurcation, and one of the fixed points
will be stable and the other unstable. Following these fixed points as 𝑐 → 0,
we observe that the plus root goes to one, which is the appropriate stable fixed
point when there is no fishing. We therefore conclude that the plus root is stable
and the minus root is unstable. As 𝑐 decreases further from this bifurcation,
the minus root collides with the fixed point 𝑥* = 0 at 𝑐 = 𝛼. This appears to
be a transcritical bifurcation and assuming an exchange of stability occurs, we
must have the fixed point 𝑥* = 0 becoming unstable for 𝑐 < 𝛼. The resulting
bifurcation diagram is shown in Fig. 7.6.
The purpose of simple mathematical models applied to complex ecological
problems is to offer some insight. Here, we have learned that overfishing (in
the model 𝑐 > 41 (1 + 𝛼)2 ) during one year can potentially result in a sudden
collapse of the fish catch in subsequent years, so that governments need to be
particularly cautious when contemplating increases in fishing quotas.
𝑟˙ = 𝜇𝑟 − 𝑟3 ,
𝜃˙ = 𝜔 + 𝑏𝑟2 ,
7.3. TWO-DIMENSIONAL BIFURCATIONS 95
where 𝑥 = 𝑟 cos 𝜃 and 𝑦 = 𝑟 sin 𝜃. The parameter 𝜇 controls the stability of the
fixed point at the origin, the parameter 𝜔 is the frequency of oscillation near
the origin, and the parameter 𝑏 determines the dependence of the oscillation
frequency at larger amplitude oscillations. Although we include 𝑏 for generality,
our qualitative analysis of these equations will be independent of 𝑏.
The equation for the radius is of the form of the supercritical pitchfork
√
bifurcation. The fixed points are 𝑟* = 0 and 𝑟* = 𝜇 (note that 𝑟 > 0), and
the former fixed point is stable for 𝜇 < 0 and the latter is stable for 𝜇 > 0. The
transition of the eigenvalues of the Jacobian from negative real part to positive
real part can be seen if we transform these equations to cartesian coordinates.
We have using 𝑟2 = 𝑥2 + 𝑦 2 ,
˙ sin 𝜃
𝑥˙ = 𝑟˙ cos 𝜃 − 𝜃𝑟
= (𝜇𝑟 − 𝑟3 ) cos 𝜃 − (𝜔 + 𝑏𝑟2 )𝑟 sin 𝜃
= 𝜇𝑥 − (𝑥2 + 𝑦 2 )𝑥 − 𝜔𝑦 − 𝑏(𝑥2 + 𝑦 2 )𝑦
= 𝜇𝑥 − 𝜔𝑦 − (𝑥2 + 𝑦 2 )(𝑥 + 𝑏𝑦);
˙ cos 𝜃
𝑦˙ = 𝑟˙ sin 𝜃 + 𝜃𝑟
= (𝜇𝑟 − 𝑟3 ) sin 𝜃 + (𝜔 + 𝑏𝑟2 )𝑟 cos 𝜃
= 𝜇𝑦 − (𝑥2 + 𝑦 2 )𝑦 + 𝜔𝑥 + 𝑏(𝑥2 + 𝑦 2 )𝑥
= 𝜔𝑥 + 𝜇𝑦 − (𝑥2 + 𝑦 2 )(𝑦 − 𝑏𝑥).
𝑟˙ = 𝜇𝑟 + 𝑟3 − 𝑟5 ,
𝜃˙ = 𝜔 + 𝑏𝑟2 .
Here, the equation for the radius is of the form of the subcritical pitchfork
bifurcation. As 𝜇 increases from negative to positive values, the origin becomes
unstable and exponentially growing oscillations increase until the radius reaches
a stable fixed point far away from the origin. In practice, it may be difficult
to tell analytically if a Hopf bifurcation is supercritical or subcritical from the
equations of motion. Computational solution, however, can quickly distinguish
between the two types.
96 CHAPTER 7. NONLINEAR DIFFERENTIAL EQUATIONS
Chapter 8
Partial differential
equations
Reference: Boyce and DiPrima, Chapter 10
97
98 CHAPTER 8. PARTIAL DIFFERENTIAL EQUATIONS
J (x2)
J (x1)
X1 X2
by
𝐽 = −𝐷𝑢𝑥 , (8.1)
where the diffusion constant 𝐷 > 0 has units [𝐷] = 𝑙2 /𝑡, and we have used
the notation 𝑢𝑥 = 𝜕𝑢/𝜕𝑥. The mass flux is opposite in sign to the gradient of
concentration. The time rate of change in the mass of dye between 𝑥1 and 𝑥2
is given by the difference between the mass flux into and the mass flux out of
the infinitesimal cross sectional volume. When 𝑢𝑥 < 0, 𝐽 > 0 and the mass
flows into the volume at position 𝑥1 and out of the volume at position 𝑥2 . On
the other hand, when 𝑢𝑥 > 0, 𝐽 < 0 and the mass flows out of the volume at
position 𝑥1 and into the volume at position 𝑥2 . In both cases, the time rate of
change of the dye mass is given by
𝑑𝑀
= 𝐽(𝑥1 , 𝑡) − 𝐽(𝑥2 , 𝑡),
𝑑𝑡
or rewriting in terms of 𝑢(𝑥, 𝑡):
𝑢𝑡 = 𝐷𝑢𝑥𝑥 .
We note that the diffusion equation is identical to the heat conduction equation,
where 𝑢 is temperature, and the constant 𝐷 (commonly written as 𝜅) is the
thermal conductivity.
constant mass density of the string, 𝑇 the tension of the string, and 𝜃 the angle
between the string and the horizonal line. We consider an infinitesimal string
element located between 𝑥1 and 𝑥2 , with Δ𝑥 = 𝑥2 − 𝑥1 , as shown in Fig. 8.2.
The governing equations are Newton’s law of motion for the horizontal and
vertical acceleration of our infinitesimal string element, and we assume that the
string element only accelerates vertically. Therefore, the horizontal forces must
balance and we have
𝑇2 cos 𝜃2 = 𝑇1 cos 𝜃1 .
The vertical forces result in a vertical √
acceleration, and with√︀ 𝑢𝑡𝑡 the vertical
acceleration of the string element and 𝜌 Δ𝑥2 + Δ𝑢2 = 𝜌Δ𝑥 1 + 𝑢2𝑥 its mass,
where we have used 𝑢𝑥 = Δ𝑢/Δ𝑥, exact as Δ𝑥 → 0, we have
√︀
𝜌Δ𝑥 1 + 𝑢2𝑥 𝑢𝑡𝑡 = 𝑇2 sin 𝜃2 − 𝑇1 sin 𝜃1 .
cos 𝜃2 = cos 𝜃1 = 1,
𝑢𝑡𝑡 = 𝑐2 𝑢𝑥𝑥 ,
where 𝑐2 = 𝑇 /𝜌. Since [𝑇 ] = 𝑚𝑙/𝑡2 and [𝜌] = 𝑚/𝑙, we have [𝑐2 ] = 𝑙2 /𝑡2 so that
𝑐 has units of velocity.
The orthogonality relations for 𝑛 and 𝑚 positive integers are then given with
compact notation as the integration formulas
∫︁ 𝐿 (︁ 𝑚𝜋𝑥 )︁ (︁ 𝑛𝜋𝑥 )︁
cos cos 𝑑𝑥 = 𝐿𝛿𝑛𝑚 , (8.3)
−𝐿 𝐿 𝐿
∫︁ 𝐿 (︁ 𝑚𝜋𝑥 )︁ (︁ 𝑛𝜋𝑥 )︁
sin sin 𝑑𝑥 = 𝐿𝛿𝑛𝑚 , (8.4)
−𝐿 𝐿 𝐿
∫︁ 𝐿 (︁ 𝑚𝜋𝑥 )︁ (︁ 𝑛𝜋𝑥 )︁
cos sin 𝑑𝑥 = 0. (8.5)
−𝐿 𝐿 𝐿
For 𝑚 = 𝑛, however,
𝐿
𝐿 𝜋
∫︁ (︁ 𝑛𝜋𝑥 )︁ ∫︁
2
sin 𝑑𝑥 = sin2 (𝑛𝜉)𝑑𝜉
−𝐿 𝐿 𝜋 −𝜋
𝐿 𝜋 (︀
∫︁
)︀
= 1 − cos (2𝑛𝜉) 𝑑𝜉
2𝜋 −𝜋
[︂ ]︂𝜋
𝐿 1
= 𝜉− sin 2𝑛𝜉
2𝜋 2𝑛 −𝜋
= 𝐿.
𝐿
𝑎0 𝐿
∫︁ ∫︁
𝑛𝜋𝑥 𝑛𝜋𝑥
𝑓 (𝑥) cos 𝑑𝑥 = cos 𝑑𝑥
−𝐿 𝐿 2 −𝐿 𝐿
∞
{︃ ∫︁ }︃
𝐿 ∫︁ 𝐿
∑︁ 𝑛𝜋𝑥 𝑚𝜋𝑥 𝑛𝜋𝑥 𝑚𝜋𝑥
+ 𝑎𝑚 cos cos 𝑑𝑥 + 𝑏𝑚 cos sin 𝑑𝑥 .
𝑚=1 −𝐿 𝐿 𝐿 −𝐿 𝐿 𝐿
If 𝑛 = 0, then the second and third integrals on the right-hand-side are zero
and the first integral is 2𝐿 so that the right-hand-side becomes 𝐿𝑎0 . If 𝑛 is
a positive integer, then the first and third integrals on the right-hand-side are
zero, and the second integral is 𝐿𝛿𝑛𝑚 . For this case, we have
∫︁ 𝐿 ∞
𝑛𝜋𝑥 ∑︁
𝑓 (𝑥) cos 𝑑𝑥 = 𝐿𝑎𝑚 𝛿𝑛𝑚
−𝐿 𝐿 𝑚=1
= 𝐿𝑎𝑛 ,
where all the terms in the summation except 𝑚 = 𝑛 are zero by virtue of the
Kronecker delta. We therefore obtain for 𝑛 = 0, 1, 2, . . .
∫︁ 𝐿
1 𝑛𝜋𝑥
𝑎𝑛 = 𝑓 (𝑥) cos 𝑑𝑥. (8.6)
𝐿 −𝐿 𝐿
𝐿
𝑎0 𝐿
∫︁ ∫︁
𝑛𝜋𝑥 𝑛𝜋𝑥
𝑓 (𝑥) sin 𝑑𝑥 = sin 𝑑𝑥
−𝐿 𝐿 2 −𝐿 𝐿
∞
{︃ ∫︁ 𝐿 ∫︁ 𝐿 }︃
∑︁ 𝑛𝜋𝑥 𝑚𝜋𝑥 𝑛𝜋𝑥 𝑚𝜋𝑥
+ 𝑎𝑚 sin cos 𝑑𝑥 + 𝑏𝑚 sin sin 𝑑𝑥 .
𝑚=1 −𝐿 𝐿 𝐿 −𝐿 𝐿 𝐿
102 CHAPTER 8. PARTIAL DIFFERENTIAL EQUATIONS
Here, the first and second integrals on the right-hand-side are zero, and the
third integral is 𝐿𝛿𝑛𝑚 so that
∫︁ 𝐿 ∞
𝑛𝜋𝑥 ∑︁
𝑓 (𝑥) sin 𝑑𝑥 = 𝐿𝑏𝑚 𝛿𝑛𝑚
−𝐿 𝐿 𝑚=1
= 𝐿𝑏𝑛 .
Hence, for 𝑛 = 1, 2, 3, . . . ,
∫︁ 𝐿
1 𝑛𝜋𝑥
𝑏𝑛 = 𝑓 (𝑥) sin 𝑑𝑥. (8.7)
𝐿 −𝐿 𝐿
Our results for the Fourier series of a function 𝑓 (𝑥) with period 2𝐿 are thus
given by (8.2), (8.6) and (8.7).
We examine in turn the Fourier series for an even or an odd function. First,
if 𝑓 (𝑥) is even, then from (8.6) and (8.7) and our facts about even and odd
functions,
2 𝐿
∫︁
𝑛𝜋𝑥
𝑎𝑛 = 𝑓 (𝑥) cos 𝑑𝑥,
𝐿 0 𝐿 (8.8)
𝑏𝑛 = 0.
The Fourier series for an even function with period 2𝐿 is thus given by the
Fourier cosine series
∞
𝑎0 ∑︁ 𝑛𝜋𝑥
𝑓 (𝑥) = + 𝑎𝑛 cos , 𝑓 (𝑥) even. (8.9)
2 𝑛=1
𝐿
−1
−2*pi −pi 0 pi 2*pi
and the Fourier series for an odd function with period 2𝐿 is given by the Fourier
sine series
∞
∑︁ 𝑛𝜋𝑥
𝑓 (𝑥) = 𝑏𝑛 sin , 𝑓 (𝑥) odd. (8.11)
𝑛=1
𝐿
2𝑥
𝑓 (𝑥) = 1 − .
𝜋
Because 𝑓 (𝑥) is even, it can be represented by the Fourier cosine series given by
(8.8) and (8.9). The coefficient 𝑎0 is
2 𝜋
∫︁
𝑎0 = 𝑓 (𝑥)𝑑𝑥
𝜋 0
2 𝜋
∫︁ (︂ )︂
2𝑥
= 1− 𝑑𝑥
𝜋 0 𝜋
]︂𝜋
𝑥2
[︂
2
= 𝑥−
𝜋 𝜋 0
= 0.
104 CHAPTER 8. PARTIAL DIFFERENTIAL EQUATIONS
𝑋𝑇 ′ = 𝐷𝑋 ′′ 𝑇,
𝑋 ′′ + 𝜆𝑋 = 0, 𝑇 ′ + 𝜆𝐷𝑇 = 0. (8.16)
Because of the boundary conditions, we must first consider the equation for
𝑋(𝑥). To solve, we need to determine the boundary conditions at 𝑥 = 0 and
𝑥 = 𝐿. Now, from (8.14) and (8.15),
Since 𝑇 (𝑡) is not identically zero for all 𝑡 (which would result in the trivial
solution for 𝑢), we must have 𝑋(0) = 0. Similarly, the boundary condition at
𝑥 = 𝐿 requires 𝑋(𝐿) = 0. We therefore consider the two-point boundary value
problem
𝑋 ′′ + 𝜆𝑋 = 0, 𝑋(0) = 𝑋(𝐿) = 0. (8.17)
106 CHAPTER 8. PARTIAL DIFFERENTIAL EQUATIONS
The equation given by (8.17) is called an ode eigenvalue problem. The allowed
values of 𝜆 and the corresponding functions 𝑋(𝑥) are called the eigenvalues and
eigenfunctions of the differential equation. Since the form of the general solution
of the ode depends on the sign of 𝜆, we consider in turn the cases 𝜆 > 0, 𝜆 < 0
and 𝜆 = 0. For 𝜆 > 0, we write 𝜆 = 𝜇2 and determine the general solution of
𝑋 ′′ + 𝜇2 𝑋 = 0
to be
𝑋(𝑥) = 𝐴 cos 𝜇𝑥 + 𝐵 sin 𝜇𝑥.
Applying the boundary condition at 𝑥 = 0, we find 𝐴 = 0. The boundary
condition at 𝑥 = 𝐿 then yields
𝐵 sin 𝜇𝐿 = 0.
The solution 𝐵 = 0 results in the trivial solution for 𝑢 and can be ruled out.
Therefore, we must have
sin 𝜇𝐿 = 0,
which is an equation for the eigenvalue 𝜇. The solutions are
𝜇 = 𝑛𝜋/𝐿,
𝑋 ′′ − 𝜇2 𝑋 = 0
to be
𝑋(𝑥) = 𝐴 cosh 𝜇𝑥 + 𝐵 sinh 𝜇𝑥,
where we have previously introduced the hyperbolic sine and cosine functions
in S3.4.1. Applying the boundary condition at 𝑥 = 0, we find 𝐴 = 0. The
boundary condition at 𝑥 = 𝐿 then yields
𝐵 sinh 𝜇𝐿 = 0,
𝑋 ′′ = 0,
We now turn to the equation for 𝑇 (𝑡). The equation corresponding to the
eigenvalue 𝜆𝑛 , using (8.18), is given by
𝑇 ′ + 𝑛2 𝜋 2 𝐷/𝐿2 𝑇 = 0,
(︀ )︀
satisfy the pde given by (8.12) and the boundary conditions given by (8.14) for
every positive integer 𝑛.
The principle of linear superposition for homogeneous linear differential
equations then states that the general solution to (8.12) and (8.14) is given
by
∞
∑︁
𝑢(𝑥, 𝑡) = 𝑏𝑛 𝑢𝑛 (𝑥, 𝑡)
𝑛=1
∞
(8.22)
−𝑛2 𝜋 2 𝐷𝑡/𝐿2
∑︁
= 𝑏𝑛 sin (𝑛𝜋𝑥/𝐿)𝑒 .
𝑛=1
The final solution step is to satisfy the initial conditions given by (8.13). At
𝑡 = 0, we have
∞
∑︁
𝑓 (𝑥) = 𝑏𝑛 sin (𝑛𝜋𝑥/𝐿). (8.23)
𝑛=1
We immediately recognize (8.23) as a Fourier sine series (8.11) for an odd func-
tion 𝑓 (𝑥) with period 2𝐿. Equation (8.23) is a periodic extension of our original
𝑓 (𝑥) defined on 0 ≤ 𝑥 ≤ 𝐿, and is an odd function because of the boundary
condition 𝑓 (0) = 0. From our solution for the coefficients of a Fourier sine series
(8.10), we determine
2 𝐿
∫︁
𝑛𝜋𝑥
𝑏𝑛 = 𝑓 (𝑥) sin 𝑑𝑥. (8.24)
𝐿 0 𝐿
Thus the solution to the diffusion equation with homogeneous Dirichlet bound-
ary conditions defined by (8.12), (8.13) and (8.14) is given by (8.22) with the
𝑏𝑛 coefficients computed from (8.24).
2 𝐿
∫︁
𝐿 𝑛𝜋𝑥
𝑏𝑛 = 𝛿(𝑥 − ) sin 𝑑𝑥
𝐿 0 2 𝐿
2
= sin (𝑛𝜋/2)
𝐿
⎧
⎨ 2/𝐿 if 𝑛 = 1, 5, 9, . . . ;
= −2/𝐿 if 𝑛 = 3, 7, 11, . . . ;
0 if 𝑛 = 2, 4, 6, . . . .
⎩
With 𝑏𝑛 determined, the solution for 𝑢(𝑥, 𝑡) given by (8.22) can be written as
∞ (︂ )︂
2 ∑︁ 𝑛 (2𝑛 + 1)𝜋𝑥 −(2𝑛+1)2 𝜋2 𝐷𝑡/𝐿2
𝑢(𝑥, 𝑡) = (−1) sin 𝑒 .
𝐿 𝑛=0 𝐿
𝑢𝑡 = 𝐷𝑢𝑥𝑥 .
so that
𝑋(0) = 𝐶1 /𝑇 (𝑡).
However, our separation of variables ansatz assumes 𝑋(𝑥) to be independent of
𝑡! We therefore say that inhomogeneous boundary conditions are not separable.
The proper way to solve a problem with inhomogeneous boundary conditions
is to transform it into another problem with homogeneous boundary conditions.
As 𝑡 → ∞, we assume that a stationary concentration distribution 𝑣(𝑥) will
attain, independent of 𝑡. Since 𝑣(𝑥) must satisfy the diffusion equation, we have
𝑣 ′′ (𝑥) = 0, 0 ≤ 𝑥 ≤ 𝐿,
Since 𝑣(𝑥) must satisfy the same boundary conditions of 𝑢(𝑥, 𝑡), we have 𝑣(0) =
𝐶1 and 𝑣(𝐿) = 𝐶2 , and we determine 𝐴 = 𝐶1 and 𝐵 = (𝐶2 − 𝐶1 )/𝐿.
We now express 𝑢(𝑥, 𝑡) as the sum of the known asymptotic stationary con-
centration distribution 𝑣(𝑥) and an unknown transient concentration distribu-
tion 𝑤(𝑥, 𝑡):
𝑢(𝑥, 𝑡) = 𝑣(𝑥) + 𝑤(𝑥, 𝑡).
Substituting into the diffusion equation, we obtain
𝜕 𝜕2
(𝑣(𝑥) + 𝑤(𝑥, 𝑡)) = 𝐷 2 (𝑣(𝑥) + 𝑤(𝑥, 𝑡))
𝜕𝑡 𝜕𝑥
or
𝑤𝑡 = 𝐷𝑤𝑥𝑥 ,
The resulting equations may then be solved for 𝑤(𝑥, 𝑡) using the technique for
homogeneous boundary conditions, and 𝑢(𝑥, 𝑡) subsequently determined.
Again, we consider in turn the cases 𝜆 > 0, 𝜆 < 0 and 𝜆 = 0. For 𝜆 > 0, we
write 𝜆 = 𝜇2 and determine the general solution of (8.26) to be
−𝜇𝐴 sin 𝜇𝐿 = 0.
The solution 𝐴 = 0 results in the trivial solution for 𝑢 and can be ruled out.
Therefore, we must have
sin 𝜇𝐿 = 0,
with solutions
𝜇 = 𝑛𝜋/𝐿,
where 𝑛 is an integer. We have thus determined the eigenvalues 𝜆 = 𝜇2 > 0 to
be
2
𝜆𝑛 = (𝑛𝜋/𝐿) , 𝑛 = 1, 2, 3, . . . , (8.27)
with corresponding eigenfunctions
For 𝜆 < 0, we write 𝜆 = −𝜇2 and determine the general solution of (8.26) to be
𝜇𝐴 sinh 𝜇𝐿 = 0,
𝑋(𝑥) = 𝐴 + 𝐵𝑥,
𝑇 ′ + 𝑛2 𝜋 2 𝐷/𝐿2 𝑇 = 0,
(︀ )︀
satisfy the pde given by (8.12) and the boundary conditions given by (8.25) for
every nonnegative integer 𝑛.
The principle of linear superposition then yields the general solution as
∞
∑︁
𝑢(𝑥, 𝑡) = 𝑐𝑛 𝑢𝑛 (𝑥, 𝑡)
𝑛=0
∞
(8.31)
𝑎0 ∑︁ 2 2 2
= + 𝑎𝑛 cos (𝑛𝜋𝑥/𝐿)𝑒−𝑛 𝜋 𝐷𝑡/𝐿 ,
2 𝑛=1
which we recognize as a Fourier cosine series (8.9) for an even function 𝑓 (𝑥)
with period 2𝐿. We have obtained a cosine series for the periodic extension of
𝑓 (𝑥) because of the boundary condition 𝑓 ′ (0) = 0, which is satisfied by an even
function. From our solution (8.8) for the coefficients of a Fourier cosine series,
we determine
2 𝐿
∫︁
𝑛𝜋𝑥
𝑎𝑛 = 𝑓 (𝑥) cos 𝑑𝑥. (8.33)
𝐿 0 𝐿
Thus the solution to the diffusion equation with homogenous Neumann bound-
ary conditions defined by (8.12), (8.13) and (8.25) is given by (8.31) with the
coefficients computed from (8.33).
2 𝐿
∫︁
𝐿 𝑛𝜋𝑥
𝑎𝑛 = 𝛿(𝑥 − ) cos 𝑑𝑥
𝐿 0 2 𝐿
2
= cos (𝑛𝜋/2)
⎧𝐿
⎨ 2/𝐿 if 𝑛 = 0, 4, 8, . . . ;
= −2/𝐿 if 𝑛 = 2, 6, 10, . . . ;
0 if 𝑛 = 1, 3, 5, . . . .
⎩
The first two terms in the series for 𝑢(𝑥, 𝑡) are given by
2 [︁ 2 2
]︁
𝑢(𝑥, 𝑡) = 1/2 − cos (2𝜋𝑥/𝐿)𝑒−4𝜋 𝐷𝑡/𝐿 + . . . .
𝐿
112 CHAPTER 8. PARTIAL DIFFERENTIAL EQUATIONS
Notice that as 𝑡 → ∞, 𝑢(𝑥, 𝑡) → 1/𝐿: the dye mass is conserved in the pipe
(since the pipe ends are sealed) and eventually diffused uniformly throughout
the pipe of length 𝐿.
Again we use the method of separation of variables and try the ansatz
Substitution of our ansatz (8.37) into the wave equation (8.34) and separating
variables results in
𝑋 ′′ 1 𝑇 ′′
= 2 = −𝜆,
𝑋 𝑐 𝑇
yielding the two ordinary differential equations
𝑋 ′′ + 𝜆𝑋 = 0, 𝑇 ′′ + 𝜆𝑐2 𝑇 = 0. (8.38)
We solve first the equation for 𝑋(𝑥). The appropriate boundary conditions for
𝑋 are given by
𝑋(0) = 0, 𝑋(𝐿) = 0, (8.39)
and we have solved this equation for 𝑋(𝑥) previously in S8.5.1 (see (8.17)).
A nontrivial solution exists only when 𝜆 > 0, and our previously determined
solution was
2
𝜆𝑛 = (𝑛𝜋/𝐿) , 𝑛 = 1, 2, 3, . . . , (8.40)
with corresponding eigenfunctions
The remaining condition to satisfy is the initial displacement of the string, the
first equation of (8.36). We have
∞
∑︁
𝑓 (𝑥) = 𝑏𝑛 sin (𝑛𝜋𝑥/𝐿),
𝑛=1
which is observed to be a Fourier sine series (8.11) for an odd function with
period 2𝐿. Therefore, the coefficients 𝑏𝑛 are given by (8.10),
2 𝐿
∫︁
𝑛𝜋𝑥
𝑏𝑛 = 𝑓 (𝑥) sin 𝑑𝑥, 𝑛 = 1, 2, 3, . . . . (8.44)
𝐿 0 𝐿
Our solution to the wave equation with plucked string is thus given by (8.43)
and (8.44). Notice that the solution is time periodic with period 2𝐿/𝑐. The
corresponding fundamental frequency is the reciprocal of the period and is given
by 𝑓 = 𝑐/2𝐿. From our derivation of the wave equation in S8.2, the velocity 𝑐
is related to the density of the string 𝜌 and tension of the string 𝑇 by 𝑐2 = 𝑇 /𝜌.
Therefore, the fundamental frequency (pitch) of our “guitar string” increases (is
raised) with increasing tension, decreasing string density, and decreasing string
length. Indeed, these are exactly the parameters used to construct, tune and
play a guitar.
The wave nature of our solution and the physical significance of the velocity
𝑐 can be made more transparent if we make use of the trigonometric identity
1 (︀ )︀
sin 𝑥 cos 𝑦 = sin (𝑥 + 𝑦) + sin (𝑥 − 𝑦) .
2
With this identity, our solution (8.43) can be rewritten as
∞ (︂ )︂
1 ∑︁ 𝑛𝜋(𝑥 + 𝑐𝑡) 𝑛𝜋(𝑥 − 𝑐𝑡)
𝑢(𝑥, 𝑡) = 𝑏𝑛 sin + sin . (8.45)
2 𝑛=1 𝐿 𝐿
114 CHAPTER 8. PARTIAL DIFFERENTIAL EQUATIONS
The first and second sine functions can be interpreted as a traveling wave moving
to the left or the right with velocity 𝑐. This can be seen by incrementing time,
𝑡 → 𝑡 + 𝛿, and observing that the value of the first sine function is unchanged
provided the position is shifted by 𝑥 → 𝑥 − 𝑐𝛿, and the second sine function is
unchanged provided 𝑥 → 𝑥 + 𝑐𝛿. Two waves travelling in opposite directions
with equal amplitude results in a standing wave.
Our solution proceeds as previously, except that now the homogeneous initial
condition on 𝑇 (𝑡) is 𝑇 (0) = 0, so that 𝐴 = 0 in (8.42). The wave equation
solution is therefore
∞
∑︁ 𝑛𝜋𝑥 𝑛𝜋𝑐𝑡
𝑢(𝑥, 𝑡) = 𝑏𝑛 sin sin . (8.47)
𝑛=1
𝐿 𝐿
The coefficient of the Fourier sine series for 𝑔(𝑥) is seen to be 𝑛𝜋𝑐𝑏𝑛 /𝐿, and we
have
2 𝐿
∫︁
𝑛𝜋𝑐𝑏𝑛 𝑛𝜋𝑥
= 𝑔(𝑥) sin 𝑑𝑥,
𝐿 𝐿 0 𝐿
or ∫︁ 𝐿
2 𝑛𝜋𝑥
𝑏𝑛 = 𝑔(𝑥) sin 𝑑𝑥.
𝑛𝜋𝑐 0 𝐿
then the solution to the wave equation can be determined using the principle of
linear superposition. Suppose 𝑣(𝑥, 𝑡) is the solution to the wave equation with
initial condition (8.36) and 𝑤(𝑥, 𝑡) is the solution to the wave equation with
initial conditions (8.46). Then we have
since 𝑢(𝑥, 𝑡) satisfies the wave equation, the boundary conditions, and the initial
conditions given by (8.48).
8.7. THE LAPLACE EQUATION 115
𝑋 ′′ 𝑌 ′′
=− = 𝜆,
𝑋 𝑌
with 𝜆 the separation constant. We thus obtain the two ordinary differential
equations
𝑋 ′′ − 𝜆𝑋 = 0, 𝑌 ′′ + 𝜆𝑌 = 0.
The homogeneous boundary conditions are 𝑋(0) = 0, 𝑌 (0) = 0 and 𝑌 (𝑏) = 0.
We have already solved the equation for 𝑌 (𝑦) in S8.5.1, and the solution yields
the eigenvalues
(︁ 𝑛𝜋 )︁2
𝜆𝑛 = , 𝑛 = 1, 2, 3, . . . ,
𝑏
with corresponding eigenfunctions
𝑛𝜋𝑦
𝑌𝑛 (𝑦) = sin .
𝑏
The remaining 𝑋 equation and homogeneous boundary condition is therefore
𝑛2 𝜋 2
𝑋 ′′ − 𝑋 = 0, 𝑋(0) = 0,
𝑏2
and the solution is the hyperbolic sine function
𝑛𝜋𝑥
𝑋𝑛 (𝑥) = sinh ,
𝑏
times a constant. Writing 𝑢𝑛 = 𝑋𝑛 𝑌𝑛 , multiplying by a constant and summing
over 𝑛, yields the general solution
∞
∑︁ 𝑛𝜋𝑥 𝑛𝜋𝑦
𝑢(𝑥, 𝑡) = 𝑐𝑛 sinh sin .
𝑛=0
𝑏 𝑏
which we recognize as a Fourier sine series for an odd function with period 2𝑏,
and coefficient 𝑐𝑛 sinh (𝑛𝜋𝑎/𝑏). The solution for the coefficient is given by
∫︁ 𝑏
𝑛𝜋𝑎 2 𝑛𝜋𝑦
𝑐𝑛 sinh = 𝑓 (𝑦) sin 𝑑𝑥.
𝑏 𝑏 0 𝑏
∇2 𝑢 = 0, (8.50)
8.7. THE LAPLACE EQUATION 117
After taking the partial derivatives of 𝑥 and 𝑦 using (8.51), we can write the
transformation (8.52) in matrix form as
(︂ )︂ (︂ )︂ (︂ )︂
𝜕𝑢/𝜕𝑟 cos 𝜃 sin 𝜃 𝜕𝑢/𝜕𝑥
= . (8.53)
𝜕𝑢/𝜕𝜃 −𝑟 sin 𝜃 𝑟 cos 𝜃 𝜕𝑢/𝜕𝑦
then (︂ )︂
1 𝑑 −𝑏
A−1 = .
det A −𝑐 𝑎
Therefore, since the determinant of the 2 × 2 matrix in (8.53) is 𝑟, we have
(︂ )︂ (︂ )︂ (︂ )︂
𝜕𝑢/𝜕𝑥 cos 𝜃 − sin 𝜃/𝑟 𝜕𝑢/𝜕𝑟
= . (8.54)
𝜕𝑢/𝜕𝑦 sin 𝜃 cos 𝜃/𝑟 𝜕𝑢/𝜕𝜃
𝜕 𝜕 sin 𝜃 𝜕 𝜕 𝜕 cos 𝜃 𝜕
= cos 𝜃 − , = sin 𝜃 + . (8.55)
𝜕𝑥 𝜕𝑟 𝑟 𝜕𝜃 𝜕𝑦 𝜕𝑟 𝑟 𝜕𝜃
To find the Laplacian in polar coordinates with minimum algebra, we combine
(8.55) using complex variables as
(︂ )︂
𝜕 𝜕 𝑖𝜃 𝜕 𝑖 𝜕
+𝑖 =𝑒 + , (8.56)
𝜕𝑥 𝜕𝑦 𝜕𝑟 𝑟 𝜕𝜃
so that the Laplacian may be found by multiplying both sides of (8.56) by its
complex conjugate, taking care with the computation of the derivatives on the
118 CHAPTER 8. PARTIAL DIFFERENTIAL EQUATIONS
right-hand-side:
𝜕2 𝜕2
(︂ )︂ (︂ )︂
𝑖𝜃 𝜕 𝑖 𝜕 −𝑖𝜃 𝜕 𝑖 𝜕
+ = 𝑒 + 𝑒 −
𝜕𝑥2 𝜕𝑦 2 𝜕𝑟 𝑟 𝜕𝜃 𝜕𝑟 𝑟 𝜕𝜃
𝜕2 1 𝜕 1 𝜕2
= 2+ + .
𝜕𝑟 𝑟 𝜕𝑟 𝑟2 𝜕𝜃2
We have therefore determined that the Laplacian in polar coordinates is given
by
𝜕2 1 𝜕 1 𝜕2
∇2 = 2 + + 2 2. (8.57)
𝜕𝑟 𝑟 𝜕𝑟 𝑟 𝜕𝜃
We now consider the solution of the Laplace equation in a circle with radius
𝑟 < 𝑎 subject to the boundary condition
𝑢(𝑟, 𝜃) = 𝑅(𝑟)Θ(𝜃),
and substitution into the Laplace equation (8.50) using (8.57) yields
1 1
𝑅′′ Θ + 𝑅′ Θ + 2 𝑅Θ′′ = 0,
𝑟 𝑟
or
𝑅′′ 𝑅′ Θ′′
𝑟2 +𝑟 =− = 𝜆,
𝑅 𝑅 Θ
where 𝜆 is the separation constant. We thus obtain the two ordinary differential
equations
𝑟2 𝑅′′ + 𝑟𝑅′ − 𝜆𝑅 = 0, Θ′′ + 𝜆Θ = 0.
The Θ equation is solved assuming periodic boundary conditions with period
2𝜋. If 𝜆 < 0, then no periodic solution exists. If 𝜆 = 0, then Θ can be constant.
If 𝜆 = 𝜇2 > 0, then
Θ(𝜃) = 𝐴 cos 𝜇𝜃 + 𝐵 sin 𝜇𝜃,
and the requirement that Θ is periodic with period 2𝜋 forces 𝜇 to be an integer.
Therefore,
𝜆𝑛 = 𝑛2 , 𝑛 = 0, 1, 2, . . . ,
with corresponding eigenfunctions
there are two real solutions when 𝑛 > 0 and degenerate solutions when 𝑛 = 0.
When 𝑛 > 0, the solution for R(r) is
The requirement that 𝑢(𝑟, 𝜃) is finite in the circle forces 𝐵 = 0 since 𝑟−𝑛 becomes
unbounded as 𝑟 → 0. When 𝑛 = 0, the solution for 𝑅(𝑟) is
𝑅𝑛 (𝑟) = 𝐴 + 𝐵 ln 𝑟,
and again finite 𝑢 in the circle forces 𝐵 = 0. Therefore, the solution for 𝑛 =
0, 1, 2, . . . is 𝑅𝑛 = 𝑟𝑛 . Thus the general solution for 𝑢(𝑟, 𝜃) may be written as
∞
𝐴0 ∑︁ 𝑛
𝑢(𝑟, 𝜃) = + 𝑟 (𝐴𝑛 cos 𝑛𝜃 + 𝐵𝑛 sin 𝑛𝜃), (8.60)
2 𝑛=1
where we have separated out the 𝑛 = 0 solution to write our solution in a form
similar to the standard Fourier series given by (8.2). The remaining boundary
condition (8.58) specifies the values of 𝑢 on the circle of radius 𝑎, and imposition
of this boundary condition results in
∞
𝐴0 ∑︁ 𝑛
𝑓 (𝜃) = + 𝑎 (𝐴𝑛 cos 𝑛𝜃 + 𝐵𝑛 sin 𝑛𝜃). (8.61)
2 𝑛=1
Equation (8.61) is a Fourier series for the periodic function 𝑓 (𝜃) with period 2𝜋,
i.e., 𝐿 = 𝜋 in (8.2). The Fourier coefficients 𝑎𝑛 𝐴𝑛 and 𝑎𝑛 𝐵𝑛 are therefore given
by (8.6) and (8.7) to be
1 2𝜋
∫︁
𝑎𝑛 𝐴𝑛 = 𝑓 (𝜑) cos 𝑛𝜑𝑑𝜑, 𝑛 = 0, 1, 2, . . . ;
𝜋 0
1 2𝜋
∫︁
𝑎𝑛 𝐵𝑛 = 𝑓 (𝜑) sin 𝑛𝜑𝑑𝜑, 𝑛 = 1, 2, 3, . . . . (8.62)
𝜋 0
A remarkable fact is that the infinite series solution for 𝑢(𝑟, 𝜃) can be summed
explicitly. Substituting (8.62) into (8.60), we obtain
∞ (︁ )︁𝑛
∫︁ 2𝜋 [︃ ]︃
1 ∑︁ 𝑟
𝑢(𝑟, 𝜃) = 𝑑𝜑𝑓 (𝜑) 1 + 2 (cos 𝑛𝜃 cos 𝑛𝜑 + sin 𝑛𝜃 sin 𝑛𝜑)
2𝜋 0 𝑛=1
𝑎
∞ (︁ )︁𝑛
∫︁ 2𝜋 [︃ ]︃
1 ∑︁ 𝑟
= 𝑑𝜑𝑓 (𝜑) 1 + 2 cos 𝑛(𝜃 − 𝜑) .
2𝜋 0 𝑛=1
𝑎
We can sum the infinite series by writing 2∑︀cos 𝑛(𝜃 − 𝜑) = 𝑒𝑖𝑛(𝜃−𝜑) + 𝑒−𝑖𝑛(𝜃−𝜑) ,
∞
and using the sum of the geometric series 𝑛=1 𝑧 𝑛 = 𝑧/(1 − 𝑧) to obtain
∞ (︁ )︁𝑛 ∞ (︂ )︂𝑛 ∑︁ ∞ (︂ )︂𝑛
∑︁ 𝑟 ∑︁ 𝑟𝑒𝑖(𝜃−𝜑) 𝑟𝑒−𝑖(𝜃−𝜑)
1+2 cos 𝑛(𝜃 − 𝜑) = 1 + +
𝑛=1
𝑎 𝑛=1
𝑎 𝑛=1
𝑎
𝑟𝑒𝑖(𝜃−𝜑)
(︂ )︂
=1+ + c.c.
𝑎 − 𝑟𝑒𝑖(𝜃−𝜑)
𝑎2 − 𝑟2
= 2 .
𝑎 − 2𝑎𝑟 cos (𝜃 − 𝜑) + 𝑟2
120 CHAPTER 8. PARTIAL DIFFERENTIAL EQUATIONS
Therefore,
2𝜋
𝑎2 − 𝑟2
∫︁
𝑓 (𝜑)
𝑢(𝑟, 𝜃) = 𝑑𝜑,
2𝜋 0 𝑎2 − 2𝑎𝑟 cos (𝜃 − 𝜑) + 𝑟2