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 - -
[EMAIL PROTECTED]
[[alternative HTML version deleted]]
______________________________________________
R-help@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.
