Goal
Build a theoretical titration curve for the phosphoric acid buffer (1M).
I provide a fully reproducible and self-contained example (of my failures ^.^).
Model equations
Acid-base equilibrium equations for phosphoric acid are:
Model implementation
Ka.1 <- 7.1 * 10^-3
Ka.2 <- 6.3 * 10^-8
Ka.3 <- 4.5 * 10^-13
Kw <- 10^-14
balance <- function(vars, Na_ca, P_ca, convert.fun=function(x) x){
# Apply positive only constraint
vars <- convert.fun(vars)
H <- vars[1]
H3A <- vars[2]
H2A <- vars[3]
HA <- vars[4]
A <- vars[5]
Na <- convert.fun(Na_ca)
eq.system <- c(H3A + H2A + HA + A - P_ca,
H + Na - Kw/H - H2A - 2*HA - 3*A,
H * H2A / Ka.1 - H3A,
H * HA / Ka.2 - H2A,
H * A / Ka.3 - HA
)
return(eq.system)
}
Notice that convert.fun is there to try different ways of forcing positive values on concentrations.
The return value is a vector of the model's equations, equated to zero (is this right?).
Iteration
I wished to solve the system for all possible Na+ concentrations, up to 3 equivalence "volumes".
I set initial conditions that woked for the lowest one: [Na]=0.
Then solved it with nleqslv and used the result to "seed" the next iteration.
And it seemed to work nicely:
But, on close inspection, the issues will become obvious.
But, before that, some code!
Setup initial conditions and results matrix:
P_ca <- 1
ci.start <- c(H=10^-1, H3A=0.9, H2A=0.1, HA=0.1, A=0.1)
Na.seq <- seq(from=0,to=3*P_ca,by=P_ca/1000)
varnames <- c("Na", "H", "H3A", "H2A", "HA", "A")
result.m <- matrix(ncol = length(varnames), nrow = length(Na.seq))
colnames(result.m) <- varnames
result.m[,1] <- Na.seq
Iteration:
convert.fun <- function(x) abs(x)
for(i in 1:length(Na.seq)){
Na_ca <- result.m[i,1]
if(i == 1){ # Si es la primera iteración,
ci <- ci.start # usar los valores "start" como C.I.
} else { # Si no,
ci <- result.m[i-1, 2:6] # usar los valores de la solución anterior
}
result <- nleqslv::nleqslv(x = ci,
fn = balance,
Na = Na_ca, P = P_ca,
convert.fun = convert.fun,
method="Newton", #method="Broyden",
global="dbldog",
control = list(allowSingular=TRUE,
maxit=1000))
result$x <- convert.fun(result$x)
result.m[i,2:6] <- result$x
stopifnot(all(result$x >= 0))
} # END LOOP
result.df <- as.data.frame(result.m)
Notice that convert.fun is now abs(x) (is this ok?).
The problem
The problem with the last plot is that the right part of it is flattened out.
The problem is even more obvious in the following plot:
The red curve is supossed to end up at the top, and the purple one at the bottom. This seems to start happening at Na~2, but after a few more iterations, the result flattens out (and becomes exactly constant).
Possible hints for the savvy
- The problem is a bit worse using
method="Broyden"instead of"Newton". nleqslv's return message changes from "Function criterion near zero" to "x-values within tolerance 'xtol'".- I also tried adding a Jacobian. That didnt change the result, but at the problematic point I get something like this:
Chkjac possible error in jacobian[2,1] = 2.7598836276240e+06
Estimated[2,1] = 1.1104869955110e+04
I am now really out of ideas! And would really appreciate some help or guidance.
