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.
There are some ideas how to speedup interalg on IP in many times or even orders, but they are not implemented yet.
Maybe some scipy or other solvers will be connected in future. Also, maybe binding IP to oofun will be performed, like
f4 = integrator(f1+2*f2+3*f3, domain, engine = 'interalg')
See also: ODE