Variant of the Secant Method

[Albert Chan]

(12-11-2020 04:34 PM)Albert Chan Wrote:  Quadratic interpolation is simply 2 linear interpolations, followed by linear interpolation.

x3 y3
x2 y2 y23
x1 y1 y13 y123

It may seems we need 3 secant root steps for 3-points fit, in a loop, we only need 2

x2 y2
x3 y3 y23
x1 y1 y12 y123

Next iteration y23 → y12, can be reused.

Secant setup code

Code:

function newx(a,fa, b,fb) return b - fb/(fb-fa)*(b-a) end

function secant_init(f,a,b,eps,verbal)
    eps = eps or 1e-9
    local f2 = verbal and function(x) print(x); return f(x) end or f
    local fa = f2(a); if fa==0 then return nil, a end
    local fb = f2(b); if fb==0 then return nil, b end
    if signbit(eps) and abs(fa)<abs(fb) then
        a,fa, b,fb = b,fb, a,fa
        if verbal then print('-- points swapped') end
    end
    return f2, a,fa, b,fb, abs(eps)
end

Code:

function S.secant2(f,a,b,eps,verbal,c)  -- 3 pt fit
    local fa,fb,fc,ab,bc
    f, a,fa, b,fb, eps = secant_init(f,a,b,eps,verbal)
    if not f then return a,a end
    ab = newx(a,fa, b,fb)
    c = c or ab
    fc = f(c); if fc==0 then return c,c end
    while abs(b-c) > eps do
        bc, ab = ab, newx(b,fb, c,fc)
        a,fa, b,fb, c = b,fb, c,fc, newx(bc,fa, ab,fc)
        fc = f(c); if fc==0 then return c,c end
    end
    return finite(newx(b,fb, c,fc)) or (b+c)/2
end

lua> f = fn'x: x^3 - 8'
lua> x = S.secant2(f, 5, 4, 1e-9, true)
5
4
3.081967213114754
2.3945608933295457
2.0980163342039635
2.0088865345029276
2.0001129292461193
2.0000000388749792
2.0000000000000164
2

Result matched post #8

---

For completeness, here is S.secant3 (OP "improved" secant)

Code:

function S.secant3(f,a,b,eps,verbal,c)  -- slope extrapolated
    local fa,fb,fc
    f, a,fa, b,fb, eps = secant_init(f,a,b,eps,verbal)
    if not f then return a,a end
    c = c or newx(a,fa, b,fb)
    fc = f(c); if fc==0 then return c,c end
    while abs(b-c) > eps do
        local ca,cb = c-a,c-b
        local slope = newx((fc-fa)/ca,ca, (fc-fb)/cb,cb)
        a,fa, b,fb, c = b,fb, c,fc, c-fc/slope
        fc = f(c); if fc==0 then return c,c end
    end
    return finite(newx(b,fb, c,fc)) or (b+c)/2
end

lua> x = S.secant3(f, 5, 4, 1e-9, true)
5
4
3.081967213114754
2.2862188297178117
2.010344209437878
1.99979593345267
2.0000000722313933
2.000000000000015
2