目标
建立磷酸缓冲液 (1M) 的理论滴定曲线。
我提供了一个完全可重现且独立的示例（我的失败^.^）。
模型方程
磷酸的酸碱平衡方程为：
模型实现
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)
}
请注意，convert.fun可以尝试不同的方法来强制浓度为正值。
返回值是模型方程的向量，等于零（对吗？）。
迭代
我希望解决系统中所有可能的 Na+ 浓度，最多 3 个当量“体积”。
我设置了唤醒最低条件的初始条件：[Na]=0.
然后解决它nleqslv并使用结果来“播种”下一次迭代。
它似乎工作得很好：
但是，仔细观察，问题将变得显而易见。
但是，在那之前，一些代码！
设置初始条件和结果矩阵：
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
迭代：
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)
请注意，convert.fun现在abs(x)是（这样可以吗？）。
问题
最后一个情节的问题在于它的右侧部分被展平了。
这个问题在下面的情节中更加明显：
红色曲线应该在顶部结束，紫色曲线在底部。这似乎从 开始发生Na~2，但经过几次迭代后，结果变平（并变得完全恒定）。
精明的可能提示
- 使用
method="Broyden"而不是"Newton". nleqslv的返回消息从“接近零的函数标准”更改为“公差 'xtol' 内的 x 值”。- 我还尝试添加雅可比行列式。这并没有改变结果，但在有问题的地方，我得到了这样的结果：
Chkjac possible error in jacobian[2,1] = 2.7598836276240e+06
Estimated[2,1] = 1.1104869955110e+04
我现在真的没有想法了！并且非常感谢一些帮助或指导。
