Weighted Essentially Non-Oscillatory Shepard

• WENO–Shepard is a nonlinear interpolation technique that integrates WENO’s adaptive weighting into the classical Shepard approach for scattered data.
• It employs least squares linear fits and local smoothness indicators to adjust weights, ensuring accurate capture of sharp transitions.
• The method achieves first-order convergence in smooth regions while effectively confining overshoot and diffusion near discontinuities.

The Weighted Essentially Non-Oscillatory Shepard (WENO–Shepard) method is a nonlinear interpolation technique designed to improve the accuracy and stability of the classical Shepard method for scattered data, particularly in the presence of discontinuities. By embedding the core ideas of the WENO interpolation framework into the Shepard approach, the method adaptively suppresses the influence of data near discontinuities and efficiently captures sharp transitions while maintaining desirable properties in smooth regions. The method, as formalized in “Weighted Essentially Non-Oscillatory Shepard method” by Levin, Ramón, Ruiz-Álvarez, and Yáñez, provides a robust interpolant for scientific computing and data analysis tasks involving high-dimensional, scattered, and possibly piecewise-smooth functions (Levin et al., 2024).

1. Classical Shepard Interpolation: Principles and Limitations

The classical Shepard interpolation constructs a smooth, partition-of-unity approximation for functions defined on scattered data. Given a dataset $\{x_i, f_i = f(x_i)\}_{i=1}^N$ in a bounded domain $\Omega \subset \mathbb{R}^n$ with fill-distance $h$, and a compactly supported, radial weight function $\omega(r)$, the unnormalized weights are $\omega_i(x) = \omega(\|x - x_i\| / h)$ and the normalized Shepard weights are

$$W_i(x) = \frac{\omega_i(x)}{\sum_{j=1}^N \omega_j(x)}.$$

The Shepard interpolant is

$I_S[f](x) = \sum_{i=1}^N W_i(x) f_i.$

This approach guarantees partition of unity, non-negativity, and exact reproduction of constants. For $C^\nu$ weights, $I_S$ inherits the same smoothness and achieves first-order convergence in the uniform norm for $C^1$ functions: $\Omega \subset \mathbb{R}^n$0 However, near sharp discontinuities (across codimension-1 surfaces), Shepard weights do not adapt, so nodes from both sides of a jump contribute, resulting in severe smearing (diffusion) of the discontinuity over a region of width $\Omega \subset \mathbb{R}^n$1.

2. Core Concepts of Weighted Essentially Non-Oscillatory (WENO) Interpolation

WENO schemes were devised to provide non-oscillatory approximations for hyperbolic partial differential equations, particularly near discontinuities. The methodology involves:

• Splitting the interpolant construction into several lower-order sub-stencils.
• Evaluating a local smoothness indicator $\Omega \subset \mathbb{R}^n$2 for each sub-stencil: $\Omega \subset \mathbb{R}^n$3 if the local data is smooth; $\Omega \subset \mathbb{R}^n$4 across jumps.
• Defining nonlinear weights: $\Omega \subset \mathbb{R}^n$5 where $\Omega \subset \mathbb{R}^n$6 are baseline (linear) weights, $\Omega \subset \mathbb{R}^n$7 prevents division by zero, and $\Omega \subset \mathbb{R}^n$8 tunes sensitivity.
• Normalizing: $\Omega \subset \mathbb{R}^n$9 This weighting scheme effectively suppresses oscillatory contributions from discontinuous substencils, localizing the support of the interpolant near smooth data.

3. Construction of the WENO–Shepard Interpolant

The WENO–Shepard method generalizes the WENO philosophy to the scattered data context. The process is as follows:

• Shepard base weights: For each data node $h$0, compute $h$1 and normalized $h$2 as above.
• Smoothness indicators: For each $h$3, select a local stencil $h$4 with $h$5 a constant. Fit a degree-1 polynomial $h$6 via least squares on $h$7, then set

$h$8

Here, $h$9 in $\omega(r)$0 regions but $\omega(r)$1 if $\omega(r)$2 straddles a jump.

• Nonlinear weights: Define the unnormalized nonlinear weights

$\omega(r)$3

with $\omega(r)$4, $\omega(r)$5 (so $\omega(r)$6 in standard WENO notation).

• Final weights and interpolant: Normalize

$\omega(r)$7

Nodes with stencils crossing a jump (large $\omega(r)$8) have $\omega(r)$9 as $\omega_i(x) = \omega(\|x - x_i\| / h)$0, thus suppressing overshoots and diffusion near discontinuities.

4. Algorithmic Steps for Implementation

The WENO–Shepard algorithm may be implemented as follows for scattered data in $\omega_i(x) = \omega(\|x - x_i\| / h)$1:

1. Set the global fill-distance $\omega_i(x) = \omega(\|x - x_i\| / h)$2 and constant $\omega_i(x) = \omega(\|x - x_i\| / h)$3. Thus, each stencil radius is $\omega_i(x) = \omega(\|x - x_i\| / h)$4.
2. Choose a radial compactly-supported weight $\omega_i(x) = \omega(\|x - x_i\| / h)$5 (examples: Wendland’s $\omega_i(x) = \omega(\|x - x_i\| / h)$6, $\omega_i(x) = \omega(\|x - x_i\| / h)$7).
3. For each data node $\omega_i(x) = \omega(\|x - x_i\| / h)$8:
   • Build the stencil $\omega_i(x) = \omega(\|x - x_i\| / h)$9.
   • Compute the least-squares linear fit $$W_i(x) = \frac{\omega_i(x)}{\sum_{j=1}^N \omega_j(x)}.$$0 and the corresponding $$W_i(x) = \frac{\omega_i(x)}{\sum_{j=1}^N \omega_j(x)}.$$1.
4. For each evaluation point $$W_i(x) = \frac{\omega_i(x)}{\sum_{j=1}^N \omega_j(x)}.$$2:
   • Compute unnormalized Shepard weights $$W_i(x) = \frac{\omega_i(x)}{\sum_{j=1}^N \omega_j(x)}.$$3.
   • Compute $$W_i(x) = \frac{\omega_i(x)}{\sum_{j=1}^N \omega_j(x)}.$$4, then $$W_i(x) = \frac{\omega_i(x)}{\sum_{j=1}^N \omega_j(x)}.$$5.
   • Normalize $$W_i(x) = \frac{\omega_i(x)}{\sum_{j=1}^N \omega_j(x)}.$$6.
   • Return $$W_i(x) = \frac{\omega_i(x)}{\sum_{j=1}^N \omega_j(x)}.$$7.

Discontinuities across hypersurfaces $$W_i(x) = \frac{\omega_i(x)}{\sum_{j=1}^N \omega_j(x)}.$$8 are handled automatically; nodes with stencils intersecting $$W_i(x) = \frac{\omega_i(x)}{\sum_{j=1}^N \omega_j(x)}.$$9 have large $I_S[f](x) = \sum_{i=1}^N W_i(x) f_i.$0, resulting in negligible weight.

5. Analytical Guarantees: Smoothness, Accuracy, and Behavior Near Discontinuities

Three main theoretical results characterize the properties of the WENO–Shepard method (Levin et al., 2024):

• Smoothness Preservation: If the selected kernel $I_S[f](x) = \sum_{i=1}^N W_i(x) f_i.$1, then $I_S[f](x) = \sum_{i=1}^N W_i(x) f_i.$2.
• Convergence in Smooth Regions: For $I_S[f](x) = \sum_{i=1}^N W_i(x) f_i.$3 and quasi-uniform data with fill-distance $I_S[f](x) = \sum_{i=1}^N W_i(x) f_i.$4,

$I_S[f](x) = \sum_{i=1}^N W_i(x) f_i.$5

matching the convergence of classical Shepard.

• Suppression of Spurious Diffusion: If $I_S[f](x) = \sum_{i=1}^N W_i(x) f_i.$6 is piecewise $I_S[f](x) = \sum_{i=1}^N W_i(x) f_i.$7 with a discontinuity across a smooth hypersurface $I_S[f](x) = \sum_{i=1}^N W_i(x) f_i.$8, for $I_S[f](x) = \sum_{i=1}^N W_i(x) f_i.$9 a distance at least $C^\nu$0 from $C^\nu$1,

$C^\nu$2

In contrast, classical Shepard diffuses the jump over a region of width $C^\nu$3; WENO–Shepard sharply confines the diffusion close to $C^\nu$4.

6. Empirical Performance: Numerical Results

Validation of the WENO–Shepard method is provided through two principal numerical experiments:

• Smooth Benchmark (Franke’s Function): Testing on both regular grids and scattered Halton nodes (using Wendland $C^\nu$5, $C^\nu$6 kernels and an appropriate shape parameter), both classical Shepard and WENO–Shepard show expected $C^\nu$7 convergence in both maximum absolute error (MAE) and root mean squared error (RMSE). WENO–Shepard matches or marginally improves classical errors in smooth zones.
• Piecewise-Smooth Benchmark (Discontinuous Functions): For test functions with a unit-height jump across a line, circle, or square boundary, constructed as

$C^\nu$8

classical Shepard interpolates severe jump smearing, while WENO–Shepard sharply demarcates the discontinuity. Graphical and tabular results display this suppression of non-physical diffusion near jumps.

7. Summary and Significance

The WENO–Shepard method unites the adaptability of WENO weights with the dimension-agnostic ease of the Shepard framework. By modulating the classical linear weights $C^\nu$9 with nonlinear terms $I_S$0, the construction maintains the classic convergence rate and smoothness on regular regions and nearly removes overshoot and shock smearing near discontinuities. The method is straightforward to implement in any dimension, requires only zero-degree (linear) least squares fits for $I_S$1, and operates robustly on scattered data (Levin et al., 2024). This approach yields a valuable interpolation tool for scientific computing and applied mathematics, especially when faithfully capturing localized features in sparsely sampled, high-dimensional, or discontinuous data is essential.