{VERSION 5 0 "SUN SPARC SOLARIS" "5.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 1 }{CSTYLE "" -1 256 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 }
{CSTYLE "" -1 257 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 
258 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 259 "" 0 1 0 0 
0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 260 "" 1 12 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Helvetica" 1 12 
0 0 0 0 1 2 2 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE 
"Heading 1" 0 3 1 {CSTYLE "" -1 -1 "" 1 12 0 0 0 0 0 1 0 0 0 0 0 0 0 
1 }1 0 0 0 6 6 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 2" 3 4 1 {CSTYLE "" 
-1 -1 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }0 0 0 -1 4 4 0 0 0 0 0 0 
-1 0 }{PSTYLE "R3 Font 0" -1 256 1 {CSTYLE "" -1 -1 "Helvetica" 1 12 
0 0 0 0 2 1 2 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE 
"R3 Font 2" -1 257 1 {CSTYLE "" -1 -1 "Courier" 1 14 0 0 0 0 2 2 2 0 
0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "R3 Font 2" -1 
258 1 {CSTYLE "" -1 -1 "Courier" 1 12 0 0 0 0 2 2 2 0 0 0 0 0 0 1 }0 
0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }}
{SECT 0 {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT 257 13 "blowflies.mws" }{TEXT -1 47 "   Maple worksheet t
o accompany Hassell, 1978, " }{TEXT 258 47 "The dynamics of arthropod \+
predator-prey systems" }{TEXT -1 35 ", Monographs in Population Biolog
y " }{TEXT 259 2 "13" }{TEXT -1 40 ", and the course text Edelstein-Ke
shet, " }{TEXT 256 30 "Mathematical Models in Biology" }{TEXT -1 61 ",
 Birkhauser Mathematics Series, McGraw Hill, New York, 1988." }}{PARA 
0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 871 "The first part
 explores the standard density independent and dependent models.  The \+
students are asked to investigate the qualitative effects on the  beha
vior of the populations under modification of relevant parameters.  Th
ey are also asked to examine the eigenvalues of the linearized system \+
at the equilibrium point, at this point purely numerically.  The secon
d part is concerned with a more general (algebraic) eigenvalue analysi
s of both the original and the desity dependent models.  After lineari
zation eigenvalues are computed and interpreted, with the goal of unde
rstanding how behaviour of the model is affected by parameter values. \+
There is a (we think rather nice!) demonstration that Maple can handle
 parameter space plots. The last part investigates the stability analy
sis from a graphical pont of view, giving an idication why the jacobia
n matrix is used." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{SECT 1 {PARA 3 
"" 0 "" {TEXT -1 48 "Section 1a. Density independent model (missing!)
" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 36 "Section 1b.  Density depende
nt model" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 307 "Enter values of growth rate r and density dependence par
ameter q.  The program will generate plots of populations vs. time and
 phase plots.  It also will calculate the eigenvalues of the linearize
d equations around the equilibrium values of host and parasitoid for t
he particular values of r and q chosen.  " }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 37 "restart;  with(plots):  with(linalg):" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 55 "r:=1.5;  q:=0.6;  c:=1.0;   a:=0.2;  # paramet
er values" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 141 "txtr:=convert
(evalf(r),string): txtq:=convert(evalf(q),string):\nlabelr:=`r =  `: l
abelq:=`,  q = `: fulltitle:= cat(labelr,txtr,labelq,txtq):" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 77 "P_hat:= r/a*(1-q); N_hat:= P
_hat/(1-exp(-a*P_hat)); K:= N_hat/q; # equilibria" }}}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 43 "N[0]:= 11; P[0]:= 4;  # initial conditions \+
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 113 "for i from 0 to 200 do
\n   N[i+1]:=N[i]* exp(r * (1 - N[i]/K) - a*P[i]);\n   P[i+1]:=c*N[i]*
(1 - exp(-a*P[i]));\nod:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 88 
"plot(\{[[t,P[t]] $t = 0.. 200], [[t,N[t]] $t = 0..200] \}, style=line
, title=`fulltitle`);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 89 "pl
ot([[N[t],P[t]] $ t=0..200], style=line, title=cat(`Phase portrait for
 `,`fulltitle`));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 237 "\nCompare the results of these si
mulations to Figure 2.7 in Hassell.  Now go back and try this again wi
th new values for r and q.  You should use the Hassell reprint as a gu
ide to pick values that yield qualitatively different behaviors. " }}}
}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 33 "Section 1c.  Equilibrium analysi
s" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 465 "\nTo calculate the eigenvalu
es of the linearized system, we recast the equations in more symbolic \+
terms, converting N(t) to N and P(t) to P.  We then calculate the part
ial derivatives of each equation with respect to both N and P.   This \+
can be done by hand, as done by calculating a11, a12, a21, a22.  Alter
natively, the jacobian command will do the calculation directly.  Once
 the matrix is assembled, the eigenvals command is used to calclulate \+
the eigenvalues.  " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 242 "Remember that since this is a difference equation model,
 stability occurs  (possibly after damped oscillation) when the eigenv
alues are less than 1 in absolute value.  Oscillations occur when ther
e are imaginary components to the eigenvalues." }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 17 "r:=1.5;  q:=0.6; " }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 41 "Host_rate:= N * exp(r * (1 - N/K) - a*P);" }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 53 "a11:= diff(Host_rate, N);   a12:= diff(Ho
st_rate, P);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "A11:= evalf(subs(\{
N=N_hat, P=P_hat\}, a11));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "A12:=
 evalf(subs(\{N=N_hat, P=P_hat\}, a12));" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 41 "Parasitoid_rate:= c * N* (1 - exp(-a*P));" }}{PARA 0 
"> " 0 "" {MPLTEXT 1 0 61 "a21:=diff(Parasitoid_rate, N); a22:=diff(Pa
rasitoid_rate, P);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "A21:= evalf(s
ubs(\{N=N_hat, P=P_hat\}, a21));  " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
43 "A22:= evalf(subs(\{N=N_hat, P=P_hat\}, a22));" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 47 "A:= matrix( 2, 2, [ [A11, A12], [A21, A22] ] );" }}
}{EXCHG {PARA 0 "" 0 "" {TEXT -1 94 "\nHere is the alternate method fo
r computing the coefficient matrix using the jacobian command." }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 96 "J:= jacobian([Host_rate, Parasitoid
_rate], [N, P]); # computes the partials and forms the matrix" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "A:= evalf(subs(\{N=N_hat, P=
P_hat\}, evalm(J)));  # same as before" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 22 "eigens:= eigenvals(A);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 97 "abs(eigens[1]); # this gives the magnitude (\"absolut
e value\") of the first of the two eigenvalues" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 220 "How do you interpret the behavior of the
 system if these are the eigenvalues?  Go back to the beginning of thi
s part of the worksheet, reset r and q, and see if the eigenvalue calc
ulations lead to a different conclusion." }}}}{SECT 1 {PARA 3 "" 0 "" 
{TEXT -1 59 "Section 2.  Eigenvalue analysis -- dependence on paramete
rs" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT 
-1 369 "We exhibit  dependence on parameters of the behaviour of the s
ystem near the equilibrium point in both the original and density depe
ndent Nicholson-Bailey models.  In the original model, host growth is \+
unconstrained in the absence of parasitoids; the rate is lambda.  We u
se the coefficient matrix A obtained by linearization at the equilibri
um (see text E-K, page 82)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 36 "restart: with(plots):  with(linalg):" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 81 "c:= 1;  # to get formal rath
er than floating point calculation of the eigenvalues" }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 34 "mu:= 1 - 1 / lambda;   b:= c * mu;" }}{PARA 0 ">
 " 0 "" {MPLTEXT 1 0 63 "A:= matrix( 2 , 2, [ [1, -ln(lambda)/b], [b, \+
ln(lambda)/mu] ]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 80 "It
 is easy to get the eigenvalues, but watch out, they could be complex \+
numbers." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "evalues:= eigenvals(A);
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 
697 "These are a big ugly mess, at least inside the square root.  You \+
can actually do a good bit of housecleaning (by hand) to see that the \+
first one at least has to have a real part that is bigger than 1,  plu
s an imaginary part.  This is already enough to demonstrate instabilty
 of the equilibrium point.  However, to get an idea what is going on, \+
I decided to substitute some values for lambda.  Remember evalf conver
ts to decimal form.  We can also compute the magnitudes of these compl
ex numbers, but it is clear that they are bigger than 1.  We begin by \+
computing the the eigenvalue magnitude that corresponds to a small rep
roductive rate lambda; then we try larger and larger reproductive rate
s." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "evalf(subs( lambda=1.001, [ev
alues] )); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 99 "abs( % [1] )
;  # take the first eigenvalue and compute its magnitude (both have th
e same magnitude)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "evalf(
 subs( lambda=1.1, [evalues] ) );   abs( % [1] );" }}}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 55 "evalf( subs( lambda=1.5, [evalues] ) );   ab
s( % [1] );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "evalf( subs(
 lambda=2.0, [evalues] ) );   abs( % [1] );" }}}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 53 "evalf( subs( lambda=10, [evalues] ) );  abs( % [1]
 );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "evalf( subs( lambda=
1000, [evalues] ) );   abs( % [1] );" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 391 "\nApparently in every case (except for the last one) we get tr
ajectories in the phase plane that spiral away from the equilibrium.  \+
Since spiralling in the phase plane (parasite population against host \+
population) represents oscillation  of both populations as functions o
f time, what this amounts to is bigger and bigger oscillations in the \+
populations--just what you observed in Figure 3.3!" }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 202 "Let's try to analyze the
 eigenvalues abstractly.  A little algebraic hacking around shows that
 if you bring the denominator inside the square root,  and do some rea
rranging and factoring, the result is:" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 86 "mess:=1-6*ln(lambda)*(lambda-2/3)/(lambda-1) + (ln(lambda))^2 \+
*lambda^2 /(lambda-1)^2;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 40 "Let's plot this as a function of lambda.
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "plot( mess, lambda = 1 .. 400 )
;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 
218 "As you can see, for large enough lambda the eigenvalues will be r
eal (since mess is positive and has a real square root), but for lambd
a in any sensible range, the situation that we observed in Figure 3.3 \+
is  typical.\n" }}{PARA 0 "" 0 "" {TEXT -1 430 "Now let's suppose the \+
host population is subject to density dependent growth.  To linearize \+
the system we need to compute partial derivatives and evaluate them at
 the equilibrium values of N and P.  We use the right hand side of equ
ations (28) on page 84 of E-K.  Recall that ultimately we want to look
 at the dependence of the eigenvalues on r and q only, so let us subst
itute the values used in the text: c = 1, K = Nequil / q ." }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "a_11:= diff( N * exp( r * (1 - N / \+
K ) - a * P) , N) ;  " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "a_12:= dif
f( N * exp( r * (1 - N / K ) - a * P) , P );  " }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 51 "a_21:= diff(c *  N * ( 1 - exp( -a * P ) ) , N ) ; " 
}}{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "a_22:= diff( c * N * ( 1 - exp( -
a * P ) ) , P ) ;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 71 "Pequil
:= r / a * ( 1 - q ); Nequil:= Pequil / (1 - exp(- a * Pequil)); " }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "K:= Nequil / q;   c:= 1;" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "A_11:= subs(\{N = Nequil, P \+
= Pequil\}, a_11);  " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 90 "A_12:= subs
(\{N = Nequil, P = Pequil\}, a_12);\nA_21:= subs(\{N = Nequil, P = Peq
uil\}, a_21); " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "A_22:= subs(\{N =
 Nequil, P = Pequil\}, a_22);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 0 "" }{TEXT -1 0 "" }{MPLTEXT 1 0 53 "A:= matrix( 2, 2, [  [ A_11, A
_12] , [A_21, A_22] ]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 179 "This gives us the linearization coeffici
ent matrix.   Note that the parameter a has disappeared.  Next we give
 an alternative approach using more of Maple's built-in capabilities.
" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 46 "host_rhs:= N * exp( r * (1 - N / K ) - a * P);" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "parasitoid_rhs:= c *  N * ( 1 - exp
( -a * P ) ) ;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "J:= jacob
ian( [host_rhs, parasitoid_rhs], [N, P] );" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 85 "B:= subs(\{N=Nequil, P=Pequil\}, evalm(J));  # eval
uate the Jacobian at the equilibrium" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 60 "The same as A
?  Use whichever syntax you prefer from now on!" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 384 "\nOnc
e again, we perform an eigenvalue analysis.  Substitution of various v
alues of r and q  should yield results consitstent with Figure 2.7 of \+
the Hassell handout.  If we look in (r, q) space the boundary curves o
f the regions should be determined by seeing where these eigenvalues h
ave imaginary part equal zero (the internal curve) or magnitude equal \+
to one (the outside boundary)." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 76 "e
values:= [eigenvals(A)]; # Note there are two of them, separated by a \+
comma" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 105 "The eigenvalues involve the parameters a, r, and q; let \+
us substitute the value of a given in the text.  " }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 36 "mu:= evalues[1];  eta:= evalues[2] ;" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 664 "\nThe following procedure examines the e
igenvalues and tests to see if the magnitude is greater than or less t
han one, and whether there is an imaginary component or not.  Dependin
g upon which of these conditions is met, the procedure returns a numbe
r between 0 and 4.  Note that it is essentially impossible to test equ
ality of floating point (decimal) numbers on a computer (you can test \+
if an integer is zero or not, but whether 10E-8 or 10E-10 or 10E-12 is
 zero will depend on the computer hardware and software).  For this re
ason, we test instead for \"close enough\" to 0 or 1.  The input argum
ents rr and qq are used to compute the value of the eigenvalue mu. " }
}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "behav:= proc(rr, qq) " }}{PARA 0 "
> " 0 "" {MPLTEXT 1 0 449 " ev:=evalf( subs( \{r=rr, q=qq\},  mu) );\n
  evi:=Im(ev);  evmag:=abs(ev);\n  if evmag > 1 and abs(evi) <10^(-8) \+
then zone:=4; fi; #  explosion (red)\n  if evmag > 1 and evi <>0  then
 zone:=3; fi; # unstable oscillation (green)\n  if abs(evmag - 1)<10^(
-8) then zone:=2; fi; # stable oscillation (yellow)\n  if evmag < 1 an
d evi <> 0 then zone:=1; fi; # damped oscillation (blue)\n  if evmag <
 1 and abs(evi) <10^(-8) then zone:=0; fi;#over-damping (purple)" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "result:= zone;" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 4 "end:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 57 "First let's test out the procedure on som
e sample values." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "behav( 6, 0.65)
;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "behav( 1, 0.3 );" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "behav( 2, 0.5);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "behav( 3 , 0.6 );" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 16 "behav( 1, 0.9 );" }}}{EXCHG {PARA 0 "" 0 
"" {TEXT -1 123 "\nNow use the procedure to reproduce Hassell Fig 2.7.
  The orientation is set so that the view of the 3D plot is from above
." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 195 "contourplot3d(behav, 0.1.. 6.
0, 0.01.. 0.9, grid=[45,45], orientation = [-90, 0], shading=ZHUE, axe
s=BOXED, style = PATCHCONTOUR, title=`Stability Regions of Density Dep
endent N-B: Hassell 2.7`);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 56 "Show the imaginary part of mu as a functi
on of r and q. " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 136 "Implot:= plot3d
(Im(mu), q=0.1.. 0.9, r=0.1.. 6.0, axes=BOXED, shading=ZHUE, style=PAT
CHCONTOUR, title=`Imaginary Parts of Eigenvalues`):" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 18 "display( Implot );" }}}{EXCHG {PARA 0 "" 
0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 51 "Show the magnitude o
f mu as a function of r and q. " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 130 
"Magplot:=plot3d(abs(mu), q=0.1..0.9, r=0.1..6.0, shading=ZHUE, axes=B
OXED, style=PATCHCONTOUR, title=`Magnitudes of Eigenvalues`):" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "display( Magplot );" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 50 "Sh
ow the real part of mu as a function of r and q." }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 129 "Replot:= plot3d(Re(mu), q=0.1..0.9, r=0.1..6.0, shad
ing=ZHUE, axes=BOXED, style=PATCHCONTOUR, title=`Real Parts of Eigenva
lues`):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "display( Replot \+
);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
-1 47 "Display real and imaginary components together." }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 67 "display( \{ Replot, Implot \} , title = `Real \+
and imaginary parts `);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 194 "Here is an alternate way to get the boun
dary curves.  Caution:  Maple's implicit plotting is impressionistic a
t best.  It is very senstitive to change in mesh fineness, and it can \+
be very slow. " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "magnitude:= abs(m
u);  " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 76 "Outer:= implicitplot( magn
itude = 1, r=0.1..6.0, q=0.01..0.9, color = plum):" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 107 "Inner:= implicitplot( magnitude = Re(mu), \+
r=0.1..6.0, q=0.01..0.9, grid=[50,50], color=green): # Patience!!" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "display( \{Outer, Inner\} );
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 4 "
" 0 "" {TEXT -1 12 "Section 3.  " }{TEXT 260 44 "Graphical illustratio
n of stability analysis" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "restart: gc(): " }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "with(plots):" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 26 "eq1:=n=n*lambda*exp(-a*p);" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 25 "eq2:=p=c*n*(1-exp(-a*p));" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 27 "parms:=a=.068,c=1,lambda=2;" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "Eq1:=subs(parms,eq1);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "Eq2:=subs(parms,eq2);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "equil:=solve(\{Eq1,Eq2\},\{n,p\});
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "plot3d(rhs(Eq1),p=0..30
,n=0..30,axes=boxed,title=`Host Growth Function`);" }}}{EXCHG {PARA 0 
"" 0 "" {TEXT -1 358 "The above graph shows the host population at tim
e t+1 as a function of the host and parasitoid populations at time t. \+
 This surface is analogous to the growth function lines used in the co
bweb analysis of the logistic equation by May and Oster.  It is a surf
ace rather than a line, because both hosts and parasitoids influence t
he growth of the host species." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 62 "plot3d(n,p=0..30,n=0..30,axes=boxed,title=`Host Zero Growth`);
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 187 "The above plane represents t
he host population at time t+1 which equals the host population at tim
e t (equivalent to the 45 degree line plotted by May and Oster in thei
r cobweb analysis)." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 79 "plot
3d(rhs(Eq2),p=0..30,n=0..30,axes=boxed,title=`Parasitoid Growth Functi
on`);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 176 "The above graph represe
nts the parasitoid population at time t+1 as a function of host and pa
rasitoid populations at time t.  Note that this growth function is ver
y non-linear." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 68 "plot3d(p,p
=0..30,n=0..30,axes=boxed,title=`Parasitoid Zero Growth`);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 152 "The above graph represents the parasitoi
d population at t+1 when it equals the population size at time t (equi
valent to May and Oster's 45 degree line)." }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 272 "The entrie
s in the jacobian matrix j below are the partial derivatives of the ho
st and parasitoid functions with respect to host and parasitoid popula
tion density.  They represent the slopes of the host and parasitoid gr
owth functions along the host and parasitoid axes.  " }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 14 "with (linalg):" }}}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 39 "j:=jacobian([rhs(Eq1),rhs(Eq2)],[n,p]);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 108 "The jacobian is evaluated at the \+
equilibrium point, and its eigenvalues describe the behavior of the sy
stem." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "jj:=subs(op(2,[equ
il]),evalm(j));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "ev:=eige
nvals(jj);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 71 "The abs function is
 used to determine the magnitude of the eigenvalues." }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 32 "abs(op(1,[ev])),abs(op(2,[ev]));" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 327 "In this case, the eigenvalues hav
e  magnitudes greater than 1, indicating growth away from the equilibr
ium.  The eigenvalues also have imaginary components indicating oscill
atory behavior.  This combination indicates expanding oscillations.  T
he parameters used are from Fig 3.3 in Edelstein-Keshet and Fig 2.3 in
 Hassell 1978." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "op(1,op(2
,[equil]));op(2,op(2,[equil]));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
65 "We will plot the functions close to the equilibrium listed above:
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 97 "h_n:=plot3d(rhs(Eq1),p=
rhs(op(1,op(2,[equil])))..rhs(op(1,op(2,[equil])))+.1,n=0..30,axes=box
ed):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 97 "h_p:=plot3d(rhs(Eq1
),p=0..30,n=rhs(op(2,op(2,[equil])))..rhs(op(2,op(2,[equil])))+.1,axes
=boxed):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "display(\{h_n,h
_p\},title=`Host Function Near Equilibrium`);" }}}{EXCHG {PARA 0 "" 0 
"" {TEXT -1 240 "The above figure plots the host growth function along
 lines parallel to the host and parasitoid axes and passing through th
e equilibrium:  The slopes of these lines are the partial derivatives \+
of the growth function with respect to h and p." }}}{EXCHG {PARA 0 "" 
0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 97 "p_n:=p
lot3d(rhs(Eq2),p=rhs(op(1,op(2,[equil])))..rhs(op(1,op(2,[equil])))+.1
,n=0..30,axes=boxed):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 97 "p_
p:=plot3d(rhs(Eq2),p=0..30,n=rhs(op(2,op(2,[equil])))..rhs(op(2,op(2,[
equil])))+.1,axes=boxed):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
64 "display(\{p_n,p_p\},title=`Parasitoid Function Near Equilibrium`);
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 246 "The above figure plots the p
arasitoid growth function along lines parallel to the host and parasit
oid axes and passing through the equilibrium:  The slopes of these lin
es are the partial derivatives of the growth function with respect to \+
h and p." }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 ">
 " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 47 "Sec
tion 4.  Incorporating a refuge for the host" }}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 37 "restart;  with(plots):  with(linalg):" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "params:= r=2.80,  b=0.10, c=1.0, a=
0.2, K=100; # parameter values" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 56 "N[i+1]:=r*N[i]* (b*K/N[i] +(1-b*K/N[i])*exp( - a*P[i]));" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "P[i+1]:=c*N[i]*(1-b*K/N[i])*(1 - ex
p(-a*P[i]));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "N[0]:= 40; \+
P[0]:= 20;  # initial conditions " }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 163 "for i from 0 to 50 do\n   N[i+1]:=subs( params,r*N[i
]* (b*K/N[i] +(1-b*K/N[i])*exp( - a*P[i])));\n   P[i+1]:=subs(params,c
*N[i]*(1-b*K/N[i])*(1 - exp(-a*P[i])));\nod:" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 33 "final_N:= N[50]; final_P:= P[50];" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "plot(\{[[t,P[t]] $t = 0.. 50], [[t,
N[t]] $t = 0..50] \}, style=line);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "plot([[
N[t],P[t]] $ t=0..50], style=line);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 34 "\nHere is the
 equilibrium analysis." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "h
ost:= r*N*(b*K/N +(1-b*K/N)*exp( - a*P));" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 43 "parasitoid:= c*N*(1-b*K/N)*(1 - exp(-a*P));" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 84 "J:= jacobian([host, parasito
id], [N, P]); # compute the partials and form the matrix" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "A:= evalf(subs(\{N=final_N, P=final
_P, params\}, evalm(J))); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
13 "eigenvals(A);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "abs(%[
1]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 
"" 0 "" {TEXT -1 233 "How do you interpret the behavior of the system \+
if these are the eigenvalues?  Go back to the beginning of this part o
f the worksheet, reset the parameter values, and see if the eigenvalue
 calculations lead to a different conclusion." }}}{EXCHG {PARA 0 "" 0 
"" {TEXT -1 0 "" }}}}}{MARK "0 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 
1 }{PAGENUMBERS 0 1 2 33 1 1 }