如何求解特定多变量非线性方程组的全部实根

给定以下用于计算平衡浓度的方程组：

X * (1 + K_AC * C + K_AB * Y) = A_0
Y * (1 + K_BC * C + K_AB * X) = Z
C * (1 + K_AC * X + K_BC * Y) = C_0

已知条件 #

• 参数 K_AB、K_AC、K_BC 为已知正参数（存在复杂相互依赖关系）
• A_0、C_0、C 为已知正数，且满足 C < C_0

求解需求 #

• 得到方程组的全部实根，关注 X>0、Y 和 Z 为实数的解（X、Y、Z 均为正的解是否存在由 A_0、C_0、C 决定）
• 推测解的数量少于100，需支持频繁求解不同参数下的方程组，要求高效且支持设置搜索边界（类似 scipy.optimize.Bounds），优先选择牛顿法、hybr法这类方法，不考虑basinhoppin或二分法

已尝试方法 #

使用过 scipy.optimize 模块的 root、fsolve 等函数，但仅能找到单个解，无法获取全部解

可行方案建议 #

1. 多初始点遍历+局部优化
   • 基于参数范围生成一批覆盖合理搜索边界的初始点集合，对每个初始点调用 scipy.optimize.root（选用hybr或lm方法，需边界约束时用trust-constr）求解，最后通过设定容差（如1e-8）去重，保留独立解
   • 初始点生成可结合第三个方程的线性约束：K_AC*X + K_BC*Y = (C_0/C) - 1，在 X>0 范围内生成满足该约束的(X,Y)组合，再代入第一个方程反推合理初始值，缩小搜索范围
2. 符号推导降维求解
   • 通过代数变换简化方程组：从第三个方程解出 K_AC*X = (C_0/C -1) - K_BC*Y，代入第一个方程消去X，得到关于Y的单变量方程；用sympy将其转化为多项式形式后，调用numpy.roots获取所有实根，再回代计算X和Z
   • 该方法能直接获取全部实根，效率极高，适合频繁求解不同参数的场景，还可直接结合X>0的条件筛选有效解
3. 全局优化器结合边界约束
   • 使用支持边界约束的全局优化工具，比如scipy.optimize.shgo（基于单纯形的全局优化方法），它能在指定边界内搜索所有局部解，无需手动生成大量初始点，适合参数频繁变动的场景

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