Lab 6:  Solving differential equations

We have already learned how to use Mathematica to iterate a recursion equation, but how do we use Mathematica to solve a differential equation?

In this lab, we will investigate the continuous-time model of population growth in the presence of habitat degradation. Habitat loss can cause the growth rate of a population to decline over time.  Assuming that this decline is linear and that the initial growth rate is r0 at t = 0, we have a modified model of exponential growth:

		dn/dt = (r0 - delta t) n

We'll learn how to use Mathematica to:

	(1)  Graph the solution to a differential equation
	(2)  Find the general solution of it.

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

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

The way that Mathematica numerically solves differential equations is a trade secret, but it is said to involve 500 pages of code.  Roughly, Mathematica determines what n must be as a function of time to have the right slope, dn/dt, at each point in time.

First, you need to tell Mathematica what range of time you want to focus on (because it can't use the above procedure for all time).  Say we want 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 habitat destruction model, we have two equations we want to hold true:

	D[n[t],t] == (r0-delta*t)*n[t]
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] == (r0 - delta*t)*n[t], n[0] == n0}, {n[t]}, {t, tmin, tmax}]

Mathematica 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[habitat, n]
habitat[r0_, delta_, n0_, tmin_, tmax_] := habitat[r0, delta, n0, tmin, tmax] = 
   NDSolve[{D[n[t], t] == (r0 - delta*t)*n[t], n[0] == n0}, {n[t]}, {t, tmin, tmax}]

[NOTE:  habitat[r0_,delta_,n0_,tmin_,tmax_] := 
	habitat[r0,delta,n0,tmin,tmax] =
notation ensures that Mathematica remembers the numerical solution that it finds.]

When we call habitat 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).

habitat[0.1, 0.01, 1000, 0, 100]

We can then plot n[t] by making a substitution of the approximate
function habitat for n[t].

For example, with r0=0.1, delta =0.01, n0=1000, tmin=0, and tmax=100:

plot1 = Plot[Evaluate[n[t] /. habitat[0.1, 0.01, 1000, 0, 100]], {t, 0, 100}, 
  PlotRange -> {{0, 100}, All}]

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

Question # 1

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

			  -0.5, -0.1,  0.1, 0.5.

How is the shape of the plot similar/different for positive and  negative values of r0?

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 DSolve solves the differential equation for all t, not just t between tmin and tmax.

The DSolve command we would use is therefore:

DSolve[{D[n[t], t] == (r0 - delta*t)*n[t]}, {n[t]}, t]

Question # 2

Enter the above equation into Mathematica and see what it gives you.  Solve for the constant of integration (C[1]) to obtain the general solution.

Assuming that delta is positive, consider what will happen to the population size as time goes to infinity (as t gets larger and larger).  Say what the ultiimate fate of the population must be.

Question # 3

Compare the numerical solution to the differential equation provided by NDSolve to a simulation of the discrete-time recursion equation:

		n[t+1] = (r0 - delta t) n[t]

Write a discrete-time simulation using code from previous weeks and plot the results.  Compare the plots to that seen from NDSolve (above) with r0=0.1, delta =0.01, n0=1000, tmin=0, and tmax=100.  You can show two plots together using the command Show[plot1,plot2]

Are the results for the continuous- and discrete-time models similar?  Would you expect this to be true if r0 were larger?  Explore how changing r0 affects the similarity of the models.