使用scipy fsolve求解耦合方程组多根遇异常，求排查思路

现有如下耦合方程组：

two_exponential = lambda x, kernel, c: np.array([x[0] - np.exp(kernel[0] * x[0] + kernel[2] * x[1] + c), x[1] - np.exp(kernel[1] * x[1] + kernel[3] * x[0] + c)])

希望使用scipy.fsolve求解该方程组的根（交点），当前实现是针对b11,b22, b12, b21的不同配置，遍历指定范围内的初始值寻找根：

b = np.array([b11, b22, b12, b21])
x_min_plot = -10
x_max_plot = 35
x_1 = np.linspace(x_min_plot, x_max_plot, 100)
x_2 = np.linspace(x_min_plot, x_max_plot, 100)
x_1, x_2 = np.meshgrid(x_1, x_2)
z_1 = -x_1 + np.exp(b[0] * x_1 + b[2] * x_2 + c)
z_2 = -x_2 + np.exp(b[1] * x_2 + b[3] * x_1 + c)
x_sols = []
x_min = 0
x_max = 35

for x in np.arange(x_min, x_max, 5):
    for y in np.arange(x_min, x_max, 5):
        initial = np.array([x, y])
        x_sol = fsolve(two_exponential, initial, args=(b, c), full_output=1)
        if x_sol[2] == 1: # if the solution converged
            x_sols.append(np.round(x_sol[0], 2))
# [x for i, x in enumerate(x_sols) if not np.isclose(x, x_sols[i-1], atol = 1e-1).all()]
x_sols = np.unique(x_sols, axis=0)


print(f'z*: {np.round(x_sols, 2)}')
if x_sol[2] != 1:
    print('no solution')

通过对结果取整来去重，仅保留唯一根，但该代码在部分参数配置下无法正确求解，请问可能的问题原因是什么？

• 初始值覆盖范围不全：当前初始值只遍历了[0,35)区间、步长为5的点，但绘图范围是[-10,35]，如果根落在负数区间或者遍历点的间隙中，fsolve根本没机会收敛到这些根；且步长5过大，容易跳过靠近根的初始点，导致迭代发散。
• 指数函数的强非线性特性：方程组包含指数项，当参数b或c使得指数项快速增长/衰减时，函数导数会剧烈变化，fsolve依赖的牛顿法容易出现Jacobian矩阵奇异、迭代步长失控的情况，直接导致收敛失败。
• 去重逻辑不合理：用np.round(x_sol[0],2)配合np.unique去重精度太粗糙——差异小于0.01的不同根会被误判为同一个，同一根的微小迭代误差又可能被当成不同根保留；注释里用np.isclose的判断逻辑更科学，却被弃用了。
• 收敛结果判断错误：最后判断“无 solution”时，用的是最后一次循环的x_sol状态，哪怕前面已经找到多个收敛的根，只要最后一次迭代没收敛，就会输出“no solution”，完全误导结果判断。
• 未调整fsolve的收敛参数：fsolve的默认收敛阈值（如xtol、ftol）可能不匹配你的方程组，比如根附近函数值变化极小时，默认阈值可能无法触发收敛判定，明明接近根却被判定为失败，可尝试手动设置更高精度参数（如xtol=1e-8）。

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