Lab 9:  Solving differential equations

In this lab, we'll compare the dynamics of the one- and two-species models with density-dependent growth.  We'll use the continuous time models and learn how to numerically and analytically solve them with Mathematica.

The one-species model with density-dependent growth is the logistic:

		dn/dt = r n (1-n/k)

The two-species model with density-dependent growth and competition is:

		dn1/dt = r1 n1 (1-(n1+a12*n2)/k1)
		dn2/dt = r2 n2 (1-(n2+a21*n1)/k2)

There are two things we would like to be able to do with such differential equations:

	(1)  Graph it.
	(2)  Find the general solution of it.

Mathematica can help us do both of these.

Numerically solving differential equations
(See Mathematica book, page 696)

First determine the parameters needed to solve the equation numerically:
	 r, 
	 k, 
	n0 (the initial population size at time 0)

Mathematica numerically solves differential equations by a series of guesses.  It starts with a rough guess of the general solution (n as a function of time) and sees if it has the right derivative (dn/dt) at various points in time.  If not, it changes the rough solution to bring its slope closer to dn/dt.

To do this, though, you need to tell Mathematica what range of time you want to focus on (because it can't use the above procedure of guessing for all times).  Say we care to know the solution for n between 

	tmin (the first time point that you want to know about), 
	tmax (the last time point that you want to know about).
	
To numerically solve a differential equation, we use NDSolve:

	NDSolve[{equations you want to hold true},
		{variables you want to solve for}, {t,tmin,tmax}]

In the case of the logistic model, we have two equations we want to hold true:

	D[n[t],t] == r*n[t]*(1-n[t]/k)
and	
	n[0] == n0
	
and we want to solve for one variable (n[t]).

[NOTE:  The == is used to denote an equation that we want to be true.  We're not setting the right hand side to the left hand side.]

The NDSolve command we would use is:

NDSolve[{D[n[t],t] == r*n[t]*(1-n[t]/k),
			n[0] == n0},
			{n[t]}, {t,tmin,tmax}]

Mathematic won't solve this yet, since it doesn't know what the parameters are.  What I do is define a new function of all of the parameters and set this function to the above NDSolve:

Clear[logistic,n]
logistic[r_,k_,n0_,tmin_,tmax_] := 
	logistic[r,k,n0,tmin,tmax] =
	NDSolve[{D[n[t],t] == r*n[t]*(1-n[t]/k),
			n[0] == n0},
			{n[t]}, {t,tmin,tmax}]

[NOTE:  logistic[r_,k_,n0_,tmin_,tmax_] := 
	logistic[r,k,n0,tmin,tmax] =
This notation ensures that Mathematica remembers the numerical solution that it finds.]

When we call logistic with all the parameters, Mathematica will return its best numerical guess for the function, n[t], between tmin and tmax (Mathematica calls this an interpolating function).

logistic[0.1,1000,100,0,100]

We can then plot n[t] by substituting this function for n[t].

For example, with r=0.2, K=1000, n0=100, tmin=0, and tmax=100:

Plot[Evaluate[n[t] /. logistic[0.2,1000,100,0,100]],{t,0,100},
	PlotRange->{{0,100},{0,1100}}]

[NOTE:  "Evaluate" tells Mathematica to numerically solve the equation now.  Otherwise it just keeps it in NDSolve notation.]

Question # 1

Plot n[t] as above for the following values of r (try others if you wish):

			  -0.1, -1,  0.1, 1.

In which cases is n=K approached and in which cases is n=0 approached?

Do these results match the local stability analysis you performed in a previous homework?

Question # 2

Write the Mathematica code to find the numerical solution to the two-species competition model.  Hint:  write a function called competition that is a function of r1, k1, a12, n10, r2, k2, a21, n20, tmin, tmax.  This competition function should contain the NDSolve command, within which are the two differential equations describing the rate of change of species 1 and 2 (see above).  Don't forget to specify two initial conditions (n10 and n20) and to ask Mathematica to solve for both n1[t] and n2[t].

Clear[competition,n1,n2]
competition[r1_,k1_,a12_,n10_,r2_,k2_,a21_,n20_,tmin_,tmax_] := 
	competition[r1,k1,a12,n10,r2,k2,a21,n20,tmin,tmax] =
	NDSolve[{D[n1[t],t] == r1*n1[t]*(1-(n1[t]+a12*n2[t])/k1),
		D[n2[t],t] == r2*n2[t]*(1-(n2[t]+a21*n1[t])/k2),
			n1[0] == n10, n2[0] == n20},
			{n1[t],n2[t]}, {t,tmin,tmax}]

Plot n1[t] and compare this to the graph of the logistic.  Use the following parameter values:

(a)  r1=0.1, k1=1000, n10=100, a12 = 0, r2=0.1, k2=1000, n20=100, a21 = 0, tmin=0, tmax =100.
(b)  r1=0.1, k1=1000, n10=100, a12 = 0.2, r2=0.1, k2=1000, n20=100, a21 = 0, tmin=0, tmax =100.
(c)  r1=0.1, k1=1000, n10=100, a12 = 2.0, r2=0.1, k2=1000, n20=100, a21 = 0, tmin=0, tmax =100.
(d)  At least one more parameter set of your choosing.

For example, I used the plot function:

Plot[Evaluate[n1[t] /. competition[0.1,1000,0,100,0.1,1000,0,100,0,100]],{t,0,100},
	PlotRange->{{0,100},{0,1100}}]

General solution of a differential equation
(See Mathematica book, page 637)

Mathematica knows the general analytical solution to some (but by no means all) differential equations.

This is helpful in many cases where solving the differential equation is laborious or difficult.

The procedure is very similar, except we use DSolve rather than NDSolve (the N stands for "numerical").

The main difference is that DSolve solves the differential equation generally, and does not like to use initial starting conditions (like n[0] == n0).  Also, DSolve solves the differential equation for all t, not just t between tmin and tmax.

The command we would use to solve the logistic one-species logistic model is:

DSolve[{D[n[t],t] == r*n[t]*(1-n[t]/k)},
			{n[t]}, t]

Question # 3

Enter the above equation into Mathematica and see what it gives you.

What should the constant of integration (C[1]) be if n[0] = n0?