(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)

The reason is not the alogrithm, but how mpmath numbers work.
Because mpmath does not do overflow / underflow, the integrand may generate noise.
(with IEEE double, generated Inf/NaN will allow the algorithm to fix it ...)

With limits of ±∞, even slight noise will make it hard to converge. (*)

>>> from mpmath import *
>>> f = lambda x: 1 - sqrt(2)*cosh(x)/sqrt(cosh(2*x))
>>> for x in range(15,20): print x, f(x)
...
15 -9.37028232783632e-14
16 -1.26565424807268e-14
17 -1.99840144432528e-15
18 -2.22044604925031e-16
19 -2.22044604925031e-16

We can explicitly remove the noise, then integrate.

>>> from double_exponential import *
>>> f2 = lambda x: 0.0 if abs(x) > 19 else f(x)

>>> double_exponential(f2, -inf, inf)
(mpf('-0.69314718055994773'), mpf('3.8011316316755028e-10'), 61, 4, 2)

This is almost as good as the "sharpened" version of f:

>>> def f3(x): y=1/cosh(2*x); return -y/(1 + sqrt(1+y))
...
>>> double_exponential(f3, -inf, inf)
(mpf('-0.69314718055994584'), mpf('3.8011260805603797e-10'), 69, 4, 2)


(*) Tracing the function calls, errors started at point #10:
old version: f(6704.2155162232693) → mpf('0.0')
new version: f(6704.215516223293) → mpf('-2.2204460492503131e-16')

The old version just get lucky ...

>>> f3(6704.2155162232693)
mpf('-6.1999991241406458e-5824')