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
parameters).
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:
out<-as.data.frame(out)
par(mfrow=c(4,1))
for (i in 2:(dim(out)[2]))plot(out[,1],out[,i], pch=".",
xlab="time",
ylab=names(out)[i])
If I try doing the same thing with lsoda:
out <- lsoda(y=y,times=times,func=func, parms=parms, rtol=1e-1, atol1e-1)
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,
Chris
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"])
-
p["p1"]*y["A"]*y["C"]
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"])
-
p["p4"]*y["B"]*p["p10"]
Cd <- (p["p1"]+p["p7"])*y["A"]*y["D"]
-
p["p1"]*y["A"]*y["C"]-p["p9"]*y["C"]
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
` ?? . , ,><(((??>
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Chris Knight <christopher.knight at plant-sciences.oxford.ac.uk> writes:> 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 > parameters). > > The thing that confuses me is that, whilst I can run the function rk4:.....> 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......> 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))I'm not going to dig deeply into your problem formulation, but my hunch is that you need to rescale your problem so that you get rid of the very small and large numbers in the parameters and starting values. Extreme values there tend to confuse convergence criteria. (The difference between the two solvers is probably exactly that lsoda has internal loops and convergence criteria, whereas rk4 just chugs along in fixed time steps.) -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
lsoda doesn't pass along the names attribute of the state vector, y. If
you want to use the names of the state vector in your code, you need to
reassign it inside your ode function. I should either fix this in the
code for lsoda, or at least document it! For example, I can run the
following modification of your code:
func2 <- function(t, y, p)
{
names(y) <- c("A","B","C","D") ###
Here put the names attribute back
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"])
-
p["p1"]*y["A"]*y["C"]
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"])
-
p["p4"]*y["B"]*p["p10"]
Cd <-
(p["p1"]+p["p7"])*y["A"]*y["D"] -
p["p1"]*y["A"]*y["C"]-p["p9"]*y["C"]
Dd <-p["p9"]*y["C"] -
p["p7"]*y["A"]*y["D"]
list(c(Ad, Bd, Cd, Dd))
}
out2 <- lsoda(y, times, func2, parms)
> out2[c(1:10,1000),]
time A B C D
[1,] 0.00 2.500000e-06 2.500000e-06 1.700000e-06 5.700000e-07
[2,] 0.05 2.432561e-06 2.500379e-06 1.671689e-06 5.312515e-07
[3,] 0.10 2.368078e-06 2.499286e-06 1.639363e-06 4.980014e-07
[4,] 0.15 2.306621e-06 2.496841e-06 1.603709e-06 4.697523e-07
[5,] 0.20 2.248381e-06 2.493370e-06 1.566611e-06 4.451401e-07
[6,] 0.25 2.193360e-06 2.488923e-06 1.528312e-06 4.239702e-07
[7,] 0.30 2.141477e-06 2.483726e-06 1.489881e-06 4.053217e-07
[8,] 0.35 2.092710e-06 2.477833e-06 1.451555e-06 3.889875e-07
[9,] 0.40 2.046809e-06 2.471378e-06 1.413729e-06 3.744582e-07
[10,] 0.45 2.003632e-06 2.464438e-06 1.376627e-06 3.614428e-07
[11,] 49.95 2.675890e-06 5.563057e-08 1.595362e-09 -7.457989e-11
In the event the results of lsoda and rk4 differ, I would trust the
lsoda result over that of rk4 (see the warnings in the help file for
rk4).
On June 10, 2004, Chris Knight wrote:
[snip ...]
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
parameters).
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:
out<-as.data.frame(out)
par(mfrow=c(4,1))
for (i in 2:(dim(out)[2]))plot(out[,1],out[,i], pch=".",
xlab="time",
ylab=names(out)[i])
If I try doing the same thing with lsoda:
out <- lsoda(y=y,times=times,func=func, parms=parms, rtol=1e-1, atol1e-1)
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,
Chris
[code deleted...]
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Dr. Christopher Knight Tel:+44 1865 275111
Dept. Plant Sciences +44 1865 275790
South Parks Road
Oxford OX1 3RB Fax:+44 1865 275074
` ?? . , ,><(((??>
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
R. Woodrow Setzer, Jr. Phone: (919) 541-0128
Experimental Toxicology Division Fax: (919) 541-4284
Pharmacokinetics Branch
NHEERL B143-01; US EPA; RTP, NC 27711