[R] uniroot problem

[Andrew Simmons]

In this line:

for (i in 1:100) {
   ATOT[i] <- uniroot(ASTA,c(-1e3,1e3),tol=1e-20,maxiter=1e4,
                      SID=SID,KA=KA[i],H=10^-pH[1])$root
}

you've used variable pH, but have no definition for this variable.


Also, I can't tell which location has a mistake, but these aren't equivalent:

ASTA <- H + SID - kw/H - ATOT*KA/KA+H
ATOT <-  (H+SID-kw/H)*(KA+H)/KA

I'm assuming the first one is incorrect, it seems redundant to
multiply then divide by KA, so change ASTA to:

ASTA <- function(SID,H,ATOT,KA) H + SID - kw/H - ATOT*KA/(KA+H)

I tried something like this:

ASTA <- function(SID,H,ATOT,KA) H + SID - kw/H - ATOT*KA/(KA+H)

SID <- 0.01199
H <- 10^-4.57
kw <- 1e-14

KA <- 10^-seq(1, 10, length.out = 100)

ATOT1 <- (H + SID - kw/H) * (KA + H)/KA

ATOT2 <- vapply(KA, function(ka) {
    uniroot(ASTA, c(-10000, 10000), SID = SID, H = H, KA = ka, tol =
1e-20, maxiter = 10000)$root
}, 0)


ASTA(SID, H, ATOT1, KA)
ASTA(SID, H, ATOT2, KA)

and it seems to work for me. The biggest absolute diff is ~1e-13

On Fri, Sep 16, 2022 at 2:05 AM Troels Ring <tring using gvdnet.dk> wrote:
>
> Dear friends - I have a problem with root finding using uniroot as goes:
>
> I have the function
>
> ASTA <- function(SID,H,ATOT,KA) {H + SID - kw/H - ATOT*KA/KA+H }
>
> I know SID = 0.01199 - and I know H = 10^-4.57 - and I know kw = 1e-14
>
> Therefore, as I make KA = 10^-seq(1,10,length=100), by letting ASTA = 0,
> I can find ATOT for each KA as
>
> ATOT <-  (H+SID-kw/H)*(KA+H)/KA
>
>
> and
>
>    summary(ATOT) =    Min.  1st Qu.   Median     Mean  3rd Qu. Max.
>                                       0.012    0.013    0.115 171.263
> 18.278 3234.407
>
> However, if I make use of uniroot by:
>
> KA <- seq(1,10,length=100)
> KA <- 10^-KA
>
> ATOT <- c()
>
> for (i in 1:100) {
>    ATOT[i] <- uniroot(ASTA,c(-1e3,1e3),tol=1e-20,maxiter=1e4,
>                       SID=SID,KA=KA[i],H=10^-pH[1])$root
> }
>
> I get 100 identical ATOT values of 0.01204383 and no warnings.
>
> I'm on Windows 10, using slightly old
>
> $version.string
> [1] "R version 4.0.5 (2021-03-31)"
>
> Best wishes
>
> Troels Ring, MD
> Aalborg, Denmark
>
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.