Rational Function Interpolation

This documentation was generated with the assistance of AI. Please report any inaccuracies.

Some functions are poorly approximated by polynomials but well approximated by rational functions — quotients of polynomials — because rational functions can model poles that would otherwise force a polynomial approximation to oscillate wildly. This package implements two rational interpolators covering two different trade-offs, described in Press et al., 2007, Section 3.4 (page 124): the classical Bulirsch-Stoer diagonal rational interpolant (RationalInterpolator), and the pole-free barycentric variant (BarycentricRationalInterpolator, Section 3.4.1, page 127).

Diagonal rational interpolation (Bulirsch-Stoer)

A rational function through points is:

with unknown coefficients. Bulirsch and Stoer found a Neville-like recurrence, directly analogous to polynomial interpolation's, that builds the diagonal rational interpolant (numerator and denominator of equal degree, or denominator one degree higher) column by column:

started from and an initial for the (empty) case. Exactly as with Neville’s algorithm, this is rewritten in terms of small correction differences and , which is what RationalInterpolator actually implements — the same tableau-building shape as Polynomial Interpolation, but with a different correction formula and, correspondingly, an error estimate that comes for free the same way.

Because this method allows poles anywhere on the real axis, an exact pole at the requested is detected and reported as a failure, rather than silently returning an infinite or nonsensical value. Plausible values for the number of local points are 4 to 7.

Barycentric rational interpolation

Constructing a global diagonal rational interpolant across the entire table (letting its order grow to ) risks introducing poles inside the interpolation interval, even when the underlying function has none. Barycentric rational interpolation avoids this entirely, at a chosen, fixed approximation order , by writing every rational interpolant in barycentric form:

with weights chosen so the approximation order is (error for point spacing ):

This formulation has no poles anywhere on the real axis, competes very favorably with cubic splines (often smaller error, and infinitely smooth rather than only twice continuously differentiable), and — unusually for the methods in this package — is worth experimenting with at higher orders ( ) rather than being kept deliberately low.

Algorithm

flowchart TD subgraph Diagonal["RationalInterpolator (diagonal, Bulirsch-Stoer)"] A1["Initialize c[i] = d[i] = y[i] for the m nearest points"] --> B1["For each tableau column,\nupdate c[i], d[i] via the rational-interpolation recurrence"] B1 --> C1{"pole detected\n(zero-over-zero)?"} C1 -->|yes| D1["Report interpolation failure"] C1 -->|no| E1["Pick c or d, add correction to y"] E1 --> F1["Return y; last correction = error estimate dy"] end subgraph Barycentric["BarycentricRationalInterpolator"] A2["Precompute weights w_k for all N nodes (once)"] --> B2["interpolate(x): sum w_i*y_i/(x-x_i)\nover sum w_i/(x-x_i), over ALL N points"] end

Library implementation

Both extend BaseInterpolator, but differ in scope: RationalInterpolator is a local method like Polynomial Interpolation, re-centering on a different subrange of points for each call; BarycentricRationalInterpolator is global — it precomputes weights once over all points and overrides interpolate(x) to skip the locate/hunt search entirely, since every call touches the whole table anyway.

Class Scope Notes

RationalInterpolator

Local, m nearest points (constructor parameter)

getDy() returns the error estimate; throws if the interpolant has a pole at the requested x.

BarycentricRationalInterpolator

Global, fixed order d (constructor parameter, must be less than the number of points)

No error estimate; cost is O(N*d) once to build the weights, O(N) per subsequent interpolate(x) call (more expensive per call than a spline, which is O(d)).
Start with a small d for BarycentricRationalInterpolator and increase it only if accuracy demands it — the library follows Press et al., 2007's own advice here.

Example

// local diagonal rational interpolation, 5 nearest points
final var rational = new RationalInterpolator(x, y, 5);
final var value1 = rational.interpolate(target);

// global barycentric rational interpolation, order 4
final var barycentric = new BarycentricRationalInterpolator(x, y, 4);
final var value2 = barycentric.interpolate(target);