Univariate Multinode Shepard Method

• The univariate multinode Shepard method is a rational blending framework that combines local polynomial interpolants using partition-of-unity weights to form a global operator.
• It approximates Caputo fractional derivatives by replacing the function with a multinode Shepard interpolant and integrating via Gauss–Jacobi quadrature.
• The method effectively addresses Bagley–Torvik boundary and initial value problems through collocation, yielding stable and accurate linear systems.

Searching arXiv for the specified paper and closely related multinode Shepard work. arXiv_search: query="(Dell'Accio et al., 11 Aug 2025) OR \"The univariate multinode Shepard method for the Caputo fractional derivatives\" OR multinode Shepard Bagley-Torvik", max_results=5 arXiv_search: query="(Dell'Accio et al., 9 Aug 2025) OR \"multinode Shepard\" quasi-histopolation", max_results=5 The univariate multinode Shepard method is a rational blending framework in which local polynomial models on small node blocks are combined by partition-of-unity weights to produce a global interpolatory operator. In the formulation developed by Dell’Accio, Di Tommaso, and Ferrara, this operator is used to approximate Caputo fractional derivatives through Gauss–Jacobi quadrature and then to solve Bagley–Torvik boundary value problems and initial value problems numerically. In that setting, the method couples polynomial reproduction, differentiability controlled by a smoothness parameter, and collocation on a prescribed node set $X\subset[0,T]$ (Dell'Accio et al., 11 Aug 2025).

1. Analytical setting and definition

For a sufficiently smooth function $f:[0,T]\to\mathbb{R}$, fractional differentiation is taken in the Caputo sense. If $\alpha>0$ and $m=\lceil \alpha\rceil$, the derivative is defined by

$$D^\alpha f(x)=\frac{1}{\Gamma(m-\alpha)}\int_0^x (x-t)^{m-\alpha-1} f^{(m)}(t)\,dt,\qquad x>0.$$

The numerical problem addressed in the cited work is to approximate $D^\alpha f$ by replacing $f$ with a univariate multinode Shepard interpolant $\mathcal M_\mu[f]$, where $\mu>m$ ensures that derivatives up to order $m$ exist (Dell'Accio et al., 11 Aug 2025).

The node set is

$f:[0,T]\to\mathbb{R}$0

together with a covering

$f:[0,T]\to\mathbb{R}$1

where each block $f:[0,T]\to\mathbb{R}$2 contains exactly $f:[0,T]\to\mathbb{R}$3 nodes. The term “multinode” refers to this use of blocks rather than single nodes as the local supports of the construction. On each block, a degree-$f:[0,T]\to\mathbb{R}$4 Lagrange interpolant is formed; these local polynomials are then blended into a global operator.

This framework separates local approximation from global assembly. The local objects encode polynomial exactness on each block, while the global weights enforce interpolation and partition of unity. A plausible implication is that the method is designed to preserve local polynomial structure without sacrificing a globally coherent approximation.

2. Multinode Shepard operator and structural properties

For $f:[0,T]\to\mathbb{R}$5, the multinode Shepard weight functions are

$f:[0,T]\to\mathbb{R}$6

These weights satisfy three structural properties: $f:[0,T]\to\mathbb{R}$7, $f:[0,T]\to\mathbb{R}$8, and each $f:[0,T]\to\mathbb{R}$9 vanishes at all nodes not in $\alpha>0$0 (Dell'Accio et al., 11 Aug 2025).

If $\alpha>0$1 denotes the unique degree-$\alpha>0$2 Lagrange interpolant on $\alpha>0$3, the global operator is

$\alpha>0$4

with coefficients

$\alpha>0$5

The operator interpolates nodal data,

$\alpha>0$6

and reproduces every polynomial of degree at most $\alpha>0$7,

$\alpha>0$8

Repeated differentiation yields

$\alpha>0$9

The paper states that closed-form but somewhat bulky formulas for $m=\lceil \alpha\rceil$0 can be written in terms of $m=\lceil \alpha\rceil$1, $m=\lceil \alpha\rceil$2, and their derivatives. This identifies the differentiated multinode Shepard operator as the immediate surrogate for $m=\lceil \alpha\rceil$3 inside the Caputo integral.

3. Approximation of Caputo fractional derivatives

The approximation strategy begins with

$m=\lceil \alpha\rceil$4

which leads to

$m=\lceil \alpha\rceil$5

The remaining integral is handled by a Gauss–Jacobi quadrature after the change of variables $m=\lceil \alpha\rceil$6, $m=\lceil \alpha\rceil$7 (Dell'Accio et al., 11 Aug 2025).

Under that transformation,

$m=\lceil \alpha\rceil$8

Gauss–Jacobi quadrature of order $m=\lceil \alpha\rceil$9 is then applied with weight $$D^\alpha f(x)=\frac{1}{\Gamma(m-\alpha)}\int_0^x (x-t)^{m-\alpha-1} f^{(m)}(t)\,dt,\qquad x>0.$$0. If $$D^\alpha f(x)=\frac{1}{\Gamma(m-\alpha)}\int_0^x (x-t)^{m-\alpha-1} f^{(m)}(t)\,dt,\qquad x>0.$$1 and $$D^\alpha f(x)=\frac{1}{\Gamma(m-\alpha)}\int_0^x (x-t)^{m-\alpha-1} f^{(m)}(t)\,dt,\qquad x>0.$$2 are the corresponding Jacobi nodes and weights for parameters $$D^\alpha f(x)=\frac{1}{\Gamma(m-\alpha)}\int_0^x (x-t)^{m-\alpha-1} f^{(m)}(t)\,dt,\qquad x>0.$$3, the resulting rule is

$$D^\alpha f(x)=\frac{1}{\Gamma(m-\alpha)}\int_0^x (x-t)^{m-\alpha-1} f^{(m)}(t)\,dt,\qquad x>0.$$4

Substituting $$D^\alpha f(x)=\frac{1}{\Gamma(m-\alpha)}\int_0^x (x-t)^{m-\alpha-1} f^{(m)}(t)\,dt,\qquad x>0.$$5 gives the fully discrete approximation

$$D^\alpha f(x)=\frac{1}{\Gamma(m-\alpha)}\int_0^x (x-t)^{m-\alpha-1} f^{(m)}(t)\,dt,\qquad x>0.$$6

This construction makes the fractional derivative approximation depend only on nodal values $$D^\alpha f(x)=\frac{1}{\Gamma(m-\alpha)}\int_0^x (x-t)^{m-\alpha-1} f^{(m)}(t)\,dt,\qquad x>0.$$7, the local polynomial structure encoded in $$D^\alpha f(x)=\frac{1}{\Gamma(m-\alpha)}\int_0^x (x-t)^{m-\alpha-1} f^{(m)}(t)\,dt,\qquad x>0.$$8, and a standard weighted quadrature rule.

4. Algorithmic realization and collocation

The computational outline consists of precomputing the local Lagrange basis on each block, evaluating multinode Shepard weights at each target point, assembling the global coefficients $$D^\alpha f(x)=\frac{1}{\Gamma(m-\alpha)}\int_0^x (x-t)^{m-\alpha-1} f^{(m)}(t)\,dt,\qquad x>0.$$9, computing the derivatives $D^\alpha f$0 at the quadrature nodes, and then applying the discrete fractional formula. The paper notes that one may store and reuse many of the intermediary weight-functions and their derivatives if multiple evaluations are needed (Dell'Accio et al., 11 Aug 2025).

For a target point $D^\alpha f$1, the implementation proceeds by computing, for each block, the numerator

$D^\alpha f$2

normalizing to obtain $D^\alpha f$3, forming

$D^\alpha f$4

and then evaluating $D^\alpha f$5 at the quadrature points

$D^\alpha f$6

The derivatives are assembled by product- and chain-rule formulas.

When the method is used for a fractional differential equation by collocation, the differential operator is written in terms of $D^\alpha f$7 and the quadrature-based fractional derivative. This yields a linear system for the unknown nodal values. The computational structure is therefore interpolation-based rather than based on a direct discretization of the Caputo integral kernel alone.

5. Application to the Bagley–Torvik equation

The principal application in the cited work is the Bagley–Torvik equation. For the boundary-value problem,

$D^\alpha f$8

with

$D^\alpha f$9

the approximation is sought in the form

$f$0

Collocation is imposed at the interior nodes $f$1, $f$2,

$f$3

Together with the boundary conditions $f$4 and $f$5, this yields an $f$6 linear system for $f$7 (Dell'Accio et al., 11 Aug 2025).

For the initial-value problem,

$f$8

the second boundary condition is replaced by

$f$9

The resulting overdetermined $\mathcal M_\mu[f]$0 linear system is solved in the least-squares sense.

The experimental results reported for these Bagley–Torvik problems confirm the method’s effectiveness, particularly in accurately approximating the equation for both BVPs and IVPs. The formulation is therefore both an approximation procedure for Caputo derivatives and a collocation mechanism for fractional differential equations.

6. Exactness, error behavior, and relation to adjacent Shepard constructions

The key theoretical statement is an exactness theorem: if the true solution $\mathcal M_\mu[f]$1 of a fractional BVP or IVP is a polynomial of degree $\mathcal M_\mu[f]$2, then the multinode-Shepard collocation solution $\mathcal M_\mu[f]$3 coincides exactly with $\mathcal M_\mu[f]$4 at all $\mathcal M_\mu[f]$5. The proof sketch relies on polynomial reproduction, exact reduction of derivatives and fractional-derivative approximations to integrals of $\mathcal M_\mu[f]$6, and the partition-of-unity property. This exactness should not be read as a statement for arbitrary smooth functions; for non-polynomial data, the paper instead appeals to standard Shepard-type error estimates (Dell'Accio et al., 11 Aug 2025).

Specifically, if $\mathcal M_\mu[f]$7 and the maximum fill of the covering is $\mathcal M_\mu[f]$8, then

$\mathcal M_\mu[f]$9

where $\mu>m$0 depends on $\mu>m$1 and $\mu>m$2. Combined with smoothness of quadrature node distributions,

$\mu>m$3

provided $\mu>m$4 is chosen large enough to guarantee $\mu>m$5 derivatives of each $\mu>m$6.

The numerical experiments reported in the same work show exact recovery to machine precision whenever the true solution is a polynomial of degree $\mu>m$7; exponential-like decrease of the maximum and mean absolute errors as the local polynomial degree increases for smooth nonpolynomial test functions such as $\mu>m$8, $\mu>m$9, and $m$0; robustness with respect to node layouts, including equispaced, blended equispaced–Chebyshev, and mock-Chebyshev nodes; and stable conditioning of the collocation matrix for moderate $m$1.

A related 2025 development applies multinode Shepard functions as blending functions in a $m$2 rational quasi-histopolation operator for bounded functions with jumps. There, local histopolation polynomials are blended on small patches by multinode Shepard weights, and the construction is presented as defeating both the Runge and Gibbs phenomena. This related use underscores that the univariate multinode Shepard idea is not restricted to fractional differential equations: it is a general rational blending architecture in which local polynomial information is assembled globally by partition-of-unity weights (Dell'Accio et al., 9 Aug 2025).