r - R nleqslv 困难 - 解决酸碱缓冲液中的 pH 值
问题描述
目标
建立磷酸缓冲液 (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
我现在真的没有想法了!并且非常感谢一些帮助或指导。
解决方案
您应该始终测试终止代码nleqslv
以确定是否找到了解决方案。并以某种方式显示终止代码和/或消息nleqslv
返回。你会看到在某些情况下没有找到更好的点。因此任何结果都是无效和无用的。
您使用了太多的值,Na.seq
以至于不可能穿过树木的木头。
我建议从一组非常有限的值开始Na.seq
。就像是
Na.seq <- seq(from=0,to=3*P_ca,by=P_ca/10)
这也是为了在结果中包含终止代码
varnames <- c("Na", "H", "H3A", "H2A", "HA", "A", "termcd")
result.m <- matrix(ncol = length(varnames), nrow = length(Na.seq))
然后将迭代循环更改为此
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
}
iter.trace <- 1
cat("Iteration ",i,"\n\n")
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,trace=iter.trace))
cat("\n\n ",result$message,"\n\n")
result$x <- convert.fun(result$x)
result.m[i,2:6] <- result$x
result.m[i,7] <- result$termcd
stopifnot(all(result$x >= 0))
} # END LOOP
并开始分析输出以找出问题所在和位置。
附录
我有理由确定解决问题的困难(部分)是由数值困难引起的。通过上述修改,我将Ka.1
、Ka.2
、Ka.3
和的值更改Kw
为
Ka.1 <- 7.1 * 10^-1
Ka.2 <- 6.3 * 10^-3
Ka.3 <- 4.5 * 10^-3
Kw <- 10^-3
然后找到解决方案没有问题(所有终止代码都是1)。我怀疑K...
常量的非常小的值是问题的原因。检查系统是否存在可能的错误或尝试更改变量的测量单位。
推荐阅读
- node.js - 如何在nodejs中使用'select if'函数和sequelize
- mysql - 加密的电子邮件 ID 不区分大小写 MySQL
- sinch - 如何在 Sinch 中结束“群呼”
- android - Why do we need LiveData and ViewModel
- lua - Lua中的“$”符号是什么意思?
- r - 试图了解附加在 R 中的使用
- python - Python 3.7 - ImportError:绘图需要matplotlib
- android - 从 JSON 显示部分数据的困难
- angular - Angular 6服务器端渲染与角度通用给出错误区域未定义
- multithreading - GNU Parallel:限制作业在使用 -n0 用于常量/固定参数时不起作用