Predator Prey System
with a stable periodic orbit
1st Session - Simple Analysis
Systems Dynamics Study Group
Ellis S. Nolley
11/7/2001
12/20/2001 Systems Dynamics Study Group 1
Topics
• Overview
– Simple Analysis – 1st session, 11/7/2001
– Rigorous Analysis – 2nd session, 11/27/2001
– Simulation Results – 3rd session, 12/11/2001
• Mathematical Model
• Fixed Points
• Stable Periodic Orbit
Reference: McGehee & Armstrong, Journal of Differential Equations,
23, 30-52 (1977)
12/20/2001 Systems Dynamics Study Group 2
Model
x = amount of prey, y = amount of predator
g(x)
dx/dt = xg(x) – yp(x)
dy/dt = y[-s + cp(x)] g(x)
x
k
g(x) is a growth function,
g(x), monotonic non-increasing, dg(x)/dx 0
p(x) is predation function
p(x), monotonic increasing, dp(x)/dx >0 , p(0)=0
12/20/2001 Systems Dynamics Study Group 3
dx/dt = xg(x) – yp(x)
dy/dt = y[-s + cp(x)]
Fixed Points
3 Fixed points: (x*,y*), (0,0), (k,0)
(x*,y*)
dy/dt = 0, dx/dt = 0 for (x*,y*)
At dy/dt=0, y>0, then p(x*) = s/c, y*=x*g(x*)/p(x*)
Assume Lim p(x) = a, as x-> inf+
1) x* > s/c, otherwise there is no fixed point
2) y* > 0, in order to have a system
3) If there is a k, g(k)=0, then x* 0, y*>0
12/20/2001 Systems Dynamics Study Group 4
dx/dt = xg(x) – yp(x)
Fixed Points (Cont’d) dy/dt = y[-s + cp(x)]
Let’s look at the slope on x=k
y
dy/dx = (dy/dt)/(dx/dt)
At x=k, g(k)=0
Recall: -s+cp(x*)=0, p’(x)>0, x*0
-p(k) denominator 0, then delta y>0 (numerator)
So, the vectors are coming in.
12/20/2001 Systems Dynamics Study Group 5
dx/dt = xg(x) – yp(x)
Analysis at Fixed Points dy/dt = y[-s + cp(x)]
(0,0)
What happens at x=0 (y axis)?
dy/dt= y(-s) 0
So, (0,0) is a saddle point.
x
(0,0)
12/20/2001 Systems Dynamics Study Group 6
Analysis at dx/dt = xg(x) – yp(x)
dy/dt = y[-s + cp(x)]
Carrying Capacity
(k,0)
At (k,0), g(k) = 0
dx/dt = xg(x) – yp(x)
= xg(x)
g(x) is monotonic non-increasing.
For x0 y
For x>k, g(x) 0
So dy/dt = y[-s + cp(x)] > 0 for y>0, x=k
x
(k,0) is a saddle point k
12/20/2001 Systems Dynamics Study Group 7
dx/dt = xg(x) – yp(x)
dy/dt = y[-s + cp(x)]
Prey Isocline
y
At the prey isocline, dx/dt = 0 y(0)=g(0)/p’(0)
y= xg(x)/p(x)
and goes through (k,0) and (x*,y*).
To find y(0): by L’Hospital’s Rule, (x*,y*)
y(0) = [x g’(0) + g(0)]/p’(0) = g(0)/p’(0) > 0 x
k
Recall: L’Hospitals Rule: if f(x) & g(x) both go to either 0 or infinity as x->a,
Then lim f(x)/g(x)] = lim [df(x)/dx]/[dg(x)/dx], as x-> a
12/20/2001 Systems Dynamics Study Group 8
The Vector Space dx/dt = xg(x) – yp(x)
dy/dt = y[-s + cp(x)]
Since delta x0 near x=k,
then vector is up near x=k
y
Vectors can only turn around at the critical pt.
At x=x* above y*, dx/dt0
Left of x*, dy/dt 0
x
(x*,y*) is unstable if tangent is positive. k
Pick a line tangent to dy/dx at (k,x**)
All vectors cross it inward
12/20/2001 Systems Dynamics Study Group 9
dx/dt = xg(x) – yp(x)
Periodic Orbit dy/dt = y[-s + cp(x)]
y
• Fixed points are unstable.
• All vectors enter the region
and move away from the boundary.
(k,x**)
• Stable periodic orbit exists around
(x*,y*)
the unstable fixed point.
x
k
12/20/2001 Systems Dynamics Study Group 10
dx/dt = xg(x) – yp(x)
dy/dt = y[-s + cp(x)]
Next Session
Simple Mathematics – 1st session
Rigorous Mathematics – 2nd session, 11/27
Simulation Results – 3rd session
12/20/2001 Systems Dynamics Study Group 11
Thank you!
12/20/2001 Systems Dynamics Study Group 12
Predator Prey System
with a stable periodic orbit
2nd Session - Rigorous Analysis
Systems Dynamics Study Group
Ellis S. Nolley
11/27/2001
12/20/2001 Systems Dynamics Study Group 13
Topics
• Overview
– Simple Analysis – 1st session, 11/7/2001
– Rigorous Analysis – 2nd session, 11/27/2001
– Simulation Results – 3rd session, 12/11/2001
• Mathematical Model
• Fixed Points & Eigenvalues
• Poincare-Bendixon Theorem
4 key slides: #23 – 26
12/20/2001 Systems Dynamics Study Group 14
References
• McGehee & Armstrong, Journal of Differential Equations, 23, 30-52 (1977)
• Morris Hirsch & Stephen Smale, Differential Equations, Dynamical Systems
and Linear Algebra, 1974, Academic Press
Ch 3-5, Linear Systems, Eigenvalues & Exponentials of Operators
Ch 9-12, Stability, Differential Equations on Electrical Systems, Poincare-Bendixon
Theorem, Ecology
• Michael Spivak, Calculus on Manifolds, 1965, W.A Benjamin
• Raghavan Narasimhan, Analysis on Real & Complex Manifolds, 1968, North-
Holland Publishing Company
12/20/2001 Systems Dynamics Study Group 15
Where to find
these References
Mathematics Library, Vincent Hall,
Vincent Hall, 206 Church Street,
3rd Floor, Mpls, MN 55455
University of MN
http://onestop.umn.edu/Maps/VinH/VinH-map.html
12/20/2001 Systems Dynamics Study Group 16
Model
x = amount of prey, y = amount of predator
g(x)
dx/dt = xg(x) – yp(x)
dy/dt = y[-s + cp(x)]
g(x)
x
g(x) is a growth function, k
g(x), monotonic non-increasing, dg(x)/dx 0
p(x) is predation function
p(x), monotonic increasing, dp(x)/dx >0 , p(0)=0
12/20/2001 Systems Dynamics Study Group 17
dx/dt = xg(x) – yp(x)
Jacobian & Eigenvalue Review dy/dt = y[-s + cp(x)]
dx/dt = xg(x) – yp(x)
dy/dt = y[-s + cp(x)]
z’(t) = f(z) = [f1(z1,…zn), …, fn(z1…,zn)]
Note above: z1=x, z2=y, f1(z)=xg(x)-yp(x), f2(z)=y[-s+cp(x)]
dF(z,t)/dt = f(z); F(n)(z)=dnF(z)/dzn,n=0, … ∞; F(0)(z)=F(z)
If z є B(z0,ε) ={z|z-z0| 0
12/20/2001 Systems Dynamics Study Group 20
Why Re λ dx/dt = xg(x) – yp(x)
dy/dt = y[-s + cp(x)]
determines stability dz/dt = f(z)
z’ = f(z); f(z0)=0, z0 fixed point, λk eigenvalues.
Suppose λj has Re λj >0. Pick z close to z0
f(z) = k=0∞Σ f(k)(z0)(z-z0)k/k! = f(z0) + f(1)(z0)(z-z0) + … Taylor Series
~ f(1)(z0) (z-z0) = (z-z0)Σckλk ; d f(z)/z ~ Σckλk dt
ln f(z) ~ Σckλkt; f(z) ~ c*eΣλkt
|f(z)| ~ |c*| |eλjt| |eΣλkt|; λ = Re λ + i Im λ ; |ei w|= |Cos(Im w) + i Sin(Im w)| = 1
lim |f(z)| ~ lim(|c*| |eRe(λj)t | |eΣλkt|) as t-> ∞
Then, |eRe(λj)t | -> large because Re λj >0 So, z0 is an unstable fixed point.
If all Re λj 0 as t-> ∞ So, z0 is a stable fixed point.
12/20/2001 Systems Dynamics Study Group 21
dx/dt = xg(x) – yp(x) det [f(1)(z0)-λI] =
dy/dt = y[-s + cp(x)] (g(x0)+x0g’(x0)-y0p’(x0)-λ) (-s+cp(x0)-λ)
dz/dt = f(z) Fixed Points – cy0p’(x0)p(x0) = 0
dx/dt = xg(x) – yp(x) = 0
dy/dt = y[-s + cp(x)] = 0
y
1. (0,0), p(0) = 0
2. (k,0), g(k) = 0 (x*,y*)
x
k
3. (x*,y*), 00
12/20/2001 Systems Dynamics Study Group 22
dx/dt = xg(x) – yp(x) det [f(1)(z0)-λI] =
dy/dt = y[-s + cp(x)] (g(x0)+x0g’(x0)-y0p’(x0)-λ) (-s+cp(x0)-λ)
dz/dt = f(z) (0,0) – cy0p’(x0)p(x0) = 0
(g(x0)+x0g’(x0)-y0p’(x0)-λ) (-s+cp(x0)-λ) – cy0p’(x0)p(x0) = 0
x0=y0=p(x0)=0
(g(0)-λ)(-s-λ)=0
λ=g(0),-s; y
λ1= g(0) > 0, corresponds to x axis x
(0,0)
λ2= -s 0, p(x*) 0 y
λ1= kg’(k) 0, corresponds to y axis
(k,0) is unstable x
k
12/20/2001 Systems Dynamics Study Group 24
dx/dt = xg(x) – yp(x) det [f(1)(z0)-λI] =
(g(x0)+x0g’(x0)-y0p’(x0)-λ) (-s+cp(x0)-λ)
dy/dt = y[-s + cp(x)]
dz/dt = f(z)
(x*,y*) – cy0p’(x0)p(x0) = 0
(g(x0)+x0g’(x0)-y0p’(x0)-λ) (-s+cp(x0)-λ) – cy0p’(x0)p(x0) = 0
Note: -s+cp(x0)=0; x0=x*, y0=y*
(g(x*)+x*g’(x*)-y*p’(x*)-λ)(-λ) – cy*p’(x*)p(x*) = 0
λ2 - [g(x*)+x*g’(x*)-y*p’(x*)]λ – cy*p’(x*)p(x*) = 0
B C>0
λ = (B +/– sqrt(B2 + 4C))/2
Note: slope of prey isocline, (dy/dt) at (x*,y*) y
= d(dx/dt)dx = g(x)+xg’(x)-yp’(x) = B
If B > 0, (x*,y*) is unstable
λ1 = [B – sqrt(B2 + 4C)]/2 0
λ2 = [B + sqrt(B2 + 4C)]/2 > 0
(x*,y*)
If B 0
C1 – has a continuous first
(x*,y*)
derivative
x
k
12/20/2001 Systems Dynamics Study Group 26
dx/dt = xg(x) – yp(x)
dy/dt = y[-s + cp(x)]
Poincaré-Bendixon Rationale
F(z+,t1)=z+1=(x1,y1)
F(z+,t2)=z+2=(x2,y2) y
Z–1
Z–2
lim z+k -> z, as k->∞ Z+1
Z+2
Z
B>0
F(z–,t1 )=z–1=(x1,y1)
F(z–,t2)=z–2=(x2,y2)
lim F(z–k) -> z-, as k->∞ x
z- 0
dx/dt = xg(x) – yp(x) (x*,y*)
dy/dt = y[-s + cp(x)], s>0, c>0 x
k
g(x) is a growth function,
g(x), monotonic non-increasing, g’(x) 0, g(k)=0
p(x) is predation function
p(x), monotonic increasing, p’(x) >0 , p(0)=0
(x*,y*) fixed point, x*>0, y*,>0 => x*g(x*)-y*p(x*)=0; -s+cp(x*)=0
B = g(x*)+ x*g’(x*)-yp’(x*) > 0
Then, the dynamical system has a stable periodic orbit.
12/20/2001 Systems Dynamics Study Group 28
Thank you!
12/20/2001 Systems Dynamics Study Group 29
Predator-Prey System
with a stable periodic orbit
Systems Dynamics Study Group
3rd Session – Simulation Results
Ellis S. Nolley
12/20/2001
12/20/2001 Systems Dynamics Study Group 30
References
• McGehee & Armstrong, Journal of Differential Equations, 23, 30-52 (1977)
• Vensim ® PLE software (free for educational use)
www.vensim.com/download.html
• Vensim Tutorial by Craig Kirkwood, Arizona State University
www.public.asu.edu/~kirkwood/sysdyn/SDRes.htm
• Vensim User Guide
www.vensim.com/ffiles/venple.pdf
12/20/2001 Systems Dynamics Study Group 31
Topics
• Overview
– Simple Analysis – 1st session, 11/7/2001
– Rigorous Analysis – 2nd session, 11/27/2001
– Simulation Results – 3rd session, 12/11/2001
• Model Parameters
• Simulation Results
• Vensim Techniques
• Bifurcation
• Extra: Mathematics of Parameter Selection
12/20/2001 Systems Dynamics Study Group 32
y
Model
B>0
dx/dt = xg(x) – yp(x)
(x*,y*)
dy/dt = y[-s + cp(x)], s>0, c>0
x
k
g(x) is a growth function,
g(x), monotonic non-increasing, g’(x) 0, g(k)=0
p(x) is predation function
p(x), monotonic increasing, p’(x) >0 , p(0)=0
(x*,y*) fixed point, x*>0, y*,>0 => x*g(x*)-y*p(x*)=0; -s+cp(x*)=0
B = g(x*)+ x*g’(x*)-yp’(x*) > 0
Then, the dynamial system has a stable periodic orbit.
12/20/2001 Systems Dynamics Study Group 33
dx/dt = xg(x) – yp(x) g(0)>0, g’(x)0, B>0
Model Parameters
g(x) = a0+a1x, x*~147.4, a0=54, a1= -0.15
p(x) = b ln(x+1), b=4, s=200, c=10
g(0) = 54>0, g’(x)= -0.150
B= g(x)+xg’(x)-xg(x)p’(x)/p(x), for x=x*, p(x*)=s/c
~ 54+2(-0.15)147.4+4(1)[54-0.15(147.4)]/(200/10)
since x/(x+1) ~ 1
~ 54 - 44.2 - 6.4 = 3.4 > 0
12/20/2001 Systems Dynamics Study Group 34
Inside Orbit
X & Y Log Inside: time 0 - 40
500
400
Y - Predator
300
200
100
0 Time Series of first 3 Periods
0 100 200 300 400
120
Amount X & Y
X - Prey 100
80
x
60
40 y
20
0
0 0.5 1 1.5 2 2.5
Time
12/20/2001 Systems Dynamics Study Group 35
Outside Orbit
X & Y Log Outside: time 0 - 40
500
400
Y - Predator
300
200
100
0
0 100 200 300 400
X - Prey Time Series of first 2 Periods
500
Amount X & Y
400
300 x
200 y
100
0
0
0.08
0.16
0.24
0.32
0.4
0.48
0.56
0.64
0.72
0.8
0.88
0.96
Time
12/20/2001 Systems Dynamics Study Group 36
Inside & Outside Orbits
X & Y Log Inside: time 0 - 40
500
X & Y Log Outside: time 0 - 40
400
Y - Predator
300
200
500
100
400 0
Y - Predator
0 100 200 300 400
X - Prey
300
200
100
0
0 100 200 300 400
X - Prey
12/20/2001 Systems Dynamics Study Group 37
Combined Inside & Outside Orbits
X & Y Log Inside: time
X & Y Log Outside: time00- -40
40
500
500
400
400
Y - Predator
300
300
200
200
100
100
0
0
0
0 100 200 300
300 400
400
X - Prey
12/20/2001 Systems Dynamics Study Group 38
dx/dt = xg(x) – yp(x)
Vensim Model Layout dy/dt = y[-s + cp(x)]
12/20/2001 Systems Dynamics Study Group 39
g(x)
12/20/2001 Systems Dynamics Study Group 40
p(x)
12/20/2001 Systems Dynamics Study Group 41
dx/dt
12/20/2001 Systems Dynamics Study Group 42
x
12/20/2001 Systems Dynamics Study Group 43
dy/dt
12/20/2001 Systems Dynamics Study Group 44
y
12/20/2001 Systems Dynamics Study Group 45
Other Vensim Techniques
• Select Runge Kutta Integration (RK4).
• Select initial points (x,y)=(1,1) for an outside orbit
and (x,y)=(125,200) for an inside orbit.
• Select 0.005 for a step size in Model/Settings
• Select a custom graph/table to export to Excel
– Control Panel, Graphs, New, Name title, select
variables x & y, click on scale between them
– Click on As Table, click on “running down”
– Click on Ok, close
12/20/2001 Systems Dynamics Study Group 46
Run and Export Text File
• Click on Run Simulation
• Click on Control Panel
• Click on graph name, click on Display
• Click on File, then Save As
12/20/2001 Systems Dynamics Study Group 47
Import Text File into Excel
• Run Excel
• Click Open, select txt type, select file, click Open,
Finish.
• Click on Chart Wizard, XY (scatter), click on Data
Source icon (to right of data range), click and drag
over x & y data, click on Data Source icon, complete
the chart.
• Create a time series chart using t,x,y data the same
way as above, dragging over several periods of data.
• Then, alter step size and initial points in Vensim
• Create other charts for parameter changes by
edit/copy sheet, run new simulation & copy/paste
simulation data onto new sheet’s data region
12/20/2001 Systems Dynamics Study Group 48
Example
Excel
Result
12/20/2001 Systems Dynamics Study Group 49
dx/dt = xg(x) – yp(x) g(0)>0, g’(x)0, B>0
12/20/2001 Systems Dynamics Study Group 50
X & Y Log Outside: time 0 - 40 a0=54 X & Y Log Outside: time 0 - 40 a0=51
500 500
400 400
Y - Predator
Y - Predator
300 300
200 200
100 100
0 0
0 100 200 300 400 0 100 200 300 400
X - Prey X - Prey
X & Y Log Outside: time 0 - 40 a0=40 X & Y Log Outside: time 0 - 40 a0=49.7
250 400
350
200
Y - Predator
300
Y - Predator
150 250
200
100 150
50 100
50
0 0
0 100 200 300 0 100 200 300 400
X - Prey X - Prey
12/20/2001 Systems Dynamics Study Group 51
Projection of
Stable Attractor onto X Axis
x
Actual
~ 147.4 boundary shape
is not described
a0
~ 22.1 ~ 49.7
12/20/2001 Systems Dynamics Study Group 52
dx/dt = xg(x) – yp(x) g(0)>0, g’(x)0
• Generalized Lotka-Volterra Predator-Prey Model
• Internal Fixed Point (x*,y*)
– Stable when B0
• Existence Proof y
• Simulation Results
• Vensim techniques
• Bifurcation B>0
(x*,y*)
x
k
12/20/2001 Systems Dynamics Study Group 53
Another Reference
• C. Neuhauser, “Mathematical Challenges in Spatial Ecology,” Notices
of the American Mathematical Society, 48, 1304-1314 (Dec 2001)
http://www.ams.org/notices/200111/fea-neuhauser.pdf
University of Minnesota, EEB dept of CBS
12/20/2001 Systems Dynamics Study Group 54
Predator-Prey models Competition
Typical Competition Beliefs
• Survival of the fittest
• Competition develops excellence
• Diversity increases stability
• Complexity decreases stability
• One competitor per niche
• Good designs stabilize desirable behavior
and destabilize undesirable behavior
What are likely outcomes of well defined systems?
What systems produce specific outcomes?
12/20/2001 Systems Dynamics Study Group 55
Thank you!
12/20/2001 Systems Dynamics Study Group 56
Extra
Mathematics of Parameter Selection
12/20/2001 Systems Dynamics Study Group 57
dx/dt = xg(x) – yp(x) g(0)>0, g’(x)0, B>0
Recall that any continuous function within a closed bounded region can
be uniformly approximated by polynomials. (Stone-Weierstrauss)
Let g(x) ε R(m), p(x) ε R(n) real polynomials of degree m & n
g(x)=0Σmakxk, a0>0, a10, g’(x)0
p(x)= 1Σnbkxk, b0=0, b1>0, bn>0, since p(0)=0, p’(x)>0 for x>0
B= g(x)+xg’(x)-xg(x)p’(x)/p(x), for x=x*
+ - -
= 0Σmakxk+x(1Σmkakxk-1)-(0Σmakxk)(1Σnkbkxk)/(1Σnbkxk)
= 0Σmakxk+1Σmkakxk -(0Σmakxk)(1Σnkbkxk)/(1Σnbkxk)
= a0[1-(1Σnkbkxk)/(1Σnbkxk)]+ 1Σmkakxk -(1Σmakxk)(1Σnkbkxk)/(1Σnbkxk)
Note: if aj0 & bk>0 for all k>0, then (1Σnkbkxk)/(1Σnbkxk)>1,
then B0 for some j:1= 3
12/20/2001 Systems Dynamics Study Group 58
dx/dt = xg(x) – yp(x) g(0)>0, g’(x)0, B>0
g(x)= 0Σmakxk , a0>0, a10, a10, p(0)=0
Find x=x*, s/c = b ln(x +1), x*= e[s/(bc)] - 1
y*= (a0x+a1x2 )(c/s)
12/20/2001 Systems Dynamics Study Group 59
dx/dt = xg(x) – yp(x) g(0)>0, g’(x)0, B>0
B= g(x)+xg’(x)-xg(x)p’(x)/p(x), for x=x*
= a0 + a1x + a1x – x(a0 + a1x)b/[x+1] /(b [ln (x+1)])
= a0 (1-(cb/s)[x/(x+1)]) + 2a1x - [(a0 + a1x)b[x/(x+1)]/(s/c)], since p(x*)=s/c
= a0 (1-(cb/s)[x/(x+1)]) + a1x(2-(cb/s)[x/(x+1)]) > 0,
a0 > - a1x[(2s-cb)[x/(x+1)]/s][s/(s-cb[x/(x+1)])]
a0 > - a1x[(2s-cb[x/(x+1)])/(s-cb[x/(x+1)])]
2s>cb[x/(x+1)] and s>cb[x/(x+1)] => s>cb[x/(x+1)]
or 2s s cb>cb[x/(x+1)], since x - a1x[(2s-cb[x/(x+1)])/(s-cb[x/(x+1)])]
4) Verify y* = (a0x+a1x2 )(c/s) >0
5) Verify that B>0
12/20/2001 Systems Dynamics Study Group 60
dx/dt = xg(x) – yp(x)
Log (Cont’d)
g(0)>0, g’(x)0, B>0
Model Parameter Selections
1) b=4, c=10, s > bc=40, Let s=200
2) x* = e[200/(4*10)] –1 = e 5 –1 ~ 148.4 – 1 = 147.4
3) a0 > - a1x(2s-cb[x/(x+1)])/(s-cb[x/(x+1])
= - a1*147.4(2*200 - 40[1])/(200 - 40[1]) , since [x/(x+1)]~1
= - a1*147.4(360/160),
= - a1(332.7)
let a1= -0.15, a0 > 49.7 , Let a0 > 54
4) y* = x*g(x*)/p(x*) = 147.4(200-0.15*147.4)/(200/10) =
~ 235.0 > 0
5) B = g(x) + xg’(x) – xg(x)p’(x)/p(x) =
= (a0+a1x) + 2a1x – b[x/(x+1)] (a0+a1x)
= [54 – 0.15(147.4)] + 4(-0.15)[54 – 0.15(147.4)] , since [x/(x+1)]~1
= 54 – 44.2 – 6.4 = 3.4 > 0
12/20/2001 Systems Dynamics Study Group 61
Thank you!
12/20/2001 Systems Dynamics Study Group 62