Mathematica Appendix

When the authors teach this course, we use a variety of user friendly specialty software such as spreadsheets and statistical packages.  This approach enables the students to spend the bulk of their time on modeling as opposed to programming and debugging.  This course could be taught entirely using a programming language such as C/C++, BASIC or TRUEBASIC, Pascal, or FORTRAN.  In the process the students would become expert programmers of mathematical algorithms, but will have much less to devote to the modeling and interpretation process.  Indeed writing a package to compute ANOVA tables alone would be a large project, when our intent is to introduce students to the ideas, not the details of the computations.  An option between many specialized programs or one programming language is a mathematical programming language such as Mathematica or Maple.  Every model in this book could be done using Mathematica, though most any model is more easily built using something else.  On the other hand, the programs which are easier to use are limited in the different types of models they can handle.  A spreadsheet program would be a good second choice for the number of types of models it can handle.

This appendix, then, is intended to illustrate how to build the models in the text using Mathematica.  While comments will be made throughout the examples, this appendix is not intended to be a Mathematica manual or tutorial.  Rather, it assumes some familiarity on the users part with Mathematica, and merely demonstrates how to put the pieces together, how to use special features, and how to avoid common problems.  To get familiar with Mathematica, the reader is directed to the book The Beginner's Guide to Mathematica by Gray and Glynn, the Tour of Mathematica in Wolfram's Mathematica book (which can be found in the electronic help in version 3.0), or the Introduction to Mathematica notebook included with the text Calculus&Mathematica by Davis, Porta, and Uhl.

In the examples in this appendix, we have tried to use Mathematica as a mathematical tool rather than a programming language whenever possible.  Thus we have tried perform calculations in a manner and order reminiscint of standard mathematics.  Reader skilled in programming will undoubtedly see ways to optimize code or use advanced programming capabilities of Mathematica.  It is our opinion that while fancy code may be more efficient and elegant, for beginning users it can obscure the mathematical and modeling ideas.

Chapter 1

Mathematica preliminaries

When getting Mathematica to solve recurrsion relations, the concept of caching is very important.  For example, consider the Fibonacci recurrence relation.  Naively this can be set up as:

Clear[x,n];
x[0]=1;
x[1]=1;
x[n_]:=x[n-1]+x[n-2];

Now suppose we want the Fibonnaci number when n=15.  We add the //Timing command to keep track of the amount of time require to find x[15].

x[15]//Timing

Notice that if we run it again, it takes about the same amount of time.

x[15]//Timing

This time we again define a recursive Fibonacci function, but this time we add a piece of code that seems redundant.  In particular, note the x[n]:=x[n]=... construction.  This is called caching and everytime it computes an x[n] value it stores it in an internal lookup table.  Notice that when x[15] is requested, x[14] and x[13] are needed.  With the original code, these are computed regardless of whether they have been computed before.  With caching, if  x[14] and x[13] are already computed they are merely looked up resulting in the great time savings observed below.

Clear[x,n];
x[0]=1;
x[1]=1;
x[n_]:=x[n]=x[n-1]+x[n-2];

x[15]//Timing

If we request x[15] again, it only needs to be looked up so it takes even less time.

x[15]//Timing

Warning!  Since the values of x are stored, it is essential to use the Clear command every time the recursive function is changed.  While clearing variables is always good practice, often it can be skipped without any problems.  This is not the case when caching is being used.  Anytime unexpected results are being obtained, clearing the variables should be a considered as a first debugging device.

Sandhill Crane

In this section we replicate the sandhill crane spreadsheet analysis.

The simple model

This code should be reminiscent of the Fibonnaci sequence discussed above.

Clear[pop,n];
pop[0]=100;
pop[n_]:=pop[n]=1.0194pop[n-1];

The Table command can be used to create a list of data.  Here the data consists of ordered pairs.  We give the list the name cranes so that it can be used later.

cranes=Table[{n,pop[n]},{n,0,14}]

The suffix //TableForm can be used to put the data in a nice table format.

cranes//TableForm

ListPlot[cranes,PlotJoined->True,AxesLabel->{"Years", "Population"}];

A model with three growth rates

We replicate our analysis for best, medium, and worst conditions, and show all three on a single graph.  One can see each part of the simple analysis replicated three times with different names.  To show all graphs on the same coordinate axes, give each one a name and then use the Show command.  Here we have made the graph of each model a different thickness so they can be identified on the combined graph.

Clear[popBest, popMedium, popWorst,n];
rBest=0.0194;
rMedium=-0.0324;
rWorst=-0.0382;

popBest[0]=100;
popMedium[0]=100;
popWorst[0]=100;

popBest[n_]:=popBest[n]=(1+rBest) popBest[n-1];
popMedium[n_]:=popMedium[n]=(1+rMedium) popMedium[n-1];
popWorst[n_]:=popWorst[n]=(1+rWorst) popWorst[n-1];

cranesBest=Table[{n,popBest[n]},{n,0,14}];
cranesMedium=Table[{n,popMedium[n]},{n,0,14}];
cranesWorst=Table[{n,popWorst[n]},{n,0,14}];

pBest=ListPlot[cranesBest,PlotJoined->True,AxesLabel->{"Years", "Population"}, PlotStyle->Thickness[0.005]];
pMedium=ListPlot[cranesMedium,PlotJoined->True,AxesLabel->{"Years", "Population"},PlotStyle->Thickness[0.01]];
pWorst=ListPlot[cranesWorst,PlotJoined->True,AxesLabel->{"Years", "Population"},PlotStyle->Thickness[0.015]];

Show[pBest,pMedium,pWorst];

Hacking Chicks

To modify a model, we can either make the change and rerun it, or, if we are interested in keeping the results cut and paste the old code and modify it appropriately.  Here we hack 5 chicks to the medium and worst cases.

Clear[popBest, popMedium, popWorst,n];
chicks=5;
yearsToRun=120;
rMedium=-0.0324;
rWorst=-0.0382;

popMedium[0]=100;
popWorst[0]=100;

popMedium[n_]:=popMedium[n]=(1+rMedium) popMedium[n-1]+chicks;
popWorst[n_]:=popWorst[n]=(1+rWorst) popWorst[n-1]+chicks;

cranesMedium=Table[{n,popMedium[n]},{n,0,yearsToRun}];
cranesWorst=Table[{n,popWorst[n]},{n,0,yearsToRun}];

pMedium=ListPlot[cranesMedium,PlotJoined->True,AxesLabel->{"Years", "Population"},PlotStyle->Thickness[0.01]];
pWorst=ListPlot[cranesWorst,PlotJoined->True,AxesLabel->{"Years", "Population"},PlotStyle->Thickness[0.015]];

Show[pMedium,pWorst];

Annual Plants

In this section, we use Mathematica to numerically analyze the annual plant model.  We first build a recursive model, then use the symbolic algebraic capabilities to find a closed form solution.

Define parameters.

Clear[p,alpha1, beta1, gamma1, sigma1, n];
alpha1=.5;
beta1=.25;
gamma1=2;
sigma1=.8;

Define the recursion relation

p[n_]:=p[n]=alpha1*sigma1*gamma1*p[n-1]+
beta1*sigma1*(1-alpha1)*sigma1*gamma1*p[n-2];

p[-1]=100;
p[0]=95;

TableForm[Table[{n, p[n]}, {n, 0, 10}]]

\(\(ListPlot[Table[{n, p[n]}, {n, 0, 10}], PlotJoined -> True, 
	
      AxesLabel -> \ {"Years", \ "Number of Plants"}]; \)\)

To find the closed form solution, we need to solve the characteristic equation.

Solve[
      (l^2) - alpha1 sigma1 gamma1 l - 
          beta1 sigma1 (1 - alpha1) sigma1 gamma1 == 0, l]

{{l -> -0.165685424949238019527`}, {
        l -> 0.965685424949238019505`}}

Use the initial conditions to find the constants.

Solve[{95 == C1 + C2, 
        100 == -(C1/0.165685424949238019527)` + 
            (C2/0.965685424949238019505)`}, {C1, C2}]

To get these numbers in useable form, we will use the Copy Output From Above menu item and then cut and paste the values into the desired equation form.

{{C1 -> -0.229707730091957898379`, 
        C2 -> 95.2297077300919578984`}}

pp[n_] = -0.229707730091957898379`*(-0.165685424949238019527`)^n + 
  95.2297077300919578984`*(0.965685424949238019505`)^n;

The following table compares the output from the recursive model and the closed form model.

TableForm[Table[{n, p[n],  pp[n]}, {n, 0, 10}]]