How to compute a definite integral with SciPy

Definite integrals turn a function and interval into an accumulated area or total quantity. In Python analysis code, SciPy handles cases where an analytic antiderivative is inconvenient, unavailable, or not worth deriving for a one-off numerical check.

The scipy.integrate.quad() function evaluates a one-variable callable between two finite or infinite bounds and returns both the estimated integral and an absolute error estimate. Extra function parameters go through args, which keeps the integration variable as the first argument and leaves model constants explicit.

A parameterized quadratic keeps the output easy to audit because the exact value is 5 / 3 over 0 to 1. The same pattern applies to probability densities, fitted model functions, and other scalar callables as long as the function returns finite values over the interval.

1. Create a Python script named integrate_definite.py.

   integrate_definite.py

   from math import isclose
    
   from scipy.integrate import quad
    
    
   def integrand(x, scale, offset):
       return scale * x**2 + offset
    
    
   value, error = quad(integrand, 0.0, 1.0, args=(2.0, 1.0))
   expected = 5.0 / 3.0
    
   print(f"integral: {value:.12f}")
   print(f"estimated_error: {error:.2e}")
   print(f"expected: {expected:.12f}")
   print(f"matches_expected: {isclose(value, expected, rel_tol=1e-12, abs_tol=error * 10)}")

   quad() calls the integrand with the integration variable first. Values in args fill the remaining function parameters.

2. Run the script and confirm the integral value.

   $ python3 integrate_definite.py
   integral: 1.666666666667
   estimated_error: 1.85e-14
   expected: 1.666666666667
   matches_expected: True

3. Set the integrand callable with the integration variable first.

   def integrand(x, scale, offset):
       return scale * x**2 + offset

   The callable must return a scalar value for each scalar x supplied by quad().

4. Call quad() with the interval and extra arguments.

   value, error = quad(integrand, 0.0, 1.0, args=(2.0, 1.0))

   The returned error is SciPy's estimate of the absolute integration error, not a comparison against a known exact result. Use numpy.inf for an unbounded integration limit.

5. Compare the result with the known value for the quadratic integral.

   expected = 5.0 / 3.0
   print(f"matches_expected: {isclose(value, expected, rel_tol=1e-12, abs_tol=error * 10)}")

6. Tighten the tolerance only when the integrand is smooth enough for that accuracy.

   quad() uses default absolute and relative tolerances near 1.49e-8. Pass epsabs or epsrel when the default error target is too loose or unnecessarily strict for the calculation.

7. Remove integrate_definite.py after adapting the calculation.

   $ rm integrate_definite.py