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

[Peter Dalgaard]

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