# 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](
https://www.gnu.org/software/gsl/doc/html/integration.html#qawf-adaptive-integration-for-fourier-integrals).
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.
Raw data
{
"_id": null,
"home_page": "https://github.com/cimatosa/tsquad",
"name": "tsquad",
"maintainer": "",
"docs_url": null,
"requires_python": ">=3.9,<4.0",
"maintainer_email": "",
"keywords": "numerical-integration,tanh-sinh-method,Fourier-transform,infinite-series",
"author": "Richard Hartmann",
"author_email": "richard.hartmann@tu-dresden.de",
"download_url": "https://files.pythonhosted.org/packages/cd/ef/7f9467d578406f9dacf3518c5d557a5b061ea50a6fbef64ad2cd9b9122dd/tsquad-0.2.5.tar.gz",
"platform": null,
"description": "# tsquad - Numeric Integration Using tanh-sinh Variable Transformation\n\nThe `tsquad` package provides general purpose integration routines.\nFor simple integration tasks this package is as good as the standard routines\navailable e.g. by `scipy`. \nAlthough, the `tsquad` routines can directly integrate handle **complex functions** \n(along the real axes). \nMost strikingly, the `tsquad` method is particularly suited to efficiently handle \n**singularities** by exploiting the *tanh-sinh* variable transformation \n(also known as *double exponential* approach), as nicely described \nhere http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/dhb-tanh-sinh.pdf.\nAs of that, integrals with lower/upper **bounds at infinity** can be treated\nwithin that scheme.\n\nFurthermore, integrals over **oscillatory functions with an infinite boundary** are\nhandled very efficiently by extrapolating the partial sums (obtained from integrating\npiecewise over half the period) using Shanks' transform (Wynn epsilon algorithm).\nThis approach (or at least similar) is also available in the [GSL-library](\nhttps://www.gnu.org/software/gsl/doc/html/integration.html#qawf-adaptive-integration-for-fourier-integrals).\n\nFor a visualization of great speedup of that extrapolation, see the \nexample in `examples_shanks.py`.\n\n# Install\n\n## poetry\n\nUse `poetry` to add `tsquad` to your dependencies with `poetry add tsquad`.\n\nOr edit your `pyproject.toml` file like that\n \n [tool.poetry.dependencies]\n tsquad = \"^0.2.0\"\n\nand run `poetry install`.\n\n## pip\n\nRun `pip install tsquad`.\n\n## git\n\nYou find the source at https://github.com/cimatosa/tsquad.\n\n# Examples\n\n## Singularities\n\nConsider the function `f(x) = 1/x^0.9`. It diverges at `x=0`, however, the integral\nover `[0, 1]` is still finite (`I=10`).\nUsing the tanh-sinh integration scheme allows to efficiently obtained a highly\naccurate numerical result.\n\n >>> import tsquad\n >>> f = lambda x: 1/x**0.9\n >>> tsq = tsquad.QuadTS(f=f)\n >>> tsq.quad(0, 1)\n QuadRes(I=10.000000000000004, err=3.868661047952931e-14, func_calls=73, rec_steps=1)\n\nNote that the function has been called only 73 times.\n\n## Infinite boundary\n\nInfinite boundary condition can be treated efficiently, too.\nThey are mapped to an integral over a finite interval.\nThe resulting singularity poses no difficulty for the tanh-sinh method.\nThe infinite boundary can be specified either as `str` (`'inf'` and `'-inf'`)\nor by `math.inf` as well as `numpy.inf`.\n\nAs example consider the integrand `1/(1 + (x+1)^2)`.\nIts indefinite integral is the `arctan`, so integrating over the whole real\naxes yield `pi`.\n\n >>> import tsquad\n >>> tsq = tsquad.QuadTS(f=f)\n >>> tsq.quad(a='-inf', b='inf')\n QuadRes(I=3.1415926535897643, err=2.3768122480443337e-12, func_calls=528, rec_steps=4)\n >>> import math\n >>> math.pi\n 3.141592653589793\n\n## Fourier integral\n\nIntegrating oscillatory functions needs special care.\nIn particular if the bounds are infinite.\n\nUsually, it is possible to integrate over single periods of the oscillation with \nhigh accuracy.\nSumming up the individual integrals yields results for larger intervals at the \nprice of summing up the errors, too.\nFor very rapidly oscillating functions this might cause some trouble.\n\nFor infinite bounds this sum of partial integrals becomes infinite, often with\nvery slow convergence.\nThe convergence can be accelerated significantly if the terms of the partial \nsum alternate in sign by using the Shanks transform \n(implemented using Wynn's epsilon algorithm).\nFor Fourier integrals (sin, cos or exp(i ...)) the alternating signs\nin the partial sum is realized by summing up finite integrals over **half** the period.\n\nAs example consider the half-sided Fourier integral of the algebraically \ndecaying function `1/(1+1j*tau)^(s+1)`.\nSo we aim to integrate\n\n int_0^inf 1/(1+1j*tau)^(s+1) * exp(1j * w * tau) * d tau .\n\nNote that an analytic expression is possible, which involves the incomplete gamma function\nwith complex-valued second argument.\n\n\n >>> import tsquad\n >>> import mpmath as mp\n >>> s = 0.5\n >>> w = 1\n >>> f = lambda tau, s: 1/(1+1j*tau)**(s+1)\n >>> tsq = tsquad.QuadTS(f=f, args=(s,))\n >>> tsq.quad_Fourier(0, 'inf', w=w)\n QuadRes(I=(1.3040986643460277+0.1523180276515212j), err=None, func_calls=2113, rec_steps=16)\n >>> (-1j*w)**s / 1j**(s+1) * mp.exp(-w) * mp.gammainc(-s, -w-1j*1e-16)\n mpc(real='1.30409866434658424220297858', imag='0.152318027651073881301097481')\n\n\n# Mathematical Details\n\nNote that the tanh-sinh approach is particularly useful when the integrand is singular at `x=0` \nand the integration goes from `a=0` to `b`, i.e.,\n\n int_0^b f(x) dx .\n\n\nTo correctly account for a singularity at a different location consider rewriting the \nintegrand such that the singularity is located at `x=0`.\n\nThe method used here is based on the variable transformation `x -> t` with\n\n x(t) = b/2(1 - g(t))\n g(t) = tanh(pi/2 * sinh(t))\n ->\n x(t) = b / 2 / (e^(pi/2*sinh(t)) cosh(pi/2 sinh(t)))\n\nwhich maps the interval `[0, b]` to `[-inf, inf]`,\n\n\n int_0^b f(x) dx \n = int_-inf^inf f(x(t)) b/2 g'(t) dt \n = int_-inf^inf I(t) dt \n = dt sum_k f_k w_k\n g'(t) = pi cosh(t)/(2 cosh^2(pi/2 sinh(t))) .\n\nIn that way the singularity \"occurs\" at `t=inf`. \nIt has been shown that the error cannot be smaller than `max(|I(t_min)|, |I(t_max)|)`. \nSo `t_min` and `t_max` are chosen such that this limitation is below the desired accuracy.\n\n# MIT-license\n\n\nCopyright 2023 Richard Hartmann\n\nPermission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the \u201cSoftware\u201d), 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:\n\nThe above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.\n\nTHE SOFTWARE IS PROVIDED \u201cAS IS\u201d, 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.\n",
"bugtrack_url": null,
"license": "MIT",
"summary": "numeric quadrature using the Tanh-Sinh variable transformation",
"version": "0.2.5",
"project_urls": {
"Homepage": "https://github.com/cimatosa/tsquad",
"Repository": "https://github.com/cimatosa/tsquad"
},
"split_keywords": [
"numerical-integration",
"tanh-sinh-method",
"fourier-transform",
"infinite-series"
],
"urls": [
{
"comment_text": "",
"digests": {
"blake2b_256": "2fe36e2629544f4677ef47d2788cd1d401c608b423aba5354ce0e463cc6cc975",
"md5": "7ab849d0dec91c80710b9ef1af9cd7ba",
"sha256": "52f1e3b2a9cdd2347d4e1013d8339060d4b4537b9afd79a37ef15ceb0f308ad9"
},
"downloads": -1,
"filename": "tsquad-0.2.5-py3-none-any.whl",
"has_sig": false,
"md5_digest": "7ab849d0dec91c80710b9ef1af9cd7ba",
"packagetype": "bdist_wheel",
"python_version": "py3",
"requires_python": ">=3.9,<4.0",
"size": 193812,
"upload_time": "2023-11-16T09:22:29",
"upload_time_iso_8601": "2023-11-16T09:22:29.846079Z",
"url": "https://files.pythonhosted.org/packages/2f/e3/6e2629544f4677ef47d2788cd1d401c608b423aba5354ce0e463cc6cc975/tsquad-0.2.5-py3-none-any.whl",
"yanked": false,
"yanked_reason": null
},
{
"comment_text": "",
"digests": {
"blake2b_256": "cdef7f9467d578406f9dacf3518c5d557a5b061ea50a6fbef64ad2cd9b9122dd",
"md5": "b631ce3b9ba4f1b7805838e561a12768",
"sha256": "cf2909fd44e3c07e5c548156d91faed487a9baa263fe4aba5999ba88423fe8b1"
},
"downloads": -1,
"filename": "tsquad-0.2.5.tar.gz",
"has_sig": false,
"md5_digest": "b631ce3b9ba4f1b7805838e561a12768",
"packagetype": "sdist",
"python_version": "source",
"requires_python": ">=3.9,<4.0",
"size": 192815,
"upload_time": "2023-11-16T09:22:32",
"upload_time_iso_8601": "2023-11-16T09:22:32.764785Z",
"url": "https://files.pythonhosted.org/packages/cd/ef/7f9467d578406f9dacf3518c5d557a5b061ea50a6fbef64ad2cd9b9122dd/tsquad-0.2.5.tar.gz",
"yanked": false,
"yanked_reason": null
}
],
"upload_time": "2023-11-16 09:22:32",
"github": true,
"gitlab": false,
"bitbucket": false,
"codeberg": false,
"github_user": "cimatosa",
"github_project": "tsquad",
"travis_ci": false,
"coveralls": false,
"github_actions": false,
"lcname": "tsquad"
}
JSON
