Blog of Andrés Aravena

Homework 8

21 April 2019. Deadline: Friday, 26 April, 9:00.

Please download the file homework8.R and write your results there. Send the your answers to my mailbox.

1. Replicate plots from last class

In the last class we simulated the water formation system under different rates and initial conditions. The system is represented by this graph: water_system

The code to make one simulation is this:

water_formation <- function(N, r1_rate,  r2_rate, H_ini, O_ini, W_ini) {
  W <- d_W <- rep(NA, N) # Water, quantity and change on each time
  H <- d_H <- rep(NA, N) # Hydrogen
  O <- d_O <- rep(NA, N) # Oxygen
  W[1] <- W_ini
  H[1] <- H_ini
  O[1] <- O_ini
  d_W[1] <- d_H[1] <- d_O[1] <- 0 # the initial change is zero

  for(i in 2:N) {
    d_W[i] <-    r1_rate*H[i-1]*H[i-1]*O[i-1] - r2_rate*W[i-1]
    d_O[i] <-   -r1_rate*H[i-1]*H[i-1]*O[i-1] + r2_rate*W[i-1]
    d_H[i] <- -2*r1_rate*H[i-1]*H[i-1]*O[i-1] + 2*r2_rate*W[i-1]

    W[i] <- W[i-1] + d_W[i]
    O[i] <- O[i-1] + d_O[i]
    H[i] <- H[i-1] + d_H[i]
  return(data.frame(Step=1:N, W, H, O))

Please write the code to draw the following plot. You will need to simulate several times for different rates and initial conditions

2. Simulate a chaotic system

In the class we also saw a chaotic system, that can have high sensitivity to initial conditions. The system is very simple


After simplification, this system can be simulated by the following code:

quad_map <- function(N, A, x_ini) {
  x <- rep(NA, N)
  x[1] <- x_ini
  for(i in 2:N) {
    x[i] <- A * x[i-1] * (1 - x[i-1])

As we said, the behavior depends on the value of A. Please write the code that simulates this system for different values of A from 2.9 to 4 by 0.001, with x_ini=0.4 and N=1500. Take the final 500 values of each simulation and draw them using pch=".". You should see this:

There are at least two ways to do this plot.

3. Many dice

In class we simulated throwing two dice at the same time. What if we throw 3 dice? Or 10? Or 100?

Please write the code to simulate N times the experiment “throw M dice and count the total of points”. Something like

many_dice <- function(N, M) {
    # repeat N times
        # throw M dice
        # sum them
        # put the result in a vector `ans`
    # return `ans` vector

If you do it right, you should be able to draw a plot like this:

Deadline: Friday, 26 April, 9:00.