library(package = "deSolve")
library(package = "viridis")

1. 1D Bifurcation Diagram

This problem set is going to be a bit shorter given that you will be giving your remodeling paper presentation this week. I really only want you to make 2 graphs—2 bifurcation plots; 1 in each of #1 and #2. You’ll be doing problem 1 with a different equation this week in class. But for problem 2, you can do with what you know now.

For the first bifurcation plot, I would like for you to recreate the figure from the Ricker model Wikipedia page. You’ll notice this is exactly the equation you learned in chapter 1, but we used \(b\) and they use \(1/K\) as the crowding/density-dependence parameter.

We’re going to cover this for the discrete-time logistic equation, so you should at least have an example to follow by Tuesday. In short, a bifurcation diagram is one that shows a qualitative change in dynamics. “Qualitative” hear, meaning the point of the graph isn’t supposed to show a quantitative change (e.g., the equilibrium increasing), but rather that there are differences in the types of dynamics.

We can see at small values of \(r < 2\), in this instance, that the populations reach ~1,000 (what we set as \(K\)) and that’s it. This is what you see what you for the behavior of our “normal” logistic equation, where populations change and approach \(K\), and just stay there.

But at slightly larger values of \(r\), \(\sim 2.5 > r > 2\), you can see, reading upwards from any point of \(r\), that there are two points. This is because the population is cycling! It’s effectively overshooting, undershooting, overshooting, etc. \(K\) and is always returning to the same two values.

At slightly larger values of \(r\), \(\sim 2.6 > r > \sim 2.5\), you can see that the population is cycling but between four points!

And finally, \(r > \sim 2.6\), you can see that there are some regular cycles, but the points otherwise appear to be random: this is chaotic dynamics.

So, this is called a 1-dimensional bifurcation because we are showing the qualitative change in dynamics—asymptotic stability, 2-point cycles, 4-point cycles, chaos, etc.—by changing a single variable, \(r\).

To run this analysis, it’s analogous to what you did in the last problem set where you create a bunch of \(r\) values, and loop over them to numerically solve them. There are two important differences, however:

  1. It’s discrete-time, so for ode(), you want to set the argument method = "iteration".
  2. For the logistic model you analyzed in problem set 3, we just saved the final value of our time series since we solved over “a long period of time” and we know that, for the \(r-\alpha\) logistic equation, that it asymptotically approaches \(r/\alpha\). If we just looked at the last value for this Ricker model, the figure below is what we’d get! Notice that we miss so much! For 2-point cycles, for instance, we just see 1/2 of the cycle! This is why, for this bifurcation diagram, we want to save more than 1 final value. For me, I usually use at least the last 10 values. What this means, is that instead of saving one value of final densities per \(r\) in a vector, we need to save the last, say, 10 values for each \(r\). A matrix is ideal for this, where you create one with, say, the number of rows being the number of \(r\) values, and the number of columns being the number of final densities you are saving.
  3. For regular plots, you plot a vector of x values and a vector of y values. Here, it’s a bit different; here, you want to plot several y values (like 10) per \(r\). To do this, you still need to give plot() vectors of the same length, but you want to repeat the same \(r\) value for every final density (e.g., 10).
  4. Last, aesthetics. I like to make the points smaller using the character expansion argument cex = to some value less than 1 (>1 makes them bigger!). Since many points overlap, I like to make them translucent. We can use hexcolors for this, and add how much transparency we want at the end from 00 (fully transparent) to FF (fully opaque) at the end. If I wanted, for instance, red in hexcolor, I would write col = "#FF00000". And if I wanted to make it medium transparent, I would write col = "#FF0000077".

Remember that this is a bifurcation plot because its purpose is to show a qualitative change in dynamics.

Again, we will go over this on Tuesday, so you’ll be able to do this problem after then.

2. 2D Bifurcation Diagram

Second, and last for this week, I want you to extent the idea of a bifurcation diagram from 1 to 2 dimensions. We aren’t going to really going to look at qualitatively different type of dynamics, but we are going to look at the output of our model qualitatively. This means that we aren’t going to identify cycles, extinction, chaos, etc., but rather just take our 2D heat maps where we change two parameters, and qualify/bin the values. You’d do the same with dynamics, but that’ll be learned later.

So, for this question, I want you to take the analysis that you did for question 3.3., from problem set 03 and run it. From there, arbitrarily choose an value within the range of your matrix of densities, assign it to a group (e.g., an arbitrary value) and the rest to another group (e.g., an arbitrary value).

For my question 3.3. in the key, for instance, I am going choose 10. We can pretend that we curious about the sets of parameter values that give us an equilibrium of 10. Here’s what the beat made and 2D bifurcation diagram would look like side-by-side:

Let’s see how to do this on an example matrix. First, let’s create a matrix with some values. I’m going to do pair-wise multiplication between 1:3 and 5:8. (You don’t need to do this, I’m just making a toy example to, if you will, f(around + find out),)

vec1 <- 1:3
vec2 <- 5:8
toy_mat <- outer(X = vec1, Y = vec2, FUN = "*")
toy_mat
##      [,1] [,2] [,3] [,4]
## [1,]    5    6    7    8
## [2,]   10   12   14   16
## [3,]   15   18   21   24

Now, if we wanted, say, to make of graph of which values are < 13: how would we do this? First, let’s try just using logical operators like <, >, ==, !=, etc. If we run toy_mat < 13, we’ll see:

##       [,1]  [,2]  [,3]  [,4]
## [1,]  TRUE  TRUE  TRUE  TRUE
## [2,]  TRUE  TRUE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE

You can correspond that to the matrix above for values that are less than < 13—i’n’it cool!?

So, now we want to assign all of those that are > 13 as something: let’s do the number 1. To do so, you index your matrix with the output of the logical operators and assign it 0. This is what toy_mat[toy_mat < 13] <- 0 would like:

##      [,1] [,2] [,3] [,4]
## [1,]    0    0    0    0
## [2,]    0    0   14   16
## [3,]   15   18   21   24

Now, do this for values ≥ 13 and assign them to 1. Your toy example should look like this when plotted with image():

##      [,1] [,2] [,3] [,4]
## [1,]    0    0    0    0
## [2,]    0    0    1    1
## [3,]    1    1    1    1

So, do this over your \(r\)-\(\alpha\) matrix you made in 3.3 for last week’s problem set with some arbitrary threshold.