-
Notifications
You must be signed in to change notification settings - Fork 0
/
simple_predator_pray.Rmd
94 lines (76 loc) · 1.93 KB
/
simple_predator_pray.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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
---
title: "Untitled"
author: "Enzo Moraes Mescall"
date: "`r Sys.Date()`"
output: html_document
---
```{r setup, include=FALSE}
library(deSolve)
library(phaseR)
library(tidyverse)
```
```{r}
predator.pray.model = function(t, State, Pars) {
with(as.list(c(State, Pars)), {
dPredator = a*x - b*x*y
dPrey = d*x*y -c*y
equations = c(dPredator, dPrey)
return(list(equations))
})
}
# Time vectors, parameters, and initial condition
t = seq(0, 100, 0.1)
Pars = c(a = 0.3,
b = 0.05,
c = 0.4,
d = 0.1)
X_0 = c(x = 4, y = 2)
```
```{r}
ode(
func = predator.pray.model,
y = X_0,
times = t,
parms = Pars
) %>% as.data.frame() -> out
```
```{r}
out %>%
gather(variable, value, -time) %>%
ggplot(aes(x = time, y = value, color = variable)) +
geom_line() +
theme_classic() +
labs(x = 'time (yr)',y = 'populations')
```
```{r}
a_vals = c(0.3, 0.7, 1)
b_vals = c(0.05, 0.1, 0.3)
expand.grid(a = a_vals, b = b_vals) %>%
group_by(a, b) %>%
do({ode(func = predator.pray.model,
y = X_0,
times = t,
parms = c(a=.$a, b=.$b, c = 0.4, d = 0.05)) %>%
as.data.frame()}) %>%
ggplot(aes(x = time)) +
geom_line(aes(y = x), color = "red") +
geom_line(aes(y = y), color = "blue") +
facet_grid(a ~ b, scales='free_y', labeller=label_both) +
labs(y = "Population")
```
```{r}
flowField = flowField(predator.pray.model,
xlim = c(0, 20),
ylim = c(0, 30),
parameters = Pars,
add = FALSE)
trajectory1 = trajectory(predator.pray.model,
y0 = c(10, 2),
tlim = c(0, 3000),
parameters = Pars,
col = "purple")
nullcline1 = nullclines(predator.pray.model,
xlim = c(0, 20),
ylim = c(0, 30),
parameters = Pars)
```