Accuracy of Integral with epsilon

[robve]

(10-12-2021 04:47 PM)Claudio L. Wrote:  So I guess the trick is as long as you don't ask the integration algorithm for an error smaller than the system precision less 9 digits, you'll be fine and get a nice result, with an error in the order of (precision-9), regardless of the tolerance you requested.

Catastrophic cancellation errors exponentially increase towards negative infinity for this integrand.

I may be going off on a tangent in this thread, but what other sources are there for noise, besides cancellation and cosh/sinh/exp roundoff/differences? Can the integration method be noisy? I'd expect that to be insignificant. It may depend on the accuracy of the weight calculations. Simple weights pose no problem, e.g. Simpson, Romberg are exact for (low order) polynomials anyway. Weights derived from table-based methods, such as Gauss-Kronrod may collect small errors in the weight calculations. I was curious to test this with a simple constant function \( \int_{-10}^{10} 2\,dx \) to integrate with eps=10^-6 tolerance (or about, what works to get comparable results w.r.t. number of evaluation points):

- Gauss-Kronrod G10K21 with IEEE 754 doubles = 40.000000000000007.
- WP-34S Tanh-Sinh C version with IEEE 754 double: 39.999999999987885 (23 points)
- qthsh Tanh-Sinh with IEEE 754 double: 40.000000000000057 (29 points)
- Boost Math Tanh-Sinh with IEEE 754 double: 40.000000000001464 (25 points)

Tanh-Sinh methods may suffer "weight issues" due to their point distributions. There are also numerical differences in the implementations. WP-34S and qthsh use the more efficient Michalski-Mosig quadrature method. Boost Math uses the standard Takahashi-Mori Tanh-Sinh quadrature method. And there are other internal differences. WP-34S computes the weights as follows (see forum thread):

Code:

      double ch = cosh(t);
      double ecs = cosh(M_PI/2 * sqrt(ch*ch - 1)); // = cosh(pi/2*sinh(t))
      double w = 1/(ecs*ecs);
      double r = sqrt(ecs*ecs - 1)/ecs; // = sinh(pi/2*sinh(t))/cosh(pi/2*sinh(t)) = tanh(pi/2*sinh(t))

Whereas qthsh computes the weights as follows:

Code:

      double u = exp(1/t-t); // = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2
      double r = 2*u/(1+u);  // = 1 - tanh(sinh(j*h))
      double w = (t+1/t)*r/(1+u); // = cosh(j*h)/cosh(sinh(j*h))^2

The noise here is still rather insignificant compared to cancellation and diminishes with smaller given tolerances toward 10^-15 with more points evaluated by Tanh-Sinh.

- Rob