Spectral Methods with Bernstein Polynomials

• The paper presents a spectral method that constructs global approximants using Bernstein polynomials, enabling robust and structure-preserving numerical solutions.
• Key features include the advantageous endpoint interpolation, nonnegativity, and tensor-product extensions which support accurate multivariate discretizations.
• Implementation optimizations such as closed-form derivative evaluation and sparse matrix assembly lead to high efficiency, exponential convergence, and reliable performance.

A spectral method based on Bernstein polynomials is a numerical framework in which global approximants are constructed using the Bernstein polynomial basis, taking advantage of its distinctive analytic, geometric, and computational properties. This approach is increasingly leveraged in the spectral discretization of differential equations, integral equations, and operator approximation, providing a robust alternative to methods relying on orthogonal bases. The Bernstein basis is defined over a finite or infinite interval and can be systematically extended to multivariate and tensor-product settings, delivering structure-preserving and boundary-respecting approximations for smooth and nonsmooth problems.

1. Mathematical Foundation and Bernstein Basis Construction

The univariate Bernstein basis of degree $n$ on $[a, b]$ is given by: $$B_{i,n}(x) = \binom{n}{i} \frac{(x-a)^i (b-x)^{n-i}}{(b-a)^n},\quad 0 \leq i \leq n,$$ with the partition of unity property $\sum_{i=0}^n B_{i,n}(x) = 1$. Any function $f \in C([a, b])$ is approximated by its Bernstein expansion: $$B_n(f; x) = \sum_{i=0}^n f\left(a + \frac{i}{n} (b-a)\right) B_{i,n}(x).$$

For multivariate/multidimensional problems and tensor product domains $\Omega \subset \mathbb{R}^d$, the basis functions take the form: $$B_{i_1,n_1}(x_1) \cdots B_{i_d,n_d}(x_d),\qquad 0 \leq i_j \leq n_j,$$ yielding a tensor-product structure suitable for high-dimensional spectral discretizations, such as in collocation and Petrov-Galerkin schemes (Mirkov et al., 2012, Jani et al., 2016).

Key algebraic properties exploited in spectral methods include:

• Endpoint interpolation: $B_{0,n}(a) = 1,\, B_{n,n}(b) = 1$, with all intermediate indices vanishing at the interval endpoints.
• Nonnegativity and convex-hull property: The range of a polynomial in the Bernstein basis is contained in the convex hull of its coefficients.

2. Collocation, Galerkin, and Petrov-Galerkin Formulations

Spectral methods using the Bernstein basis typically proceed by representing the solution of a PDE or integral equation as a sum of Bernstein polynomials multiplied by unknown coefficients. In the collocation framework (Mirkov et al., 2012), the governing equation (e.g., a differential or integral operator $L$ acting on $u$) is enforced at a discrete set of collocation nodes. Substituting the Bernstein interpolant into $L[u](x)$ and evaluating at nodes yields a linear algebraic system: $A \mathbf{b} = \mathbf{c}.$ Boundary conditions, often Dirichlet or Neumann, are directly imposed by taking advantage of the endpoint interpolation property; coefficients corresponding to boundary basis functions can be eliminated or set to prescribed values, simplifying enforcement and reducing the number of unknowns.

For time-dependent or fractional PDEs, dual and modal bases constructed from linear combinations of Bernstein and dual Bernstein polynomials lead to well-conditioned Petrov-Galerkin schemes. Specifically, the dual-Petrov-Galerkin method (Jani et al., 2016), and modal basis construction (Jani et al., 2016), build test functions that ensure sparse (banded) matrices and compatibility with (fractional) boundary conditions. The operational matrices of differentiation exploit biorthogonality and recurrence relations for computational efficiency.

3. Numerical Efficiency, Derivatives, and Matrix Assembly

Bernstein-based spectral methods incorporate several algorithmic optimizations to improve computational efficiency:

• Non-recursive closed-form derivative evaluation: Higher-order derivatives of the Bernstein basis functions are calculated via explicit, non-recursive formulas, such as

$$D^p B_{i,n}(x) = \frac{n!}{(n-p)! (b-a)^p} \sum_{k=\max(0,i+p-n)}^{\min(i,p)} \binom{p}{k} B_{i-k, n-p}(x),$$

which eliminates the need for repeated differentiation (Mirkov et al., 2012).

• Efficient binomial coefficient calculation using a multiplicative formula rather than factorials, stabilizing computations for large polynomial degree $n$.
• Sparse and banded matrix structures: Operational (differentiation) matrices and mass matrices are constructed to be banded, exploiting the local support and recurrence of the dual/modal bases (Jani et al., 2016, Jani et al., 2016). This yields favorable sparsity for fast direct solvers (e.g., LU decomposition).

4. Accuracy, Convergence, and Robustness

Bernstein spectral methods exhibit high regularity and exponential or spectral convergence for smooth analytic solutions:

• Exponential decay of $L^2$ or maximum error norms with increasing polynomial degree is observed in a range of elliptic, Helmholtz, and biharmonic problems (Mirkov et al., 2012).
• In fractional PDEs, the fully discrete scheme is analytically shown to be unconditionally stable, with spatial errors decaying spectrally and temporal errors decaying at rate $O(\tau^{2-\alpha})$ for Caputo derivatives (Jani et al., 2016, Jani et al., 2016).
• For moderate degrees ($n\lesssim 20$), machine precision can be achieved, with error stagnation or deterioration only at very high degrees due to floating point limitations (large factorial growth).
• Bernstein polynomials do not exhibit the Gibbs phenomenon for functions with jumps; their uniform convergence property controls spurious oscillations and ensures monotonicity and total variation diminishing properties, which is further leveraged for shock capturing (Glaubitz, 2019).

5. Specialized Techniques: Bounds, Nonnegativity, and Integer Constraints

Variational extensions of the Bernstein spectral method enforce a variety of qualitative constraints directly on the expansion coefficients:

• Imposition of positivity or upper/lower bounds on the approximant is achieved by constraining the Bernstein coefficients (due to the convex hull property). Bounds-preserving spectral approximants are constructed via quadratic programming with linear inequality constraints (Allen et al., 2021).
• Mass conservation, or other linear equality constraints, is incorporated by additional Lagrange multipliers, resulting in KKT systems solvable by direct linear algebra (Allen et al., 2021).
• In integer-constrained spectral approximation, the so-called Bernstein lattice uses integer combinations of the basis. Approximants from the Bernstein lattice are shown to be uniformly dense in the space of continuous functions with integer values at endpoints, providing explicit sup-norm error bounds in terms of the modulus of continuity (Güntürk et al., 2023).

6. Implementation Strategies and Notable Software

The computational pipeline for Bernstein spectral methods typically includes:

• Assembly of collocation or Galerkin matrices with efficient handling of basis function evaluations and derivatives,
• Mapping between solution vectors and coefficient arrays,
• Modular codebases for problem set-up, including flexible specification of polynomial degree, spatial domain, boundary conditions, and right-hand sides.

Publicly available libraries, such as the Python package “bernstein-poly” (Mirkov et al., 2012), encapsulate these routines and offer rapid prototyping for high-order elliptic PDEs and extendibility to Navier–Stokes or more complex systems.

A typical implementation snippet:

n = 12 # polynomial order nd = linspace(x1, x2, n+1) ... def rhs(x,y): return ... # right hand side for i in range(0, n+1): x = nd[i] for j in range(0, n+1): y = nd[j] node = (n+1)*i + j if on_boundary(x, y): f[node] = boundary_value(x, y) else: f[node] = rhs(x, y)

This modular approach enables adaptation to a variety of boundary conditions and PDE operators with only minor adjustments.

7. Applications and Future Directions

Bernstein spectral frameworks have been successfully applied to:

• Elliptic boundary value problems (including Poisson, Helmholtz, biharmonic),
• Fractional diffusion and advection-dispersion equations (with Caputo derivatives),
• Shock capturing in conservation laws via structure-preserving reconstructions and TV-diminishing operators (Glaubitz, 2019),
• Optimal control (allowing direct transcription of state and costate trajectories), where convex hull and non-oscillatory properties are critical for safe and feasible controls (Cichella et al., 2018),
• Operator approximation (e.g., data-driven spectral approximation of the Koopman operator) (Yadav et al., 4 Mar 2024),
• Integer-constrained digital approximation schemes (Güntürk et al., 2023).

The research trajectory includes extensions to multidimensional and simplex geometries, adaptivity and domain decomposition (for very high accuracy and large domains), efficient inversion/preconditioning of mass matrices (Allen et al., 2019), and synthesis with convex optimization for enforcing sophisticated constraints (Allen et al., 2021, Lasserre, 2023).

Bernstein-based spectral methods thus deliver a comprehensive toolbox for high-order, robust, and structure-preserving numerical simulation across a broad class of operator equations and computational problems.