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.
Examples
Example 1: Demo case 1
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: Demo case 2
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.5707999999999998
Example 3: Demo case 3
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.4636250000000002
Example 4: Demo case 4
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
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):
"""Convert and validate that a value is a finite float."""
try:
numeric = float(value)
except Exception:
return f"Invalid input: {name} must be a number."
if math.isnan(numeric) or math.isinf(numeric):
return f"Invalid input: {name} must be finite."
return numeric
def validate_tolerance(value, name, allow_zero=False):
"""Validate that a tolerance value is non-negative (or positive if allow_zero=False)."""
converted = validate_finite_float(value, name)
if isinstance(converted, str):
return converted
if allow_zero:
if converted < 0.0:
return f"Invalid input: {name} must be non-negative."
else:
if converted <= 0.0:
return f"Invalid input: {name} must be positive."
return converted
# Validate function_table structure
if not isinstance(function_table, list) or len(function_table) < 2:
return "Invalid input: 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 "Invalid input: 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 "Invalid input: function_table must contain numeric values."
# Validate that values are finite
if math.isnan(x_val) or math.isinf(x_val):
return f"Invalid input: x value at row {i} must be finite."
if math.isnan(y_val) or math.isinf(y_val):
return f"Invalid input: 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 "Invalid input: 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 "Invalid input: 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"Invalid input: 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):
"""Perform piecewise-linear interpolation at 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"Integration error: {exc}"
# Validate result is finite
if not isinstance(result, (int, float)):
return "Integration error: non-numeric result returned."
numeric_result = float(result)
if math.isnan(numeric_result) or math.isinf(numeric_result):
return "Integration error: result is not finite."
return numeric_result