sicmutils.numerical.interpolate.rational — sicmutils 0.14.0

This namespace contains a discussion of rational function interpolation, and different methods for fitting rational functions N points and evaluating them at some value x.

bs-merge clj/s

(bs-merge x)

Implements the recursion rules described in Press's Numerical Recipes, section 3.2 http://phys.uri.edu/nigh/NumRec/bookfpdf/f3-2.pdf to generate x_l, x_r, C and D for a tableau node, given the usual left and left-up tableau entries.

This merge function ALSO includes a 'zero denominator fix used by Bulirsch and Stoer and Henon', in the words of Sussman from rational.scm in the scmutils package.

If the denominator is 0, we pass along C from the up-left node and d from the previous entry in the row. Otherwise, we use the algorithm to calculate.

TODO understand why this works, or where it helps!

bs-prepare clj/s

(bs-prepare [x fx])

Processes an initial point [x (f x)] into the required state:

[x_l, x_r, C, D]

The recursion starts with $C = D = f(x)$.

bulirsch-stoer clj/s

(bulirsch-stoer 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.

P(x) is rational function fit (using the Bulirsch-Stoer algorithm, of similar style to Neville's algorithm described in polynomial.cljc) to every point in the supplied sequence points.

"The Bulirsch-Stoer algorithm produces the so-called diagonal rational function, with the degrees of numerator and denominator equal (if m is even) or with the degree of the denominator larger by one if m is odd." ~ Press, Numerical Recipes, p105

The implementation follows Equation 3.2.3 on on page 105 of Press: http://phys.uri.edu/nigh/NumRec/bookfpdf/f3-2.pdf.

Column

If you supply an integer for the third (optional) column argument, bulirsch-stoer will return that /column/ offset 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 2-point approximations between pairs of points; 2 gives 3-point approximations between successive triplets, etc.

References:

• Stoer & Bulirsch, 'Introduction to Numerical Analysis': https://www.amazon.com/Introduction-Numerical-Analysis-Applied-Mathematics/dp/144193006X
• PDF of the same: http://www.math.uni.wroc.pl/~olech/metnum2/Podreczniki/(eBook)%20Introduction%20to%20Numerical%20Analysis%20-%20J.Stoer,R.Bulirsch.pdf
• Press's Numerical Recipes (p105), Section 3.2 http://phys.uri.edu/nigh/NumRec/bookfpdf/f3-2.pdf

bulirsch-stoer-recursive clj/s

(bulirsch-stoer-recursive points x)

Returns the value of P(x), where P is rational function fit (using the Bulirsch-Stoer algorithm, of similar style to Neville's algorithm described in polynomial.cljc) to every point in the supplied sequence points.

points: is a sequence of pairs of the form [x (f x)]

"The Bulirsch-Stoer algorithm produces the so-called diagonal rational function, with the degrees of numerator and denominator equal (if m is even) or with the degree of the denominator larger by one if m is odd." ~ Press, Numerical Recipes, p105

The implementation follows Equation 3.2.3 on on page 105 of Press: http://phys.uri.edu/nigh/NumRec/bookfpdf/f3-2.pdf.

References:

• Stoer & Bulirsch, 'Introduction to Numerical Analysis': https://www.amazon.com/Introduction-Numerical-Analysis-Applied-Mathematics/dp/144193006X
• PDF of the same: http://www.math.uni.wroc.pl/~olech/metnum2/Podreczniki/(eBook)%20Introduction%20to%20Numerical%20Analysis%20-%20J.Stoer,R.Bulirsch.pdf
• Press's Numerical Recipes (p105), Section 3.2 http://phys.uri.edu/nigh/NumRec/bookfpdf/f3-2.pdf

modified-bulirsch-stoer clj/s

(modified-bulirsch-stoer points x)

Similar to bulirsch-stoer (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 division, making the algorithm slightly more efficient.

See the bulirsch-stoer docstring for usage information, and info about the required structure of the arguments.

References:

• Press's Numerical Recipes (p105), Section 3.2 http://phys.uri.edu/nigh/NumRec/bookfpdf/f3-2.pdf

modified-bulirsch-stoer-fold clj/s

(modified-bulirsch-stoer-fold x)

Returns a function that consumes an entire sequence xs of points, and returns a sequence of successive approximations of x using rational functions fitted to the points in reverse order.

modified-bulirsch-stoer-fold-fn clj/s

(modified-bulirsch-stoer-fold-fn x)

Returns a function that accepts:

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

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

modified-bulirsch-stoer-scan clj/s

(modified-bulirsch-stoer-scan x)

Returns a function that consumes an entire sequence xs of points, and returns a sequence of SEQUENCES of successive rational function approximations of x; one for each of the supplied points.

For a sequence a, b, c... you'll see:

[(modified-bulirsch-stoer [a] x) (modified-bulirsch-stoer [b a] x) (modified-bulirsch-stoer [c b a] x) ...]