[R] Re: using LSODA in R

[Stephen Ellner]

I think an additional problem is that the derivatives function 
must be 'told' where to 'look' in y for the values of R, C, and P 
that are used to compute Rprime,Cprime, Pprime - it has no way of
'knowing' that the model's state vector y consists of
the variables (R,C,P). Appended below is an example of code
that works to solve a similar system of equations. Perhaps 
it might be useful to add something like that to the 
?lsoda page?  

> function(times,y,p)
> {
> Rprime <-
> (R*(1-R))-((xc*yc*C*R)/(R+R0))-((w*xp*ypr*P*R)/(R02+((1-w)*C)+(w*R)))
> Cprime <-
> (-1*(xc*C)*(1-(yc*R)/(R+R0)))-(((1-w)*xp*ypc*P*C)/((w*R)+((1-w)*C)+C0))
> Pprime <-
> 
>(-1*P)-(((1-w)*xp*ypc*C*P)/((w*R)+((1-w)*C)+C0))+((w*xp*ypr*P*R)/((w*R)+((1-w)*C)+R02))
>list(c(Rprime, Cprime, Pprime))
>}

require(odesolve);

# Set parameter values
bmaxc=3.3; bmaxb=2.25; Kc=4.3; Kb=15; ni=80;
m=0.05; lambda=0.4; epsilon=0.25; delta=0.95;

f2k=function(t,y,parms){
        dy=rep(0,4); 
        n=y[1]; c=y[2]; r=y[3]; b=y[4];
        fcn=bmaxc*n/(Kc+n);  
        fbc=bmaxb*c/(Kb+c);
        dy[1]=delta*(ni-n)-fcn*c;
        dy[2]=fcn*c-fbc*b/epsilon-delta*c;
        dy[3]=fbc*r-(delta+m+lambda)*r;
        dy[4]=fbc*r-(delta+m)*b;
        list(dy);
}   

y0=c(1,70,.01,.01); times=(0:240)/4; parms=0; 
out=lsoda(y0, times, f2k, parms, rtol=1e-4, atol=1e-4);
matplot(times,cbind(out[,3]/5,out[,5]),type="l",col=c("green","red"), lty=1,
xlab="Time (d)", ylab="Prey (green) and Predators (red)"); 


Stephen P. Ellner (spe2 at cornell.edu)
Department of Ecology and Evolutionary Biology
Corson Hall, Cornell University, Ithaca NY 14853-2701
Phone (607) 254-4221    FAX (607) 255-8088