The attached worksheet shows how to evaluate and graphically analyze an autonomous first-order nonlinear recurrence with two dependent variables and multiple symbolic parameters. 

This worksheet shows how a small module that simply encapsulates the given information of a problem combined with some use statements can greatly facilitate the organization of one's work, can encapsulate the setting of parameter values, and can allow one to work with symbolic parameters.

Edit: In the first version of this Post, I forgot to include the qualifier "autonomous".  The system being autonomous substantially simplifies its treatment.

Autonomous first-order nonlinear recurrences with parameters and multiple dependent variables

Author: Carl Love <> 20-Oct-2018


The techniques used in this worksheet can be applied to most autonomous first-order nonlinear recurrences with multiple dependent variables and parameters.


This worksheet shows how a small module that simply encapsulates the given information of a problem combined with some use statements


can greatly facilitate the organization of one's work,


can encapsulate the setting of parameter values,


can allow one to work with symbolic parameters.


A Problem from MaplePrimes: A discrete Lottka-Volterra population model is applied to an isolated island with a population of predators (foxes), R, and prey (rabbits), K. [Note that R is the foxes, not the rabbits! Perhaps this problem statement originated in another language.] The change over one time period is given by

K[n+1]:= K[n]*(-b*R[n]+a+1);  R[n+1]:= R[n]*(b*e*K[n]-c+1),

where a, b, c, e are parameters of the model. In this problem we will use a= 0.15, b= 0.01, c= 0.02, e= 0.01, when numeric values are needed.


a) Show that there exists an equilibrium (values of K[n] and R[n] such that K[n+1] = K[n] and R[n+1] = R[n]).


b) Write Maple code that solves the recurrence numerically. Assume that if any population is less than 0.5 then it has gone extinct and set the value to 0. Check that your program is idempotent at the equilibrium.



We begin by collecting all the given information (except for specific numeric values) into a module. The ModuleApply lets the user set the numeric values later.


For all two-element vectors used in this worksheet, K is the first value and R is the second value.

KandR:= module()
   a, b, c, e, #parameters

   #procedure that lets user set parameter values:
   ModuleApply:= proc({
       a::algebraic:= KandR:-a, b::algebraic:= KandR:-b,
       c::algebraic:= KandR:-c, e::algebraic:= KandR:-e
   local k;
      for k to _noptions do thismodule[lhs(_options[k])]:= rhs(_options[k]) od;
   end proc,

   Extinct:= (x::realcons)-> `if`(x < 0.5, 0, x) #force small, insignificant values to 0
   #Procedure that does one symbolic iteration
   #(Note that this procedure uses Vector input and output.)
   iter_symb:= KR-> KR *~ <-b*KR[2]+a+1, b*e*KR[1]-c+1>, 

   #Such simple treatment as above is only possible for autonomous

   iter_num:= Extinct~@iter_symb #one numeric iteration
end module:

#The following expression is the discrete equivalent of the derivative (or gradient).
#It represents the change over one time period.
P:= <K,R>:  
OneStep:= KandR:-iter_symb(P) - P

Vector(2, {(1) = K*(-R*b+a+1)-K, (2) = R*(K*b*e-c+1)-R})

#An equilibrium occurs when the gradient is 0.
Eq:= <K__e, R__e>:
Eqs:= solve({seq(eval(OneStep=~ 0, [seq(P=~ Eq)]))}, [seq(Eq)]);

[[K__e = 0, R__e = 0], [K__e = c/(b*e), R__e = a/b]]

#We're only interested here in nonzero solutions.
EqSol:= remove(S-> 0 in rhs~(S), Eqs)[];

[K__e = c/(b*e), R__e = a/b]

#Set parameters:
KandR(a= 0.15, b= 0.01, c= 0.02, e= 0.01);

#Show idempotency at equilibrium:
use KandR in Eq0:= eval(Eq, EqSol); print(Eq0 = iter_num(Eq0)) end use:

(Vector(2, {(1) = 200.0000000, (2) = 15.00000000})) = (Vector(2, {(1) = 200.0000000, (2) = 15.00000000}))

#procedure that fills a Matrix with computed values of a 1st-order recurrence.
#(A more-efficient method than this can be used for linear recurrences.)
#This procedure has no dependence on the module.
Iterate:= proc(n::nonnegint, iter, init::Vector[column])
local M:= Matrix((n+1, numelems(init)), init^+, datatype= hfloat), i;
   for i to n do M[i+1,..]:= iter(M[i,..]) od;
end proc:

We want to see what happens if the initial conditions deviate slightly from the equilibrium. It turns out that any deviation (as long as the
initial values are still nonnegative!) will cause the same effect. I simply chose the deviation <7,2> because it was the smallest for which

the plot clearly showed what happens using the scale that I wanted to show the plot at. By using a finer scale, it is possible to see the

"outward spiral" efffect from even the tiniest deviation.

dev:= <7,2>:
use KandR in KR:= Iterate(1000, iter_num, Eq0 + dev) end use:

       KR, #trajectory of population
       KR[[1,1],..], #1st point
       KR[-[1,1],..], #last point,
       <Eq0|Eq0>^+, #equilibrium
       #every 100th point (helps show time scale):
       KR[100*[$1..iquo(numelems(KR[..,1]), 100)-1], ..]
   #This group of options are all lists, each element of which corresponds
   #to one of the above components of the plot:
   style= [line, point$4],
   symbol= [solidcircle$4, soliddiamond],
   symbolsize= [18$4, 12],
   color= [black, green, red, brown, blue],
   thickness= [0$5],
   legend= [`pop.`, init, final, equilibrium, `100 periods`],

   #This group of options are lists, each element of which corresponds to one
   #coordinate axis (horizontal, then vertical).
   view= [0..max(KR[..,1]), 0..max(KR(..,2))],
   labels= [rabbits, foxes],
   labeldirections= [horizontal,vertical],
   size= [700,700], #measured in pixels

   #options applied to whole plot:
   labelfont= [TIMES, BOLDITALIC, 14],
   title= "Population of foxes and rabbits over time" "\n", titlefont= [TIMES,16],
      "\n" "Choosing an initial point near the equilibrium causes"
      "\n" "outward spiraling divergence." "\n",
   gridlines= false

A fieldplot helps show what happens for any starting values. An arrow is drawn from each of a 2-D grid of point. The magnitude and direction of the arrow show the gradient (as a vector) in this case.

   K= 0..max(KR[..,1]),  R= 0..max(KR[..,2]), grid= [16,16],

   #arrow-specific options:
   anchor= tail, fieldstrength= log, arrows= slim, color= "DarkGreen",

   #other options (same as any 2D plot):
   labels= [rabbits, foxes], labeldirections= [horizontal,vertical],
   labelfont= [TIMES, BOLDITALIC, 14],
   title= "One-step population changes from any point" "\n", titlefont= [TIMES,16],
   caption= "\n" "All trajectories spiral outward from the equilibrium." "\n",
   size= [700,700],
   gridlines= false

The above plot is computed only from the symbolic discrete gradient expression OneStep; it does not use the computed population values from the first plot. It only uses the maxima of those computed values to determine the length of the axes.


Conclusion: While this is interesting stuff mathematically, and makes for great plots, divergence from the equilibrium doesn't seem realistic to me.




Please Wait...