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 at 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 at 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.