-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy path04_ODEs_SIR_code.qmd
55 lines (46 loc) · 1.45 KB
/
04_ODEs_SIR_code.qmd
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
---
title: "04. Ordinary differential equations (ODEs): SIR model"
---
## Practical 1: Susceptible-Infectious-Recovered model implementation
```{r, eval = TRUE, message = FALSE}
library(deSolve) # For solving systems of ODEs
# Define model function
SIR_model <- function(times, state, parms)
{
# Get variables
S <- state["S"]
I <- state["I"]
R <- state["R"]
N <- S + I + R
# Get parameters
beta <- parms["beta"]
gamma <- parms["gamma"]
# Define differential equations
dS <- -(beta * I / N) * S
dI <- (beta * I / N) * S - gamma * I
dR <- gamma * I
res <- list(c(dS, dI, dR))
return (res)
}
# Define parameter values
parms <- c(beta = 0.4, gamma = 0.2)
# Define time to solve equations
times <- seq(from = 0, to = 100, by = 1)
# Define initial conditions
N <- 100
I_0 <- 1
R_0 <- 0
S_0 <- N - I_0 - R_0
y <- c(S = S_0, I = I_0, R = R_0)
# Solve equations
output_raw <- ode(y = y, times = times, func = SIR_model, parms = parms)
# Convert matrix to data frame for easier manipulation
output <- as.data.frame(output_raw)
# Plot model output
plot(output$time, output$S, type = "l", col = "blue", lwd = 2, ylim = c(0, N),
xlab = "Time", ylab = "Number")
lines(output$time, output$I, lwd = 2, col = "red", type = "l")
lines(output$time, output$R, lwd = 2, col = "green", type = "l")
legend("topright", legend = c("Susceptible", "Infectious", "Recovered"),
col = c("blue", "red", "green"), lwd = 2, bty = "n")
```