QUAD
Overview
The QUAD function computes the definite integral of a function defined by a table of (x, y) data points over a specified interval using adaptive quadrature. This is essential for calculating areas, cumulative quantities, or any scenario where the integral of discrete or sampled data is needed.
This implementation leverages the scipy.integrate.quad function from the SciPy library, which is built on the QUADPACK Fortran library—a well-established collection of routines for numerical integration developed in the 1980s. QUADPACK remains widely used due to its reliability and efficiency for general-purpose integration tasks.
The function accepts tabular data (x, y pairs) and constructs a piecewise-linear interpolation to define a continuous function from the discrete points. It then applies adaptive quadrature to integrate this interpolated function between the specified limits a and b. Adaptive quadrature works by recursively subdividing the integration interval and concentrating computational effort where the integrand varies most rapidly, achieving high accuracy while minimizing function evaluations.
The underlying algorithm uses a 21-point Gauss-Kronrod quadrature rule within each subinterval, combined with global adaptive subdivision and extrapolation to handle potential singularities. The integration aims to satisfy the error criterion:
|I - \text{result}| \leq \max(\text{epsabs}, \text{epsrel} \times |I|)
where I is the true integral value and \text{result} is the numerical approximation. The default tolerances (epsabs = 1.49e-08 and epsrel = 1.49e-08) provide near machine-precision accuracy for well-behaved functions.
For more details on the mathematical foundations, see the SciPy integration documentation and the QUADPACK reference.
This example function is provided as-is without any representation of accuracy.
Excel Usage
=QUAD(function_table, a, b, epsabs, epsrel)function_table(list[list], required): Table of [x, y] pairs defining the function to integrate, with strictly increasing x valuesa(float, required): Lower limit of integrationb(float, required): Upper limit of integrationepsabs(float, optional, default: 1.49e-8): Absolute error tolerance for integrationepsrel(float, optional, default: 1.49e-8): Relative error tolerance for integration
Returns (float): The estimated value of the definite integral. str: An error message if input validation fails or integration errors occur.
Example 1: Parabola integral over [0, 4]
Inputs:
| function_table | a | b | |
|---|---|---|---|
| 0 | 0 | 0 | 4 |
| 1 | 1 | ||
| 2 | 4 | ||
| 3 | 9 | ||
| 4 | 16 |
Excel formula:
=QUAD({0,0;1,1;2,4;3,9;4,16}, 0, 4)Expected output:
22
Example 2: Sine integral over half period
Inputs:
| function_table | a | b | epsabs | |
|---|---|---|---|---|
| 0 | 0 | 0 | 3.1416 | 1e-9 |
| 1.5708 | 1 | |||
| 3.1416 | 0 |
Excel formula:
=QUAD({0,0;1.5708,1;3.1416,0}, 0, 3.1416, 1e-9)Expected output:
1.5708
Example 3: Exponential decay with custom tolerances
Inputs:
| function_table | a | b | epsabs | epsrel | |
|---|---|---|---|---|---|
| 0 | 1 | 0 | 5 | 1e-9 | 1e-9 |
| 2.5 | 0.0821 | ||||
| 5 | 0.0067 |
Excel formula:
=QUAD({0,1;2.5,0.0821;5,0.0067}, 0, 5, 1e-9, 1e-9)Expected output:
1.46363
Example 4: Linear function integral over [1, 3]
Inputs:
| function_table | a | b | epsrel | |
|---|---|---|---|---|
| 1 | 3 | 1 | 3 | 1e-9 |
| 2 | 5 | |||
| 3 | 7 |
Excel formula:
=QUAD({1,3;2,5;3,7}, 1, 3, 1e-9)Expected output:
10
Python Code
Show Code
from bisect import bisect_left
import math
from scipy.integrate import quad as scipy_integrate_quad
def quad(function_table, a, b, epsabs=1.49e-08, epsrel=1.49e-08):
"""
Numerically integrate a function defined by a table of x, y values over [a, b] using adaptive quadrature.
See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html
This example function is provided as-is without any representation of accuracy.
Args:
function_table (list[list]): Table of [x, y] pairs defining the function to integrate, with strictly increasing x values
a (float): Lower limit of integration
b (float): Upper limit of integration
epsabs (float, optional): Absolute error tolerance for integration Default is 1.49e-08.
epsrel (float, optional): Relative error tolerance for integration Default is 1.49e-08.
Returns:
float: The estimated value of the definite integral. str: An error message if input validation fails or integration errors occur.
"""
def validate_finite_float(value, name):
try:
numeric = float(value)
except Exception:
return f"Error: {name} must be a number."
if math.isnan(numeric) or math.isinf(numeric):
return f"Error: {name} must be finite."
return numeric
def validate_tolerance(value, name, allow_zero=False):
converted = validate_finite_float(value, name)
if isinstance(converted, str):
return converted
if allow_zero:
if converted < 0.0:
return f"Error: {name} must be non-negative."
else:
if converted <= 0.0:
return f"Error: {name} must be positive."
return converted
try:
# Validate function_table structure
if not isinstance(function_table, list) or len(function_table) < 2:
return "Error: function_table must be a 2D list with at least two rows."
processed_rows = []
for i, row in enumerate(function_table):
if not isinstance(row, list) or len(row) < 2:
return "Error: each row in function_table must contain at least two values."
try:
x_val = float(row[0])
y_val = float(row[1])
except Exception:
return "Error: function_table must contain numeric values."
# Validate that values are finite
if math.isnan(x_val) or math.isinf(x_val):
return f"Error: x value at row {i} must be finite."
if math.isnan(y_val) or math.isinf(y_val):
return f"Error: y value at row {i} must be finite."
processed_rows.append([x_val, y_val])
xs = [row[0] for row in processed_rows]
ys = [row[1] for row in processed_rows]
if any(x2 <= x1 for x1, x2 in zip(xs, xs[1:])):
return "Error: x values in function_table must be strictly increasing."
converted_a = validate_finite_float(a, "a")
if isinstance(converted_a, str):
return converted_a
converted_b = validate_finite_float(b, "b")
if isinstance(converted_b, str):
return converted_b
if converted_b < converted_a:
return "Error: b must be greater than or equal to a."
x_min, x_max = xs[0], xs[-1]
if converted_a < x_min or converted_b > x_max:
return f"Error: integration limits [{converted_a}, {converted_b}] must fall within function_table x range [{x_min}, {x_max}]."
converted_epsabs = validate_tolerance(epsabs, "epsabs", allow_zero=True)
if isinstance(converted_epsabs, str):
return converted_epsabs
converted_epsrel = validate_tolerance(epsrel, "epsrel", allow_zero=False)
if isinstance(converted_epsrel, str):
return converted_epsrel
def interpolated_function(x_value):
# Find the right interval using binary search
index = bisect_left(xs, x_value)
# Handle edge cases
if index == 0:
return ys[0]
if index >= len(xs):
return ys[-1]
# Linear interpolation between adjacent points
x0, x1 = xs[index - 1], xs[index]
y0, y1 = ys[index - 1], ys[index]
weight = (x_value - x0) / (x1 - x0)
return y0 + weight * (y1 - y0)
try:
result, _ = scipy_integrate_quad(
interpolated_function,
converted_a,
converted_b,
epsabs=converted_epsabs,
epsrel=converted_epsrel,
)
except Exception as exc: # noqa: BLE001 - repackage SciPy errors as strings
return f"Error: integration failed: {exc}"
# Validate result is finite
if not isinstance(result, (int, float)):
return "Error: integration returned a non-numeric result."
numeric_result = float(result)
if math.isnan(numeric_result) or math.isinf(numeric_result):
return "Error: integration returned a non-finite result."
return numeric_result
except Exception as e:
return f"Error: {str(e)}"