R中使用nleqslv求解1M磷酸缓冲液滴定曲线的迭代收敛异常

构建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){
  # 施加仅取正值的约束
  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用于尝试不同方式强制浓度为正值，返回值为令模型方程右端为0的向量（请确认该实现是否正确？）

迭代逻辑 #

需要对最高到3倍当量点的所有Na+浓度序列求解该方程组，初始条件适配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){                   # 第一次迭代
    ci <- ci.start              # 采用设定的初始条件
  } else {                      # 其余迭代
    ci <- result.m[i-1, 2:6]    # 采用前一次的解作为初值
  }
  
  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))
  
} # 循环结束

result.df <- as.data.frame(result.m)

此处convert.fun采用abs(x)实现（请确认该实现是否合理？）

最终生成的曲线右侧被拉平，在以下物种浓度分布图中问题更明显：

物种浓度分布示意图

理论上红色曲线（HA）应最终趋近于0、紫色曲线（A）应趋近于1，二者在Na浓度接近2时开始出现该趋势，但迭代数次后结果完全平坦化、变为恒定值。

• 采用method="Broyden"替代"Newton"时问题更严重
• nleqslv的返回信息从“函数准则接近0”变为“x值处于公差'xtol'范围内”
• 尝试提供自定义雅可比矩阵后结果未发生变化，但在问题出现点得到如下提示：

Chkjac  possible error in jacobian[2,1] =  2.7598836276240e+06
                             Estimated[2,1] =  1.1104869955110e+04

目前暂无解决思路，恳请相关指导与帮助。

内容的提问来源于stack exchange，提问作者Naiky