*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:

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])
}
return(x)
}
```

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.

- Make a very long vector with 500 values from each simulation, and plot all in one command
- Draw first an empty plot using
`type="n"`

, and then use the function`points()`

to add the 500 values of each simulation separately

# 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: