lebesgue

The Lebesgue function

recursivenodes.lebesgue.lebesgue(x, d, n, nodes, Vinv=None)[source]

The Lebesgue function in the biunit \(d\)-simplex with a choice of nodes for degree-\(n\) polynomials.

Notes

Let \(X_n = \{\boldsymbol{x}_i\}\) be the set of \(\binom{n+d}{d}\) nodes. The associated shape functions \(\{\phi_i\}\) are polynomials with degree at most \(n\) such that \(\phi_i(\boldsymbol{x}_j) = \delta_{ij}\).

The Lebesgue function is the sum of the absolute value of the shape functions,

\[l(\boldsymbol{x}) = \sum_{i=1}^{\binom{n+d}{d}} | \phi_i(\boldsymbol{x}) |.\]
Parameters:
  • x (ndarray) – Shape (\(N\), \(d\)), points at which to evaluate the Lebesgue function.

  • d (int) – The spatial dimension.

  • n (int) – The polynomial degree that the nodes are unisolvent for.

  • nodes (ndarray) – Shape (\(\binom{n+d}{d}\), \(d\)), the nodes defining the shape functions.

  • Vinv (ndarray, optional) – Shape (\(\binom{n+d}{d}\), \(\binom{n+d}{d}\)), the inverse of the Vandermonde matrix for the Proriol-Koornwinder-Dubiner basis evaluated at the nodes. These polynomials are used as a basis for stable evluation in the biunit \(d\)-simplex, so the inverse of the Vandermonde matrix must be computed if it is not already known.

Returns:

Shape (\(N\),), the value of the Lebesgue function for the nodes at the points x.

Return type:

ndarray

recursivenodes.lebesgue.lebesguemin(x, d, n, nodes, Vinv=None)[source]

\(-\sum_j l(\boldsymbol{x}_j)\) and its gradient for the Lebesgue function \(l\) defined by a set of nodes.

Parameters:
  • x (ndarray) – Shape (\(N\), \(d\)), points at which to evaluate the Lebesgue function.

  • d (int) – The spatial dimension.

  • n (int) – The polynomial degree that the nodes are unisolvent for.

  • nodes (ndarray) – Shape (\(\binom{n+d}{d}\), \(d\)), the nodes defining the shape functions.

  • Vinv (ndarray, optional) – See documentation for lebesgue().

Returns:

\(-\sum_j l(\boldsymbol{x}_j)\) and its gradient with respect to x, reshaped to a vector with shape (\(Nd\),). (This function and return format are intended for use by scipy.optimize.minimize.)

Return type:

(float, ndarray)

recursivenodes.lebesgue.lebesgueminhess(x, d, n, nodes, Vinv=None)[source]

The Hessian of \(-\sum_j l(\boldsymbol{x}_j)\) for the Lebesgue function \(l\) defined by a set of nodes.

Parameters:
  • x (ndarray) – Shape (\(N\), \(d\)), points at which to evaluate the Lebesgue function.

  • d (int) – The spatial dimension.

  • n (int) – The polynomial degree that the nodes are unisolvent for.

  • nodes (ndarray) – Shape (\(\binom{n+d}{d}\), \(d\)), the nodes defining the shape functions.

  • Vinv (ndarray, optional) – See documentation for lebesgue().

Returns:

A block diagonal matrix of shape \((Nd,Nd)\), the Hessian of the Lebesgue function with respect to x. (This function and return format are intended for use by scipy.optimize.minimize.)

Return type:

scipy.sparse.block_diag

The Lebesgue constant

recursivenodes.lebesgue.lebesguemax(d, n, nodes, Vinv=None, solver_options={'barrier_tol': 0.001, 'gtol': 0.001}, ex_nodes=None, logger=None)[source]

An estimate of the Lebesgue constant \(\Lambda_n^{\max}(\Delta^d)\), the maximum of the Lebesgue function on the biunit \(d\)-simplex for a set of nodes for polynomials up to degree \(n\).

Notes

If \(l\) is the Lebesgue function defined by the nodes (see lebesgue()), then

(1)\[\Lambda_n^{\max}(\Delta^d) = \max_{\boldsymbol{x}\in\Delta^d} l(\boldsymbol{x}).\]

This constant is significant in quantifying the interpolation quality of a set of nodes because of the error estimate

\[\| f - p^n(f) \|_{\infty} \leq (1 + \Lambda_n^{\max}) \inf_{q\in\mathcal{P}_n(\Delta^d)} \| f - q \|_{\infty},\]

where \(p^n(f)\) is the interpolant of \(f\) defined by the nodes.

\(\Lambda_n^{\max}(\Delta^d)\) is nonconvex and only piecewise smooth, with many local maxima. This function is designed to estimate \(\Lambda_n^{\max}(\Delta^d)\) when the nodes have the same lattice-like layout as equispaced nodes (as all but one of the node sets implemented in this package do, see nodes).

When this lattice-like structure is present, the maxima always occur in one of the “void”s between nodes in the lattice:

>>> import matplotlib.pyplot as plt
>>> from recursivenodes import recursive_nodes
>>> from recursivenodes.nodes import equispaced
>>> from recursivenodes.utils import coord_map
>>> from recursivenodes.lebesgue import lebesgue
>>> x = equispaced(2, 200, domain='biunit')
>>> nodes = recursive_nodes(2, 6, domain='biunit')
>>> l = lebesgue(x, 2, 6, nodes)
>>> x = coord_map(x, 'biunit', 'equilateral')
>>> ax = plt.gca()
>>> ax.tricontour(x[:,0], x[:, 1], l, 41)
<matplotlib.tri._tricontour.TriContourSet object at ...>
>>> ax.set_aspect('equal')
>>> plt.show()
_images/lebesgue-1.png

Also, though the Lebesgue function is only piecewise smooth, it is guaranteed to be smooth at the maxima. An efficient way to estimate the \(\Lambda_n^{\max}(\Delta^d)\) for nodes with lattice-like structure is to try to optimize the Lebesgue function from one initial guess at the centroid of each void. The node sets in this package are also fully symmetric, which means that initial guesses can be placed in only a fraction of the voids. Here, for instance, are the 3D degree 7 equispaced nodes and the initial guesses used by lebesguemax:

_images/lebesgue-2.png

This function uses scipy.optimize.minimize (method=’trust-constr’) to find the local maxima from the initial guesses. While it is not theoretically guaranteed that the global maximum is always found in this way, the values of \(\Lambda_n^{\max}(\Delta^d)\) that have been compute with lebesguemax agree with previously computed values in [BP06, LP06, War06].

Parameters:
  • d (int) – The spatial dimension.

  • n (int) – The polynomial degree that the nodes are unisolvent for.

  • nodes (ndarray) – Shape (\(\binom{n+d}{d}\), \(d\)), the nodes defining the shape functions.

  • Vinv (ndarray, optional) – See documentation for lebesgue().

  • solver_options (dict, optional) – Options to control the behavior of scipy.optimize.minimize.

  • ex_nodes (ndarray, optional) – if nodes includes nodes that do not touch the boundary of the simplex, ex_nodes is a set of nodes that includes extra boundary locations that can be used to help define initial guesses.

  • logger (logging.logger, optional) – Logger for progress information.

Returns:

The estimate and estimated argmax of (1).

Return type:

(float, ndarray)

recursivenodes.lebesgue.lebesguemax_reference(d, n, nodes='recursive')[source]

Look up previously computed Lebesgue constants (\(\Lambda_n^{\max}(\Delta^d)\)) for node sets implemented in this package and other node sets in the literature.

Parameters:
  • d (int) – The spatial dimension.

  • n (int) – The polynomial degree that the nodes are unisolvent for.

  • nodes (str, optional) –

    The name of the node set to query. Choices include:

    • 'recursive': recursivenodes.recursive_nodes().

    • 'equispaced': recursivenodes.nodes.equispaced().

    • 'warburton': recursivenodes.nodes.warburton().

    • 'blyth_luo_pozrikidis': recursivenodes.nodes.blyth_luo_pozrikidis().

    • 'Chen-Babuska': Chen-Babuška type nodes [CB95, CB96], computed in [Rot05].

    • 'Fekete': Fekete nodes, computed in [Rot05].

    • 'Lebesgue': Lebesgue-constant-minimizing nodes, computed in [Rot05].

    • 'Heinrichs-asymmetric': Lebesgue-constant-minimizing nodes [Hei05].

    • 'Heinrichs-symmetric': Symmetric Lebesgue-constant-minimizing nodes [Hei05].

    • 'Hesthaven': Hesthaven electrostatic nodes [Hes98].

    • 'Hesthaven-Teng': Hesthaven-Teng electrostatic nodes [Hes98].

    • 'Rapetti-Sommariva-Vianello-leb': Lebesgue-constant-minimizing nodes, computed in [RSV12].

    • 'Rapetti-Sommariva-Vianello-lebgl': Lebesgue-constant-minimizing nodes wih Lobatto-Gauss-Legendre edges, computed in [RSV12].

    • 'Rapetti-Sommariva-Vianello-lebgls': Symmetric Lebesgue-constant-minimizing nodes wih Lobatto-Gauss-Legendre edges, computed in [RSV12].

Returns:

\(\Lambda_n^{\max}(\Delta^d)\) for that node set.

Return type:

float

Raises:

KeyError – The node set is unknown, or \(\Lambda_n^{\max}(\Delta^d)\) has not been precomputed for it.