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

[robve]

(03-28-2021 06:59 PM)Albert Chan Wrote:  Some issues I noticed for Tanh-Sinh quadrature.

Thanks for sharing your findings! It's great to play with these methods to see their practical strengths and weaknesses.

There are limitations to the applicability of the method that are generally known (see e.g. Boost Math Tanh-Sinh for a decent explanation) and issues that are caused in practice using this method. But I am sure I don't need to tell you that

The implementation in C and BASIC I wrote are based on several other implementations (see links). I have not optimized the code or tried to fix problems that I already saw in the code early on (more about this further below).

(03-28-2021 06:59 PM)Albert Chan Wrote:  >>> from mpmath import *
>>> f = lambda x: 1/sqrt(-ln(x))

The crash with this function does not happen in the C and BASIC versions I've worked on (here output with k, error, and points for debugging):

k=3 err=2.29154e-08 pts=32
qthsh = 1.77245384129412

However, it happens to the same function after simply reversing the x range over [0,1]:

double f(double x) { return 1/sqrt(-log(1-x)); }

This crashes at k=2:

k=2 err=nan pts=17
int qthsh = -inf

When prepping my implementation, I noticed the following conditions in some existing implementations of the method I compared to look for certain details such as this, which I've adopted in the exact same way to avoid hitting the endpoints:

if (a+x > a)
  fp = f(a+x);
if (b-x < b)
  fm = f(b-x);

Do you spot the problem too?





If a=0 or b=0 then x can be arbitrarily small.

On the other hand, if a=1 or b=1 then this "fix" works the way it was intended.

This happens in the "best Tanh-Sinh implementation on the planet".

Further, perhaps in the Python and Lua implementations there is no guard at all, since your case failed either way?

Fortunately, we can simply modify the code to something like this:

if (1+x > 1) {
  fp = f(a+x);
  fm = f(b-x);
}

Hopefully the C/C++ compiler does not optimize the condition to x > 0, which is algebraically identical, but numerically destroys the benefits of having this test (because x > 0 is always true). The type of optimizations were a serious problem in the past. Compiler optimizations wrongly assumed that floating point arithmetic is associative and commutative. Modern C/C++ compilers should not mess this up.

With this simple change I propose, it works as expected:

k=3 err=2.29154e-08 pts=29
int qthsh = 1.77245384129412

I can't think of a problem with this change, though the value 1 could be picked smaller such as eps or MachEps for example. It depends on how close you allow x to approach the endpoints.

I also noticed in the public implementations of the method with this incomplete fix that fp and fm are simply the old f(a+x) and f(b-x) values i.e. at neighboring points closer to the center of the x range. I believe that should be fine, though I doubt this was researched thoroughly. It has an ad-hoc stink to it.

I wonder how many of these public Tanh-Sinh implementations are numerically flawed.

- Rob