Multivariate polynomial interpolation in Newton forms

Abstract

Techniques of univariate Newton interpolating polynomials are extended to multivariate data points by different generalizations and practical algorithms. The Newton basis format, with divided-difference algorithm for coefficients, generalizes in a straightforward way when interpolating at nodes on a grid within certain regions. While this idea, known as the tensor product case, has been around for over a century, this exposition uses modern multi-index notation to provide generality, concise algorithms, and rigor, appropriate as a topic in introductory numerical analysis. Node configurations where the tensor product case applies are called lower sets, and they include n-dimensional rectangles (boxes) and triangles (corners), the latter having unique interpolation over polynomials of fixed degree. Arbitrary distinct nodes do not ensure unique interpolating polynomials but, when possible, a different basis generalization, which we call Newton-Sauer, results in a nice triangular linear system for the coefficients. For the right number of distinct nodes, an algorithm based on Gaussian elimination produces either this Newton-Sauer basis or a polynomial that is zero on all nodes, showing that unique interpolation is impossible. The algorithm may also be used on any distinct nodes to produce a polynomial subspace of minimal degree where unique interpolation is possible. A variation of this algorithm using matrix block operations produces another basis generalization that we call Newton-Olver, which exactly agrees with the classic one for a corner of nodes and computes an interpolating polynomial using formulas similar to divided differences.