Lab 3:  Writing simulations in Mathematica

Mathematica can be used as a programming language in which you can run
simulations.  However, there are a few tricks you need to learn about writing
simulations in Mathematica.

Say that you have the discrete generation recursion equation:

	x[t] = a x[t-1]			(exponential growth)
	
x[t] is a function of one parameter (a) and x[t-1], so you could simply write
the function:

	x[t_,a_] = a*x[t-1,a] 
	
THIS WON'T WORK THOUGH!!

It sends Mathematica into an infinite loop because it immediately tries to find
x at time t as a function of x at time t-1, but it doesn't know that either, so it
tries to find x at time t-2, etc.

It would be OK to use

	x[t_,a_] := a*x[t-1,a] 
	
The := tells Mathematica not to calculate anything until you tell it to.

The second step would be to set the initial conditions.  Say:

	x[0,a_] = 10

Remember, the underline (_) tells Mathematica to make the definition for all
values of a.

You could then ask Mathematica what x is at any time t for particular values of a:

	x[10, 1.1]
	x[100,1.1]

THIS TURNS OUT TO BE REALLY SLOW!!

The problem is that Mathematica does not remember previous values that it
has calculated. For every time t, Mathematica goes back to the recursion to
figure out x all the way up from 0 to t.

To make Mathematica remember, we use a trick:

	x[t_,a_] := x[t,a] = a*x[t-1,a] 
	
How does this help?  The first time Mathematica sees x[10, 1.1] it goes to
the recursion to calculate it, but the second part (x[t,a] =) then sets x[10, 1.1]
forever more to the value on the far right.  Every time x[10, 1.1] is called in
the future, Mathematica now remembers it!

Example

Clear[x,a,t]
x[t_,a_] := x[t,a] = a*x[t-1,a]
x[0,a_] := 1

(I always Clear all the parameters and variables first.)

x[10,1.1]

x[100,1.1]

ListPlot[Table[x[t,1.1],{t,0,100}]]

Three hints about graphing:

(1)  You can't use Plot, which only works for continuous functions
(we only know x at the discrete time points of each generation).
You must use ListPlot, which plots a list of points.  In this case
the list of points is generated by a table giving x for generations 0 to 100.
	
(2)  If you want to see the points connected, type before the last ]:

	PlotJoined->True
	
(3)  If you want to change the range of the plot, type before the last ]:

	PlotRange->{{0,110},{0,10000}}

ListPlot[Table[x[t,1.1],{t,0,100}],
	PlotJoined->True,
	PlotRange->{{0,110},{0,10000}}]

Question for the lab

The melanic form of the British peppered moth was observed to increase from
a frequency of 0.00001 to 0.8 in 50 generations.

The melanic form is dominant and we will call the fitnesses of the three genotypes:

	WAA = 1,
	WAa = 1+s
	Waa = 1+s
	
The recursion equation for diploid selection developed in class is then:

	p[t+1] = ((1+s)*p[t]^2+(1+s)*p[t]*(1-p[t])) / ((1+s)*p[t]^2+(1+s)*2*p[t]*(1-p[t])+(1-p[t])^2)

Rewrite this equation in the form of a simulation in Mathematica.

Use Mathematica to simulate this equation for various values of s to estimate how
strong selection must have been to take p from 0.00001 to 0.8 in 50 generations.

Solution for the lab

Clear[p,t,s]
p[t_,s_]:=p[t,s]=(
	((1+s)*p[t-1,s]^2+(1+s)*p[t-1,s]*(1-p[t-1,s]))/
	((1+s)*p[t-1,s]^2+(1+s)*2*p[t-1,s]*(1-p[t-1,s])+(1-p[t-1,s])^2)
	)

p[0,s_]:=0.00001

ListPlot[Table[p[t,0.3],{t,0,100}],
	PlotJoined->True]