3 Modelling with Ordinary Differential Equations

Question 1

When are evolutionary processes modelled with ODEs?

Answer 1

Evo process assumed to happen continously in time, when one single species involved modelled with an ** autonomous ODE**, referred to as scalar ODE
with the form
dx/dt ≡ x* = f(x(t))

Given f(x)∈ C¹
usually nonlinear
x(t)#relative # individuals in a pop
initial condition x(0)=x₀

Question 2

species interacting

Answer 2

May be predator,prey models we will consider a set of coupled ODEs

evolution of two interacting species can be modelled by
x₁ = f₁(x₁(t), x₂(t)) and
x₂ = f₂(x₁(t), x₂(t)),

where x₁, x₂ represent #individuals of species 1,2
(e.g. fractions of individuals, or densities of each in the pop)
given C¹ functs f₁(x₁, x₂) and f2(x₁, x₂)
initial condition (x₁(0), x₂(0)) is assumed to be known

Question 3

assumptions for modelling with ODES

Answer 3

(i) evolution is asynchronous, with reproduction/death occurring at different times and generations can therefore overlap;
(ii) populations are assumed to be large enough and
spatially-homogeneous so that fluctuations and space can be ignored;
(iii) genetics and mating are ignored.

*sol exists and is unique, only accept physical sols (>0)
* Will a given species go extinct?
* Under which circumstances do interacting species coexist?

Question 4

Verhulst equation

Answer 4

Continuous time example of the logistic map (Now we can solve)
the exponential growth of the pop size N(t) at time t is limited by a carrying capacity K
N*=rN(1- N/k)
where r>0 is the rare of growth limited by carrying capacity K>0.

with 0 < N(0) < K

Growth limiting term~ -rN²/k

Question 5

Verhulst equation

dimensions

N*=rN(1- N/k)
where r>0 is the rare of growth limited by carrying capacity K>0.

Answer 5

[r]=1/[T]=inverse of time;

N and K are numbers=>
[N]=[K]=1

Question 6

N*=rN(1- N/k)

without growth limitation:

Answer 6

N*=rN

⇒ N(t) = N(0)eʳᵗ → ∞as t→∞

exponential growth as t increases it explodes
(unrealistic!)

Question 7

N * =rN(1- N/k)
Useful also to non-dimensionalise the equation by changing time:

Answer 7

t → 𝜏 ≡ rt

where 𝜏 is dimensionless
: [𝜏 ] = [r] [t]
= (1/time) . time = 1
as r is a rate

so dx/dt
= rx(t)(1-x(t))
becomes

dx/d𝜏 =
dx/dt . dt/d𝜏
= x(𝜏)(1-x(𝜏))

Question 8

verhulst eq x = N/K ≥ 0 giving the density of
the population relative to the carrying capacity

Answer 8

Diving both sides of by K readily
gives
dx/dt = rx (1 − x),

where the left- (LHS) and right-hand-side (RHS) have dimension
1/T = [r].

Question 9

Solution of the Verhulst equation
Q4 Example sheet 1

Answer 9

dx/dt=x*=f(x(t)) is separable

rearrange and integrate both sides,partial fractions
∫ₓ_₀ ˣ dx~/f(x~)= ∫₀ ᵗ dt~= t

solving for tau then t with IC

x(t) = x₀eʳᵗ/(1 + x₀(eʳᵗ − 1)

sol x(t) s shaped approaches 1 as t tends to infinity, initially exp growth
growth slows
carrying capacity

Question 10

Linear stability analysis of Verhulst’s model

Answer 10

two equilibria
x∗₀=0 extinction
x∗₁=1 N=Kx reaches carrying capacity K
look at stability

f′(x) = r(1 − 2x)
⇒ f′(x∗0) = r > 0,
f′(x∗1) = −r < 0

⇒ confirms that x∗
0 = 0 is unstable whereas x∗
1 = 1 is asymptotically stable

Question 11

equilibria sols of ODES

Answer 11

equilibrium points of ODEs, denoted by x∗, are the points where the solution does neither increase nor decrease. For scalar ODEs, equilibria are thus
obtained by demanding x* = 0,

usually only look at FP not sols as usually cant solve

Question 12

Linear stability analysis: General Method

Answer 12

To determine the stability of a fixed point x∗ and the dynamics around it, we write
x = x∗ + y,
where y is small deviation that changes in time.
subbed and taylor expand about x∗
y* = f(x∗ + y) = f(x∗) + yf′(x∗) + higher order terms,
( f(x∗)=0 since SS)

where f′ = (d/dy)f = (d/dx)f.

To linear order we find that the deviation y around x∗ changes in time according to
y* = f′(x∗)y
This has sol by integrating
y(t)=y(0)e⁽ᶠ′⁽ˣ∗⁾⁾ᵗ
so when t → ∞
(if f′(x∗) < 0) :
y(t) → 0 vanishes in this case x(t) → x∗. STABLE EQ

(If f′(x∗) > 0:
y(t) → ∞ (grows/explodes) In this case x(t) deviates greatly from x∗ UNSTABLE EQ

Question 13

linear stability analysis conclusions

Answer 13

**hyperbolic equilibrium **
x∗ is asymptotically stable and approached exponentially if f′(x∗) < 0,

and

unstable when f′(x∗) > 0.

When f′(x∗) = 0 the fixed point is non-hyperbolic and no conclusion can be drawn from the linear stability analysis (other methods are necessary).

Question 14

differences with linear stability analysis of ODEs and scalar maps

Answer 14

Methodology is similar:
*generally nonlinear model cannot be solved => linearization about
physical fixed points/equilibria.
*in both cases y either increases or decreases exponentially,

but stability of x∗ is here in ODEs determined by the sign of f′(x∗) (whereas
in Scalar maps it is whether MODULUS|f′(x∗)| > 1 or |f′(x∗)| < 1.

Question 15

summary
differences with linear stability analysis of ODEs and scalar maps

Answer 15

• Fixed point x* for a scalar (1D) map is given by x=f(x).
• Equilibrium x* for a scalar (1D) ODE is f(x)=0
  For both maps and ODEs, stability of hyperbolic fixed point/equilibrium x is determined by
  the eigenvalues of the Jacobian (at x), which is just the value of f ‘(x) in the 1D case.

However stability conditions differ:
- Hyperbolic fixed point x* of a scalar map is asymptotically stable if |f ‘(x)|<1 and unstable if
|f ‘(x)|>1
- Hyperbolic equilibrium x* of a scalar ODE is asymptotically stable if f ‘(x)<0 and unstable if
f ‘(x)>0.

These notions, with their similarities and differences, are generalized to planar maps and ODEs
For both maps and ODEs, linear stability analysis doesn’t say anything about stability of non-hyperbolic fixed points/ equilibria.
However:
- x* is a non-hyperbolic fixed point of a scalar map if |f ‘(x)|=1
- x is a non-hyperbolic fixed point of a scalar ODE if f ‘(x*)=0

Question 16

Hartman-Grobman theorem (Thm. 3.1) e

Answer 16

linear stability analysis not limited to linear regime

hyperbolic
fixed point x∗ is asymptotically stable if f′(x∗) < 0 and is unstable if f′(x∗) > 0

Question 17

Direction field and phase portrait:

Answer 17

Useful qualitative information by noting that f(x(t)) gives rate of change of x at each t =>
- Evaluating f(x(t)) at each t and gives the slope of solution x(t) at each point (t,x)
- Graph of x(t) vs t forms the “direction field”
- Connecting vectors of direction field from => phase portrait giving solutions of x*=f(x(t))

DIRECTION FIELD obtained by evaluating slope dx/dt=f(x) of so, x (t) at each point (t,x)
idea is to evaluate f(x(t)) at each point of your phase space and see how sol changes
using the fact that f(x(t)) gives the slope of sol at any point x(t)~ the rate of change of x at each t

PHASE PORTRAIT
A family of solution curves from the direction field forms the phase portrait by connecting them and choosing an IC

Question 18

PHASE PORTAIT FOR
Verhulst model:
f(x(t)) = rx(t)(1 − x(t))
With 0< x₀<1, f(x)> 0 for 0<x<1 and f(x=0)=f(x=1)=0

Answer 18

we have
f(x) ≥ 0 with f(0) = f(1) = 0
hence
dx/dt>0 for all x ̸= 0,

DIRECTION FIELD
0< x₀<1 and 0 < x < 1
tells us f(x) always positive so
slope points up and right for 0 < x < 1 and becomes flat near x = 0 and x = 1, hence yielding the S-shape curve

it is clear from the direction field that solutions from different initial conditions never cross (otherwise the solution would no longer be unique

DESCRIPTION
Hence, the solution x(t) starting from 0< x₀<1 grows as f(x(t))>0, and keeps growing until it approaches
due to the growth limiting term

The direction field consists of vectors (x,f(x)) pointing up right and becoming horizontal near x=1 (stable equilibrium)
and x=0 (unstable). The solution x(t) hence connects x(0) to x=1 by forming an S-shaped function.

Question 19

J
∗ ≡ (∂fi/∂xj )|x=x∗ is the Jacobian of f at x
∗

hyperbolic?

Answer 19

The equilibrium x∗is hyperbolic if J ∗has no eigenvalues with zero real parts. An equilibrium point x∗
is non-hyperbolic when the Jacobian matrix J∗
has at least one eigenvalue
that is purely imaginary

Question 20

Definition 3.1 (Equilibrium, Jacobian, hyperbolic / non-hyperbolic equilibria). A

Answer 20

An equilibrium point (or “equilibrium”, “rest point”, “fixed point” or “stationary point”) of
x˙ = f(x(t)), denoted by x∗ = (x∗_1, . . . , x∗_n)^T, is a point where x does neither increase
nor decrease. An equilibrium point therefore satisfies f(x∗) = 0. Equilibria are found by
solving simultaneously the n equations: f_1(x∗) = · · · = f_n(x∗) = 0.
J∗ ≡ (∂fi/∂xj )|x=x∗ is the Jacobian of f at x∗

Question 21

n-dimensional ODE systems

Answer 21

dx(t)/dt ≡ x*(t) = f(x(t))

where f = (f_1, . . . , f_n)
T : R^n → R^n is a C^1 vector field and
x = (x_1, . . . , x_n)^T

systems of ODEs are often nonlinear and cannot be solved

Question 22

Modelling with scalar ODEs: the Verhulst logistic model

Answer 22

with n = 1, f(x(t)) = rx(t){1 − x(t)}

1D ODE

Question 23

For an n-dimensional ODE systems with suitable given initial condition x(0)

how to find x∗

Answer 23

Cauchy Lipschitz Theorem => existence and uniqueness of solution under general conditions

How to find?

Look for physical equilibria x∗
linearise about ** x* ** = f(x(t))

Equilibria x∗ =(x₁,…,xₙ)^T need to be physical
if each x_i describe a density of a population= then all x_i ≥ 0

obtained by solving
f(x∗)=0
f₁( x∗)= ….,=fₙ( x∗)=0

. For hyperbolic equilibria, describe by linear system
y* = J∗ y

where J∗ = (∂fᵢ/∂xⱼ )|ₓ₌ₓ∗ is the Jacobian (n × n matrix) of f at x∗
and y=x-x∗ small deviation

Question 24

hyperbolic equilibrium

Answer 24

x∗ is an hyperbolic equilibrium if has no purely imaginary eigenvalues, and non-hyperbolic otherwise

eigienvalues of J∗

meaning we can predict the dynamics about equilibrium from the jacobian

so linearisation is helpful to study dynamics around a hyperbolic equilibrium only

Question 25

Hartman-Grobman theorem (Thm. 3.1)

Answer 25

If x∗ is a hyperbolic fixed point of x˙ = f(x(t)),
the phase portrait of the nonlinear system is equivalent to that of the linearized system
y˙ = J ∗y in a neighbourhood of x∗(with y = x − x∗)

in the sense that there is a homeomorphism between the solutions of x˙ = f(x(t)) and y˙ = J∗y. Hence, there is a one-to-one continuous mapping between the solutions of both systems

Question 26

Why do we linearise

Answer 26

> useful to linearize around hyperbolic equilibria to understand the stability of :intuitively, it is stable if, starting close to it, the solution doesn’t deviate from it. Linearization is not helpful if is non-hyperbolic (does not allow to conclude on stability).

Question 27

planar ODE systems

Poincare-Bendixson theorem

Answer 27

n=2
result in Poincare-Bendixson theorem
It says that for n=2, the solutions of form cycles or terminate on an equilibrium

=> no chaos in 2D ODE systems!

(cyclic solution or terminate in equilibrium)
(We saw chaotic behaviour in the system for maps: simple scalar logistic map, when r>3.57…. not the case for planar ODEs)

Question 28

Theorem 3.2 (Poincar´e-Bendixson).

glossed over

Answer 28

. For x˙ = f(x(t)) and f = (f_1, f_2)^T, and x(t) =(x1(t), x2(t))^T. We consider a closed bounded region D. Suppose a solution H from a given point of the plane lies within D.

Then, one of the following is true: (i) H is a closed
trajectory, e.g. a limit cycle or a closed orbit; (ii) H asymptotically tends to a closed trajectory, e.g. a limit cycle; (iii) H terminates on a equilibrium point.

Proof: see, e.g., in P. Glendinning, “Stability, Instability & Chaos: An Introduction to the Theory
of Nonlinear Differential Equations” (Cambridge, 1994

The Poincar´e-Bendixson theorem therefore says that the dynamics predicted by (3.2) in
the plane is limited: Since different solutions of (3.2) cannot cross (uniqueness) the only
possibilities are either to approach an asymptotically stable fixed point or to settle on a periodic orbit (persisting oscillatory behaviour). Remarkably, Thm. 3.2 implies that the behaviour of (3.2) is never chaotic.)

Question 29

Theorem 3.4 (Local stability of nonlinear planar ODEs).

Question 30

3.3 Linear planar (two-dimensional) ODEs
-used for
-general form

Answer 30

models dynamics of two interacting species with ODES:
x˙₁ = f₁(x₁(t), x₂(t)) and
x˙₂ = f₂(x₁(t), x₂(t)),

LINEAR SO:
x˙₁ = ax₁+bx₂ and
x˙₂ = cx₁+dx₂
usually not solved- dynamics modelled
where x= (x₁ x₂)^T
A=
[a b]
[c d] is a non-singular matrix (det not 0 hence A-1 well defined)
x˙(t) = Ax(t), with x(0) given

A has eigenvalues λ±

clearly origin x∗ =(0,0) unique sol of Ax(t)=0

Question 31

Stability of x∗ = 0 and dynamics about it?

Answer 31

Formally: x(t) = eᴬᵗx(0) = Σₖ₌₀∞(At)ᵏ/k! x(0), which is difficult to compute directly

Using prev techniques: Choose matrix P
use similarity transformation
A and B are similar matrices with the same eigenvalues λ±
ξ(t) = P⁻¹x(t) with ξ(0) = P⁻¹x(0)
B=P⁻¹AP
where P is chosen so that computing powers of B is easy
ξ(t) *˙ = Bξ(t)

ξ(t)= eᴮᵗξ(0)
= Σₖ₌₀∞(Bt)ᵏ/k! ξ(0) easier to compute

x(t)= eᴬᵗx(0)
= P eᴮᵗ ξ(0)
= P eᴮᵗ P⁻¹ x(0)

thus
When t → ∞, |x(t)| and |ξ(t)| have same behaviour: either both vanish (0 asymptotically stable),or both → ∞ (0 unstable)

Question 32

If A has two distinct real eigenvalues

Answer 32

can choose s.t diagonal

If A has two distinct real eigenvalues, λ±, we can always choose P to diagonalize
the matrix B =
[λ₊ 0]
[0 λ₋]
In this case eᴮᵗ =
[exp(λ₊t) 0]
[0 exp(λ₋t)]
thus
eᴬᵗ =
[exp(λ₊t) 0]
P [0 exp(λ₋t)]P⁻¹

Question 33

if A has one real double eigenvalues

Answer 33

λ₊ =λ₋=λ
two eigenvectors are associated with it
if we can diagonalise:
B=
[λ 0]
[0 λ]

and thus
eᴮᵗ =
[exp(λt) 0]
[0 exp(λt)]

If A can’t be diagonalised choose P to obtain Jordan-Block form
B=
[λ 1]
[0 λ]

eᴮᵗ =
[exp(λt) texp(λt)]
[0 exp(λt)]

Question 34

if A has has two complex conjugate eigenvalues

Answer 34

λ±=α ± iβ (wlog β>0)
eᴮᵗ
=exp(αt) *
[cos βt sin βt
[− sin βt cos βt]

exp(αt) correspond to amplitude to either
grow in time α >0
decrease in time α <0
giving a spiral/focus

Question 35

Stability of x∗=0 and dynamics about it depends on
λ±

Answer 35

x∗ = (0, 0) is

• locally asymptotically stable when Re(λ+) < 0 and Re(λ−) < 0, since in this case |ξ(t)| → 0 and |x(t)| → 0 as t → ∞;
• locally unstable when Re(λ+) > 0 or/and Re(λ−) > 0, since in this case |ξ(t)| → ∞
  and |x(t)| → ∞ as t → ∞;
• a center when Re(λ±) = 0. In this case, the linear system (3.8) is characterized by closed orbits surrounding the origin, and |ξ(t)| = |ξ(0)|.

Question 36

reminder non hyperbolic when

Answer 36

x∗ is non-hyperbolic when Re(λ+)Re(λ−) = 0, and, as discussed below, the behaviour of a nonlinear system around a non-hyperbolic equilibrium cannot be inferred from its linearization (3.8). This means that a linear center around x ∗ may not result in a center with closed orbits in a system exhibiting nonlinearities (but it can also result in closed orbits).

predicts its a center but

Question 37

For planar systems, eigenvalues of A are given by

Answer 37

λ± =
[tr(A) ± sqrt((tr(A))2 − 4 det(A)) ]/2

which applies to planar systsme

Question 38

Theorem 3.3 (Stability of linear planar ODEs)

When is the origin asymptotically stable?

x˙(t) = Ax(t), with x(0) given

Answer 38

For the linear system , the following
statements are equivalent:
* The origin x∗ = (x∗1, x∗2)^T = (0, 0)T is an asymptotically stable equilibrium point.

• The eigenvalues λ+ and λ− have both a negative real part.
• det(A) = λ+λ− = ad − bc > 0 and tr(A) = λ+ + λ− = a + d < 0
  POSITIVE DET and NEGATIVE TRACE

det is non zero

Question 39

Dynamics about x*=(0,0): 6 possibilities summarized by the Jury’s conditions for planar ODEs:

Answer 39

not to be confused with the Jury’s conditions
for planar maps!

(a) The equilibrium x∗ is an attractor (asymptotically stable) if det(A) > 0 and tr(A) < 0,
i.e. if Re(λ±) < 0. We thus distinguish two types of behaviour:
(a.i) x∗ is a stable spiral (focus) if 4 det(A) > (tr(A))²
, in this case: Im(λ±) ̸= 0;
(a.ii) otherwise x∗is a stable node (λ± have no imaginary part: Im(λ±) = 0).

(b) The equilibrium x∗ is a repeller (unstable) if det(A) > 0 and tr(A) > 0, i.e. if Re(λ±) > 0. We thus distinguish two types of behaviour:
(b.i) x∗ is an unstable spiral (focus) if 4 det(A) > (tr(A))² , i.e. when Im(λ±) ̸= 0;
(b.ii) x∗is an unstable node otherwise (when Im(λ±) = 0).

(c) When det(A) < 0, we have λ+λ− < 0, the equilibrium x∗ is a saddle (unstable): one eigenvalue is negative and the other is positive; solutions deviate from x∗ along one
eigenvector and approach x∗along the other; resulting in deviating from x∗

(d) x∗ is a centre if tr(A) = 0 and det(A) > 0, when λ± are purely imaginary (Re(λ±) =0, Im(λ+) = −Im(λ−) ̸= 0). In this case, x∗ is non-hyperbolic. λ± are purely imaginary (complex conjugate) ⇒ Re(λ±) = 0

It has to be stressed that x∗ is a centre only for the dynamics of the linear ODEs, (3.8), but it may be a spiral (focus) or a center in nonlinear systems (depending on nonlinear terms).

Question 40

Similarities with Planar ODEs and linear planar maps

Answer 40

We note striking similarities with the analysis of linear planar maps of Sec. 2.5. These
stem from the similarity of the expressions of exp(Bt)
and B^t.

However, there are also important differences: the stability of ODEs equilibria depends on the sign of λ±, whereas the stability of fixed points of planar maps depends on whether |λ±| are greater or less than one

Remember that linear
stability analysis is
useful only for
hyperbolic equilibria,
but it does not say
ANYTHING about the
stability of non-hyperbolic
equilibria

Question 41

3.4 Nonlinear two-dimensional ODE systems: Linear stability analysis

x˙₁ = f₁(x₁(t), x₂(t)) and
x˙₂ = f₂(x₁(t), x₂(t)),

with given initial condition x₁(0), x₂(0) and where f₁, f₂ ∈ C¹

Answer 41

Cauchy-Lipschitz Thm ensures existence
and uniqueness. Poincare-Bendixson=>
reach an equilibrium or cyclic behaviour

Equilibirum: x∗ = (x∗₁, x∗₂)^T
and satisfies simultaneously
f₁(x∗)= f₁((x∗₁, x∗₂) = 0 and
f₂(x∗)= f₂((x∗₁, x∗₂) = 0
typically describe fractions of individuals of species:
x₁≥0, x₂≥0 for all t≥0
PHYSICAL

system can have 0,1 or many physical equilibria

LINEAR STABILITY ANALYSIS

Question 42

LINEAR STABILITY ANALYSIS for

3.4 Nonlinear two-dimensional ODE systems: Linear stability analysis

x˙₁ = f₁(x₁(t), x₂(t)) and
x˙₂ = f₂(x₁(t), x₂(t)),

with given initial condition x₁(0), x₂(0) and where f₁, f₂ ∈ C¹

Answer 42

Introduce:
y ≡ x − x∗
Look at deviations
y = (y₁, y₂)^T about equilibrium xc = (x∗₁, x∗₂)^T

How do y₁,=x₁- x∗₁ and y₂= x₂-x₂∗
change in time?
Substituting x ≡ y + x∗ into ODEs and Taylor-expand to linear order in y₁, y₂

y₁= x₁ = f₁(y₁ + x∗₁, y₂ + x∗₂)
= f₁(x∗₁,x∗₂) + y₁∂ᵧ₁ f₁(x∗₁,x∗₂) + y₂∂ᵧ₂ f₁(x∗₁,x∗₂) + higher order
=0 + …

y₂= x₂ = f₂(y₁ + x∗₁, y₂ + x∗₂)
= f₂(x∗₁,x∗₂) + y₁∂ᵧ₁ f₂(x∗₁,x∗₂) + y₂∂ᵧ₂ f₂(x∗₁,x∗₂) + higher order y’s
=0 + …

keeping only linear terms we have jacobian to describe system:

y˙(t) = J∗y(t)
where
J∗ =
[ (∂f1/∂x1 (x∗) ∂f1/∂x2(x∗) ]
[ (∂f2/∂x1(x∗) ∂f2/∂x2(x∗) ]

is the Jacobian evaluated at x∗

System of planar linear ODEs
saying if (y_1,y_2) grows or decays
=> depends on the eigenvalues of J∗

Question 43

Nonlinear two-dimensional ODE systems: Linear stability analysis

remarks

Answer 43

the stability of a hyperbolic equilibrium and dynamics
is given by the linear stability analysis and provided by the Jury’s conditions on and for planar ODEs.

x∗ is non-hyperbolic if J∗ has at least one purely imaginary eigenvalue

The linear stability analysis and Jury conditions don’t say anything about the stability of non-hyperbolic equilibria. Linear stability does not say anything also about the case

where 4det(J∗) = (tr(J∗))²: borderline case between complex and real eigenvalues (diagram)

Question 44

Phase plane analysis

Answer 44

*studying solution curve or trajectory (x₁(t),x₂(t))
*useful for qualitative understanding of global dynamics in x₁-x₂ plane
*Solution with given IC exists and is unique by Cauchy-Lipschitz thm there is only ONE TRAJECTORY passing though a given (x₁(0),x₂(0))
trajectories of sols therefore do not cross each other
* Trajectories (x₁(t),x₂(t)) are parametric equations with t as parameter

Question 45

Phase plane
slope of trajectory

Answer 45

Slope at (x₁,x₂) of trajectory in x₁-x₂ plane given by

x˙₂ /x˙₁ = dx₂/dx₁=
f₂(x₁, x₂)/f₁(x₁, x₂)

ratio gives slope of trajectory in the phase plane

Question 46

Phase plane
tangent vector

Answer 46

The tangent vector
(f₁(x₁, x₂), f₂(x₁, x₂))^T

gives the direction of the trajectory everywhere in the plane, except at the equilibria
(where the flow stops).

Question 47

Phase plane
Nullcline

Answer 47

The x₁-nullcline are the set of all points in the phase plane s.t
f₁(x₁, x₂)=0
The x₂-nullcline are the set of all points in the phase plane s.t
f₂(x₁, x₂)=0

The tangent vector along the x₁ -nullcline is
(0)
(f₂) parallel to x₂-axis

The tangent vector along the x₂ -nullcline is
(0)
(f₁) parallel to x₁-axis

• Solutions along x₁-nullclines are parallel to x₂ axis;
  solutions along x₂-nullclines are parallel to x₁ axis
  =>nullclines intersect at equilibria

Nullclines separate phase space where change sign of x * _i
used to sketch phase portrait

Question 48

Two-species interaction models:
consider quadratic
interactions between species 1 and 2 in an autonomous model of ODEs

Answer 48

x₁(t), x₂(t) represent densities of 2 interacting pops:
requires x₁(t)≥ 0 , x₂(t)≥ 0 for all t≥ 0
system:
f₁(x₁, x₂)= αx₁(t)+ βx₁(t)x₂(t) and
f₂(x₁, x₂) =γx₁(t) +δx₁(t)x₂(t), where |α|, |γ| are the birth/death rates of each species,
|β| and |δ|are the rates of predation (non-linear contributions).

(α, β, γ, δ : parameters (positive or negative), allowing us to model “mutualism”,”competition” and “predator-prey interactions”)

NONLINEAR SYSTEM

model is:
d/dt (x)=
d/dt(
(x₁)
(x₂) =
(αx₁(t)+ βx₁(t)x₂(t) )
(γx₁(t) +δx₁(t)x₂(t))
=
(f₁(x₁, x₂))
(f₂(x₁, x₂))

Question 49

Two species interaction models:
model is:
d/dt (x)=
d/dt(
(x₁)
( x₂) =
(αx₁(t)+ βx₁(t)x₂(t) )
(γx₁(t) +δx₁(t)x₂(t))
=
(f₁(x₁, x₂))
(f₂(x₁, x₂))

NULLCLINES

EQUILIBRIA

Answer 49

x₁-nullcline:
f₁(x₁, x₂)=αx₁(t)+ βx₁(t)x₂(t)=0
x₁=0 and
x₂=-α/β

x₂-nullcline
f₂(x₁, x₂)= γx₁(t) +δx₁(t)x₂(t) =0
x₁= - γ/δ
x₂=0

Nullclines are straight
horizontal/vertical lines

x∗_{a,b} (when physical) given by f_1(x∗a,b) = f_2(x∗a,b) = 0 =

Physical where cross:
x∗_a = (0, 0)^T [extinction, always physical],
x∗_b = (−γ/δ, −α/β)^T [coexistence when γ/δ < 0 and α/β < 0]

Question 50

Two species interaction models:
model is:
d/dt (x)=
d/dt(
(x₁)
( x₂) =
(αx₁(t)+ βx₁(t)x₂(t) )
(γx₁(t) +δx₁(t)x₂(t))
=
(f₁(x₁, x₂))
(f₂(x₁, x₂))

STABILITY

Answer 50

Linear stability => Jacobian
J=
[α + βx₂ βx₁]
[ δx₂ γ + δx₁]
@x∗_a = (0, 0)^T

J∗_a=
[α 0]
[0 γ]
λ− = α and λ+ = γ
tr(J∗_a) = α + γ, det(J∗_a) = αγ

@x∗_b = (−γ/δ, −α/β)^T

J∗_b=
[0 −βγ/δ]
[−αδ/β 0]

λ± = ±√αγ
tr(J∗_b) = 0, det(J∗_b) = −αγ
JURY CONDITIONS FOR ODES

Question 51

Jury’s conditions for ODEs:

Two species interaction models:
model is:
d/dt (x)=
d/dt(
(x₁)
( x₂) =
(αx₁(t)+ βx₁(t)x₂(t) )
(γx₁(t) +δx₁(t)x₂(t))
=
(f₁(x₁, x₂))
(f₂(x₁, x₂))

Answer 51

first:
if α, γ have the same sign, opposite to sign of β, δ=>x∗_b is a saddle (unstable)

if α, γ have opposite sign and =>x∗_b is physical+non-hyperbolic, no conclusion

second:
if α<0, γ<0 x∗_a is asymptotically stable (stable node)
if αγ<0 x∗_a is a saddle
if α>0 γ>0 x∗_a is an unstable node

Question 52

Models of mutualism:
α, β, γ, δ : parameters

equilibria

Answer 52

α<0 γ<0 : linear terms
if no interactions species would die
δ>0 β>0 : nonlinear terms
mutualism; species help each other

equilibria:
x∗_a = (0, 0)^T [extinction, always physical],
STABLE NODE
tr(J∗_a) = α + γ < 0, det(J∗_a) = αγ > 0

x∗_b = (−γ/δ, −α/β)^T
= (|γ/δ|, |α/β|)^T [coexistence when γ/δ < 0 and α/β < 0]
UNSTABLE SADDLE
tr(J∗_b) = 0, det(J∗_b) = −αγ < 0

Question 53

Models of mutualism:
Nullclines

Slope of solution trajectories in x1-x2 plane given by

Answer 53

nullclines:
x₁-nullcline:
x₁=0 and x₂=-α/β = |α/β|

x₂-nullclines
x₁= - γ/δ = |γ/δ| and x₂=0

along the nullclines in this example the trajectories will be perpendicular
x˙ _1 = f_1 = βx_1(x_2 − |α/β|)
x˙ _2 = f_2 = δx_2(x_1 − |γ/δ|)

In the physical quadrant, seperate into 4 regions:
x˙ _1>0 if x_2 >|α/β|

x˙ _2>0 if x_1 > |γ/δ|
DIAGRAM
sketch trajectories
either both extinct or both thrive (trajectories at origin or grow without bound)

Model of mutualism: if populations are large enough they thrive and grow without bound.
Otherwise, both go extinct.

Question 54

Models of Competition
α, β, γ, δ : parameters

equilibria

Answer 54

Species 1 can grow while species 2 can only die, and competition is detrimental for both species.

α>0 γ<0 : linear terms
birth rate of 1, death rate of 2
δ<0 β<0 : nonlinear terms
interaction detriment for both species
(|β| and |δ| are the rate at individuals of species 1 and 2 die due to harsh competition for resources)
|δ| rate at which species 1 kills species 2

equilibria:
only physical equilibrium
x∗_a = (0, 0)^T [extinction, always physical],

det(J∗_a) = α γ < 0,

x∗_a is a saddle (unstable). Species 2 goes extinct, while species 1 grows without bound

Question 55

Predator prey interactions model
α, β, γ, δ : parameters

equilibria

When x_1 → 0 :
When x_2 → 0 :

Answer 55

α<0 |α| death rate of species 1
left alone will only die
β>0 rate of predation of species 1 on species 2

γ>0 birth rate of species 2
if left alone reproduces

|δ|rate at which species 1 kills species 2
species 1: predators
species 2: prey

2 physical equilibria=>
x∗_a = (0, 0)^T EXTINCTION (UNSTABLE SADDLE)
x∗_b=(γ/|δ|,|α|/β)^T COEXISTENCE (non-hyperbolic)

When x_1 → 0 :x_2 ≈ γx_2 so x_2~ exp(γt) converges to infinity as t tends to infinity PREY EXPLODE WIHOUT PRED
When x_2 → 0 :x_1 ≈ -|α|x_1 so x_1~ exp(-|α|t) converges to 0 as t tends to infinity PRED WIPED OUT WITHOUT PREY

Less predators => more prey => more predators => less prey => less predators => …
Coexistence equilibrium is non-hyperbolic- cant use linearity analysis. Can predators and prey coexist?

Question 56

Lotka-Volterra model
is a predator-prey model

Answer 56

x~#prey in a reference volume V
y~ #predators in V

Assumptions:
prey are limited only by predators;
predator’s functional response is linear (predation term ∝ y);
no competition between predators in finding prey, predation term ∝ xy; when x → 0 (no prey left), the population of predator vanishes exponentially quickly (y → 0);
every prey’s death contributes identically to the predators’ growth.

Question 57

Lotka-Volterra predator prey model

FORM

FISHING EFFORT

Answer 57

dx/dτ = [a − pϵ]x(τ ) − gx(τ )y(τ ),
dy/dτ = fx(τ )y(τ ) − [b + qϵ]y(τ ),

a (prey birth rate),
b (predator death rate),
g and f (predation rates)
all positive parameters, and τ denotes the time

x decreases if interaction between
y increases if interactions between
FISHING EFFORT 0<pϵ<a>0
enhanced/contribution to death rate caused by fishing</a>

Question 58

dx/dτ = [a − pϵ]x(τ ) − gx(τ )y(τ ),
dy/dτ = fx(τ )y(τ ) − [b + qϵ]y(τ ),

Lotka-Volterra predator prey model
Equilibria

Ratio of predator/prey in coexistence equilibrium:

Answer 58

Physical SS
found by solving simultaneously dx/dτ =dy/dτ =0

x∗a ≡ (x∗a, y∗a)^T = (0, 0)^T => extinction (saddle, unstable)
and
coexistence equilibrium: x∗_b = (x∗_b , y∗_b )^T =
((qϵ + b)/f , (a − pϵ)/g)^T

Question 59

Volterra principle

Answer 59

Ratio of pred/prey in coexistence eq:
y∗_b/x∗_b = [f (a − pϵ)]/g[(qϵ + b)]

This ratio increases as fishing effort decreases

Volterra principle
decrease epsilon, decrease fishing effort means increase predator proportion
(less fishing=> more predators)

Question 60

Lotka-Volterra predator prey model

Does the LVM settle in the coexistence equilibrium ?
Does the Volterra principle hold?

Answer 60

Jacobian at coexistence eq. is
J∗ᵦ =
[0 −gxᵦ]
[fy∗ᵦ 0]

eigenvalues
λ± = ± i sqrt(fgx∗_b y∗_b)
purely imaginary
so x∗_b is non-hyperbolic
=> no conclusion from linear stability
so use phase plane portait

Question 61

Lotka-Volterra predator prey model
PHASE PORTRAIT

Answer 61

x−nullclines are x = 0 and y = (a − pϵ)/g
y−nullclines are y = 0 and x = (b + qϵ)/f

x >˙ 0 where 0 < y < (a − pϵ)/g and
x <˙ 0 where y > (a − pϵ)/g
y >˙ 0 where x > (b + qϵ)/f and
y <˙ 0 where 0 < x < (b + qϵ)/f

suggests orbits surrounding x∗_b

Question 62

Lotka-Volterra predator prey model
Further analysis: constant of motion

STEP 1 NON DIMENSIONALISATION
introduce u and v
x → u = x/x∗_b ,
y → v = y/y∗_b ,
τ → t = (a − pϵ)τ

Answer 62

Further analysis: constant of motion => closed orbits and oscillatory behaviour
STEP 1 NON DIMENSIONALISATION
introduce u and v
x → u = x/x∗_b ,
y → v = y/y∗_b ,
τ → t = (a − pϵ)τ

STEP 2
du/dt = (1/x∗_b) (dx/dτ) (dτ/dt) = u(1 − v)

dv/dt = (1/y∗_b) ( dy/dτ) (dτ/dt)
= (( qϵ + b)/(a − pϵ)) v(u − 1)

so
u* = du/dt = u(t)[1 − v(t)]
v* = dv/dt = κv(t)[u(t) − 1

with κ= (qϵ + b)/(a − pϵ)

The nonhyperbolic coexistence equilibrium point is now
(u, v) = (u∗, v∗) = (1, 1)

CONSTANT OF MOTION
Φ(u(t), v(t)) = (u − ln u) + κ⁻¹(v − ln v)

Question 63

CONSTANT OF MOTION
Φ(u(t), v(t)) = (u − ln u) + κ⁻¹(v − ln v)
Lotka-Volterra predator prey model

Answer 63

If you look at time derivative of this quantity it is 0:
means this is constant
means that the eq for u and v preserves this quanitity Φ
Φ does not change:
Φ=dΦ/dt
= u - (u)/u + (1/κ)(v - (v/v))
=(u/u_(u-1) + (v/κv)(v-1) using u/u =1-v and v*/κv = u-1
=(1-v)(u-1)+(u-1)(v-1)=0

since
Φ dot =0
Φ(t)=Φ(0) is conserved
so we have closed orbits around (u∗, v∗) (x∗_b, y∗_b)
closed orbits in the phase plane arounf point
oscillations around point in the LVM time series

orbit forever: Trajectories not started at IC SS
Solutions of the LVM are closed orbits surrounding the coexistence equilibrium and set by the conservation of for all regular oscillations around coexistence equilibrium that is marginally stable-stays at that particular orbit, period and oscillation

Question 64

We now know that is not asymptotically stable => Does the Volterra principle actually hold?

Answer 64

In Example Sheet 2, we show that the time average of (u, v) over one period T
orbit is
u¯ = integral over [0,T] of u(t).dt/T=1
v¯ = integral over [0,T] of v(t).dt/T=1
(¯u, v¯) = (u∗, v∗) = (1, 1), and (¯x, y¯) = (x∗_b , y∗_b )

The Volterra principle therefore holds for the time-averaged variables

Question 65

LVM
closed orbits are set by the initial value of
ϕ(u(0), v(0))

Answer 65

Two slightly different initial conditions lead to different oscillatory behaviour.

LVM closed orbits are therefore “marginally stable” or “neutrally stable”: they are not robust to any small perturbations => a small initial change modifies the period/frequency and amplitude entirely

Oscillations in the Lotka-Volterra model are not robust (only marginally stable) as they result from the constant of motion ϕ (changing IC changes dynamics)

Various attempts to generalise the model and obtain robust (stable) oscillatory dynamics, e.g. Holling response => stable limit cycle and sustained oscillations

Question 66

Extra MATH5567M topic: Hopf bifurcation in predator-prey models & Holling II response

Answer 66

dynamics can change when an equilibrium point becomes non hyperbolic

Pred-prey models characterised by stable limit cyclces can predict robust and persistent oscillations

the kind of oscillation may dep on IC or not

Question 67

Supercritical Hopf bifurcation

Answer 67

In planar 2D ODEs
HB occurs when an equilibrium loses its stability as a pair of complex conjugate eigenvalues of J∗crosses the imaginary axis of the complex plane.

HB accompanied by the emergence of a limit cycle yielding a regular oscillatory behaviour around x∗

model that admits a SHB can be recast

e.g in the example spiral towards origin, µ<0 ; HB µ=0;µ>0 spiral from origin to outside stable limit cycle or from out to stable limit cycle

Application to predator-prey model: LVM with growth limiting term on prey and “Holling II
response” => Hopf bifurcation, limit cycle, stable oscillations.

Question 68

Supercritical Hopf bifurcation
normal form

Answer 68

Normal form of supercritical HB:
dx₁/dt = x˙₁ = −x₂ + x₁(µ − x₁² − x₂²);
dx˙/dt=x˙ ₂ = x₁ + x₂(µ − x₁² − x₂²);

EQP obtained by setting
x˙₁ = x˙₂ = 0
⇒ (x₁ , x₂ ) = (x∗1, x∗2) = (0, 0)
only one
stability from:
Jacobian J(x₁ , x₂ )=
[µ − 3x₁² − x₂² −1 − 2x₁x₂]
[1 − 2x₁x₂ µ − x₁² − 3x₂²]

J∗ = J(0, 0) =
[µ −1]
[1 µ]

eigenvalues λ± = µ ± i

(x∗_1, x∗_2) = (0, 0) is
asymptotically stable when µ < 0 (SPIRAL INWARDS)
and unstable when µ > 0 (SPIRAL OUTWARDS)

Useful to use polar coordinates defined by
A supercritical Hopf bifurcation arises when µ = 0

Question 69

A supercritical Hopf bifurcation arises when µ = 0

in predator-prey models

working out

Answer 69

CHANGE OF BEHAVIOUR/ BIFURCATION at µ = 0

using polar coords
x_1 = r cos θ, x_2 = r sin θ
r^2= x_1^2+x_2^2
differentiating gives
2rr= 2(x₁x₁ + x₂x₂*)
divide both sides by 2r

r= (x₁/r)x₁ + (x₂/r)x*₂
= cosθ[ -rsin θ + (r/(µ-r²)cos θ] + sinθ[ rcosθ + (r/(µ-r²)sinθ] +
=r(µ-r²)(cos²θ+sin²θ)= r(µ-r²)

d/dt ( cosθ)= - θ* sinθ = (d/dt)(x₁/r)= [x₁r- r x₁]/r²
gives
θ= -x₁/rsinθ)+ rx₁/(r²sinθ) = (rx₁)/(rx₂)- (x₁*/x₂)
=

[(µ-r²)rcosθ+rsinθ -(µ-r²)rcosθ]/rsinθ = 1
gives…

Question 70

A supercritical Hopf bifurcation arises when µ = 0

in predator-prey models
NORMAL FORM IN POLAR COORDINATES
RADIAL EQ

Answer 70

r˙ = r(µ − r²) radial dynamics
˙θ = 1 angular dynamics
thus
∂r˙/∂r = µ − 3r²
RADIAL EQ

Question 71

A supercritical Hopf bifurcation arises when µ = 0

in predator-prey models

radial eq equilibria

r˙ = r(µ − r²) radial dynamics
˙θ = 1 angular dynamics
thus
∂r˙/∂r = µ − 3r²
RADIAL EQ

Answer 71

This is a scalar ODE
Radial equation has 2 equilibria r˙ = 0:
r = 0 (i.e. x1 = x2 = 0) and
r = r∗ = √ µ

STABILITY:
∂r˙/∂r|ᵣ₌₀= µ
∂r˙/∂r|ᵣ₌√µ = −2µ

µ < 0,
r = 0 is asymptotically stable,
r∗ = √µ is unstable

µ > 0,
r∗ = √µ is asymptotically stable
r = 0 is unstable

**a stable limit cycle of radius **
r∗ = √µ emerges when µ > 0
that is a stable periodic orbit of
radius r∗ = √µ and
period 2π arises when µ > 0

Question 72

cycles on LVM vs PRED PREY

Answer 72

*different nature to LVM; here r∗ asymptotically stable no matter IC you end up on the cycle whereas in the LVM you dep on IC as constant of motion will set the orbit

*pred-prey lead to limit cycle; robust oscillations around equilibrium

*LVM imaginary stable oscillation; changing IC changes orbits

Question 73

3.7.2 Predator-prey model with Holling (type II) response

Criticisms to LVM
reasoning

Answer 73

prey growing exponentially in absence of predators

structurally unstable nature of the closed orbits

we address these by introducing a growth-limiting term (logistic type) in the ODE for the prey

we assume a Holling type 2 response obtained by replacing the LVM predation terms by terms in the form xy/(1+x) which will give rise to robust oscilattions associated with a stable limit cycle

Question 74

3.7.2 Predator-prey model with Holling (type II) response

model

Answer 74

x: number of prey, y: number of predators =>
x ≥ 0 and y ≥ 0

x˙ = x (1 − x/γ) − xy/(a+x) = f₁(x, y)
(x*growth-lim. term - predation/Holling II funct)

y˙ = xy/2(1 + x)−y/4 = f₂(x, y)
+ (predation) Holling II num- linear term death of pred

γ > 1 is the bifurcation parameter and will affect behaviour

Question 75

3.7.2 Predator-prey model with Holling (type II) response

equilibria
stability

x˙ = x (1 − x/γ) − xy/(a+x) = f₁(x, y)
(x*growth-lim. term - predation/Holling II funct)

y˙ = xy/2(1 + x)−y/4 = f₂(x, y)
+ (predation) Holling II num- linear term death of pred

Answer 75

f₁(x, y) =f₂(x, y)=0
physical SS gamma >1

x∗_0 = (0, 0)T (complete extinction)
x∗_1 = (γ, 0)T (predators die out
x∗_2 = (1, 2(1 − γ−1))T
Jacobian at (x,y) is
[1 − 2x/γ − (y)/((1+x)²) − x/(1+x)]
[(y/2) (1 + x)⁻² x/(2(1+x)) − (1/4]

x∗_0 is always a saddle (unstable)
x∗_1 is always a saddle (unstable)

J(x∗_2) = (1/2)
( 1 − (3/γ) −1)
( (γ−1)/2γ 0)

tr(J∗_2) = (1/2)(1 − 3/γ)< 0 if γ < 3 and > 0 if γ > 3,
det(J∗2) = (γ − 1)/8γ> 0

x∗_2 is asymptotically stable if 1< γ <3
(stable spiral if √13 − 2 < γ < 3, o/w node)
and is unstable when γ >3

*γ >3 all equilibria unstable

Question 76

when γ >3 all equilibria unstable

x˙ = x (1 − x/γ) − xy/(a+x) = f₁(x, y)
(x*growth-lim. term - predation/Holling II funct)

y˙ = xy/2(1 + x)−y/4 = f₂(x, y)
+ (predation) Holling II num- linear term death of pred

Answer 76

all equilibria unstable
species will not go extinct
pred not die out
and not coexist at the given density

Question 77

When γ = 3

x˙ = x (1 − x/γ) − xy/(a+x) = f₁(x, y)
(x*growth-lim. term - predation/Holling II funct)

y˙ = xy/2(1 + x)−y/4 = f₂(x, y)
+ (predation) Holling II num- linear term death of pred

Answer 77

tr(J∗_2) = 0 ⇒
purely imaginary eigenvalues (conjugate pair): Hopf bifurcation

thus when γ > 3
there is a stable limit cycle and regular/robust oscillations around x∗_2

Linear stability does not help to see whether the limit cycle is stable.
Phase plane analysis + Poincar´e-Bendixson theorem

Question 78

PHASE PLANE

x˙ = x (1 − x/γ) − xy/(a+x) = f₁(x, y)
(x*growth-lim. term - predation/Holling II funct)

y˙ = xy/2(1 + x)−y/4 = f₂(x, y)
+ (predation) Holling II num- linear term death of pred

Answer 78

x-nullclines: f_1=0
x=0
y = (1 − x/γ)(1 + x) (parabola)

y-nullcline:f_2=0
y=0
x=1

x˙ <0 above parabola
x˙ >0 below parabola
y˙ >0 if x>1
y˙ <0 if x<1

stable limit cycle generated (diagram)

Question 79

Poincar´e-Bendixson theorem in example

x˙ = x (1 − x/γ) − xy/(a+x) = f₁(x, y)
(x*growth-lim. term - predation/Holling II funct)

y˙ = xy/2(1 + x)−y/4 = f₂(x, y)
+ (predation) Holling II num- linear term death of pred

Answer 79

we consider a closed bounded region (within dotted
rectangle) punctured at x∗_2 (repeller when γ > 3).

By thm since the trapping region contains neither saddle nor attractor, it contains a closed orbit
=> limit cycle is stable
=> regular and robust oscillations

Here, different initial conditions x(0), y(0) lead to the same
regular oscillations (after a short transient), with
frequency/period independent of the initial condition
=> ROBUST OSCILLATORY BEHAVIOUR
(IC transient behaviour initially different but then settles into same cycle)
That’s a striking difference with the oscillations in the LVM
which dep on the IC