pH平衡计算器数值不稳定：算法问题与约束实现咨询

背景 #

开发pH平衡计算器时，已将问题转化为包含核心约束H_plus*OH_minus=10^-14的联立方程组：使用fsolve可行但对初始条件敏感且无法施加浓度非负约束；使用Nelder-Mead算法收敛后pKw值偏离14，存在三个技术问题待解决。

1. Nelder-Mead算法表现不佳的原因 #

Nelder-Mead是无导数的单纯形优化算法，失效核心原因是变量尺度差异过大：

• 体系中HA、A⁻、Na⁺浓度在0.1mol/L量级，H⁺、OH⁻在1e-7mol/L量级，两者相差10^6倍。目标函数为各方程平方和，大尺度变量对应的方程（如质量平衡、电荷平衡）的平方误差占总误差主导，算法会优先满足这些大项，对Kw约束这类小尺度项的偏差容忍度极高——即使H⁺*OH⁻偏离1e-14几个数量级，其平方误差在总误差中占比仍可忽略，导致pKw偏离14。
• 收紧容差无效的原因：双精度浮点数固有精度极限约为1e-16，设置xatol=1e-40或fatol=1e-40远超硬件能达到的精度，算法会在达到浮点数精度极限后提前收敛，无法进一步优化小尺度项的偏差。

2. 在fsolve中施加浓度非负约束 #

fsolve本身是无约束根求解器，需通过以下方式间接实现非负约束：

方法1：变量替换（最可靠）

将所有非负变量替换为平方项，确保变量取值天然非负：

import numpy as np
from scipy.optimize import fsolve

def equations(variables, Ka, Kw, initial_concentrations):
    # 替换变量：HA = x1², A_minus=x2², H_plus=x3², OH_minus=x4², ion_Na=x5²
    x1, x2, x3, x4, x5 = variables
    HA = x1**2
    A_minus = x2**2
    H_plus = x3**2
    OH_minus = x4**2
    ion_Na = x5**2

    HA_init, A_minus_init, _, _, ion_Na_init = initial_concentrations
    eq1 = HA + A_minus - HA_init - A_minus_init
    eq2 = H_plus + ion_Na - OH_minus - A_minus
    eq3 = H_plus * OH_minus - Kw
    eq4 = Ka * HA - A_minus * H_plus
    eq5 = ion_Na - ion_Na_init
    return [eq1, eq2, eq3, eq4, eq5]

def equilibrate(initial_concentrations, initial_guess):
    Ka = 1.8e-5
    Kw = 1e-14
    # 初始猜测替换为平方根，避免0的平方根报错
    guess_sqrt = np.sqrt(np.maximum(initial_guess, 1e-10))
    solution_sqrt = fsolve(equations, guess_sqrt, args=(Ka, Kw, initial_concentrations), xtol=1e-12)
    # 还原为原变量
    solution = solution_sqrt**2
    return solution

initial_concentrations = [0.0, 0.1, 1e-7, 1e-7, 0.1]
initial_guess = initial_concentrations
HA, A_minus, H_plus, OH_minus, ion_Na = equilibrate(initial_concentrations, initial_guess)

pH = -np.log10(H_plus)
pKw = -np.log10(H_plus) - np.log10(OH_minus)
print(f"HA: {HA:.8f}")
print(f"A-: {A_minus:.8f}")
print(f"H+: {H_plus:.8f}")
print(f"OH-: {OH_minus:.8f}")
print(f"pH: {pH:.2f}")
print(f"pKw: {pKw:.2f}")

方法2：用带约束的根求解器替代fsolve

使用scipy.optimize.root的trust-constr方法，直接设置变量 bounds：

import numpy as np
from scipy.optimize import root

def equations(variables, Ka, Kw, initial_concentrations):
    HA, A_minus, H_plus, OH_minus, ion_Na = variables
    HA_init, A_minus_init, _, _, ion_Na_init = initial_concentrations
    eq1 = HA + A_minus - HA_init - A_minus_init
    eq2 = H_plus + ion_Na - OH_minus - A_minus
    eq3 = H_plus * OH_minus - Kw
    eq4 = Ka * HA - A_minus * H_plus
    eq5 = ion_Na - ion_Na_init
    return [eq1, eq2, eq3, eq4, eq5]

initial_concentrations = [0.0, 0.1, 1e-7, 1e-7, 0.1]
initial_guess = initial_concentrations
Ka = 1.8e-5
Kw = 1e-14

# 设置变量非负约束
bounds = [(0, None) for _ in initial_guess]
result = root(equations, initial_guess, args=(Ka, Kw, initial_concentrations), 
              method='trust-constr', bounds=bounds, options={'xtol':1e-12})

HA, A_minus, H_plus, OH_minus, ion_Na = result.x
pH = -np.log10(H_plus)
pKw = -np.log10(H_plus) - np.log10(OH_minus)
print(f"HA: {HA:.8f}")
print(f"A-: {A_minus:.8f}")
print(f"H+: {H_plus:.8f}")
print(f"OH-: {OH_minus:.8f}")
print(f"pH: {pH:.2f}")
print(f"pKw: {pKw:.2f}")

3. 更优的求解方案 #

针对酸碱平衡问题，推荐以下两种高效可靠的方案：

方案1：简化方程组+对数变换

先通过惰性离子约束减少变量数，再对乘性方程做对数变换，统一变量尺度：

import numpy as np
from scipy.optimize import minimize

def objective(variables, Ka, Kw, initial_concentrations):
    # 变量简化：ion_Na已知为initial_concentrations[4]，只保留HA, A_minus, H_plus
    HA, A_minus, H_plus = variables
    HA_init, A_minus_init, _, _, ion_Na = initial_concentrations
    OH_minus = Kw / H_plus  # 直接用Kw约束计算OH-，减少变量
    
    # 对数变换优化尺度
    eq1 = np.log(HA + A_minus) - np.log(HA_init + A_minus_init)
    eq2 = np.log(H_plus + ion_Na) - np.log(OH_minus + A_minus)
    eq3 = np.log(Ka * HA) - np.log(A_minus * H_plus)
    return eq1**2 + eq2**2 + eq3**2

initial_concentrations = [0.0, 0.1, 1e-7, 1e-7, 0.1]
initial_guess = [0.005, 0.095, 1e-9]  # 更合理的初始猜测
Ka = 1.8e-5
Kw = 1e-14

# 设置非负约束，H+设下限避免log(0)
bounds = [(0, None), (0, None), (1e-15, None)]
result = minimize(objective, initial_guess, args=(Ka, Kw, initial_concentrations),
                  bounds=bounds, method='L-BFGS-B', options={'ftol':1e-12})

HA, A_minus, H_plus = result.x
OH_minus = Kw / H_plus
pH = -np.log10(H_plus)
pKw = -np.log10(H_plus) - np.log10(OH_minus)
print(f"HA: {HA:.8f}")
print(f"A-: {A_minus:.8f}")
print(f"H+: {H_plus:.8f}")
print(f"OH-: {OH_minus:.8f}")
print(f"pH: {pH:.2f}")
print(f"pKw: {pKw:.2f}")

方案2：专用酸碱平衡求解框架

对于复杂酸碱体系，可使用专门的化学计算库（如pyneqsys），内置针对平衡问题的数值优化策略，自动处理尺度问题和约束：

# 需先安装：pip install pyneqsys
from pyneqsys import ChemicalEquilibriumSystem

# 定义平衡反应
reactions = [
    {'reactants': {'HA': 1}, 'products': {'A-': 1, 'H+': 1}, 'K': 1.8e-5},
    {'reactants': {'H2O': 1}, 'products': {'H+': 1, 'OH-': 1}, 'K': 1e-14}
]
# 初始浓度
init_conc = {'HA': 0.0, 'A-': 0.1, 'Na+': 0.1, 'H+': 1e-7, 'OH-': 1e-7}
# 电荷平衡约束
charge_balance = {'Na+': 1, 'H+': 1, 'A-': -1, 'OH-': -1}

sys = ChemicalEquilibriumSystem(reactions, charge_balance=charge_balance)
sol = sys.solve_equilibrium(init_conc)

print(f"HA: {sol.concentrations['HA']:.8f}")
print(f"A-: {sol.concentrations['A-']:.8f}")
print(f"H+: {sol.concentrations['H+']:.8f}")
print(f"OH-: {sol.concentrations['OH-']:.8f}")
print(f"pH: {sol.pH:.2f}")
print(f"pKw: {sol.pKw:.2f}")

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