Hello there,
Suppose you want to solve the following system of ODE's (a simple
Lotka-Volterra predator prey model)
dP/dt = beta*P*V - mu*P
dV/dt = r*V - beta*P*V
where P and V are the numbers of predators and prey. Now, this is
easy to do, but suppose you have a system of equations like this,
dP1/dt = beta1*P1*V1 - mu1*P1
dP2/dt = beta2*P2*V2 - mu2*P2
dV1/dt = r1*V1 - beta1*P1*V1
dV2/dt = r1*V2 - beta1*P2*V2
System 1 and system 2 are independent but that doesn't need to be the
case. Now can you specify this system in lsoda() as:
dP/dt = beta*P*V - mu*P
dV/dt = r*V - beta*P*V
but now the initial state variables are 1x2 vectors and the
parameters beta, mu and r are also vectors of size 1x2:
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) {
y=matrix(y,nc=2,byrow=T)
beta=parms[,1]; mu = parms[,2], r = parms[,3]
dP/dt = beta*P*V - mu*P
dV/dt = r*V - beta*P*V
list(c(dP/dt,dV/dt)
}
I tried this but it does not give me the right answer. lsoda() gives
warnings and produces weird results..
thanks, Jorge
Jorge, If you'll send me details of the error messages, I'll see what I can do to help. I notice in your posting there is a missing ')' in the next to last line of model, but no error message; I don't suppose that could have anything to do with it? (send the details to my work email: setzer.woodrow AT epa DOT gov). In general, I'm willing to help people trying to use lsoda, but it is generally better to email me directly (as the package author, as recommended in the posting guide) rather than through r-help. Woody Setzer Woodrow Setzer National Center for Computational Toxicology US Environmental Protection Agency Research Triangle Park, NC 27711
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