R nleqslv 困难 - 解决酸碱缓冲液中的 pH 值
R nleqslv difficulties - solving for pH in an acid-base buffer
目标
建立磷酸缓冲液 (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 可以尝试不同的方法来强制浓度为正值。
return 值是模型方程的向量，等于零（对吗？）。
迭代
我希望求解所有可能的 Na+ 浓度的系统，最多 3 个当量“体积”。
我设置了唤醒最低的初始条件：[Na]=0.
然后用 nleqslv 求解并使用结果为下一次迭代“播种”。
它似乎运行良好：
但是，仔细检查后，问题就会变得很明显。
But, before that, some code!
设置初始条件和结果矩阵：
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 的 return 消息从“函数标准接近零”更改为“x 值在公差范围内 'xtol'”。- 我也试过添加雅可比行列式。那并没有改变结果，但是在有问题的地方我得到了这样的东西：
Chkjac possible error in jacobian[2,1] = 2.7598836276240e+06
Estimated[2,1] = 1.1104869955110e+04
我现在真是没脑子了！非常感谢您的帮助或指导。
您应该始终测试 nleqslv 的终止代码以确定是否已找到解决方案。并以某种方式显示终止代码 and/or 消息 nleqslv returns。您会发现在某些情况下找不到更好的点。因此任何结果都是无效和无用的。
您为 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... 常量的非常小的值是问题的原因。检查系统是否存在错误或尝试更改变量的测量单位。
解决方案详情
查找详细信息和完整代码 at this repo。
数值方法奏效了，提供的分析答案at chemistry stackexchange 恰到好处:)
遗憾的是，它与 Julia Martín 等人 (DOI 10.20431/2349-0403.0409002) 的实验数据不符。也许我会 post 在 chemistry stackexchange 上问一个关于它的问题。
感谢所有帮助过的人 <3
最后，数值模拟的重要图：
目标
建立磷酸缓冲液 (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 可以尝试不同的方法来强制浓度为正值。
return 值是模型方程的向量，等于零（对吗？）。
迭代
我希望求解所有可能的 Na+ 浓度的系统，最多 3 个当量“体积”。
我设置了唤醒最低的初始条件：[Na]=0.
然后用 nleqslv 求解并使用结果为下一次迭代“播种”。
它似乎运行良好：
但是，仔细检查后，问题就会变得很明显。
But, before that, some code!
设置初始条件和结果矩阵：
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的 return 消息从“函数标准接近零”更改为“x 值在公差范围内 'xtol'”。- 我也试过添加雅可比行列式。那并没有改变结果，但是在有问题的地方我得到了这样的东西：
Chkjac possible error in jacobian[2,1] = 2.7598836276240e+06
Estimated[2,1] = 1.1104869955110e+04
我现在真是没脑子了！非常感谢您的帮助或指导。
您应该始终测试 nleqslv 的终止代码以确定是否已找到解决方案。并以某种方式显示终止代码 and/or 消息 nleqslv returns。您会发现在某些情况下找不到更好的点。因此任何结果都是无效和无用的。
您为 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... 常量的非常小的值是问题的原因。检查系统是否存在错误或尝试更改变量的测量单位。
解决方案详情
查找详细信息和完整代码 at this repo。
数值方法奏效了，提供的分析答案at chemistry stackexchange 恰到好处:)
遗憾的是，它与 Julia Martín 等人 (DOI 10.20431/2349-0403.0409002) 的实验数据不符。也许我会 post 在 chemistry stackexchange 上问一个关于它的问题。
感谢所有帮助过的人 <3
最后，数值模拟的重要图：