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 \, dx
This 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.
**Example 1: Rectangular region with x*y^2**
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.666667
Example 2: Curved bounds between sin(x) and cos(x)
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.414214
Example 3: Gaussian integral over [-1,1] x [-1,1]
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.23099
Example 4: Linear function with variable y bounds
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
Show 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"Error: {name} must be a number."
if math.isnan(converted) or math.isinf(converted):
return f"Error: {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 "Error: expression arguments must be strings."
processed = re.sub(r"\^", "**", expr.strip())
if not processed:
return "Error: 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 "Error: 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
try:
# 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 "Error: 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"Error: {error}"
except Exception as exc: # noqa: BLE001 - want to surface SciPy errors
return f"Error: scipy.integrate.dblquad failed: {exc}"
if not isinstance(result, (int, float)):
return "Error: scipy.integrate.dblquad returned a non-numeric result."
numeric_result = float(result)
if math.isnan(numeric_result) or math.isinf(numeric_result):
return "Error: scipy.integrate.dblquad returned a non-finite result."
return numeric_result
except Exception as e:
return f"Error: {str(e)}"