R help list subscribers, I am a new user of R. I am attempting to use R to explore a set of equations specifying the dynamics of a three trophic level food chain. I have put together this code for the function that is to be evaluted by LSODA. My equations Rprime, Cprime, and Pprime are meant to describe the actual equation of the derivative. When I run LSODA, I do not get the output that these equations should be giving. Can someone tell me if I have set this function up correctly to use with LSODA when the user is specifying the equation of the derivative or offer some advice for using LSODA in R? An example of how to code for user specified differential equations would be great. 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)) } The above is the function yprime which the documentation for the odesolve says that I may specify. Thanks for any help that anyone can provide. Ivan Kautter _________________________________________________________________ Compare high-speed Internet plans, starting at $26.95. https://broadband.msn.com (Prices may vary by service area.)
Try returning list(c(Rprime,Cprime,Pprime),NULL) -- the first element in the returned list should be a numeric *vector* of the derivatives. Ben On Tue, 4 Nov 2003, Ivan Kautter wrote:> R help list subscribers, > > I am a new user of R. I am attempting to use R to explore a set of > equations specifying the dynamics of a three trophic level food chain. I > have put together this code for the function that is to be evaluted by > LSODA. My equations Rprime, Cprime, and Pprime are meant to describe the > actual equation of the derivative. When I run LSODA, I do not get the > output that these equations should be giving. Can someone tell me if I have > set this function up correctly to use with LSODA when the user is specifying > the equation of the derivative or offer some advice for using LSODA in R? > An example of how to code for user specified differential equations would be > great. > > 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)) > } > > The above is the function yprime which the documentation for the odesolve > says that I may specify. > > Thanks for any help that anyone can provide. > > Ivan Kautter > > _________________________________________________________________ > Compare high-speed Internet plans, starting at $26.95. > https://broadband.msn.com (Prices may vary by service area.) > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-help >-- 620B Bartram Hall bolker at zoo.ufl.edu Zoology Department, University of Florida http://www.zoo.ufl.edu/bolker Box 118525 (ph) 352-392-5697 Gainesville, FL 32611-8525 (fax) 352-392-3704
Thanks for the response, Ben. It's not a case that the list isn't coming out correctly. It's that the numbers that are coming out are not the numbers that these equations should be producing if I have specified the equations correctly in R code for use with LSODA. So the question is more if I have the code right when the user specifies the differential equations for LSODA.>From: Ben Bolker <bolker at zoo.ufl.edu> >Reply-To: bolker at zoo.ufl.edu >To: Ivan Kautter <ivankautter at hotmail.com> >CC: R help list <r-help at stat.math.ethz.ch> >Subject: Re: [R] using LSODA in R >Date: Wed, 5 Nov 2003 08:43:53 -0500 (EST) > > > Try returning list(c(Rprime,Cprime,Pprime),NULL) -- the first element in >the returned list should be a numeric *vector* of the derivatives. > > Ben > >On Tue, 4 Nov 2003, Ivan Kautter wrote: > > > R help list subscribers, > > > > I am a new user of R. I am attempting to use R to explore a set of > > equations specifying the dynamics of a three trophic level food chain. >I > > have put together this code for the function that is to be evaluted by > > LSODA. My equations Rprime, Cprime, and Pprime are meant to describe >the > > actual equation of the derivative. When I run LSODA, I do not get the > > output that these equations should be giving. Can someone tell me if I >have > > set this function up correctly to use with LSODA when the user is >specifying > > the equation of the derivative or offer some advice for using LSODA in >R? > > An example of how to code for user specified differential equations >would be > > great. > > > > 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)) > > } > > > > The above is the function yprime which the documentation for the >odesolve > > says that I may specify. > > > > Thanks for any help that anyone can provide. > > > > Ivan Kautter > > > > _________________________________________________________________ > > Compare high-speed Internet plans, starting at $26.95. > > https://broadband.msn.com (Prices may vary by service area.) > > > > ______________________________________________ > > R-help at stat.math.ethz.ch mailing list > > https://www.stat.math.ethz.ch/mailman/listinfo/r-help > > > >-- >620B Bartram Hall bolker at zoo.ufl.edu >Zoology Department, University of Florida http://www.zoo.ufl.edu/bolker >Box 118525 (ph) 352-392-5697 >Gainesville, FL 32611-8525 (fax) 352-392-3704 >_________________________________________________________________ Compare high-speed Internet plans, starting at $26.95. https://broadband.msn.com (Prices may vary by service area.)
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