-
Notifications
You must be signed in to change notification settings - Fork 0
/
density-dependence.Rmd
78 lines (60 loc) · 1.37 KB
/
density-dependence.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
---
title: "Density dependence"
output: html_document
editor_options:
chunk_output_type: console
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Packages
```{r}
library(purrr)
library(deSolve)
library(tibble)
```
## Functions
```{r}
seq2 <- function(...) seq(..., le = 512)
plot2 <- function(..., col = 4) plot(..., type = "l", lwd = 2, col = col)
lines2 <- function(..., col = 4) lines(..., lwd = 2, col = col)
ode2 <- function(...) as_tibble(as.data.frame(deSolve::ode(...)))
```
## Models
Model 1:
$$
\frac{dN}{dt} = \frac{1}{(1 + N)^\alpha} b N
$$
Model 2:
$$
\frac{dN}{dt} = \left(1 - \frac{N}{K}\right) b N
$$
```{r}
xs <- seq2(0, 10)
plot(NA, xlim = range(xs), ylim = 0:1, xlab = "density", ylab = "per capita rate")
walk(c(.1, .5, 1, 2, 3), ~ lines2(xs, 1 / (1 + xs)^.x))
title("model 1")
```
```{r}
xs <- seq2(0, 10)
plot(NA, xlim = range(xs), ylim = 0:1, xlab = "density", ylab = "per capita rate")
walk(c(1, 3, 5, 8, 10), ~ lines2(xs, (1 - xs / .x)))
title("model 2")
```
```{r}
model <- function(N0, b, alpha, times) {
ode2(c(N = N0),
times,
function(time, state, pars) {
with(as.list(c(state, pars)), {
dN <- b * N / (1 + N) ^ alpha
list(dN)
})
},
c(b = b,
alpha = alpha))
}
```
```{r}
with(model(1, 1, 1, seq2(0, 500)), plot2(time, N, ylim = c(0, 1000)))
```