CHAPTER 5

1

The model with b* = 45 can be simulated for 100 generations starting with both types at densities of 0.5 as follows:

f[{x_,y_}]:=(P0=.02 Exp[-.01(100 x + bstar y)];
			{100 P0 x,(bstar P0 + 0.5)y})

bstar=45;

res=NestList[f,{0.5,0.5},100]

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.)

2

The model with b* = 70 can be simulated for 100 generations starting with both types at densities of 0.5 as follows:

Clear[bstar]

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})

bstar=70;

res=NestList[f,{0.5,0.5},100]

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.

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:

Clear[bstar]

g[y_]=.02 Exp[-.01 bstar y] bstar + .5 Exp[-y];

bstar=70;

equil=FindRoot[g[y]==1,{y,.5}]

Now complete the calculations:

P1=.5 Exp[-y]/.equil

bstar/(1-P1)

Repeat these calculations with b* = 77 and 85.

3

This exercise is best done with pencil and paper.

4

4a

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:

func[t_]:={w[t+1]->.35 u[t] w[t] + .84 w[t],
q[t+1]->.35(1-u[t]) w[t] + q[t]}

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:

j=q[22];
j=Expand[j/.func[21]]

j is maximized by setting u[21]=0. We do this and evaluate j in terms of values in week 20:

u[21]=0;
j=Expand[j/.func[20]]

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]:

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

The optimal strategy is to make workers only until the last four weeks and then to make queens only.

4b

From the values in the text, the values of R and μ 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:

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]}

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

The optimal strategy is to make workers only until the last 27 days and then to make queens only.

5

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:

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]}

In week 24 the objective function is j=q[24]:

j=q[24];

In week 23 the objective function is q[23]*q[24], which can be expanded in terms of values in week 23:

j=Expand[q[23] j/.func[23]]

This is maximized by setting u[23]=0; we then repeat the calculation for week 22, after updating the objective function:

u[23]=0;
j=Expand[q[22] j/.func[22]]

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:

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}]

Clear the trial values for u and determine the optimal values for u assuming the numbers of queens and workers just calculated:

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}

The value of u which maximizes this function can be seen by plotting jj against u; it is at u=0:

Plot[jj,{u,0,1}]

Equate u[22] to zero and go back to week 21:

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}]

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.

6

Calculate the l[x] and m[x] functions as in Exercise 4.2a

temp={.25,.46,.77,.65,.67,.64,.88};
l=Table[Product[temp[[i]],{i,x}],{x,7}];

m={1.28,2.28,2.28,2.28,2.28,2.28,2.28};

The age-specific sensitivities from Equations 50 are:

Table[Sum[l[[x]] m[[x]],{x,a+1,7}],{a,0,6}]

The age-specific sensitivities from Equations 51 are:

l

7

I first define a matrix of the fitnesses of the four genotypes in the two environments:

w={{1.00,0.6,0.785,0.776},{0.58,1.0,0.785,0.815}};

I now calculate the arithmetic and geometric mean fitnesses for the four genotypes. (This calculation can easily be done by hand.)

Apply[Plus,w]*.5

Apply[Times,w]^.5

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.

f[x_List]:=(rand=Random[Integer,{1,2}];
							100 Table[x[[i]] w[[rand,i]],{i,4}]/x.w[[rand]])

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.

res=NestList[f,{25,25,25,25},1000];

8

8a

I define the log geometric mean fitness from Equation 54 and then find the germination rate that maximizes it.

f=(1-p) Log[(1-g) (1-d)] + p Log[(1-g) (1-d) + g s];

Solve[D[f,g]==0,g]

8b

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.

g=(s p+d-1)/(s+d-1);

func[x_,y_]:=(d=x;p=y;
			{d,p,g,s}/.FindRoot[f==0,{s,1/p}])

dlist={.01,.1,.5};plist={.1,.25,.75};

Table[func[dlist[[i]],plist[[j]]],{i,3},{j,3}]//TableForm