[R] odesolve: lsoda vs rk4

About this list Date view Thread view Subject view Author view Attachment view

From: Chris Knight (christopher.knight@plant-sciences.oxford.ac.uk)
Date: Fri 11 Jun 2004 - 01:46:56 EST

Message-id: <1086882415.1328.168.camel@dops7118.plants.ox.ac.uk>

I'm trying to use odesolve for integrating various series of coupled 1st
order differential equations (derived from a system of enzymatic
catalysis and copied below, apologies for the excessively long set of

The thing that confuses me is that, whilst I can run the function rk4:

out <- rk4(y=y,times=times,func=func, parms=parms)

and the results look not unreasonable:

for (i in 2:(dim(out)[2]))plot(out[,1],out[,i], pch=".", xlab="time",

If I try doing the same thing with lsoda:

out <- lsoda(y=y,times=times,func=func, parms=parms, rtol=1e-1, atol=

I run into problems with a series of 'Excessive precision requested'
warnings with no output beyond the first time point.

Fiddling with rtol and atol doesn't seem to do very much.

What is likely to be causing this (I'm guessing the wide range of the
absolute values of the parameters can't be helping), is there anything I
can sensibly do about it and, failing that, can I reasonably take the
rk4 results as being meaningful?

Any help much appreciated,
Thanks in advance,


func <- function(t, y, p)
Ad <- p["p2"]*(p["p1"]*y["A"]*y["D"])/(p["p2"]+p["p3"]) +
p["p6"]*(p["p4"]*y["B"]*p["p10"])/(p["p5"]+p["p6"]) -
Bd <- p["p3"]*(p["p1"]*y["A"]*y["D"])/(p["p2"]+p["p3"]) +
p["p5"]*(p["p4"]*y["B"]*p["p10"])/(p["p5"]+p["p6"]) -
Cd <- (p["p1"]+p["p7"])*y["A"]*y["D"] -
Dd <-p["p9"]*y["C"] - p["p7"]*y["A"]*y["D"]
list(c(Ad, Bd, Cd, Dd))

parms<-c(p1=4.8e5, p2=1.25, p3=1.3, p4=1e6, p5=1, p6=1.25, p7=1e6,
p8=16, p9=0.35, p10=0.235e-6)
y<-c(A=2.5e-6,B=2.5e-6, C=1.7e-6, D=0.57e-6)
times <- c(0.05 * (0:999))

Dr. Christopher Knight                               Tel:+44 1865 275111
Dept. Plant Sciences                                     +44 1865 275790
South Parks Road
Oxford     OX1 3RB                                   Fax:+44 1865 275074
`  . , ,><(((> 

______________________________________________ R-help@stat.math.ethz.ch mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

About this list Date view Thread view Subject view Author view Attachment view

This archive was generated by hypermail 2.1.3 : Wed 30 Jun 2004 - 23:05:05 EST