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

[robve]

(04-18-2021 01:30 AM)robve Wrote:  In your example the left side's magnitude f(a+d/v)/v is larger than f(a+d*v)*v which means that the sign changes for a different reason. In my test this shows that d should be small (because the test also takes the magnitudes of the difference into account). So your example is not a counter example of my earlier post on why this might work.

Here are the numbers.

For your example exp(-0.01/x)/x^2 the change is not due to the right side's magnitude, but the left side:
d diff (left-right)
2 diff=2.92578 (3.92079 - 0.995012)
4 diff=14.3751 (15.3726 - 0.997503)
8 diff=58.0807 (59.0794 - 0.998751)
16 diff=217.149 (218.149 - 0.999375)
32 diff=742.577 (743.577 - 0.999688)
64 diff=2158.79 (2159.79 - 0.999844)
128 diff=4554.36 (4555.36 - 0.999922)
256 diff=5065.24 (5066.24 - 0.999961)
512 diff=1565.58 (1566.58 - 0.99998)
1024 diff=36.4476 (37.4476 - 0.99999)
2048 diff=-0.994646 (0.00534945 - 0.999995)
4096 diff=-0.999998 (2.72909e-11 - 0.999998)
8192 diff=-0.999999 (1.77573e-28 - 0.999999)
16384 diff=-0.999999 (1.87946e-63 - 0.999999)
32768 diff=-1 (5.26361e-134 - 1)
65536 diff=-1 (1.03212e-275 - 1)

For exp(-.01*x) the change is due to the right side's magnitude:
d diff (left-right)
2 diff=-2.92578 (0.995012 - 3.92079)
4 diff=-14.3751 (0.997503 - 15.3726)
8 diff=-58.0807 (0.998751 - 59.0794)
16 diff=-217.149 (0.999375 - 218.149)
32 diff=-742.577 (0.999688 - 743.577)
64 diff=-2158.79 (0.999844 - 2159.79)
128 diff=-4554.36 (0.999922 - 4555.36)
256 diff=-5065.24 (0.999961 - 5066.24)
512 diff=-1565.58 (0.99998 - 1566.58)
1024 diff=-36.4476 (0.99999 - 37.4476)
2048 diff=0.994646 (0.999995 - 0.00534945)
4096 diff=0.999998 (0.999998 - 2.72909e-11)
8192 diff=0.999999 (0.999999 - 1.77573e-28)
16384 diff=0.999999 (0.999999 - 1.87946e-63)
32768 diff=1 (1 - 5.26361e-134)
65536 diff=1 (1 - 1.03212e-275)

Part of tweaking the test to get an optimal d is to observe the cliff: there is a sharp drop in magnitude in the right part, from 1566.58 to 37.4476 while the left side is still smooth, followed by a sign change. This suggests an optimal d is somewhere near 512-1024.

Edit: I just now realized I had changed the test earlier today with respect to the magnitude check included in the test as a filter. I had not published this change or said anything about it. To produce the numbers shown above:

Code:

  else if (isfinite(a)) {
    int i, dh;
    double ph = 0;
    dh = d = 1;
    mode = 1; // Exp-Sinh
    for (i = 1; i <= 16; ++i) {
      double r = 1 << i;
      double fl = f(a + d/r)/r, fr = f(a + d*r)*r;
      double h = r*(fl - fr);
      if (dh == d && sign(ph) != sign(h) && fabs(fl) > fabs(fr))
        dh = 1 << (i-1);
      ph = h;
      printf("%g diff=%g (%g - %g)\n", r, h, r*fl, r*fr);
    }
    if (dh != d)
      printf("optimal d=%d\n", dh);

Note that this scales the left fl and right fr by r just to make the value's magnitudes comparable to display, otherwise this is unnecessary. This is not the final version, since this can be further improved.

- Rob