ode(y, times, func, parms,
method=c("lsoda","lsode","lsodes","lsodar","vode","daspk", "euler", "rk4",
"ode23", "ode45"), ...)y has a name attribute, the names will be used to label the output matrix.times must be the initial time If func is an R
funcrkMethod.y plus the number of "global" values returned
in the second element of the return from func, plus an additional column (the first) for the time value.
There will be one row for each element in times unless the integrator returns with an unrecoverable error.
If y has a names attribute, it will be used to label the columns of the output value.
The output will have the attributes istate, and rstate, two vectors with several useful elements.
The first element of istate returns the conditions under which the last call to the integrator returned. Normal is istate = 2.
If verbose = TRUE, the settings of istate and rstate will be written to the screen. See the help for the selected integrator for details.ode.band for solving models with a banded Jacobian
ode.1D for integrating 1-D models
ode.2D for integrating 2-D models
aquaphy, ccl4model, where ode is used lsoda, lsode, lsodes, lsodar, vode, daspk,
rk, rkMethod
LVmod <- function(Time,State,Pars) {
with(as.list(c(State,Pars)),
{ Ingestion <- rIng * Prey*Predator GrowthPrey <- rGrow * Prey*(1-Prey/K) MortPredator <- rMort * Predator
dPrey <- GrowthPrey - Ingestion dPredator <- Ingestion*assEff -MortPredator
return(list(c( dPrey, dPredator)))
}) }
pars <- c(rIng =0.2, # /day, rate of ingestion rGrow =1.0, # /day, growth rate of prey rMort =0.2 , # /day, mortality rate of predator assEff =0.5, # -, assimilation efficiency K =10 ) # mmol/m3, carrying capacity
yini <- c(Prey=1,Predator=2) times <- seq(0,200,by=1) out <- as.data.frame(lsoda(func= LVmod, y=yini, parms=pars, times=times))
matplot(out$time,out[,2:3],type="l",xlab="time",ylab="Conc", main="Lotka-Volterra",lwd=2) legend("topright",c("prey", "predator"),col=1:2, lty=1:2)
######################################### ## Example2: Resource-producer-consumer Lotka-Volterra model #########################################
## Note: ## 1. parameter and state variable names made ## accessible via "with" statement ## 2. function sigimp passed as an argument (input) to model ## (see also lsoda and rk examples)
lvmodel <- function(t, x, parms, input) {
with(as.list(c(parms,x)), { import <- input(t) dS <- import - b*S*P + g*K #substrate dP <- c*S*P - d*K*P #producer dK <- e*P*K - f*K #consumer res<-c(dS, dP, dK) list(res) }) } ## The parameters parms <- c(b=0.0, c=0.1, d=0.1, e=0.1, f=0.1, g=0.0)
## vector of timesteps times <- seq(0, 100, length=101) ## external signal with rectangle impulse signal <- as.data.frame(list(times = times, import = rep(0,length(times)))) signal$import[signal$times >= 10 & signal$times <=11] <-="" 0.2="" sigimp="" approxfun(signal$times,="" signal$import,="" rule="2)" ##="" start="" values="" for="" steady="" state="" xstart="" c(s="1," p="1," k="1)" solve="" model="" out="" as.data.frame(ode(y="xstart,times=" times,="" func="lvmodel," parms,="" input="sigimp))" plot(out$p,out$k,type="l" ,lwd="2,xlab="producer",ylab="consumer")