Multinode Shepard Interpolant

• Multinode Shepard interpolant is a partition-of-unity operator that blends local polynomial models from multiple node subsets using normalized inverse-distance weights.
• It achieves accurate reconstruction and derivative consistency by ensuring exact polynomial reproduction and leveraging overlapping local patches.
• Its versatility is demonstrated in applications such as scattered data interpolation, DEM reconstruction, PDE collocation, and fractional derivative approximation.

The multinode Shepard interpolant is a generalized Shepard-type approximation operator in which a global approximant is formed by blending local interpolants built on subsets of multiple nodes rather than raw nodal values. In its standard form, the construction starts from a node set $X\subset\Omega$, a family of local subsets or patches that are unisolvent for a prescribed polynomial space, and normalized rational weights obtained from inverse products of distances to the nodes of each patch. The resulting operator is simultaneously local, because each constituent interpolant depends on a small node subset, and global, because the normalized weights form a partition of unity. Recent work has developed multinode Shepard constructions for scattered-data interpolation, Cartesian-grid interpolation, digital elevation model reconstruction, elliptic and parabolic PDE collocation, fractional-derivative approximation, and quasi-histopolation; related kernel-based multilevel work has also produced a nodal/cardinal representation that is structurally similar at the level of “sum of nodal values times basis functions,” while remaining distinct from classical Shepard normalization (Dell'Accio et al., 11 Aug 2025).

1. Canonical operator and its local-to-global structure

In the scattered-data formulation used for the two-dimensional Black–Scholes equation, one assumes a node set

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$

together with a covering $\mathcal{T}=\{t_1,\dots,t_m\}$, where each subset

$t_j=\{\mathbf{x}_{j_k}\mid k=1,\dots,\tau\}$

contains exactly

$\tau=\binom{s+p}{s}$

points and is unisolvent for polynomial interpolation in $\mathbb{P}_p(\mathbb{R}^s)$. For each patch $t_j$, the local interpolant $p_j[f]\in\mathbb{P}_p(\mathbb{R}^s)$ is uniquely determined by the interpolation conditions

$$p_j[f](\mathbf{x}_{j_k})=f(\mathbf{x}_{j_k}),\qquad k=1,\dots,\tau.$$

The multinode Shepard operator is then

$$\mathcal{M}_{\mu}[f](\mathbf{x}) =\sum_{j=1}^{m}B_{\mu,j}(\mathbf{x})\,p_j[f](\mathbf{x}),\qquad \mu>0,$$

with multinode Shepard functions

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$0

These are normalized rational weights: each local factor is the inverse product of distances from the evaluation point to all nodes in the corresponding tuple, and the denominator enforces global normalization. The same structural pattern appears in the regular-grid DEM formulation, where local tensor-product interpolants are blended by

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$1

and in the univariate formulation, where blocks $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$2 of $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$3 nodes are blended by

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$4

In all of these cases, the essential mechanism is the same: local polynomial information is assembled through Shepard-type rational partitioning rather than through a single global interpolant (Dell'Accio et al., 11 Aug 2025).

2. Basis functions, cardinality, and polynomial reproduction

A central feature of multinode Shepard methods is that the local interpolants can be rewritten in barycentric or Lagrange form and then regrouped into a global nodal expansion. In the scattered-data Black–Scholes formulation,

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$5

where the local Lagrange polynomials satisfy

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$6

After rearrangement, the operator becomes

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$7

with

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$8

The functions $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$9 form the multinode cardinal Shepard basis. An analogous regrouping appears in the univariate setting,

$\mathcal{T}=\{t_1,\dots,t_m\}$0

where

$\mathcal{T}=\{t_1,\dots,t_m\}$1

The structural properties repeatedly emphasized in the literature are the following. First, the raw patch weights form a partition of unity: $\mathcal{T}=\{t_1,\dots,t_m\}$2 Second, they vanish at nodes that do not belong to the associated patch: $\mathcal{T}=\{t_1,\dots,t_m\}$3 Third, at nodal points, sums over the relevant patch family restore the interpolation conditions: $\mathcal{T}=\{t_1,\dots,t_m\}$4 For $\mathcal{T}=\{t_1,\dots,t_m\}$5, the gradients vanish at nodes,

$\mathcal{T}=\{t_1,\dots,t_m\}$6

and for $\mathcal{T}=\{t_1,\dots,t_m\}$7, so do the Hessians,

$\mathcal{T}=\{t_1,\dots,t_m\}$8

Polynomial reproduction is obtained because each local polynomial interpolant reproduces the target polynomial on its patch and the weights sum to one. In the univariate case, the operator is designed to reproduce polynomials up to degree $\mathcal{T}=\{t_1,\dots,t_m\}$9. In the grid-based tensor-product variant, if

$t_j=\{\mathbf{x}_{j_k}\mid k=1,\dots,\tau\}$0

then each local tensor-product interpolant is exact on $t_j=\{\mathbf{x}_{j_k}\mid k=1,\dots,\tau\}$1, and consequently

$t_j=\{\mathbf{x}_{j_k}\mid k=1,\dots,\tau\}$2

The same paper also states exact reproduction for first and second derivatives of such polynomials: $t_j=\{\mathbf{x}_{j_k}\mid k=1,\dots,\tau\}$3 and

$t_j=\{\mathbf{x}_{j_k}\mid k=1,\dots,\tau\}$4

This exactness underlies the use of multinode Shepard bases in collocation and reconstruction settings (Harrak et al., 12 Jun 2026).

3. Geometric variants and major formulations

The terminology “multinode Shepard interpolant” now covers several geometrically distinct constructions.

Formulation	Local approximation	Weight construction
Scattered-data multinode Shepard	Polynomial interpolants on unisolvent tuples $t_j=\{\mathbf{x}_{j_k}\mid k=1,\dots,\tau\}$5	Inverse products of distances over tuple nodes
Rectangular-grid multinode Shepard for DEMs	Tensor-product polynomial interpolants on overlapping $t_j=\{\mathbf{x}_{j_k}\mid k=1,\dots,\tau\}$6 blocks	Normalized products over block-node distances
Grid-Based Multinode Shepard Collocation Method	Tensor-product Lagrange interpolation on overlapping Cartesian subgrids	Shepard-type partition functions $t_j=\{\mathbf{x}_{j_k}\mid k=1,\dots,\tau\}$7 from inverse-distance products
Quasi-histopolation with multinode Shepard functions	Local histopolation polynomials on continuity-preserving patches	Rational multinode weights $t_j=\{\mathbf{x}_{j_k}\mid k=1,\dots,\tau\}$8 built from local node sets
Univariate multinode Shepard for fractional derivatives	Local polynomial interpolants $t_j=\{\mathbf{x}_{j_k}\mid k=1,\dots,\tau\}$9 on overlapping blocks $\tau=\binom{s+p}{s}$0	Multinode basis functions $\tau=\binom{s+p}{s}$1

On regular rectangular grids used for digital elevation models, the node set is

$\tau=\binom{s+p}{s}$2

with local overlapping blocks

$\tau=\binom{s+p}{s}$3

each containing

$\tau=\binom{s+p}{s}$4

nodes. The global interpolant is

$\tau=\binom{s+p}{s}$5

This specializes multinode Shepard interpolation to a Cartesian setting while retaining the characteristic partition-of-unity assembly (Barrera et al., 11 Aug 2025).

For elliptic PDEs on rectangles, the Grid-Based Multinode Shepard Collocation Method constructs local Cartesian subgrids

$\tau=\binom{s+p}{s}$6

with geometric supports

$\tau=\binom{s+p}{s}$7

and forms the global approximant

$\tau=\binom{s+p}{s}$8

where

$\tau=\binom{s+p}{s}$9

and

$\mathbb{P}_p(\mathbb{R}^s)$0

A stabilized logarithmic evaluation of the weights is also described in order to avoid overflow and underflow (Harrak et al., 12 Jun 2026).

For bounded integrable functions with jumps, the quasi-histopolation construction first splits the interval into continuity intervals and then blends local histopolation polynomials $\mathbb{P}_p(\mathbb{R}^s)$1 through multinode Shepard functions

$\mathbb{P}_p(\mathbb{R}^s)$2

yielding

$\mathbb{P}_p(\mathbb{R}^s)$3

The node sets in overlapping patches must be chosen carefully: the paper explicitly notes that if they do not share consistent points on overlaps, the weights may behave poorly there (Dell'Accio et al., 9 Aug 2025).

A nearby but distinct development appears in multilevel kernel interpolation. For nested data, the multilevel operator admits the nodal form

$\mathbb{P}_p(\mathbb{R}^s)$4

with a cardinal basis

$\mathbb{P}_p(\mathbb{R}^s)$5

This is not Shepard interpolation in the standard normalization-by-distance sense, but it is explicitly described as the closest analog in that work to a weighted cardinal interpolant written as a sum of nodal values times nodal basis functions. This suggests that the multinode Shepard paradigm belongs to a broader family of generalized nodal representations, while remaining distinguished by its specific rational inverse-distance normalization (Gollwitzer et al., 15 Jun 2026).

4. Approximation properties, smoothness, and error behavior

The DEM formulation provides an explicit approximation theorem. Assuming

$\mathbb{P}_p(\mathbb{R}^s)$6

defining

$\mathbb{P}_p(\mathbb{R}^s)$7

and requiring

$\mathbb{P}_p(\mathbb{R}^s)$8

the error satisfies

$\mathbb{P}_p(\mathbb{R}^s)$9

where

$t_j$0

The associated approximation order is therefore

$t_j$1

The proof is described as relying on the partition-of-unity property, the Stancu remainder formula for tensor-product interpolation, bounds on products of coordinate distances within each block, and a summability argument over blocks at increasing annular distance from the evaluation point (Barrera et al., 11 Aug 2025).

For the Cartesian-grid collocation formulation, the local tensor-product interpolant on a subgrid of size $t_j$2, with

$t_j$3

satisfies

$t_j$4

where

$t_j$5

If $t_j$6 and $t_j$7, this simplifies to

$t_j$8

The paper emphasizes that the global discretization inherits locality while preserving exact tensor-product polynomial reproduction (Harrak et al., 12 Jun 2026).

The quasi-histopolation formulation stresses smoothness and suppression of oscillatory artifacts. Its operator is described as $t_j$9, and the paper states that it reconstructs discontinuous functions without Gibbs-type oscillations while also defeating the Runge phenomenon. The global pointwise error bound has the form

$p_j[f]\in\mathbb{P}_p(\mathbb{R}^s)$0

with $p_j[f]\in\mathbb{P}_p(\mathbb{R}^s)$1 under grid refinement and $p_j[f]\in\mathbb{P}_p(\mathbb{R}^s)$2 controllable through the number $p_j[f]\in\mathbb{P}_p(\mathbb{R}^s)$3 of local nodes, provided

$p_j[f]\in\mathbb{P}_p(\mathbb{R}^s)$4

A corresponding segmentwise integral bound is also given: $p_j[f]\in\mathbb{P}_p(\mathbb{R}^s)$5 The paper also derives explicit upper bounds for neighboring multinode Shepard weights in terms of $p_j[f]\in\mathbb{P}_p(\mathbb{R}^s)$6, $p_j[f]\in\mathbb{P}_p(\mathbb{R}^s)$7, and patch separation, thereby quantifying off-patch decay (Dell'Accio et al., 9 Aug 2025).

In practical numerical studies, the reported error behavior is strongly problem-dependent but consistently tied to local polynomial degree, patch geometry, and node layout. This suggests that the interpolant’s performance is governed less by a single global smoothness mechanism than by the interaction among local unisolvency, overlap, and rational weight localization.

5. Computational and numerical performance

The trivariate partition-of-unity RBF method, although not a multinode Shepard scheme in the narrow polynomial-patch sense, is explicitly described as a Shepard-type interpolation with higher-order data because the blended objects are local approximants $p_j[f]\in\mathbb{P}_p(\mathbb{R}^s)$8 rather than raw point values. Its global interpolant is

$p_j[f]\in\mathbb{P}_p(\mathbb{R}^s)$9

with a possible normalized Shepard choice

$$p_j[f](\mathbf{x}_{j_k})=f(\mathbf{x}_{j_k}),\qquad k=1,\dots,\tau.$$0

where $$p_j[f](\mathbf{x}_{j_k})=f(\mathbf{x}_{j_k}),\qquad k=1,\dots,\tau.$$1 is the inverse Euclidean norm. The main computational contribution is a cube-partition search in $$p_j[f](\mathbf{x}_{j_k})=f(\mathbf{x}_{j_k}),\qquad k=1,\dots,\tau.$$2 that restricts neighbor inspection to at most

$$p_j[f](\mathbf{x}_{j_k})=f(\mathbf{x}_{j_k}),\qquad k=1,\dots,\tau.$$3

cubes for each local search. Sorting the $$p_j[f](\mathbf{x}_{j_k})=f(\mathbf{x}_{j_k}),\qquad k=1,\dots,\tau.$$4 nodes costs $$p_j[f](\mathbf{x}_{j_k})=f(\mathbf{x}_{j_k}),\qquad k=1,\dots,\tau.$$5, sorting the $$p_j[f](\mathbf{x}_{j_k})=f(\mathbf{x}_{j_k}),\qquad k=1,\dots,\tau.$$6 evaluation points costs $$p_j[f](\mathbf{x}_{j_k})=f(\mathbf{x}_{j_k}),\qquad k=1,\dots,\tau.$$7, and local RBF systems are treated as constant-sized under regular covering assumptions. Reported timings show the acceleration clearly:

• $$p_j[f](\mathbf{x}_{j_k})=f(\mathbf{x}_{j_k}),\qquad k=1,\dots,\tau.$$8: $$p_j[f](\mathbf{x}_{j_k})=f(\mathbf{x}_{j_k}),\qquad k=1,\dots,\tau.$$9 s vs. $$\mathcal{M}_{\mu}[f](\mathbf{x}) =\sum_{j=1}^{m}B_{\mu,j}(\mathbf{x})\,p_j[f](\mathbf{x}),\qquad \mu>0,$$0 s,
• $$\mathcal{M}_{\mu}[f](\mathbf{x}) =\sum_{j=1}^{m}B_{\mu,j}(\mathbf{x})\,p_j[f](\mathbf{x}),\qquad \mu>0,$$1: $$\mathcal{M}_{\mu}[f](\mathbf{x}) =\sum_{j=1}^{m}B_{\mu,j}(\mathbf{x})\,p_j[f](\mathbf{x}),\qquad \mu>0,$$2 s vs. $$\mathcal{M}_{\mu}[f](\mathbf{x}) =\sum_{j=1}^{m}B_{\mu,j}(\mathbf{x})\,p_j[f](\mathbf{x}),\qquad \mu>0,$$3 s,
• $$\mathcal{M}_{\mu}[f](\mathbf{x}) =\sum_{j=1}^{m}B_{\mu,j}(\mathbf{x})\,p_j[f](\mathbf{x}),\qquad \mu>0,$$4: $$\mathcal{M}_{\mu}[f](\mathbf{x}) =\sum_{j=1}^{m}B_{\mu,j}(\mathbf{x})\,p_j[f](\mathbf{x}),\qquad \mu>0,$$5 s vs. $$\mathcal{M}_{\mu}[f](\mathbf{x}) =\sum_{j=1}^{m}B_{\mu,j}(\mathbf{x})\,p_j[f](\mathbf{x}),\qquad \mu>0,$$6 s. For Franke-type test functions, RMSE values as low as $$\mathcal{M}_{\mu}[f](\mathbf{x}) =\sum_{j=1}^{m}B_{\mu,j}(\mathbf{x})\,p_j[f](\mathbf{x}),\qquad \mu>0,$$7 and $$\mathcal{M}_{\mu}[f](\mathbf{x}) =\sum_{j=1}^{m}B_{\mu,j}(\mathbf{x})\,p_j[f](\mathbf{x}),\qquad \mu>0,$$8 are reported at $$\mathcal{M}_{\mu}[f](\mathbf{x}) =\sum_{j=1}^{m}B_{\mu,j}(\mathbf{x})\,p_j[f](\mathbf{x}),\qquad \mu>0,$$9, depending on the kernel and test function (Cavoretto et al., 2014).

The DEM-oriented multinode Shepard method is evaluated on a real mountainous DTM from Sierra Nevada, Granada, Spain, using a $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$00 reference DEM and coarsened $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$01 and $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$02 grids. Against QGIS cubic B-spline interpolation, the reported mean absolute vertical discrepancies are $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$03 m versus $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$04 m for $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$05, and $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$06 m versus $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$07 m for $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$08. Against bilinear interpolation, the reported values are $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$09 m versus $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$10 m and $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$11 m versus $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$12 m, respectively. Mean horizontal discrepancies are also lower for the multinode Shepard reconstruction, approximately $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$13 m and $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$14 m in the two scenarios, compared with $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$15 m and $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$16 m for cubic B-splines and $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$17 m and $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$18 m for bilinear interpolation (Barrera et al., 11 Aug 2025).

For elliptic PDE collocation on Cartesian grids, sparsity and conditioning are major reported advantages. The derivative matrices produced by the Grid-Based Multinode Shepard Collocation Method are described as highly sparse, with reported sparsity levels around $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$19–$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$20 in many tests. Condition numbers are reported as typically $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$21–$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$22 for GBMSC, compared with often $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$23–$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$24 for Kansa RBF collocation; MSC is usually better than Kansa but often less well-conditioned than GBMSC (Harrak et al., 12 Jun 2026).

6. PDE, fractional, and operator-theoretic applications

In time-dependent PDEs, the multinode Shepard interpolant functions as a spatial discretization basis. For the two-dimensional Black–Scholes equation,

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$25

with

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$26

the approximate solution is written as

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$27

Spatial derivatives of the PDE operator are applied directly to the basis $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$28, producing collocation matrices $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$29 and $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$30 from values of $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$31 and their first and second derivatives. Time stepping uses Backward Euler for the first step and BDF2 thereafter: $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$32

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$33

The reported numerical findings include mean errors around $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$34 for MS-FD on Halton nodes, with RBF-FD mean errors around $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$35, and condition numbers about $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$36 at the first time step and $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$37 thereafter (Dell'Accio et al., 11 Aug 2025).

For elliptic boundary value problems on rectangles, GBMSC uses the global expansion

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$38

then collocates the differential operator at grid nodes. For Poisson’s equation,

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$39

the stiffness matrix is

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$40

with boundary rows replaced to enforce Dirichlet conditions. The paper identifies the cardinal property

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$41

as crucial for collocation because the nodal unknowns are exactly the coefficients of the global expansion (Harrak et al., 12 Jun 2026).

For fractional calculus, the univariate multinode Shepard operator is inserted into the Caputo derivative

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$42

After approximating $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$43 by $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$44, the resulting integral is evaluated by Gauss–Jacobi quadrature with order

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$45

The full approximation becomes

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$46

This is then used in collocation schemes for Bagley–Torvik boundary and initial value problems. The paper reports pointwise errors typically between $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$47 and $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$48 in fractional-derivative tests, and mean errors often around $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$49 to $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$50 for polynomial-solution Bagley–Torvik BVP examples, numerically confirming an exactness theorem for polynomial solutions (Dell'Accio et al., 11 Aug 2025).

Across these applications, the multinode Shepard interpolant is not merely a reconstruction formula. It acts as a basis-generating mechanism whose interpolation, differentiation, and polynomial-reproduction properties can be transferred directly into discretizations for PDEs, fractional operators, and inverse problems.

7. Relation to classical Shepard interpolation and recurring misconceptions

Classical Shepard interpolation is usually understood as a direct normalized inverse-distance average of nodal values. Multinode Shepard interpolation differs in a precise and now standard way: it replaces constants or point values by local approximants. In the scattered-data formulation, the global approximation blends local polynomials $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$51; in the trivariate PU-RBF formulation, it blends local RBF interpolants $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$52; in quasi-histopolation, it blends local histopolation polynomials; and in Cartesian-grid variants, it blends tensor-product Lagrange or tensor-product polynomial interpolants. For this reason, one paper explicitly remarks that the PU-RBF construction can be viewed as a Shepard-type interpolation with higher-order data, while also stating that it is more precisely a partition-of-unity RBF interpolant (Cavoretto et al., 2014).

A common misconception is that any nodal/cardinal representation of the form

$$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$53

is automatically a Shepard method. The multilevel kernel formulation shows why this is not correct. It produces a cardinal nodal representation with basis functions $$X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\subset \Omega\subset \mathbb{R}^s$$54, and in the nested case it is exactly a weighted nodal interpolant, but the basis is generated by kernel interpolation on multiple scales through a residual-correction algorithm rather than by normalized inverse-distance factors. The same source therefore states that it is not a Shepard interpolant in the standard normalization-by-distance sense, even though it is structurally comparable at the level of nodal expansion (Gollwitzer et al., 15 Jun 2026).

Another recurring distinction concerns scattered versus grid-based constructions. On scattered data, local unisolvency of each patch must be ensured explicitly. On Cartesian grids, tensor-product structure makes local unisolvency automatic, simplifies algebra, and yields sparse matrices with moderate conditioning. This is why recent grid-based papers present their methods as retaining the “good features of multinode Shepard constructions—locality, interpolation, and partition-of-unity assembly” while replacing scattered-node patch selection by overlapping Cartesian subgrids (Harrak et al., 12 Jun 2026).

In current usage, then, the multinode Shepard interpolant is best understood as a family of partition-of-unity rational operators built from multi-node local models. What remains invariant across formulations is the normalized blending of local approximants through inverse-distance products; what varies is the geometry of the patches, the local approximation space, and the downstream objective, ranging from surface reconstruction to PDE and fractional-operator discretization.