如何无需预知正确解即可用scipy fsolve求解非线性方程组

局部非线性求解器（如fsolve）基于梯度迭代思路，对初始值高度敏感，容易收敛到局部伪解，可通过以下方案稳定求解：

方案1：优先化简方程组，降低求解维度 #

你给出的方程组可直接通过变量替换降维：根据eq2和eq3可得z = y = 100000/x，代入eq1后原三维方程组直接简化为单变量方程：

x - 0.095 * np.log(100000 / x) - 1.2022 = 0

单变量方程求解难度极低，可使用不需要初始值、仅需给定取值区间的区间求根法（如scipy.optimize.brentq）稳定求解，只要保证区间两个端点的函数值符号相反即可：

from scipy.optimize import brentq
import numpy as np

def single_eq(x):
    return x - 0.095 * np.log(100000 / x) - 1.2022

# 仅需给定大致取值区间，不需要知道精确解位置
x_sol = brentq(single_eq, a=1, b=10)
y_sol = 100000 / x_sol
z_sol = y_sol
print(x_sol, y_sol, z_sol)

该方案优先级最高，降维后求解稳定性会提升数个量级。

方案2：无法化简时，先用全局优化算法定位初始值 #

如果方程组复杂度高无法做代数化简，可以先构造残差平方和目标函数，用全局优化算法找到近似解，再将其作为fsolve的初始值做局部精修：

from scipy.optimize import fsolve, differential_evolution
import numpy as np

def equations(vars):
    x, y, z = vars
    eq1 = x - 0.095*np.log(z) - 1.2022
    eq2 = y - 100000/x
    eq3 = z - y
    return [eq1, eq2, eq3]

# 构造残差平方和函数，作为全局优化的目标
def residual_sum(vars):
    return np.sum(np.square(equations(vars)))

# 给变量设置合理的取值边界，不需要知道精确解
bounds = [(0.1, 10), (1, 1e6), (1, 1e6)]
# 全局优化找近似解
res = differential_evolution(residual_sum, bounds=bounds)
# 用近似解作为初始值调用fsolve精修
x, y, z = fsolve(equations, res.x)
print(x, y, z)

全局优化算法不会被局部极小值困住，能在你给定的合理边界内找到接近真实解的位置。

方案3：使用带约束的求解器替代fsolve #

fsolve本身不支持变量边界约束，很容易跑到不合理的取值区间（比如对数的参数接近0的位置），可以改用scipy.optimize.root中支持边界约束的求解方法，比如trust-constr，只要给定合理的变量上下界，就能避免收敛到伪解。

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