Awesome

The rather modern tanh-sinh quadrature is different from classical Gaussian integration methods in that it doesn't integrate any function exactly, not even polynomials of low degree. Its tremendous usefulness rather comes from the fact that a wide variety of functions, even seemingly difficult ones with (integrable) singularities, can be integrated with arbitrary precision.

Install with

pip install tanh-sinh

and use it like

import tanh_sinh
import numpy as np

val, error_estimate = tanh_sinh.integrate(
    lambda x: np.exp(x) * np.cos(x),
    0,
    np.pi / 2,
    1.0e-14,
    # Optional: Specify first and second derivative for better error estimation
    # f_derivatives={
    #     1: lambda x: np.exp(x) * (np.cos(x) - np.sin(x)),
    #     2: lambda x: -2 * np.exp(x) * np.sin(x),
    # },
)

If you want more digits, use mpmath for arbitrary precision arithmetic:

import tanh_sinh
from mpmath import mp
import sympy

mp.dps = 50

val, error_estimate = tanh_sinh.integrate(
    lambda x: mp.exp(x) * sympy.cos(x),
    0,
    mp.pi / 2,
    1.0e-50,  # !
    mode="mpmath",
)

If the function has a singularity at a boundary, it needs to be shifted such that the singularity is at 0. (This is to avoid round-off errors for points that are very close to the singularity.) If there are singularities at both ends, the function can be shifted both ways and be handed off to integrate_lr; For example, for the function 1 / sqrt(1 - x**2), this gives

import numpy
import tanh_sinh

# def f(x):
#    return 1 / numpy.sqrt(1 - x ** 2)

val, error_estimate = tanh_sinh.integrate_lr(
    lambda x: 1 / numpy.sqrt(-(x**2) + 2 * x),  # = 1 / sqrt(1 - (x-1)**2)
    lambda x: 1 / numpy.sqrt(-(x**2) + 2 * x),  # = 1 / sqrt(1 - (-(x-1))**2)
    2,  # length of the interval
    1.0e-10,
)
print(numpy.pi)
print(val)

3.141592653589793
3.1415926533203944

Awesome

Relevant publications