我在做一个酸碱平衡计算器。从酸/碱的组成开始,计算器应该根据物理常数得出每种组分的平衡浓度。我把它表述成一组需要同时求解的方程。一个关键的约束是H_plus*OH_minus
必须是10^-14
。我把这个也算在方程式里了。
Fsolve.这个方法有效,但是对初始条件很敏感,我不能设置约束条件(我需要浓度>=0):
import numpy as np
from scipy.optimize import fsolve
def equations(variables, Ka, Kw, initial_concentrations):
HA, A_minus, H_plus, OH_minus, ion_Na = variables
HA_init, A_minus_init, H_plus_init, OH_minus_init, ion_Na_init = initial_concentrations
# Define the equations
eq1 = HA + A_minus - HA_init - A_minus_init # Mass balance on acetate
eq2 = H_plus + ion_Na - OH_minus - A_minus # charge balance
eq3 = H_plus * OH_minus - Kw # Kw for water
eq4 = Ka * HA - A_minus * H_plus # HA acid dissociation
eq5 = ion_Na - ion_Na_init # inert ion
return [eq1, eq2, eq3, eq4, eq5]
def equilibrate(initial_concentrations, initial_guess):
# User-supplied equilibrium constant (Ka) and initial concentration of acetic acid (HA_initial)
Ka = 1.8e-5 # mol/L, acetic acid dissociation constant
Kw = 1e-14 # mol/L, ionic product of water
# Use fsolve from SciPy to find the equilibrium concentrations
solution = fsolve(equations, initial_guess, args=(Ka, Kw, initial_concentrations), xtol=1e-18)
return solution
initial_concentrations = [0.0, 0.1, 1e-7, 1e-7, 0.1]
initial_guess = initial_concentrations
# Extract the concentrations
HA, A_minus, H_plus, OH_minus, ion_Na = equilibrate(initial_concentrations, initial_guess)
# Calculate pH
pH = -np.log10(H_plus)
pKw = -np.log10(H_plus) -np.log10(OH_minus)
# Display the results
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}") # THIS IS THE CONDIRMATORY CHECK - Should be 14.
Nelder-Mead。尽管我设置了严格的容差,但这会收敛,但会给出无意义的值(pKw):
import numpy as np
from scipy.optimize import minimize
def equations(variables, Ka, Kw, initial_concentrations):
HA, A_minus, H_plus, OH_minus, ion_Na = variables
HA_init, A_minus_init, H_plus_init, OH_minus_init, ion_Na_init = initial_concentrations
# Define the equations
eq1 = HA + A_minus - HA_init - A_minus_init # Mass balance on acetate
eq2 = H_plus + ion_Na - OH_minus - A_minus # Charge balance
eq3 = H_plus * OH_minus - Kw # Kw for water
eq4 = Ka * HA - A_minus * H_plus # HA acid dissociation
eq5 = ion_Na - ion_Na_init # Inert ion
return eq1**2 + eq2**2 + eq3**2 + eq4**2 + eq5**2 # Use squared error for the minimization
def equilibrate(initial_concentrations, initial_guess, Ka, Kw):
# Define bounds for variables (concentrations should be >= 0)
bounds = [[0, 1], [0, 1], [0, 1], [0, 1], [0, 1]]
# Use minimize from SciPy with BFGS method to find the equilibrium concentrations
result = minimize(equations, initial_guess, args=(Ka, Kw, initial_concentrations),
bounds=bounds,method='Nelder-Mead', options={'disp': True, 'xatol':1E-40, 'fatol':1E-40, 'maxfev':1000000})
return result.x
# User-supplied equilibrium constant (Ka) and initial concentration of acetic acid (HA_initial)
Ka = 1.8e-5 # mol/L, acetic acid dissociation constant
Kw = 1e-14 # mol/L, ionic product of water
initial_concentrations = [0.00, 0.1, 1e-7, 1e-7, 0.1]
initial_guess = initial_concentrations
# Extract the concentrations
HA, A_minus, H_plus, OH_minus, ion_Na = equilibrate(initial_concentrations, initial_guess, Ka, Kw)
# Calculate pH
pH = -np.log10(H_plus)
pKw = -np.log10(H_plus) -np.log10(OH_minus)
# Display the results
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}") # THIS IS THE CONDIRMATORY CHECK - Should be 14.
确认检查值到处都是,而不是14。
1.为什么Nelder-Mead算法做得这么差?是因为绝对值接近误差吗?我想我可以通过收紧公差来解决这个问题。
1.如何在fsolve中施加正约束?
1.有没有更好的方法来做到这一点?
1条答案
按热度按时间lmyy7pcs1#
不要强迫Nelder-Mead,除非你有一个很好的理由,我在这里没有看到一个。
ion_Na
实际上不是一个自由度,所以要么通过边界将其固定到其初始值,要么根本不将其包含在变量中。包括(和测试)成本函数的雅可比矩阵。
最小化显式变量
最小化,隐式变量
还有一种不同的方法,根据pKw和ion_Na的约束条件隐式设置浓度变量;它更快,但可能不正确:
根本的问题是,你的系统在不同的数量级上运行。归一化缩放和/或对数域中的计算可能会使其表现得更好。
最小化,仅限约束
这种方法当然是最准确的,但有点慢,而且-如果输入是畸形的,会产生负浓度,而不是给你一个近似的解决方案-会简单地拒绝问题是不可行的。