emmy.polynomial.interpolate — org.mentat/emmy 0.32.0

This namespace contains a discussion of polynomial interpolation, and different methods for fitting a polynomial of degree N-1 to N points and evaluating that polynomial at some different x.

lagrange clj/s

(lagrange points x)

Generates a lagrange interpolating polynomial that fits every point in the supplied sequence points (of form [x (f x)]) and returns the value of the polynomial evaluated at x.

The Lagrange polynomial has this form:

g(x) =  (f(a) * [(x-b)(x-c)...] / [(a-b)(a-c)...])
      + (f(b) * [(x-a)(x-c)...] / [(b-a)(b-c)...])
      + ...

for points [a f(a)], [b f(b)], [c f(c)] etc.

This particular method of interpolating x into the polynomial is inefficient; any new calculation requires fully recomputing. Takes O(n^2) operations in the number of points.

modified-neville clj/s

(modified-neville points x)

Similar to neville (the interface is identical) but slightly more efficient. Internally this builds up its estimates by tracking the delta from the previous estimate.

This non-obvious change lets us swap an addition in for a multiplication, making the algorithm slightly more efficient.

See neville for usage information, and info about the required structure of the arguments.

The structure of the modified-neville algorithm makes it difficult to select a particular column. See neville if you'd like to generate polynomial approximations between successive sequences of points.

References:

• "A comparison of algorithms for polynomial interpolation", A. Macleod
• Press's Numerical Recipes (p103), chapter 3

modified-neville-fold clj/s

(modified-neville-fold x)

Given some point x, returns a fold that accumulates rows of an interpolation tableau providing successively better estimates (at the value x) of a polynomial interpolated to all seen points.

The 2-arity aggregation step takes:

• previous-row: previous row of an interpolation tableau
• a new point of the form [x_new (f x_new)]

and returns the next row of the tableau using the algorithm described in modified-neville.

modified-neville-scan clj/s

(modified-neville-scan x)

Returns a function that consumes an entire sequence xs of points of the form [x_i, f(x_i)] and returns a lazy sequence of successive approximations of x using polynomials fitted to the first point, then the first and second points, etc. using the algorithm described in modified-neville.

Equivalent to ([[modified-neville]] xs x).

modified-neville-sum clj/s

(modified-neville-sum x)

Returns a function that consumes an entire sequence xs of points of the form [x_i, f(x_i)] and returns the best approximation of x using a polynomial fitted to all points in xs using the algorithm described in modified-neville.

Faster than, but equivalent to, (last ([[modified-neville]] xs x))

neville clj/s

(neville points x)

(neville points x column)

Takes:

• a (potentially lazy) sequence of points of the form [x (f x)] and
• a point x to interpolate

and generates a lazy sequence of approximations of P(x). Each entry in the return sequence incorporates one more point from points into the P(x) estimate.

Said another way: the Nth in the returned sequence is the estimate using a polynomial generated from the first N points of the input sequence:

p0 p01 p012 p0123 p01234

This function generates each estimate using Neville's algorithm:

$$P(x) = [(x - x_r) P_l(x) - (x - x_l) P_r(x)] / [x_l - x_r]$$

Column

If you supply an integer for the third column argument, neville will return that column of the interpolation tableau instead of the first row. This will give you a sequence of nth-order polynomial approximations taken between point i and the next n points.

As a reminder, this is the shape of the tableau:

 p0 p01 p012 p0123 p01234
 p1 p12 p123 p1234 .
 p2 p23 p234 .     .
 p3 p34 .    .     .
 p4 .   .    .     .

So supplying a column of 1 gives a sequence of linear approximations between pairs of points; 2 gives quadratic approximations between successive triplets, etc.

References:

• Press's Numerical Recipes (p103), chapter 3
• Wikipedia, Neville's Algorithm

neville-fold clj/s

(neville-fold x)

Given some point x, returns a fold that accumulates rows of an interpolation tableau providing successively better estimates (at the value x) of a polynomial interpolated to all seen points.

The 2-arity aggregation step takes:

• previous-row: previous row of an interpolation tableau

• a new point of the form [x_new (f x_new)]

  and returns the next row of the tableau using the algorithm described in neville.

neville-incremental clj/s

(neville-incremental points x)

Takes a potentially lazy sequence of points and a point x and generates a lazy sequence of approximations of P(x).

entry N in the returned sequence is the estimate using a polynomial generated from the first N points of the input sequence.

neville-recursive clj/s

(neville-recursive points x)

Top-down implementation of Neville's algorithm

Returns the value of P(x), where P is a polynomial fit (using Neville's algorithm) to every point in the supplied sequence points (of form [x (f x)])

The efficiency and results should be identical to [[emmy.numerical.interpolate/lagrange]]. This function represents a step on the journey toward more incremental methods of polynomial interpolation.

References:

• Press's Numerical Recipes (p103), chapter 3
• Wikipedia, Neville's Algorithm

neville-scan clj/s

(neville-scan x)

Returns a function that consumes an entire sequence xs of points of the form [x_i, f(x_i)] and returns a lazy sequence of successive approximations of x using polynomials fitted to the first point, then the first and second points, etc. using the algorithm described in neville.

Equivalent to ([[neville]] xs x).

neville-sum clj/s

(neville-sum x)

Returns a function that consumes an entire sequence xs of points of the form [x_i, f(x_i)] and returns the best approximation of x using a polynomial fitted to all points in xs using the algorithm described in neville.

Faster than, but equivalent to, (last ([[neville]] xs x))

tableau-fn clj/s

(tableau-fn prepare merge points)

Returns a Newton-style approximation tableau, given:

• prepare: a fn that processes each element of the supplied points into the state necessary to calculate future tableau entries.

• merge: a fn of land r the tableau entries:

l -- return
   /
  /
 /
r

the inputs are of the same form returned by prepare. merge should return a new structure of the same form.

• points: the (potentially lazy) sequence of points used to generate the first column of the tableau.

tableau-fold-fn clj/s

(tableau-fold-fn prepare merge present)

Given prepare and merge and present functions, returns a fold capable of aggregating a point of the form [x, f(x)] into an accumulating tableau row (generating the next tableau row).

The 0-arity of the returned function returns an empty row, [].

The 1-arity calls the supplied present on the accumulated tableau row.

The 2-arity scans the supplied merge across all entries in the accumulating row, producing a new row.

More detail on the arguments:

• prepare: a fn that processes each element of the supplied points into the state necessary to calculate future tableau entries.

• merge: a fn of land r the tableau entries:

l -- return / / / r

the inputs are of the same form returned by prepare. merge should return a new structure of the same form.

• present: Transforms a tableau row into an estimate at some value x of the polynomial interpolated to hit all supplied points.