(*^ ::[ 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 ^*)