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

[Albert Chan]

(04-24-2021 02:31 AM)robve Wrote:  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.

Max. of 12 points ... impressive
For illustration purpose, I tried 1 from Appendix C, where suggested d does most good.
(there were 1 that does better, but integrand is oscillatory, so it might be just luck)

lua> f = function(x) local y=sqrt(x*x-1); return exp(-y*11/1000)*x/y end
lua> Q.exp_sinh_opt_d(f, 1) -- default eps=1e-9, d=1.
512

Code:

N   x         f(1+x)                   f(1+x)*x                sign(top-bot)
1   1/2       1.3252418334225524       0.6626209167112762
2   2         1.028168248457613        2.056336496915226            -    
3   1/131072  255.99046499382536       0.0019530522536760357
4   131072    0                        0                            +
5   1/512     16.012410275011778       0.03127423881838238
6   512       0.003542270480289993     1.8136424859084763           -
7   1/8912    63.9948588286239         0.007180751663894064
8   8192      7.245672531983614e-040   5.9356549382009764e-036      +
9   1/2048    32.000714078579314       0.015625348671181306
10  2048      1.6271885884174952e-010  3.33248222907903e-007        +
11  1/1024    22.632978376948273       0.022102517946238548
12  1024      1.2686220440383643e-005  0.01299068973095285          +

I also made a version based on peak, and extrapolate for x-intercept.

I've noticed after the peak, the curve mostly go straight down.
(code use tangent line of last point, of the fitted quadratic)

Of course, with peak, there are really 2 valleys.
Code assumed optimal d is to the right of peak.

Code:

function Q.peak(f, a, d)
    a = a or 0
    d = d or 1
    local L, R, k, P = f(a+d), f(a+d*2)*2, 2, huge-huge
    if R/L < 1.2 then return Q.peak(f, a, d/100) end
    while L/R < 1.2 do      -- stop if a bit passes peak
        P, L, k = L, R, k+k
        R = f(a+d*k)*k
    end
    d = d*k                 -- (0,R),(1,L),(2,P), 0 as reference
    local d0, d1 = L-R, P-L -- finite difference, fit quadratic p(x)
    P = d0/(d1-d0) - 1/2    -- interpolate for peak, p'(-x)=0, for x
    R = R/(1.5*d0-0.5*d1)   -- extrapolate x-intercept, -x = R/p'(0)
    return d*2^P, d*2^R     -- peak, x-intercept
end

lua> Q.peak(f, 1)
89.13375894266683       399.8238469921945

Note that code evaluated an extra point passes the peak (when points are too close)
This is to get a more reasonable slope for the tangent line.
Without the extra point, we would get peak, x-intercept ≈ 88.25, 9106.

Code:

N   x         f(1+x)                   f(1+x)*x
1   1         1.1329087918426253       1.1329087918426253
2   2         1.028168248457613        2.056336496915226
3   4         0.9670764022530607       3.8683056090122427
4   8         0.9119448784679965       7.295559027743972
5   16        0.8311515864637837       13.29842538342054
6   32        0.6960220436099958       22.272705395519864
7   64        0.4892914169125088       31.314650682400565
8   128       0.2419734386421159       30.972600146190835
9   256       0.05919187288015089      15.153119457318628

Let's see effect of suggested d, vs default d=1

lua> function quad(f,...) local s,e = Q.quad(Q.count(f),...); return s,e,Q.n end

lua> quad(f, 1, huge) -- default d=1
90.90909088146049       3.999609164190213e-009       273
lua> quad(f, 1, huge, 512) -- sign change based d
90.90909089915341       6.360683358493127e-011       143
lua> quad(f, 1, huge, 400) -- x-intercept based d
90.90909089845148       2.0323842611238905e-011     143
lua> quad(f, 1, huge, 89) -- peak based d
90.90909089656193       4.928270670014641e-011       141

The peak may be used for d.
In that case, exp-sinh transformed plot, when folded, turn into half a bell-shaped curve.