Edited to take into account suggestions and observations by @FrédéricC and @nbruin.
Sage's numerical_integral
function uses a heuristic to guess an error bound,
provide no certified error bound.
The good way to go is to use Arb, which is easy from Sage.
Example (from http://fredrikj.net/math/scan2018.pdf):
sage: C = ComplexBallField(100)
sage: C.integral(lambda x, _: cos(x) * sin(x), 0, 1)
[0.35403670913678559674939205737 +/- 8.68e-30]
Side observation: the file gauss_legendre.pyx
in src/sage/numerical/
strangely does not result in a gauss_legendre
entry in the reference manual
under sage/numerical
.
Did you look at the doc: numerical_integral?