deSolve ODE Output Question
On May 20, 2015, at 6:45 AM, walke554 wrote:
Hello,
I am working on a simple ODE problem with the deSolve package, and I was
hoping that someone could answer a question about how the deSolve package
does integration.
Here is my program:
#The function
STDMod<-function(t,y,p){
IH = y[1];
IL = y[2];
with(as.list(p), {
dIH.dt = Beta[1,1]*(nH-IH)*IH + Beta[1,2]*(nH-IH)*IL - gamma*IH;
dIL.dt = Beta[2,1]*(nL-IL)*IH + Beta[2,2]*(nL-IL)*IL - gamma*IL;
return(list(c(dIH.dt,dIL.dt)));
})
}
#giving the parameters
Beta = matrix(data=c(10, 0.1, 0.1, 1.0), ncol=2, nrow=2)
nH = 0.2
nL = 0.8
IH0 = 1e-5
IL0 = 0
gamma = 1
p = list(Beta=Beta, gamma=gamma, nH=nH, nL=nL)
y0 = c(IH0, IL0)
#Running the ode integrator
steps= 10;
t = seq(from=0, to=30, by=.01);
out = ode(y=y0, times=t, func=STDMod, parms=p);
My understanding is that the 'out' matrix would be the values of the STDmod
function at each time step, given the initial values of the state parameters
IH and IL. However, for the very first time step, I am getting a value
different than what the derivative and the initial values add up to.
My reading of the returned value from a call to ode is that it is not the values returned by the 'func'-function which only contain the derivative vectors, but rather the integrated values. That is what you are seeing in any case. You need to remember that integrators calculate the derivative and then further multiply by the "dt" increment which in your case will make the incremental values one-hundredth of the value you calculated by hand below. (Basic calculus.)
My output: <http://r.789695.n4.nabble.com/file/n4707448/ROutput.png> but 'IH0' + 'the value of dIH.dt = Beta[1,1]*(nH-IH)*IH + Beta[1,2]*(nH-IH)*IL - gamma*IH for timestep one' = 1.999e-5. Is there something that I am missing? I am hoping to use this for more complicated derivatives, but I want to make sure I am using the package correctly first. Thanks! View this message in context: http://r.789695.n4.nabble.com/deSolve-ODE-Output-Question-tp4707448.html Sent from the R help mailing list archive at Nabble.com.
I normally also respond to the list and the questioner but have gotten an annoying number of bounces from Nabble users so am not responding directly.
David Winsemius Alameda, CA, USA