[R] R and Excel disagreement - Goal Seek versus uniroot

From: Troels Ring <tring_at_gvdnet.dk>
Date: Sat, 12 Apr 2008 10:01:51 +0200


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_at_gvdnet.dk


	[[alternative HTML version deleted]]

______________________________________________
R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Received on Sat 12 Apr 2008 - 08:05:15 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sat 12 Apr 2008 - 11:30:28 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive