From: Peter Dalgaard <p.dalgaard_at_biostat.ku.dk>

Date: Sat, 12 Apr 2008 11:13:12 +0200

$root

[1] -1.466662866313071e-12

$minimum

[1] -1.466662866319826e-12

Date: Sat, 12 Apr 2008 11:13:12 +0200

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 ______________________________________________ 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 - 09:19:47 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 - 09:30:27 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.
*