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

[robve]

Some progress. The following updated and improved method predicts an optimal splitting point d for Exp-Sinh quadratures \( \int_a^\infty f(x)\,dx = \int_a^d f(x)\,dx + \int_d^\infty f(x)\,dx \). The method uses bisection with an exponential scale. The method inspects 2 or 4 points to check if an optimization is applicable. Bisection uses just 8 points to find a close-to-optimal d. This reduces the overhead to only 5.95 evaluation points on average for the 208 integrals tested, requiring 179.38 points each on average with eps=1e-9.

Code:

// return optimized Exp-Sinh integral split point d
double exp_sinh_opt_d(double (*f)(double), double a, double eps, double d) {
  double h2 = f(a + d/2) - f(a + d*2)*4;
  int i = 1, j = 32;     // j=32 is optimal to find r
  if (isfinite(h2) && fabs(h2) > 1e-5) { // if |h2| > 2^-16
    double r, fl, fr, h, s = 0, lfl, lfr, lr = 2;
    do {                 // find max j such that fl and fr are finite
      j /= 2;
      r = 1 << (i + j);
      fl = f(a + d/r);
      fr = f(a + d*r)*r*r;
      h = fl - fr;
    } while (j > 1 && !isfinite(h));
    if (j > 1 && isfinite(h) && sign(h) != sign(h2)) {
      lfl = fl;          // last fl=f(a+d/r)
      lfr = fr;          // last fr=f(a+d*r)*r*r
      do {               // bisect in 4 iterations
        j /= 2;
        r = 1 << (i + j);
        fl = f(a + d/r);
        fr = f(a + d*r)*r*r;
        h = fl - fr;
        if (isfinite(h)) {
          s += fabs(h);  // sum |h| to remove noisy cases
          if (sign(h) == sign(h2)) {
            i += j;      // search right half
          }
          else {         // search left half
            lfl = fl;    // record last fl=f(a+d/r)
            lfr = fr;    // record last fl=f(a+d*r)*r*r
            lr = r;      // record last r
          }
        }
      } while (j > 1);
      if (s > eps) {     // if sum of |h| > eps
        h = lfl - lfr;   // use last fl and fr before the sign change
        r = lr;          // use last r before the sign change
        if (h != 0)      // if last diff != 0, back up r by one step
          r /= 2;
        if (fabs(lfl) < fabs(lfr))
          d /= r;        // move d closer to the finite endpoint
        else
          d *= r;        // move d closer to the infinite endpoint
      }
    }
  }
  return d;
}

// the quad() routine is updated in two branches to use this method:

  else if (isfinite(a)) {
    mode = 1; // Exp-Sinh
    d = exp_sinh_opt_d(f, a, eps, d);
    c = a;
    v = a + d;
  }
  else if (isfinite(b)) {
    mode = 1; // Exp-Sinh
    d = exp_sinh_opt_d(f, b, eps, -d);
    c = b;
    v = b + d;
    sign = -sign;
  }

This new method of optimization is generic and not specifically tailored to the 208 integrands. Integrands that are suitable for Exp-Sinh will benefit, but those that do not converge or very slowly converge such as periodic and oscillating functions will not benefit from an optimized splitting point d. Of the 208 integrals tested, 29 are improved significantly in accuracy and/or number of points evaluated (63 fewer points on average). Of these 208 integrands, 14 are slightly worse (11 more points on average). However, practically there is almost no impact because 8 of these 14 integrands are not suitable for Exp-Sinh quadrature and return a large error with and without the optimized d. The full Exp-Sinh quadrature results with eps=1e-9 are listed in Appendix C.

- Rob