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

From: Woodrow Setzer <wsetzer_at_mindspring.com>
Date: Thu 27 Oct 2005 - 00:18:54 EST


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



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Thu Oct 27 00:35:01 2005

This archive was generated by hypermail 2.1.8 : Thu 27 Oct 2005 - 03:32:08 EST