---
title: Univariate Multinode Shepard Method
url: https://www.emergentmind.com/topics/univariate-multinode-shepard-method
type: topic
---

# Univariate Multinode Shepard Method

Searching arXiv for the specified paper and closely related multinode Shepard work.
arXiv_search: query="2508.08067 OR \"The univariate multinode Shepard method for the Caputo fractional derivatives\" OR multinode Shepard Bagley-Torvik", max_results=5
arXiv_search: query="2508.07070 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]\) [2508.08067].

## 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 [2508.08067].

The node set is
\[
X=\{x_1,\dots,x_n\}\subset[0,T],\qquad x_1<\cdots<x_n,
\]
together with a covering
\[
\mathcal F=\{F_1,\dots,F_s\},
\]
where each block \(F_k\subset X\) contains exactly \(p\) 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-\((p-1)\) 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 \(\mu>0\), the multinode Shepard weight functions are
\[
B_{\mu,k}(x)=
\frac{\displaystyle \prod_{i=1}^p |x-x_{k_i}|^{-\mu}}
{\displaystyle \sum_{\ell=1}^s \prod_{i=1}^p |x-x_{\ell_i}|^{-\mu}},
\qquad k=1,\ldots,s.
\]
These weights satisfy three structural properties: \(\sum_{k=1}^s B_{\mu,k}(x)\equiv 1\), \(B_{\mu,k}(x)\ge 0\), and each \(B_{\mu,k}\) vanishes at all nodes not in \(F_k\) [2508.08067].

If \(P_k[f]\) denotes the unique degree-\((p-1)\) Lagrange interpolant on \(F_k\), the global operator is
\[
\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(x_i),
\]
with coefficients
\[
g_i(x)=\sum_{k:\,x_i\in F_k} B_{\mu,k}(x)\,\ell_{k,i}(x).
\]
The operator interpolates nodal data,
\[
\mathcal M_\mu[f](x_j)=f(x_j),
\]
and reproduces every polynomial of degree at most \(p-1\),
\[
\mathcal M_\mu[p](x)=p(x).
\]

Repeated differentiation yields
\[
\bigl(\mathcal M_\mu[f]\bigr)^{(m)}(x)=\sum_{i=1}^n g_i^{(m)}(x)\,f(x_i).
\]
The paper states that closed-form but somewhat bulky formulas for \(g_i^{(m)}\) can be written in terms of \(B_{\mu,k}\), \(\ell_{k,i}\), and their derivatives. This identifies the differentiated multinode Shepard operator as the immediate surrogate for \(f^{(m)}\) inside the Caputo integral.

## 3. Approximation of Caputo fractional derivatives

The approximation strategy begins with
\[
f^{(m)}(t)\approx \bigl(\mathcal M_\mu[f]\bigr)^{(m)}(t)
=\sum_{i=1}^n g_i^{(m)}(t)\,f(x_i),
\]
which leads to
\[
D^\alpha f(x)\approx \frac{1}{\Gamma(m-\alpha)}
\sum_{i=1}^n f(x_i)\int_0^x (x-t)^{m-\alpha-1} g_i^{(m)}(t)\,dt.
\]
The remaining integral is handled by a Gauss–Jacobi quadrature after the change of variables \(t=\tfrac{x}{2}(1+\xi)\), \(\xi\in[-1,1]\) [2508.08067].

Under that transformation,
\[
\int_0^x (x-t)^{m-\alpha-1}u(t)\,dt
=
\Bigl(\tfrac{x}{2}\Bigr)^{m-\alpha}
\int_{-1}^1 (1-\xi)^{m-\alpha-1}
u\!\Bigl(\tfrac{x}{2}(1+\xi)\Bigr)\,d\xi.
\]
Gauss–Jacobi quadrature of order \(N\) is then applied with weight \((1-\xi)^{m-\alpha-1}(1+\xi)^0\). If \(\{\xi_k\}_{k=1}^N\) and \(\{w_k\}_{k=1}^N\) are the corresponding Jacobi nodes and weights for parameters \((m-\alpha-1,0)\), the resulting rule is
\[
\int_0^x (x-t)^{m-\alpha-1}u(t)\,dt
\approx
\Bigl(\tfrac{x}{2}\Bigr)^{m-\alpha}
\sum_{k=1}^N w_k^{(m-\alpha-1,0)}
u\!\Bigl(\tfrac{x}{2}(1+\xi_k^{(m-\alpha-1,0)})\Bigr).
\]

Substituting \(u(t)=g_i^{(m)}(t)\) gives the fully discrete approximation
\[
D^\alpha f(x)\approx
\frac{1}{\Gamma(m-\alpha)}
\Bigl(\tfrac{x}{2}\Bigr)^{m-\alpha}
\sum_{i=1}^n f(x_i)
\sum_{k=1}^N w_k^{(m-\alpha-1,0)}
\,g_i^{(m)}\!\Bigl(\tfrac{x}{2}(1+\xi_k)\Bigr).
\]
This construction makes the fractional derivative approximation depend only on nodal values \(f(x_i)\), the local polynomial structure encoded in \(g_i^{(m)}\), 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 \(g_i\), computing the derivatives \(g_i^{(m)}\) 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 [2508.08067].

For a target point \(x^\ast\), the implementation proceeds by computing, for each block, the numerator
\[
\prod_{x_{k_i}\in F_k}|x^\ast-x_{k_i}|^{-\mu},
\]
normalizing to obtain \(B_{\mu,k}(x^\ast)\), forming
\[
g_i(x^\ast)=\sum_{k:\,x_i\in F_k} B_{\mu,k}(x^\ast)\,\ell_{k,i}(x^\ast),
\]
and then evaluating \(g_i^{(m)}\) at the quadrature points
\[
t_k=\frac{x^\ast}{2}(1+\xi_k).
\]
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 \(\mathcal M_\mu\) 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,
\[
\rho\,y''(x)+\lambda\,D^\alpha y(x)+\sigma\,y(x)=h(x),\qquad x\in(0,T),
\]
with
\[
y(0)=\gamma_1,\qquad y(T)=\gamma_2,
\]
the approximation is sought in the form
\[
\tilde y(x)=\mathcal M_\mu[y](x)=\sum_{i=1}^n g_i(x)\,y_i,
\qquad y_i=y(x_i).
\]
Collocation is imposed at the interior nodes \(x_j\), \(j=2,\dots,n-1\),
\[
\rho\,\tilde y''(x_j)+\lambda\,D^\alpha \tilde y(x_j)+\sigma\,\tilde y(x_j)=h(x_j).
\]
Together with the boundary conditions \(\tilde y(x_1)=\gamma_1\) and \(\tilde y(x_n)=\gamma_2\), this yields an \((n-2)\times(n-2)\) linear system for \(y_2,\dots,y_{n-1}\) [2508.08067].

For the initial-value problem,
\[
\rho\,y''(x)+\lambda\,D^\alpha y(x)+\sigma\,y(x)=h(x),\qquad
y(0)=\gamma_1,\quad y'(0)=\gamma_2,
\]
the second boundary condition is replaced by
\[
\sum_{i=1}^n g_i'(0)\,y_i=\gamma_2-g_1'(0)\,\gamma_1.
\]
The resulting overdetermined \((n\times n)\) 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 \(\phi(x)\) of a fractional BVP or IVP is a polynomial of degree \(\le p-1\), then the multinode-Shepard collocation solution \(\tilde y(x)=\mathcal M_\mu[y]\) coincides exactly with \(\phi(x)\) at all \(x\). The proof sketch relies on polynomial reproduction, exact reduction of derivatives and fractional-derivative approximations to integrals of \(\phi\), 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 [2508.08067].

Specifically, if \(f\in C^p\) and the maximum fill of the covering is \(h\), then
\[
\|f-\mathcal M_\mu[f]\|_{L^\infty}\le C\,h^p\|f^{(p)}\|_{L^\infty},
\]
where \(C\) depends on \(\mu\) and \(p\). Combined with smoothness of quadrature node distributions,
\[
\bigl\|D^\alpha f-D^\alpha \mathcal M_\mu[f]\bigr\|_{L^\infty}=O(h^p),
\]
provided \(\mu\) is chosen large enough to guarantee \(p\) derivatives of each \(B_{\mu,k}\).

The numerical experiments reported in the same work show exact recovery to machine precision whenever the true solution is a polynomial of degree \(\le p-1\); exponential-like decrease of the maximum and mean absolute errors as the local polynomial degree increases for smooth nonpolynomial test functions such as \(\sin x\), \(e^{2x}\), and \(x^{9/2}\); robustness with respect to node layouts, including equispaced, blended equispaced–Chebyshev, and mock-Chebyshev nodes; and stable conditioning of the collocation matrix for moderate \(n\).

A related 2025 development applies multinode Shepard functions as blending functions in a \(C^\infty\) 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 [2508.07070].

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