Biophysics 101, Fall 2002
Problem Set 1, part 3

Warm-up Notebook

warmup_answers

This notebook is for beginners that have done the Mathematica tutorial. Most of it is no harder than a coloring book  It starts with cell manipulation and formatting, and it works its way up to discrete recursive functions that operate on Mathematica tables and arrays. Many of us are new to Mathematica, so if you discover anything useful, please include them in your finished notebook.

1:  Saving notebooks

Your first challenge is to save a copy of this notebook, and you get to name it after you.  On the menu bar at the top of the screen its "File-> Save As.... ". If your name is Jack Flapjack then name it Jack_Flapjack.nb. If your name is not Jack Flapjack then name it after your own name instead. Put your answers in your self-named copy and turn it in to your TF.

2: Inserting and deleting Mathematica cells

If you click in the thin horizontal space between cells, a thin black line appears. Hit the enter key, and voila, a new cell. Insert a cell between the second and third cells in this notebook, and type your name in it. If you are feeling stylish, format your name text with the many options available under the format menu. On the menu, check out "Format->Show Toolbar" Though short on features, the formatting toolbar may be worth it's screen real estate. To delete a cell, click to highlight it's blue bracket on the right, and press the delete key. To copy/paste a cell - highlight the blue bracket of the cell you want to copy, <right-click>->copy,then  make the thin horizontal black line appear where you want it to go, and <right-click>->paste.

3: Mathematica cell types

Notice the subtle differences in the cell bracket shapes. They tell you what type of cell it is. This type, with the double line at the top, means its a non-evaluatable text cell. If you press <shift-enter> while the cursor is in this type of cell, Mathematica will beep at you. 
(2 pts) Try it on the cell below. Now to make it evaluatable, highlight the cell bracket, and then its "Cell->Cell Properties->Cell Evaluatable"

2+2

The default cell type, for example the one you just typed your name in,  is called "Input Form", and its cell bracket has a small triangle at the top. Change the property of the cell containing your name to non-evaluatable, and convert its display type to Text thru the menu via "Cell->Display As->Text".

This is an input type cell in text display mode. If you press <shift-enter> in this cell, you will see what this sentence means to the Mathematica kernel (not much).

This is a traditional form cell .   It has the coolest bracket shape,  
        but it' s a relic of old Mathematica versions that' 
          s  still  compatable with version  4.2.

4: Cell grouping

Notice how most of these cells are nested, but the two cells in this section are not. To group them, first go on the menu to "Cell->Cell Grouping" and make sure "Manual Grouping" is checked. Then select both brackets by dragging the mouse along the right margin, and its "Cell->Cell Grouping->Group Cells". Double click on the outermost bracket to make the nest expand and collapse.

5: Named lists

In the next cell there are some named lists and a name based list operation. (2 pts)  Evaluate the expressions by placing the cursor to the right of the last line, and <shift-enter>.  To get the total oil consumption for all five years, replace the "*" with a "." and evaluate again.

census  =  {7,  12,  11,  17,  28}
perCapitaOil  =  {300,  445,  425,  625,  
        850}
census  .  perCapitaOil

6: Calculated lists using the Table method

The first example creates a length b array of a. Evaluate the expression.

a = 3
b = 4
Table[a,  {b}]

Evaluate these expressions.

Prime[5]
Table[Prime[n],  {n,  1,  10}]
Random[Real,  {0, 1}]

Table[
      Random[Real,  {1,  50}],  {20}]

(2 pts) Now insert a cell below this one that creates a list of 20 random numbers. Choose your own range (in the above case the range is between 0 and 1), and your number type choices are (Integer, Real, or Complex).

Evaluate these different ways to count to ten.

Table[x,  {x,  10}]
Table[x,  {x,  1,  10}]
Table[x/2,  {x, 2, 20, 2}]
Range[10]

Inset a formula in a cell at the bottom of this group that counts backwards from ten.

Table[x,  {x,  10,  1,  -1}]

7: List generation using functions with the Array method

This example introduces functions and arrays. Notice the underscore (x_) after the variable definition in the function. Arg 1 in the Array specification is the creation method that takes an integer, arg 2 is the array length, arg 3 is the starting point. (2 pts) Evaluate these expressions.

f[x_]  :=  x^2;
qT =  
      Array[f,  20,  -12]
qT
ListPlot[qT];
ListPlot[qT, PlotJoined -> True]

Note 1) Notice how the arrays are printed out when you evaluate the cell. This can be annoying when you want to look at a complicated graph but don't want to see all of it's underlying points printed out. Put a ";" at the end of the expressions that you don't want to display the results of.

Note 2) It's considered good Mathematica practice to begin your custom variable names with a lower case letter, which will prevent name collisions with the Mathematicas names, all of which begin with uppercase letters.

Note 3) Don't forget to put the colon in front of the equals sign ":=" when you define a function. This can be very frustrating, because leaving it off won't cause an error. So if your function is giving you unexpected results, don't forget to re-read this note!

8: Using lists to graph a built in function

You can compose your own custom functions with Mathematica built in functions. Remember N from the tutorial? - evaluate these expressions.

k  =  3;
f[x_]  = k*Tanh[x];
f[3];
N[f[3]]

(2 pts) Make a plot of Tanh[x], with the y values between -0.9 and 0.9. Bonus: Make it exactly between -0.9 and 0.9 (hint: look in Help->Master Index->Built-in Functions->Mathematical Functions).

Its possible to do this using functions with plot or using a lists with listplot.  In order to get the y range exactly between -0.9 and 0.9, use the ArcTanh function.

Plot[ 
      Tanh[x],  {x,  ArcTanh[-0.9],  
        ArcTanh[0.9]} ]
numPoints  =  50;
stepSize  =  (ArcTanh[0.9]  - ArcTanh[-0.9])/
        numPoints
dataPoints  =  
        Table[Tanh[ArcTanh[-0.9] +  x*stepSize],  {x,  1 ,  
            numPoints  -  1}];
ListPlot[dataPoints,  PlotJoined -> True]

9: The Table method

Similar to the Array method, these two have a lot of overlapping functionality. Evaluate and notice we didn't have to redefine f, because this variable has notebook level scope. (2 pts) Create a second list populated with decimal approximations of f.

myTableGeneratedList  = 
        Table[f[x],  {x, 2, 20, 2}];
myTableGeneratedList

myList  =  
        Table[N[f[x]],  {x,  2,  20,  2}];
myList

10: Recursive functions

Recursive functions take whole numbers as inputs,  while their output variables can be as simple as true/false, whole numbers, or as complicated as you can imagine. A recursive function, rF, is made out of a regular function, f, and a seed value, s. s is usually the same type of thing as the output of f. 
rF(0) = s
rF(1) = f(s)
rF(2) = f(f(s)), or generally speaking, rF(n) = f(rF(n-1)).
(2 pts) In the cell below is the definition of a Mathematica function that doubles at every step. Insert some expressions that evaluate it at a few more steps.  What is the function "f" in this case?

rF[t_] :=  rF[t]  =  rF[t - 1]*2;
rF[0]  =  7;
rF[3]

rF[4]

    rF[5]

    rF[6]

In this case, f(x) = 2x.

11: Modeling population growth

Recursive functions are perfect for modeling systems that change in a step by step way at regular time points. Dynamic systems such as populations of organisms don't change their numbers on a fixed schedule, but recursive functions can be a good rough model, since the population at year n has a big influence on the population at year n+1.  (2 pts) Use the methods in steps 6-9 to create a graph of F defined above. Can you think of a real biological system that F is a good model for?

ListPlot[Table[rF[t],  {t,  0, 10}], PlotJoined -> True]

This population chart models pools of organisims with no predators and insignificant resource constraints, like bacterial cultures in "mid log phase", and people until recently.

12: A more complex model

The logistic model is presented below. It is considered to be the simplest dynamical system that exhibits chaotic behavior. It is implemented here with three parameters: r, k, and g. Set these three numbers to different values, evaluate, and observe how the graph changes. In Mathematica it is necessary to clear a function before you change the parameters that are part of it's definition, which explains the first command: Clear[pop]. Note: use values for r, 0 < r < 4, or you will get overflow errors. Also, when setting r to a whole number, such as 2, put a decimal point after the number: write "r = 2." instead of "r=2" . The latter will cause Mathematica to hang.  
(12 pts) Match the four  parameters to the following definitions:
1) The population at time zero.
2) The number of generations.
3) The capacity of the enviroment to sustain the population.
4) The fertility of the organisims.
Create three copies of this cell. In one cell choose parameters that produce a graph that shows unchecked exponential growth, similar to the graph you created in (11). 
In the second cell, set parameters that produce a graph that shows periodic behavior (i.e. the population changes from generation to generation in a regular way.
In the thrid cell, set parameters that result in chaotic behavior. Which parameter controls the chaotic/non-chaotic behavior of the system?

r - A measure of the fecundity of the population.
k - The carrying capacity of the enviroment: (1-pop[t]/k) becomes negative when the population exceeds k.
g - The number of generations plotted.
s - The population at time 0.

The parameters below show unrestrained population growth.   A large carrying capacity was chosen, and saturation was avoided by choosing small enough k, g, and s.

\(\(\(
\)\(Clear\ [pop]

    \(r\  = \ 1.2;\)

    \(k\  = \ 10000;\)

    \(g\  = \ 30;\)

    \(s\  = \ 2;\)

    \(pop[t_]\  := \ \(pop[t]\  = \ 
          r\ pop[t - 1] \((1 - \ pop[t - 1]/k)\)\);\)

    \(pop[0]\  = \ 5;\)

    \(dt\  = Table[pop[x], \ {x, 0, g, 1}];\)

    ListPlot[dt, \ PlotJoined\  -> True, \ 
      PlotRange -> {0, \ k/16}]\)\)\)

The parameters below show a period 4 cycle. Substituting in r = 3.4 will show a period two,  and r = 3.562 will show a period 8 cycle.  These cycles show up better by selecting large g and by setting PlotJoined->False. Starting with s near k eliminates pre-convergence noise as well.

\(\(\(
\)\(Clear\ [pop]

    \(r\  = \ 3.51;\)

    \(k\  = \ 1000;\)

    \(g\  = \ 1000;\)

    \(s\  = \ 1000;\)

    \(pop[t_]\  := \ \(pop[t]\  = \ 
          r\ pop[t - 1] \((1 - \ pop[t - 1]/k)\)\);\)

    \(pop[0]\  = \ 5;\)

    \(dt\  = Table[pop[x], \ {x, 0, g, 1}];\)

    ListPlot[dt, \ PlotJoined\  -> False, \ 
      PlotRange -> {k/4, \ k}]\)\)\)

Its r values just under four show chaos. There are actually a few values that have limit cycles that are not powers of two. If you pick an r value close enough to one of them, you can easily see them.

\(\(\(
\)\(Clear\ [pop]

    \(r\  = \ 3.97;\)

    \(k\  = \ 1000;\)

    \(g\  = \ 1000;\)

    \(s\  = \ 1000;\)

    \(pop[t_]\  := \ \(pop[t]\  = \ 
          r\ pop[t - 1] \((1 - \ pop[t - 1]/k)\)\);\)

    \(pop[0]\  = \ 5;\)

    \(dt\  = Table[pop[x], \ {x, 0, g, 1}];\)

    ListPlot[dt, \ PlotJoined\  -> False, \ 
      PlotRange -> {0, \ 3  k/2}]\)\)\)

Extra Sections

13: Building recursive lists with NestList

The Mathematica NestList method provides another way to model recursiveness. The NestList second arg is the inital value of the recursive variable, and the third arg is the number of recursive iterations. Notice how NestList can be used with multivariate functions, provided the variable to be treated recursively is listed last.

Clear[g1];
g1[z_]  :=  N[z^2];
NestList[g1,  3,  4]
Clear[g2];
g2[x_][y_][z_]  :=  x*y*z;
NestList[g2[b][c],  a,   5]

14: Sensitivity analysis with list differences

The code in the cell below offers a technique on how to compare two alternative histories. Apply this technique to the Logistic map. Suggested steps: 1) Define a non-recursive version of the pop[] function above, with all of the appropriate parameters.  2) Use NestList to populate two lists, that are each seeded with nearly identical starting populations. 3) Plot the history difference.  Compare results of the history differences in the chaotic and non-chaotic regimes.

Clear[browny];
browny[x_]  :=  
        x  +  Random[Real,  {0, 1}]  - 
          0.5 ;
hist1  =  NestList[browny,  0,  10];
hist2  =  NestList[browny,  0,  10];
ListPlot[hist2 - hist1,  PlotJoined  -> True]

The function fg below is the stripped down logistic map, where the carrying capacity is normalized to 1. Here, if we start with initial conditions far enough apart, in this case 0.5 and 0.75 each history will be at a different step in the period 4 limit cycle. If you choose values that are close enough together, the graph will become very boring, since each trajectory will move thru the same limit cycle in paralel.

Clear[fg];

    r  =  3.51

    fg[x_]  :=  r x (1 - x);

    h1  =  NestList[fg,  0.5,  1000];

    h2  =  NestList[fg,  0.75,  1000];

    ListPlot[h2 - h1,  PlotJoined  -> False];

When we are in the 8 cycle regime, we need to start with nearly the same initial condition to make the histories coincide. Try starting with 0.55, and 0.56 for h1 and h2.  In this case the convergence is slow enough to generate an interesting plot.

Clear[fg];

    r  =  3.562

    fg[x_]  :=  r x (1 - x);

    h1  =  NestList[fg,  0.50,  1000];

    h2  =  NestList[fg,  0.55,  1000];

    ListPlot[h2 - h1,  PlotJoined  -> False];

When we are in the 8 cycle regime, we need to start with nearly the same initial condition to make the histories coincide. Try starting with 0.55, and 0.56 for h1 and h2.  In this case the convergence is slow enough to generate an interesting plot.

Clear[fg];

    r  =  3.562

    fg[x_]  :=  r x (1 - x);

    h1  =  NestList[fg,  0.50,  1000];

    h2  =  NestList[fg,  0.55,  1000];

    ListPlot[h2 - h1,  PlotJoined  -> False];

In the chaotic regime of the logistic function, any two points eventualy have divergent histories. You may not always see it though, due to calculation precision limits. One way to get Mathematica to use, say 400 digits of precision is to put `400 in front of the numbers. If you increase precision from 400 digits to say, 800, the graph will remain chaotic longer.

Clear[fg];
r  =  3.989`400;
fg[x_]  :=  r x (1 - x);
h1  =  
        NestList[fg,  0.700000000000000000000`400,  
          1000];
h2  =  
        NestList[fg,  0.700000000000000010000`400,  
          1000];
ListPlot[h2 - h1,  PlotJoined  -> False];