Scientific Computing I Part I Module 3: Population Modelling - - PDF document

Scientific Computing I

Module 3: Population Modelling – Continuous Models Michael Bader

Lehrstuhl Informatik V

Winter 2007/2008

Part I ODE Models Discrete vs. Contiuous Models

dp dt = F(p,t,...) p(t) =? discrete model: continuous model: p(t) ∈ N individuals p : R → R,p(t) =? Move to Continuous Models: easier(?) type of mathematical problem: differential equations, calculus analytical solutions available(?)

Model of Malthus (1798)

Only one species:

1

birth rate γ (number of births per time interval) proportional to size of population

2

death rate δ proportional to size of population

3

thus: constant growth (or decay) rate: r = γ −δ Modelling: constant growth rate dp dt = r ·p growth within a time interval p(t+δt) = p(t)+r ·p(t)·δt

Model of Malthus – Differential Equation

Written as an ordinary differential equation: ˙ p(t) = r ·p(t) Requires initial condition (population at start) p(0) = p0 Analytical solution: p(t) = p0 ert

Model of Malthus – Solutions

The model of Malthus describes exponential growth or decay of a population:

1,5 1 0,5 8 6 4 2 p 3 2,5 2

Model of Verhulst (19th century)

Objective:

model populations that approach saturation value

Assumptions:

growth/death rate depend on population size; assume linear dependency: g(t) = g0 −g1 ·p(t) d(t) = d0 +d1 ·p(t) leads to differential equation: ˙ p(t) = g(t)−d(t) = (g0 −d0)

• =:α

−(g1 +d1)

• =:β

p(t)

Model of Verhulst – Saturation

solve initial value problem: ˙ p(t) = α −βp(t), p(0) = p0 solution: p(t) = p∞ +e−βt (p0 −p∞), p∞ = α β

1,5 1 0,5 8 6 4 2 p 3 2,5 2

Model of Verhulst – Logistic Growth

saturation model does no longer model exponential growth idea: let growth/death rate decrease linearly with size of population but keep growth/death rate proportional to population size leads to differential equation: ˙ p(t) = (α −βp(t))p(t)

Logistic Growth

• ther formulation

˙ p(t) = α

• 1− p(t)

β

• p(t)

solution: p(t) = β (1−e−αt)+ β

p0e−αt

1,5 1 0,5 8 6 4 2 p 3 2,5 2

Logistic Growth with Threshold

extended version of Verhulst’s model: ˙ p(t) = α

• 1− p(t)

β

• 1− p(t)

δ

• p(t)

solutions (β = 2,δ = 4):

p 5 4 3 2 1 8 6 4 2

Example – The Passenger Pigeon

beginning of the 19th century, estimated population in North America: four billion hunting diminished its number below a critical threshold (late 1880s) The last passenger pigeon died on September, 1st 1914.

Scientific Computing I

Module 3: Population Modelling – Continuous Models (Parts II and III) Michael Bader

Lehrstuhl Informatik V

Winter 2007/2008

Part II More Than One Species – Systems of ODE A Linear Model

similar to Verhulst’s saturation model additional growth term proportional to other species leads to system of differential equations: ˙ p(t) = b1 +a11p(t)+a12q(t) ˙ q(t) = b2 +a21p(t)+a22q(t) typically:

b1 > 0,b2 > 0 (growth term) a11 < 0,a22 < 0 (saturation) a12,a21?

First Example: Arms Race

armament of two (hostile) countries

• ur suspicion: a12 > 0, a21 > 0

Observation:

long-time behaviour depends on size of parameters steady-state solutions exist solutions exist that show unlimited growth

Second Example: Competition

two species sharing a common natural habitat competition: a12 < 0, a21 < 0

Observation:

long-time behaviour depends on size of parameters steady-state solutions exist some scenarios are physically incorrect! (negative population size)

A Non-Linear Model

similar to Verhulst’s logistic growth model additional growth term proportional to other species leads to system of differential equations: ˙ p(t) = (b1 +a11p(t)+a12q(t))p(t) ˙ q(t) = (b2 +a21p(t)+a22q(t))q(t) typically:

b1 > 0,b2 > 0 (growth term) a11 < 0,a22 < 0 (saturation) a12,a21?

The Non-Linear Competition Model

two species sharing a common natural habitat competition: a12 < 0, a21 < 0

Possible Scenarios:

steady-state

• ne species dies out (extinction)

no obvious nonsense

Competition – Steady State

system of differential equations: ˙ p(t) =

• 5

2 + √ 3 24 − 5 8p(t)− √ 3 24 q(t)

• p(t)

˙ q(t) =

• 7

8 + 3 √ 3 2 − 3 √ 3 8 p(t)− 7 8q(t)

• q(t)

solution for p0 = 1

4, q0 = 3:

2 1 t 10 8 6 4 p 2 4 3

Competition – Extinction

system of differential equations: ˙ p(t) =

• 71

8 − 23 12p(t)− 25 12q(t)

• p(t)

˙ q(t) =

• 73

8 − 25 12p(t)− 23 12q(t)

• q(t)

solution for p0 = 1

4, q0 = 1 4:

t 10 p 8 4 6 3 2 4 1 2

Predator-Prey

two species: predator p and prey q predator eats prey: a12 > 0 prey is eaten by predator: a21 < 0

Possible Scenarios:

stable oscillations

• ne species dies out (what happens with the other,

then?)

Predator-Prey by Lotka & Volterra

system of differential equations: ˙ p(t) =

• −1

2 + 1 200q(t)

• p(t)

˙ q(t) = 1

5 − 1 50p(t)

• q(t)

solution for p0 = 6, q0 = 50:

t 100 80 p 60 160 40 120 80 20 40

Open Questions

Methods to Analyse a Given Model?

predict approximate solution or shape of solution? predict possible steady states? predict critical points? (species on edge of extiction?)

Methods to Improve Modeling?

predict failure of the model? tune parameters to model a specific situation?

Part III Discussion and Analysis of ODE Models Analysing the Slope of a Solution

Example: Model of Malthus

˙ p(t) = αp(t) for a sensible solution: p(t) > 0 α decides slope of solution:

α > 0: growing population (accelerated growth) α < 0: receding population (decelerated reduction)

Points of Equilibrium

Example: Model of Verhulst (saturation)

˙ p(t) = α −βp(t) equilibrium: ˙ p(t) = 0

• nly, if p(t) = α

β

Example: Logistic Growth

˙ p(t) = α

• 1− p(t)

β

• p(t)

constant solution, if p(t) = β or p(t) = 0

Critical Points

Observation on Logistic Growth:

constant solution p(t) = β, if p(0) = β constant solution p(t) = 0, if p(0) = 0 equilibrium at p = β is reached for nearly all initial conditions ⇒ attractive (stable) equilibrium equilibrium at p = 0 is not reached for any other initial conditions (“repulsive”) ⇒ unstable equilibrium

Critical Points – Derivatives

Examine derivatives:

critical point p = ¯ p attractive equilibrium (asymptotically stable): ˙ p < 0 for p = ¯ p+ε ˙ p > 0 for p = ¯ p−ε unstable equilibrium: ˙ p > 0 for p = ¯ p+ε ˙ p < 0 for p = ¯ p−ε

• therwise: saddle point

Direction Field

plot derivatives vs. time and size of population:

Example: Logistic Growth

˙ p(t) = α

• 1− p(t)

β

• p(t)

p(t) 3 2,5 2 1,5 1 t 0,5 10 8 6 4 2

Direction Field (2)

Example: Logistic Growth with Threshold

˙ p(t) = α

• 1− p(t)

β

• 1− p(t)

δ

• p(t)

p(t) t 5 10 4 3 8 2 1 6 4 2

Identifying Critical Points

attractive equilibrium: unstable equilibrium saddle point

Critical Points in 2D

Example: Arms Race

system of differential equations equilibrium: ˙ p = 0, ˙ q = 0 ˙ p(t) = b1 +a11p(t)+a12q(t) = 0 ˙ q(t) = b2 +a21p(t)+a22q(t) = 0 solution of a linear system of equations: a11p(t)+a12q(t) = −b1 a21p(t)+a22q(t) = −b2 in most cases one critical point critical line, if system matrix is singular

Direction Field for a System of ODE

example: 2D system of differential equations: ˙ p(t) = b1 +a11p(t)+a12q(t) ˙ q(t) = b2 +a21p(t)+a22q(t) natural exension: 3D plot: t vs. p vs. q 1D direction field for p vs. t or q vs. t not sufficient: what values to chose for q (or p resp.)? but: stationary problem ⇒ independent of t thus: plot directions depending on p and q

2D Direction Field – Arms Race

system of differential equations: ˙ p(t) =

3 2 −p(t)+ 1 2q(t)

˙ q(t) = 0+ 1

2p(t)−q(t)

direction field – with critical point at (2,1):

1 4 q 3 p 0,5 2 1,5 2 1

Arms Race – unlimited growth

system of differential equations: ˙ p(t) =

1 2 − 3 4p(t)+q(t)

˙ q(t) = −5

4 +p(t)− 3 4q(t)

direction field – with critical point at (2,1):

q 4 p 2 3 5 5 2 1 3 1 4

Arms race – the peaceful neighbour

system of differential equations: ˙ p(t) = 0− 3

4p(t)+q(t)

˙ q(t) =

5 2 −p(t)− 3 4q(t)

direction field – with critical point at

• 8

5, 6 5

• :

1 0,5 2,5 p 0,5 q 3 1 1,5 3 2,5 2 1,5 2

Nonlinear System – Competition

system of differential equations: ˙ p(t) =

• 5

2 + √ 3 24 − 5 8p(t)− √ 3 24 q(t)

• p(t)

˙ q(t) =

• 7

8 + 3 √ 3 2 − 3 √ 3 8 p(t)− 7 8q(t)

• q(t)

direction field – critical points at (4,1),...:

4 4 5 2 3 1 2 1 3 p q

Nonlinear System – Extinction

system of differential equations: ˙ p(t) =

• 71

8 − 23 12p(t)− 25 12q(t)

• p(t)

˙ q(t) =

• 73

8 − 25 12p(t)− 23 12q(t)

• q(t)

critical points at (0,4.76...),(4.63...,0),...:

0,5 3 q 1,5 p 3 4 1 2 2 1 5 2,5

Lotka & Volterra

system of differential equations: ˙ p(t) =

• −1

2 + 1 200q(t)

• p(t)

˙ q(t) = 1

5 − 1 50p(t)

• q(t)

direction field – with critical point at (10,100):

q 40 30 150 200 100 50 20 10 p

2D Critical Points – Summary

Different types of critical points in 2D: attractive/stable equilibrium (arms race – steady state) unstable equilibrium saddle point (arms race – unlimited growth) attractive “spiral point” (“peaceful neighbour”) unstable “spiral point” centre of “rotation” (Lotka-Volterra) ⇒ How to discriminate between these types?

Homogeneous Systems of ODE

Homogeneous System in matrix-vector-notation: ˙ x = Ax x : R → Rn, A ∈ Rn×n example: x(t) = (p(t),q(t)) Solutions: let xλ be an eigenvector: Axλ = λxλ then xλeλt is a solution: Axλeλt = λxλeλt = d dt

• xλeλt

q.e.d.

Eigenvectors and Eigenvalues

Corollaries: the solutions of the homogeneous system ˙ x = Ax are linear combinations of the respective eigen-solutions: xhom(t) = ∑

λ

aλxλeλt, aλ ∈ R the solutions of the inhomogeneous system ˙ x = Ax+b are x(t) = −A−1b+xhom(t)

• bservation: xc = −A−1b is a critical point!

Eigenvalues and Critical Points

the ODE system ˙ x = Ax+b is solved by x(t) = xc +∑

λ

aλxλeλt xc attractive equilibrium, lim

t→∞x(t) = xc,

• nly if eλt → 0 for all eigenvalues λ

λ ∈ R ⇒ λ < 0 λ = µ +iν ⇒ µ < 0 (eiνt = cosνt+isinνt)

Stability of Linear Systems

Overview:

• eigenval. (λj = µj +iνj)

critical point stability real, all λ < 0 node stable, attr. real, all λ > 0 node unstable real, λk > 0,λl < 0 saddle point unstable complex, all µ < 0 spiral point stable, attr. complex, all µ > 0 spiral point unstable complex, all µ = 0 centre stable

Stability of 2D Systems

Real Eigenvalues:

λ1 < 0, λ2 < 0, attractive equilibrium

x2 x1 eig2 eig1

Stability of 2D Systems

Real Eigenvalues:

λ1 > 0, λ2 > 0, unstable equilibrium

x2 x1 eig2 eig1

Stability of 2D Systems

Real Eigenvalues:

λ1 > 0, λ2 < 0, saddle point

x2 x1 eig2 eig1

Stability of 2D Systems

Complex Eigenvalues:

µ1 < 0, µ2 < 0, spiral point (asympt. stable)

x2 x1 eig2 eig1

Stability of 2D Systems

Complex Eigenvalues:

µ1 > 0, µ2 > 0, spiral point (unstable)

x2 x1 eig2 eig1

Stability of 2D Systems

Complex Eigenvalues:

µ1 = µ2 = 0, centre of oscillation

x2 x1 eig2 eig1

Stability of Non-Linear Systems

2D system of ODE: ˙ x(t) = f(x(t)), f : Rn → Rn nonlinear critical point at xc: f(xc) = 0 for analysis of critical points: linearization ˙ x(t) = f(x(t)) ≈ f(xc)

=0

+Jf(xc)(x(t)−xc) examine eigenvalues of Jf(xc)