Dear friends - occurring in Windows R2.6.2
I am modeling physical chemistry in collaboration with a friend who has
preferred working in Excel. I used uniroot, and find a solution to a two
buffer problem in acid-base chemistry which I believe is physiologically
sensible. Using "goal seek" in Excel my friend found another plausible
root, quite close to zero, and a plot of the function fd below shows
that it may be about OK - but Excel fails to find the plausible root
indicated by uniroot.
fd <- function(H,SID,KA1,KA2,ATOT1,ATOT2){
H^4+H^3*(SID+KA1+KA2)+H^2*(SID*(KA1+KA2)+KA1*KA2-ATOT1*KA1-ATOT2*KA2-kw)+
H*(SID*KA1*KA2-KA1*KA2*(ATOT1+ATOT2)-kw*(KA1+KA2))-kw*KA1*KA2}
KA1 <- 1e-6;KA2 <- 1e-5;ATOT1 <- 0.05; ATOT2 <- 0.03; SID <-
0.05;kw <- 4.4e-14
options(digits=20)
uniroot(fd,c(0,1),tol
.Machine$double.eps,maxiter=100000,KA1=KA1,KA2=KA2,SID=SID,
ATOT1=ATOT1,ATOT2=ATOT2)$root #1.16219184891115e-06
And here is the plot indicating to me that the root I found at 1.16e-06
might be right and also that the root at 0 might be OK too - but less
interesting.
H <- seq(0,2e-6,length=10000)
plot(H,fd(H,SID,KA1,KA2,ATOT1,ATOT2))
abline(v=1.16219184891115e-06)
abline (h=0)
text(9.0e-7,3e-19,"H=1.1622e-6")
However, using optimize another solution is found:
xmin <- optimize(fd,c(0,1),tol
.Machine$double.eps,KA1=KA1,KA2=KA2,SID=SID,
ATOT1=ATOT1,ATOT2=ATOT2)
xmin
#$minimum
#[1] 6.102746508205e-07
#$objective
#[1] -9.72253673562293e-20
I would very much appreciate some help on doing this modeling with very
small numbers in R. Comments on the capability of Excel's goal seek
would also be welcome.
Best wishes
Troels
--
Troels Ring - -
Department of nephrology - -
Aalborg Hospital 9100 Aalborg, Denmark - -
+45 99326629 - -
tring@gvdnet.dk
[[alternative HTML version deleted]]
Peter Dalgaard
2008-Apr-12 09:13 UTC
[R] R and Excel disagreement - Goal Seek versus uniroot
Troels Ring wrote:> Dear friends - occurring in Windows R2.6.2 > I am modeling physical chemistry in collaboration with a friend who has > preferred working in Excel. I used uniroot, and find a solution to a two > buffer problem in acid-base chemistry which I believe is physiologically > sensible. Using "goal seek" in Excel my friend found another plausible > root, quite close to zero, and a plot of the function fd below shows > that it may be about OK - but Excel fails to find the plausible root > indicated by uniroot. > > fd <- function(H,SID,KA1,KA2,ATOT1,ATOT2){ > H^4+H^3*(SID+KA1+KA2)+H^2*(SID*(KA1+KA2)+KA1*KA2-ATOT1*KA1-ATOT2*KA2-kw)+ > H*(SID*KA1*KA2-KA1*KA2*(ATOT1+ATOT2)-kw*(KA1+KA2))-kw*KA1*KA2} > > KA1 <- 1e-6;KA2 <- 1e-5;ATOT1 <- 0.05; ATOT2 <- 0.03; SID <- > 0.05;kw <- 4.4e-14 > options(digits=20) > uniroot(fd,c(0,1),tol > .Machine$double.eps,maxiter=100000,KA1=KA1,KA2=KA2,SID=SID, > ATOT1=ATOT1,ATOT2=ATOT2)$root #1.16219184891115e-06 > > > And here is the plot indicating to me that the root I found at 1.16e-06 > might be right and also that the root at 0 might be OK too - but less > interesting. > > H <- seq(0,2e-6,length=10000) > plot(H,fd(H,SID,KA1,KA2,ATOT1,ATOT2)) > abline(v=1.16219184891115e-06) > abline (h=0) > text(9.0e-7,3e-19,"H=1.1622e-6") > > > However, using optimize another solution is found: > > xmin <- optimize(fd,c(0,1),tol > .Machine$double.eps,KA1=KA1,KA2=KA2,SID=SID, > ATOT1=ATOT1,ATOT2=ATOT2) > xmin > #$minimum > #[1] 6.102746508205e-07 > #$objective > #[1] -9.72253673562293e-20 > > > I would very much appreciate some help on doing this modeling with very > small numbers in R. Comments on the capability of Excel's goal seek > would also be welcome. > > Best wishes > Troels > >A number of points to be made here: (1) uniroot finds a root, optimize a minimum. In this case, your function is fairly close to a quadratic so the plot looks like a parabola with two root and a minimum about midways. (And a good starting point could be to drop off the 3rd and 4th order terms in H and analyze the quadratic using your highschool arithmetic.) (2) Since the value at zero is -4.4e-25 and the function is decreasing there, the root on the left is not zero, but slightly to the left of it (at -1.47e-12, it would seem), so the only root in the interval (0,1) is indeed the one at 1.16e-6. Presumably, Excel is using local linearization around zero and finding the nearest root. (3) you _can_ use optimize (or optim) to find roots, but you need to square the objective function first, and you need rather better starting values, because the squared function is not convex (plus you have the risk of finding an optimum which is not at zero). (4) you have to be very careful with optimizers and root finders if the scale of the objective function and its parameters is not close to 1, because they can (and will) misinterpret the small values as indications of convergence. For uniroot and optimize, you may need to set the tolerance well below .Machine$double.eps (this is floating point, so double.eps is only relevant for values that are not themselves small), and for optim. you can play with control parameters fnscale and parscale, e.g. > uniroot(fd,c(-1e-6,0),tol + 1e-70,maxiter=100000,KA1=KA1,KA2=KA2,SID=SID, + ATOT1=ATOT1,ATOT2=ATOT2) $root [1] -1.466662866313071e-12 $f.root [1] 0 $iter [1] 4 $estim.prec [1] 3.187027620317954e-19 > optimize(fq,c(-1e-6,0),tol + 1e-70,KA1=KA1,KA2=KA2,SID=SID, + ATOT1=ATOT1,ATOT2=ATOT2) $minimum [1] -1.466662866319826e-12 $objective [1] 4.10609522514202e-72 > optim(0, fq, control=list(fnscale=1e-70, parscale=1e-12),method="BFGS", + KA1=KA1,KA2=KA2,SID=SID, + ATOT1=ATOT1,ATOT2=ATOT2) $par [1] -1.466662866312404e-12 $value [1] 4.000708456650843e-74 $counts function gradient 510 20 $convergence [1] 0 $message NULL (5) This is a fourth degree polynomial in H. R has a polyroot() function... -- O__ ---- Peter Dalgaard ?ster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907