Speed of ODE solving in R and Python

# Load the deSolve package
library(deSolve)

# Define the SIR model
sir_model <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    dS <- -beta * S * I
    dI <- beta * S * I - gamma * I
    dR <- gamma * I
    return(list(c(dS, dI, dR)))
  })
}

# Initial state values
initial_state <- c(S = 0.99, I = 0.01, R = 0.0)

# Parameters
parameters <- c(beta = 0.3, gamma = 0.1)

# Time points
times <- seq(0, 1600, by = 1)

# Solve the system of ODEs
start_time <- Sys.time()
for(i in 1:100) {
  output <- ode(y = initial_state, times = times, func = sir_model, parms = parameters)
}
end_time <- Sys.time()
r_runtime <- end_time - start_time

print("Runtime for deSolve in R:")

[1] "Runtime for deSolve in R:"

print(r_runtime)

Time difference of 1.270186 secs

import numpy as np
from scipy.integrate import odeint
import time

# Define the SIR model
def sir_model(y, t, beta, gamma):
    S, I, R = y
    dS = -beta * S * I
    dI = beta * S * I - gamma * I
    dR = gamma * I
    return [dS, dI, dR]

# Initial state values
initial_state = [0.99, 0.01, 0.0]

# Parameters
beta = 0.3
gamma = 0.1

# Time points
times = np.linspace(0, 1600, 1601)

# Solve the system of ODEs
start_time = time.time()
for i in range(0, 100):
  output = odeint(sir_model, initial_state, times, args=(beta, gamma))
end_time = time.time()
python_runtime = end_time - start_time

print("Runtime for scipy.integrate.odeint in Python:")

Runtime for scipy.integrate.odeint in Python:

print(python_runtime)

0.035394906997680664