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

[robve]

(04-08-2021 11:20 AM)emece67 Wrote:  With these small changes in the Python code (in fact I copied it from your post above) I do see differences in the results. Usually only in the last digit, but enough to get (a really small) decrement on correct digits and a small but noticeable increment in TNFE (this time mainly in the sin(40πx) case). In 32 out of 160 cases the true absolute error get (slightly) worse. Of them, 14 have mutated from 0 true absolute error, to != 0.

Looks like this falls within machine precision MachEps? Is this just noise or is this significant? Perhaps Python mpmath mp.dps is too small?

One of the reasons why the current WP34S code must be worse compared to the optimized version is simple: the current WP34S code uses two square roots to compute sinh and tanh from cosh. This approach is almost always less accurate compared to computing sinh and tanh directly, which the optimizations perform with the usual hyperbolic formulas based on exp. Furthermore, the fewer arithmetic operations to produce these, the better, as in the optimized code.

Using exp(t) in the inner loop instead of the strength-reduced expt recurrence should be the most accurate. Try it.

However, the expt strength-reduced recurrence is accurate enough when the number of inner loop iterations are limited, which we know is the case. The expt error to exp(t) is at most \( (1+\epsilon)^m-1 \) for \( m \) inner loop iterations. Consider for example \( \epsilon=10^{-15} \) and \( m=100 \) then the error of expt to exp(t) is at most \( 10^{-13} \). But IEEE 754 MachEps is smaller and the integration results are better than that because the inner loop does not exceed 100 iterations (requiring 200 function evaluations) and is typically much lower.

From the 818 integrals tested, for those that do not produce exactly the same integration result, 754 produce smaller errors after optimization of the WP34S code. Almost all of the remaining differ in the last 15'th digit or within MachEps which is more noise so the difference can be discarded as insignificant.

Note that Tanh-Sinh is not suitable for some types of integrals. Integrands can be ill-conditioned, i.e. a small change in x can cause a very large change in y. One can find examples with the "right conditions" when the optimized version differs from the non-optimized version, but this is not at all related to the numerical stability of the (optimized) algorithm.

I would suggest to consider the Michalski & Mosig rule to accelerate the WP34S Tanh-Sinh further with qthsh. The spreadsheet I shared also shows the qthsh has far fewer ALERT than WP34S. For example x^(-3/4)*(1-x)^(-0.25)/(3-2*x) completely fails with the WP34S algorithm after 395 evaluations, but is accurate with qthsh after just 69 evaluations (integral 113 in the list).

Therefore, I am very confident that the suggested optimizations are more than worthwhile to reduce the computational cost significantly (especially on a calculator) and to produce more accurate results in (almost) every case.

Just my 2c.

- Rob