SciPy - Integrate

SciPy - Numerical Integration (scipy.integrate)

Numerical integration is essential when a function cannot be integrated analytically or when the analytical solution is too complex. SciPy’s scipy.integrate module provides a suite of powerful and versatile functions for numerical integration (also known as quadrature).

Key functions in scipy.integrate include:

Function	Description
quad	General-purpose single integration (adaptive quadrature).
dblquad	General-purpose double integration.
tplquad	General-purpose triple integration.
nquad	General-purpose n-dimensional multiple integration.
fixed_quad	Integrate func(x) using Gaussian quadrature of order n.
quadrature	Integrate with Gaussian quadrature to a specified tolerance.
romberg	Romberg integration.
trapz or trapezoid	Use trapezoidal rule on discrete data points.
simps or simpson	Use Simpson’s rule on discrete data points.
solve_ivp	Solve initial value problems for systems of ordinary differential equations (ODEs).
solve_bvp	Solve boundary value problems for systems of ODEs.

(Note: trapz and simps are now also available as trapezoid and simpson for clarity. polyint and poly1d are NumPy functions often used in conjunction with integration.)

Single Integrals with quad

The quad function is the primary tool for computing a definite integral of a function of one variable: $\int_{a}^{b} f(x)dx$.

The basic syntax is scipy.integrate.quad(func, a, b, args=(), ...), where:

• func is the Python function or callable to integrate.
• a is the lower limit of integration.
• b is the upper limit of integration.
• args is an optional tuple of extra arguments to pass to func.

quad returns a tuple (integral_value, absolute_error_estimate).

Let’s integrate the Gaussian function $f(x) = e^{-x^2}$ from 0 to 1:

import numpy as np
import scipy.integrate

# Define the function to integrate
def gaussian_function(x):
    return np.exp(-x**2)

# Alternatively, using a lambda function:
# gaussian_function_lambda = lambda x: np.exp(-x**2)

# Set integration limits
lower_limit = 0
upper_limit = 1

# Perform the integration
integral_value, error_estimate = scipy.integrate.quad(gaussian_function, lower_limit, upper_limit)

print(f"Integral of e^(-x^2) from {lower_limit} to {upper_limit}:")
print(f"Value: {integral_value}")
print(f"Estimated Absolute Error: {error_estimate}")

Output:

Integral of e^(-x^2) from 0 to 1:
Value: 0.7468241328124271
Estimated Absolute Error: 8.291413475940725e-15

The quad function can handle infinite limits using np.inf and -np.inf. It can also integrate many standard NumPy functions directly, provided they take a single array argument and return an array.

Example with parameters using the args argument:

import numpy as np
import scipy.integrate

# Define a function f(x, a, b) = a*x + b
def linear_func(x, a, b):
    return a * x + b

# Integrate f(x, 2, 1) from 0 to 1
integral_value, _ = scipy.integrate.quad(linear_func, 0, 1, args=(2, 1)) # a=2, b=1
print(f"Integral of 2x + 1 from 0 to 1: {integral_value}") # Expected: [x^2 + x]_0^1 = 2

Output:

Integral of 2x + 1 from 0 to 1: 2.0

Multiple Integrals

For multiple integrals, SciPy provides dblquad, tplquad, and nquad.

Double Integrals with dblquad

The dblquad function computes double integrals: $\int_{a}^{b} \int_{g(x)}^{h(x)} f(x,y) ,dy ,dx$.

Syntax: scipy.integrate.dblquad(func, a, b, gfun, hfun, args=())

• func(y, x, ...): The function to integrate (note the order of arguments: y first, then x).
• a, b: Lower and upper limits for x (outer integral).
• gfun(x), hfun(x): Functions defining the lower and upper limits for y (inner integral), which can depend on x.

Let’s compute $\int_{0}^{1/2} \int_{0}^{\sqrt{1-4y^2}} 16xy ,dx ,dy$. Wait, the original example was $\int_{0}^{1/2} dy \int_{0}^{\sqrt{1-4y^2}} 16xy :dx$. This means $x$ is the inner integral, $y$ is the outer. SciPy’s dblquad integrates $dy dx$ where $y$ is inner, $x$ is outer, if func(y,x). If the problem is $\int_{y_low}^{y_high} dy \int_{x_low(y)}^{x_high(y)} f(x,y) ,dx$, we need to match dblquad’s expectation or swap variables. The original example f = lambda x, y: 16 * x * y with dblquad(f, 0, 0.5, g, h) implies f is f(inner_var, outer_var). Given limits g and h are for the inner variable, and 0, 0.5 are for outer. If g and h depend on y (as in h = lambda y: sqrt(1-4*y**2)), then y must be the outer variable for dblquad. Let’s re-verify dblquad’s argument order for func and limits.

The documentation for dblquad says: func(y, x) is integrated with y from gfun(x) to hfun(x) and x from a to b. So $x$ is outer, $y$ is inner. Our integral $\int_{0}^{1/2} dy \int_{0}^{\sqrt{1-4y^2}} 16xy :dx$ needs to be rephrased if $y$ is outer. Let’s stick to the original integral formulation for the example, assuming $y$ is the outer variable for the user, and $x$ is inner. This means limits for $x$ can depend on $y$. dblquad requires func(inner_var, outer_var). So if $x$ is inner and $y$ is outer, it’s func(x,y). The limits for $x$ are $g(y)$ and $h(y)$, and limits for $y$ are $a$ and $b$. This means: dblquad(lambda x, y: 16*x*y, y_lower, y_upper, lambda y_val: x_lower_limit(y_val), lambda y_val: x_upper_limit(y_val)).

Let’s evaluate $\int_{y=0}^{y=0.5} \int_{x=0}^{x=\sqrt{1-4y^2}} 16xy ,dx ,dy$. Here, $y$ is the outer variable, $x$ is the inner variable.

import numpy as np
import scipy.integrate

# Function to integrate: f(x, y) = 16xy
# dblquad expects func(inner_variable, outer_variable)
# Here, x is inner, y is outer.
integrand = lambda x, y: 16 * x * y

# Limits for the outer variable y
y_lower = 0
y_upper = 0.5

# Limits for the inner variable x (can be functions of y)
x_lower_limit_func = lambda y_val: 0
x_upper_limit_func = lambda y_val: np.sqrt(1 - 4 * y_val**2)

# Perform double integration
integral_value, error_estimate = scipy.integrate.dblquad(
    integrand,
    y_lower,          # Outer variable (y) lower limit
    y_upper,          # Outer variable (y) upper limit
    x_lower_limit_func, # Inner variable (x) lower limit function (of y)
    x_upper_limit_func  # Inner variable (x) upper limit function (of y)
)

print(f"Double integral value: {integral_value}")
print(f"Estimated error: {error_estimate}")

Output:

Double integral value: 0.5000000000000001
Estimated error: 1.7092350012594845e-14

tplquad and nquad

tplquad works similarly for triple integrals, and nquad generalizes this to N-dimensional integrals. For nquad, you provide a list of limit functions or ranges for each dimension. These are powerful but can be computationally intensive for high dimensions.

Integrating Sampled Data

If you have data points $(x_i, y_i)$ rather than a function, you can use methods like the trapezoidal rule (trapz or trapezoid) or Simpson’s rule (simps or simpson).

import numpy as np
from scipy.integrate import trapezoid, simpson

# Sample data points
x_samples = np.linspace(0, np.pi, 10)
y_samples = np.sin(x_samples)

# Integrate using trapezoidal rule
integral_trapz = trapezoid(y_samples, x_samples)
print(f"Integral of sin(x) from 0 to pi (Trapezoidal with 10 points): {integral_trapz}")

# Integrate using Simpson's rule (requires an odd number of samples, or an even number of intervals)
x_samples_simps = np.linspace(0, np.pi, 11) # 11 points = 10 intervals
y_samples_simps = np.sin(x_samples_simps)
integral_simps = simpson(y_samples_simps, x_samples_simps)
print(f"Integral of sin(x) from 0 to pi (Simpson with 11 points): {integral_simps}")

# Analytical result: integral of sin(x) from 0 to pi is [-cos(x)]_0^pi = -(-1) - (-1) = 2
print(f"Analytical result: {2.0}")

Output:

Integral of sin(x) from 0 to pi (Trapezoidal with 10 points): 1.9796508989590044
Integral of sin(x) from 0 to pi (Simpson with 11 points): 2.0001095173150043
Analytical result: 2.0

The scipy.integrate module is a cornerstone of scientific computing in Python. For more complex problems, such as solving differential equations, refer to functions like solve_ivp (Initial Value Problems) and solve_bvp (Boundary Value Problems). Always consult the official SciPy documentation for detailed usage and advanced options.