DBLQUAD
Overview
The DBLQUAD function computes the double integral of a function f(y, x) over a two-dimensional region, where the outer integral runs from x = a to x = b, and the inner integral runs from y = g(x) to y = h(x). Mathematically, it evaluates:
\int_{a}^{b} \int_{g(x)}^{h(x)} f(y, x) \, dy \, dxThis implementation wraps the scipy.integrate.dblquad function from the SciPy scientific computing library. Under the hood, dblquad performs two nested one-dimensional integrations using routines from QUADPACK, a well-established FORTRAN library for numerical quadrature originally developed in 1983.
The integration algorithm uses adaptive quadrature with 21-point Gauss-Kronrod rules within each subinterval, combined with Wynn’s epsilon algorithm for series acceleration. This approach can handle a variety of integrand behaviors, including those with singularities, by adaptively refining the integration mesh where needed. For infinite bounds, the algorithm maps the infinite range onto a finite interval before applying the same adaptive strategy.
The function accepts both constant boundary values and variable boundaries specified as expressions of x. The epsabs and epsrel parameters control the absolute and relative error tolerances respectively, with default values of 1.49 \times 10^{-8}. The algorithm attempts to achieve an error satisfying:
|I - \text{result}| \leq \max(\text{epsabs}, \text{epsrel} \cdot |I|)where I is the true value of the inner integral.
Double integrals are fundamental in applications such as computing areas, volumes, center of mass, probability distributions over two-dimensional regions, and evaluating physical quantities in fields like electromagnetics and fluid dynamics. For related integration functions, see SciPy’s quad for single integrals, tplquad for triple integrals, and nquad for N-dimensional integrals.
This example function is provided as-is without any representation of accuracy.
Excel Usage
=DBLQUAD(func_expr, a, b, gfun_expr, hfun_expr, epsabs, epsrel)func_expr(str, required): Expression for the function f(y, x) to integratea(float, required): Lower limit of integration for xb(float, required): Upper limit of integration for xgfun_expr(str, required): Expression or constant for the lower boundary y = g(x)hfun_expr(str, required): Expression or constant for the upper boundary y = h(x)epsabs(float, optional, default: 1.49e-8): Absolute tolerance for the integrationepsrel(float, optional, default: 1.49e-8): Relative tolerance for the integration
Returns (float): The computed double integral, or an error message (str) if invalid.
Examples
Example 1: Demo case 1
Inputs:
| func_expr | a | b | gfun_expr | hfun_expr |
|---|---|---|---|---|
| x * y^2 | 0 | 2 | 0 | 1 |
Excel formula:
=DBLQUAD("x * y^2", 0, 2, "0", "1")Expected output:
0.6667
Example 2: Demo case 2
Inputs:
| func_expr | a | b | gfun_expr | hfun_expr | epsabs |
|---|---|---|---|---|---|
| 1 | 0 | 0.7853981633974483 | sin(x) | cos(x) | 1e-9 |
Excel formula:
=DBLQUAD("1", 0, 0.7853981633974483, "sin(x)", "cos(x)", 1e-9)Expected output:
0.4142
Example 3: Demo case 3
Inputs:
| func_expr | a | b | gfun_expr | hfun_expr | epsabs | epsrel |
|---|---|---|---|---|---|---|
| exp(-(x2 + y2)) | -1 | 1 | -1 | 1 | 1e-9 | 1e-9 |
Excel formula:
=DBLQUAD("exp(-(x**2 + y**2))", -1, 1, "-1", "1", 1e-9, 1e-9)Expected output:
2.231
Example 4: Demo case 4
Inputs:
| func_expr | a | b | gfun_expr | hfun_expr | epsrel |
|---|---|---|---|---|---|
| 3 * x * y | 0 | 1 | x | 2 - x | 1e-9 |
Excel formula:
=DBLQUAD("3 * x * y", 0, 1, "x", "2 - x", 1e-9)Expected output:
1
Python Code
import math
import re
from scipy.integrate import dblquad as scipy_dblquad
def dblquad(func_expr, a, b, gfun_expr, hfun_expr, epsabs=1.49e-08, epsrel=1.49e-08):
"""
Compute the double integral of a function over a two-dimensional region.
See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.dblquad.html
This example function is provided as-is without any representation of accuracy.
Args:
func_expr (str): Expression for the function f(y, x) to integrate
a (float): Lower limit of integration for x
b (float): Upper limit of integration for x
gfun_expr (str): Expression or constant for the lower boundary y = g(x)
hfun_expr (str): Expression or constant for the upper boundary y = h(x)
epsabs (float, optional): Absolute tolerance for the integration Default is 1.49e-08.
epsrel (float, optional): Relative tolerance for the integration Default is 1.49e-08.
Returns:
float: The computed double integral, or an error message (str) if invalid.
"""
def _convert_float(value, name):
try:
converted = float(value)
except Exception:
return f"Invalid input: {name} must be a number."
if math.isnan(converted) or math.isinf(converted):
return f"Invalid input: {name} must be finite."
return converted
def _parse_callable(expr, varnames):
# expr must be a string expression that can be evaluated to a numeric value
if not isinstance(expr, str):
return "Invalid input: expression arguments must be strings."
processed = re.sub(r"\^", "**", expr.strip())
if not processed:
return "Invalid input: expression arguments must not be empty."
allowed = {"__builtins__": {}}
for _name in dir(math):
if _name.startswith("_"):
continue
func = getattr(math, _name)
allowed[_name] = func
allowed[_name.upper()] = func
for varname in varnames:
allowed[varname] = 0.0
allowed[varname.upper()] = 0.0
try:
code = compile(processed, "<dblquad_expr>", "eval")
except Exception:
return "Invalid input: unable to parse expression."
def _callable(*args):
local_env = {}
for var, value in zip(varnames, args):
numeric = float(value)
local_env[var] = numeric
local_env[var.upper()] = numeric
try:
result = eval(code, allowed, local_env)
except Exception as exc: # noqa: BLE001 - need to surface parse errors cleanly
raise ValueError(f"Expression evaluation failed: {exc}") from exc
if isinstance(result, bool):
return float(result)
if isinstance(result, (int, float)):
numeric = float(result)
if math.isnan(numeric) or math.isinf(numeric):
raise ValueError("Expression produced non-finite value.")
return numeric
raise ValueError("Expression must evaluate to a numeric value.")
return _callable
# Validate scalar numeric inputs
converted_a = _convert_float(a, "a")
if isinstance(converted_a, str):
return converted_a
converted_b = _convert_float(b, "b")
if isinstance(converted_b, str):
return converted_b
if converted_b < converted_a:
return "Invalid input: lower limit must be less than or equal to upper limit."
converted_epsabs = _convert_float(epsabs, "epsabs")
if isinstance(converted_epsabs, str):
return converted_epsabs
converted_epsrel = _convert_float(epsrel, "epsrel")
if isinstance(converted_epsrel, str):
return converted_epsrel
parsed_func = _parse_callable(func_expr, ["y", "x"])
if isinstance(parsed_func, str):
return parsed_func
if isinstance(gfun_expr, str):
parsed_gfun = _parse_callable(gfun_expr, ["x"])
if isinstance(parsed_gfun, str):
return parsed_gfun
else:
g_value = _convert_float(gfun_expr, "gfun_expr")
if isinstance(g_value, str):
return g_value
def parsed_gfun(x: float) -> float:
return g_value
if isinstance(hfun_expr, str):
parsed_hfun = _parse_callable(hfun_expr, ["x"])
if isinstance(parsed_hfun, str):
return parsed_hfun
else:
h_value = _convert_float(hfun_expr, "hfun_expr")
if isinstance(h_value, str):
return h_value
def parsed_hfun(x: float) -> float:
return h_value
try:
result, _ = scipy_dblquad(
parsed_func,
converted_a,
converted_b,
parsed_gfun,
parsed_hfun,
epsabs=converted_epsabs,
epsrel=converted_epsrel,
)
except ValueError as error:
return f"Invalid input: {error}"
except Exception as exc: # noqa: BLE001 - want to surface SciPy errors
return f"scipy.integrate.dblquad error: {exc}"
if not isinstance(result, (int, float)):
return "scipy.integrate.dblquad error: non-numeric result."
numeric_result = float(result)
if math.isnan(numeric_result) or math.isinf(numeric_result):
return "scipy.integrate.dblquad error: result is not finite."
return numeric_result