Single Integration with SciPy's quad() Function

The quad() function from the scipy.integrate module provides a robust and versatile method for performing single-variable numerical integration, commonly known as definite integration. It is a cornerstone in scientific computing for evaluating integrals of continuous functions over specified intervals. quad() leverages adaptive quadrature techniques, ensuring both high accuracy and computational efficiency.

What is SciPy's quad() Function?

The scipy.integrate.quad() function is designed to compute the definite integral of a given function over a specified range. It returns two values: the numerical result of the integration and an estimate of the absolute error associated with the computation.

Mathematical Representation:

The quad() function evaluates integrals of the form:

$$ \int_{a}^{b} f(x) , dx $$

It is particularly well-suited for smooth functions with finite or infinite integration limits and can handle both real and complex-valued integrands.

Syntax of quad() Function

scipy.integrate.quad(func, a, b, args=(), full_output=0, epsabs=1.49e-08, epsrel=1.49e-08,
                     limit=50, points=None, weight=None, wvar=None, wopts=None,
                     maxp1=50, limlst=50, complex_func=False)

Important Parameters:

• func: The function to be integrated. This function must accept a single argument (the variable of integration) and return a scalar value.
• a: The lower limit of integration (a float).
• b: The upper limit of integration (a float).
• args (optional): A tuple of extra arguments to be passed to the func.
• epsabs: The absolute error tolerance. The algorithm stops when the absolute error is less than this value (default: 1.49e-08).
• epsrel: The relative error tolerance. The algorithm stops when the relative error is less than this value (default: 1.49e-08).
• limit: The maximum number of subintervals to use in the adaptive quadrature (default: 50).
• full_output: If set to 1, the function returns additional information about the convergence of the integration process.
• complex_func: Set to True if the function (func) returns complex numbers. This allows for the integration of complex-valued functions.

Examples of Using quad()

1. Single Integration of sin(x) from 0 to π

This example demonstrates the basic usage of quad() to integrate the sine function.

import numpy as np
from scipy import integrate

def f(x):
    """The function to integrate: sin(x)."""
    return np.sin(x)

# Integrate sin(x) from 0 to pi
result, error = integrate.quad(f, 0, np.pi)

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

Output:

Integral Result: 2.0
Estimated Error: 2.220446049250313e-14

2. Handling Infinite Limits in Integration

The quad() function can elegantly handle improper integrals with infinite limits by using np.inf or -np.inf for the integration bounds.

Example: Integrating e^(-x) from 0 to ∞

import numpy as np
from scipy import integrate

def g(x):
    """The function to integrate: e^(-x)."""
    return np.exp(-x)

# Integrate e^(-x) from 0 to infinity
result, error = integrate.quad(g, 0, np.inf)

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

Output:

Integral Result: 1.0
Estimated Error: 1.1102230246251565e-14

3. Controlling Error Tolerance in Integration

You can refine the accuracy of the integration or potentially reduce computation time by adjusting the epsabs and epsrel parameters.

Example: Custom Tolerances for e^(-x)

import numpy as np
from scipy import integrate

def g(x):
    """The function to integrate: e^(-x)."""
    return np.exp(-x)

# Integrate e^(-x) from 0 to infinity with custom tolerances
result, error = integrate.quad(g, 0, np.inf, epsabs=1e-9, epsrel=1e-9)

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

Output:

Integral Result: 1.0
Estimated Error: 7.105427357601002e-10

Note: The specific error value may vary slightly based on the execution environment and internal algorithm details, but it should be within the specified tolerance.

4. Integrating Complex-Valued Functions

To integrate functions that return complex numbers, set the complex_func parameter to True. SciPy will then handle the integration of both the real and imaginary parts separately.

Example: Integrating exp(i·x) from 0 to π

Here, we integrate $e^{ix}$ from $0$ to $\pi$. The analytical solution is $\frac{e^{i\pi} - e^{i0}}{i} = \frac{-1 - 1}{i} = \frac{-2}{i} = 2i$.

import numpy as np
from scipy import integrate

def complex_function(x):
    """The function to integrate: exp(i*x)."""
    return np.exp(1j * x)

# Integrate exp(i*x) from 0 to pi, specifying it's a complex function
result, error = integrate.quad(complex_function, 0, np.pi, complex_func=True)

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

Output:

Integral Result: (3.6775933888827275e-17+2j)
Estimated Error: (2.2102239425853306e-14+2.220446049250313e-14j)

Note: The output is very close to $2i$. The small real part is due to numerical precision, and the error estimate reflects this.

Conclusion

The scipy.integrate.quad() function is an indispensable tool for one-dimensional numerical integration. Its key strengths include:

• Versatile Limits: Supports integration over both finite and infinite intervals.
• Adjustable Precision: Allows users to control the accuracy of the results through epsabs and epsrel.
• Complex Number Support: Capable of integrating complex-valued functions.
• Comprehensive Output: Returns both the integrated value and an estimate of the error.
• Adaptive Quadrature: Employs advanced algorithms for efficient and accurate computations.

Whether you are tackling problems in physics, engineering, statistics, economics, or any field requiring numerical integration, quad() provides a reliable and powerful solution, especially when analytical integration is not feasible or too complex.