Double Integration with SciPy's dblquad()

Double integration is a cornerstone of calculus and numerical analysis, essential for calculating areas, volumes, and solving various scientific and engineering problems. In Python, the scipy.integrate module offers the powerful dblquad() function to perform double integrals over specified regions in the XY-plane. This guide provides a comprehensive overview of its usage.

What is Double Integration?

Double integration involves computing the integral of a function with two variables, denoted as $f(x, y)$, over a defined two-dimensional region $R$. Mathematically, it is represented as:

$$ \iint_R f(x, y) ,dA $$

When integrating over a rectangular region defined by $a \le x \le b$ and $c \le y \le d$, or a more general region, it is often expressed as iterated integrals. For a region where $y$ varies from $g(x)$ to $h(x)$ and $x$ varies from $a$ to $b$, the integral becomes:

$$ \int_a^b \left( \int_{g(x)}^{h(x)} f(x, y) ,dy \right) ,dx $$

The process involves evaluating the inner integral (with respect to $y$) first, treating $x$ as a constant. The result of this inner integral, which will be a function of $x$, is then integrated with respect to $x$ over its specified bounds.

Overview of scipy.integrate.dblquad()

The scipy.integrate.dblquad() function in SciPy is designed to numerically compute double integrals. It implements an adaptive quadrature algorithm that automatically adjusts to the integrand's behavior, aiming to meet specified error tolerances. The function evaluates the integral by performing nested integrations.

Syntax

scipy.integrate.dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)

Parameters

• func: The function to be integrated. It must accept two arguments, conventionally y and x (in that order), representing the variables of integration. The order is important for how dblquad passes arguments.
• a, b: The lower and upper bounds for the outer integral (typically with respect to $x$). These can be constants or functions of other independent variables if needed, but for standard double integration, they are constants.
• gfun, hfun: Functions that define the inner integral's bounds (typically with respect to $y$). gfun represents the lower bound for $y$, and hfun represents the upper bound for $y$. These functions can depend on the outer integration variable ($x$).
• args (optional): A tuple of extra arguments to pass to func. This is useful if your integrand requires additional parameters.
• epsabs (optional): The absolute error tolerance for the integration. Defaults to 1.49e-08.
• epsrel (optional): The relative error tolerance for the integration. Defaults to 1.49e-08.

Return Values

The function returns a tuple containing:

1. The estimated value of the integral.
2. An estimate of the absolute error in the result.

Examples

Example 1: Basic Double Integration

Integrate the function $f(x,y) = x \cdot y$ over the region defined by $0 \le x \le 1$ and $0 \le y \le 2x$.

import scipy.integrate as integrate
import numpy as np

# Define the integrand. Note the order of arguments: y, x
def integrand(y, x):
    return x * y

# Define the bounds for y as functions of x
lower_y_bound = lambda x: 0
upper_y_bound = lambda x: 2 * x

# Perform the double integration
# Outer integral bounds (x): 0 to 1
# Inner integral bounds (y): 0 to 2*x
result, error = integrate.dblquad(
    integrand,
    0,          # a: lower bound for x
    1,          # b: upper bound for x
    lower_y_bound, # gfun: lower bound for y (as a function of x)
    upper_y_bound  # hfun: upper bound for y (as a function of x)
)

print(f"Result: {result}")
print(f"Estimated Error: {error}")

Output:

Result: 0.5
Estimated Error: 2.2060128823111155e-14

Example 2: Handling Infinite Limits

SciPy's dblquad() can handle integration over infinite bounds by using np.inf or -np.inf. Let's integrate $f(x,y) = e^{-(x^2 + y^2)}$ over the first quadrant ($x \ge 0, y \ge 0$).

import scipy.integrate as integrate
import numpy as np

# Define the integrand
def integrand_gaussian(y, x):
    return np.exp(-x**2 - y**2)

# Perform the double integration over infinite bounds
# x from 0 to infinity
# y from 0 to infinity
result_inf, error_inf = integrate.dblquad(
    integrand_gaussian,
    0,          # a: lower bound for x
    np.inf,     # b: upper bound for x
    lambda x: 0, # gfun: lower bound for y
    lambda x: np.inf # hfun: upper bound for y
)

print(f"Result (infinite bounds): {result_inf}")
print(f"Estimated Error (infinite bounds): {error_inf}")

Output:

Result (infinite bounds): 0.7853981633973343
Estimated Error (infinite bounds): 1.4647640380321503e-08

(Note: The theoretical result for this integral is $\pi/4 \approx 0.785398$)

Example 3: Custom Error Tolerance

You can adjust the accuracy of the integration by modifying epsabs and epsrel. Let's use tighter tolerances for the previous Gaussian integral.

import scipy.integrate as integrate
import numpy as np

def integrand_gaussian(y, x):
    return np.exp(-x**2 - y**2)

# Perform integration with custom tolerances
result_custom_tol, error_custom_tol = integrate.dblquad(
    integrand_gaussian,
    0, np.inf,
    lambda x: 0,
    lambda x: np.inf,
    epsabs=1e-10, # Tighter absolute tolerance
    epsrel=1e-10  # Tighter relative tolerance
)

print(f"Result (custom tolerance): {result_custom_tol}")
print(f"Estimated Error (custom tolerance): {error_custom_tol}")

Output:

Result (custom tolerance): 0.785398163397448
Estimated Error (custom tolerance): 9.692990732633506e-11

Example 4: Integrating Complex Functions

SciPy's dblquad() does not directly support complex-valued integrands. To integrate a complex function $f(x,y) = u(x,y) + i \cdot v(x,y)$, you need to integrate its real and imaginary parts separately and then combine the results.

Let's integrate $f(x,y) = e^{-(x^2+y^2)} + i \cdot \sin(x+y)$ over a large finite region to approximate an infinite one (e.g., $0 \le x \le 100, 0 \le y \le 100$).

import scipy.integrate as integrate
import numpy as np

# Define the real part of the complex function
def real_part_integrand(y, x):
    return np.exp(-x**2 - y**2)

# Define the imaginary part of the complex function
def imag_part_integrand(y, x):
    return np.sin(x + y)

# Define a large finite bound to approximate infinite region
bound = 100

# Integrate the real part
real_result, real_error = integrate.dblquad(
    real_part_integrand,
    0, bound,
    lambda x: 0,
    lambda x: bound
)

# Integrate the imaginary part
imag_result, imag_error = integrate.dblquad(
    imag_part_integrand,
    0, bound,
    lambda x: 0,
    lambda x: bound
)

# Combine the results to get the complex integral
complex_result = real_result + 1j * imag_result

print(f"Complex Result: {complex_result}")
print(f"Real Part Error: {real_error}")
print(f"Imaginary Part Error: {imag_error}")

Output:

Complex Result: (0.7853981633971309-0.13943398500584472j)
Real Part Error: 1.1078396381028211e-08
Imaginary Part Error: 1.4892850109719667e-08

Conclusion

scipy.integrate.dblquad() is a versatile and robust tool for numerical double integration in Python. It efficiently handles various scenarios, from simple finite regions to integrals over infinite domains, and allows fine-grained control over accuracy through error tolerances. By understanding its parameters and proper usage, you can effectively leverage it for your scientific and mathematical computations.