This answer to this question works only for situations in which the desired solution to the coupled functions is not restricted to a certain range.

But what if, for example, we wanted a solution such that 0 < x < 10 and 0 < y < 10?

There are functions within scipy.optimize that find roots to a function within a given interval (e.g., brentq), but these work only for functions of one variable.

Why does scipy fall short of providing a root solver that works for multi-variable functions within specific ranges? How might such a solver be implemented?

  • 1
    $\begingroup$ The param to be optimised, x0 is a ndarray . . . so your "one variable" is an array - i.e. any number of scalar or concatenated array variables. In what way is this not working for multi-variable functions? $\endgroup$ Sep 14, 2015 at 6:24
  • $\begingroup$ I upvoted @NeilSlater's comment, but I shouldn't have. The question is about multivariate (rectangular) constraints, not about multivariate argument. $\endgroup$
    – Valentas
    Sep 14, 2015 at 6:55
  • $\begingroup$ @NeilSlater when you say, "in what way is this not working," what is the "this" you're referring to? i.e., which scipy function are you referring to? the different functions aren't suitable to my case for different reasons. $\endgroup$
    – abcd
    Sep 14, 2015 at 7:24
  • $\begingroup$ Sorry, I missed the constraints as being the main blocker to using fsolve, and thought you had in turn missed the first param being an array. Perhaps making the question more self-contained would have helped me spot the difference. $\endgroup$ Sep 14, 2015 at 7:55

1 Answer 1


As a workaround, you could minimize another function that includes both the objective and the constraints, then check if sol.fun is (numerically) equal to zero.

from scipy.optimize import minimize
import numpy as np

f = lambda x: np.sin(x).sum() #the function to find roots of
L = np.array([-1,-1]) #lower bound for each coordinate
U = np.array([1, 1])  #upper bound for each coordinate
g = lambda x: f(x) **2 +  max(0, (L - x).max(), (x - U).max())
sol = minimize(g, [0.5,0.5])

Also, scipy.optimize seems to have some optimisers that support rectangular bounds, i.e. differential_evolution (since version 0.15.0).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.