[R] solving ODE's in matrix form with lsoda()

[Woodrow Setzer]

Just a followup.  I suppose you meant something like this:

library(odesolve)

y <- c(10, 20, 10, 20)
parms <- matrix(c(0.05, 0.1, 0.2, 0.05, 0.1, 0.2), nc=3, byrow=T)
 model <- function(times, y, parms) {
   P <- y[1:2]
   V <- y[3:4]
   beta <- parms[,1]
   mu <- parms[,2]
   r <- parms[,3]
   dPdT <- beta*P*V - mu*P
   dVdt <- r*V - beta*P*V
   list(c(dPdT, dVdt))
 }

out <- lsoda(y, times=(0:100), model, parms)

which gives the following output (for the first 8 time units) and no errors or warnings:

> source("C:/home/testode.R")
> out
       time          1            2           3            4
  [1,]    0 10.0000000 20.000000000 10.00000000 2.000000e+01
  [2,]    1 13.7603737 33.615668238  6.72208454 6.079882e+00
  [3,]    2 16.1560654 35.516702438  3.85780113 1.275331e+00
  [4,]    3 16.8730079 33.195852918  2.05089685 2.778086e-01
  [5,]    4 16.4682671 30.260712324  1.08485347 6.942984e-02
  [6,]    5 15.5186874 27.434895458  0.59474861 2.006117e-02
  [7,]    6 14.3655789 24.839030740  0.34399530 6.638867e-03
  [8,]    7 13.1757074 22.479990536  0.21106101 2.486824e-03
  [9,]    8 12.0240882 20.342411249  0.13733214 1.042244e-03

Perhaps you have some errors in your model code?

Woody

Woodrow Setzer
National Center for Computational Toxicology
US Environmental Protection Agency
Research Triangle Park, NC 27711