Multinode Shepard Method: Theory & Applications

• Multinode Shepard method is an interpolation scheme that uses local polynomial approximants blended with normalized inverse-distance weights to achieve higher polynomial reproduction.
• It is effectively applied in digital elevation model reconstruction, financial PDE collocation, and fractional derivative approximation to deliver improved accuracy compared to classical methods.
• The approach forms unisolvent clusters of nodes, ensuring enhanced derivative behavior and maintaining a partition-of-unity structure for reliable global approximation.

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 (Barrera et al., 11 Aug 2025), to meshfree collocation for two-dimensional Black–Scholes equations (Dell'Accio et al., 11 Aug 2025), to quasi-histopolation with jump discontinuities (Dell'Accio et al., 9 Aug 2025), and to univariate approximation of Caputo fractional derivatives and Bagley–Torvik equations (Dell'Accio et al., 11 Aug 2025). Related developments include triangular Shepard–Bernoulli constructions for scattered bivariate data (Dell'Accio et al., 2014), scaled Shepard approximation on quasi-uniform sets (Senger et al., 2017), grid-based tensor-product collocation on Cartesian grids (Harrak et al., 12 Jun 2026), and data-dependent adaptive weighting for reduced smearing near discontinuities (Kuruc et al., 18 Jun 2026).

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 (Barrera et al., 11 Aug 2025): $$\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 (Barrera et al., 11 Aug 2025). 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 (Dell'Accio et al., 11 Aug 2025): $$\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 (Dell'Accio et al., 11 Aug 2025).

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 (Barrera et al., 11 Aug 2025).

2. Relation to classical Shepard and modified Shepard families

Classical Shepard interpolation in $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$0 has the form

$$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$1

or, more generally, a normalized weighted average with a radial kernel (Barrera et al., 11 Aug 2025, Senger et al., 2017). 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 (Dell'Accio et al., 2014).

A useful historical branch is the Shepard–Bernoulli construction, where compactly supported normalized Shepard weights are combined with local generalized Taylor polynomials on triangles (Dell'Accio et al., 2014). 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

$$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$2

interpolates scattered data, reproduces bivariate polynomials up to degree $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$3, and has regularity controlled by the local weight smoothness (Dell'Accio et al., 2014).

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 (Senger et al., 2017). 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 (Kuruc et al., 18 Jun 2026). 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 (Kuruc et al., 18 Jun 2026).

3. Approximation properties and polynomial reproduction

The most important theoretical advantage of multinode Shepard methods is polynomial reproduction beyond constants. If each local approximant $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$4 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 (Barrera et al., 11 Aug 2025, Dell'Accio et al., 11 Aug 2025). In the two-dimensional rectangular-grid specialization for tensor-product spaces,

$$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$5

the local interpolants are tensor-product polynomials and the global method reproduces all monomials $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$6 with $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$7, $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$8 (Barrera et al., 11 Aug 2025).

For digital elevation models on regular Cartesian grids, the approximation analysis uses Stancu’s remainder formula for tensor-product Lagrange interpolation. Under the assumptions

$$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$9

the rectangular-grid multinode Shepard interpolant satisfies the uniform estimate (Barrera et al., 11 Aug 2025)

$r$0

where

$r$1

For uniform refinement, this becomes

$r$2

matching the expected tensor-product polynomial order (Barrera et al., 11 Aug 2025).

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 $r$3, $r$4 nodes and supports tensor-product Lagrange interpolation, then on each patch (Harrak et al., 12 Jun 2026)

$r$5

and in the symmetric case $r$6, $r$7,

$r$8

The same work proves exact global reproduction of $r$9 and of its first and second derivatives under the grid-based blended operator (Harrak et al., 12 Jun 2026). 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 $t=t(d,r)=\binom{d+r}{r}.$0, and the local interpolants have degree $t=t(d,r)=\binom{d+r}{r}.$1, the multinode Shepard collocation scheme recovers the exact polynomial solution under a nonsingularity or full-rank assumption on the assembled linear system (Dell'Accio et al., 11 Aug 2025). 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

$t=t(d,r)=\binom{d+r}{r}.$2

the data are naturally arranged in an elevation matrix

$t=t(d,r)=\binom{d+r}{r}.$3

and the grid is covered by overlapping rectangular blocks of size $t=t(d,r)=\binom{d+r}{r}.$4 (Barrera et al., 11 Aug 2025). Under the divisibility assumption

$t=t(d,r)=\binom{d+r}{r}.$5

the blocks are

$t=t(d,r)=\binom{d+r}{r}.$6

Each block supports a unique tensor-product interpolant $t=t(d,r)=\binom{d+r}{r}.$7, 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 (Barrera et al., 11 Aug 2025).

The resulting grid-based multinode Shepard interpolant is

$t=t(d,r)=\binom{d+r}{r}.$8

where

$t=t(d,r)=\binom{d+r}{r}.$9

These weights are nonnegative, sum to one, and vanish at grid nodes outside the associated block (Barrera et al., 11 Aug 2025).

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 (Harrak et al., 12 Jun 2026). In that setting, the global approximation has the form

$p_j[f]$0

with

$p_j[f]$1

This formulation yields a cardinal basis $p_j[f]$2, exact nodal interpolation, exact tensor-product polynomial reproduction, and sparse derivative matrices well suited to elliptic PDE discretization (Harrak et al., 12 Jun 2026).

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 (Dell'Accio et al., 11 Aug 2025). The operator is written as

$p_j[f]$3

with

$p_j[f]$4

For $p_j[f]$5, $p_j[f]$6, one has $p_j[f]$7, so each stencil contains six nodes and supports a quadratic bivariate polynomial (Dell'Accio et al., 11 Aug 2025).

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

$p_j[f]$8

followed by BDF time discretization, with sparse matrices induced by locality of the basis (Dell'Accio et al., 11 Aug 2025). Numerical experiments on several node configurations show that the method, denoted MS-FD in space and BDF2 in time, achieves mean errors around $p_j[f]$9 and maximum errors around $$\mathcal{MS}_u[f](\boldsymbol{x})=\sum_{j=1}^q W_{u,j}(\boldsymbol{x})\,p_j[f](\boldsymbol{x}),$$0 in the reported tests, with condition numbers from about $$\mathcal{MS}_u[f](\boldsymbol{x})=\sum_{j=1}^q W_{u,j}(\boldsymbol{x})\,p_j[f](\boldsymbol{x}),$$1 to $$\mathcal{MS}_u[f](\boldsymbol{x})=\sum_{j=1}^q W_{u,j}(\boldsymbol{x})\,p_j[f](\boldsymbol{x}),$$2 depending on the node distribution (Dell'Accio et al., 11 Aug 2025).

The univariate multinode Shepard method for Caputo fractional derivatives uses overlapping subsets $$\mathcal{MS}_u[f](\boldsymbol{x})=\sum_{j=1}^q W_{u,j}(\boldsymbol{x})\,p_j[f](\boldsymbol{x}),$$3, local Lagrange polynomials $$\mathcal{MS}_u[f](\boldsymbol{x})=\sum_{j=1}^q W_{u,j}(\boldsymbol{x})\,p_j[f](\boldsymbol{x}),$$4, and subset-level weights

$$\mathcal{MS}_u[f](\boldsymbol{x})=\sum_{j=1}^q W_{u,j}(\boldsymbol{x})\,p_j[f](\boldsymbol{x}),$$5

leading to

$$\mathcal{MS}_u[f](\boldsymbol{x})=\sum_{j=1}^q W_{u,j}(\boldsymbol{x})\,p_j[f](\boldsymbol{x}),$$6

This nodal representation is then differentiated and inserted into a Gauss–Jacobi approximation of the Caputo derivative, enabling collocation for Bagley–Torvik equations (Dell'Accio et al., 11 Aug 2025). 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 (Dell'Accio et al., 11 Aug 2025).

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 $$\mathcal{MS}_u[f](\boldsymbol{x})=\sum_{j=1}^q W_{u,j}(\boldsymbol{x})\,p_j[f](\boldsymbol{x}),$$7 m, domain $$\mathcal{MS}_u[f](\boldsymbol{x})=\sum_{j=1}^q W_{u,j}(\boldsymbol{x})\,p_j[f](\boldsymbol{x}),$$8, elevation range $$\mathcal{MS}_u[f](\boldsymbol{x})=\sum_{j=1}^q W_{u,j}(\boldsymbol{x})\,p_j[f](\boldsymbol{x}),$$9 m, and average slope $$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.$$0 with range from $$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.$$1 to $$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.$$2 (Barrera et al., 11 Aug 2025). Coarsened $$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.$$3 m and $$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.$$4 m grids are reconstructed by a biquadratic multinode Shepard operator with $$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.$$5 and $$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.$$6, then evaluated back on the fine $$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.$$7 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, $$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.$$8 m to $$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.$$9 m	$p_j[f]$0 m	$p_j[f]$1 m
Vertical mean error, $p_j[f]$2 m to $p_j[f]$3 m	$p_j[f]$4 m	$p_j[f]$5 m
Horizontal mean discrepancy $p_j[f]$6, $p_j[f]$7 m to $p_j[f]$8 m	$p_j[f]$9 m	$$\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.$$0 m
Horizontal mean discrepancy $$\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.$$1, $$\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.$$2 m to $$\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.$$3 m	$$\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.$$4 m	$$\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.$$5 m

The same study also reports lower mean and standard deviation than bilinear interpolation in both vertical and horizontal accuracy, especially in the $$\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.$$6 m case (Barrera et al., 11 Aug 2025). Visual assessment through hill-shade and 3D perspective indicates that the reconstructed surfaces are almost indistinguishable from the original at the tested resolutions (Barrera et al., 11 Aug 2025).

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 (Dell'Accio et al., 11 Aug 2025). 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 $$\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.$$7 mean and $$\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.$$8 max for MS-FD, versus $$\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.$$9 mean and $$W_{\mu,i}(\mathbf{x}_k)=\delta_{ik},\qquad \sum_{i=1}^n W_{\mu,i}(\mathbf{x})=1,$$0 max for RBF-FD (Dell'Accio et al., 11 Aug 2025).

In fractional calculus, the univariate multinode Shepard method is used to approximate Caputo derivatives of functions such as $$W_{\mu,i}(\mathbf{x}_k)=\delta_{ik},\qquad \sum_{i=1}^n W_{\mu,i}(\mathbf{x})=1,$$1, $$W_{\mu,i}(\mathbf{x}_k)=\delta_{ik},\qquad \sum_{i=1}^n W_{\mu,i}(\mathbf{x})=1,$$2, and $$W_{\mu,i}(\mathbf{x}_k)=\delta_{ik},\qquad \sum_{i=1}^n W_{\mu,i}(\mathbf{x})=1,$$3, and then to solve Bagley–Torvik equations. The paper states that pointwise errors are typically between $$W_{\mu,i}(\mathbf{x}_k)=\delta_{ik},\qquad \sum_{i=1}^n W_{\mu,i}(\mathbf{x})=1,$$4 and $$W_{\mu,i}(\mathbf{x}_k)=\delta_{ik},\qquad \sum_{i=1}^n W_{\mu,i}(\mathbf{x})=1,$$5 for local polynomial degree $$W_{\mu,i}(\mathbf{x}_k)=\delta_{ik},\qquad \sum_{i=1}^n W_{\mu,i}(\mathbf{x})=1,$$6, with improved accuracy as $$W_{\mu,i}(\mathbf{x}_k)=\delta_{ik},\qquad \sum_{i=1}^n W_{\mu,i}(\mathbf{x})=1,$$7 increases and competitive conditioning for mixed node distributions (Dell'Accio et al., 11 Aug 2025).

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 (Dell'Accio et al., 9 Aug 2025). The global operator is

$$W_{\mu,i}(\mathbf{x}_k)=\delta_{ik},\qquad \sum_{i=1}^n W_{\mu,i}(\mathbf{x})=1,$$8

where each $$W_{\mu,i}(\mathbf{x}_k)=\delta_{ik},\qquad \sum_{i=1}^n W_{\mu,i}(\mathbf{x})=1,$$9 satisfies local integral constraints and each $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$00 is a multinode rational weight based on $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$01 subnodes inside an interval $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$02 (Dell'Accio et al., 9 Aug 2025). The method is exact on polynomials up to degree $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$03, is $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$04 on the designated open set $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$05, and numerically defeats the oscillations that plague classical global histopolation on equispaced grids (Dell'Accio et al., 9 Aug 2025).

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 $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$06 and $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$07, with the theoretical requirement

$$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$08

for the approximation theorem, जबकि the experiments use the biquadratic choice $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$09, $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$10, and $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$11, well above the theoretical threshold $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$12 (Barrera et al., 11 Aug 2025). In the Black–Scholes application, the fixed choices are polynomial degree $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$13, shape parameter $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$14, BDF2 time stepping, and Leja-selected local stencils of size six drawn from $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$15 nearest neighbors, typically with $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$16 (Dell'Accio et al., 11 Aug 2025). In the fractional setting, $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$17 is again a common choice, and local polynomial degree is varied to control accuracy and conditioning (Dell'Accio et al., 11 Aug 2025).

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,

$$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$18

subtracting the local maximum, optionally clipping, exponentiating, and then normalizing (Harrak et al., 12 Jun 2026). 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 (Barrera et al., 11 Aug 2025, Harrak et al., 12 Jun 2026). 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 $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$19, a regular block covering satisfying the divisibility condition, and a sufficiently large weight exponent (Barrera et al., 11 Aug 2025). The Black–Scholes collocation work is empirically strong but does not provide a full PDE convergence theorem in the same level of detail (Dell'Accio et al., 11 Aug 2025). 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 (Harrak et al., 12 Jun 2026). The quasi-histopolation method assumes knowledge of jump locations or reliable smoothness indicators for their detection (Dell'Accio et al., 9 Aug 2025).

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 (Barrera et al., 11 Aug 2025). 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 (Kuruc et al., 18 Jun 2026).

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 $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$20, 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 (Harrak et al., 12 Jun 2026).

Several generalizations are already explicit in the literature. The original multivariate multinode formulation encompasses triangular and hexagonal variants in two dimensions for scattered data (Barrera et al., 11 Aug 2025). 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 (Dell'Accio et al., 2014). Fractional collocation demonstrates that the local polynomially reproducing rational basis can be differentiated and integrated inside nonlocal operators (Dell'Accio et al., 11 Aug 2025). Quasi-histopolation shows that the local objects being blended need not even be nodal interpolants; they may instead satisfy segment integral conditions (Dell'Accio et al., 9 Aug 2025).

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 $$\{\sigma_j\}_{j=1}^q,\qquad \sigma_j=\{\boldsymbol{x}_{j_1},\dots,\boldsymbol{x}_{j_t}\},$$21 under the stated hypotheses (Kuruc et al., 18 Jun 2026). 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 (Barrera et al., 11 Aug 2025, Dell'Accio et al., 11 Aug 2025, Dell'Accio et al., 9 Aug 2025, Dell'Accio et al., 11 Aug 2025, Harrak et al., 12 Jun 2026, Kuruc et al., 18 Jun 2026). 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.