IP

OpenOpt has solver interalg based on interval analysis, that is capable of getting solution with specifiable accuracy.

Function evaluations based solvers (like SciPy quad, dblquad, tplquad) can miss some essential regions, like in the case from the picture (and situations in multidimensional spaces can be even more complex). Unlike them, interalg is based on interval analysis and is capable of getting solution with required precision (that is mathematically proved). It operates on pieces of space instead of function evaluation in points and never lies about obtained residual (except of small rounding errors due to finite lattice numbers representation effects, e.g. for float64 type it's about 10^-16, i.e. with the type 1 + 10^-16 = 1).

Also, SciPy routines work in space with dimension up to 3, while interalg can handle any dimensional space.

interalg drawbacks:

• it may work slower than his function evaluation based competitors
• only functions with overloaded interval arithmetic operations are supported (for now it's +, -, *, /, pow (**), sin, cos, arcsin, arccos, arctan, sinh, cosh, exp, sqrt, abs, log, log2, log10, floor, ceil, sum, min, max, sorted monotone 1-d splines (R->R) with order 1 or 3; since v. 0.38: arcsinh, arccosh, tanh, arctanh; future plans: general 1-d splines, N-dimensional splines, erf, tan, cotan, maybe others (some monotone, unimodal functions or functions with all known local extrema could be connected on demand))
• currently only box domains D are handled; FuturePlans: add constraints like lb_i <= c_i(x) <= ub_i to integration domain (D) definition.