Tanh-sinh quadrature(The Double Exponential Formulas)

本文主要是介绍Tanh-sinh quadrature(The Double Exponential Formulas)，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

Tanh-sinh quadrature(The Double Exponential Formulas for Numerical Integration)

见维基百科
https://en.m.wikipedia.org/wiki/Tanh-sinh_quadrature

见Tanh-sinh quadrature numerical integration method converging to wrong value
https://www.e-learn.cn/topic/2848897

见Stackoverflow
https://stackoverflow.com/questions/24986588/tanh-sinh-quadrature-numerical-integration-method-converging-to-wrong-value

参考文献
<T Hidetosi，M Masatake, Double Exponential Formulas for Numerical Integration,DOI:10.2977/prims/1195192451>

参考文献
<TakuyaOoura, An IMT-type quadrature formula with the same asymptotic performance as the DE formula,DOI:10.1016/j.cam.2007.01.002>

测试如下（采用Python高精度浮点算法包 mpmath）

# -*- coding: utf-8 -*-
"""
Created on Fri Oct 22 20:14:43 2021I=Int_{a}^{b} f(z)dz其中a<b为任意实数,z属于[a,b]为任意区间首先将区间[a,b]变换到 t属于 [-1,1]z=(b-a)*t/2 +(b+a)/2则I=((b-a)/2)*Int_{-1}^{1} g(t)dt其中g(t)=f((b-a)*t/2 +(b+a)/2)采用双指数数值积分或Tanh-sinh 数值积分现将积分区间变换到[-inf,inf],则t=tanh((pi/2)*sinh(nh))w=cosh(nh)/cosh((pi/2)*sinh(nh))**2I_h=h*pi/2*sum_{n=-N}^{N} f(t)*w具体可以参考文献<T Hidetosi，M Masatake, Double Exponential Formulas for Numerical          Integration,DOI:10.2977/prims/1195192451>"""
from mpmath import mp,mpf,cosh,sinh,pi,tanh,sqrt
if(1):mp.dps = 100h = mpf(2**-12);#def weights(k):num = mpf(0.5)*h*pi*cosh(k*h)den = cosh(mpf(0.5)*pi*sinh(k*h))**2return (num/den)def abscissas(k,a,b):f=((b-a)/2)*(tanh(mpf(0.5)*pi*sinh(k*h)))+(b+a)/2return fdef f(x):return 1/sqrt(1 - mpf(x)**2)N = 20000#通过作线性变换，求任意区间的数值积分a=0b=1result = 0for k in range(-N, N+1):result = result + weights(k)*f(abscissas(k,a,b))*(b-a)/2print(result)print(result - pi)
'''```
#可以发现对N的要求很高，且很慢
#N=20000
3.141592653589793238462643383279502884197169395623305210054471789091557055855748154170301951289498386
-3.751800610920472803216259350430460844457732874052618682441090144344372471319795201134275503228835472e-45
#N=2000
1.449489205637340355829215405535223522203982454439832629403540056142008755610237611478998353099179717
-1.692103447952452882633427977744279361993186944935273191571404536165807650675971387149036472242937351#N=200
0.1536915064631849189177884604412165536645462018136925450754802115004847789524209471365251505235660057
-2.987901147126608319544854922838286330532623197561413275899464380807331627333788051491509674818551062
'''