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

[robve]

Albert, thanks for sharing additional examples and details.

Let me first say that I'm still "playing around" with this method using several implementations to compare (but I do not have much time beyond weekends.)

Tanh-Sinh also crashes in Boost Math library for \( \int_0^1 1/\sqrt{-\log(1-x)} \,dx \) with the Abort message:

libc++abi.dylib: terminating with uncaught exception of type boost::wrapexcept<boost::math::evaluation_error>: Error in function tanh_sinh<double>::integrate: The tanh_sinh quadrature evaluated your function at a singular point and got -inf. Please narrow the bounds of integration or check your function for singularities.
Abort

but does not crash for \( \int_0^1 1/\sqrt{-\log x} \) which is the symmetric equivalent(!), i.e. x approaching 0 is more problematic than approaching 1. But the reason is probably different, since there is no explicit guard against evaluating at a+x or b-x when x is too small (IMHO a practice that can have problems of its own, i.e. to ignore function values may produce inaccurate error estimates and less precise integrals).

(03-29-2021 02:28 AM)Albert Chan Wrote:  The problem is uf = lambda u: 6*u*(1-u) * f(u*u*(3-2*u))

∫(uf(u), u=0..1) will touch f(x=1) = ∞, much earlier than expected.

This function works in
Boost Math: Tanh-sinh=0.5 err=3.66374e-14 points=74
C implementation in this post: qthsh=0.5 k=4 err=1.55431e-15 points=57

The Boost Math Tanh-Sinh is a rather straightforward implementation. By comparison, the WP 43S code is more clever, but both are not always robust to functions that have problems in close proximity to the endpoints.

The Tanh-Sinh transformation does produce more points closer to the endpoints and performs particularly well for U-shaped functions or L-shaped on the domain of integration. Functions that are sensitive to errors in close proximity to the endpoints may not be usable with this method, depending on how many levels are (internally) used and if the implementation has guards that attempt to prevent this (if there are any at all). For the latter case, the "fix" in the code I looked at is likely a response to this problem (but still fails).

(03-29-2021 02:28 AM)Albert Chan Wrote:  ∫(f(x), x=0..1) did not fail in Python mpmath either.

A quick look at mpmath's quadrature internals reveals some interesting details. The implementation produces a list of (w,x) with weights w and nodes x to sum, see
260     class TanhSinh(QuadratureRule):

The mpf precision is set extra high:

311     def calc_nodes(self, degree, prec, verbose=False):
...
329        extra = 20
330        ctx.prec += extra

After producing the Tanh-Sinh points in calc_nodes, the function to integrate is evaluated at all x in the node list, regardless how close x is to the endpoints:

308        S += self.ctx.fdot((w,f(x)) for (x,w) in nodes)

where fdot ultiplies the weight w with f(x) represented in mpf, e.g. when w and x are mpf reals:

892     def fdot(ctx, A, B=None, conjugate=False):
944        real.append(mpf_mul(a._mpf_, b._mpf_))
...
972        s = mpf_sum(real, prec, rnd)
...
976        s = ctx.make_mpf(s)
977        return s

It looks like there is no error in this case due to higher precision. More luck than wisdom perhaps?

Let me add that Boost Math Tanh-Sinh does appear to offer a fix by passing the distance to the endpoint to the integrator:

sqrt(x / (1 - x * x))

This an example of an integral that has all its area close to a non-zero endpoint, the problem here is that the function being integrated returns "garbage" values for x very close to 1. We can easily fix this issue by passing a 2 argument functor to the integrator: the second argument gives the distance to the nearest endpoint, and we can use that information to return accurate values, and thus fix the integral calculation.

As stated in the table with example cases spanning simple to problematic.

And, of course:

* As always, there are caveats: For instance, if the function you want to integrate is not holomorphic on the unit disk,
* then the rapid convergence will be spoiled. In this case, a more appropriate quadrature is (say) Romberg, which does not
* require the function to be holomorphic, only differentiable up to some order.