# recursivenodes: Recursive, parameter-free, explicitly defined interpolation nodes for simplices¶

Contents:

## What is this?¶

recursivenodes.recursive_nodes(d, n, family='lgl', domain='barycentric')

Recursively defined nodes for $$\mathcal{P}_n(\Delta^d)$$, the polynomials with degree at most $$n$$ on the $$d$$-simplex, based on a 1D node family.

Notes

Some definitions:

• A 1D node set $$X_k=(x_{k,0}, \dots, x_{k,k})$$ is a sorted list of $$k+1$$ points in $$[0,1]$$ that is symmetric about $$1/2$$.

• A 1D node family $$\boldsymbol{X} = \{X_k\}$$ is a collection of 1D node sets for every degree $$k$$.

• The barycentric triangle is a canonical domain for the $$d$$-simplex in $$\mathbb{R}^{d+1}$$: positive coordinates $$\boldsymbol{b} = (b_0, b_1, \dots b_d)$$ such that $$\sum_i b_i = 1$$.

• For a multi-index $$\boldsymbol{\alpha}$$, let $$\#\boldsymbol{\alpha}$$ be its length, $$|\boldsymbol{\alpha}|$$ its sum, and $$\boldsymbol{\alpha}_{\backslash i}$$ the multi-index created by removing $$\boldsymbol{\alpha}_i$$: nodes that define a basis of $$\mathcal{P}_n(\Delta^d)$$ can be indexed by $$\boldsymbol{\alpha}$$ such that $$\#\boldsymbol{\alpha} = d+1$$ and $$|\boldsymbol{\alpha}| = n$$.

• For a vector $$\boldsymbol{b}$$, let $$\boldsymbol{b}_{+i}$$ be the vector created by inserting a zero for the $$i$$-th coordinate.

The recursive definition of barycentric node coordinates, $$\boldsymbol{b}_{\boldsymbol{X}}(\boldsymbol{\alpha}) \in \mathbb{R}^{\#\boldsymbol{\alpha}}$$ has the base case

(1)$\boldsymbol{b}_{\boldsymbol{X}}(\boldsymbol{\alpha}) = (x_{|\boldsymbol{\alpha}|,\boldsymbol{\alpha}_0}, x_{|\boldsymbol{\alpha}|,\boldsymbol{\alpha}_1}), \quad \#\boldsymbol{\alpha} = 2,$

and the recursion

(2)$\boldsymbol{b}_{\boldsymbol{X}}(\boldsymbol{\alpha}) = \frac{\sum_i x_{|\boldsymbol{\alpha}|,|\boldsymbol{\alpha}_{\backslash i}|} \boldsymbol{b}_{\boldsymbol{X}}(\boldsymbol{\alpha}_{\backslash i})_{+i}} {\sum_i x_{|\boldsymbol{\alpha}|,|\boldsymbol{\alpha}_{\backslash i}|}}, \quad \#\boldsymbol{\alpha} > 2.$

The full set of nodes is

(3)$R^d_{\boldsymbol{X},n} = \{\boldsymbol{b}_{\boldsymbol{X}}(\boldsymbol{\alpha}): \#\boldsymbol{\alpha} = d+1, |\boldsymbol{\alpha}| = n\}.$
Parameters
• d (int) – The dimension of the simplex

• n (int) – The maximum degree of the polynomials

• family (optional) – The 1D node family used to define the coordinates in the barycentric $$d$$-simplex. The default family='lgl' corresponds to the shifted Lobatto-Gauss-Legendre nodes. See “The family optional argument” for using other node families.

• domain (str, optional) – The domain for the $$d$$-simplex where the returned coordinates will be defined. See “The domain optional argument” for the choices and their formal definitions.

Returns

The nodes $$R^d_{\boldsymbol{X},n}$$ defined in (3), as a 2D array with $$\binom{n+d}{d}$$ rows. If domain='barycentric', it has $$d+1$$ columns, otherwise it has $$d$$ columns.

Return type

ndarray

Examples

Nodes for $$\mathcal{P}^4(\Delta^2)$$ in barycentric coordinates:

>>> recursive_nodes(2, 4)
array([[0.        , 0.        , 1.        ],
[0.        , 0.17267316, 0.82732684],
[0.        , 0.5       , 0.5       ],
[0.        , 0.82732684, 0.17267316],
[0.        , 1.        , 0.        ],
[0.17267316, 0.        , 0.82732684],
[0.2221552 , 0.2221552 , 0.5556896 ],
[0.2221552 , 0.5556896 , 0.2221552 ],
[0.17267316, 0.82732684, 0.        ],
[0.5       , 0.        , 0.5       ],
[0.5556896 , 0.2221552 , 0.2221552 ],
[0.5       , 0.5       , 0.        ],
[0.82732684, 0.        , 0.17267316],
[0.82732684, 0.17267316, 0.        ],
[1.        , 0.        , 0.        ]])


The same nodes on the unit triangle:

>>> recursive_nodes(2, 4, domain='unit')
array([[0.        , 0.        ],
[0.        , 0.17267316],
[0.        , 0.5       ],
[0.        , 0.82732684],
[0.        , 1.        ],
[0.17267316, 0.        ],
[0.2221552 , 0.2221552 ],
[0.2221552 , 0.5556896 ],
[0.17267316, 0.82732684],
[0.5       , 0.        ],
[0.5556896 , 0.2221552 ],
[0.5       , 0.5       ],
[0.82732684, 0.        ],
[0.82732684, 0.17267316],
[1.        , 0.        ]])


If we construct the node set not from the Lobatto-Gauss-Legendre 1D node family, but from the equispaced 1D node family, we get equispaced 2D nodes:

>>> recursive_nodes(2, 4, family='equi', domain='unit')
array([[0.  , 0.  ],
[0.  , 0.25],
[0.  , 0.5 ],
[0.  , 0.75],
[0.  , 1.  ],
[0.25, 0.  ],
[0.25, 0.25],
[0.25, 0.5 ],
[0.25, 0.75],
[0.5 , 0.  ],
[0.5 , 0.25],
[0.5 , 0.5 ],
[0.75, 0.  ],
[0.75, 0.25],
[1.  , 0.  ]])


This is what they look like mapped to the equilateral triangle:

>>> import matplotlib.pyplot as plt
>>> from recursivenodes import recursive_nodes
>>> nodes_equi = recursive_nodes(2, 4, family='equi', domain='equilateral')
>>> nodes_lgl = recursive_nodes(2, 4, domain='equilateral')
>>> plt.scatter(nodes_lgl[:,0], nodes_lgl[:,1], marker='o', label='recursive LGL')
<matplotlib.collections.PathCollection object at ...>
>>> plt.scatter(nodes_equi[:,0], nodes_equi[:,1], marker='^', label='equispaced')
<matplotlib.collections.PathCollection object at ...>
>>> plt.gca().set_aspect('equal')
>>> plt.legend()
<matplotlib.legend.Legend object at ...>
>>> plt.show()


## Why does it matter?¶

Runge’s phenomenon: when constructing a polynomial interpolant $$p_f^n \in \mathcal{P}_n(\Delta^d)$$ that interpolates $$f$$ at given points $$X_n = \{\boldsymbol{x}_{n,i}\}$$ (“interpolation nodes”),

$p_f^n(\boldsymbol{x}_{n,i}) = f(\boldsymbol{x}_{n,i}), \quad i=1,\dots, \binom{n+d}{d},$

the locations of the $$\boldsymbol{x}_{n,i}$$’s has a big effect on how different $$p_f^n$$ is from $$f$$ on $$\Delta^d$$. Equispaced nodes can define bad interpolants, most famously for the Witch of Agnesi:

Here is the degree 15 equispaced polynomial interpolant:

Now here is the degree 15 interpolant constructed from recursive_nodes():

The most accepted way to characterize the interpolation quality of a set of nodes $$X_n$$ on the $$d$$-simplex is “The Lebesgue constant”, $$\Lambda_n^{\max}(X_n):\mathbb{R}^{\binom{n+d}{d}d} \to \mathbb{R}$$.

## Are recursive_nodes() the best?¶

No.

But in some cases they are the best that have yet been computed.

To find the best interpolation nodes requires minimizing $$\Lambda_n^{\max}(X_n)$$: it is nonconvex and there is no known explicit formula for nodes that minimize it.

### 2D¶

For the triangle, authors have optimized related functions [CB95][Hes98][Rot05], directly optimized $$\Lambda_n^{\max}(X_n)$$ over all potential nodes [Rot05][Hei05][RSV12], and directly optimized $$\Lambda_n^{\max}(X_n)$$ while preserving desirable properties such as symmetry [Hei05][RSV12].

In this figure, recursive_nodes() (blue) are on par with the best known node sets until about $$n = 10$$, and then perform worse. Both Roth [Rot05] (orange), and Rapetti, Sommariva, and Vianello [RSV12] (green) have calculated node sets that try to minimize the Lebesgue constant directly. That they disagree for some degrees is an indicator of the difficulty of finding the global minimum of $$\Lambda_n^{\max}(X_n)$$.

Rapetti, Sommariva, and Vianello in the same work also computed constrained minima of $$\Lambda_n^{\max}(X_n)$$, restricting the node sets to:

• node sets whose trace nodes on the edges are the LGL nodes (red), and

• node sets with full $$S_3$$ symmetry whose trace nodes on the edges are the LGL nodes (purple).

In contrast to node sets implicitly defined by a minimization problem are node families with explicitly defined construction rules. Equispaced and recursive_nodes() have already been seen: this package also implements the nodes of Blyth, Luo and Pozrikidis [BP06][LP06] and the widely used warp & blend nodes of Warburton [War06] (although these nodes do have a single tuning parameter that has been optimized separately for degrees up to 15).

Among explicitly defined node families, recursive_nodes() are usually within ~2% of Warburton’s nodes, which perform best. This is because, despite the different ways the node families are defined and computed, they are quite similar:

>>> from recursivenodes import recursive_nodes
>>> from recursivenodes.nodes import warburton
>>> import matplotlib.pyplot as plt
>>> nr = recursive_nodes(2, 15, domain='equilateral')
>>> nw = warburton(2, 15, domain='equilateral')
>>> ax = plt.gca()
>>> ax.scatter(nr[:,0], nr[:,1], marker='x', label='recursive_nodes')
<matplotlib.collections.PathCollection object at ...>
>>> ax.scatter(nw[:,0], nw[:,1], marker='+', label='warburton')
<matplotlib.collections.PathCollection object at ...>
>>> ax.set_aspect('equal')
>>> ax.legend()
<matplotlib.legend.Legend object at ...>
>>> plt.show()


### 3D¶

In 3D, directly optimizing $$\Lambda_n^{\max}(X_n)$$ becomes more challenging, because the size of $$X_n$$ grows cubically with $$n$$. It seems none of the authors who directly optimized for nodes for the triangle have yet followed up with tetrahedral nodes. Only approaches that optimize simpler, related functions have been extended to the tetrahedron [CB96][HT00].

In this figure, indirect optimization approaches are not noticeably better or worse than recursive_nodes() for the degrees that have been computed.

Among explicitly defined node families, recursive_nodes() are barely worse than the best for $$n < 7$$ and are best by an increasing margin for larger degrees.

## So why use recursive_nodes()?¶

def _recursive(d, n, alpha, family):
'''The barycentric d-simplex coordinates for a
multiindex alpha with sum n, based on a 1D node family.'''
xn = family[n]
b = np.zeros((d+1,))
if d == 1:
b[:] = xn[[alpha[0], alpha[1]]]
return b
weight = 0.
for i in range(d+1):
alpha_noti = alpha[:i] + alpha[i+1:]
n_noti = n - alpha[i]
w = xn[n_noti]
br = _recursive(d-1, n_noti, alpha_noti, family)
b[:i] += w * br[:i]
b[i+1:] += w * br[i:]
weight += w
b /= weight
return b

• They have full $$S_{d+1}$$ symmetry in every dimension $$d$$.

• The boundary traces of $$d$$-simplex recursive nodes are $$(d-1)$$-simplex recursive nodes, and in the base case the edges always match the 1D node set used in construction.

• When built from the Lobatto-Gauss-Legendre family, they perform about as well in 2D as other explicitly constructed node families, and much better in 3D.

• Until optimal nodes are computed for the tetrahedron, they have the best known Lebesgue constants in 3D for $$n \geq 7$$.

## How does it work?¶

Equispaced triangle nodes project onto 1D equispaced nodes on the edges of the triangle:

Blyth and Pozrikidis observed in [BP06] that Fekete points of the triangle, which for low degrees have good Lebesgue constants, nearly project onto Lobatto-Gauss-Legendre nodes on the edges of the triangle:

Some intriguing observations can be made regarding the location of some of the Fekete points in a given set. […] If an imaginary line is drawn through the nodes […] the two Fekete nodes sit on this line, close to the two zeros of the second Lobatto polynomial, Lo2, scaled by the length of the imaginary line.

If we try to reverse engineer a point that projects onto Lobatto-Gauss-Legendre nodes, the lines nearly meet, but not quite:

Going back to equispaced nodes, we find that the intersection point of the three projection lines can be written in barycentric coordinates relative to the projection points, and that these relative barycentric coordinates are themselves points from a 1D equispaced node set:

(Note that the relative barycentric coordinates have not been scaled to sum to 1.)

Using the same idea by analogy with Lobatto-Gauss-Legendre points, we have the recursive node definition:

Or, put into the form given in the definition:

## Where is there more detail?¶

A preprint has been submitted to arXiv. You can read a version that includes source code for all tables and figures here.