# [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
# 6.102746508205e-07
#\$objective
# -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