项目作者: nschloe

项目描述 :
tanh-sinh quadrature for Python
高级语言: Python
项目地址: git://github.com/nschloe/tanh_sinh.git
创建时间: 2020-05-16T09:28:58Z
项目社区:https://github.com/nschloe/tanh_sinh

开源协议:GNU General Public License v3.0

下载



logo

PyPi Version
PyPI pyversions
GitHub stars
PyPi downloads

Discord

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

  1. pip install tanh-sinh

and use it like

  1. import tanh_sinh
  2. import numpy as np
  3. val, error_estimate = tanh_sinh.integrate(
  4. lambda x: np.exp(x) * np.cos(x),
  5. # Optional: Specify the function with its first and second derivative for
  6. # better error estimation
  7. # (
  8. # lambda x: np.exp(x) * np.cos(x),
  9. # lambda x: np.exp(x) * (np.cos(x) - np.sin(x)),
  10. # lambda x: -2 * np.exp(x) * np.sin(x),
  11. # )
  12. 0,
  13. np.pi / 2,
  14. 1.0e-14,
  15. )

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

  1. import tanh_sinh
  2. from mpmath import mp
  3. import sympy
  4. mp.dps = 50
  5. val, error_estimate = tanh_sinh.integrate(
  6. lambda x: mp.exp(x) * sympy.cos(x),
  7. 0,
  8. mp.pi / 2,
  9. 1.0e-50, # !
  10. mode="mpmath",
  11. )

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

  1. import numpy
  2. import tanh_sinh
  3. # def f(x):
  4. # return 1 / numpy.sqrt(1 - x ** 2)
  5. val, error_estimate = tanh_sinh.integrate_lr(
  6. lambda x: 1 / numpy.sqrt(-(x**2) + 2 * x), # = 1 / sqrt(1 - (x-1)**2)
  7. lambda x: 1 / numpy.sqrt(-(x**2) + 2 * x), # = 1 / sqrt(1 - (-(x-1))**2)
  8. 2, # length of the interval
  9. 1.0e-10,
  10. )
  11. print(numpy.pi)
  12. print(val)
  1. 3.141592653589793
  2. 3.1415926533203944

Relevant publications