CHAPTER 3

1

1a

Use pencil and paper.

1b

The procedure sol[a,b,t] calculates the numerical solution of the differential equation from time 0 to t, starting with n1=a, n2=b. ParametricPlot gives a parametric plot (phase plane diagram) of the solution. You should calculate some more examples and plot them on the same graph. You will need to add the arrows manually.

sol[a_,b_,time_]:=
	NDSolve[{n1'[t]==.5(1-(n1[t]+3n2[t]/2)/1000)n1[t],
			n2'[t]==.5(1-(n2[t]+3n1[t]/2)/1000)n2[t],
			n1[0]==a,n2[0]==b},
			{n1,n2},{t,time}]

example=sol[200,600,20];

ParametricPlot[Evaluate[{n1[t],n2[t]}/.example],
					{t,0,20},
				Axes->Automatic,
				PlotRange->{{0,1500},{0,1500}}]

2

We first define the two functions representing the differential equations for the Lotka-Volterra competition model. It is convenient to write them as indexed functions:

f[1][n1_,n2_]:=r1(1-(n1+g12 n2)/k1)n1
f[2][n1_,n2_]:=r2(1-(n2+g21 n1)/k2)n2

In principle we could find the four sets of equilibria using "Solve" but in practice this takes an inordinate time. However, it is clear that they are the equilibria (i) with both species absent, (ii) with species 1 present and 2 absent, (iii) with species 2 present and 1 absent, and (iv) with both species present. The first three are easy to find. The fourth can be found using "Solve" on the pair of linear equations:

equil=Solve[{1-(n[1]+g12 n[2])/k1 == 0,
		1-(n[2]+g21 n[1])/k2 == 0},{n[1],n[2]}]

(You will see shortly why I am now using indexed variables for n.) The solution has two braces round it, allowing for the possibility of several solutions. We shall use Flatten to remove one of the braces, and at the same time try to simplify the result:

equil=Simplify[Flatten[equil]]

The matrix of partial derivatives of the functions f[i] with respect to the variables n[j] can be calculated as follows:

jac=Table[D[f[i][n[1],n[2]],n[j]],{i,2},{j,2}]

It remains to evaluate this matrix at an appropriate equilibrium:

jac/.{n[1]->k1,n[2]->0}

If necessary use Simplify[] on the result. To put in numerical values for the parameters you may use a second replacement rule:

ans=jac/.equil/.{r1->.5,r2->.5,g12->2/3,
							g21->2/3,k1->1000,k2->1000}

To find the eigenvalues and eigenvectors use Eigensystem[]:

Eigensystem[ans]

3

Mathematica is useful to simplify the algebra. I illustrate the most complicated case, that of the internal equilibrium, using jac and equil from the previous exercise. From there on, use pencil and paper.

mat=jac/.equil;

d=Simplify[Det[mat]]

t=Simplify[Sum[mat[[i,i]],{i,2}]]

4

See Exercise 2.

5

5a

This is also similar to Exercise 2. It is first necessary to clear the old definition of f and then to redefine the functions.

Clear[f]

f[1][n1_,n2_]:=r1 n1(1-n1/k1) - a1 n1 n2/(1+b n1)
f[2][n1_,n2_]:=-r2 n2 + a2 n1 n2/(1+b n1)

r1=r2=1;a1=a2=.01;b=.005;

equil={n[1]->200,n[2]->200(1-200/k1)};

jac=Table[D[f[i][n[1],n[2]],n[j]],{i,2},{j,2}];

eigen=Eigenvalues[jac/.equil]

Plot[Max[Re[eigen]],{k1,201,1000}]

5b

See Exercise 1b.

6

Start, as before, by evaluating the Jacobian matrix at the two equilibria. (I write nn for N and dd for D since capital letters are reserved for built-in functions.)

Clear[f]

f[1][x_,y_]:=d(nn - x)+(dd - d)y - beta x y
f[2][x_,y_]:=beta x y -(dd + gamma)y

jac=Table[D[f[i][n[1],n[2]],n[j]],{i,2},{j,2}];

equil={{n[1]->nn,n[2]->0},
			{n[1]->(dd+gamma)/beta,
		n[2]->(nn-(dd+gamma)/beta)/(1+gamma/d)}};

Simplify[jac/.equil[[1]]]//TableForm

The determinant is positive when (dd + gamma) > beta nn, in which case the trace is negative.

Simplify[jac/.equil[[2]]]//TableForm

The determinant is positive when (dd + gamma) > beta nn, in which case the trace is negative.

7

We may use the model of Exercise 6 if we take nn=1, d=dd=1/70, gamma=50, beta=500; I write these in decimal form to force numerical evaluation.

nn=1.;
d=dd=1/70.;
gamma=50.;
beta=500.;

equil[[2]]

Eigenvalues[jac/.%]

N[2Pi/Im[%[[1]]]]

8

8a

Use pencil and paper.

8b

Stability analysis is similar for discrete-time models, but depends on the eigenvalues being less than 1 in absolute value rather than on their real part being negative.

Clear[f]

f[1][n1_,n2_]:=r n1 Exp[-n2]
f[2][n1_,n2_]:=n1(1-Exp[-n2])

equil={n[1]->r Log[r]/(r-1),n[2]->Log[r]};

eigen=Eigenvalues[Table[D[f[i][n[1],n[2]],n[j]],
		{i,2},{j,2}]/.equil];

Plot[Max[Abs[eigen]],{r,1.01,5}]

Plot[Im[eigen[[1]]],{r,1.01,5}]

8c

The following program will iterate the equations using NestList acting on a two-dimensional function. Change the value of R, or the starting values, to obtain different simulations.

g[{n1_,n2_}]:={f[1][n1,n2],f[2][n1,n2]}/.r->2.0

res=NestList[g,{1.0,1.0},50];

{h,p}=Transpose[res];

host=ListPlot[h,PlotJoined->True,
		PlotRange->All];

parasitoid=ListPlot[p,PlotJoined->True,
		PlotRange->All];

Show[host,parasitoid]

9

See previous exercise.