---
title: 'Multinode Shepard Method: Theory & Applications'
url: https://www.emergentmind.com/topics/multinode-shepard-method
type: topic
---

# Multinode Shepard Method: Theory & Applications

Searching arXiv for recent and foundational papers on multinode Shepard methods and closely related variants.
The multinode Shepard method is a family of Shepard-type approximation and interpolation schemes in which the elementary contribution associated with a single node is replaced by a contribution associated with a local cluster of nodes, typically through a local polynomial interpolant blended by normalized inverse-distance-like weights. In contrast with classical Shepard interpolation, which reproduces only constants, multinode formulations are designed to inherit higher polynomial reproduction, higher approximation order, and improved derivative behavior while retaining partition-of-unity structure and rational weighting. Recent work has specialized the method to regular rectangular grids for digital elevation model reconstruction [2508.07764], to meshfree collocation for two-dimensional Black–Scholes equations [2508.08023], to quasi-histopolation with jump discontinuities [2508.07070], and to univariate approximation of Caputo fractional derivatives and Bagley–Torvik equations [2508.08067]. Related developments include triangular Shepard–Bernoulli constructions for scattered bivariate data [1406.5962], scaled Shepard approximation on quasi-uniform sets [1702.04764], grid-based tensor-product collocation on Cartesian grids [2606.14503], and data-dependent adaptive weighting for reduced smearing near discontinuities [2606.20332].

## 1. Definition and core construction

The defining idea of the multinode Shepard method is to replace pointwise data values in a Shepard average by local approximants attached to subsets of nodes. In the general multivariate setting, one starts from a node set
\[
\mathcal{X}=\{\boldsymbol{x}_1,\dots,\boldsymbol{x}_n\}\subset\mathbb{R}^d,
\]
chooses a family of subsets or clusters
\[
\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},
\]
and requires each cluster to be unisolvent for a polynomial space of prescribed degree. For total degree at most \(r\), the cluster size is
\[
t=t(d,r)=\binom{d+r}{r}.
\]
Each cluster supports a local polynomial interpolant \(p_j[f]\), and the global approximant is formed as a normalized weighted blend [2508.07764]:
\[
\mathcal{MS}_u[f](\boldsymbol{x})=\sum_{j=1}^q W_{u,j}(\boldsymbol{x})\,p_j[f](\boldsymbol{x}),
\]
with multinode weights
\[
W_{u,j}(\boldsymbol{x})
= \frac{\displaystyle\prod_{\iota=1}^{t} \|\boldsymbol{x}-\boldsymbol{x}_{j_\iota}\|_2^{-u}}
{\displaystyle\sum_{l=1}^{q}\prod_{\lambda=1}^{t} \|\boldsymbol{x}-\boldsymbol{x}_{l_\lambda}\|_2^{-u}},
\qquad u>0.
\]

These weights satisfy nonnegativity and partition of unity, and vanish at nodes outside the associated cluster [2508.07764]. This cluster-level weighting is the source of the term “multinode”: one basis function depends on a group of nodes rather than a single one. A closely related formulation appears in the collocation setting, where the local polynomial \(p_j[f]\) is written in barycentric Lagrange form and the blended operator is recast into a purely nodal cardinal basis [2508.08023]:
\[
\mathcal{M}_\mu[f](\mathbf{x}) = \sum_{j=1}^m B_{\mu,j}(\mathbf{x})\,p_j[f](\mathbf{x})
= \sum_{i=1}^n W_{\mu,i}(\mathbf{x})\,f_i.
\]
In that formulation, the multinode cardinal basis functions satisfy
\[
W_{\mu,i}(\mathbf{x}_k)=\delta_{ik},\qquad \sum_{i=1}^n W_{\mu,i}(\mathbf{x})=1,
\]
so interpolation at the original nodes is recovered [2508.08023].

The principal distinction from classical Shepard interpolation is therefore structural rather than merely parametric. Classical Shepard assigns one normalized inverse-distance weight to each data value. Multinode Shepard assigns one normalized inverse-distance-product weight to each local stencil or patch, and each patch contributes a polynomial rather than a constant. This immediately raises the degree of exactness from zero to the degree of the local interpolants, provided the covering and unisolvency assumptions are satisfied [2508.07764].

## 2. Relation to classical Shepard and modified Shepard families

Classical Shepard interpolation in \(\mathbb{R}^d\) has the form
\[
S_u[f](\boldsymbol{x})
=\frac{\sum_{i=1}^n \|\boldsymbol{x}-\boldsymbol{x}_i\|^{-u} f(\boldsymbol{x}_i)}
{\sum_{i=1}^n \|\boldsymbol{x}-\boldsymbol{x}_i\|^{-u}},
\qquad u>0,
\]
or, more generally, a normalized weighted average with a radial kernel [2508.07764, 1702.04764]. Its most basic algebraic property is constant reproduction. In the scattered-data literature, this low degree of exactness and the appearance of flat spots for certain parameter regimes motivated combined, modified, and local Shepard variants, in which each weight is paired with a local polynomial, spline, or generalized Taylor approximant rather than a single function value [1406.5962].

A useful historical branch is the Shepard–Bernoulli construction, where compactly supported normalized Shepard weights are combined with local generalized Taylor polynomials on triangles [1406.5962]. Although not identical to the inverse-distance-product formulation of recent multinode schemes, it belongs to the same conceptual lineage: local polynomial exactness is lifted to the global interpolant through partition of unity. The bivariate Shepard–Bernoulli operator
\[
S_{B_m}[f](\boldsymbol{x})
= \sum_{i=1}^N \widetilde{W}_{\mu,i}(\boldsymbol{x})\, P_m^{\Delta_2(i)}[f](\boldsymbol{x})
\]
interpolates scattered data, reproduces bivariate polynomials up to degree \(m\), and has regularity controlled by the local weight smoothness [1406.5962].

Another related line is scaled Shepard approximation, where the base kernel is dilated in accordance with the fill distance of a quasi-uniform node set. That scheme is not interpolatory in general because the base kernel is bounded at the origin, but it yields an optimal-order Jackson-type inequality for bounded continuous functions on convex domains [1702.04764]. This suggests that the broader Shepard paradigm encompasses both interpolatory multinode constructions and quasi-interpolatory scaled variants, distinguished primarily by singular versus bounded kernels and by whether local polynomial enrichment is used.

Recent data-dependent Shepard interpolation modifies the nodewise shape parameter according to smoothness indicators in order to reduce smearing near jump discontinuities [2606.20332]. That work does not introduce local polynomial stencils into the final interpolant, so it is not a multinode Shepard method in the strict recent sense. Nevertheless, it is directly relevant conceptually because it shows that weight adaptation alone can improve interface behavior, and it suggests a plausible extension in which multinode Shepard weights inherit local smoothness-dependent shape parameters [2606.20332].

## 3. Approximation properties and polynomial reproduction

The most important theoretical advantage of multinode Shepard methods is polynomial reproduction beyond constants. If each local approximant \(p_j[f]\) belongs to a polynomial space reproduced exactly on its cluster, and the multinode weights form a partition of unity, then the global interpolant reproduces the same polynomial space [2508.07764, 2508.08023]. In the two-dimensional rectangular-grid specialization for tensor-product spaces,
\[
\mathbb{P}_r(\mathbb{R})\otimes\mathbb{P}_s(\mathbb{R}),
\]
the local interpolants are tensor-product polynomials and the global method reproduces all monomials \(x^i y^j\) with \(0\le i\le r\), \(0\le j\le s\) [2508.07764].

For digital elevation models on regular Cartesian grids, the approximation analysis uses Stancu’s remainder formula for tensor-product Lagrange interpolation. Under the assumptions
\[
f\in C^{r+1,s+1}(\Omega),\qquad
u>\frac{3+r+s}{t},\qquad t=(r+1)(s+1),
\]
the rectangular-grid multinode Shepard interpolant satisfies the uniform estimate [2508.07764]
\[
\|f-\mathcal{MS}_u[f]\|_\infty
\le \Delta\, C\, l_{\max}^{\delta_{\min}},
\]
where
\[
\Delta=\max\{\|f^{(r+1,0)}\|_\infty,\|f^{(0,s+1)}\|_\infty,\|f^{(r+1,s+1)}\|_\infty\},
\qquad
\delta_{\min}=\min\{r+1,s+1\}.
\]
For uniform refinement, this becomes
\[
\|f-\mathcal{MS}_u[f]\|_\infty = O(h^{\min\{r+1,s+1\}}),
\]
matching the expected tensor-product polynomial order [2508.07764].

The grid-based tensor-product collocation variant developed later for elliptic problems on Cartesian grids establishes analogous local interpolation estimates. If each local subgrid has \(n_x=m_x+1\), \(n_y=m_y+1\) nodes and supports tensor-product Lagrange interpolation, then on each patch [2606.14503]
\[
\|u-\mathcal I_k[u]\|_{\Omega_k}
\le C_k\left(h^{m_x+1}+h^{m_y+1}+h^{m_x+m_y+2}\right),
\]
and in the symmetric case \(m_x=m_y=m\), \(h<1\),
\[
\|u-\mathcal I_k[u]\|_{\Omega_k}\le \mathcal C_k\,h^{m+1}.
\]
The same work proves exact global reproduction of \(\mathbb P_{m_x,m_y}\) and of its first and second derivatives under the grid-based blended operator [2606.14503]. This suggests a strong structural continuity between interpolation-oriented multinode Shepard schemes and PDE-oriented collocation variants.

In the univariate fractional-derivative setting, exactness is formulated differently. When the exact solution of a Bagley–Torvik boundary or initial value problem is a polynomial of degree \(r\le d\), and the local interpolants have degree \(d\), the multinode Shepard collocation scheme recovers the exact polynomial solution under a nonsingularity or full-rank assumption on the assembled linear system [2508.08067]. This is a direct consequence of polynomial reproduction by the underlying interpolant.

## 4. Structured-grid and tensor-product formulations

A major recent specialization concerns regular rectangular grids, motivated by DEM reconstruction and later by elliptic PDE collocation. On a Cartesian grid
\[
\mathcal{X}_m\times\mathcal{Y}_n=\{x_1<\dots<x_m\}\times\{y_1<\dots<y_n\},
\]
the data are naturally arranged in an elevation matrix
\[
\mathcal Z=[z_{ij}],\qquad z_{ij}=f(x_i,y_j),
\]
and the grid is covered by overlapping rectangular blocks of size \((r+1)\times (s+1)\) [2508.07764]. Under the divisibility assumption
\[
\bmod(m-1,r)=\bmod(n-1,s)=0,
\]
the blocks are
\[
\sigma_{k,\ell}
= \{x_{(k-1)r+1},\dots,x_{kr+1}\}
\times
\{y_{(\ell-1)s+1},\dots,y_{\ell s+1}\}.
\]
Each block supports a unique tensor-product interpolant \(p_{k,\ell}\in \mathbb P_r\otimes \mathbb P_s\), often represented in a barycenter-centered Taylor-like form, with coefficients computed by a Vandermonde-type linear system solved by Gaussian elimination with partial pivoting [2508.07764].

The resulting grid-based multinode Shepard interpolant is
\[
\mathcal{MS}_u[\cdot](x,y)
=
\sum_{k,\ell} W_{u,k,\ell}(x,y)\, p_{k,\ell}[\cdot](x,y),
\]
where
\[
W_{u,k,\ell}(x,y)=
\frac{
\displaystyle\prod_{i=1}^{r+1}\prod_{j=1}^{s+1}
\bigl((x-x_{(k-1)r+i})^2+(y-y_{(\ell-1)s+j})^2\bigr)^{-u/2}
}{
\displaystyle\sum_{\kappa,\lambda}
\prod_{i=1}^{r+1}\prod_{j=1}^{s+1}
\bigl((x-x_{(\kappa-1)r+i})^2+(y-y_{(\lambda-1)s+j})^2\bigr)^{-u/2}
}.
\]
These weights are nonnegative, sum to one, and vanish at grid nodes outside the associated block [2508.07764].

A further development, the Grid-Based Multinode Shepard Collocation method, retains the same conceptual structure but uses overlapping sliding Cartesian subgrids, tensor-product Lagrange interpolation, and global basis functions explicitly assembled from local tensor-product basis functions and multinode weights [2606.14503]. In that setting, the global approximation has the form
\[
\mathcal G[u,\mathcal X](\mathbf x)
=
\sum_{k=1}^M \mathcal W_k(\mathbf x)\,\mathcal I_k[u](\mathbf x)
=
\sum_{\lambda=1}^N u_\lambda\,B_\lambda(\mathbf x),
\]
with
\[
\mathcal W_k(\mathbf x)
=
\frac{\omega_k(\mathbf x)}{\sum_{l=1}^M \omega_l(\mathbf x)},
\qquad
\omega_k(\mathbf x)=\prod_{\mathbf x'\in\mathcal X_k}\|\mathbf x-\mathbf x'\|^{-\mu},
\qquad \mu>2.
\]
This formulation yields a cardinal basis \(B_\lambda\), exact nodal interpolation, exact tensor-product polynomial reproduction, and sparse derivative matrices well suited to elliptic PDE discretization [2606.14503].

A plausible implication is that tensor-product grid structure removes one of the central implementation burdens of scattered-data multinode Shepard methods, namely the local unisolvency search. On Cartesian grids, the polynomial spaces and local node layouts are fixed a priori, so the method trades geometric flexibility for algorithmic regularity.

## 5. Scattered-data, meshfree, and collocation variants

In scattered-data settings, the multinode Shepard method is often interpreted as a meshfree partition-of-unity technique. For the two-dimensional Black–Scholes equation, the spatial approximation uses local quadratic polynomials on six-node stencils, blended by multinode Shepard basis functions [2508.08023]. The operator is written as
\[
\mathcal M_\mu[f](\mathbf x)=\sum_{j=1}^m B_{\mu,j}(\mathbf x)\,p_j[f](\mathbf x),
\]
with
\[
\tilde B_{\mu,j}(\mathbf x)=\prod_{\ell=1}^{\tau}\|\mathbf x-\mathbf x_{j_\ell}\|^{-\mu},
\qquad
B_{\mu,j}(\mathbf x)=
\frac{\tilde B_{\mu,j}(\mathbf x)}{\sum_{k=1}^m \tilde B_{\mu,k}(\mathbf x)},
\qquad
\tau=\binom{s+p}{s}.
\]
For \(s=2\), \(p=2\), one has \(\tau=6\), so each stencil contains six nodes and supports a quadratic bivariate polynomial [2508.08023].

The differentiability of the rational weights and local polynomials allows the direct computation of first and second derivatives of the approximation, making the method suitable for collocation of second-order PDE operators such as the two-dimensional Black–Scholes operator. The paper assembles a semi-discrete system
\[
\frac{d}{dt}\mathbf P_I(t)=\mathcal A\,\mathbf P_I(t)+\mathcal B\,\mathbf P_B(t),
\]
followed by BDF time discretization, with sparse matrices induced by locality of the basis [2508.08023]. Numerical experiments on several node configurations show that the method, denoted MS-FD in space and BDF2 in time, achieves mean errors around \(10^{-4}\) and maximum errors around \(10^{-2}\) in the reported tests, with condition numbers from about \(2.8\) to \(33.3\) depending on the node distribution [2508.08023].

The univariate multinode Shepard method for Caputo fractional derivatives uses overlapping subsets \(F_k\subset X\), local Lagrange polynomials \(P_k[f]\), and subset-level weights
\[
B_{\mu,k}(x)=
\frac{\prod_{i=1}^{p}|x-x_{k_i}|^{-\mu}}
{\sum_{\ell=1}^{s}\prod_{i=1}^{p}|x-x_{\ell_i}|^{-\mu}},
\]
leading to
\[
\mathcal M_\mu[f](x)=\sum_{k=1}^s B_{\mu,k}(x)\,P_k[f](x)
=\sum_{i=1}^n g_i(x)\,f_i.
\]
This nodal representation is then differentiated and inserted into a Gauss–Jacobi approximation of the Caputo derivative, enabling collocation for Bagley–Torvik equations [2508.08067]. The paper reports high accuracy for fractional derivative approximation and exactness for polynomial solutions of the associated boundary and initial value problems under the stated rank conditions [2508.08067].

These PDE-oriented formulations make clear that the multinode Shepard method is not limited to interpolation in the narrow sense. It also functions as a rational, polynomially reproducing basis-generation mechanism for derivative approximation, collocation, and operator discretization.

## 6. Applications: digital elevation models, finance, fractional equations, and discontinuities

The rectangular-grid DEM application is one of the most detailed empirical demonstrations of the method. The test case is a real digital terrain model of a mountainous area in the Sierra Nevada, Granada, Spain, with reference cell size \(2\) m, domain \(5420\ \text{m}\times 4886\ \text{m}\), elevation range \(1055\) m, and average slope \(28\%\) with range from \(0\%\) to \(81.5\%\) [2508.07764]. Coarsened \(8\) m and \(16\) m grids are reconstructed by a biquadratic multinode Shepard operator with \(r=s=2\) and \(u=4\), then evaluated back on the fine \(2\) m grid.

The reported vertical and horizontal discrepancies against cubic B-spline resampling are summarized below.

| Scenario | Multinode Shepard | Cubic B-spline |
|---|---:|---:|
| Vertical mean error, \(8\) m to \(2\) m | \(0.277\) m | \(0.596\) m |
| Vertical mean error, \(16\) m to \(2\) m | \(0.549\) m | \(0.886\) m |
| Horizontal mean discrepancy \(H_d\), \(8\) m to \(2\) m | \(0.227\) m | \(0.521\) m |
| Horizontal mean discrepancy \(H_d\), \(16\) m to \(2\) m | \(0.461\) m | \(0.790\) m |

The same study also reports lower mean and standard deviation than bilinear interpolation in both vertical and horizontal accuracy, especially in the \(8\) m case [2508.07764]. Visual assessment through hill-shade and 3D perspective indicates that the reconstructed surfaces are almost indistinguishable from the original at the tested resolutions [2508.07764].

In finance, the multinode Shepard collocation method is applied to the two-dimensional Black–Scholes equation on a triangular domain with near-field and far-field boundary conditions [2508.08023]. Across several node distributions, the multinode Shepard scheme outperforms an RBF-FD comparison in mean error and often yields comparable or slightly smaller maximum error. For example, on a Halton-based node set at the final time step, the reported errors are approximately \(1.97\times 10^{-4}\) mean and \(1.15\times 10^{-2}\) max for MS-FD, versus \(1.59\times 10^{-3}\) mean and \(3.42\times 10^{-2}\) max for RBF-FD [2508.08023].

In fractional calculus, the univariate multinode Shepard method is used to approximate Caputo derivatives of functions such as \(\sin x\), \(x^{9/2}\), and \(e^{2x}\), and then to solve Bagley–Torvik equations. The paper states that pointwise errors are typically between \(10^{-20}\) and \(10^{-6}\) for local polynomial degree \(d=8\), with improved accuracy as \(d\) increases and competitive conditioning for mixed node distributions [2508.08067].

A different application domain is quasi-histopolation of functions with jumps. There, multinode Shepard functions are used to blend local histopolation polynomials on overlapping intervals in a way that suppresses both Runge and Gibbs phenomena [2508.07070]. The global operator is
\[
\mathcal Q_{\mathcal H,\mu}[f](x)
=
\sum_{\iota=1}^{M} W_{\mu,\iota}(x)\,p_\iota(x),
\qquad
\mu\in 2\mathbb Z_+,
\]
where each \(p_\iota\) satisfies local integral constraints and each \(W_{\mu,\iota}\) is a multinode rational weight based on \(K\) subnodes inside an interval \(U_\iota\) [2508.07070]. The method is exact on polynomials up to degree \(d\), is \(C^\infty\) on the designated open set \(\Omega\), and numerically defeats the oscillations that plague classical global histopolation on equispaced grids [2508.07070].

## 7. Implementation issues, parameter choices, and limitations

Multinode Shepard methods are governed by a small number of recurring design choices: local polynomial degree, stencil or cluster geometry, overlap pattern, and weight exponent. In the DEM setting, the main parameters are \(r,s\) and \(u\), with the theoretical requirement
\[
u>\frac{3+r+s}{t},\qquad t=(r+1)(s+1),
\]
for the approximation theorem, जबकि the experiments use the biquadratic choice \(r=s=2\), \(t=9\), and \(u=4\), well above the theoretical threshold \(7/9\) [2508.07764]. In the Black–Scholes application, the fixed choices are polynomial degree \(p=2\), shape parameter \(\mu=4\), BDF2 time stepping, and Leja-selected local stencils of size six drawn from \(6+q\) nearest neighbors, typically with \(q=10\) [2508.08023]. In the fractional setting, \(\mu=4\) is again a common choice, and local polynomial degree is varied to control accuracy and conditioning [2508.08067].

A recurrent computational issue is numerical stabilization of the weights. Because the raw weights are products of many inverse powers of distances, overflow and underflow are natural risks. The grid-based collocation variant explicitly addresses this by working with logarithms,
\[
\lambda_k(\mathbf x)=\log \omega_k(\mathbf x),
\]
subtracting the local maximum, optionally clipping, exponentiating, and then normalizing [2606.14503]. This suggests that any implementation on large stencils or dense grids benefits from a log-sum-exp style normalization.

Another practical issue is locality. The raw definitions are global because every weight is normalized over all clusters or blocks. Yet the inverse-distance-product structure implies fast decay of remote contributions when the exponent is sufficiently large. Several papers note that distant blocks have negligible influence in practice and can be ignored to accelerate evaluation [2508.07764, 2606.14503]. A plausible implication is that scalable implementations on large DEMs or dense collocation grids should use localized block screening or neighborhood truncation, even though the analytic formulas are globally normalized.

The main limitations are structural. The rigorous DEM error estimate assumes smoothness \(f\in C^{r+1,s+1}(\Omega)\), a regular block covering satisfying the divisibility condition, and a sufficiently large weight exponent [2508.07764]. The Black–Scholes collocation work is empirically strong but does not provide a full PDE convergence theorem in the same level of detail [2508.08023]. The grid-based elliptic collocation method is presently restricted to two-dimensional rectangular domains with Cartesian grids, although it is presented as conceptually extendable to higher dimensions [2606.14503]. The quasi-histopolation method assumes knowledge of jump locations or reliable smoothness indicators for their detection [2508.07070].

A separate but relevant limitation concerns noisy data. The DEM paper notes that its method does not explicitly address regularization, whereas the cubic B-spline comparison in QGIS uses Tikhonov regularization, which may be advantageous for noisy measurements [2508.07764]. This suggests that classical multinode Shepard formulations are primarily interpolatory or exact-constraint methods; for noisy inverse problems, regularized local models or adaptive weight mechanisms such as data-dependent shape parameters may be desirable [2606.20332].

## 8. Generalizations and conceptual position

The multinode Shepard method now occupies a broad methodological space between classical inverse-distance interpolation, partition-of-unity meshfree approximation, local polynomial reconstruction, and collocation-based PDE discretization. The general framework is inherently multivariate and can be defined on scattered data in \(\mathbb R^d\), but recent work has shown that strong specializations emerge on structured grids, where tensor-product interpolation removes unisolvency search and yields sparse, moderately conditioned collocation matrices [2606.14503].

Several generalizations are already explicit in the literature. The original multivariate multinode formulation encompasses triangular and hexagonal variants in two dimensions for scattered data [2508.07764]. Shepard–Bernoulli operators show that the local approximants need not be ordinary Lagrange polynomials; generalized Taylor polynomials on triangles also fit naturally within the same partition-of-unity framework [1406.5962]. Fractional collocation demonstrates that the local polynomially reproducing rational basis can be differentiated and integrated inside nonlocal operators [2508.08067]. Quasi-histopolation shows that the local objects being blended need not even be nodal interpolants; they may instead satisfy segment integral conditions [2508.07070].

The most interesting emerging direction is adaptation. Data-dependent Shepard interpolation modifies nodewise shape parameters according to smoothness indicators and proves that the smearing belt near discontinuities can be narrowed relative to classical Shepard by approximately a factor \(1/\kappa\) under the stated hypotheses [2606.20332]. Although that work is not itself a multinode polynomial scheme, it strongly suggests that adaptive shape-parameter control could be combined with multinode local polynomial patches. This suggests a hybrid research direction: retain polynomial reproduction in smooth regions while sharpening locality near discontinuities.

In encyclopedic terms, the multinode Shepard method is best understood not as a single operator but as a design principle. Its invariant ingredients are a partition-of-unity Shepard-type rational weighting and local approximants constructed from multiple nodes. Around that principle, current research branches into structured-grid interpolation, scattered-data approximation, rational quasi-histopolation, financial PDE collocation, fractional differential equations, and adaptive discontinuity-aware weighting [2508.07764, 2508.08023, 2508.07070, 2508.08067, 2606.14503, 2606.20332]. The common outcome is higher algebraic exactness and stronger locality than classical Shepard interpolation, achieved without abandoning the normalized rational blending that defines the Shepard family.

Source: https://www.emergentmind.com/topics/multinode-shepard-method