Numerical solving Lotka-Volterra ODE in R

Question

Aim: I am trying to numerically solve a Lotka-Volterra ODE in R, using de sde.sim() function in the sde package. I would like to use the sde.sim() function in order to eventually transform this system into an SDE. So initially, I started with a simple ODE system (Lotka-Volterra model) without a noise term.

The Lotka-Volterra ODE system: $$ \left\{\begin{matrix}\frac{dx}{dt}= \alpha x -\beta xy. \\ \frac{dy}{dt}= \delta xy -\gamma y \end{matrix}\right. $$

with initial values for x = 10 and y = 10.

The parameter values for alpha, beta, delta and gamma are 1.1, 0.4, 0.1 and 0.4 respectively (mimicking this example).

Attempt to solve the problem:

library(sde)
d <- expression((1.1 * x[0] - 0.4 * x[0] * x[1]), (0.1 * x[0] * x[1] - 0.4 * x[1]))
s <- expression(0, 0)
X <- sde.sim(X0=c(10,10), T = 10, drift=d, sigma=s) 
plot(X)

However, this does not seem to generate a nice cyclic behavior of the predator and prey population.

Question:

What is going wrong in the attempt to solve this system of ODEs in sde.sim()?

Answer 1

I did not yet succeed to use the sde.sim() function of the SDE package, however, I succeeded to solve the system (with and without noise) using the suggestions of Chris and the diffeqr package in R.

library(plotly)

# Lotka-Volterra Model (SDE without noise)
f <- function(u,p,t) {
  du1 = p[1]*u[1]-p[2]*u[1]*u[2]
  du2 = p[3]*u[1]*u[2]-p[4]*u[2]
  return(c(du1,du2))
}

g <- function(u,p,t) {
  return(c(0*u[1],0*u[2]))
}

u0 = c(10.0, 10.0)
tspan <- list(0.0,100.0)
p = c(1.1, 0.4, 0.1, 0.4)
sol = diffeqr::sde.solve(f,g,u0,tspan,p=p,saveat=0.005)
udf = as.data.frame(sol$u)

plot_ly(udf, x = sol$t, y = ~V1, name = 'Prey', type = 'scatter', mode = 'lines') %>%
  add_trace(y = ~V2, name = 'Predator', mode = 'lines')

So, not the solution to the problem, but a step in the right direction.