Variant of the Secant Method

[Albert Chan]

Quadratic interpolation is simply 2 linear interpolations, followed by linear interpolation.

x3 y3
x2 y2 y23
x1 y1 y13 y123

Similarly, cubic interpolation is 2 quadratic interpolations, followed by linear interpolation.

x4 y4
x3 y3
x2 y2 y234
x1 y1 y134 y1234

Using this recursive nature of polynomial interpolation, I get this neat Lua code

Code:

function interp(x1,y1,x2,y2,x3,y3,x4,y4)    -- estimate y = p(0)
    if x3 then
        y2 = interp(x2,y2,x3,y3,x4,y4)
        y1 = interp(x1,y1,x3,y3,x4,y4)
    end
    local y, p = y2, x2/(x2-x1)
    if p > 0.5 then y, p = y1, x1/(x2-x1) end
    return y - p*(y2-y1)
end

Extending interp(...) is trivial, by adding (x5,y5), ...
The down-side is each bump up in polynomial degree, doubled the work.
For high degree fit, divided-difference style is better, because it can cache previous results.

OTTH, assuming we have convergence, old points are very, very far away.
Fitting those probably may not help much ...

---

For slope of improved secant, we can also consider this an interpolation problem.
Example, for a cubic fit, we generate divided-difference table:

x4 y4
x3 y3 y34
x2 y2 y24 y234
x1 y1 y14 y123 y1234

p(x) = y4 + (x-x4)*(y34 + (x-x3)*(y234 + (x-x2)*y1234)) = y4 + (x-x4)*M

p'(x) = M + (x-x4)*M'       → p'(x4) = M

Basing all values relative to (x4, y4). Example X3 = x3-x4, Y3 = y3-y4 ...
M is just interpolation of individual slopes.

X3      Y3/X3
X2      Y2/X2
X1      Y1/X1      M

Again, extending this is simply adding more slopes to fit

Code:

function secant(x1,y1,x2,y2,x3,y3,x4,y4)    -- estimate root of p(x)
    if not x3 then return interp(y1,x1,y2,x2) end
    if x4 then x3,y3,x4,y4 = x4,y4,x3-x4,(y3-y4)/(x3-x4) end
    x1, y1, x2, y2 = x1-x3, y1-y3, x2-x3, y2-y3
    return x3 - y3 / interp(x1,y1/x1, x2,y2/x2, x4,y4)    
end

Example:

lua> function f(x) return sin(x)/x - 1/6.5 end
lua> x1, x2 = 2.711, 2.712
lua> y1, y2 = f(x1), f(x2)
lua> x3 = secant(x1,y1, x2,y2)
lua> y3 = f(x3)

lua> interp(y1,x1, y2,x2, y3,x3) -- inverse interpolate, x = g(0)
2.711312910356982
lua> secant(x1,y1, x2,y2, x3,y3) -- improved secant, x = x3 - y3/slope
2.7113129103569826