问题描述

目标

建立磷酸缓冲液 (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 的返回消息从“函数标准接近零”更改为“x 值在容差 'xtol' 内”。
我也尝试添加雅可比行列式。这并没有改变结果，但在有问题的地方，我得到了这样的东西：

    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... 常量的非常小的值是问题的原因。检查系统是否存在可能的错误或尝试更改变量的度量单位。

解决方案详情

查找详细信息和完整代码 at this repo。

数值方法奏效了，提供的解析答案at chemistry stackexchange 很高兴：)

遗憾的是，它与 Julia Martín et al (DOI 10.20431/2349-0403.0409002) 的实验数据不符。也许我会在 chemistry stackexchange 上发布一个关于它的问题。

感谢所有帮助过我的人

最后，数值模拟的重要图：

chemistry nonlinear-equation nonlinear-optimization numerical-methods r r

R nleqslv 困难 - 在酸碱缓冲液中求解 pH