# [R] Re: using LSODA in R

Stephen Ellner spe2 at cornell.edu
Thu Nov 6 15:50:55 CET 2003

```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; c=y; r=y; b=y;
fcn=bmaxc*n/(Kc+n);
fbc=bmaxb*c/(Kb+c);
dy=delta*(ni-n)-fcn*c;
dy=fcn*c-fbc*b/epsilon-delta*c;
dy=fbc*r-(delta+m+lambda)*r;
dy=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

```