(*^
::[ Information =
"This is a Mathematica Notebook file. It contains ASCII text, and can be
transferred by email, ftp, or other text-file transfer utility. It should
be read or edited using a copy of Mathematica or MathReader. If you
received this as email, use your mail application or copy/paste to save
everything from the line containing (*^ down to the line containing ^*)
into a plain text file. On some systems you may have to give the file a
name ending with ".ma" to allow Mathematica to recognize it as a Notebook.
The line below identifies what version of Mathematica created this file,
but it can be opened using any other version as well.";
FrontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.2";
MacintoshStandardFontEncoding;
fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e8, 24, "Times";
fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 18, "Times";
fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, e6, 14, "Times";
fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20, 18, "Times";
fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15, 14, "Times";
fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12, 12, "Times";
fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times";
fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times";
fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-5, 12, "Courier";
fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier";
fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, L-5, 12, "Courier";
fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier";
fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B65535, L-5, 12, "Courier";
fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, 12, "Courier";
fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, 10, "Geneva";
fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7, 12, "Times";
fontset = leftheader, inactive, L2, 12, "Times";
fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7, 12, "Times";
fontset = leftfooter, inactive, L2, 12, "Times";
fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times";
fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times";
fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times";
fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times";
fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times";
fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times";
fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times";
fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times";
paletteColors = 128; automaticGrouping; currentKernel;
]
:[font = section; inactive; preserveAspect; startGroup]
CHAPTER 5
:[font = subsection; inactive; Cclosed; preserveAspect; startGroup]
1
:[font = text; inactive; preserveAspect]
The model with b* = 45 can be simulated for 100 generations starting with both types at densities of 0.5 as follows:
:[font = input; preserveAspect]
f[{x_,y_}]:=(P0=.02 Exp[-.01(100 x + bstar y)];
{100 P0 x,(bstar P0 + 0.5)y})
:[font = input; preserveAspect]
bstar=45;
:[font = input; preserveAspect]
res=NestList[f,{0.5,0.5},100]
:[font = text; inactive; preserveAspect; endGroup]
The model can be simulated with a different value of b* by reassigning the value of this variable. (For this to work it is important that bstar should not have a value when the function f[] is defined, otherwise this value will be incorporated in the definition and cannot be changed.)
;[s]
5:0,0;138,1;143,0;186,1;189,0;286,-1;
2:3,13,9,Times,0,12,0,0,0;2,13,9,Times,1,12,0,0,0;
:[font = subsection; inactive; Cclosed; preserveAspect; startGroup]
2
:[font = text; inactive; preserveAspect]
The model with b* = 70 can be simulated for 100 generations starting with both types at densities of 0.5 as follows:
:[font = input; preserveAspect]
Clear[bstar]
:[font = input; preserveAspect]
f[{x_,y_}]:=(P0=.02 Exp[-.01(100 x + bstar y)];
P1=.5 Exp[-1.0(x+y)];
{100 P0 x,(bstar P0 + P1 )y})
:[font = input; preserveAspect]
bstar=70;
:[font = input; preserveAspect]
res=NestList[f,{0.5,0.5},100]
:[font = text; inactive; preserveAspect]
To determine the criterion for perennials to invade a population of annuals, first determine the equilibrium number of annuals, then P1 and then b*/(1ŠP1). This is easily done with pencil and paper: the equilibrium number of annuals is log 2, P1 = 0.5 exp(Šlog 2) = 0.25, b*/(1ŠP1) = 4b*/3 = 93.33 when b* = 70. Since b = 100 > 93.33, annuals can resist invasion by a few perennials.
;[s]
9:0,0;134,1;135,0;152,1;153,0;244,1;245,0;279,1;280,0;385,-1;
2:5,13,9,Times,0,12,0,0,0;4,18,12,Times,64,10,0,0,0;
:[font = text; inactive; preserveAspect]
To determine the criterion for annuals to invade a population of perennials, first determine the equilibrium number of perennials, then P1 and then b*/(1ŠP1). This is a little more complicated to calculate. First determine the equilibrium number of perennials:
;[s]
5:0,0;137,1;138,0;155,1;156,0;261,-1;
2:3,13,9,Times,0,12,0,0,0;2,18,12,Times,64,10,0,0,0;
:[font = input; preserveAspect]
Clear[bstar]
:[font = input; preserveAspect]
g[y_]=.02 Exp[-.01 bstar y] bstar + .5 Exp[-y];
:[font = input; preserveAspect]
bstar=70;
:[font = input; preserveAspect]
equil=FindRoot[g[y]==1,{y,.5}]
:[font = text; inactive; preserveAspect]
Now complete the calculations:
:[font = input; preserveAspect]
P1=.5 Exp[-y]/.equil
:[font = input; preserveAspect]
bstar/(1-P1)
:[font = text; inactive; preserveAspect; endGroup]
Repeat these calculations with b* = 77 and 85.
:[font = subsection; inactive; Cclosed; preserveAspect; startGroup]
3
:[font = text; inactive; preserveAspect; endGroup]
This exercise is best done with pencil and paper.
:[font = subsection; inactive; Cclosed; preserveAspect; startGroup]
4
:[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup]
4a
:[font = text; inactive; preserveAspect]
First define a function which replaces w[t+1] and q[t+1] in terms of their values in week t according to the recurrence relation:
:[font = input; preserveAspect]
func[t_]:={w[t+1]->.35 u[t] w[t] + .84 w[t],
q[t+1]->.35(1-u[t]) w[t] + q[t]}
:[font = text; inactive; preserveAspect]
We want to find the values of the control variable u[t] (between 0 and 1) which maximize q[22], the number of queens in week 22 at the end of the season, assumed to begin in week 0. Define the objective function as j=q[22], and then evaluate it in terms of values in week 21:
:[font = input; preserveAspect]
j=q[22];
j=Expand[j/.func[21]]
:[font = text; inactive; preserveAspect]
j is maximized by setting u[21]=0. We do this and evaluate j in terms of values in week 20:
:[font = input; preserveAspect]
u[21]=0;
j=Expand[j/.func[20]]
:[font = text; inactive; preserveAspect]
This is maximized by setting u[20]=0. In general it is clear that j, expressed in terms of values of week t, depends on u[t] through a single term which is proportional to u[t] w[t] and is maximized by equating u[t] to 0 or 1 according as the coefficient of this term is negative or positive. We now write a general program to find u[t]:
:[font = input; preserveAspect]
u[21]=.;
j=q[22];
Do[j=Expand[j/.func[t-1]];
u[t-1]=If[Negative[Coefficient[j,u[t-1] w[t-1]]],0,1],
{t,22,1,-1}]
Table[{t,u[t]},{t,0,21}]//MatrixForm
:[font = text; inactive; preserveAspect; endGroup]
The optimal strategy is to make workers only until the last four weeks and then to make queens only.
:[font = subsubsection; inactive; preserveAspect; startGroup]
4b
:[font = text; inactive; preserveAspect]
From the values in the text, the values of R and m on a daily basis are 0.05 and 0.023, respectively. The discrete-time analogue of Equation 41 with time measured in days is
w(t+1) = 0.05u(t)w(t) + 0.977w(t)
q(t+1) = 0.05[1 - u(t)]w(t) + q(t)
We now repeat the previous calculations with the season ending on day 154:
;[s]
3:0,0;49,1;50,0;328,-1;
2:2,13,9,Times,0,12,0,0,0;1,14,10,Symbol,0,12,0,0,0;
:[font = input; preserveAspect]
Clear[u]
func[t_]:={w[t+1]->.05 u[t] w[t] + .977 w[t],
q[t+1]->.05(1-u[t]) w[t] + q[t]}
:[font = input; preserveAspect]
j=q[154];
Do[j=Expand[j/.func[t-1]];
u[t-1]=If[Negative[Coefficient[j,u[t-1] w[t-1]]],0,1],
{t,184,1,-1}]
Table[{t,u[t]},{t,0,153}]//MatrixForm
:[font = input; preserveAspect; endGroup; endGroup]
The optimal strategy is to make workers only until the last 27 days and then to make queens only.
;[s]
1:0,1;98,-1;
2:0,12,10,Courier,1,12,0,0,0;1,12,9,Times,0,12,0,0,0;
:[font = subsection; inactive; Cclosed; preserveAspect; startGroup]
5
:[font = text; inactive; preserveAspect]
As before we define a function which replaces w[t] and q[t] in terms of their values in the previous week according to the recurrence relation:
:[font = input; preserveAspect]
Clear[u]
func[t_]:={w[t+1]->.35 u[t] w[t] + .84 w[t],
q[t+1]->.35(1-u[t]) w[t] + q[t]}
:[font = text; inactive; preserveAspect]
In week 24 the objective function is j=q[24]:
:[font = input; preserveAspect]
j=q[24];
:[font = text; inactive; preserveAspect]
In week 23 the objective function is q[23]*q[24], which can be expanded in terms of values in week 23:
:[font = input; preserveAspect]
j=Expand[q[23] j/.func[23]]
:[font = text; inactive; preserveAspect]
This is maximized by setting u[23]=0; we then repeat the calculation for week 22, after updating the objective function:
:[font = input; preserveAspect]
u[23]=0;
j=Expand[q[22] j/.func[22]]
:[font = text; inactive; preserveAspect]
It is not possible to find the value of u which maximizes this expression without knowing the relative numbers of queens and workers in week 22, which cannot be determined without knowing the values of u in previous weeks. This circularity can be broken by an iterative procedure. We first assume a trial solution for u[t] and calculate the numbers of queens and workers in different weeks:
:[font = input; preserveAspect]
Do[u[t]=If[t<18,1,0],{t,0,23}]
ww[t_]:=ww[t]=.35 u[t-1] ww[t-1] + .84 ww[t-1];
qq[t_]:=qq[t]=.35(1-u[t-1]) ww[t-1] + qq[t-1];
ww[0]=10;qq[0]=0;
Do[ww[t];qq[t],{t,1,24}]
:[font = text; inactive; preserveAspect]
Clear the trial values for u and determine the optimal values for u assuming the numbers of queens and workers just calculated:
:[font = input; preserveAspect]
Clear[u]
j=q[24];
j=Expand[q[23] j/.func[23]];
u[23]=0;
j=Expand[q[22] j/.func[22]];
jj=j/.{w[22]->ww[22],q[22]->qq[22],u[22]->u}
:[font = text; inactive; preserveAspect]
The value of u which maximizes this function can be seen by plotting jj against u; it is at u=0:
:[font = input; preserveAspect]
Plot[jj,{u,0,1}]
:[font = text; inactive; preserveAspect]
Equate u[22] to zero and go back to week 21:
:[font = input; preserveAspect]
u[22]=0;
j=Expand[q[21] j/.func[21]];
jj=j/.{w[21]->ww[21],q[21]->qq[21],u[21]->u};
Plot[jj,{u,0,1}]
:[font = text; inactive; preserveAspect; endGroup]
The optimal value of u[21] is also zero. Continue backwards in this way. Remember that if j is the objective function at time t+1, the objective function at time t is q[t] j when t>=20 but is j when t<20. At t=17 you will find that there is an internal maximum near u=0.4; its exact value can be found by defining d=D[jj,u] and then solving for d=0 near u=0.4 using FindRoot. For t>17 the optimal strategy is u=0 and for t<17 it is u=1.
Now recalculate the numbers of queens and workers using the above values of u, and then recalculate the optimal values of u. Continue this process until the values of u have converged to 3 decimal places.
:[font = subsection; inactive; Cclosed; preserveAspect; startGroup]
6
:[font = text; inactive; preserveAspect]
Calculate the l[x] and m[x] functions as in Exercise 4.2a
:[font = input; preserveAspect]
temp={.25,.46,.77,.65,.67,.64,.88};
l=Table[Product[temp[[i]],{i,x}],{x,7}];
:[font = input; preserveAspect]
m={1.28,2.28,2.28,2.28,2.28,2.28,2.28};
:[font = text; inactive; preserveAspect]
The age-specific sensitivities from Equations 50 are:
:[font = input; preserveAspect]
Table[Sum[l[[x]] m[[x]],{x,a+1,7}],{a,0,6}]
:[font = input; preserveAspect]
The age-specific sensitivities from Equations 51 are:
;[s]
1:0,1;54,-1;
2:0,12,10,Courier,1,12,0,0,0;1,12,9,Times,0,12,0,0,0;
:[font = input; preserveAspect; endGroup]
l
:[font = subsection; inactive; Cclosed; preserveAspect; startGroup]
7
:[font = text; inactive; preserveAspect]
I first define a matrix of the fitnesses of the four genotypes in the two environments:
:[font = input; preserveAspect]
w={{1.00,0.6,0.785,0.776},{0.58,1.0,0.785,0.815}};
:[font = text; inactive; preserveAspect]
I now calculate the arithmetic and geometric mean fitnesses for the four genotypes. (This calculation can easily be done by hand.)
:[font = input; preserveAspect]
Apply[Plus,w]*.5
:[font = input; preserveAspect]
Apply[Times,w]^.5
:[font = text; inactive; preserveAspect]
I now write a function which will calculate the numbers of the four genotypes next year given a list of their numbers this year. Random[Integer,{1,2}] returns a random integer which is equally likely to be 1 or 2, used to code for wet and dry years.
;[s]
3:0,0;129,1;151,0;250,-1;
2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0;
:[font = input; preserveAspect]
f[x_List]:=(rand=Random[Integer,{1,2}];
100 Table[x[[i]] w[[rand,i]],{i,4}]/x.w[[rand]])
;[s]
2:0,1;1,0;98,-1;
2:1,12,10,Courier,1,12,0,0,0;1,13,10,Times,2,12,0,0,0;
:[font = text; inactive; preserveAspect]
Finally, I use NestList to simulate the history of the population for 1000 years. You may extract appropriate elements from the list res. The idea of using a random number generator to simulatre a stochastic process is a very useful tool.
;[s]
5:0,0;15,1;23,0;133,1;136,0;239,-1;
2:3,13,9,Times,0,12,0,0,0;2,13,9,Times,1,12,0,0,0;
:[font = input; preserveAspect]
res=NestList[f,{25,25,25,25},1000];
:[font = text; inactive; preserveAspect; endGroup]
:[font = subsection; inactive; Cclosed; preserveAspect; startGroup]
8
:[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup]
8a
:[font = text; inactive; preserveAspect]
I define the log geometric mean fitness from Equation 54 and then find the germination rate that maximizes it.
:[font = input; preserveAspect]
f=(1-p) Log[(1-g) (1-d)] + p Log[(1-g) (1-d) + g s];
:[font = input; preserveAspect; endGroup]
Solve[D[f,g]==0,g]
:[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup]
8b
:[font = text; inactive; preserveAspect]
I define the optimal germination rate, g, given s, p, and d; I write a procedure to find the value of s by equating the log geometric mean fitness to zero; and I use this procedure to evaluate g and s for different values of d and p.
;[s]
19:0,0;39,1;40,0;48,1;49,0;51,1;52,0;58,1;59,0;102,1;103,0;193,1;194,0;199,1;200,0;225,1;226,0;231,1;232,0;234,-1;
2:10,13,9,Times,0,12,0,0,0;9,13,9,Times,1,12,0,0,0;
:[font = input; preserveAspect]
g=(s p+d-1)/(s+d-1);
:[font = input; preserveAspect]
func[x_,y_]:=(d=x;p=y;
{d,p,g,s}/.FindRoot[f==0,{s,1/p}])
:[font = input; preserveAspect]
dlist={.01,.1,.5};plist={.1,.25,.75};
:[font = input; preserveAspect; endGroup; endGroup; endGroup]
Table[func[dlist[[i]],plist[[j]]],{i,3},{j,3}]//TableForm
^*)