library(package = "deSolve")

Writining math in RMarkdown

  1. There are two ways of writing math in R, either inline or displayed. In-line is used when you want to write math that is simple and short, like if I wanted to write, \(C = \tau r\), as the circumference of a circle in relation to its radius. To write inline, use two dollar signs and write the math in the middle; e.g., $C = \tau r$ gives you what I just typed. In display mode, you will create a line break, which you use for more complex and important math; e.g., \[a^2 + b^2 = c^2\]. To write math in display mode, you use two dollar signs; e.g., $$a^2 + b^2 = c^2$$. For different symbols, letters, etc. see the Resources page on the course website and follow the link to \(\LaTeX\).
    1. Tell me about your three favorite equations in-line and make sure that you include them.

    My three favorite equations are the logistic equation, \(\mathrm{d}N/\mathrm{d}t = rN - \alpha N^2\), entropy of a random variable as \(H(X) = - \sum _{x \in \chi} p(x) \log p(x)\), and Hardy-Weinberg equilibrium \(p^2 + 2pq + q^2 = 1\).

    1. Using the \(\LaTeX\) guide as a reference, write the three logistic models from the quiz. Separate lines using \\ and, to align the equations, use \begin{aligned} after the first pair of dollar signed and \end{aligned} before the second pair of dollar signed and put an & before each equals sign to align them, so they look something like this: \[ \begin{aligned} \frac{\mathrm{d}N}{\mathrm{d}t} & = N\left(r - \alpha N\right) \\ \frac{\mathrm{d}N}{\mathrm{d}t} & = rN\left(\frac{K - N}{K}\right) \\ \frac{\mathrm{d}N}{\mathrm{d}t} & = rN\left(1 - \alpha N\right) \end{aligned} \] So now you have it—you can write math using RMarkdown!

Looks like this:
$$
\begin{aligned} \frac{\mathrm{d}N}{\mathrm{d}t} & = N\left(r - \alpha N\right) \\ \frac{\mathrm{d}N}{\mathrm{d}t} & = rN\left(\frac{K - N}{K}\right) \\ \frac{\mathrm{d}N}{\mathrm{d}t} & = rN\left(1 - \alpha N\right)
\end{aligned}
$$

Making graphs in R

  1. Before we proceed in making graphs, I want to note that you are free to use any method that you choose, but I would prefer that you use plot(), etc. to make graphs in this course because (1) what you learn in base R is transferable but specialized packages like ggplot() are not and (2) although there is no plotting in base R that I cannot do (feel free to test this out) I cannot say the same about ggplot() and the rest of the tidy universe.
    Onward: making plots in R is fundamentally simple. The three steps you need for 2-dimensional graphs (i.e., have a horizontal and vertical axis) are: (1) create a vector of numbers representing your \(x\) variable, (2) plug \(x\) into an equation, and (3) plot it using plot(x = , y = ) where you plug in the names of your \(x\) and \(y\) variables. That’s really it.
    Let’s try this by plotting a quadratic function. If I create a sequence of values of a variable \(x\), like x_vec <- seq(from = -5, to = 5, length.out = 10) and make \(y\) a function of \(x\), like y_vec <- 0.5*x_vec^2 + 0.1*x_vec - 2. When I plot, plot(x = x_vec, y = y_vec, type = "l", las = 1), I get
    Make this plot.

Here is what the chunk would look like:
```{r echo = F, fig.dim = c(4,4)}
par(mar = c(4.5, 4.5, 0.1, 0.1)) x_vec <- seq(from = -5, to = 5, length.out = 10) y_vec <- 0.5*x_vec^2 + 0.1*x_vec - 2 plot(x = x_vec, y = y_vec, type = "l", las = 1)
```

  1. I want you to make a couple more plots. First, choose one of your favorite equations from above and plot it in this way. Make sure you make the length.out argument within the sequence of \(x\) values longer than just 10: try whatever makes it sufficiently smooth like 100, 1,000, or 10,000. Second, I want you to add 3 additional curves onto a new graph of your favorite equation. For each of the three curves, vary a parameter. For instance, if I use the quadratic function from above, I would plot it and change, say, the coefficient of \(x^2\). To add this to your plot, use the lines() function. Change the color of the curve just to tell the difference between them; e.g.,
x_vec <- seq(from = -5, to = 5, length.out = 100)

y_vec <- 0.5*x_vec^2 + 0.1*x_vec - 2
y_vec1 <- 1*x_vec^2 + 0.1*x_vec - 2
y_vec2 <- 2*x_vec^2 + 0.1*x_vec - 2
y_vec3 <- 4*x_vec^2 + 0.1*x_vec - 2

par(mar = c(4.5, 4.5, 0.1, 0.1)) # sets the whitespace on the bottom, left, top, right
plot(x = x_vec, y = y_vec, type = "l", las = 1, xlab = "x", ylab = "y")
  lines(x = x_vec, y = y_vec1, col = "red")
  lines(x = x_vec, y = y_vec2, col = "#0000FF") # You can use hex codes
  lines(x = x_vec, y = y_vec3, col = "gold")

Here I’ll use Einstein’s mass-energy equivalance, \(E = mc^2\). I am going to use different speeds of light, \(c\), as the parameter I vary.

m_vec <- seq(from = 0, to = 1, length.out = 10000)
c_constant <- 3*10^9
  E_vec_1 <- m_vec*c_constant^2
c_constant_1 <- 3*10^8.9
  E_vec_2 <- m_vec*c_constant_1^2
c_constant_2 <- 3*10^8.8
  E_vec_3 <- m_vec*c_constant_2^2
c_constant_3 <- 3*10^8.7
  E_vec_4 <- m_vec*c_constant_3^2


par(mar = c(4.5, 4.5, 0.1, 0.1)) # sets the whitespace on the bottom, left, top, right
plot(x = m_vec, y = E_vec_1/10^18, type = "l", las = 1, xlab = "Mass (kg)", ylab = "E (joules) / 10^18")
  lines(x = m_vec, y = E_vec_1/10^18, col = "red")
  lines(x = m_vec, y = E_vec_2/10^18, col = "#0000FF") # You can use hex codes
  lines(x = m_vec, y = E_vec_3/10^18, col = "gold")

Visualizing density dependence

  1. Since the beginning of this course includes learning a programming language (R), then I hope you will excuse this problem set for being very R/skills forward. For this last question, it will be doing things you’ve done above, but with logistic equations.

    1. First, plot a time series of the logistic equation. Use any parameter values that you like and use either the \(r\)-\(K\) or \(r\)-\(\alpha\). In lab this week you have already created and plotted a time series, so modify your code (hint). > To detemine what the cumulative populaiton size will be over some time, then use numerical integration to find the solution.
  logisitc <- function(time, N, parms) {
    with(as.list(c(N, parms)), {
    dNdt <- r*N - alpha*N*N
    return(list(c(dNdt)))
    })
  }

  parm_vec <- c(r = 0.4, alpha = 0.005)
  N0 <- c(N = 0.1)
  time_vec <- seq(from = 0, to = 40, length.out = 1000)
  ode_out <- ode(y = N0, times = time_vec, func = logisitc, parms = parm_vec)
  par(mar = c(4.5, 4.5, 0.5, 0.5))
  plot(x = ode_out[,1], y = ode_out[,2], las = 1, xlab = "Time", ylab = "Density (N)", type = "l")

  1. Second, make two plots of your same logistic equation (see last week’s key on how to make plots side-by-side), but where you have the population growth rate on a vertical axis and the per-capita growth rate on the other axis, with both having density on the horizontal axis. This is just like the quiz. But also remember you are not integrating and just plotting some function like what you did in the previous section.

Just like the problems above. There is no need to numerically integrate because I am not asking to to find \(N\) at some time. Instead, I am asking you to plot \(x\) and \(f(x)\):

  parm_vec <- c(r = 0.4, alpha = 0.005)
  N_vec <- seq(from = 0, to = 90, length.out = 1000)
  pgr <- parm_vec["r"]*N_vec - parm_vec["alpha"]*N_vec*N_vec
  pcgr <- parm_vec["r"] - parm_vec["alpha"]*N_vec

  par(mar = c(4.5, 4.5, 0.5, 0.5), mfrow = c(1, 2))
  plot(x = N_vec, y = pgr, las = 1, type = "l", lwd = 2, xlab = "Density (N)", ylab = "Popuation growth rate (dN/dt)")
  plot(x = N_vec, y = pcgr, las = 1, type = "l", lwd = 2, xlab = "Density (N)", ylab = "Per-capita growth rate (dN/dt/N)")

  1. Last, recreate parts of the Figure 5.12 from Stevens link. Choose a value of \(r_d\) that gives you asymptotic stability, oscillations (damped or perpetual), and chaos. Use the numerical solver you used in a. above, but in the ode() function, write the argument method = "iteration".
 logisitc_diff <- function(time, N, parms) {
    with(as.list(c(N, parms)), {
    Ntp1 <-  rd*N*(1 - alpha*N) + N
    return(list(c(Ntp1)))
    })
  }

  parm_vec_1 <- c(rd = 1.3, alpha = 0.01)
  parm_vec_2 <- c(rd = 1.9, alpha = 0.01)
  parm_vec_3 <- c(rd = 2.2, alpha = 0.01)
  parm_vec_4 <- c(rd = 2.8, alpha = 0.01)
  N0 <- c(N = 10)
  time_vec <- 0:20
  ode_out_rd1 <- ode(y = N0, times = time_vec, func = logisitc_diff, parms = parm_vec_1, method = "iteration")
  ode_out_rd2 <- ode(y = N0, times = time_vec, func = logisitc_diff, parms = parm_vec_2, method = "iteration")
  ode_out_rd3 <- ode(y = N0, times = time_vec, func = logisitc_diff, parms = parm_vec_3, method = "iteration")
  ode_out_rd4 <- ode(y = N0, times = time_vec, func = logisitc_diff, parms = parm_vec_4, method = "iteration")

    par(mar = c(4.5, 4.5, 0.5, 0.5), mfrow = c(2,2))
  plot(x = ode_out_rd1[,1], y = ode_out_rd1[,2], las = 1, xlab = "Time", ylab = "Density (N)", type = "l")
  plot(x = ode_out_rd2[,1], y = ode_out_rd2[,2], las = 1, xlab = "Time", ylab = "Density (N)", type = "l")
  plot(x = ode_out_rd3[,1], y = ode_out_rd3[,2], las = 1, xlab = "Time", ylab = "Density (N)", type = "l")
  plot(x = ode_out_rd4[,1], y = ode_out_rd4[,2], las = 1, xlab = "Time", ylab = "Density (N)", type = "l")