(PC-12xx~14xx) qthsh Tanh-Sinh quadrature

[Albert Chan]

(04-08-2021 11:20 AM)emece67 Wrote:  integrals (1 - sqrt(2)*cosh(x)/sqrt(cosh(2*x)), from -inf to +inf) needed one more level.
This particular case worsened dramatically, from 13 correct digits with strength reduction to only 4 with these changes)

I believe there is nothing to worry about.
Integrand suffered from subtraction cancellations.
A slight shift to the nodes/weights (even if more precise) does not translate to better result.

Even for 4 good digits, consider yourself lucky. mpmath quadts only get 3.

>>> f = lambda x: 1-sqrt(2)*cosh(x)/sqrt(cosh(2*x))
>>> quadts(f, [-inf, inf], error=True)
(mpf('-0.69306562029329821'), mpf('0.001'))

Actually, mpmath "cheated", calcuating in higher precision.
See what happen if we pre-calculated sqrt(2), under mp.dps=15 setting.

>>> f2 = lambda x, k=sqrt(2): 1-k*cosh(x)/sqrt(cosh(2*x))
>>> quadts(f2, [-inf, inf], error=True) # failed convergence
(mpf('-1886.3975082035695'), mpf('1.0'))

This is a weakness of tanh-sinh quadrature. It relied on accurate sample points.
Gauss-Legendre quadrature, using curve-fitting of samples, handle "fuzzy" more easily.

>>> quadgl(f2, [-inf, inf], error=True)
(mpf('-0.69314718056249158'), mpf('1.0e-14'))

We could rewrite integrand, to perform better. Let c = cosh(x), d = cosh(2x)
With half-angle formula, we have c*c = (1+d)/2

1 - √(2)*c/√(d) = 1 - √(1+1/d) = -(1/d) / (1 + sqrt(1+1/d))

>>> def f3(x): y=1/cosh(2*x); return -y/(1 + sqrt(1+y))
...
>>> quadts(f3, [-inf, inf], error=1)
(mpf('-0.69314718055994529'), mpf('1.0e-29'))

(04-06-2021 11:36 PM)robve Wrote:  When reusing previous points is better as default for f on [a,b]:

Code:

1/(1+COS(x)^x)        [-1.570796327,1.570796327]    lim x->pi/2+ f(x)=1       lim x->-pi/2- f(x)=0
x^2/(1-COS(x))        [0,1.570796327]        lim x->0- f(x)=0
x*SIN(x)/(1-COS(x))    [0,3.141592654]        lim x->0- f(x)=0

Same issue with all of these. A simple rewrite "sharpen" integrand greatly.

First one had the form 1/(1+u). With identity 1/(1+u) + 1/(1+1/u) = 1/(1+u) + u/(1+u) = 1

∫(1/(1+cos(x)^x), x=-pi/2 .. pi/2) = ∫(1, x=0 .. pi/2) = pi/2

With identity sin(x/2)^2 = (1-cos(x))/2, we can rewrite the others:

x^2/(1-cos(x)) = x^2 / (2*sin(x/2)^2) = (x/sin(x/2))^2 / 2
x*sin(x)/(1-cos(x)) = (x*2*sin(x/2)*cos(x/2) / (2*sin(x/2)^2)) = x / tan(x/2)