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