D. Shepard introduced in 1968 in [13] a very powerful method for approximating a given function $f$ on a set of scattered data, method that nowadays is named after him. The procedure has an easy implementation and it is expressed as a combination between some basis functions and the values of the function $f$ on a given set of interpolation nodes. However, two of its major drawbacks are the high computational cost and the low degree of exactness. Several authors have studied them and proposed different solutions to overcome them, such as modifying the basis functions or combining the Shepard operator with other interpolation operators for an increased degree of exactness (see, e.g., [1]; [2]; [3]; [4]; [5]; [6]; [7]; [8]; [9]; [10]).

In the univariate case, when $f$ is a real-valued function defined on a subset $X$ of $ℝ$, for a given set of $K$ interpolation nodes, $x_i\in X,i=1,\mathrm{\dots },K$, the Shepard operator is defined as

with the basis functions $A_{i,\mu }$ given by

$x\in X$ and $\mu >0$ an arbitrary parameter.

This paper focuses on introducing a new univariate Shepard operator, combined with polynomials constructed based on the least squares approach. In Section 2, after we construct these polynomials, we study several properties of them and of the combined Shepard operators (interpolation property, degree of exactness, linearity). Finally, we study the errors, based on Peano’s Theorem. Section 3 is dedicated to numerical examples that show the benefits of these Shepard operators.

R. J. Renka introduced in 1988 in [11] an algorithm for improving the bivariate Shepard operator, considering a quadratic polynomial that interpolates the function $f$ on a set of given nodes and also approximates the data in a weighted least squares way. Later on, in 1999 in [12] he improved this method by replacing the quadratic polynomial with a cubic one. In 2010 in [14], W. I. Thacker et al. emphasized the main disadvantages of these two methods and suggested the combination of the Shepard operator with a linear polynomial that still fits the data in a weighted least squares sense.

Using some ideas for the bivariate case presented in the above mentioned papers, we are going to consider an improvement for the classical Shepard operator in the univariate case, by combining it with polynomials of degree $n,n\in \mathbb{N}$, constructed following the weighted least squares approximation technique.

Consider $X\subset ℝ$, $f:X\to ℝ$ and $K$ given real nodes, denoted by $x_j$, $j=1,\mathrm{\dots },K$. The values of the function $f$ on the given nodes are known and denoted by $f_j=f(x_j),j=1,\mathrm{\dots },K$.

Under these assumptions, for a point $x\in X$, let us define the $n$th degree polynomial function $C_j^n[f]$, $j=1,\mathrm{\dots },K$, $n\in \mathbb{N}$, as

where the coefficients $a_{j,k}$ are found such that they minimize the sum of the weighted squared residuals

where

for $i,j=1,\mathrm{\dots },K$ and $\mu >0$.

In order to find the coefficients $a_{j,k}$ (i.e, obtain the minimum of expression (4)), we follow the weighted least squares reasoning, take the partial derivatives of $E_j$ with respect to each unknown, set them to zero and solve the resulting system:

Further, for every $j=1,\mathrm{\dots },K$ one obtains

Let us make the notation

Then, the system of normal equations that has to be solved in order to find the coefficients $a_{j,k},k=1,\mathrm{\dots },n,$ has the form

for each $j=1,\mathrm{\dots },K$.

For every $j=1,\mathrm{\dots },K,$ we can write the normal equations that appear above in a matricial form as

where $M_j$ is a $n\times n$ matrix having on the entry $(r,s)$ the element $\sum _{\begin{array}{c}i=1\\ i\ne j\end{array}}^K{\lambda }_{i,j}\cdot {(x_i-x_j)}^{r+s}$, $b_j$ is a vector of $n$ elements with $\sum _{\begin{array}{c}i=1\\ i\ne j\end{array}}^K{\lambda }_{i,j}\cdot {(x_i-x_j)}^k\cdot (f_i-f_j)$ on the $k$th entry and $a_j={(a_{j,1},a_{j,2},\mathrm{\dots },a_{j,n})}^T$ is the vector of unknowns.

We now introduce the interpolation formula for the univariate Shepard combined with a polynomial, that is given by

with $R_n[f]$ denoting the remainder.

Considering the space $H^m[a,b]$, $m\in \mathbb{N}\setminus \{0\}$ of functions $f\in C^{m-1}[a,b]$ (continuously differentiable up to order $m-1$, inclusively) with $f^{(m-1)}$ absolutely continuous on $[a,b]$, we obtain the following result for the remainder of the formula:

For the numerical experiments we consider four well-known real-valued test functions (see, e.g., [9]):

For each function $f_i,i=1,\mathrm{\dots },4$, we compare the test results obtained by considering the linear, quadratic and cubic interpolants $S{P}_j[f_i],j=1,2,3,$ with the ones obtained for some other combined Shepard operators, well-known in the literature. We study the maximum approximation errors for the Shepard operators of Lagrange, Taylor and Bernoulli type, for all of them considering the first, second and third order. They are denoted by $S{L}_k,S{T}_k$ and $S{B}_k$, $k=1,2,3$, respectively. We recall their definitions in the sequel.