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

[Albert Chan]

(04-29-2021 11:38 PM)Albert Chan Wrote:  Without knowing the shape of curve, peak for d is the safe choice.

This is work in progress, but if we got peak and both x-intercepts, we have an idea of curve shape.

Code:

function Q.peak2(f, a, d)       -- find peak of f(a+x)*x, guess=d
    a = a or 0
    d = d or 1
    local k, L, C, R = 2, f(a+d), f(a+d*2)
    while L/C > 2.353 do d=d/2; R=C; C=L; L=f(a+d) end
    if R then d=d*4 else d=d*2; R=C; C=L end
    while R/C > 0.425 do        -- stop if f(a+x)*x passes peak
        d = d+d
        L, C, R = C, R, f(a+d)
    end
    L, C = L/4, C/2             -- 3 points: (0,R),(1,C),(2,L)
    local d0, d1 = C-R, L-C     -- finite difference, fit quadratic p(x)
    local P = d0/(d1-d0) - 0.5  -- interpolate for peak, p'(x)=0, for -x
    R = R/(1.5*d0-0.5*d1)       -- extrapolate x-intercept, -x = R/p'(0)    
    local shift = (L/C > 0.85)  -- extra point for left-side x-intercept   
    if shift then C=L; L=f(a+d/8)/8; d0=d1; d1=L-C end
    L = L/(1.5*d1-0.5*d0) - (shift and 3 or 2)
    return d*2^P, L-P, R-P  -- both x-intercepts, log2 scale, rel. to peak
end

lua> f = function(x) return log(1+9*(x*x))/(1+16*(x*x)) end -- previous post example
lua> Q.peak2(f1)
0.9456442962396107       -2.580895653494933       3.5331822306621703

Peak for f(x)*x, at x ≈ 0.9456
The other 2 numbers are x-intercepts, in log2 scale, relative to peak:
Left-side x-intercept ≈ 0.9456 * 2^-2.581 ≈ 0.158
Right-side x-intercept ≈ 0.9456 * 2^3.533 ≈ 10.9

We should pick the steeper side (x-intercept with minimum log2 absolute value)
If there is no steeper side, curve is relatively symmetrical, just use peak for d.

Recommended d = left-side = 0.158

---

Both of these has peak at 1.0 ... which side to use for d ?

lua> f1 = function(x) return exp(-x) end
lua> f2 = function(y) return f1(1/y)/y^2 end
lua> Q.peak2(f1)
0.932573172731174       -2.983219913988427         2.6203047422648424
lua> Q.peak2(f2)
1.072301916075236       -2.6203047422648424       2.9832199139884263

Noted the symmetry of f1, f2 intercepts. This is expected, f2 ≡ f1, with x=1/y, dx=-dy/y²
Again, picking the steeper side, symmetry is carried to optimal d:

f1 optimal d ≈ 2^2.62 ≈ 6.15
f2 optimal d ≈ 2^-2.62 ≈ 0.163