2

I have interest in using the R language and environment to numerically solve a system of linear ordinary differential equations. The numerical solver, deSolve, handles this just fine when I write and equation for each equation. For example:

library(package = "deSolve")

logistic.growth <- function(t, state, parms) {
   with(as.list(c(state, parms)), {
   dx <- a*x - b*x*x - c*x*y
   dy <- d*y - e*y*x - f*y*y
   return(list(c(dx, dy)))   
   })
}

parameters <- c(a = 0.5, b = 1, c = 0.75, d = 0.1, e = 1.25, f = 0.25)
init <- c(x = 0.1, y = 0.2)

time.series <- ode(y = init, func = logistic.growth, times = seq(from = 0, to = 10, by = 0.1), parms = parameters)

If I wanted to solve for two equations, I would simply add another equation, etc. That's no problem.

With a few equations this is simple enough and not a burden. My problem arises from wanting to create a large system of linear equations (tens or hundreds) and wanting to use a matrix of linear coefficients for all of the pairwise interactions terms.

What would be the most efficient and elegant ways of doing this?

Chris Moore
  • 211
  • 2
  • 6

1 Answers1

2

I found this question after a couple of months, so I assume a solution was already found. Nevertheless, here an example, just for the record:

library(deSolve)

## matrix multiplication makes this model compact
multi <- function(t, n, parms) {
  with(parms, {
    dn <- r * n  + n * (A %*% n)
    return(list(c(dn)))
  })
}
times <- seq(from = 0, to = 500, by = 0.1)
init  <- c(prey1 = 1, prey2 = 1, pred1 = 2, pred2 = 2)

## parms is a list, containing a vector and a matrix
parms <- list(
  r = c(r1 = 0.1, r2 = 0.1, r3 = -0.1, r4 = -0.1),
  ## only pairwise interactions:
  A = matrix(c(0.0, 0.0, -0.2,  0.0,     # prey 1
               0.0, 0.0,  0.0, -0.1,     # prey 2
               0.2, 0.0,  0.0,  0.0,     # predator 1; eats prey 1
               0.0, 0.1,  0.0,  0.0),    # predator 2; eats prey 2
             nrow = 4, ncol = 4, byrow = TRUE)
)
out <- ode(init, times, multi, parms)
plot(out)

More material about this: http://desolve.r-forge.r-project.org

tpetzoldt
  • 121
  • 5