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]]

