IA Differential Equations (Full)

Part IA — Differential Equations

Based on lectures by M. G. Worster

Notes taken by Dexter Chua

Michaelmas 2014

These notes are not endorsed by the lecturers, and I have modified them (often

significantly) after lectures. They are nowhere near accurate representations of what

was actually lectured, and in particular, all errors are almost surely mine.

Basic calculus

Informal treatment of differentiation as a limit, the chain rule, Leibnitz’s rule, Taylor

series, informal treatment of

and

notation and l’Hˆopital’s rule; integration as an

area, fundamental theorem of calculus, integration by substitution and parts. [3]

Informal treatment of partial derivatives, geometrical interpretation, statement (only)

of symmetry of mixed partial derivatives, chain rule, implicit differentiation. Informal

treatment of differentials, including exact differentials. Differentiation of an integral

with resp ect to a parameter. [2]

First-order linear differential equations

Equations with constant coefficients: exp onential growth, comparison with discrete

equations, series solution; modelling examples including radioactive decay.

Equations with non-constant coefficients: solution by integrating factor. [2]

Nonlinear first-order equations

Separable equations. Exact equations. Sketching solution trajectories. Equilibrium

solutions, stability by perturbation; examples, including logistic equation and chemical

kinetics. Discrete equations: equilibrium solutions, stability; examples including the

logistic map. [4]

Higher-order linear differential equations

Complementary function and particular integral, linear independence, Wronskian (for

second-order equations), Abel’s theorem. Equations with constant coefficients and

examples including radioactive sequences, comparison in simple cases with difference

equations, reduction of order, resonance, transients, damping. Homogeneous equations.

Response to step and impulse function inputs; introduction to the notions of the Heav-

iside step-function and the Dirac delta-function. Series solutions including statement

only of the need for the logarithmic solution. [8]

Multivariate functions: applications

Directional derivatives and the gradient vector. Statement of Taylor series for functions

. Local extrema of real functions, classification using the Hessian matrix. Coupled

first order systems: equivalence to single higher order equations; solution by matrix

metho ds. Non-degenerate phase portraits local to equilibrium points; stability.

Simple examples of first- and second-order partial differential equations, solution of

the wave equation in the form f (x + ct) + g(x − ct). [5]

Contents

0 Introduction

1 Differentiation

1.1 Differentiation

1.2 Small o and big O notations

1.3 Methods of differentiation

1.4 Taylor’s theorem

1.5 L’Hopital’s rule

2 Integration

2.1 Integration

2.2 Methods of integration

3 Partial differentiation

3.1 Partial differentiation

3.2 Chain rule

3.3 Implicit differentiation

3.4 Differentiation of an integral wrt parameter in the integrand

4 First-order differential equations

4.1 The exponential function

4.2 Homogeneous linear ordinary differential equations

4.3 Forced (inhomogeneous) equations

4.3.1 Constant forcing

4.3.2 Eigenfunction forcing

4.4 Non-constant coefficients

4.5 Non-linear equations

4.5.1 Separable equations

4.5.2 Exact equations

4.6 Solution curves (trajectories)

4.7 Fixed (equilibrium) points and stability

4.7.1 Perturbation analysis

4.7.2 Autonomous systems

4.7.3 Logistic Equation

4.8 Discrete equations (Difference equations)

5 Second-order differential equations

5.1 Constant coefficients

5.1.1 Complementary functions

5.1.2 Second complementary function

5.1.3 Phase space

5.2 Particular integrals

5.2.1 Guessing

5.2.2 Resonance

5.2.3 Variation of parameters

5.3 Linear equidimensional equations

5.4 Difference equations

5.5 Transients and damping

5.6 Impulses and point forces

5.6.1 Dirac delta function

5.7 Heaviside step function

6 Series solutions

7 Directional derivative

7.1 Gradient vector

7.2 Stationary points

7.3 Taylor series for multi-variable functions

7.4 Classification of stationary points

7.5 Contours of f (x, y)

8 Systems of differential equations

8.1 Linear equations

8.2 Nonlinear dynamical systems

9 Partial differential equations (PDEs)

9.1 First-order wave equation

9.2 Second-order wave equation

9.3 The diffusion equation

0 Introduction

In this course, it is assumed that students already know how to do calculus.

While we will define all of calculus from scratch, it is there mostly to introduce

the big and small

notation which will be used extensively in this and future

courses (as well as for the sake of completeness). It is impossible for a person

who hasn’t seen calculus before to learn calculus from those few pages.

Calculus is often used to model physical systems. For example, if we know

that the force

m¨x

on a particle at any time

is given by

−

1, then we

can write this as

m¨x = t

− 1.

We can easily integrate this twice with respect to

, and find the position

as a

function of time.

However, often the rules governing a physical system are not like this. Instead,

the force on the particle is more likely to depend on the position of the particle,

instead of what time it is. Hence the actual equation of motion might be

m¨x = x

− 1.

This is an example of a differential equation. We are given an equation that

a function

obeys, often involving derivatives of

, and we have to find all

functions that satisfy this equation (of course, the first equation is also a

differential equation, but a rather boring one).

A closely related notion is difference equations. These are discrete analogues

of differential equations. A famous example is the Fibonacci sequence, which

states that

n+2

− F

n+1

− F

= 0.

This specifies a relationship between terms in a sequence (

), and we want to

find an explicit solution to this equation.

In this course, we will develop numerous techniques to solve different differ-

ential equations and difference equations. Often, this involves guessing of some

sort.

1 Differentiation

We will first quickly go through basic notions of differentiation and integration.

You should already be familiar with these from A levels or equivalent.

1.1 Differentiation

Definition (Derivative of function). The derivative of a function

(

) with

respect to x, interpreted as the rate of change of f(x) with x, is

= lim

h→0

f(x + h) − f(x)

A function

(

) is differentiable at

if the limit exists (i.e. the left-hand and

right-hand limits are equal).

Example.

(

) =

|x|

is not differentiable at

= 0 as

lim

h→0

|h|−|0|

= 1 and

lim

h→0

−

|h|−|0|

= −1.

Notation. We write

(

) =

(

). Also,



f(x)



(

) =

(x).

Note that the notation

represents the derivative with respect to the

argument. For example, f

(2x) =

d(2x)

1.2 Small o and big O notations

Definition (O and o notations).

(i)

“

(

) =

(

)) as

x → x

” if

lim

x→x

f(x)

g(x)

= 0. Intuitively,

(

) is much

smaller than g(x).

(ii)

“

(

) =

(

)) as

x → x

” if

f(x)

g(x)

is bounded as

x → x

. Intuitively,

f(x) is about as big as g(x).

Note that for f(x) = O(g(x)) to be true, lim

x→x

f(x)

g(x)

need not exist.

Usually,

is either 0 or infinity. Clearly, we have

(

) =

(

)) implies

f(x) = O(g(x)).

Note that this is an abuse of notation. We are not really saying that

(

)

is “equal” to

(

)), since

(

)) itself is not a function. Instead,

(

))

represents a class of functions (namely functions that satisfy the property above),

and we are saying that

is in this class. Technically, a better notation might

(

)

∈ o

(

)), but in practice, writing

(

) =

(

)) is more common and

more convenient.

Example.

– x = o(

√

x) as x → 0 and

√

x = o(x) as x → ∞.

– sin 2x = O(x) as x → 0 as sin θ ≈ θ for small θ.

– sin 2x = O(1) as x → ∞ even though the limit does not exist.

This notation will frequently be used in calculus. For example, if we want to

ignore all terms second order in

in an expression, we can write out the first

order terms and then append +

(

). In particular, we can use it to characterize

derivatives in a different way.

Proposition.

f(x

+ h) = f(x

) + f

)h + o(h)

Proof. We have

) =

f(x

+ h) − f(x

)

o(h)

by the definition of the derivative and the small

notation. The result follows.

1.3 Methods of differentiation

Theorem (Chain rule). Given f(x) = F (g(x)), then

Proof. Assuming that

exists and is therefore finite, we have

= lim

h→0

F (g(x + h)) − F (g(x))

= lim

h→0

F [g(x) + hg

(x) + o(h)] − F (g(x))

= lim

h→0

F (g(x)) + (hg

(x) + o(h))F

(g(x)) + o(hg

(x) + o(h)) − F (g(x))

= lim

h→0

(x)F

(g(x)) +

o(h)

= g

(x)F

(g(x))

Theorem (Product Rule). Give f(x) = u(x)v(x). Then

(x) = u

(x)v(x) + u(x)v

(x).

Theorem (Leibniz’s Rule). Given f = uv, then

(n)

(x) =

r=0





(r)

(n−r)

where f

(n)

is the n-th derivative of f.

1.4 Taylor’s theorem

Theorem (Taylor’s Theorem). For n-times differentiable f, we have

f(x + h) = f(x) + hf

(x) +

(x) + ··· +

(n)

(x) + E

where E

= o(h

) as h → 0. If f

(n+1)

exists, then E

= O(h

n+1

Note that this only gives a local approximation around

. This does not

necessarily tell anything about values of f far from x (but sometimes does).

An alternative form of the sum above is:

f(x) = f(x

) + (x − x

) + ··· +

(x − x

)

(n)

) + E

When the limit as

n → ∞

is taken, the Taylor series of

(

) about the point

x = x

is obtained.

1.5 L’Hopital’s rule

Theorem (L’Hopital’s Rule). Let

(

) and

(

) be differentiable at

, and

lim

x→x

f(x) = lim

x→x

g(x) = 0. Then

lim

x→x

f(x)

g(x)

= lim

x→x

(x)

Proof.

From the Taylor’s Theorem, we have

(

) =

(

) + (

x − x

)

(

) +

o(x − x

), and similarly for g(x). Thus

lim

x→x

f(x)

g(x)

= lim

x→x

f(x

) + (x − x

) + o(x − x

)

g(x

) + (x − x

) + o(x − x

)

= lim

x→x

) +

o(x−x

)

x−x

) +

o(x−x

)

x−x

= lim

x→x

(x)

2 Integration

2.1 Integration

Definition (Integral). An integral is the limit of a sum, e.g.

f(x) dx = lim

∆x→0

n=0

f(x

)∆x.

For example, we can take ∆

b−a

and

∆

. Note that an integral

need not be defined with this particular ∆

and

. The term “integral” simply

refers to any limit of a sum (The usual integrals we use are a special kind known

as Riemann integral, which we will study formally in Analysis I). Pictorially, we

have

a x

···

n+1

···

The area under the graph from

n+1

(

)∆

(∆

). Provided

that f is differentiable, the total area under the graph from a to b is

lim

N→∞

N−1

n=0

(f(x

)∆x)+N·O(∆x

) = lim

N→∞

N−1

n=0

(f(x

)∆x)+O(∆x) =

f(x) dx

Theorem (Fundamental Theorem of Calculus). Let

(

) =

(

) d

. Then

(x) = f(x).

Proof.

F (x) = lim

h→0

x+h

f(t) dt −

f(t) dt

= lim

h→0

x+h

f(t) dt

= lim

h→0

[f(x)h + O(h

)]

= f(x)

Similarly, we have

f(t) dt = −f (x)

and

g(x)

f(t) dt = f (g(x))g

(x).

Notation. We write

(

) d

(

) d

, where the unspecified lower limit

gives rise to the constant of integration.

2.2 Methods of integration

Integration is substantially harder than differentiation. Given an expression,

the product rule and chain rule combined with a few standard derivatives can

usually differentiate it with ease. However, many seemingly-innocent expressions

can turn out to be surprisingly difficult, if not impossible, to integrate. Hence

we have to come up with a lot of tricks to do integrals.

Example (Integration by substitution). Consider

1−2x

√

x−x

. Write

x−x

and du = (1 − 2x) dx. Then the integral becomes

√

= 2

√

u + C = 2

x − x

+ C.

Trigonometric substitution can be performed with reference to the following

table: if the function in the 2nd column is found in the integrand, perform

the substitution in the 3rd column and simplify using the identity in the first

column:

Useful identity Part of integrand Substitution

cos

θ + sin

θ = 1

√

1 − x

x = sin θ

1 + tan

θ = sec

θ 1 + x

x = tan θ

cosh

u − sinh

u = 1

√

− 1 x = cosh u

cosh

u − sinh

u = 1

√

1 + x

x = sinh u

1 − tanh

u = sech

u 1 −x

x = tanh u

Example. Consider

√

2x − x

1 − (x − 1)

. Let

x −

1 =

sin θ

and thus dx = cos θ dθ. The expression becomes

cos

θ dθ =

cos 2θ + 1

dθ

sin 2θ +

θ + C

sin

−1

(x − 1) +

(x − 1)

2x − x

+ C.

Theorem (Integration by parts).

dx = uv −

dx.

Proof.

From the product rule, we have (

)

. Integrating the whole

expression and rearranging gives the formula above.

Example. Consider

∞

−x

. Let

−x

. Then

= 1 and

v = −e

−x

. We have

∞

−x

dx = [−xe

−x

]

∞

−x

= 0 + [−e

−x

]

∞

= 1

Example (Integration by Parts). Consider

log x

. Let

log x

and

= 1.

Then u

and v = x. So we have

log x dx = x log x −

= x log x − x + C

3 Partial differentiation

3.1 Partial differentiation

So far, we have only considered functions of one variable. If we have a function

of multiple variables, say

(

x, y

), we can either differentiate it with respect to

or with respect to y.

Definition (Partial derivative). Given a function of several variables

(

x, y

the partial derivative of

with respect to

is the rate of change of

varies,

keeping y constant. It is given by

∂f

∂x



= lim

δx→0

f(x + δx, y) − f(x, y)

δx

Example. Consider

(

x, y

) =

. Computing the partial derivative

is equivalent to computing the regular derivative with the other variables treated

as constants. e.g.

∂f

∂x



= 2x + y

Second and mixed partial derivatives can also be computed:

∂

∂x

= 2 + y

∂

∂y∂x

∂

∂y



∂f

∂x



= 2ye

+ 2xy

It is often cumbersome to write out the variables that are kept constant. If

all other variables are being held constant, we simply don’t write them out, and

just say

∂f

∂x

Another convenient notation is

Notation.

∂f

∂x

, f

∂

∂y∂x

It is important to know how the order works in the second case. The left

hand side should be interpreted as (

)

. So we first differentiate with respect

to x, and then y. However, in most cases, this is not important, since we have

Theorem. If f has continuous second partial derivatives, then f

= f

We will not prove this statement and just assume it to be true (since this is

an applied course).

3.2 Chain rule

Consider an arbitrary displacement in any direction (

x, y

)

→

(

δx, y

δy

We have

δf = f(x + δx, y + δy) − f(x, y)

= f(x + δx, y + δy) − f(x + δx, y) + f(x + δx, y) − f(x, y)

= f

(x + δx, y)δy + o(δy) + f

(x, y)δx + o(δx)

= (f

(x, y) + o(1))δy + o(δy) + f

(x, y)δx + o(δx)

δf =

∂f

∂x

δx +

∂f

∂y

δy + o(δx, δy)

Take the limit as δx, δy → 0, we have

Theorem (Chain rule for partial derivatives).

df =

∂f

∂x

dx +

∂f

∂y

dy.

Given this form, we can integrate the differentials to obtain the integral form:

df =

∂f

∂x

dx +

∂f

∂y

dy,

or divide by another small quantity. e.g. to find the slope along the path

(x(t), y(t)), we can divide by dt to obtain

∂f

∂x

∂f

∂y

If we pick the parameter t to be the arclength s, we have

∂f

∂x

∂f

∂y







∂f

∂x

∂f

∂y



= ˆs · ∇f,

which is known as the directional derivative (cf. Chapter 7).

Alternatively, the path may also be given by

(

). So

(

x, y

(

)).

Then the slope along the path is

∂f

∂x



∂f

∂y

The chain rule can also be used for the change of independent variables, e.g.

change to polar coordinates x = x(r, θ), y = y(r, θ). Then

∂f

∂θ



∂f

∂x



∂x

∂θ



∂f

∂y



∂y

∂θ



3.3 Implicit differentiation

Consider the contour surface of a function

(

x, y, z

) given by

(

x, y, z

) = const.

This implicitly defines

(

x, y

). e.g. If

(

x, y, z

) =

= 5,

then we can have

5−yz

. Even though

(

x, y

) cannot be found explicitly

(involves solving quintic equation), the derivatives of

(

x, y

) can still be found

by differentiating F (x, y, z) = const w.r.t. x holding y constant. e.g.

∂

∂x

(xy

+ yz

+ z

x) =

∂

∂x

+ 2yz

∂z

∂x

+ z

+ 5z

∂z

∂x

= 0

∂z

∂x

= −

+ z

2yz + 5z

In general, we can derive the following formula:

Theorem (Multi-variable implicit differentiation). Given an equation

F (x, y, z) = c

for some constant c, we have

∂z

∂x



= −

(∂F )/(∂x)

(∂F )/(∂z)

Proof.

dF =

∂F

∂x

dx +

∂F

∂y

dy +

∂F

∂z

∂F

∂x



∂F

∂x



∂F

∂y

∂x



∂F

∂z

∂x



= 0

∂F

∂x

∂F

∂z

∂x



= 0

∂z

∂x



= −

(∂F )/(∂x)

(∂F )/(∂z)

3.4

Differentiation of an integral wrt parameter in the

integrand

Consider a family of functions

(

x, c

). Define

(

b, c

) =

(

x, c

. Then by

the fundamental theorem of calculus, we have

∂I

∂b

(

b, c

). On the other hand,

we have

∂I

∂c

= lim

δc→0

δc

f(x, c + δc)dx −

f(x, c) dx

= lim

δc→0

f(x, c + δc) − f(x, c)

δc

lim

δc→0

f(x, c + δc) − f(x, c)

δc

∂f

∂c

In general, if I(b(x), c(x)) =

b(x)

f(y, c(x))dy, then by the chain rule, we have

∂I

∂b

∂I

∂c

= f(b, c)b

(x) + c

(x)

∂f

∂c

dy.

So we obtain

Theorem (Differentiation under the integral sign).

b(x)

f(y, c(x)) dy = f(b, c)b

(x) + c

(x)

∂f

∂c

This is sometimes a useful technique that allows us to perform certain

otherwise difficult integrals, as you will see in the example sheets. However, here

we will only have time for a boring example.

Example. Let I =

−λx

dx. Then

dλ

−x

−λx

dx.

If I =

−λx

dx. Then

dλ

= e

−λ

−x

−λx

dx.

4 First-order differential equations

A differential equation is an equation that involves derivatives, such as

= 0. Unlike regular equations where the solution is a number, the solution to

a differential equation is a function that satisfies the equation. In this chapter,

we will look at first-order differential equations, where only first derivatives are

involved.

4.1 The exponential function

Often, the solutions to differential equations involve the exponential function.

Consider a function f(x) = a

, where a > 0 is constant.

The derivative of this function is given by

= lim

h→0

x+h

− a

= a

lim

h→0

− 1

= λa

= λf(x)

where

lim

h→0

− 1

(0) = const. So the derivative of an exponential

function is a multiple of the original function. In particular,

Definition (Exponential function).

exp

(

) =

is the unique function

satisfying f

(x) = f(x) and f(0) = 1.

We write the inverse function as ln x or log x.

Then if y = a

= e

x ln a

, then y

= e

x ln a

ln a = a

ln a. So λ = ln a.

Using this property, we find that the value of e is given by

e = lim

k→∞



1 +



≈ 2.718281828 ··· .

The importance of the exponential function lies in it being an eigenfunction of

the differential operator.

Definition (Eigenfunction). An eigenfunction under the differential operator

is a function whose functional form is unchanged by the operator. Only its

magnitude is changed. i.e.

= λf

Example. e

is an eigenfunction since

= me

4.2 Homogeneous linear ordinary differential equations

Before we start, we will first define the terms used in the title. These are useful

criteria for categorizing differential equations.

Definition (Linear differential equation). A differential equation is linear if the

dependent variable (y, y

, y

etc.) appears only linearly.

Definition (Homogeneous differential equation). A differential equation is

homogeneous if y = 0 is a solution.

Definition (Differential equation with constant coefficients). A differential

equation has constant coefficients if the independent variable

does not appear

explicitly.

Definition (First-order differential equation). A differential equation is first-

order if only first derivatives are involved.

Theorem. Any linear, homogeneous, ordinary differential equation with constant

coefficients has solutions of the form e

This theorem is evident when we consider an example.

Example. Given 5

− 3y = 0. Try y = e

. Then

= me

5me

− 3e

= 0

Since this must hold for all values of

, there exists some value

for which

= 0 and we can divide by

(note in this case this justification is not

necessary, because

is never equal to 0. However, we should justify as above

if we are dividing, say,

). Thus 5

m −

3 = 0 and

= 3

5. So

3x/5

is a

solution.

Because the equation is linear and homogeneous, any multiple of a solution

is also a solution. Therefore y = Ae

3x/5

is a solution for any value of A.

But is this the most general form of the solution? One can show that an

-order linear differential equation has

and only

independent solutions. So

y = Ae

3x/5

is indeed the most general solution.

We can determine A by applying a given boundary condition.

Discrete equations

Suppose that we are given the equation 5

−

= 0 with boundary condition

= 0. This gives a unique function

. We can approximate this by

considering discrete steps of length

between

and

n+1

. (Using the simple

Euler numerical scheme,) we can approximate the equation by

n+1

− y

− 3y

≈ 0.

Rearranging the terms, we have

n+1

≈

(1 +

)

. Applying the relation

successively, we have



1 +



n−1



1 +



1 +



n−2



1 +



For each given value of x, choose h = x/n, so

= y



1 +

(x/n)



Taking the limit as n → ∞, we have

y(x) = lim

n→∞



1 +

3x/5



= y

3x/5

in agreement with the solution of the differential equation.

solution of diff. eq.

discrete approximation

Series solution

We can also try to find a solution in the form of a Taylor Series

∞

n=0

We have y

n−1

. Substituting these into our equation

− 3y = 0

5(xy

) − 3x(y) = 0

(5n − 3x)x

= 0

Consider the coefficient of

: 5

−

n−1

= 0. Since this holds for all values

, when

= 0, we get 0

= 0. This tells that

is arbitrary. If

n >

0, then

n−1

n(n − 1)

n−2

= ··· =





Therefore we have

y = a

∞

n=0





= a

3x/5

4.3 Forced (inhomogeneous) equations

Recall that a homogeneous equation is an equation like 5

−

= 0, with no

or constant terms floating around. A forced, or inhomogeneous, equation is

one that is not homogeneous. For example, 5

−

= 10 or 5

−

are

forced equations. We call the “extra” terms 10 and x

the forcing terms.

4.3.1 Constant forcing

Example. Consider 5

−

= 10. We can spot that there is a equilibrium

(constant) solution y = y

= −

with y

= 0.

The particular solution

is a solution of the ODE. Now suppose the general

solution is

. We see that 5

−

= 0. So

satisfies the homogeneous

equation we already solved, and

y = −

+ Ae

3x/5

Note that any boundary conditions to determine

must be applied to the full

solution y and not the complementary function y

This is the general method of solving forced equations. We first find one

particular solution to the problem, often via educated guesses. Then we solve

the homogeneous equation to get a general complementary solution. Then the

general solution to the full equation is the sum of the particular solution and

the general complementary solution.

4.3.2 Eigenfunction forcing

This is the case when the forcing term is an eigenfunction of the differential

operator.

Example. In a radioactive rock, isotope A decays into isotope B at a rate

proportional to the number

of remaining nuclei A, and B also decays at a rate

proportional to the number b of remaining nuclei B. Determine b(t).

We have

= −k

= k

a − k

Solving the first equation, we obtain a = a

−k

. Then we have

+ k

b = k

−k

We usually put the variables involving

on the left hand side and the others on

the right. The right-hand term k

−k

is the forcing term.

Note that the forcing term is an eigenfunction of the differential operator

on the LHS. So that suggests that we can try a particular integral

−k

Substituting it in, we obtain

−k

C + k

C = k

C =

− k

Then write

. We get

= 0 and

−k

. All together, we

have the general solution

b =

− k

−k

+ De

−k

Assume the following boundary condition:

= 0 when

= 0, in which case we

can find

b =

− k



−k

− e

−k



The isotope ratio is

− k

1 − e

−k

So given the current ratio

b/a

, with laboratory determined rates

and

, we

can determine the value of t, i.e. the age of the rock.

4.4 Non-constant coefficients

Consider the general form of equation

a(x)y

+ b(x)y = c(x).

Divide by a(x) to get the standard form

+ p(x)y = f(x).

We solve this by multiplying an integrating factor

(

) to obtain (

)

+ (

µp

)

µf.

We want to choose a

such that the left hand side is equal to (

µy

)

. By the

product rule, we want µp = µ

, i.e.

p =

dµ

p dx =

dµ

= ln µ(+C)

µ = exp



p dx



Then by construction, we have (µy)

= µf and thus

y =

µf dx

, where µ = exp



p dx



Example. Consider

+ (1

−x

)

= 1. To obtain it in standard form, we have

1−x

y =

. We have µ = exp



(

− 1) dx



= e

ln x−x

= xe

−x

. Then

y =

−x

−e

−x

+ C

−x

−1

Suppose that we have a boundary condition

is finite at

= 0. Since we have

−1

, we have to ensure that

−

→

0 as

x →

0. Thus

= 1, and by

L’Hopital’s rule, y → 1 as x → 0.

4.5 Non-linear equations

In general, a first-order equation has the form

Q(x, y)

+ P (x, y) = 0.

(this is not exactly the most general form, since theoretically we can have powers

or even other complicated functions of

, but we will not consider that

case here).

4.5.1 Separable equations

Definition (Separable equation). A first-order differential equation is separable

if it can be manipulated into the following form:

q(y) dy = p(x) dx.

in which case the solution can be found by integration

q(y) dy =

p(x) dx.

This is the easy kind.

Example.

y − 3y)

− 2xy

= 4x

4x + 2xy

y − 3y

2x(2 + y

)

y(x

− 3)

2 + y

dy =

− 3

2 + y

dy =

− 3

ln(2 + y

) = ln(x

− 3) + C

2 + y

= ln A(x

− 3)

+ 2 = A(x

− 3)

4.5.2 Exact equations

Definition (Exact equation).

(

x, y

)

(

x, y

) = 0 is an exact equation iff

the differential form

(

x, y

) d

(

x, y

) d

is exact, i.e. there exists a function

f(x, y) for which

df = Q(x, y) dy + P (x, y) dx

(

x, y

) d

(

x, y

) d

is an exact differential of

, then d

(

x, y

) d

(

x, y

) d

. But by the chain rule, d

∂f

∂x

∂f

∂y

and this equality holds

for any displacements dx, dy. So

∂f

∂x

= P,

∂f

∂y

= Q.

From this we have

∂

∂y∂x

∂P

∂y

∂

∂x∂y

∂Q

∂x

We know that the two mixed 2nd derivatives are equal. So

∂P

∂y

∂Q

∂x

The converse is not necessarily true. Even if this equation holds, the differential

need not be exact. However, it is true if the domain is simply-connected.

Definition (Simply-connected domain). A domain

is simply-connected if it

is connected and any closed curve in

can be shrunk to a point in

without

leaving D.

Example. A disc in 2D is simply-connected. A disc with a “hole” in the middle

is not simply-connected because a loop around the hole cannot be shrunk into a

point. Similarly, a sphere in 3D is simply-connected but a torus is not.

Theorem. If

∂P

∂y

∂Q

∂x

through a simply-connected domain

, then

is an exact differential of a single-valued function in D.

If the equation is exact, then the solution is simply

= constant, and we

can find f by integrating

∂f

∂x

= P and

∂f

∂y

= Q.

Example.

6y(y − x)

+ (2x − 3y

) = 0.

We have

P = 2x − 3y

, Q = 6y(y − x).

Then

∂P

∂y

∂Q

∂x

= −6y. So the differential form is exact. We now have

∂f

∂x

= 2x − 3y

∂f

∂y

= 6y

− 6xy.

Integrating the first equation, we have

f = x

− 3xy

+ h(y).

Note that since it was a partial derivative w.r.t.

holding

constant, the

“constant” term can be any function of

. Differentiating the derived

w.r.t

we have

∂f

∂y

= −6xy + h

(y).

Thus h

(y) = 6y

and h(y) = 2y

+ C, and

f = x

− 3xy

+ 2y

+ C.

Since the original equation was d

= 0, we have

= constant. Thus the final

solution is

− 3xy

+ 2y

= C.

4.6 Solution curves (trajectories)

Example. Consider the first-order equation

= t(1 − y

We can solve it to obtain

1 − y

= t dt

1 + y

1 − y

+ C

1 + y

1 − y

= Ae

y =

A − e

−t

A + e

−t

We can plot the solution for different values of

and obtain the following graph:

−1

But can we understand the nature of the family of solutions without solving

the equation? Can we sketch the graph without solving it?

We can spot that

1 are two constant solutions and we can plot them

first. We also note (and mark on our graph) that y

= 0 at t = 0 for any y.

Then notice that for t > 0, y

> 0 if −1 < y < 1. Otherwise, y

< 0.

Now we can find isoclines, which are curves along which

(i.e.

) is constant:

− y

) =

for some constant

. Then

= 1

− D/t

. After marking a few

isoclines, can sketch the approximate form of our solution:

−1

In general, we sketch the graph of a differential equation

= f(t, y)

by locating constant solutions (and determining stability: see below) and iso-

clines.

4.7 Fixed (equilibrium) points and stability

Definition (Equilibrium/fixed point). An equilibrium point or a fixed point of

a differential equation is a constant solution

. This corresponds to

= 0

for all t.

These are usually the easiest solutions to find, since we are usually given an

expression for

and we can simply equate it to zero.

Definition (Stability of fixed point). An equilibrium is stable if whenever

is deviated slightly from the constant solution

y → c

t → ∞

. An

equilibrium is unstable if the deviation grows as t → ∞.

The objective of this section is to study how to find equilibrium points and

their stability.

Example. Referring to the differential equation above (

˙y

− y

)), we see

that the solutions converge towards

= 1, and this is a stable fixed point. They

diverge from y = −1, and this is an unstable fixed point.

4.7.1 Perturbation analysis

Perturbation analysis is used to determine stability. Suppose

is a fixed

point of

(

y, t

), so

(

a, t

) = 0. Write

(

), where

(

) is a small

perturbation from

. We will later assume that

is arbitrarily small. Putting

this into the differential equation, we have

dε

= f(a + ε, t)

= f(a, t) + ε

∂f

∂y

(a, t) + O(ε

)

= ε

∂f

∂y

(a, t) + O(ε

)

Note that this is a Taylor expansion valid when

ε 

1. Thus

(

) can be

neglected and

dε

∼

∂f

∂y

We can use this to study how the perturbation grows with time.

This approximation is called a linearization of the differential equation.

Example. Using the example ˙y = t(1 − y

) above, we have

∂f

∂y

= −2yt =

(

−2t at y = 1

2t at y = −1

= 1,

˙ε

−

tε

and

−t

. Since

ε →

0 as

t → ∞

y →

1. So

= 1 is a

stable fixed point.

On the other hand, if we consider

−

1, then

˙ε

= 2

tε

and

. Since

ε → ∞ as t → ∞, y = −1 is unstable.

Technically

ε → ∞

is not a correct statement, since the approximation used

is only valid for small

. But we can be sure that the perturbation grows (even

if not → ∞) as t increases.

4.7.2 Autonomous systems

Often, the mechanics of a system does not change with time. So

˙y

is only a

function of y itself. We call this an autonomous system.

Definition (Autonomous system). An autonomous system is a system in the

form ˙y = f(y), where the derivative is only (explicitly) dependent on y.

Near a fixed point

, where

(

) = 0, write

(

). Then

˙ε

(

) =

kε

for some constant

. Then

. The stability of the

system then depends solely on sign of k.

Example. Consider a chemical reaction NaOH + HCl

→

O + NaCl. We

have

NaOH + HCl → H

O + NaCl

Number of molecules a b c c

Initial number of molecules a

0 0

If the reaction is in dilute solution, then the reaction rate is proportional to

Thus

= λab

= λ(a

− c)(b

− c)

= f(c)

We can plot

as a function of c, and wlog a

< b

˙c

We can also plot a phase portrait, which is a plot of the dependent variable only,

where arrows show the evolution with time,

We can see that the fixed point c = a

is stable while c = b

is unstable.

We can also solve the equation explicitly to obtain

c =

[1 − e

−(b

−a

)λt

]

− a

−λ(b

−a

Of course, this example makes no sense physically whatsoever because it assumes

that we can have negative values of

and

. Clearly in reality we cannot have,

say,

−

1 mol of NaOH, and only solutions for

c ≤ a

are physically attainable,

and in this case, any solution will tend towards c = a

4.7.3 Logistic Equation

The logistic equation is a simple model of population dynamics. Suppose we

have a population of size

. It has a birth rate

αy

and a death rate

βy

. With

this model, we obtain

= (α − β)y

y = y

(α−β)t

Our population increases or decreases exponentially depending on whether the

birth rate exceeds death rate or vice versa.

However, in reality, there is fighting for limited resources. The probability of

some piece of food (resource) being found is

∝ y

. The probability of the same

piece of food being found by two individuals is

∝ y

. If food is scarce, they fight

(to the death), so death rate due to fighting (competing) is γy

for some γ. So

= (α − β)y − γy

= ry



1 −



where

α − β

and

r/γ

. This is the differential logistic equation. Note

that it is separable and can be solved explicitly.

However, we find the phase portrait instead. If

α − β >

0, then the

graph is a quadratic parabola, and we see that Y is a stable fixed point.

O Y

Now when the population is small, we have

˙y ' ry

So the population grows exponentially. Eventually, the stable equilibrium

reached.

4.8 Discrete equations (Difference equations)

Since differential equations are approximated numerically by computers with

discrete equations, it is important to study the behaviour of discrete equations,

and their difference with continuous counterparts.

In the logistic equation, the evolution of species may occur discretely (e.g.

births in spring, deaths in winter), and we can consider the population at certain

time intervals (e.g. consider the population at the end of each month). We might

have a model in the form

n+1

= λx

(1 − x

This can be derived from the continuous equation with a discrete approximation

n+1

− y

∆t

= ry



1 −



n+1

= y

+ r∆ty



1 −



= (1 + r∆t)y

−

r∆t

= (1 + r∆t)y



1 −



r∆t

1 + r∆t





Write

λ = 1 + r∆t, x



r∆t

1 + r∆t



then

n+1

= λx

(1 − x

This is the discrete logistic equation or logistic map.

If λ < 1, then deaths exceed births and the population decays to zero.

n+1

= x

We see that x = 0 is a fixed point.

In general, to find fixed points, we solve for

n+1

, i.e.

(

) =

. For

the logistic map, we have

λx

(1 − x

) = x

[1 − λ(1 − x

)] = 0

= 0 or x

= 1 −

When 1 < λ < 2, we have

n+1

= x

We see that

= 0 is an unstable fixed point and

= 1

−

is a stable fixed

point.

When 2 < λ < 3, we have

n+1

= x

There is an oscillatory convergence to x

= 1 −

When

λ >

3, we have a limit cycle, in which

oscillates between 2 values,

i.e. x

n+2

= x

. When λ = 1 +

√

6 ≈ 3.449, we have a 4-cycle, and so on.

n+1

= x

We have the following plot of the stable solutions for different values of

(plotted

as r in the horizontal axis)

Credits: Wikimedia Commons: Jordan Pierce - Public Domain CC0 1.0

Note that the fixed point still exists after

= 3, but is no longer stable. Similarly,

the 2-cycles still exist after λ = 1 +

√

6, but it is not stable.

5 Second-order differential equations

We now move on to second-order differential equations. While we will only look

at second-order equations, most of the methods in this section apply to higher

order differential equations as well.

5.1 Constant coefficients

The general form of an equation with constant coefficients is

+ by

+ cy = f(x).

We solve this in two steps:

(i)

Find the complementary functions which satisfy the homogeneous equation

+ by

+ cy = 0.

(ii) Find a particular solution that satisfies the full equation.

5.1.1 Complementary functions

Recall that

λx

is an eigenfunction of the differential operator

. Hence it is

also an eigenfunction of the second derivative





If the complementary function has the form

λx

, then

λe

λx

and

= λ

λx

. Substituting into the differential equation gives

Definition (Characteristic equation). The characteristic equation of a (second-

order) differential equation ay

+ by

+ c = 0 is

aλ

+ bλ + c = 0.

In this case there are two solutions to the characteristic equation, giving (in

principle) two complementary functions y

= e

and y

= e

and

are distinct, then

and

are linearly independent and complete

— they form a basis of the solution space. The (most) general complementary

function is

= Ae

+ Be

Example.

−

+ 6

= 0. Try

λx

. The characteristic equation is

− 5λ + 6 = 0. Then λ = 2 or 3. So the general solution is y = Ae

+ Be

Note that A and B can be complex constants.

Example (Simple harmonic motion).

+ 4

= 0. Try

λx

. The character-

istic equation is

+ 4 = 0, with solutions

. Then our general solution

2ix

−2ix

. However, if this is in a case of simple harmonic motion

in physics, we want the function to be real (or look real). We can write

y = A(cos 2x + i sin 2x) + B(cos 2x −i sin 2x)

= (A + B) cos 2x + i(A − B) sin 2x

= α cos 2x + β sin 2x

where α = A + B and β = i(A − B), and α and β are independent constants.

In effect, we have changed the basis from {e

2ix

, e

−2ix

} to {cos 2x, sin 2x}.

Example (Degeneracy). y

− 4y

+ 4y = 0.

Try

λx

. We have

−

+ 4 = 0 and (

λ −

= 0. So

= 2 or 2. But

and

are clearly not linearly independent. We have only managed to find

one basis function of the solution space, but a second order equation has a 2

dimensional solution space. We need to find a second solution.

We can perform detuning. We can separate the two functions found above

from each other by considering

−

+ (4

− ε

)

= 0. This turns into the

equation we want to solve as

ε →

0. Try

λx

. We obtain

−

+ 4

− ε

The two roots are λ = 2 ± ε. Then

y = Ae

(2+ε)x

+ Be

(2−ε)x

= e

[Ae

εx

+ Be

−εx

]

Taking the limit as ε → 0, we use the Taylor expansion of e

εx

to obtain

y = e

[(A + B) + εx(A − B) + O(Aε

, Bε

)]

We let (

) =

and

(

A − B

) =

. This is perfectly valid for any non-zero

ε. Then A =

(α +

) and B =

(α −

). So we have

y = e

[α + βx + O(Aε

, Bε

)]

We know for any

, we have a solution of this form. Now we turn the procedure

around. We fix some

and

. Then given any

, we can find some constants

(depending on

) such that the above holds. As we decrease the size of

, we

have A, B = O(

). So O(Aε

) = O(Bε

) = O(ε). So our solution becomes

y = e

[α + βx + O(ε)]

→ e

(α + βx)

In this way, we have derived two separate basis functions. In general, if

(

)

is a degenerate complementary function of a linear differential equation with

constant coefficients, then

(

) =

(

) is an independent complementary

function.

5.1.2 Second complementary function

In general (i.e. if we don’t have constant coefficients), we can find a second com-

plementary function associated with a degenerate solution of the homogeneous

equation by looking for a solution in the form

(

) =

(

)

(

), where

(

) is

the degenerate solution we found.

Example. Consider

−

+ 4

= 0. We have

. We try

Then

= (v

+ 2v)e

= (v

+ 4v

+ 4v)e

Substituting into the original equation gives

+ 4v

+ 4v) − 4(v

+ 2v) + 4v = 0.

Simplifying, this tells us v

= 0, which forces v to be a linear function of x. So

= (Ax + B)e

for some A, B ∈ R.

5.1.3 Phase space

If we are given a general nth order differential equation of the form

(x)y

(n)

+ a

n−1

(n−1)

+ ··· + a

(x)y

+ a

(x)y = f(x),

and we have a solution

, then we can plot a graph of

versus

, and see how

evolves with x.

However, one problem is that for such an equation, the solution is not just

determined by the initial condition

(

), but also

(

) etc. So if we

just have a snapshot of the value of

at a particular point

, we have completely

no idea how it would evolve in the future.

So how much information do we actually need? At any point

, if we

are given the first

n −

1 derivatives, i.e.

(

···

(n−1)

(

), we can

then get the

th derivative and also any higher derivatives from the differential

equation. This means that we know the Taylor series of

about

, and it follows

that the solution is uniquely determined by these conditions (note that it takes

considerably more extra work to actually prove rigorously that the solution is

uniquely determined by these initial conditions, but it can be done for sufficiently

sensible f , as will be done in IB Analysis II).

Thus we are led to consider the solution vector

Y(x) = (y(x), y

(x), ··· , y

n−1

(x)).

We say such a vector lies in the phase space, which is an

-dimensional space. So

for each

, we thus get a point Y(

) lying in the

-dimensional space. Moreover,

given any point in the phase space, if we view it as the initial conditions for our

differential equation, we get a unique trajectory in the phase space.

Example. Consider

+ 4

= 0. Suppose we have an initial condition of

(0) = 1

, y

(0) = 0. Then we can solve the equation to obtain

(

) =

cos

Thus the initial solution vector is Y

(0) = (1

0), and the trajectory as

varies is

given by Y

(

) = (

cos

x, −

sin

). Thus as

changes, we trace out an ellipse

in the clockwise direction:

(x)

Another possible initial condition is

(0) = 0

, y

(0) = 2. In this case, we obtain

the solution y(x) = sin 2x, with a solution vector Y

(x) = (sin 2x, 2 cos 2x).

Note that as vectors, our two initial conditions (1

0) and (0

2) are indepen-

dent. Moreover, as

changes, the two solution vectors Y

(

)

(

) remain

independent. This is an important observation that allows the method of varia-

tion of parameters later on.

In general, for a 2nd order equation, the phase space is a 2-dimensional space,

and we can take the two complementary functions Y

and Y

as basis vectors

for the phase space at each particular value of

. Of course, we need the two

solutions to be linearly independent.

Definition (Wronskian). Given a differential equation with solutions

, y

, the

Wronskian is the determinant

W (x) =



Definition (Independent solutions). Two solutions

(

) and

(

) are indepen-

dent solutions of the differential equation if and only if Y

and Y

are linearly

independent as vectors in the phase space for some

, i.e. iff the Wronskian is

non-zero for some x.

In our example, we have W (x) = 2 cos

2x + 2 sin

2x = 2 6= 0 for all x.

Example. In our earlier example, y

= e

and y

= xe

. We have

W =



+ 2xe



= e

(1 + 2x − 2x) = e

6= 0.

In both cases, the Wronskian is never zero. Is it possible that it is zero for

some x while non-zero for others? The answer is no.

Theorem (Abel’s Theorem). Given an equation

(

)

(

)

= 0, either

= 0 for all

, or

W 6

= 0 for all

. i.e. iff two solutions are independent for

some particular x, then they are independent for all x.

Proof. If y

and y

are both solutions, then

+ py

+ qy

) = 0

+ py

+ qy

) = 0

Subtracting the two equations, we have

− y

+ p(y

− y

) = 0

Note that

−y

and

−

(

) =

−y

+ P (x)W = 0

W (x) = W

−

P dx

Where

= const. Since the exponential function is never zero, either

= 0,

in which case W = 0, or W

6= 0 and W 6= 0 for any value of x.

In general, any linear

th-order homogeneous differential equation can be

written in the form Y

Y = 0, a system of first-order equations. It can then

be shown that

(

)

= 0, and

−

tr A dx

. So Abel’s theorem

holds.

5.2 Particular integrals

We now consider equations of the form

(

). We will come up

with several ways to find a particular integral.

5.2.1 Guessing

If the forcing terms are simple, we can easily “guess” the form of the particular

integral, as we’ve previously done for first-order equations.

f(x) y

(x)

sin kx

A sin kx + B cos kx

cos kx

polynomial p

(x) q

(x) = a

+ ··· + a

x + a

It is important to remember that the equation is linear, so we can superpose

solutions and consider each forcing term separately.

Example. Consider

−

+ 6

= 2

. To obtain the forcing term 2

we need a first order polynomial

, and to get

we need

. Thus we

can guess

= ax + b + ce

= a + 4ce

= 16ce

Substituting in, we get

16ce

− 5(a + 4ce

) + 6(ax + b + ce

) = 2x + e

Comparing coefficients of similar functions, we have

16c − 20c + 6c = 1 ⇒ c =

6a = 2 ⇒ a =

−5a + 6b = 0 ⇒ b =

Since the complementary function is

, the general solution is

y = Ae

+ Be

x +

Note that any boundary condition to determine

and

must be applied to

the full solution, not the complementary function

5.2.2 Resonance

Consider

¨y

sin ω

. The complementary solution is

A sin ω

B cos w

. We notice that the forcing is itself a complementary function. So

if we guess a particular integral

C sin ω

D cos ω

, we’ll simply find

¨y

+ ω

= 0, so we can’t balance the forcing.

This is an example of a simple harmonic oscillator being forced at its natural

frequency.

We can detune our forcing away from the natural frequency, and consider

¨y + ω

y = sin ωt with ω 6= ω

. Try

= C(sin ωt − sin ω

t).

We have

¨y

= C(−ω

sin ωt + ω

sin ω

t).

Substituting into the differential equation, we have C(ω

− ω

) = 1. Then

sin ωt − sin ω

− ω

We can simplify this to

− ω

cos



+ ω



sin



ω − ω



We let ω

− ω = ∆ω. Then

−2

(2ω + ∆ω)∆ω

cos



ω +

∆ω





sin



∆ω



sin



∆ω



−sin



∆ω





∆ω



This oscillation in the amplitude of the

cos

wave is known as beating. This

happens when the forcing frequency is close to the natural frequency. The

wavelength of the

sin

function has order

(

∆ω

) and

cos

has wavelength

(

As ∆

ω →

0, the wavelength of the beating envelope

→ ∞

and we just have the

initial linear growth.

Mathematically, since sin θ ≈ θ as θ → 0, as ∆ω → 0, we have

→

−t

2ω

cos ω

In general, if the forcing is a linear combination of complementary functions,

then the particular integral is proportional to

(the independent variable) times

the non-resonant guess.

5.2.3 Variation of parameters

So far, we have been finding particular integrals by guessing. Here we are going

to come up with a method that can systematically help us find the particular

integral. Of course, this is substantially more complicated than guessing. So if

the form of the particular integral is obvious, we should go for guessing.

Suppose we have a second order differential equation

+ p(x)y

+ q(x)y = f(x).

We then know any particular solution vector can be written in the form

Y(x) =



y(x)

(x)



and our job is to find one solution vector that satisfies the equation. We

presuppose such a solution actually exists, and we will try to find out what it is.

The trick is to pick a convenient basis for this space. Let

(

) and

(

) be

linearly independent complementary functions of the ODE. Then for each

, the

solution vectors Y

(

) = (

(

)

, y

(

)) and Y

(

) = (

(

)

, y

(

)) form a basis

of the solution space. So for each particular

, we can find some constants

u, v

(depending on x) such that the following equality holds:

Y(x) = u(x)Y

(x) + v(x)Y

(x),

since Y

and Y

are a basis.

Component-wise, we have

= uy

+ vy

(a)

= uy

+ vy

(b)

Differentiating the second equation, we obtain

= (uy

+ u

) + (vy

+ v

) (c)

If we consider (c) + p(b) + q(a), we have y

+ y

= f.

Now note that we derived the equation of

from the vector equation. This

must be equal to what we get if we differentiate (a). By (a)

- (b), we obtain

= 0. Now we have two simultaneous equations for

and

, which

we should be able to solve.

We can, for example, write them in matrix form as











Inverting the left matrix, we have







−y





So u

= −

f and v

Example.

+ 4

sin

. We know that

sin

and

cos

W = −2. We write

= u sin 2x + v cos 2x.

Using the formulae above, we obtain

cos 2x sin 2x

sin 4x

, v

−sin

cos 4x − 1

u = −

cos 4x

, v =

sin 4x

−

Therefore

(−cos 4x sin 2x + sin 4x cos 2x −

cos 2x) =

sin 2x −

cos 2x

Note that

−

x cos

is what we found previously by detuning, and

sin

a complementary function, so the results agree.

It is generally not a good idea to remember the exact formula for the results

we’ve obtained above. Instead, whenever faced with such questions, you should

be able to re-derive the results instead.

5.3 Linear equidimensional equations

Equidimensional equations are often called homogeneous equations, but this is

confusing as it has the same name as those with no forcing term. So we prefer

this name instead.

Definition (Equidimensional equation). An equation is equidimensional if it

has the form

+ bxy

+ cy = f(x),

where a, b, c are constants.

To understand the name “equidimensional”, suppose we are doing physics

and variables have dimensions. Say

has dimensions

and

has dimensions

Then

has dimensions

−1

and

has dimensions

−2

. So all terms

and y have the same dimensions.

Solving by eigenfunctions

Note that

is an eigenfunction of

. We can try an eigenfunction

We have

k−1

and thus

; and

(

k −

k−2

and

= k(k − 1)x

Substituting in, we have

ak(k − 1) + bk + c = 0,

which we can solve, in general, to give two roots

and

, and

Solving by substitution

Alternatively, we can make a substitution z = ln x. Then we can show that

+ (b − a)

+ cy = f(e

This turns an equidimensional equation into an equation with constant coeffi-

cients. The characteristic equation for solutions in the form

λz

is of form

+ (

b −a

)

= 0, which we can rearrange to become

aλ

(

λ −

1) +

bλ

= 0.

So λ = k

, k

Then the complementary function is y

= Ae

+ Be

= Ax

+ Bx

Degenerate solutions

If the roots of the characteristic equation are equal, then

, ze

}

, x

ln x}

. Similarly, if there is a resonant forcing proportional to

(or

then there is a particular integral of the form x

ln x.

These results can be easily obtained by considering the substitution method

of solving, and then applying our results from homogeneous linear equations

with constant coefficients.

5.4 Difference equations

Consider an equation of the form

n+2

+ by

n+1

+ cy

= f

We can solve in a similar way to differential equations, by exploiting linearity

and eigenfunctions.

We can think of the difference operator

[

] =

n+1

. This has an eigen-

function y

= k

. We have D[y

] = D[k

] = k

n+1

= k · k

= ky

To solve the difference equation, we first look for complementary functions

satisfying

n+2

+ by

n+1

+ cy

= 0

We try y

= k

to obtain

n+2

+ bk

n+1

+ ck

= 0

+ bk + c = 0

from which we can determine

. So the general complementary function is

= Ak

+ Bk

if k

6= k

. If they are equal, then y

= (A + Bn)k

To find the particular integral, we guess.

if k 6= k

, k

+ Bn

p−1

+ ··· + Cn + D

Example (Fibonacci sequence). The Fibonacci sequence is defined by

= y

n−1

+ y

n−2

with y

= y

= 1.

We can write this as

n+2

− y

n+1

− y

= 0

We try y

= k

. Then k

− k − 1 = 0. Then

− k − 1 = 0

k =

1 ±

√

We write k = ϕ

, ϕ

. Then y

= Aϕ

+ Bϕ

. Our initial conditions give

A + B = 1

Aϕ

+ Bϕ

= 1

We get A =

√

and B =

−ϕ

√

. So

n+1

− ϕ

n+1

√

n+1

−



−1



n+1

√

5.5 Transients and damping

In many physical systems, there is some sort of restoring force and some damping,

e.g. car suspension system.

Consider a car of mass

with a vertical force

(

) acting on it (e.g. mouse

jumping on the car). We can consider the wheels to be springs (

) with a

“shock absorber” (F = l ˙x). Then the equation of motion can be given by

M ¨x = F (t) − kx − l ˙x.

So we have

¨x +

˙x +

x =

F (t).

Note that if we don’t have the damping and the forcing, we end up with simple

harmonic motion of angular frequency

k/M

. Write

M/k

, where

is dimensionless. The timescale

M/k

is proportional to the period of the

undamped, unforced system (or 1 over its natural frequency). Then we obtain

¨x + 2κ ˙x + x = f(τ)

where, ˙x means

dτ

, κ =

√

and f =

By this substitution, we are now left with only one parameter

instead of

the original three (M, l, k).

We will consider different possible cases.

Free (natural) response f = 0

¨x + 2κ ˙x + x = 0

We try x = e

λτ

+ 2κλ + 1 = 0

λ = −κ ±

− 1

= −λ

, −λ

where λ

and λ

have positive real parts.

Underdamping

If κ < 1, we have x = e

−κτ

(A sin

√

1 − κ

τ + B cos

√

1 − κ

τ).

The period is

2π

√

1−κ

and its amplitude decays in a characteristic of

(

Note that the damping increases the period. As

κ →

1, the oscillation period

→ ∞.

Critically damping

If κ = 1, then x = (A + Bτ)e

−κτ

The rise time and decay time are both

(

) =

(1). So the dimensional rise

and decay times are O(

M/k).

Overdamping

κ >

1, then

−λ

with

< λ

. Then the decay time is

O(1/λ

) and the rise time is O(1/λ

Note that in all cases, it is possible to get a large initial increase in amplitude.

Forcing

In a forced system, the complementary functions typically determine the short-

time transient response, while the particular integral determines the long-time

(asymptotic) response. For example, if

(

) =

sin τ

, then we can guess

C sin τ + D cos τ . In this case, it turns out that x

= −

2κ

cos τ.

The general solution is thus

−λ

−λτ

−

2κ

cos τ ∼ −

2κ

cos τ

τ → ∞ since Re(λ

1,2

) > 0.

It is important to note that the forcing response is out of phase with the

forcing.

5.6 Impulses and point forces

5.6.1 Dirac delta function

Consider a ball bouncing on the ground. When the ball hits the ground at some

time

, it experiences a force from the ground for some short period of time.

The force on the ball exerted by the ground

(

) is 0 for most of the time, except

during the short period (T −ε, T + ε).

Often we don’t know (or we don’t wish to know) the details of

(

) but we

can note that it only acts for a short time of

(

) that is much shorter than the

overall time

(

−t

) of the system. It is convenient mathematically to imagine

the force acting instantaneously at time t = T , i.e. consider the limit ε → 0.

Newton’s second law gives

m¨x

(

)

− mg

. While we cannot solve it, we

can integrate the equation from T − ε to T + ε. So

T +ε

T −ε

dt =

T +ε

T −ε

F (t) dt −

T +ε

T −ε

mg dt





T +ε

T −ε

= I − 2εmg

∆p = I − O(ε)

Where ∆

is the change in momentum and the impulse

T +ε

T −ε

(

) d

is the

area under the force curve. Note that the impulse

is the only property of

that influences the macroscopic behaviour of the system. If the contact time 2

is small, we’ll neglect it and write

∆p = I

Assuming that

only acts on a negligible amount of time

, all that matters to

us is its integral I, i.e. the area under the force curve.

wlog, assume

= 0 for easier mathematical treatment. We can consider a

family of functions D(t; ε) such that

lim

ε→0

D(t; ε) = 0 for all t 6= 0;

lim

ε→0

∞

−∞

D(t; ε) dt = 1.

So we can replace the force in our example by

(

;

), and then take the limit

as ε → 0.

For example, we can choose

D(t; ε) =

√

−t

/ε

ε = 1

ε = 0.5

This has height O(1/ε) and width O(ε).

It can be checked that this satisfies the properties listed above. Note that as

ε → 0, D(0; ε) → ∞. Therefore lim

ε→0

D(0; ε) does not exist.

Definition (Dirac delta function). The Dirac delta function is defined by

δ(x) = lim

ε→0

D(x; ε)

on the understanding that we can only use its integral properties. For example,

when we write

∞

−∞

g(x)δ(x) dx,

we actually mean

lim

ε→0

∞

−∞

g(x)D(x; ε) dx.

In fact, this is equal to g(0).

More generally,

(

)

(

x − c

) d

(

) if

c ∈

(

a, b

) and 0 otherwise,

provided g is continuous at x = c.

This gives a convenient way of representing and making calculations involving

impulsive or point forces. For example, in the previous example, we can write

m¨x = −mg + Iδ(t − T ).

Example.

−y

= 3

(

x −

) with

= 0 at

= 0

, π

. Note that our function

is split into two parts by x =

First consider the region 0

≤ x <

. Here the delta function is 0, and we have

−y

= 0 and

= 0 at

= 0. Then

−x

A sinh x

B cosh x

and

obtain

= 0 from the boundary condition. In the region

< x ≤ π

, we again

obtain

C sinh

(

π − x

) +

D cosh

(

π − x

) and (from the boundary condition),

D = 0.

When

, first insist that

is continuous at

. So

. Then

note that we have to solve

− y = 3δ



x −



But remember that the delta function makes sense only in an integral. So we

integrate both sides from

−

. Then we obtain

]

−

y dx = 3

Since we assume that y is well behaved, the second integral is 0. So we are left

with

]

−

= 3

So we have

−C cosh

− A cosh

= 3

A = C =

−3

2 cosh

Then we have

y =

(

−3 sinh x

2 cosh

0 ≤ x <

−3 sinh(π−x)

2 cosh

< x ≤ π

Note that at

, our final function has continuous

, discontinuous

and

infinite

. In general, differentiating a function makes it less continuous. This

is why we insisted at first that

has to be continuous. Otherwise,

would look

like a delta function, and y

would be something completely unrecognizable.

Hence the discontinuity is always addressed by the highest order derivative

since differentiation increases the discontinuity.

5.7 Heaviside step function

Definition (Heaviside step function). Define the Heaviside step function as:

H(x) =

−∞

δ(t) dt

We have

H(x) =











0 x < 0

1 x > 0

undefined x = 0

By the fundamental theorem of calculus,

= δ(x)

But remember that these functions and relationships can only be used inside

integrals.

6 Series solutions

Often, it is difficult to solve a differential equation directly. However, we can

attempt to find a Taylor series for the solution.

We will consider equations of the form

p(x)y

+ q(x)y

+ r(x)y = 0.

Definition (Ordinary and singular points). The point

is an ordinary

point of the differential equation if

and

have Taylor series about

(i.e. are

“analytic”, cf. Complex Analysis). Otherwise, x

is a singular point.

If x

is a singular point but the equation can be written as

P (x)(x − x

)

+ Q(x)(x − x

+ R(x)y = 0,

where

and

have Taylor series about

, then

is a regular singular point.

Example.

(i)

− x

)

−

+ 2

= 0.

= 0 is an ordinary point. However,

are (regular) singular points since p(±1) = 0.

(ii) sin xy

cos xy

+ 2

= 0.

nπ

are regular singular points while all

others are ordinary.

(iii)

(1 +

√

)

−

+ 2

= 0.

= 0 is an irregular singular point because

√

x is not differentiable at x = 0.

It is possible to show that if

is an ordinary point, then the equation is

guaranteed to have two linearly independent solutions of the form

y =

∞

n=0

(x − x

)

i.e. Taylor series about

. The solution must be convergent in some neighbour-

hood of x

is a regular singular point, then there is at least one solution of the form

y =

∞

n=0

(x − x

)

n+σ

with

= 0 (to ensure

is unique). The index

can be any complex number.

This is called a Frobenius series.

Alternatively, it can be nice to think of the Frobenius series as

y = (x − x

)

∞

n=0

(x − x

)

= (x − x

)

f(x)

where f (x) is analytic and has a Taylor series.

We will not prove these results, but merely apply them.

Ordinary points

Example. Consider (1

− x

)

−

+ 2

= 0. Find a series solution about

x = 0 (which is an ordinary point).

We try

∞

n=0

. First, we write the equation in the form of an

equidimensional equation with polynomial coefficients by multiplying both sides

. This little trick will make subsequent calculations slightly nicer. We

obtain

(1 − x)

) − 2x

(xy

) + 2x

y = 0

[(1 − x

)n(n − 1) − 2x

n + 2x

= 0

[n(n − 1) + (−n

− n + 2)x

= 0

We look at the coefficient of

and obtain the following general recurrence

relation:

n(n − 1)a

+ [−(n − 2)

− (n − 2) + 2]a

n−2

= 0

n(n − 1)a

= (n

− 3n)a

n−2

Here we do not divide by anything since they might be zero.

First consider the case

= 0. The left hand side gives 0

· a

= 0 (the right

hand side is 0 since

n−2

= 0). So any value of

satisfies the recurrence

relationship, and it can take any arbitrary value. This corresponds to a constant

of integration. Similarly, by considering n = 1, a

is arbitrary.

For n > 1, n and n − 1 are non-zero. So we have

n − 3

n − 1

n−2

In this case (but generally not), we can further simplify it to obtain:

n − 3

n − 1

n − 5

n − 3

n−4

n − 5

n − 1

n−4

−1

2k − 1

2k+1

= 0.

So we obtain

y = a

[1 −

−

− ···] + a

= a



1 −



1 + x

1 − x



+ a

Notice the logarithmic behaviour near

1 which are regular singular points.

Regular singular points

Example. Consider 4

+ 2(1

− x

)

− xy

= 0. Note that

= 0 is a singular

point. However, if we multiply throughout by

to obtain an equidimensional

equation, we obtain

4(x

) + 2(1 − x

)xy

− x

y = 0.

Since

1−x

and

−

both have Taylor series,

= 0 is a regular singular

point. Try

y =

∞

n=0

n+σ

with a

6= 0.

Substituting in, we have

n+σ

[4(n + σ)(n + σ − 1) + 2(1 − x

)(n + σ) − x

]

By considering the coefficient of

n+σ

, we obtain the general recurrence relation

[4(n + σ)(n + σ − 1) + 2(n + σ)]a

− [2(n − 2 + σ) + 1]a

n−2

= 0.

Simplifying the equation gives

2(n + σ)(2n + 2σ − 1)a

= (2n + 2σ − 3)a

n−2

The n = 0 case gives the indicial equation for the index σ:

2σ(2σ − 1)a

= 0.

Since

= 0, we must have

= 0 or

. The

= 0 solution corresponds to an

analytic (“Taylor series”) solution, while

corresponds to a non-analytic

one.

When σ = 0, the recurrence relation becomes

2n(2n − 1)a

= (2n − 3)a

n−2

When

= 0, this gives 0

· a

= 0. So

is arbitrary. For

n >

0, we can divide

and obtain

2n − 3

2n(2n − 1)

n−2

We can see that a

= 0 and so are subsequent odd terms.

If n = 2k, i.e. n is even, then

4k − 3

4k(4k − 1)

2k−2

y = a



1 +

4 · 3

8 · 7 · 4 ·3

+ ···



Note that we have only found one solution in this case.

Now when σ =

, we obtain

(2n + 1)(2n)a

= (2n − 2)a

n−2

When

= 0, we obtain 0

· a

= 0, so

is arbitrary. To avoid confusion with

the a

above, call it b

instead.

When

= 1, we obtain 6

= 0 and

= 0 and so are subsequent odd terms.

For even n,

n − 1

n(2n + 1)

n−2

y = b

1/2



1 +

2 · 5

2 · 5 · 4 ·9

+ ···



Resonance of solutions

Note that the indicial equation has two roots

, σ

. Consider the two different

cases:

(i)

− σ

is not an integer, then there are two linearly independent

Frobenius solutions

y =

(x − x

)

∞

n=0

(x − x

)

(x − x

)

∞

n=0

(x − x

)

As x → x

, y ∼ (x − x

)

, where Re(σ

) ≤ Re(σ

)

(ii)

−σ

is an integer (including when they are equal), there is one solution

of the form

= (x − x

)

∞

n=0

(x − x

)

with σ

≥ σ

In this case,

will not give a valid solution, as we will later see.

Instead, the other solution is (usually) in the form

= ln(x − x

∞

n=0

(x − x

)

n+σ

This form arises from resonance between the two solutions. But if the

resonance somehow avoids itself, we can possibly end up with two regular

Frobenius series solutions.

We can substitute this form of solution into the differential equation to

determine b

Example. Consider

− xy

= 0.

= 0 is a regular singular point. It is

already in equidimensional form (x

) − x(y) = 0. Try

y =

∞

n=0

n+σ

with a

6= 0. We obtain

n+σ

[(n + σ)(n + σ − 1) − x] = 0.

The general recurrence relation is

(n + σ)(n + σ − 1)a

= a

n−1

n = 0 gives the indicial equation

σ(σ − 1) = 0.

Then

= 0

1. We are guaranteed to have a solution in the form

= 1. When

σ = 1, the recurrence relation becomes

(n + 1)na

= a

n−1

When n = 0, 0 · a

= 0 so a

is arbitrary. When n > 0, we obtain

n(n + 1)

n−1

(n + 1)(n!)

= a



1 +

144

+ ···



When σ = 0, we obtain

n(n − 1)a

= a

n−1

When

= 0, 0

· a

= 0 and

is arbitrary. When

= 1, 0

· a

. However,

= 0 by our initial constraint. Contradiction. So there is no solution in this

form (If we ignore the constraint that

= 0, we know that

is arbitrary. But

this gives exactly the same solution we found previously with σ = 1)

The other solution is thus in the form

= y

ln x +

∞

n=0

7 Directional derivative

7.1 Gradient vector

Consider a function

(

x, y

) and a displacement ds = (d

). The change in

f(x, y) during that displacement is

df =

∂f

∂x

dx +

∂f

∂y

We can also write this as

df = (dx, dy) ·



∂f

∂x

∂f

∂y



= ds · ∇f

where

∇f

gradf



∂f

∂x

∂f

∂y



are the Cartesian components of the gradient of

We write ds =

s ds, where |

s| = 1. Then

Definition (Directional derivative). The directional derivative of

in the

direction of

s is

s · ∇f.

Definition (Gradient vector). The gradient vector

∇f

is defined as the vector

that satisfies

= ˆs · ∇f.

Officially, we take this to be the definition of

∇f

. Then

∇f



∂f

∂x

∂f

∂y



is a

theorem that can be proved from this definition.

We know that the directional derivative is given by

= ˆs · ∇f = |∇f|cos θ

where

is the angle between the displacement and

∇f

. Then when

cos θ

maximized,

= |∇f|. So we know that

(i) ∇f

has magnitude equal to the maximum rate of change of

(

x, y

) in the

xy plane.

(ii) It has direction in which f increases most rapidly.

(iii)

If ds is a displacement along a contour of

(i.e. along a line in which

constant), then

= 0.

So ˆs · ∇f = 0, i.e. ∇f is orthogonal to the contour.

7.2 Stationary points

There is always (at least) one direction in which

= 0, namely the direction

parallel to the contour of f. However, local maxima and minima have

= 0

for all directions, i.e. ˆs · ∇f = 0 for all ˆs, i.e. ∇f = 0. Then we know that

∂f

∂x

∂f

∂y

= 0.

However, apart from maxima and minima, in 3 dimensions, we also have saddle

points:

In general, we define a saddle point to be a point at which

∇f

= 0 but is not a

maximum or minimum.

When we plot out contours of functions, near maxima and minima, the

contours are locally elliptical. Near saddle points, they are locally hyperbolic.

Also, the contour lines cross at and only at saddle points.

7.3 Taylor series for multi-variable functions

Suppose we have a function

(

x, y

) and a point x

. Now consider a finite

displacement δs along a straight line in the x, y plane. Then

δs

= δs · ∇

The Taylor series along the line is

f(s) = f(s

+ δs)

= f(s

) + δs

(δs)

= f(s

) + δs · ∇f +

δs

(ˆs · ∇)(ˆs ·∇)f.

We get

δs · ∇f = (δx)

∂f

∂x

+ (δy)

∂f

∂y

= (x − x

)

∂f

∂x

+ (y − y

)

∂f

∂y

and

δs

(ˆs · ∇)(ˆs ·∇)f = (δs · ∇)(δs · ∇)f



δx

∂

∂x

+ δy

∂

∂y



δx

∂f

∂x

+ δy

∂f

∂y



= δx

∂

∂x

+ 2δxδy

∂

∂x∂y

+ δy

∂

∂y



δx δy







δx

δy



Definition (Hessian matrix). The Hessian matrix is the matrix

∇∇f =





In conclusion, we have

f(x, y) = f(x

, y

) + (x − x

+ (y − y

[(x − x

)

+ 2(x − x

)(y − y

+ (y − y

)

]

In general, the coordinate-free form is

f(x) = f(x

) + δx · ∇f(x

) +

δx · ∇∇f · δx

where the dot in the second term represents a matrix product. Alternatively, in

terms of the gradient operator (and real dot products), we have

f(x) = f(x

) + δx · ∇f(x

) +

[∇(∇f · δx)] · δx

7.4 Classification of stationary points

At a stationary point x

, we know that

∇f

) = 0. So at a point near the

stationary point,

f(x) ≈ f (x

) +

δx · H · δx,

where H = ∇∇f(x

) is the Hessian matrix.

At a minimum, Every point near x

has

(x)

> f

), i.e.

·H ·δ

0 for

all δx. We say δx ·H · δx is positive definite.

Similarly, at a maximum,

· H · δ

0 for all

x. We say

· H · δ

x is

negative definite.

At a saddle,

· H · δ

x is indefinite, i.e. it can be positive, negative or zero

depending on the direction.

This, so far, is not helpful, since we do not have an easy way to know what

sign

· H · δ

x could be. Now note that

∇∇f

is symmetric (because

). So

can be diagonalized (cf. Vectors and Matrices). With respect

to these axes in which H is diagonal (principal axes), we have

δx · H · δx = (δx, δy, ··· , δz)



















δx

δy

δz







= λ

(δx)

+ λ

(δy)

+ ··· + λ

(δz)

where λ

, λ

, ···λ

are the eigenvalues of H.

So for

· H · δ

x to be positive-definite, we need

0 for all

. Similarly,

it is negative-definite iff

0 for all

. If eigenvalues have mixed sign, then it

is a saddle point.

Finally, if there is at least one zero eigenvalue, then we need further analysis

to determine the nature of the stationary point.

Apart from finding eigenvalues, another way to determine the definiteness is

using the signature.

Definition (Signature of Hessian matrix). The signature of

is the pattern of

the signs of the subdeterminants:

|{z}



| {z }

, ··· ,



··· f



| {z }

|=|H|

Proposition.

is positive definite if and only if the signature is +

, ··· ,

is negative definite if and only if the signature is

−,

, ··· ,

(

−

. Otherwise,

H is indefinite.

7.5 Contours of f(x, y)

Consider

in 2 dimensions, and axes in which

is diagonal. So



0 λ



Write x − x

= (X, Y ).

Then near x

= constant

⇒

x = constant, i.e.

= constant.

At a maximum or minimum,

and

have the same sign. So these contours

are locally ellipses. At a saddle point, they have different signs and the contours

are locally hyperbolae.

Example. Find and classify the stationary points of

(

x, y

) = 4

−

+ 10y + 6. We have

= 12x

− 12y

= −12x + 2y + 10

= 24x

= −12

= 2

At stationary points, f

= f

= 0. So we have

12x

− 12y = 0, −12x + 2y + 10 = 0.

The first equation gives

. Substituting into the second equation, we obtain

= 1

5 and

= 1

25 respectively. So the stationary points are (1

1) and (5

25)

To classify them, first consider (1

1). Our Hessian matrix



24 −12

−12 2



Our signature is

= 24 and

−

96. Since we have a +

, −

signature, this

an indefinite case and it is a saddle point.

At (5

25),



120 −12

−12 2



= 120 and

= 240

−

144 = 96.

Since the signature is +, +, it is a minimum.

To draw the contours, we draw what the contours look like near the stationary

points, and then try to join them together, noting that contours cross only at

saddles.

(5,25)

(1,1)

-2 0 2 4 6

-10

8 Systems of differential equations

8.1 Linear equations

Consider two dependent variables y

(t), y

(t) related by

˙y

= ay

+ by

+ f

(t)

˙y

= cy

+ dy

+ f

(t)

We can write this in vector notation by



˙y





a b

c d









Y = M Y + F. We can convert this to a higher-order equation by

¨y

= a ˙y

+ b ˙y

= a ˙y

+ b(cy

+ dy

+ f

) +

= a ˙y

+ bcy

+ d( ˙y

− ay

− f

) + bf

¨y

− (a + d) ˙y

+ (ad − bc)y

= bf

− df

and we know how to solve this. However, this actually complicates the equation.

So what we usually do is the other way round: if we have a high-order

equation, we can do this to change it to a system of first-order equations:

If ¨y + a ˙y + by = f , write y

= y and y

= ˙y. Then let Y =



˙y



Our system of equations becomes

˙y

= y

˙y

= f − ay

− by

Y =



0 1

−b −a





Now consider the general equation

Y = M Y + F

Y − MY = F

We first look for a complementary solution Y

= v

λt

, where v is a constant

vector. So we get

λv − Mv = 0.

We can write this as

Mv = λv.

So λ is the eigenvalue of M and v is the corresponding eigenvector.

We can solve this by solving the characteristic equation

det

(

M − λI

) = 0.

Then for each λ, we find the corresponding v.

Example.

Y =



−4 24

1 −2



Y +





The characteristic equation of



−4 − λ 24

1 −2 − λ



= 0, which gives (

8)(λ − 2) = 0 and λ = 2, −8.

When λ = 2, v satisfies



−6 24

1 −4





= 0,

and we obtain v





When λ = −8, we have



4 24

1 6





= 0,

and v



−6



. So the complementary solution is

Y = A





+ B



−6



−8t

To plot the phase-space trajectories, we first consider the cases where Y is

an eigenvector. If Y is an integer multiple of





, then it will keep moving

outwards in the same direction. Similarly, if it is an integer multiple of



−6



it will move towards the origin.

= (4, 1)

= (−6, 1)

We can now add more (curved) lines based on these two:

= (4, 1)

= (−6, 1)

To find the particular integral, we try Y

= ue

. Then





−



−4 24

1 −2











5 −24

−1 3













= −



3 24

1 5







−4

−1



So the general solution is

Y = A





+ B



−6



−8t

−





In general, there are three possible cases of

Y corresponding to three

different possible eigenvalues of M:

(i)

If both

, λ

are real with opposite sign (

0). wlog assume

Then there is a saddle as above:

= (4, 1)

= (−6, 1)

(ii)

, λ

are real with

0. wlog assume

|λ

| ≥ |λ

. Then the phase

portrait is

If both

, λ

0, then the arrows point towards the intersection and we

say there is a stable node. If both are positive, they point outwards and

there is an unstable node.

(iii) If λ

, λ

are complex conjugates, then we obtain a spiral

(

1,2

)

0, then it spirals inwards. If

(

1,2

)

0, then it spirals

outwards. If

(

1,2

) = 0, then we have ellipses with common centers

instead of spirals. We can determine whether the spiral is positive (as

shown above), or negative (mirror image of the spiral above) by considering

the eigenvectors. An example of this is given below.

8.2 Nonlinear dynamical systems

Consider the second-order autonomous system (i.e.

does not explicitly appear

in the forcing terms on the right)

˙x = f(x, y)

˙y = g(x, y)

It can be difficult to solve the equations, but we can learn a lot about phase-space

trajectories of these solutions by studying the equilibria and their stability.

Definition (Equilibrium point). An equilibrium point is a point in which

˙x = ˙y = 0 at x

= (x

, y

Clearly this occurs when

(

, y

) =

(

, y

) = 0. We solve these simultane-

ously for x

, y

To determine the stability, write x = x

+ ξ, y = y

+ η. Then

ξ = f(x

+ ξ, y

+ η)

= f(x

, y

) + ξ

∂f

∂x

) + η

∂f

∂y

) + O(ξ

, η

)

So if ξ, η  1,



˙η









This is a linear system, and we can determine its character from the eigensolu-

tions.

Example. (Population dynamics - predator-prey system) Suppose that there

are x prey and y predators. Then we have the following for the prey:

˙x = αx

|{z}

births - deaths

− βx

|{z}

natural competition

− γxy

|{z}

killed by predators

and the following for the predators:

˙y = εxy

|{z}

birth/survival rate

− δy

|{z}

natural death rate

For example, let

˙x = 8x − 2x

− 2xy

˙y = xy − y

We find the fixed points: x(8 − 2x − 2y) = 0 gives x = 0 or y = 4 − x.

We also want y(x − 1) = 0 which gives y = 0 or x = 1.

So the fixed points are (0, 0), (4, 0), (1, 3).

Near (0, 0), we have



˙η





8 0

0 −1





We clearly have eigenvalues 8, −1 with the standard basis as the eigenvectors.

Near (4

0), we have

= 4 +

. Instead of expanding partial derivatives,

we can obtain from the equations directly:

ξ = (4 + ξ)(8 − 8 − 2ξ − 2η)

= −8ξ − 8η − 2ξ

− 2ξη

˙η = η(4 + ξ − 1)

= 3η + ξη

Ignoring the second-order terms, we have



˙η





−8 −8

0 3





The eigenvalues are −8 and 3, with associated eigenvectors (1, 0), (8, −11).

Near (1, 3), we have x = 1 + ξ, y = 3 + η. So

ξ = (1 + ξ)(8 − 2 − 2ξ − 6 − 2η)

≈ −2ξ − 2η

˙η = (3 + η)(1 + ξ − 1)

≈ 3η



˙η





−2 −2

3 0





The eigenvalues are

−

± i

√

. Since it is complex with a negative real part, it

is a stable spiral.

We can determine the chirality of the spirals by considering what happens to

a small perturbation to the right





with

ξ >

0. We have



˙η





−2ξ

3ξ



. So

x will head top-left, and the spiral is counter-clockwise (“positive”).

We can patch the three equilibrium points together and obtain the following

phase portrait:

0 1 2 3 4 5

We see that (1, 3) is a stable solution which almost all solutions spiral towards.

9 Partial differential equations (PDEs)

9.1 First-order wave equation

Consider the equation of the form

∂y

∂t

= c

∂y

∂x

with

a constant and

a function of

and

. This is known as the (first-order)

wave equation. We will later see that solutions correspond to waves travelling in

one direction.

We write this as

∂y

∂t

− c

∂y

∂x

= 0.

Recall that along a path x = x(t) so that y = y(x(t), t),

∂y

∂x

∂y

∂t

∂y

∂x

+ c

∂y

∂x

by the chain rule. Now we choose a path along which

= −c. (1)

Along such paths,

= 0 (2)

So we have replaced the original partial differential equations with a pair of

ordinary differential equations.

Each path that satisfies (1) can be described by

− ct

, where

is a

constant. We can write this as x

= x + ct.

From (2), along each of these paths,

is constant. So suppose for each

, the

value of

along the path

is given by

(

), where

is an arbitrary

function. Then the solution to the wave equation is

y = f(x + ct),

By differentiating this directly, we can easily check that every function of this

form is a solution to the wave equation.

The contours of y look rather boring.

contours of y

Note that as we move up the time axis, we are simply taking the

= 0 solution

and translating it to the left.

The paths we’ve identified are called the “characteristics” of the wave equation.

In this particular example,

is constant along the characteristics (because the

equation is unforced).

We usually have initial conditions e.g.

y(x, 0) = x

− 3

Since we know that y = f(x + ct) and f(x) = x

− 3, y must be given by

y = (x + ct)

− 3.

We can plot the xy curve for different values of t to obtain this:

We see that each solution is just a translation of the t = 0 version.

We can also solve forced equations, such as

∂y

∂t

+ 5

∂y

∂x

= e

−t

, y(x, 0) = e

−x

Along each path

x −

, we have

−t

. So

(

)

− e

−t

for some

function f .

Using the boundary conditions provided, At

= 0,

(

)

−

1 and

So f (x

) − 1 = e

−x

, i.e. f(x

) = 1 + e

−x

. So

y = 1 + e

−(x−5t)

− e

−t

9.2 Second-order wave equation

We consider equations in the following form:

∂

∂t

= c

∂

∂x

This is the second-order wave equation, and is often known as the “hyperbolic

equation” because the form resembles that of a hyperbola (which has the form

− b

= 0). However, the differential equation has no connections to

hyperbolae whatsoever.

This equation models an actual wave in one dimension. Consider a horizontal

string, along the

axis. We let

(

x, t

) be the vertical displacement of the string

at the point x at time t.

Suppose that

(

) is the mass per unit length of a string. Then the restoring

force on the string

∂

∂t

is proportional to the second derivative

∂

∂x

. So

we obtain this wave equation.

(Why is it proportional to the second derivative? It certainly cannot be

proportional to

, because we get no force if we just move the whole string

upwards. It also cannot be proportional to

∂y/∂x

: if we have a straight slope,

then the force pulling upwards is the same as the force pulling downwards, and

we should have no force. We have a force only if the string is curved, and

curvature is measured by the second derivative)

To solve the equation, suppose that c is constant. Then we can write

∂

∂t

− c

∂

∂x

= 0



∂

∂t

+ c

∂

∂x



∂

∂t

− c

∂

∂x



y = 0

(

), then the first operator differentiates it to give a constant (as in

the first-order wave equation). Then applying the second operator differentiates

it to 0. So y = f (x + ct) is a solution.

Since the operators are commutative,

(

x − ct

) is also a solution. Since

the equation is linear, the general solution is

y = f(x + ct) + g(x − ct).

This shows that the solution composes of superpositions of waves travelling to

the left and waves travelling to the right.

We can show that this is indeed the most general solution by substituting

and

x −ct

. We can show, using the chain rule, that

−c

≡

−4c

ηξ

= 0. Integrating twice gives y = f (ξ) + g(η).

How many boundary conditions do we need to have a unique solution? In

ODEs, we simply count the order of the equation. In PDEs, we have to count

over all variables. In this case, we need 2 boundary conditions and 2 initial

conditions. For example, we can have:

– Initial conditions: at t = 0,

y =

1 + x

∂y

∂t

= 0

– Boundary conditions: y → 0 as x → ±∞.

We know that the solution has the form

y = f(x + ct) + g(x − ct).

The first initial condition give

f(x) + g(x) =

1 + x

(1)

The second initial condition gives

∂y

∂t

= cf

(x) − cg

(x) = 0 (2)

From (2), we know that

. So

and

differ by a constant. wlog, we can

assume that they are indeed equal, since if we had, say,

+ 2, then we can

let

= (

(

) + 1) + (

(

x − ct

) + 1) instead, and the new

and

are equal.

From (1), we must have

f(x) = g(x) =

2(1 + x

)

So, overall,

y =



1 + (x + ct)

1 + (x − ct)



Where we substituted x for x + ct and x − ct in f and g respectively.

9.3 The diffusion equation

Heat conduction in a solid in one dimension is modelled by the diffusion equation

∂T

∂t

= κ

∂

∂x

This is known as a parabolic PDE (parabolic because it resembles y = ax

Here T (x, t) is the temperature and the constant κ is the diffusivity.

Example. Consider an infinitely long bar heated at one end (

= 0). Note

in general the “velocity”

∂T /∂t

is proportional to curvature

∂

T/∂x

, while in

the wave equation, it is the “acceleration” that is proportional to curvature.

In this case, instead of oscillatory, the diffusion equation is dissipative and all

unevenness simply decays away.

Suppose

(

0) = 0,

, t

) =

(

) =

(

0 t < 0

1 t > 0

, and

(

x, t

)

→

0 as

x → ∞

. In words, this says that the rod is initially cool (at temperature 0), and

then the end is heated up on one end after t = 0.

There is a similarity solution of the diffusion equation valid on an infinite

domain (or our semi-infinite domain) in which

(

x, t

) =

(

), where

√

κt

Applying the chain rule, we have

∂T

∂x

dθ

dη

∂η

∂x

√

κt

(η)

∂

∂x

√

κt

dθ

dη

∂η

∂x

4κt

(η)

∂T

∂t

dθ

dη

∂η

∂t

= −

√

3/2

(η)

= −

(η)

Putting this into the diffusion equation yields

−

= κ

4κt

+ 2ηθ

= 0

This is an ordinary differential equation for

(

). This can be seen as a first-

order equation for

with non-constant coefficients. Use the integrating factor

µ = exp(

2η dη) = e

. So

)

= 0

= Ae

−η

θ = A

−u

du + B

= α erf(η) + B

where erf(η) =

√

−u

du from statistics, and erf(η) → 1 as η → ∞.

Now look at the boundary and initial conditions, (recall

√

κt

)) and

express them in terms of η. As x → 0, we have η → 0. So θ = 1 at η = 0.

Also, if x → ∞, t → 0

, then η → ∞. So θ → 0 as η → ∞.

(0) = 1

⇒ B

= 1. Colloquially,

(

∞

) = 0 gives

−

1. So

= 1

−erf

(

This is also sometimes written as

erfc

(

), the error function complement of

. So

T = erfc



√

κt



In general, at any particular fixed time t

, T (x) looks like

with decay length

(

√

κt

). So if we actually have a finite bar of length

, we

can treat it as infinite if

√

κt  L, or t  L

/κ