[SciPy-user] Polynomial interpolation

[Anne Archibald]

Hi,

It appears that scipy does not have a facility for using the Lagrange
polynomial to interpolate data. (Or did I miss it?) I am well aware of
the numerical difficulties this poses, but it (and its generalization,
the Hermite polynomial) does prove useful on occasion. I have written
prototype implementations of two algorithms for evaluating this
polynomial, and I'd like comments before submitting them for inclusion
in scipy.

One implementation is based on Krogh 1970, "Efficient Algorithms for
Polynomial Interpolation and Numerical Differentiation"; it allows the
construction of Hermite polynomials (where some derivatives at each
point may also be specified) and the evaluation of arbitrary
derivatives. It is based on a Neville-like algorithm, so that it does
O(n^2) work when constructed but only O(n) per point evaluated, or
O(n^2) per point for which all derivatives must be evaluated. (n here
is the degree of the polynomial.)

The other implementation is based on Berrut and Trefethen 2004,
"Barycentric Lagrange Interpolation". This implementation uses a
rewriting of the Lagrange polynomial as a rational function, and is
efficient and numerically stable. It also allows efficient updating of
the y values at which interpolation occurs, as well as addition of new
x values. It does not support evaluation of derivatives or
construction of Hermite polynomials.

Finally, I have written a "PiecewisePolynomial" class for constructing
splines in which each piece may have arbitrary degree, and for which
the function values and some derivatives are specified at each knot.
The intent is that this be used to represent solutions obtained from
ODE solvers, using one polynomial for each solver step, with the same
order as the solver's internal polynomial solution, and with (some)
derivatives matching at the ends. Such a Trajectory object would
contain additional information about the solution (for example
stiffness or error estimates) beyond what is in PiecewisePolynomial.
(I tried implementing Trajectory using splines, which are evaluated in
compiled code, but their maximum degree is 5 while the solvers will go
up to degree 12.)

They all need work, in particular, efficiency would be improved by
making the y values vectors, error checking needs to be more robust,
and documentation is not in reST form. Ultimately, too, the evaluation
functions at least should be written in a compiled language (cython?).
But I thought I'd solicit comments on the code first - is the
object-oriented design cumbersome? Is including the algorithm in the
name confusing to users? Is the calling convention for Hermite
polynomials too confusing or error-prone?

Thanks,
Anne
-------------- next part --------------
A non-text attachment was scrubbed...
Name: polyint.py
Type: text/x-python
Size: 9798 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20080419/de485866/attachment.py>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test_polyint.py
Type: text/x-python
Size: 2859 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20080419/de485866/attachment-0001.py>