Math Insight

Exploring a discrete dynamical system with R

Math 2241, Spring 2023
Name:
ID #:
Due date: Jan. 25, 2023, 11:59 p.m.
Table/group #:
Group members:
Total points: 1

Background

Mountain Nyala
Mountain nyala, Tragelaphus buxtoni
(Image: © KrisMaes, CC BY-SA 3.0 wikimedia)

The mountain nyala (Tragelaphus buxtoni) is a type of antelope found in central Ethiopia. At this point, there are only about 2,500 mountain nyalas around and their population is decreasing, due to factors such as habitat loss and hunting.

Let's develop a simple discrete dynamical system of the mountain nyala population and explore the behavior of the model using the programming language R.

  1. The first step is to introduce notation for the model. We are just going to keep track of one quantity, the size of the mountain nyala population. Hence, we just need one state variable. We'll call it $m$ for mountain nyala, and to remind ourselves that mountain nyalas are distinct from nyalas, which aren't endangered. The population size $m$ is going to vary with time, so we need notation for time as well. Let let $t$ be time in years. Since we are going to create a discrete dynamical system, we are only going to look at snapshots of the mountain nyala population. We'll imagine we count the number of mountain nyalas each year. We number our year $t=0,1,2, \ldots$, where year $0$ is the initial year. The mountain nyala population size in year $t$ is represented by the number $m_t$. Therefore, our initial population size is the number $m_0$, the population size in year 1 is $m_1$, and the population size in year 5 is
    .
    1. To form a dynamical system, we need a rule that gives the population size in one year as a function of the previous year's population size. We will assume that environmental or other factors are not changing from year to year, so that the only thing we need to keep track of is the population size $m_t$. We will create a rule so that if we know the population size $m_7$ in year 7, we'll have a formula to compute the population size $m_8$ in the next year. In general, if we know the population size $m_t$, we want a formula for the next year's population size $m_{t+1}$.

      Let's imagine that every year, 20% of the mountain nyalas give birth to a baby mountain nyala and that half of those babies survive to the next year. Therefore, the increase in the population due to new births is
      %. However, if we imagine that each year, 6% of the remaining mountain nyalas (not including the new babies) die each year, the population is actually increasing by only
      %. But, the situation looks good for the mountain nyalas, as the population is increasing (unlike in reality, but we'll get back to that later).

    2. Let's say we start with $m_0=900$ mountain nyalas in year $0$. Above, you determined that the population should increase by $_$% each year. This means, to calculate the predicted population for year $1$, you need to multiply the number $900$ by how much?
      . Similarly, you should multiply this year $1$ population by $_$ to calculate the year $2$ population, and similarly for future years.

      Therefore, according to this model, we should have $m_1 =$
      in year 1, $m_2=$
      in year 2, and
      in year 3.

    3. The next step is to transform these calculations into a dynamical rule. If you didn't know what the mountain population was in year $0$ but just represented that number by $m_0$, we can still determine an expression for the year $1$ population size $m_1$; it will just be in terms of the quantity $m_0$. Given that the population size is growing by _% each year, we can calculate that $m_1 = $
      .
    4. Similarly, if the population size in year $8$ is represented by $m_8$, the formula for the next year's population size must be $m_9 = $
      .
    5. We could write that same formula to get $m_{12}$ from $m_{11}$ and $m_{49}$ from $m_{48}$. But to save us from writing all these different formulas, we can summarize them all by writing a general formula to get from year $t$ to the next year $t+1$. We just need a formula for $m_{t+1}$ in terms of the previous year's population size $m_t$. That formula is
      $m_{t+1} = $
      , for $t=0,1,2, \ldots$

      If we know this dynamical rule combined with the initial condition $m_0$, we can calculate the population size for any $m_t$ just by applying the rule first for $t=0$ to calculate $m_1$, then applying it for $t=1$ to calculate $m_2$, and repeating the process until we get to the year $t$ we want.

      To summarize, the initial dynamical system we developed (which includes the dynamical rule and the initial condition is): \begin{align*} m_{t+1} &= _\\ m_0 &= 900. \end{align*}

  2. Let's use the programming language R to simulate the mountain nyala population. We assume that you have installed R as well as RStudio on your computer and have opened up RStudio. You should have one window labeled “Console” with a > prompt where you can type commands. We'll start by using the console directly.

    1. First, let's define a variable $r$ for the growth rate of the population. Since, at this point, we are thinking of the population as growing by 4% per year, define the variable $r$ by typing r= 0.04 in the console. (In the Environment window of RStudio, you should see that a variable $r$ has been defined.) Define the initial condition by typing m0=900. (To make it easier to type, we don't need to use underscores as in m_0 in R.)

      (In R, people often type the assignment commands as

      r <- 0.04
      m0 <- 900
      
      using “<-” rather than “=” for the assignment operator. Either way works.)

      Above, though, you didn't multiply by the number $r$. Instead, you multiplied by a different number to go from year $t$ to year $t+1$. Let's call this other number $R$ (not to be confused with the name of of the programming language R). What is the relationship between $R$ and $r$?
      $R = $

      Type this equation into the R console to define the variable $R$.

    2. Above, you wrote the dynamical rule giving $m_{t+1}$ in terms of $m_t$, where in each step you multiplied by a number. Rewrite that dynamical rule in terms of $R$ rather than using a number. The new dynamical rule is
      $m_{t+1} = $
      for $t=0,1,2, \ldots$
    3. Use R to calculate $m_1, m_2, \ldots, m_{10}$.

      To calculate $m_1$, you should type in the Console a formula in terms of R and m0. This formula is
      m1 =

      (In R, you will need to use * for multiplication, though the above answer blank will accept your answer even without an *.)

      Similarly, the equation for $m_2$ is
      m2 =

      Continue this process to calculate m3 all the way until m10, which will be our R variable for $m_{10}$, the population size in year 10. (Hint, you can press the up arrow in the Console to recall the previous command, which you can edit to create a new command to enter.)

      When you finish this procedure, what value do you get for $m_{10}$?
      (In the console, you can just enter m10 to see the current value of the variable.)

      In 10 years, the mountain nyala population size increased by what percent?
      %. This value is
      the value of $4 \times 10 = 40$, which you might expect for 10 years' growth at 4%.

    4. In R, the command c combines its arguments into a vector (see vector creation in R). You can create a vector of all the population sizes from m0 to m10 by entering:
      ms = c(m0,m1,m2,m3,m4,m5,m6,m7,m8,m9,m10)
      into the console. We can also create a vector of all times $0,1,2, \ldots, 10$ by typing ts = 0:10 into the console. Then, to plot the population sizes ms versus the times ts, enter plot(ts, ms) into the console. Or, to be more fancy and have a graph with readable labels (see plotting in R), enter the following plot command (all as one line).
      plot(ts, ms, xlab="time", ylab="population", main="Mountain nyala population versus time")

      What trend do you observe about the change in the mountain nyala population over 10 years?

    5. According to the model, will this trend continue for a long time? Let's simulate the model for 100 years. It'd be too much of a pain to continue typing the commands in by hand to calculate m11 to m100. Instead, we'll create a for-loop to do this.

      First, just in case we change our mind about 100 years, let's set a variable $t_{max}$ to 100 by typing tmax=100 in the console. That's good practice so that it will be easy to play around with different times by simply changing $t_{max}$ (and we won't have to change any other commands).

      It'll be much easier to switch to using the vector ms rather than the individual variables m0, m1, etc. for our population sizes. One note of caution, though. R starts numbering vectors with the index 1 (rather than starting with 0 like some other languages). So, to get the value of $m_0$, we need the first entry of ms, which we access by ms[1]. This means that the vector component ms[2] represents the value of
      , the component ms[6] contains the value of
      , and that ms[t+1] is for
      . To summarize, we must remember that we need to add one to the value of $t$ to get the proper index of ms.

      Using the vector ms, write a for-loop to simulate our model $$m_{t+1} = Rm_t \quad \text{for } t=0,1,2, \ldots, t_{max}-1.$$ (We only need to go to $t=t_{max}-1=99$, as then $t+1=t_{max}$.) You can read about for-loops in R for guidance on the syntax. The idea is to make a loop that is equivalent to typing in all the commands ms[2]=R*ms[1], ms[3]=R*ms[2], all the way until ms[101]=R*ms[100], where the last command calculates $m_{100}$ and stores it in ms[101]. (It's OK that the loop will overwrite the values of ms[2] through ms[11] that we calculated before, and R will automatically extend the vector ms to its final length of 101 entries.)

    6. Create a longer version of the ts vector to include all times to $t=100$ with the command ts=0:tmax. Now reenter the above plotting command to plot the population size $m_t$ for all years from $t=0$ to $t=100$.

      What does the model predict for the population size $m_{100}$ in year 100?
      (Remember, you can just type in a variable in the console to see its value. What component of the vector ms stores $m_{100}$?)

      In 100 years, the mountain nyala population size increased by what percent?
      %. This value is
      the value of $4 \times 100 = 400$, which you might expect for 100 years' growth at 4%.

      What trend do you observe about the change in the mountain nyala population over 100 years?

  3. As a reminder, if the mountain nyala population size 900 in year 0 and is increasing by 4% per year, then the dynamical system is
    $m_{t+1} =$

    $m_0 =$
    .
    1. To obtain $m_1$ from $m_0$, we multiply by
      . To obtain $m_2$ from $m_1$ we multiply by that same number
      . Therefore, if we want to obtain $m_2$ directly from $m_0$ without going through the intermediate step of calculating $m_1$ first, we should multiply by _ how many times?
      In other we should multiply by _ to the power of
      , which is
      .

      In other words, we can write $m_2$ as a function $m_0$ without any reference to $m_1$ as $m_2=$
      .

    2. To obtain $m_3$ from $m_2$, we multiply by _ one more time. Therefore, we need to multiply $m_0$ by _ to the power of
      to obtain $m_3$.

      We can write $m_3$ as a function $m_0$ as $m_3=$
      .

    3. In general, for any number $t$, to obtain $m_t$ from $m_0$, we need to multiply by _ to the power of
      .

      The resulting formula is

      $m_t = $
      ,

      which we call the solution to the dynamical system. It's called the solution, as it is a formula the population size $m_t$ that is just a function of $t$ and $m_0$ (and not a function of the population at any other time). With the solution in hand, we can compute $m_t$ for any $t$ directly.

      The solution is an exponential function. We say that the population size is growing exponentially, i.e., we have exponential growth.

    4. For a general initial condition $m_0$, calculate the population size after 100 years. $m_{100} = $
      . Now, substitute $m_0=900$ to obtain a number for this population size $m_{100} = $
      . This answer should agree with the one you computed by simulating it in R, above.

  4. So far, we have assumed that the mountain nyala population has been increasing by 4%, or a factor of $r=0.04$, each year. This means, each year, we multiply the population size by $R=$
    . We came up with the 4% by assuming a particular birth rate and survival rate.

    Solving the dynamical system isn't that much harder if we don't assume that $R=_$. Let's just leave $R$ as the unknown value $R$ and say that the mountain nyala population size is multiplied by the number $R$ each year. In other words, let's use the dynamical rule

    $m_{t+1} =$
    1. We can still solve this dynamical system. If we start with an initial population size of $m_0$, to calculate $m_t$, we need to multiply by $R$ how many times?
      times. In other words, our solution is

      $m_{t} =$
      .
    2. Do we still get exponential growth for the solution? That depend on the value of $R$. We are going to assume that $R$ is always positive. For what values of $R$ does the quantity $R^t$ get larger as $t$ increases?
      (We are expecting an inequality involving $R$.) For these values of $R$, we get exponential growth.

      On the other hand, if $0 \lt R \lt$
      , then the quantity $R^t$ gets smaller as $t$ increases. For these values of $R$, we get exponential decay.

      There is one special case where we get neither growth nor decay. That is for the value $R=$
      , in which case $R^t$ is always equal to
      .

    3. Let's assume, as above, that every year, 20% of the mountain nyalas give birth to a baby mountain nyala and that half of those babies survive to the next year. Let $d$ equal the fraction of the non-babies that die each year. If we want to mountain population to be increasing, what condition do we need on $d$ (which we assume is a positive number)?

  5. Up to this point, we've been assuming the simple dynamical rule $m_{t+1} = R m_t.$ In terms of the growth rate $r$, where $R=r+1$, we could also write the rule as $m_{t+1} = (r+1) m_t$ or even as $m_{t+1} = m_t + r m_t.$

    Keeping $r$ equal to a constant assumes that the mountain nyalas have plenty of space and food so that their growth rate stays at $r$ even as the population size increases. That's why, if $R > 1$ (equivalently, if $r > 0$), we get an exponentially growing population that eventually explodes to huge numbers.

    Let's take a look at what happens if we introduce limited resources into the model.

    1. We are going to modify the dynamical rule to include a carrying capacity that will slow down the growth when the population size gets large.

      Let $K$ be the carrying capacity, which we can think of being the number of mountain nyalas that their habitat can sustainably support (for example, due to a limited space or supply of food). If the current population size is $m_t$, what is the fraction of the carrying capacity that is being used?
      Therefore, what is the fraction of the carrying capacity that is still available?

      To modify the dynamical system to include the effects of the carrying capacity, we simply multiply the growth rate $r$ by the fraction of the carrying capacity available: $_$. The resulting model is called the logistic model, which we can write as $$m_{t+1} = m_{t} + r m_{t} _.$$

    2. Let's use R to explore the behavior of the logistic model. This time, we'll use a script file to make it easier to keep track and reuse our command. The first step is to create a new directory (or folder) on your computer for the scripts, for example, an “RWork” directory in your Documents directory. Then, in RStudio, set the working directory to that new directory by selecting Set Working Directory from the Session menu. Select Choose Directory and navigate to your new directory and click Open. In the R console, you will see a setwd (for set working directory) command being automatically entered. If you like, you can just type that setwd command in the future.

      Next, open up a new script file. You can select New File from the File menu, then click R Script. Save the file by selecting Save from the File menu. It should automatically show your new directory, as it is the working directory. Save the script with an appropriate name (maybe one expressing your affection for mountain nyalas).

      Let's start by setting the carrying capacity to 5000. First, enter a comment to remind you what is the carry capacity. In R, everything after a pound sign # is ignored. Then set K to 5000. For example, you can enter the two lines:

      # carrying capacity
      K = 5000
      

      To actually run the second line in R, press Ctrl/Cmd-Enter (or click the Run icon) while your cursor is on that line. You'll see that K=5000 is automatically entered in the console. Similarly, set the growth rate $r$ to 0.04 and $t_{max}$ to 100 by typing the lines

      # growth rate
      r = 0.04
      
      # simulation time
      tmax = 100
      

      To run the code, you can press Ctrl/Cmd-Enter on each line or just highlight the lines with your mouse before pressing Ctrl/Cmd-Enter.

      As before, we'll use the vector ms to store the values of $m_t$. We can start by just setting it to our initial condition $m_0=900$ mountain nyalas

      # initial condition
      ms = 900
      

      Be sure to run the second line. (In R, a single number like 900 is the same as a vector of one dimension, so we'll subsequently be able to add components to ms by assigning values to ms[2], ms[3], etc., as we did before.)

      Write a for-loop to simulate the evolution of our logistic model. The only difference will be to change the dynamical rule to the new logistic model you derived above. See if you can figure out how to modify the above for-loop so that R will calculate the solution $m_t$ for $t=1, 2, \ldots, t_{max}$ and store the result in the vector ms. If you get stuck and there's no one around you who can help, you can peak at the code in the below hint.

      Be sure to make a comment of what the loop is doing. Assuming you wrote the loop on more than one line, you can select the lines with the mouse before pressing Ctrl/Cmd-Enter to run all the lines together.

      Once you've run these commands, the vector ms is defined. You can examine ms or one of its components by entering them in the console. (Remember, ms[17], the 17th component of the vector ms actually stores $m_{16}$, the population in year 16.) According to the model, what will be the nyala population size in year 100?

      Complete your script by adding a line to plot the results. The same plot command as before will work.

      # plot the results
      plot(ts, ms, xlab="time", ylab="population size", main="Mountain nyala population versus time")
      

      What trend do you observe in the mountain nyala population size over 100 years?

      You've done a lot work. Be sure to save your script!

    3. To the nearest integer, what does the model predict for the nyala population size after 300 years?

      This number is exactly one of the parameters from out model. Which one?

    4. Lastly, imagine that the mountain nyala have been around for a long time, so their population size is equal to the ending population size you calculated above, i.e., _. Let's use this as an initial condition, i.e., set $$m_0 = _$$

      Imagine that, habitat loss due to human activity has reduced the carrying capacity to $K=2000$. Using your R program, estimate the mountain nyala population size after twenty years. $m_{20} = $