Over the time Shepard method, introduced in 1968 in [21], has been improved in order to get better reproduction qualities, higher accuracy and lower computational cost (see, e.g.,[2]-[9], [22], [23]).

Let $f$ be a real-valued function defined on $X\subset {ℝ}^2,$ and $(x_i,y_i)\in X,i=1,\mathrm{\dots },N$ some distinct points. The bivariate Shepard operator is defined by

where

with the parameter $\mu >0$ and $r_i(x,y)$ denoting the distances between a given point $(x,y)\in X$ and the points $(x_i,y_i),i=1,\mathrm{\dots },N$.

In [11], Franke and Nielson introduced a method for improving the accuracy in reproducing a surface with the bivariate Shepard approximation. This method has been further improved in [10], [19], [20], and it is given by:

with

where $R_w$ is a radius of influence about the node $(x_i,y_i)$ and it is varying with $i.$ $R_w$ is taken as the distance from node $i$ to the $j$th closest node to $(x_i,y_i)$ for $j>N_w$ ($N_w$ is a fixed value) and $j$ as small as possible within the constraint that the $j$th closest node is significantly more distant than the $(j-1)$st closest node (see, e.g. [20]). As it is mentioned in [14], this modified Shepard method is one of the most powerful software tools for the multivariate approximation of large scattered data sets.

A.V. Masjukov and V.V. Masjukov introduced in [15] an iterative modification for the Shepard operator that requires no artificial parameter, such as a radius of influence or number of nodes. So, they defined the iterative Shepard operator as

where $w$ is the weight function, continuously differentiable, with the properties that

and $u_j^{(k)}$ denotes the interpolation residuals at the $k$th step, with $u_j^{(0)}\equiv u_j$.

Let $f$ be a real-valued function defined on $X\subset {ℝ}^2.$ We denote by $𝐱$ the point $(x,y)\in X$ and we assume that ${𝐱}_{𝐢}=(x_i,y_i)\in X,$ $i=1,\mathrm{\dots },N^{\prime }$, are some given interpolation nodes.

The radial basis functions (RBF) are some modern and very efficient tools for interpolating scattered data, thus they are intensively used (see, e.g., [1], [12] – [14], [18]). In the sequel we use two radial basis functions that are positive definite, the inverse quadratic RBF and the inverse multiquadric RBF.

Consider the two radial basis functions as

with $ϵ$ being a shape parameter and $r_j(x,y)=\sqrt{{(x-x_j)}^2+{(y-y_j)}^2}.$

For $\beta =-1$, ${\varphi }_{𝐢}^{-\mathrm{𝟏}}$ is the inverse quadratic RBF and for $\beta =-1/2$, ${\varphi }_{𝐢}^{-\mathrm{𝟏}/\mathrm{𝟐}}$ is the inverse multiquadric RBF.

The coefficients ${\alpha }_j,a,b,c$ are obtained as solutions of systems of the form

with $r_{ij}=\sqrt{{(x_i-x_j)}^2+{(y_i-y_j)}^2}$ and $f_i=f({𝐱}_i)$.

Shortly, this system can be written as

considering the following notations:

• •

  $A\in {ℳ}_{N^{\prime }\times N^{\prime }}(ℝ)$, with the element on the entry $(i,j)$ being
  $a_{ij}={\left[1+{\left(ϵr_{ij}\right)}^2\right]}^{\beta }$, where $r_{ij}=\sqrt{{(x_i-x_j)}^2+{(y_i-y_j)}^2}$,
  $i,j=1,\mathrm{\dots },N^{\prime }$ and $\beta \in \{-1,-1/2\}$;

• •

  $O_3$ is the zero square matrix of order 3;

• •

  $𝐮={(a,b,c)}^T$, $𝜶={({\alpha }_1,\mathrm{\dots },{\alpha }_{N^{\prime }})}^T$, $\mathrm{𝟎}={(0,\mathrm{\hspace{0.33em}0},\mathrm{\hspace{0.33em}0})}^T$ ;

• •

  $𝐟={(f_1,\mathrm{\dots },f_{N^{\prime }})}^T$, with $f_i=f({𝐱}_i)$.

First, consider the classical Shepard operator given in (1.1).

Furthermore, we consider the improved form of the Shepard operator, given in (1.3).

Finally, we follow the idea proposed in [15], which consists of using an iterative procedure that requires no artificial parameters.

The functions ${\varphi }_i^{\beta }$ are given in (2.1). We follow ideas from [15] for the parameters’ choice. As an example, the sequence $\{{\tau }_k\}$ of scale factors is defined as

The setup parameter ${\tau }_k$ can be chosen such that it decreases from an initial value ${\tau }_0$, which is given for instance as

to the final value ${\tau }_K$ such that

The behaviour of $u_{\varphi }^{\beta }$ does not change very much for $\gamma$ between $0.6$ and $0.95$, as shown in [15]. One can also choose smaller values for $\gamma$ if the nodes are sparse and a decreased computational time is desired.

Finally, the weight function $w$ is given by

with

We apply the three operators on two sets of points. For the first way, we consider a set of $N$ initial interpolation nodes ${𝐱}_{𝐢},i=1,\mathrm{\dots },N,$ and for the second way, we consider a smaller set of $k\in {\mathbb{N}}^{\ast }$ knot points ${\widehat{𝐱}}_{𝐣},$  $j=1,\mathrm{\dots },k$ that will be representative for the original set. This set is obtained following the next steps (see, e.g., [16] and [17]):

We consider the following test functions (see, e.g., [10], [19], [20]):

Tables 1 - 3 contain the maximum errors for approximating the functions (3.1) by the classical, the modified and the iterative Shepard operators given, respectively, by (1.1), (1.3) and (1.5), and the errors of approximating by the operators introduced in (2.2), (2.3) and (2.4). We construct the operators for both radial basis functions - the inverse quadratic and the inverse multiquadric. For each function we consider a set of $N=100$ random points in $[0,1]\times [0,1]$, a subset of $k=25$ representative knots, $\mu =3$, $N_w=19$, $K=20$, ${\tau }_0=3$ and $\gamma =0.66,\mathrm{\hspace{0.33em}0.84},\mathrm{\hspace{0.33em}0.91}$.

In Figures 1 - 4 we plot the graphs of $f_1,f_2,f_3$ and of the corresponding Shepard operators $S_{\mu }^{\beta }f$, $S_W^{\beta }f$ and $u_{{\varphi }^{\beta }}$, combined with the inverse quadratic ($\beta =-1$) and the inverse multiquadric ($\beta =-1/2$) radial basis functions. We consider the sets of the $k=25$ representative knot points.

We remark that $S_W^{\beta }f$ and $u_{{\varphi }^{\beta }}$ have better approximation properties than the classical Shepard operator $S_{\mu }^{\beta }f$, the results for $u_{{\varphi }^{\beta }}$ depending on the values of $\gamma$. Also, we notice better approximation errors for the lower number of knots obtained using the Algorithm 2.4.