12 Thurs Oct 8: Solving the logistic growth equation using a computer 2

DUE DATE: Thursday Oct 15

Last class we learned how to make a custom function. Today we would like to write a custom function that satisfies the requirements of the func argument that is supplied to ode() so that numerical integration can be performed.

12.1 Writing the func argument for ode()

The func argument of ode() has the following constraints:

  • The arglist must include t, y, and parms.
  • The value returned must be a list of the rates of change in the dependent variable(s).

Note that if your ODE comprises of more than one variable that is changing over time, then y will be a vector. Below is an example of how func might be coded for exponential growth. Remember that the ODE for exponential growth is:

\[\begin{equation} \frac{dN(t)}{dt} = rN(t) \tag{12.1} \end{equation}\]

Below equation (12.1) is coded in the appropriate syntax for eventual use in the ode() function:

# Write a function that returns a list of the rate of change in your dependent variables.
model = function(t,y,parms){
  # It's a personal perference of mine to switch the symbols
  # to be consistent with the equation.
  N <- y[1]
  # This is this equation for dN/dt
  dN <- r*N
  # The function returns the value of the change in N.
return(list(c(dN)))
}

Let's try out the function to make sure it works. Copy and paste the function into your Console and then type:

r<-1
model(1,2,parms=NULL)

Did you get 2? Does it make sense to you that the returned value is 2? Try some other values. Are the values making sense?

Now that the func argument of ode() has been demonstrated, the remaining arguments for the ode() function are not so bad.

# The initial value of N
yini<- c(N = 1)

# the times for the numerical integration
times <- seq(0, 1, by = .1)

Rembember to assign a value for r:

# Asign a value to the parameter:
r <- 2

And then once all the arguments are set, call the function, and plot the results:

# performing the numerical integration
out <- ode(y = yini, parms = NULL, times = times, func = model)
out <- data.frame(out)
head(out)

# this line overrides a multipanel plot. I had tried plot(out) # and was having conflicts with a small plotting window on my
# laptop. Generally, the below line isn't necessary.
par(mfrow=c(1,1))

# Make the graph
plot(out$time, out$N, typ = "l")

12.2 Questions

  1. For the code model(1,2,parms=NULL) (where the function model is defined above), what do 1 and 2 correspond to? One of these quantities does not affect the value that the model() function returns. Which quantity does not affect the value the function returns? [2 marks]

  2. Write an R Script that numerically integrates the logistic growth function (equation (11.3)). Choose a value of \(r\) such that the population size increases over time. Remember to use good practices for writing your script - see 1.7.2. Your script should produce a graph with \(N(t)\) on the y-axis and time on the x-axis [10 marks]

  3. Add to your script from 2. such that another graph is produced, this time where the population size decreases over time. Note that this should just require adding new lines to the bottom of your question 1 script such that you: a) set a new value of r; b) run the ode() command; and c) make the graph. In other words, it is not necessary to define the logistic growth function, model, again because that has not changed. Further, if you are happy to use the same values of times as for question 1, you do not need to restate this either. Answer this question by adding a few lines of code to your answer to 1. You are to hand in your R script. You should answer both questions 1. and 2. in the same R script. [5 marks]