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

[robve]

I made two notable improvements to qthsh and quad that are documented in the updated article. The first post of this thread is also updated.

The first improvement is a minor code optimization. The previous qthsh and quad code fragment with the convergence test is:

Code:

    v = 2*s;
    s += p;
    ++k;
  } while (fabs(v-s) > tol*fabs(s) && k <= n);
  e = fabs(v-s)/(fabs(s)+eps);

That is, we compute \( s_{k+1} = s_k + p \) and the loop terminates when \( |2s_k-s_{k+1}| \le {\it tol}|s_{k+1}| \) .

Note that \( 2s_k-s_{k+1} = 2s_k-s_k-p = s_k-p \) and we can make the code more elegant and slightly more optimal to avoid cancellation:

Code:

    v = s-p;
    s += p;
    ++k;
  } while (fabs(v) > 10*eps*fabs(s) && k <= n);
  e = fabs(v)/(fabs(s)+eps);

where variable v is now used to store the estimated absolute error and e is the estimated relative error with +eps in the denominator to avoid catastrophic blow up when the integral value s is close to zero.

The second improvement is an addition to the article on pages 20 and 21 in section "Improved accuracy by adjusting the tolerance threshold". I felt that this section needed more data to justify the tol=10*eps choice and also needed quantitative comparisons of the estimated and actual errors for the 21 showcased integrals and the 818 test integrals.

A new table is added to the article showing the estimated relative error e to the actual relative error. The table also shows the estimated absolute error v and the actual absulute error. Notable in this table (not shown here) is that the estimated relative error e almost always lies within one order of magnitude of the actual relative error, i.e. e estimates the accuracy of the result which may be off by one digit, but typically not more.

I verified the accuracy of the estimated relative error e for a range of eps values 1e-5 to 1e-10 for the 818 test integrals:

\[
\begin{array}{|r|c|c|c|c|c|c|}
\hline
{\it eps}= & 10^{-5} & 10^{-6} & 10^{-7} & 10^{-8} & 10^{-9} & 10^{-10} \\
\hline
\mbox{average number of evaluations} & 57 & 69.4 & 81.7 & 93.4 & 106 & 132 \\
\hline
\mbox{actual relative error} \le {\it eps} & 771 & 766 & 755 & 707 & 559 & 182 \\
\hline
\mbox{estimated rel.err. is within 10}\times & 809 & 807 & 799 & 780 & 690 & 385 \\
\hline
\end{array}
\]

This table shows the number of cases when the actual relative error≤eps on the second row and the number of cases when the estimated relative error is within 10× of the actual relative error on the third row, that is, e≤10eps if the actual relative error≤eps (the method converged) and e≥eps/10 if the actual relative error≥eps (the method failed to converge). The table suggests that e can be a reliable estimate of the relative error for eps up to 1e-9. There is a sharp drop-off in the number of correct cases at eps=1e-10, i.e. not many results meet the desired 1e-10 accuracy, but rather meet a 1e-9 accuracy.

This makes sense when comparing qthsh to the VB code and the WP-34S code:

\[
\begin{array}{|r|c|c|c|c|c|c|}
\hline
& \mbox{qthsh} & \mbox{VB} & \mbox{WP-34S} \\
\hline
\mbox{average number of evaluations} & 106 & 97.8 & 62 \\
\hline
\mbox{actual relative error} \le {\it eps} & 559 & 549 & 233 \\
\hline
\mbox{estimated rel.err. is within 10}\times & 690 & 634 & 593 \\
\hline
\end{array}
\]

It should be noted that the VB code succeeds for 796 integrands instead of 814 for the others. This is due to NaN issues that aren’t caught in the VB code. Overall, qthsh is more accurate, but at a cost of a slightly higher average number of function evaluations.

To reduce the number of function evaluations in qthsh to 93.9 on average (for the 818 test cases), we can increase the tolerance, much like the VB code uses a higher tolerance to get the fewer 97.8 evaluations on average. However, for qthsh we use an optimal tol=160*eps that I empirically found to be optimal and is also largely independent of eps the integral. We also adjust the estimated relative error e accordingly:

Code:


double tol = 160*eps;
...
e = fabs(v)/(16*fabs(s)+eps);

This reduces the function evaluations by over 10% with a minimal impact on accuracy:

\[
\begin{array}{|r|c|c|c|c|c|c|}
\hline
{\it eps}= & 10^{-5} & 10^{-6} & 10^{-7} & 10^{-8} & 10^{-9} & 10^{-10} \\
\hline
\mbox{average number of evaluations} & 47.5 & 59.2 & 68.9 & 82 & 93.9 & 106 \\
\hline
\mbox{actual relative error} \le {\it eps} & 754 & 767 & 753 & 715 & 559 & 183 \\
\hline
\mbox{estimated rel.err. is within 10}\times & 805 & 803 & 796 & 777 & 682 & 352 \\
\hline
\end{array}
\]

By comparison, the 93.9 average number of evaluations of qthsh for eps=1e-9 is lower than the VB code. Also qthsh has a higher number of accurate results within the actual and estimated relative errors.

PS. feel free to rerun these tests with the code and integrands in qthsh2.zip

- Rob