numeric quadrature using the Tanh-Sinh variable transformation
Project description
tsquad - Numeric Integration Using tanh-sinh Variable Transformation
The tsquad
package provides general purpose integration routines.
For simple integration tasks this package is as good as the standard routines
available e.g. by scipy
.
Although, the tsquad
routines can directly integrate handle complex functions
(along the real axes).
Most strikingly, the tsquad
method is particularly suited to efficiently handle
singularities by exploiting the tanh-sinh variable transformation
(also known as double exponential approach), as nicely described
here http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/dhb-tanh-sinh.pdf.
As of that, integrals with lower/upper bounds at infinity can be treated
within that scheme.
Furthermore, integrals over oscillatory functions with an infinite boundary are handled very efficiently by extrapolating the partial sums (obtained from integrating piecewise over half the period) using Shanks' transform (Wynn epsilon algorithm). This approach (or at least similar) is also available in the GSL-library.
For a visualization of great speedup of that extrapolation, see the
example in examples_shanks.py
.
Install
poetry
Use poetry
to add tsquad
to your dependencies with poetry add tsquad
.
Or edit your pyproject.toml
file like that
[tool.poetry.dependencies]
tsquad = "^0.2.0"
and run poetry install
.
pip
Run pip install tsquad
.
git
You find the source at https://github.com/cimatosa/tsquad.
Examples
Singularities
Consider the function f(x) = 1/x^0.9
. It diverges at x=0
, however, the integral
over [0, 1]
is still finite (I=10
).
Using the tanh-sinh integration scheme allows to efficiently obtained a highly
accurate numerical result.
>>> import tsquad
>>> f = lambda x: 1/x**0.9
>>> tsq = tsquad.QuadTS(f=f)
>>> tsq.quad(0, 1)
QuadRes(I=10.000000000000004, err=3.868661047952931e-14, func_calls=73, rec_steps=1)
Note that the function has been called only 73 times.
Infinite boundary
Infinite boundary condition can be treated efficiently, too.
They are mapped to an integral over a finite interval.
The resulting singularity poses no difficulty for the tanh-sinh method.
The infinite boundary can be specified either as str
('inf'
and '-inf'
)
or by math.inf
as well as numpy.inf
.
As example consider the integrand 1/(1 + (x+1)^2)
.
Its indefinite integral is the arctan
, so integrating over the whole real
axes yield pi
.
>>> import tsquad
>>> tsq = tsquad.QuadTS(f=f)
>>> tsq.quad(a='-inf', b='inf')
QuadRes(I=3.1415926535897643, err=2.3768122480443337e-12, func_calls=528, rec_steps=4)
>>> import math
>>> math.pi
3.141592653589793
Fourier integral
Integrating oscillatory functions needs special care. In particular if the bounds are infinite.
Usually, it is possible to integrate over single periods of the oscillation with high accuracy. Summing up the individual integrals yields results for larger intervals at the price of summing up the errors, too. For very rapidly oscillating functions this might cause some trouble.
For infinite bounds this sum of partial integrals becomes infinite, often with very slow convergence. The convergence can be accelerated significantly if the terms of the partial sum alternate in sign by using the Shanks transform (implemented using Wynn's epsilon algorithm). For Fourier integrals (sin, cos or exp(i ...)) the alternating signs in the partial sum is realized by summing up finite integrals over half the period.
As example consider the half-sided Fourier integral of the algebraically
decaying function 1/(1+1j*tau)^(s+1)
.
So we aim to integrate
int_0^inf 1/(1+1j*tau)^(s+1) * exp(1j * w * tau) * d tau .
Note that an analytic expression is possible, which involves the incomplete gamma function with complex-valued second argument.
>>> import tsquad
>>> import mpmath as mp
>>> s = 0.5
>>> w = 1
>>> f = lambda tau, s: 1/(1+1j*tau)**(s+1)
>>> tsq = tsquad.QuadTS(f=f, args=(s,))
>>> tsq.quad_Fourier(0, 'inf', w=w)
QuadRes(I=(1.3040986643460277+0.1523180276515212j), err=None, func_calls=2113, rec_steps=16)
>>> (-1j*w)**s / 1j**(s+1) * mp.exp(-w) * mp.gammainc(-s, -w-1j*1e-16)
mpc(real='1.30409866434658424220297858', imag='0.152318027651073881301097481')
Mathematical Details
Note that the tanh-sinh approach is particularly useful when the integrand is singular at x=0
and the integration goes from a=0
to b
, i.e.,
int_0^b f(x) dx .
To correctly account for a singularity at a different location consider rewriting the
integrand such that the singularity is located at x=0
.
The method used here is based on the variable transformation x -> t
with
x(t) = b/2(1 - g(t))
g(t) = tanh(pi/2 * sinh(t))
->
x(t) = b / 2 / (e^(pi/2*sinh(t)) cosh(pi/2 sinh(t)))
which maps the interval [0, b]
to [-inf, inf]
,
int_0^b f(x) dx
= int_-inf^inf f(x(t)) b/2 g'(t) dt
= int_-inf^inf I(t) dt
= dt sum_k f_k w_k
g'(t) = pi cosh(t)/(2 cosh^2(pi/2 sinh(t))) .
In that way the singularity "occurs" at t=inf
.
It has been shown that the error cannot be smaller than max(|I(t_min)|, |I(t_max)|)
.
So t_min
and t_max
are chosen such that this limitation is below the desired accuracy.
MIT-license
Copyright 2023 Richard Hartmann
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
Project details
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.