Discontinuous Piecewise Polynomial Approximation

• Discontinuous piecewise polynomial approximation is a method that assigns separate polynomials to distinct subdomains, enabling accurate modeling of functions with abrupt changes.
• It underpins discontinuous Galerkin methods, moving least squares techniques, and specialized neural network architectures for effectively handling non-smooth data.
• Error estimates and convergence rates are established in broken Sobolev spaces, ensuring robust performance for shock-capturing and complex geometry problems.

Discontinuous piecewise polynomial approximation is a class of numerical and analytic methods in which the domain is partitioned into subregions, and on each subregion, a distinct polynomial (possibly of varying order) is assigned. Across the interfaces between subregions, these polynomials are not generally constrained to be continuous, allowing for the accurate representation of functions with discontinuities. Such constructions are central to discontinuous Galerkin (dG) finite element methods, moving least squares with discontinuous weights, discontinuous neural network architectures, and various optimization and interpolation strategies in analysis and PDEs. The theory and practice of discontinuous piecewise polynomial approximation address challenges arising from discontinuities in the solution or the domain, including those associated with non-Lipschitz geometries, heterogeneous data, multi-phase materials, and shock-capturing in conservation laws.

1. Fundamental Constructions and Function Spaces

The canonical setting for discontinuous piecewise polynomial approximation consists of a domain Ω, possibly with complex or non-Lipschitz boundary, decomposed into a collection of open subdomains or mesh elements {K}. On each K, a polynomial p_K ∈ P_p(K) of total degree ≤ p is defined. The global approximation space takes the form

$$V_{h,p} = \{ v \in L^2(\Omega) : v|_K \in P_p(K) \ \forall K \},$$

where no inter-element continuity is enforced (Hewett, 27 Nov 2025, Angermann et al., 2011). This flexibility allows for the representation of functions with jumps or non-matching derivatives across interfaces.

On non-Lipschitz or even fractal domains, the approximation theory admits the use of meshes and subdomains whose boundaries are highly irregular. The only essential requirement is that each element is contained in some shape-regular, overlapping cover (Hewett, 27 Nov 2025). The natural norm setting is in broken Sobolev spaces, with local approximation and projection theorems (see §3) holding with the same rates as in the classical case, independent of geometric boundary regularity.

2. Discontinuous Polynomial Approximation in Galerkin and Mesh-free Methods

Discontinuous Galerkin (dG) methods employ spaces of discontinuous piecewise polynomials as a core ingredient. On each element T, a local polynomial of tensor-product or total degree ≤ k is defined, with the global space assembled without continuity constraints (Angermann et al., 2011). The lack of inter-element constraints permits exact recovery of sharp interfaces, adaptivity, and local refinement.

Advanced mesh-free methods also utilize discontinuous piecewise polynomial techniques to address discontinuities in the underlying function. In moving least squares (MLS) with variably scaled discontinuous weight functions, the discontinuity structure is explicitly encoded into the weights, ensuring that the approximation is locally polynomial on smooth patches while suppressing the influence of data across jump interfaces. The MLS-VSDK approximant is constructed via local weighted polynomial least squares, with the weight function defined using a partition of Ω and an auxiliary scaling function ψ(x) that jumps across discontinuities (Esfahani et al., 2023).

Methods based on polynomials defined in terms of an argmin over auxiliary variables further enable the approximation of multivariate discontinuous functions. In this framework, the function $f$ is recovered as the $y$-minimizer of a bivariate polynomial $p(x, y)$, constructed via sum-of-squares semidefinite programming over sample data, achieving exactness in the case of piecewise-polynomial targets and obviating the need for explicit partitioning or discontinuity location (Henrion et al., 2023).

3. Approximation Theory: Error Estimates and Convergence Rates

The approximation properties of discontinuous piecewise polynomials are governed by local best-approximation results, global projection theorems, and inverse inequalities in fractional Sobolev and broken polynomial spaces. For any s ∈ [0, m], and polynomial degree p ≥ m − 1, the $L^2$-orthogonal projector Π onto $V_{h,p}$ satisfies,

$$\|u - \Pi u\|_{L^2(\Omega)} \leq C \left( \frac{h}{p+1} \right)^s \|u\|_{H^s(\Omega)},$$

for all $u \in H^s(\Omega)$, with C independent of element shape or smoothness of the boundary (Hewett, 27 Nov 2025, Angermann et al., 2011). These estimates generalize in a standard way to $W^{s,p}$-norms.

When extended to functions with discontinuity interfaces, the convergence rates in the $L^2$-norm remain optimal (of order $h^{m+1}$ or $h^p$) on smooth subdomains, while the global error is controlled in a piecewise Sobolev norm $\mathcal{W}_p^{k}(\Omega)$ (Esfahani et al., 2023). In the presence of sufficiently regular discontinuity curves or surfaces, e.g., $C^{m+1}$ boundaries, discontinuous spline and MLS-VSDK methods guarantee robust piecewise rates, with loss only in a localized neighborhood of the interface (Levin, 10 May 2024).

In approximation of piecewise analytic functions on quasi-smooth arcs in $\mathbb C$, near-best polynomial approximants achieve error

$$\|f - P_n\|_{\mathcal{L}} \leq C E_n(f;\mathcal{L}),$$

globally, and at analytic points away from junctions, stretched-exponential convergence $O(E_n e^{-c n^\sigma})$ for any $0 < \sigma < 1$ (Kryvonos, 2021). Local rates are determined by the order of matching derivatives at the discontinuity and the arc’s conformal geometry.

4. Algorithms and Problem Classes: Segmentation, Optimization, and Neural Networks

Discontinuous piecewise polynomial approximation is broadly employed in computational problems involving segmentation, denoising, and constrained optimization where the solution is inherently discontinuous or defined on a domain with disjoint feasible regions.

Optimization with Disjoint Feasible Regions

The Piecewise Polynomial Interpolation (PPI) function approach replaces disjoint-interval constraints (prohibited operating zones) in nonlinear programming by a set of differentiable equality and inequality constraints in higher dimensions. For a variable $p_k$ restricted to a union of intervals, the PPI function $m(p_k, \alpha^k)$ is constructed as a sum of piecewise quadratics activated by Heaviside functions. The transformation $m=0$ together with $0 \leq \alpha_{k,i} \leq 1$ replaces the original, non-smooth constraints, rendering the problem amenable to standard gradient-based solvers—while guaranteeing $C^1$ smoothness and normalized gradients in the admissible set (Pinheiro et al., 2020).

Piecewise Polynomial Neural Architectures

Discontinuous neural networks employ connection functions that are polynomials of varying (typically low) degree, defined on compact intervals (sub-links) and identically zero outside these intervals. For each input, only the polynomial associated with the active interval is evaluated, yielding a highly sparse, discontinuous response function. Such architectures efficiently represent highly nonlinear and discontinuous mappings, with empirical results demonstrating that higher-order sub-link polynomials significantly reduce error on standard regression and classification benchmarks (Loverich, 2015).

Segmentation and Data Fitting

In the context of image segmentation and denoising, mixed integer programming (MIP) formulations fit a discontinuous piecewise-affine function to grid data by introducing binary variables to indicate discontinuities across neighboring data points. The resulting function is enforced to be affine within each segment by big-M constraints on second derivatives, and global multicuts ensure a valid, label-consistent partition. The framework generalizes to higher-degree polynomials by constraining higher-order derivatives (Shen et al., 2020).

A summary of representative algorithms:

Approach	Domain Class	Discontinuity Handling
dG finite elements (Angermann et al., 2011)	Polyhedral/Fractal	Mesh faces, no continuity
MLS-VSDK (Esfahani et al., 2023)	ℝ^d, partitions	Discontinuous weight function
Polynomial argmin (Henrion et al., 2023)	Compact, sample set	Argmin eliminates jumps
PPI for NLP (Pinheiro et al., 2020)	Disjoint intervals	Piecewise quadratic mapping
Discontinuous NN (Loverich, 2015)	Euclidean	Interval-polynomial links
MIP fitting (Shen et al., 2020)	Grid data	Binary edge variables

5. Numerical and Analytical Performance

In dG and broken element methods on non-Lipschitz domains, optimal convergence is observed in both $L^2$ and broken Sobolev norms. Explicit estimates show

$$\left( \sum_{K\in T_h} \|u - v\|_{W^j(K)}^2 \right)^{1/2} \leq C h^{t - j} (p+1)^{j - m} \|u\|_{W^m(\Omega)}$$

with $t = \min(m, p+1)$ and $j = 0, ..., m$ (Hewett, 27 Nov 2025). In practical implementations, the reconstructed one-DOF-per-element discontinuous polynomial spaces achieve the same error as classical dG spaces with reduced degrees of freedom, and optimal multigrid-preconditioned solvers exhibit iteration counts and condition numbers independent of mesh size (Li et al., 2023).

Mesh-free MLS-VSDK achieves second-order or higher rates O($h^{\ell+1}$) on each smooth patch, as demonstrated numerically on benchmark 1D and 2D interface problems, even under interface perturbations or noise. Classical MLS stagnates at first order near jumps, while VSDK-based routines restore full convergence in regular subdomains (Esfahani et al., 2023).

In the argmin SOS polynomial approach, the method is mesh-free, requires only isolated sample values, and achieves exact recovery with low-degree bivariate polynomials for genuinely piecewise-polynomial targets, without suffering Gibbs oscillations. Empirical results on discontinuous indicators demonstrate this absence of ringing, in contrast to classical polynomial fits (Henrion et al., 2023).

6. Extensions, Generalizations, and Open Problems

Recent advances in discontinuous piecewise polynomial approximation have focused on several directions:

• Fractal and rough domains: Theoretical estimates have been extended to arbitrary non-Lipschitz, even fractal, geometries with only minimal covering assumptions (Hewett, 27 Nov 2025).
• Fractional and negative Sobolev bounds: All standard projection error results adapt via interpolation to $H^s$ with $s$ non-integer and negative $s$ under weak geometric conditions.
• Mesh adaptation and hp-refinement: Covering meshes and local projection strategies enable hp-optimal rates even when element shapes are degenerate or non-standard.
• Algorithmic extensions: Generalization to higher-order piecewise polynomials in combinatorial algorithms (e.g., MIP) is achieved via analogous higher-order smoothness constraints (Shen et al., 2020).

Open directions include sharp quantification of constants in approximation theorems relative to the Hausdorff or Assouad dimension of domain and interface, precise characterization of interpolation scales in $W^{s,p}$ for fully fractal boundaries, and geometric multigrid solvers exploiting self-similarity in mesh refinement (Hewett, 27 Nov 2025).

Discontinuous piecewise polynomial approximation, through its theoretical robustness and algorithmic flexibility, remains central to the numerical treatment of singular, heterogeneous, or sharply varying solutions across computational mathematics, scientific computing, signal processing, and learning theory.