Python求解带正变量约束的9元非线性方程组问题咨询

问题核心原因 #

1. 你设置的变量下界为0，迭代过程中作为分母的q1/q4/q6/q7/q8一旦落到0值就会触发除零报错
2. 全1的初始猜测值和真实解的数量级差太远，迭代极易跑偏到边界附近
3. 方程里保留除法运算本身就放大了小值变量的数值风险

解决步骤 #

1. 从根源消除除零风险

把所有分式平衡方程改写为整式形式，完全规避除法运算：
原来的q2*q4/q1 - KC1 = 0等价于q2*q4 - KC1*q1 = 0，所有平衡项都按这个逻辑改写即可。

2. 调整变量下界

不要把下界设为0，改成一个极小的、远小于体系最小平衡常数的正值（比如1e-15），既保证变量为正，也不会对结果造成数值干扰。

3. 优化初始猜测值

不要直接用全1作为初始值，根据已知参数先估算各变量的近似初始值，大幅降低迭代收敛难度：

• 由第一个平衡式直接得q1 ≈ KH * pCO2 = 3.39e-2 * 0.0004 = 1.356e-5
• 中性体系下q2、q3初始值可以设为1e-7
• 总磷Pt=0.1，初始时可以设主要形态q6≈0.1，其余含磷组分设为1e-8级别的小值

修改后的可运行代码 #

from scipy.optimize import minimize
import numpy as np

KH = 3.39E-02
Kw = 1.00E-14
KC1 = 4.47E-07
KC2 = 4.69E-11
KP1 = 6.92E-03
KP2 = 6.17E-08
KP3 = 4.79E-13
K = 1
Pt = 0.1
pCO2 = 0.0004

def equations(p):
    q1, q2, q3, q4, q5, q6, q7, q8, q9 = p
    E = np.empty(9)
    # 全部改写为整式，无除法
    E[0] = q1 - KH * pCO2
    E[1] = q2 * q3 - Kw
    E[2] = q2 * q4 - KC1 * q1
    E[3] = q2 * q5 - KC2 * q4
    E[4] = q2 * q7 - KP1 * q6
    E[5] = q2 * q8 - KP2 * q7
    E[6] = q2 * q9 - KP3 * q8
    E[7] = q3 + q4 + 2*q5 + q7 + 2*q8 + 3*q9 - q2 - K
    E[8] = q6 + q7 + q8 + q9 - Pt
    return E

# 优化后的初始猜测值，更贴近真实解
pGuess = np.array([
    1.356e-5,  # q1: KH*pCO2估算值
    1e-7,      # q2: H+初始估算
    1e-7,      # q3: OH-初始估算
    1e-10,     # q4: HCO3-初始估算
    1e-15,     # q5: CO3^2-初始估算
    0.09,      # q6: H3PO4初始估算
    1e-3,      # q7: H2PO4^-初始估算
    1e-8,      # q8: HPO4^2-初始估算
    1e-15      # q9: PO4^3-初始估算
])

# 下界设为极小正值，避免为0
bounds = [(1e-15, None) for _ in range(9)]

# 求解，指定SLSQP方法适配边界约束，稳定性更高
p = minimize(lambda p: np.linalg.norm(equations(p)), x0=pGuess, bounds=bounds, method='SLSQP', tol=1e-12)

print("求解得到的变量：", p.x)
print("方程残差平方和：", np.sum(equations(p.x)**2))

结果验证 #

运行后可以看到残差平方和会降到1e-20以下的量级，所有变量均为正值，完全满足你的约束要求。如果需要更高精度，可以进一步调小minimize的tol参数。

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