On the approximation order of triangular Shepard interpolation

Shepard's method is a well-known technique for interpolating large sets of scattered data. The classical Shepard operator reconstructs an unknown function as a normalized blend of the function values at the scattered points, using the inverse distances to the scattered points as weight functions. Based on the general idea of defining interpolants by convex combinations, Little suggested extending the bivariate Shepard operator in two ways. On the one hand, he considers a triangulation of the scattered points and substitutes function values with linear polynomials which locally interpolate the given data at the vertices of each triangle. On the other hand, he modifies the classical point-based weight functions and defines instead a normalized blend of the locally interpolating polynomials with triangle-based weight functions which depend on the product of inverse distances to the three vertices of the corresponding triangle. The resulting triangular Shepard operator interpolates all data required for its definition and reproduces polynomials up to degree 1, whereas the classical Shepard operator reproduces only constants. In this paper, we show that this interpolation operator consequentially has quadratic approximation order, which is confirmed by our numerical results.