A geometric approach for the addition of nodes to an interpolatory quadrature rule with positive weights

TL;DR

This paper introduces a geometric framework for adding nodes to interpolatory quadrature rules while maintaining positive weights, enabling the design of new nested quadrature rules with improved properties.

Contribution

A novel geometric method is developed to determine how nodes can be added or replaced in quadrature rules without losing positive weights, including explicit inequalities for single node addition.

Findings

Explicit inequalities for single node addition are derived.

Adding multiple nodes while preserving positivity is always possible.

New nested quadrature rules outperform Gaussian and Clenshaw-Curtis rules.

Abstract

A novel mathematical framework is derived for the addition of nodes to univariate and interpolatory quadrature rules. The framework is based on the geometrical interpretation of the Vandermonde matrix describing the relation between the nodes and the weights and can be used to determine all nodes that can be added to an interpolatory quadrature rule with positive weights such that the positive weights are preserved. In the case of addition of a single node, the derived inequalities that describe the regions where nodes can be added are explicit. Besides addition of nodes these inequalities also yield an algorithmic description of the replacement and removal of nodes. It is shown that it is not always possible to add a single node while preserving positive weights. On the other hand, addition of multiple nodes and preservation of positive weights is always possible, although the minimum number of nodes that need to be added can be as large as the number of nodes of the quadrature rule. In case of addition of multiple nodes the inequalities describing the regions where nodes can be added become implicit. It is shown that the well-known Patterson extension of quadrature rules is a special case that forms the boundary of these regions and various examples of the applicability of the framework are discussed. By exploiting the framework, two new sets of quadrature rules are proposed. Their performance is compared with the well-known Gaussian and Clenshaw-Curtis quadrature rules, demonstrating the advantages of our proposed nested quadrature rules with positive weights and fine granularity.

Figures

Tables

Table 1. Table 1 : The test functions from Genz  [ 12 ] , which depend on the shape and translation parameters a 𝑎 a and b 𝑏 b .

Integrand Family	Attribute
$u_1\left(x\right)=\cos \left(2\pi b+ax\right)$	Oscillatory
$u_2\left(x\right)={\left(a^{-2}+{\left(x-b\right)}^2\right)}^{-1}$	Product Peak
$u_3\left(x\right)={\left(1+ax\right)}^{-2}$	Corner Peak
$u_4\left(x\right)=\exp \left(-a^2{\left(x-b\right)}^2\right)$	Gaussian
$u_5\left(x\right)=\exp \left(-a\left|x-b\right|\right)$	$C^0$ function
	Discontinuous

Equations

$$\mathcal{I}u={\int }_{\Omega }u(x)\rho (x)\text{d}x={\int }_a^{b}u(x)\rho (x)\text{d}x.$$

$${\mathcal{A}}_{N}u:=\sum _{k=0}^{N}u(x_k)w_k\approx \mathcal{I}u.$$

$${\mathcal{A}}_N\phi =\mathcal{I}\phi ,\text{ for all φ∈P(K)},$$

$$\underset{V(X_N)}{\underbrace{\left(\begin{array}{ccc}{x}_0^0& \cdots & x_N^0\\ ⋮& \ddots & ⋮\\ x_0^N& \cdots & x_N^N\end{array}\right)}}\left(\begin{array}{c}{w}_0\\ ⋮\\ w_N\end{array}\right)=\left(\begin{array}{c}{\mu }_0\\ ⋮\\ {\mu }_N\end{array}\right),$$

$${\mu }_k={\int }_{\Omega }{x}^k\rho (x)\text{d}x.$$

$$\det V(X_N)=\text{leavevmode}\prod _{0\le i<j\le N}\text{leavevmode}(x_j-x_i),$$

$$|{\mathcal{A}}_{N}u-\mathcal{I}u|$$

$$=(\Vert {\mathcal{A}}_N{\Vert }_{\infty }+{\mu }_0)\Vert u-{\phi }_N{\Vert }_{\infty }.$$

$$\Vert {\mathcal{A}}_N{\Vert }_{\infty }=\text{leavevmode}\underset{\Vert u{\Vert }_{\infty }=1}{\sup }\text{leavevmode}|{\mathcal{A}}_{N}u|=\sum _{k=0}^N|w_k|=\sum _{k=0}^N{w}_k={\mu }_0,$$

$$|{\mathcal{A}}_{N}u-\mathcal{I}u|\le 2{\mu }_0\text{leavevmode}\underset{\phi \in \mathbb{P}(N)}{\inf }\text{leavevmode}\Vert u-\phi {\Vert }_{\infty }.$$

$$x_k=\cos \left(\frac{k}{N}\pi \right),\text{ for k=0,…,N}.$$

$${\int }_{\Omega }{x}^j\left[\prod _{k=0}^{N+M}(x-x_k)\right]\rho (x)\text{d}x=0,\text{ for j=0,…,M−1}.$$

$$\mathbf{v}=\sum _{\begin{array}{c}k=0\\ k\ne k_0\end{array}}^N{b}_k{\mathbf{v}}_k.$$

$$\sum _{k=0}^N{c}_k{\mathbf{v}}_k=0.$$

$$\mathbf{v}=\sum _{k=0}^N(a_k-\alpha c_k){\mathbf{v}}_k.$$

$$\alpha =\min \left(\frac{a_k}{c_k}\mid c_k>0\right)=:\frac{a_{k_0}}{c_{k_0}}.$$

$$\mathbf{v}$$

$${\mu }_j=\sum _{\begin{array}{c}k=0\\ k\ne k_0\end{array}}^N{x}_k^j(w_k-\alpha c_k),\text{ for all j=0,…,N−1}.$$

$$\sum _{k=0}^{N+1}{x}_k^j{w}_k^{(N+1)}={\mu }_j,\text{ for j=0,…,N+1}.$$

$$W_N=\{w_0^{(N)},\dots ,w_N^{(N)}\}.$$

$${\mu }_j=\sum _{k=0}^N{x}_k^j{w}_k^{(N)}=\sum _{k=0}^{N+1}{x}_k^j{w}_k^{(N+1)},\text{ for j=0,…,N}.$$

$$0=\sum _{k=0}^{N+1}{x}_k^j{w}_k^{(N+1)}-\sum _{k=0}^{N+1}{x}_k^j{w}_k^{(N)}=\left(\sum _{k=0}^{N+1}{x}_k^j{w}_k^{(N)}+\sum _{k=0}^{N+1}{x}_k^j{c}_k\right)-\sum _{k=0}^{N+1}{x}_k^j{w}_k^{(N)}=\sum _{k=0}^{N+1}{x}_k^j{c}_k.$$

$$\sum _{k=0}^{N+1}{x}_k^{N+1}{w}_k^{(N+1)}={\mu }_{N+1},$$

$${\epsilon }_{N+1}:={\mu }_{N+1}-\sum _{k=0}^N{x}_k^{N+1}{w}_k^{(N)}=\sum _{k=0}^{N+1}{x}_k^{N+1}{c}_k.$$

$$\left(\begin{array}{cccc}{x}_0^0& \cdots & x_N^0& x_{N+1}^0\\ ⋮& \ddots & ⋮& ⋮\\ x_0^N& \cdots & x_N^N& x_{N+1}^N\\ x_0^{N+1}& \cdots & x_N^{N+1}& x_{N+1}^{N+1}\end{array}\right)\left(\begin{array}{c}{c}_0\\ ⋮\\ c_N\\ c_{N+1}\end{array}\right)=\left(\begin{array}{c}0\\ ⋮\\ 0\\ {\epsilon }_{N+1}\end{array}\right),$$

$$V(X_{N+1})\mathbf{c}=\mathbf{\epsilon },$$

$$c_k=\frac{\det V_k(X_{N+1})}{\det V(X_{N+1})},$$

$$\det V_k(X_{N+1})=(-1{)}^{(N+2)+(k+1)}{\epsilon }_{N+1}\det V(X_{N+1}\setminus \{x_k\})=(-1{)}^{N+k+1}{\epsilon }_{N+1}\det V(X_{N+1}\setminus \{x_k\}),$$

$$c_k=\frac{\det V_k(X_{N+1})}{\det V(X_{N+1})}$$

$$=(-1{)}^{N+k+1}{\epsilon }_{N+1}\left(\prod _{\begin{array}{c}0\le i<j\le N+1\\ i,j\ne k\end{array}}(x_j-x_i))/(\prod _{0\le i<j\le N+1}(x_j-x_i)\right)$$

$$={\epsilon }_{N+1}/(\prod _{\begin{array}{c}j=0\\ j\ne k\end{array}}^{N+1}(x_k-x_j)).$$

$$-{\epsilon }_{N+1}/(\prod _{\begin{array}{c}j=0\\ j\ne k\end{array}}^{N+1}(x_k-x_j))\le w_k^{(N)},\text{ for k=0,…,N+1}.$$

Full text

Keywords: Quadrature rules (65D32), Numerical integration (65D30), Interpolation (65D05)

1 Introduction

This article is concerned with the addition of nodes to univariate and interpolatory quadrature rules with positive weights. If such a quadrature rule is given, the goal is to determine all sequences of nodes such that, upon adding all nodes from such a sequence to the rule, an interpolatory quadrature rule with positive weights is again obtained. The motivation of this problem is twofold. Firstly, approximations of integrals computed using interpolatory quadrature rules with positive weights converge for any absolute continuous function [27, 7, 5]. Secondly, nested quadrature rules allow for straightforward refinements of the quadrature rule approximation, which is especially relevant if the integrand is computationally expensive.

Possibly the best-known interpolatory quadrature rule is the Gaussian quadrature rule [14], which exists for virtually any probability distribution with finite moments. It has positive weights and maximal polynomial degree. However, the nodes are not nested. The Gauss–Kronrod quadrature rule is an extension of a Gaussian quadrature rule, such that two nested rules with positive weights are obtained [20, 31]. The Gauss–Kronrod–Patterson quadrature rule [23, 24] further extends this idea by repeatedly applying the same algorithm, such that a sequence of nested rules is obtained. However, it does not exist for any distribution [17, 18]. Even though many other extensions have been proposed over the years [19, 21, 13], in general it is difficult to obtain a series of nested quadrature rules with positive weights. Moreover often the smallest possible granularity between two consecutive nested quadrature rules can only be found by exhaustive search [4].

An other large group of well-known quadrature rules is formed by the Clenshaw–Curtis quadrature rules [6], or simply those quadrature rules that are based on Chebyshev approximations (the Clenshaw–Curtis rule is formed by the Chebyshev extrema). Besides having excellent interpolation properties [16], it is well-known that these quadrature rules have positive weights if the distribution under consideration is uniform (explicit expressions are known [30]). Moreover for non-uniform distributions, the condition number of the quadrature rule converges to unity [5]. However, the vanilla Clenshaw–Curtis nodes are only nested for exponentially growing numbers of nodes [15].

Both the Gaussian and Clenshaw–Curtis quadrature rules have explicitly predefined nodes based on the roots of orthogonal polynomials. This results into accurate quadrature rules, but the construction of an accurate nested quadrature rule with fine granularity based on these rules remains notoriously difficult.

In this article the goal is to propose a geometrical framework for the addition of nodes to an interpolatory quadrature rule with positive weights and use this framework to determine all interpolatory quadrature rules with positive weights that extend a rule based on predefined nodes. It will be demonstrated rigorously that the boundary of the set that contains all nodes that can be added is equivalent to the Patterson extension of quadrature rules, such that a special case of the framework is an extension of the aforementioned Gaussian quadrature rule families.

The approach taken is based on the geometrical interpretation of the linear system describing the nodes and the weights [2, 26, 9], which yields a necessary and sufficient condition for a quadrature rule to have positive weights. The framework embeds previous results on the removal of nodes from quadrature rules [2, 8, 33] and describes, besides a geometrical description of all nodes that can be added to a quadrature rule, algorithms that can be used to construct and modify interpolatory quadrature rules with positive weights.

The addition and replacement of a single node can be determined analytically, whereas numerical methods are required to determine the bounds on the sets describing multiple nodes. The focus of this article is mainly on the geometrical and mathematical aspects and not on the numerical aspects of the proposed algorithms. However, to illustrate the potential of the framework, two straightforward examples of quadrature rules with positive weights that can be constructed by exploiting the proposed techniques are discussed.

In Section 2 the nomenclature used in this article is discussed, including the motivation behind enforcing positive weights. In Section 3 the problem of adding a single node to a quadrature rule is considered, which can be solved analytically. It is not always possible to add a node to a quadrature rule such that the resulting rule has positive weights. Therefore the theory is extended to adding multiple nodes in Section 4, where the results developed for adding a single node will be used extensively. It is always possible to add multiple nodes to a quadrature rule, provided that any number of nodes may be added to the rule. To demonstrate the advantages of nested quadrature rules with positive weights, two quadrature rules that are derived in this work are compared with the well-known Gaussian and Clenshaw–Curtis quadrature rule. The details and results of this numerical experiment are discussed in Section 5. Conclusions and suggestions for future work are discussed in Section 6.

2 Preliminaries

The quadrature rule nomenclature relevant for this article is discussed in Section 2.1. The relevance of positive weights and the relation between positive weights and accuracy of a quadrature rule is briefly reviewed in Section 2.2. The mathematical notion of adding nodes to a quadrature rule can be interpreted as a non-trivial extension of the removal of nodes [2, 8, 33], which is briefly discussed in Section 2.3. Finally, the problem setting of this article and the main results obtained from this work are summarized mathematically in Section 2.4.

2.1 Nomenclature

A quadrature rule is a well-known approach to approximate a weighted integral in the interval $\Omega=[a,b]\subset\mathbb{R}$ with $-\infty\leq a<b\leq\infty$. The weighting function is a positive density function $\rho\colon\Omega\to[0,\infty)$. The main interest is to approximate the integral over a given continuous function $u\colon\Omega\to\mathbb{R}$, i.e. to approximate the following operator:

[TABLE]

A quadrature rule approximates this integral by means of a weighted average, consisting of nodes and weights, which we denote by $X_{N}=\{x_{0},\dots,x_{N}\}\subset\Omega$ and $W_{N}=\{w_{0},\dots,w_{N}\}\subset\mathbb{R}$ respectively. The quadrature rule is the following operator $\mathcal{A}_{N}$:

[TABLE]

It is common to measure the consistency of this construction by means of polynomial degree. The polynomial degree of a quadrature rule is defined as the maximum polynomial degree the quadrature rule integrates exactly, or equivalently: a quadrature rule of degree $K$ has the property

[TABLE]

where $\mathbb{P}(K)$ denotes the space of all univariate polynomials of degree $K$ or less. This definition is only meaningful if $\rho$ has finite moments, so that is assumed to be the case throughout this article.

A quadrature rule is called interpolatory if the dimension of $\mathbb{P}(K)$ is larger than or equal to the number of nodes, or in other words, if $K\geq N$. Such quadrature rules can be formed by integrating the polynomial interpolant of $u$ using the nodes $X_{N}$. As the title of this article suggests, these quadrature rules are the main focus of this work: throughout this article the interest is mainly in rules with $K=N$ (though quadrature rules with $K>N$, such as the Gaussian quadrature rules, will also be considered).

The operators $\mathcal{A}_{N}$ and $\mathcal{I}$ and the space $\mathbb{P}(K)$ are linear, so if $K=N$, (2.3) defines a linear system that can be used to determine the weights, given the nodes and the moments of the distribution. Throughout this article a monomial basis of $\mathbb{P}(K)$ is considered. In this case, the matrix of the linear system is the well-known Vandermonde matrix, denoted as follows:

[TABLE]

with $\mu_{k}$ the raw moments of $\rho$:

[TABLE]

Throughout this article it is assumed that $\mu_{k}$ is known exactly for all $k$. The notation $V(X_{N})$ is used for the matrix of this linear system. It is well-known that

[TABLE]

such that, given the nodes, (2.4) defines a unique solution of the weights provided that all nodes are distinct.

2.2 Accuracy of quadrature rules

In this article the focus is on constructing interpolatory quadrature rules with non-negative weights (which we will call with a little abuse of nomenclature a positive quadrature rule). An approximation of an integral by means of such a quadrature rule converges if the integrand is sufficiently smooth [27], which can among others be demonstrated by applying the Lebesgue inequality [5], provided that $\Omega$ is bounded. To this end, let $u$ be given and let $\varphi_{N}$ be the best approximation polynomial [32] of degree $N$ of $u$, i.e. $\varphi_{N}=\operatorname*{arg\,min}_{\varphi\in\mathbb{P}(N)}\|\varphi-u\|_{\infty}$. Then

[TABLE]

Here, it holds that

[TABLE]

where it is used that $|w_{k}|=w_{k}$. Hence the following inequality is obtained:

[TABLE]

This shows many similarities with the classical Lebesgue inequality [16, 32] and demonstrates that if $u$ can be approximated well using a polynomial, it can be integrated using a quadrature rule with positive weights. Similar results exist for unbounded $\Omega$ [5, 34].

Two well-known interpolatory quadrature rules with positive weights are the Clenshaw–Curtis and Gaussian quadrature rules.

The Clenshaw–Curtis quadrature rule [6] has nodes $X_{N}$ that are defined as follows for $\Omega=[-1,1]$:

[TABLE]

The Clenshaw–Curtis quadrature rule has positive weights if the uniform distribution is considered and for any other distribution with bounded support the sum of the absolute weights becomes arbitrary close to $\mu_{0}$ for large $N$ [5]. The quadrature rule is nested for specific levels: it holds that $X_{N_{L}}\subset X_{N_{L+1}}$ with $N_{L}=2^{L}$ (for $L=1,2,\dots$).

The nodes of the Gaussian quadrature rule [14] are defined as the roots of the orthogonal polynomials with respect to the distribution $\rho$ under consideration, e.g. Legendre polynomials for the uniform distribution, Jacobi polynomials for the Beta distribution, etc. The uniquely defined rules always have positive weights and with $N+1$ nodes the rule has degree $2N+1$, however the rules are not nested.

The Gauss–Kronrod and Gauss–Patterson quadrature rules are extensions of Gaussian quadrature rules such that upon adding $M$ nodes (with $M=N+2$ for the Gauss–Kronrod rule) to a rule of $N+1$ nodes, a (not necessarily positive) rule of degree $N+2M$ is obtained [24, 20]. The Patterson extension is also applicable to non-Gaussian quadrature rules, though possibly complex-valued nodes can be obtained. The idea is to solve the following problem for $x_{N+1},\dots,x_{N+M}$, given quadrature rule nodes $X_{N}$:

[TABLE]

Then the obtained rule has degree $N+2M$ [5, Theorem 5.1.3], is defined uniquely, and possibly has complex-valued nodes. By construction, a Gaussian quadrature rule is obtained if $M=N+1$ (the weights of the nodes in $X_{N}$ become zero). These rules are reobtained as a special case in the framework discussed in this work.

2.3 Removal of nodes

The primary focus of this article is on the addition of nodes, but the obtained mathematical expressions can be interpreted as reversing the removal of nodes from an existing quadrature rule. Using Carathéodory’s theorem, it can be shown that for each positive interpolatory quadrature rule $X_{N}$, $W_{N}$ there exist two nodes $x_{k_{0}}$ and $x_{k_{1}}$ such that $X_{N}\setminus\{x_{k_{0}}\}$ and $X_{N}\setminus\{x_{k_{1}}\}$ both form the nodes of interpolatory quadrature rules with positive weights [2, 8, 33, 28]. The details are discussed in the constructive proof of the following theorem.

Theorem 1** (Carathéodory’s theorem).**

Let $\mathbf{v}_{0},\dots,\mathbf{v}_{N}$ be $N+1$ vectors spanning an $N$-dimensional space $V$. Let $\mathbf{v}\in V$ be such that $\mathbf{v}=\sum_{k=0}^{N}a_{k}\mathbf{v}_{k}$ with all $a_{k}\geq 0$. Then there exist non-negative $b_{k}$ and a $k_{0}\in\{0,\dots,N\}$ such that

[TABLE]

Proof.

The vectors $\mathbf{v}_{0},\dots,\mathbf{v}_{N}$ are linearly dependent, since these are $N+1$ vectors spanning an $N$-dimensional space. Hence there exists a vector $\mathbf{c}={(c_{0},\dots,c_{N})}^{\text{T}}\neq\mathbf{0}$ such that

[TABLE]

Hence for any $\alpha\in\mathbb{R}$, we have that

[TABLE]

In particular, consider the following $\alpha$ and $k_{0}$:

[TABLE]

With these choices it holds that $a_{k}-\alpha c_{k}\geq 0$ for all $k$ and $a_{k_{0}}-\alpha c_{k_{0}}=0$, concluding the proof as follows:

[TABLE]

The theorem can be used straightforwardly to remove nodes from a quadrature rule. To this end, let the positive interpolatory quadrature rule $X_{N}$ and $W_{N}$ be given. The goal is to construct an interpolatory quadrature rule using $N$ nodes from $X_{N}$ (which consists of $N+1$ nodes). Therefore, let $\mathbf{v}_{0},\dots,\mathbf{v}_{N}$ be the columns of the Vandermonde matrix of degree $N-1$, i.e. $\mathbf{v}_{k}={(x_{k}^{0},\dots,x_{k}^{N-1})}^{\text{T}}$. Then $\mathbf{v}_{k}$ are $N+1$ vectors spanning an $N$-dimensional space. The proof of Carathéodory’s theorem yields that there exists a vector $\mathbf{c}$, scalar $\alpha$, and index $k_{0}$ such that

[TABLE]

Moreover, we have that $w_{k_{0}}-\alpha c_{k_{0}}=0$, so by using $X_{N-1}=\{x_{k}\in X_{N}\mid k\neq k_{0}\}$ and $W_{N-1}=\{w_{k}-\alpha c_{k}\mid k\neq k_{0}\}$ a positive interpolatory quadrature rule is obtained. Notice that the vector $\mathbf{c}$ is computable, since it is a null vector of the Vandermonde matrix of degree $N-1$ (which is a $N\times(N+1)$-matrix).

This approach can be used to compute nested quadrature rules, but limits the accuracy of those quadrature rules to the initial rule of which nodes are removed. It is of less use if this rule is inadequately accurate or if no such rule is available. A possible approach to alleviate this is to use random samples as initial quadrature rule [3], though such samples do not accurately integrate higher order moments. The necessity of an existing quadrature rule is one of the main motivations to consider the addition of nodes, since that does not require the computation of an initial quadrature rule of sufficient accuracy.

2.4 Problem setting and main results

The problem studied in this article is how to add nodes to a positive interpolatory quadrature rule such that it remains positive and interpolatory. To formulate this mathematically, let a positive interpolatory quadrature rule $X_{N}$, $W_{N}$ be given. Then the goal is to determine, for given $M$, all nodes such that the set $X_{N+M}$ contains the nodes of a positive interpolatory quadrature rule and such that the rules are nested, i.e. $X_{N}\subset X_{N+M}$. To keep the nomenclature concise, we will refer to this problem as adding $M$ nodes to a positive interpolatory quadrature rule, where by “adding” we always mean addition such that the resulting quadrature rule has positive weights. Moreover the number of nodes added to a quadrature rule should be minimal, so we are also interested in the minimal value of $M$ (with $M>0$) such that a positive interpolatory quadrature rule with nodal set $X_{N+M}$ exists.

If a positive quadrature rule is given that is not interpolatory, i.e. a quadrature rule such that $\mathcal{A}_{N}\varphi=\mathcal{I}\varphi$ for all $\varphi\in\mathbb{P}(K)$ with $K<N$, a positive interpolatory quadrature rule can be deduced from this rule by repeatedly applying Theorem 1. Therefore we assume in this article without loss of generality that all quadrature rules are interpolatory.

The approach is to formulate, for given $M$, a necessary and sufficient condition for all $M$ nodes that can be added. This condition can be used firstly to determine whether such nodes exist for a specific $M$ and secondly to determine the nodes themselves. Moreover the derived theory allows for specific adjustments of quadrature rules. These adjustments consist of replacing and removing nodes from the quadrature rule, in such a way that the degree of the rule is not affected.

The analysis is split into two sections. The addition of a single node ($M=1$) can be solved analytically and is discussed in Section 3. The addition of multiple nodes ($M>1$) can only be done analytically for special cases. Based on the theory for $M=1$, this problem is analyzed in Section 4.

3 Addition of one node

Let $X_{N}$, $W_{N}$ be a positive interpolatory quadrature rule. The goal is to determine all $x_{N+1}$ such that $X_{N+1}=X_{N}\cup\{x_{N+1}\}$ forms the nodes of a positive interpolatory quadrature rule, i.e. there exists a set of non-negative weights $W_{N+1}$ such that

[TABLE]

Here, $w^{(N+1)}_{k}$ are the weights in the set $W_{N+1}$ and $\mu_{j}$ is assumed to be known. Notice that in general $W_{N}$ and $W_{N+1}$ will completely differ, so we use the following notation for any $N$:

[TABLE]

Moreover, with a little abuse of notation we will use $w^{(N)}_{k}=0$ for all $k>N$.

In Section 3.1 we derive a necessary and sufficient condition for such an $x_{N+1}$ to exist, which depends on the current nodes, weights, and moment $\mu_{N+1}$. As such, the developed theory provides practical adjustments of a quadrature rule. These constitute addition and replacement of a node, without reducing the degree of the interpolatory quadrature rule. The details are discussed in Section 3.2 and will be very useful in the remainder of this article. In Section 3.3 the Patterson extension is discussed in light of the derived adjustments and some basic applications of the derived procedures are discussed, including the construction of a (partially) nested quadrature rule with positive weights.

3.1 Positive weight criterion

The key notion is that if the node $x_{N+1}$ is given, a vector $\mathbf{c}={(c_{0},\dots,c_{N+1})}^{\text{T}}$ can be constructed such that $w^{(N+1)}_{k}=w^{(N)}_{k}+c_{k}$ (for $k=0,\dots,N+1$). This is the vector used in Section 2.3 to remove nodes from a rule. If this vector is such that $c_{k}\geq-w_{k}^{(N)}$, then $w_{k}^{(N+1)}\geq 0$, which is the primary goal. In this section, these properties are translated to conditions on $x_{N+1}$ that describe in which cases a node can be added to a quadrature rule.

The interpolatory quadrature rule $X_{N}$, $W_{N}$ has degree $N$, so after adding $x_{N+1}$ the following should hold to ensure that the new rule is interpolatory:

[TABLE]

From this, it follows for $j=0,\dots,N$ that (using $w_{N+1}^{(N)}=0$):

[TABLE]

The goal is to construct $X_{N+1}$ and $W_{N+1}$ such that they form a quadrature rule of degree $N+1$. Hence with $\mu_{N+1}=\int_{\Omega}x^{N+1}\,\rho(x)\,\text{d}x$ given, it should hold that

[TABLE]

which can be expressed in terms of the vector $\mathbf{c}$ as

[TABLE]

The value of $\varepsilon_{N+1}$ can be interpreted as the approximation error of the quadrature rule with nodes $X_{N}$ and weights $W_{N}$ with respect to $\mu_{N+1}$. Combining (3.4) and (3.6) yields the following system of linear equations for the vector $\mathbf{c}$:

[TABLE]

or more compactly:

[TABLE]

with ${\boldsymbol{\varepsilon}}={(0,\dots,0,\varepsilon_{N+1})}^{\text{T}}$. The vector ${\boldsymbol{\varepsilon}}$ has a large number of zeros so it is convenient to apply Cramer’s rule to this linear system, which yields

[TABLE]

where $V_{k}(X_{N+1})$ is equal to $V(X_{N+1})$ with the $k$-th column replaced by ${\boldsymbol{\varepsilon}}$, where the indexing of columns is started with 0. This expression can be simplified by noticing that

[TABLE]

with $V(X_{N+1}\setminus\{x_{k}\})$ the $(N+1)\times(N+1)$ Vandermonde matrix constructed with the nodal set $X_{N+1}\setminus\{x_{k}\}$. By using (2.6), the following is obtained for $k=0,\dots,N+1$:

[TABLE]

The denominator of this expression can be written as $\ell^{\prime}_{N}(x_{k})$, where $\ell_{N}(x)=\prod_{j=0}^{N}(x-x_{j})$ is the nodal polynomial. To keep the dependence on $x_{N+1}$ clear, this notation is used sparingly in this article.

The goal is to have positive weights, i.e. $w^{(N+1)}_{k}=w^{(N)}_{k}+c_{k}\geq 0$, which can be used to prove the following theorem.

Theorem 2**.**

Let $X_{N}$, $W_{N}$ form an interpolatory quadrature rule. Then $X_{N+1}=X_{N}\cup\{x_{N+1}\}$ forms the nodal set of a positive interpolatory quadrature rule if and only if

[TABLE]

Proof.

If $X_{N+1}$ forms the nodal set of a positive interpolatory quadrature rule, then

[TABLE]

Subtracting $w_{k}^{(N)}$ from both sides of the inequality yields (3.12). Vice versa, if (3.12) holds, it follows that

[TABLE]

If ${\varepsilon}_{N+1}=0$, i.e. $\mathcal{A}_{N}x^{N+1}=\mu_{N+1}$, then the theorem yields that the new rule has positive weights if and only if the current rule has positive weights. This is not surprising: any node $x_{N+1}$ can be added to such a rule with $w^{(N+1)}_{N+1}=0$ (and with $w^{(N+1)}_{k}=w^{(N)}_{k}$ for $k=0,\dots,N$).

From a computational point of view (3.12) might not be a numerically stable way of computing the bounds that describe all nodes that can be added. In the context of quadrature rules, numerical instabilities are usually alleviated by changing the basis of the Vandermonde matrix, but this is not applicable in this case since the determinant is up to a scaling factor independent from the basis used to construct the Vandermonde matrix (and this factor cancels out in (3.9)). Nonetheless, (3.12) can be evaluated in a numerical stable way using the well-known barycentric formulation of the interpolating polynomial. The interested reader is referred to [1].

3.2 Quadrature rule adjustments

Theorem 2 describes a necessary and sufficient condition for a quadrature rule to have positive weights if both $x_{N+1}$ and ${\varepsilon}_{N+1}$ are known. A main novelty of this work is to employ a geometrical interpretation of (3.12), from which several possible adjustments of quadrature rules can be derived. The most straightforward one is that all nodes $x_{N+1}$ can be determined that yield a positive interpolatory quadrature rule upon adding one of them to an existing quadrature rule. Moreover the formula also yields procedures to replace nodes in a quadrature rule, keeping the weights positive. The latter adjustment will be useful in Section 4, where it can be used to determine all possible $M$ nodes that can be added to a rule.

In Section 3.2.1 we further consider (3.12) and discuss the geometrical relation between the new node $x_{N+1}$ and the quadrature error ${\varepsilon}_{N+1}$. In Section 3.2.2 and 3.2.3 we discuss the addition and replacement of nodes in a positive interpolatory quadrature rule such that positivity of the weights is preserved. These operations follow directly from the geometrical interpretation of Theorem 2 derived in Section 3.2.1. The removal of a node, as outlined in Section 2.3, can also be formulated as a consequence of Theorem 2, which is not done here since the removal of nodes has been considered extensively in previous work [2, 8, 33, 28].

3.2.1 Geometry of nodal addition

The inequalities from (3.12) are $N+2$ linear inequalities in $x_{N+1}$ and ${\varepsilon}_{N+1}$. This can be seen by rewriting (3.11) as follows:

[TABLE]

If two values of $x_{N+1}$, $c_{k}$ (for $k=0,\dots,N+1$), or ${\varepsilon}_{N+1}$ are known, all other values can be determined from these expressions, which enforces that the obtained quadrature rule is again interpolatory. To incorporate positive weights, we use that for $k=0,\dots,N$ it holds that

[TABLE]

By combining this with (3.12) and requiring $w^{(N)}_{k}+c_{k}\geq 0$ inequalities of the following form are obtained:

[TABLE]

These are linear inequalities describing the relation between $x_{N+1}$ and ${\varepsilon}_{N+1}$ such that $w^{(N+1)}_{k}\geq 0$ for $k=0,\dots,N$. For $k=N+1$ it holds that $w_{k}^{(N)}=0$, so by using that $c_{k}=w_{k}^{(N+1)}$, (3.15) translates to:

[TABLE]

Even though the rightmost inequalities are non-linear, their sign solely depends on the location of $x_{N+1}$ with respect to the other nodes. Hence the exact value of the product is not of importance.

Example 1**.**

The inequalities from (3.18) are visualized as functions from $x_{N+1}$ to ${\varepsilon}_{N+1}$ in Figure 1 for the quadrature rule with $X_{N}$ and $W_{N}$ as follows:

[TABLE]

This is an (obviously positive) interpolatory quadrature rule with $\Omega=[-1,1]$ and $\rho\equiv 1/2$. The solid lines in the figures depict all $(x_{N+1},{\varepsilon}_{N+1})$ pairs such that one weight becomes equal to zero (i.e. where equality is attained in inequality (3.17) or (3.18)). The region where individual weights are positive are shaded in subfigures 1, 1, 1, and 1. Subfigure 1 is the intersection of these figures and therefore depicts regions where all weights are positive. Any $(x_{N+1},{\varepsilon}_{N+1})$ pair in the shaded region describes a positive interpolatory quadrature rule that contains the original three nodes.

The left subfigures demonstrate some key properties of the derived inequalities. The inequalities are linear and switch sign at the node, which is the rightmost condition of (3.17). The characteristics of the last inequality (subfigure 1) solely depend on the location of $x_{N+1}$ with respect to the other nodes. A combination of all inequalities (subfigure 1) has varying characteristics between different nodes, but it is always a system of linear inequalities. The line ${\varepsilon}_{N+1}=0$ is contained in all shaded regions, because any node with weight equal to zero can be added to the rule if the next moment $\mu_{N+1}$ is already correctly integrated by the quadrature rule.

The relation between $x_{N+1}$ and ${\varepsilon}_{N+1}$ from (3.18) can be interpreted in two ways. Firstly, if a new node $x_{N+1}$ is given, an upper bound and lower bound on ${\varepsilon}_{N+1}$ can be determined such that upon adding $x_{N+1}$ to the quadrature rule, a positive interpolatory quadrature rule is obtained. Geometrically these are the bounds of the shaded area with the $x=x_{N+1}$ line. This interval is never empty (as ${\varepsilon}_{N+1}=0$ is always in the shaded region). Secondly, if ${\varepsilon}_{N+1}$ is given, a (possibly empty) set can be determined such that a positive interpolatory quadrature rule is obtained upon adding a node from such a set. Geometrically this is equivalent to determining the bounds of the shaded area with the $y={\varepsilon}_{N+1}$ line.

The second interpretation can be used to add nodes to a quadrature rule, i.e. ${\varepsilon}_{N+1}$ is known and the goal is to determine $x_{N+1}$ (this is discussed in Section 3.2.2). The first interpretation can be used to replace nodes within a quadrature rule: $x_{N+1}$ is added to the node and an existing node can be removed by setting its weight to zero (this is discussed in Section 3.2.3).

3.2.2 Addition of a node

A direct consequence of (3.17) is that all nodes that can be added to a quadrature rule can be defined by means of intervals, obtained via a linear inequality. The results are discussed in the following lemmas. The first focuses on keeping the existing weights of the quadrature rule positive, the second focuses on ensuring that the additional weight (i.e. of the added node) is positive.

Lemma 3**.**

Let $X_{N}$, $W_{N}$ form the nodes and the weights of a positive interpolatory quadrature rule, let ${\varepsilon}_{N+1}$ from (3.6) be given, and let index $k$ of node $x_{k}$ be given. Let $x_{N+1}^{[k]}$ be as follows:

[TABLE]

Then $w_{k}^{(N+1)}\geq 0$ if and only if $x_{N+1}\in I_{k}$ with

[TABLE]

Or in other words, if and only if $x_{N+1}$ is not between $x_{k}$ and $x_{N+1}^{[k]}$.

Proof.

Adding a node is determining an $x_{N+1}$ that solves (3.17) if ${\varepsilon}_{N+1}$ is known. Hence, to keep the $k$-th weight positive, this is equivalent to computing the solution $x^{[k]}_{N+1}$ of the following problem:

[TABLE]

Here we used $\ell_{N}^{\prime}$ to make the notation more compact. Hence if $w_{k}^{(N)}\neq 0$:

[TABLE]

The node $x^{[k]}_{N+1}$ is such that, if added to the quadrature rule, an interpolatory quadrature rule is obtained with $w^{(N+1)}_{k}=0$ (the other weights may be negative). Assume $x_{k}<x_{N+1}^{[k]}$, without loss of generality. Then any node $x_{N+1}$ with $x_{N+1}\geq x_{N+1}^{[k]}$ or $x_{N+1}<x_{k}$ solves (3.17) for a single $k$. This is equivalent to stating that $x_{N+1}\in I_{k}=\mathbb{R}\setminus[x_{k},x_{N+1}^{[k]})$. ∎

The proof of this lemma can also be stated geometrically, using one of the Figures 1, 1, or 1. If ${\varepsilon}_{N+1}$ is known, those $x_{N+1}$ that are such that $(x_{N+1},{\varepsilon}_{N+1})$ is not part of a gray region form the interval as stated in the theorem. Here, $x_{N+1}^{[k]}$ is the intersection of the line passing through $x_{k}$ and the constant line ${\varepsilon}_{N+1}$. All intervals $I_{k}$ are bounded, so there always exists a node $x_{N+1}\in(I_{0}\cap\cdots\cap I_{N})$, or in other words, there always exists a node that keeps the existing $N+1$ weights of a quadrature rule positive upon addition.

Obviously, the goal is also to ensure that the weight of the added node is positive, which can be described by means of a series of intervals. The details of this are discussed in the following lemma.

Lemma 4**.**

Let $X_{N}$, $W_{N}$ form the nodes and the weights of a positive interpolatory quadrature rule, let ${\varepsilon}_{N+1}$ from (3.6) be given. Without loss of generality, assume that $x_{0}<x_{1}<\cdots<x_{N}$. Then $w_{N+1}^{(N+1)}\geq 0$ upon addition of $x_{N+1}$ to the quadrature rule if and only if one of the following holds for all $k=0,\dots,N$:

• •

$x_{N+1}\in[x_{k-1},x_{k}]$* if the signs of $\ell_{N}^{\prime}(x_{k})$ and ${\varepsilon}_{N+1}$ are equal (e.g. both are negative);*

• •

$x_{N+1}\in[x_{k},x_{k+1}]$* if the signs of $\ell_{N}^{\prime}(x_{k})$ and ${\varepsilon}_{N+1}$ differ.*

For $k=N$, use $x_{k+1}=\infty$ and for $k=0$, use $x_{k-1}=-\infty$ (with a little abuse of notation).

Proof.

Recall the derivation of (3.18), i.e. the relation between $x_{N+1}$, $w_{N+1}^{(N+1)}$, and ${\varepsilon}_{N+1}$:

[TABLE]

It holds that $w^{(N+1)}_{N+1}>0$ if $\prod_{j=0}^{N}(x_{N+1}-x_{j})$ and ${\varepsilon}_{N+1}$ have different sign. The first term flips sign only at $x_{N+1}=x_{k}$ (for any $k=0,\dots,N$), hence if, for given $k$,

[TABLE]

it is necessary that $x_{k-1}<x_{N+1}<x_{k}$ to ensure that $\prod_{j=0}^{N}(x_{N+1}-x_{j})$ is negative and $x_{k}<x_{N+1}<x_{k+1}$ to ensure that $\prod_{j=0}^{N}(x_{N+1}-x_{j})$ is positive. A similar result holds if

[TABLE]

Combining this with the sign of ${\varepsilon}_{N+1}$ results in the statement of the lemma. ∎

Geometrically, Lemma 4 describes the intervals of Figure 1. Notice that Lemma 4 can also straightforwardly be applied to cases where $w_{k}^{(N)}=0$ (for any $k=0,\dots,N$), i.e. if the quadrature rule has weights equal to zero.

Using Lemma 3 and Lemma 4 the set $I$ can be computed such that any $x_{N+1}\in I$ can be added to a quadrature rule $X_{N}$ and $W_{N}$ such that positive weights are obtained (and adding any $x_{N+1}\notin I$ yields a rule with at least one negative weight). The procedure is to firstly compute all intervals $I_{0},\dots,I_{N}$ from Lemma 3 and construct $I=I_{0}\cup\cdots\cup I_{N}$. Secondly, Lemma 4 is used to remove intervals of the form $[x_{k-1},x_{k}]$ from $I$.

The exact details of this procedure are outlined in Algorithm 1. No advanced interval arithmetic is necessary to implement this algorithm, only a procedure that implements the removal of an interval from a series of intervals is needed.

Example 2**.**

Reconsider the quadrature rule from Example 1. Then the bounds of the intervals containing nodes that can be added, i.e. the solutions of (3.22), are depicted in Figure 2 as open circles. Here, $\mu_{N+1}=0$, so from a straightforward computation it follows that ${\varepsilon}_{N+1}=-1/9$. A constant $\rho$ is considered here. In this case, the values of $x^{[k]}_{N+1}$ are (from left to right) $-5/3$, [math], and $7/9$, of which the first is not visible in the figure. Adding any of these nodes yields a quadrature rule with positive weights, but we emphasize that this is generally not the case for other quadrature rules. Hence adding any node from the set $I=(-\infty,-5/3]\cup[0,7/9]$ yields a positive interpolatory quadrature rule. Restricting $x_{N+1}$ to the set $\Omega$ further reduces the number of possible intervals.

Notice that $I=\emptyset$ if ${\varepsilon}_{N+1}\neq 0$ and $w_{k}^{(N)}=0$. This can be derived mathematically, but it also follows from the mere fact that all weights change (see (3.22)) upon addition of a node to a quadrature rule, so $w^{(N+1)}_{k}=w^{(N)}_{k}=0$ is not possible. If ${\varepsilon}_{N+1}=0$, no node can be added to enforce that $w_{k}^{(N)}=0$. However, any node with weight equal to zero can be added, hence the formula yields $x^{[k]}_{N+1}=x_{k}$ with $w^{(N+1)}_{N+1}=0$. Technically, the quadrature rule now has a node equal to $x_{k}$ with weight equal to zero. Nonetheless, this results into a singular Vandermonde matrix (which contradicts the theory developed so far), so we do not further study this specific case.

If $\Omega=\mathbb{R}$ and the number of nodes is odd, it is always possible to add a single node to a quadrature rule: in this case the result from Lemma 4 either states that $x_{N+1}\in(-\infty,x_{0}]$ or $x_{N+1}\in[x_{N},\infty)$, but never both. Geometrically this means that the leftmost and rightmost shaded region grow to infinity and minus infinity respectively (or vice versa). Similarly, if $\Omega=\mathbb{R}$ and the number of nodes is even, it is always possible to add a single node if ${\varepsilon}_{N+1}\geq 0$.

However, in any other case (i.e. that of a bounded $\Omega$ or even number of nodes with ${\varepsilon}_{N+1}<0$) adding a single node to a quadrature rule is not always possible, as shown in the following example.

Example 3**.**

Adding a single node to the following interpolatory quadrature rule is not possible when requiring positive weights:

[TABLE]

Note that this example can be obtained straightforwardly by adding the node $1/11$ to the quadrature rule from Example 1 and redetermining the weights likewise.

3.2.3 Replacement of a node

Replacing a node is equivalent to adding a node, with the difference that the goal is to determine this node such that the weight of an existing node in the obtained quadrature rule becomes zero, i.e. $w_{k}^{(N+1)}=0$ for a $k\leq N$. This is equivalent to determining a specific $(x_{N+1},{\varepsilon}_{N+1})$ pair that yields $w_{k}^{(N+1)}=0$, which was used to determine all possible additions in Section 3.2.2. The main difference with addition is that the next moment $\mu_{N+1}$ is not used, as the number of nodes and the degree of the rule do not change. This makes ${\varepsilon}_{N+1}$ a free variable.

The relation between ${\varepsilon}_{N+1}$ and $x_{N+1}$ is already derived, so by reconsidering (3.22) with the goal to determine both $x_{N+1}$ and all ${\varepsilon}^{[k]}_{N+1}$ (indexed by $[k]$ with $k=0,\dots,N$) that make $w_{k}^{(N+1)}=0$ the following expressions are obtained:

[TABLE]

We will interpret this expression as a function of $x_{N+1}$, denoted by ${\varepsilon}^{[k]}_{N+1}\colon\Omega\to\mathbb{R}$. By using ${\varepsilon}_{N+1}={\varepsilon}^{[k]}_{N+1}(x_{N+1})$, a positive interpolatory quadrature rule with $w^{(N+1)}_{k}=0$ is obtained upon adding $x_{N+1}$ to the rule.

It follows that for every $x_{N+1}\in\Omega$ there is an $x_{k}\in X_{N}$ such that the quadrature rule with nodes $(X_{N}\cup\{x_{N+1}\})\setminus\{x_{k}\}$ is positive and interpolatory. The details are discussed in the following lemma.

Lemma 5**.**

Let $X_{N}$, $W_{N}$ form the nodes and the weights of a positive interpolatory quadrature rule and let $x_{N+1}\in\Omega$ be given. Then there exists an $x_{k}$ such that $(X_{N}\cup\{x_{N+1}\})\setminus\{x_{k}\}$ forms the nodal set of a positive and interpolatory quadrature rule.

Proof.

Let ${\varepsilon}^{[k]}_{N+1}$ be defined by (3.28). Consider ${\varepsilon}_{-}$ and ${\varepsilon}_{+}$, defined as follows:

[TABLE]

Hence ${\varepsilon}_{-}<0<{\varepsilon}_{+}$. Using Lemma 3, it follows that using either ${\varepsilon}_{-}$ or ${\varepsilon}_{+}$ to add $x_{N+1}$ results in a quadrature rule with $w_{k}^{(N+1)}\geq 0$ for $k=0,\dots,N$. Moreover, by definition of ${\varepsilon}_{-}$ and ${\varepsilon}_{+}$ these rules have one (or more) weights equal to 0. From Lemma 4 it follows that either the rule constructed using ${\varepsilon}_{-}$ or ${\varepsilon}_{+}$ has $w_{N+1}^{(N+1)}\geq 0$ (and the other has $w_{N+1}^{(N+1)}\leq 0$).

Concluding, either ${\varepsilon}_{-}$ or ${\varepsilon}_{+}$ can be used to construct a positive interpolatory quadrature rule with at least one weight equal to zero. Nodes with weights equal to zero can be removed without affecting the quadrature rules. This is equivalent to having added a node $x_{N+1}$ and removed one, say $x_{k}$, which is the statement of the theorem. ∎

The proof of the lemma is constructive, and therefore describes a straightforward method to replace nodes in a quadrature rule. Given $x_{N+1}\in\Omega$, the procedure is to compute ${\varepsilon}_{+}$ and ${\varepsilon}_{-}$, figure out whether ${\varepsilon}_{N+1}={\varepsilon}_{+}$ or ${\varepsilon}_{N+1}={\varepsilon}_{-}$ yields $w^{(N+1)}_{N+1}\geq 0$ by using using Lemma 4, and finally compute the quadrature rule after replacement. These steps are outlined in detail in Algorithm 2. Geometrically, the approach computes the two lines closest to the ${\varepsilon}_{N+1}=0$ line, i.e. the boundary of the gray region, and determines which of these lines corresponds to obtaining a quadrature rule with only positive weights (see Figure 2).

Consequently, the domain of a quadrature rule, depicted by $\Omega\subset\mathbb{R}$, can be decomposed in subsets $\Omega_{0},\dots,\Omega_{N}$ that indicate which node can be replaced. If $x_{N+1}\in\Omega_{k}$, $(X_{N}\cup\{x_{N+1}\})\setminus\{x_{k}\}$ forms the nodes of a positive interpolatory quadrature rule. Combining the results of Lemma 4 and Lemma 5, these sets can be denoted in the following way:

[TABLE]

The sets $\Omega_{k}$ have been depicted in Figure 2. Notice that the boundaries of these sets correspond to positions where two lines intersect, or in other words, those $x_{N+1}\in\Omega_{k}$ that result into two weights equal to zero, if used for replacement. One of these weights is, by construction, $w_{k}^{(N+1)}$. If the other weight is $w_{l}^{(N+1)}$, we also have $x_{N+1}\in\Omega_{l}$. This geometrical observation can be made explicit, which can be used to actually compute $\Omega_{k}$: these $x_{N+1}$ have ${\varepsilon}_{N+1}^{[k]}(x_{N+1})={\varepsilon}_{N+1}^{[l]}(x_{N+1})$, or equivalently:

[TABLE]

Hence we have proved the following lemma.

Lemma 6**.**

Let $k$ be given and let $\partial\Omega_{k}$ denote the boundary of $\Omega_{k}$. Then, for any $x_{N+1}\in\partial\Omega_{k}$, we have that

[TABLE]

for an $l\in 0,\dots,N$.

The result is a procedure to compute the boundaries of a specific $\Omega_{k}$. Firstly, for $l=0,\dots,N$, compute $x_{(k,l)}$ such that

[TABLE]

Those $x_{(k,l)}$ that yield a positive interpolatory quadrature rule upon replacement (e.g. computed using Algorithm 2), form the boundary of the interval $\Omega_{k}$. If $x_{l}<x_{(k,l)}$, it follows that $[x_{l},x_{(k,l)})\notin\Omega_{k}$, since a replacement with $x_{N+1}\in[x_{l},x_{(k,l)}]$ results in a negative $w^{(N+1)}_{l}$ (similarly for $x_{l}>x_{(k,l)}$). The procedure to determine $\Omega_{k}$ explicitly is outlined in Algorithm 3. Here, the indexing is slightly changed to be able to reuse parts of Algorithm 1, since we still need to ensure that the weight of the added node (which replaces $x_{k}$) is positive.

Equation (3.34) does not necessarily have a solution for any $l$. Geometrically this is the case if the lines through $x_{k}$ and $x_{l}$ are parallel. In such a case, one should use $x_{(k,l)}\leftarrow\infty$ or $x_{(k,l)}\leftarrow-\infty$ in Algorithm 3, depending on the sign of the nominator when computing $x_{(k,l)}$ (usually, this happens automatically when using floating point arithmetic).

The values of $x_{N+1}$ that solve (3.34) form a special case. Since $x_{N+1}\in\Omega_{k}\cap\Omega_{l}$, the quadrature rule $(X_{N}\cup\{x_{N+1}\})\setminus\{x_{k},x_{l}\}$ is positive, interpolatory, and has degree $N$, even though it consists only of $N$ nodes. The latter result is remarkable: two nodes are removed and one is added, but the degree of the quadrature rule is not affected. Such rules have a non-trivial high degree and are therefore more accurate than interpolatory quadrature rules without this property.

Example 4**.**

An example of an interpolatory quadrature rule with non-trivial high degree is $X_{N}=\{-1,1/3\}$, obtained by adding $1/3$ to the quadrature rule of Example 1 (and removing all nodes with zero weight). All nodes that can be added to obtain such a rule are the intersection of two lines in Figure 2.

More generally, all nodes $x_{k}$ that can be added to a rule can be found by determining the bounds on the shaded region and observing which node belongs to the obtained bound. Consequently, the fact that $\Omega=\bigcup_{k=0}^{N}\Omega_{k}$ follows visually from Figure 2. Hence the relation between ${\varepsilon}^{[k]}_{N+1}$ and $x_{N+1}$, as described by (3.28), are the solid lines in Figure 2.

The node $x_{(k,l)}$ only depends on the nodes $x_{j}$ with $j\neq k$ and $j\neq l$, i.e. its value is independent from $x_{k}$ and $x_{l}$. This is not evident, as (3.34) depends on these nodes. However, it can be demonstrated by using that the rule interpolatory, which yields:

[TABLE]

Here, $L_{k}(x)$ is the $k$-th Lagrange basis polynomial. Replacing this expression in (3.34) and using that $\ell^{\prime}_{N}(x_{k})=\prod_{j\neq k}(x_{k}-x_{j})$ yields an equality that can be simplified to the following:

[TABLE]

This expression is in fact a Patterson extension (consider (2.12) with $j=0$). The tight relation between the Patterson extension and the framework discussed in this article is further discussed in Section 3.3.1.

3.3 Constructing quadrature rules

In the previous section the theoretical foundation for extending a positive interpolatory quadrature rule with a single node is derived. In this section, firstly it is discussed how addition relates naturally to the Patterson extension [24, 25] of (non-Gaussian) quadrature rules. Secondly, due to the simplicity of addition and replacement of a node, quadrature rules based on these procedures can be derived numerically fast and accurately, and an example is discussed.

As discussed previously, there does not always exist a single node that can be added such that positive weights are obtained, so it is non-trivial to construct a sequence of positive interpolatory quadrature rules by consecutively adding a single node to the rule. There are various possibilities to alleviate this, e.g. by allowing negative weights, relaxing the strict requirement that all nodes of the quadrature rule have to be preserved, or by adding multiple nodes instead of one. In this article, the second and third options are further considered. For this purpose, a quadrature rule is presented based on the replacement of nodes. The rule has positive weights and is interpolatory, but is strictly speaking not fully nested. The details are considered in Section 3.3.2. The addition of multiple nodes is further discussed in Section 4.

3.3.1 Patterson extension

Remarkably, both the addition and replacement of a node can yield a Patterson extension of a quadrature rule. In both cases, the focus is on the nodes that yield a zero weight upon addition to the quadrature rule.

In Section 3.2.2 it was noticed that any weight from a quadrature rule can be made equal to zero by exploiting the relation between ${\varepsilon}_{N+1}$ and $x_{N+1}$. In Example 2 the quadrature rule $X_{N}=\{-1,-1/6,1\}$ was considered, where the nodes $-5/3$, [math], and $7/9$ are such that upon adding one of these to the rule, a rule of only three nodes with non-zero weights of degree three is obtained. Notice that these nodes are Patterson extensions of quadrature rules (as discussed in Section 2.2), as they can be interpreted as adding one node ($M=1$) to a quadrature rule of two nodes ($N=1$), obtaining a rule of degree three ($N+2M=3$). This also holds in general: for given $k$, adding one node $x^{[k]}_{N+1}$ from (3.22) (so $M=1$) to the interpolatory quadrature rule $X_{N}\setminus\{x_{k}\}$ (with degree $N-1$) yields a quadrature rule with $N+1$ nodes and degree $N+1$ (which equals $(N-1)+2M$).

In Section 3.2.3 the notation $x_{(k,l)}$ was introduced to denote nodes that, upon adding them to the rule, yield a (possibly negative) interpolatory quadrature rule with $w^{(N+1)}_{k}=w^{(N+1)}_{l}=0$. These nodes also form a Patterson extension. To see this, notice that the replacement is adding a single node to the quadrature rule $X_{N-2}=X_{N}\setminus\{x_{k},x_{l}\}$. The Patterson extension of a single node of this quadrature rule is a quadrature rule consisting of $N$ nodes of degree $(N-2)+2M=N$ (adding one node means $M=1$). By construction, this rule has the nodes $X_{N-2}\cup\{x_{(k,l)}\}$.

Example 5**.**

Reconsider for example the quadrature rule with the nodes $X_{N}=\{-1,-1/6,1\}$ and $\rho\equiv 1/2$. Then it is straightforward to determine using (3.34) that $x_{(0,1)}=-1/3$, $x_{(0,2)}=2$, and $x_{(1,2)}=1/3$. Hence these are three Patterson extensions of the quadrature rule nodes $\{1\}$, $\{-1/6\}$, and $\{-1\}$. Indeed, the quadrature rules with the nodes $\{-1/3,1\}$, $\{-1/6,2\}$, or $\{-1,1/3\}$ have degree equal to 2.

Notice that $x_{(k,l)}$ is not a Patterson extension of the quadrature rule that has been used to determine it, i.e. $X_{N}$, $W_{N}$ in (3.34). However, its definition allows for a straightforward way to determine this extension. First, add (randomly) two nodes to the quadrature rule $X_{N}$, $W_{N}$, obtaining a possibly negative interpolatory quadrature rule $X_{N+2}$, $W_{N+2}$. Then the node $x_{(N+1,N+2)}$ is the Patterson extension of the quadrature rule with nodes $X_{N}$, because upon adding this node to $X_{N+2}$, the weights of the randomly added nodes become zero. As the Patterson extension is unique, this construction is well-defined. Naturally, this is not the preferred approach to construct a Patterson extension, but it embeds such extensions into the framework discussed here.

The Patterson extension is also obtained as a special case if multiple nodes are added to a quadrature rule. This will be discussed in Section 4.3.

3.3.2 Partially nested, positive, and interpolatory quadrature rule

The addition and replacement of a single node are straightforward procedures described as the solutions of linear inequalities. However, there does not always exist a single node that can be added such that all weights remain positive. In this section, this is alleviated by relaxing the requirement that $X_{N}\subset X_{N+M}$.

To this end, let $X_{N}$ and $\hat{X}_{N+1}$ be the nodes of two positive interpolatory quadrature rules, possibly with $X_{N}\not\subset\hat{X}_{N+1}$. The nodes $\hat{X}_{N+1}$ can for example form a Gaussian quadrature rule. The idea is to iteratively replace nodes in $\hat{X}_{N+1}$ with nodes from $X_{N}$, i.e. removing $x_{k}\in\hat{X}_{N+1}$ and adding $x_{k}\in X_{N}$. Ideally, all nodes $x_{k}\in X_{N}$ can be added to $x_{k}\in\hat{X}_{N+1}$, which would yield a rule that reuses all nodes in $X_{N}$.

In other words, if $X_{N}=\{x_{0},\dots,x_{N}\}$ and $\hat{X}_{N+1}=\{\hat{x}_{0},\dots,\hat{x}_{N+1}\}$, for each node $x_{k}\in X_{N}\setminus\hat{X}_{N+1}$ the set $\Omega_{j}$ is identified (see Section 3.2.3) such that $(\hat{X}_{N+1}\cup\{x_{k}\})\setminus\{\hat{x}_{j}\}$ is the nodal set of a positive and interpolatory quadrature rule. If there is an $x_{k}$ such that $\hat{x}_{j}\not\in X_{N}$, we set $\hat{X}_{N+1}\leftarrow(\hat{X}_{N+1}\cup\{x_{k}\})\setminus\{\hat{x}_{j}\}$ and keep repeating this procedure until no such $x_{k}$ exists anymore. If there are multiple $x_{k}$ that could possibly be used to trigger a replacement in $\hat{X}_{N+1}$, the smallest one is selected in the example presented in this article.

The nodes from $X_{N}$ that cannot be added to $\hat{X}_{N+1}$ are reconsidered in consecutive iterations and added again if possible. It is difficult to theoretically quantify the number of nodes from $X_{N}$ that can be “added” this way to $\hat{X}_{N+1}$, though it is straightforward to see that there exists at least a single $x_{k}\in X_{N}$ that can be reused.

To demonstrate this procedure numerically, let $X_{1}$ and $W_{1}$ form a Gaussian quadrature rule of two nodes. If the uniform distribution is considered, evaluating all quadrature rules up to $N=19$ requires in total $22$ unique evaluations of $u$, which is two more than optimally possible considering the limitations of the framework as discussed in this work. The obtained sequence is depicted in Figure 3 (the two additional evaluations of $u$ can be found at $N=15$ and $N=18$). This result seems to be somewhat independent from the distribution, since applying the same approach to construct a sequence of quadrature rules with respect to a $\operatorname{Beta}(10,10)$ distribution requires in total $23$ function evaluations, which is three more than optimally possible (the obtained rules are depicted in Figure 3).

The main advantage of this approach compared to the previously discussed Patterson extension is that it always has positive weights. Moreover the expressions to compute the nodes contained in the quadrature rule are straightforward. However, the approach has the same disadvantage as the removal of nodes (see Section 2.3), since it requires a sequence of existing quadrature rules.

4 Addition of multiple nodes

In the previous section a counterexample of a positive interpolatory quadrature rule is discussed that can not be extended by adding a single node. In this section we will therefore study the addition of multiple nodes to a quadrature rule. The problem setting is that of Section 2.4: given a positive interpolatory quadrature rule $X_{N}$, $W_{N}$, determine $M$ as small as possible and nodes $X_{N+M}$ with $X_{N}\subset X_{N+M}$ such that $X_{N+M}$ forms the nodes of a positive interpolatory quadrature rule.

The first step is to extend the derivation of Section 3.1 for the addition of multiple nodes. The derivation is again based on Cramer’s rule. With the theory that is derived in the upcoming Section 4.1 it is not obvious how nodes can be added to the quadrature rule, but it provides geometrical insight in the location of such nodes with respect to the existing nodes. Again we can derive some non-trivial adjustments one can apply to a quadrature rule. These are discussed in Section 4.2. Similar to the case of a single node, there is a tight relation with the Patterson extension. In this case, the Patterson extension for general $M$ is recovered. This is discussed in Section 4.3, including some examples of nested quadrature rules obtained with the theory derived in this section.

4.1 Positive weight criterion

The idea is similar to the derivation of the addition of single node. Let $X_{N}$ be the initial nodal set and let $M$ be given. The goal is to determine $X_{N+M}$ with $X_{N}\subset X_{N+M}$ such that it forms the nodal set of a positive interpolatory quadrature rule.

Let $w^{(N)}_{k}$ for $k=0,\dots,N$ be the weights of $W_{N}$ and likewise let $w^{(N+M)}_{k}$ be the (unknown) weights of $W_{N+M}$. Then there exists a vector $\mathbf{c}={(c_{0},\dots,c_{N},c_{N+1},\dots,c_{N+M})}^{\text{T}}$ such that $w_{k}^{(N+M)}=w_{k}^{(N)}+c_{k}$. The goal is to construct $\mathbf{c}$ such that the obtained rule is interpolatory and positive.

With a similar reasoning as before it is straightforward to observe that the following should hold for such a vector to ensure that the obtained quadrature rule is interpolatory:

[TABLE]

and

[TABLE]

where $\varepsilon_{j}$ is as previously introduced, i.e. $\varepsilon_{j}\coloneqq\mu_{j}-\sum_{k=0}^{N}x_{k}^{j}w^{(N)}_{k}$. This can be written in the form of a linear system as follows:

[TABLE]

Applying Cramer’s rule to this system requires more bookkeeping, as the right hand side contains multiple non-zero entries. Let ${\boldsymbol{\varepsilon}}={(0,\dots,0,\varepsilon_{N+1},\dots,\varepsilon_{N+M})}^{\text{T}}$, then Cramer’s rule prescribes

[TABLE]

where $V_{k}(X_{N+M})$ is equal to $V(X_{N+M})$ with the $k$-th column (indexed from 0) replaced by ${\boldsymbol{\varepsilon}}$. The numerator can be further expanded as follows:

[TABLE]

where $V_{(j,k)}(X_{N+M})$ is the $(j,k)$-minor of $V(X_{N+M})$ (i.e. the matrix without its $j$-th row and $k$-th column, where both indices start at 0). Hence for $c_{k}$ the following expression is obtained:

[TABLE]

The same derivation is commonly used to derive the determinant of a Vandermonde matrix [11, 22], and it is well-known that the ratio of determinants obtained in this expression is an elementary symmetric polynomial. The $k$-th elementary symmetric polynomial is generally defined as the sum of all monomial permutations of length $k$, that is as follows:

[TABLE]

The elementary symmetric polynomials are only defined for $k\leq N+1$ and by convention $e_{0}\equiv 1$. Concluding, the following expression is obtained for $c_{k}$:

[TABLE]

Here, $e_{k}$ is the $k$-th elementary symmetric polynomial as defined above. With a little abuse of notation, we used:

[TABLE]

We are now in a position to formulate a theorem in similar form as Theorem 2, but then for multiple nodes. The proof is omitted, since it is equivalent to that of Theorem 2, but then with the equalities derived in this section.

Theorem 7**.**

Let $X_{N}$, $W_{N}$ form an interpolatory quadrature rule. Then $X_{N+M}=X_{N}\cup\{x_{N+1},\dots,x_{N+M}\}$ forms the nodal set of a positive interpolatory quadrature rule if and only if

[TABLE]

For $M=1$, we have that the summation only incorporates $j=N+1$, hence $e_{N+M-j}(X_{N+M}\setminus\{x_{k}\})=e_{0}(X_{N+M}\setminus\{x_{k}\})=1$, recovering Theorem 2. So Theorem 7 is indeed a strict generalization of Theorem 2.

4.2 Quadrature rule adjustments

Theorem 7 presents a necessary and sufficient condition for a quadrature rule extended with $M$ nodes to have positive weights. Contrary to the addition of a single node, it cannot be used directly to determine possible nodes that can be added to the quadrature rule. This can be seen by rewriting it in a similar form as (3.17), i.e. for $k=0,\dots,N+M$:

[TABLE]

Notice that, if $x_{N+1},\dots,x_{N+M}$ are unknowns, an $M$-variate system of $N+M+1$ polynomial inequalities is obtained. In general these systems are very difficult to solve, so we do not directly pursue a solution of the system above. Nonetheless, the system still provides a geometrical interpretation about where solutions reside, similar to the case of single node addition (though less intuitive). This is discussed in Section 4.2.1. Based on these geometrical insights, procedures to replace nodes and to add nodes, which extend those explained previously, can be derived. These procedures are discussed in Section 4.2.2 and 4.2.3 respectively.

4.2.1 Geometry of nodal addition

The type of the inequalities (4.13) (i.e. “greater than” versus “less than”) does not change between two nodes and if this type is fixed, the system consists of polynomial inequalities. Hence the region where $M$ nodes can be added is described by a continuous boundary, bounded by the polynomial inequalities of (4.13), consisting of lines, surfaces, or “hypersurfaces” through the nodes.

If one of the right hand sides of (4.13) changes sign, there is an addition of $M$ nodes such that the inequality forms an equality for a specific $k$. In such cases, there is an addition such that one of the nodes obtains a weight equal to zero. This is equivalent to the case discussed in Section 3.2.3, where a single node is added in order to set the weights of another node equal to zero.

It is difficult to visualize the addition of $M$ nodes in a similar way as we visualized the addition of one node, as there are $M$ nodes $x_{N+1},\dots,x_{N+M}$ and $M$ quadrature rule errors ${\varepsilon}_{N+1},\dots,{\varepsilon}_{N+M}$. Plotting the errors with respect to the nodes (as in Figure 2) is therefore not viable, as this is a plot from $\mathbb{R}^{M}$ to $\mathbb{R}^{M}$.

On the other hand, if the distribution $\rho(x)$ is fixed beforehand, the values of ${\varepsilon}_{N+1},\dots,{\varepsilon}_{N+M}$ are known and contour plots of the regions encompassing all $M$ nodes that can be added can be made (provided that $M$ is small enough).

Example 6**.**

Let $\rho\equiv 1/2$ with $\Omega=[-1,1]$ and reconsider the quadrature rule from Example 1. In Figure 4 lines are depicted where the inequalities from (4.13) are equalities. The shaded area depicts regions where all inequalities are valid, i.e. any coordinate $(x_{N+1},x_{N+2})$ in the shaded region can be added to the respective quadrature rule in order to obtain a positive interpolatory rule. The figure is obviously symmetric around $x_{N+1}=x_{N+2}$, as the order of addition (i.e. first adding $x_{N+1}$ and then $x_{N+2}$ or vice versa) yields equivalent quadrature rules. Selecting a coordinate $(x_{N+1},x_{N+2})$ on one of the boundaries results into one weight equal to zero. Adding the coordinates on the corners, depicted by the open circles (i.e. “the boundary of the boundary”), results into two weights equal to zero.

The dashed lines indicate where the inequalities (4.13) with $k=N+1$ and $k=N+2$ change sign. If this happens, one of the new nodes $x_{N+1}$ or $x_{N+2}$ has weight equal to zero. This line forms everywhere a boundary of the shaded area: the node with weight equal to zero can be replaced by any other node, while still resulting into an interpolatory quadrature rule with positive weights. This situation is equivalent to adding a single node $x_{N+1}$ to the quadrature rule, but gaining two degrees, as discussed in Section 3.2.2.

The addition and replacement of multiple nodes follow readily from this example. Notice that if any coordinate $(x_{N+1},\dots,x_{N+M})$ is known, the replacement for $M=1$ can be used to reach any other coordinate $(x_{N+1},\dots,x_{N+M})$ in the same region (shaded in Figure 4). Hence if all corners of those regions are determined (depicted as open circles in Figure 4), the full region can be explored straightforwardly using Algorithm 3. As these corner cases form a replacement of nodes, we start by discussing replacement of $M$ nodes. Moreover, it will be shown that these corners are a Patterson extension. Based on the algorithm to determine all these corners, addition of $M$ nodes follows straightforwardly.

4.2.2 Replacement of multiple nodes

Let $X_{N}$, $W_{N}$ be an interpolatory quadrature rule and let indices $k_{1},\dots,k_{M}$ be given such that $0\leq k_{i}\leq N$ and $k_{i}\neq k_{j}$ for $i\neq j$. In this section the goal is to determine the interpolatory quadrature rule $X_{N+M}$, $W_{N+M}$ such that $w^{(N+M)}_{k_{i}}=0$ for all $k_{i}$. Notice that this is equivalent to replacing the nodes $x_{k_{1}},\dots,x_{k_{M}}$ in the quadrature rule $X_{N}$ by the nodes $x_{N+1},\dots,x_{N+M}$. The nodes with this property are the intersections of the polynomials of (4.13) and they are depicted as open circles in Figure 4. Moreover, they describe the boundary of the set of nodes that can be added to the quadrature rule.

The desired nodes $x_{N+1},\dots,x_{N+M}$ can be determined by calculating the Patterson extension of the interpolatory quadrature rule with the nodes $X_{N}\setminus\{x_{k_{1}},\dots,x_{k_{M}}\}$, for which efficient techniques exist [25, 20, 19]. Such techniques require that $M$ must be known a priori and they do not provide a simple geometrical interpretation. Therefore we proceed by embedding the Patterson extension in the framework discussed here. This yields an alternative, new algorithm to determine these nodes, which is mainly of theoretical and geometrical interest, since it requires the computation of large numbers of roots of polynomials.

We start by solving a slightly easier problem. Assume ${\varepsilon}_{N+1}=\cdots={\varepsilon}_{N+M-1}=0$ and ${\varepsilon}_{N+M}\neq 0$. Notice that, if ${\varepsilon}_{N+M}$ is neglected, any addition of $M-1$ nodes yields a valid quadrature rule (as these nodes have zero weight). Geometrically, a fully shaded figure (if drawn as Figure 4) is obtained. This can be exploited to determine the desired nodes, as only the value of ${\varepsilon}_{N+M}$ imposes a condition on the nodes $x_{N+1},\dots,x_{N+M}$.

The nodes that yield $w^{(N+M)}_{k_{1}}=\dots=w^{(N+M)}_{k_{M}}=0$ can be found by applying Theorem 7 with $c_{k_{i}}=-w^{(N)}_{k_{i}}$ for all $i$ or by consecutively applying Theorem 2. In both cases, the following is obtained:

[TABLE]

In principle this system of polynomial equalities is difficult to solve, but it has a certain structure that can be exploited. To see this, let $\hat{\ell}_{M}(x)$ be the nodal polynomial of the nodes $x_{N+1},\dots,x_{N+M}$:

[TABLE]

which translates the system above to

[TABLE]

If the nodal polynomial $\hat{\ell}_{M}$ is known, its roots equal $x_{N+1},\dots,x_{N+M}$. The nodal polynomial has degree $M$ and it is known that its leading order coefficient equals 1. Therefore it is useful to introduce the polynomial $q_{M}(x)\coloneqq\hat{\ell}_{M}(x)-x^{M}$, which has degree $M-1$. Then (4.16) can be rewritten as follows:

[TABLE]

These are $M$ values of a polynomial of degree $M-1$, which is a well-known interpolation problem and can be solved with various well-known methods (such as barycentric interpolation [1]). If $q_{M}$ is determined, the roots of the polynomial $\hat{\ell}_{M}(x)=q_{M}(x)+x^{M}$ are the nodes $x_{N+1},\dots,x_{N+M}$. By construction these nodes are such that $w^{(N+M)}_{k_{i}}=0$ for $i=1,\dots,M$.

Even though assuming ${\varepsilon}_{N+1}=\cdots={\varepsilon}_{N+M-1}=0$ is not realistic in practical cases, this procedure can readily be extended to the general case. For this we reuse the replacement step from Section 3.2.3. If ${\varepsilon}_{N+1}\neq 0$, then a single node is added to the quadrature rule such that $w^{(N+1)}_{k_{1}}=0$. This is equivalent to applying Algorithm 2 with $x_{N+1}=x_{(k,l)}$, as discussed in Section 3.2.3. Then the obtained quadrature rule $X_{N+1}\setminus\{x_{k_{1}}\}$ has ${\varepsilon}_{N+1}=0$. By applying the procedure discussed above to these $N+1$ nodes, the nodes $x_{N+2}$ and $x_{N+3}$ can be determined such that $w^{(N+2)}_{k_{2}}=0$ and $w^{(N+2)}_{N+1}=0$, i.e. we enforce that the weight of $x_{k_{2}}$ is zero and the weight of the previously added node becomes zero. The obtained rule has $N+3$ nodes, where two nodes have weight equal to zero. This is again a replacement, but here two nodes get weight equal to zero, which is a generalization of the replacement discussed in Section 3.2.3. Those nodes are removed to reobtain a quadrature rule of $N+1$ nodes and this process is repeated iteratively until $X_{N+M}$ is obtained. The obtained rule can be interpreted as a replacement of $M$ nodes, and yields the open circles from Figure 4. It is an iterative description: a replacement of $M$ nodes is determined using a replacement of $M-1$ nodes. Geometrically, we iterate over the dimension of the figure and iteratively determine a set of nodes that can be used as a replacement.

The obtained nodes form by definition a Patterson extension of the nodal set $X_{N}\setminus\{x_{k_{1}},\dots,x_{k_{M}}\}$, since it holds that $(X_{N}\setminus\{x_{k_{1}},\dots,x_{k_{M}}\})\cup\{x_{N+1},\dots,x_{N+M}\}$ has degree $N+M$. The existence of such a Patterson extension is directly coupled to the existence of $M$ nodes that can possibly be added to $X_{N}$ in the hope of obtaining an interpolatory quadrature rule with positive weights: if $M$ nodes can be added to the quadrature rule, the Patterson extension has positive weights, since it forms the boundary of the set that describes all additions. Moreover, if all Patterson extensions of all sets $X_{N}\setminus\{x_{k_{1}},\dots,x_{k_{M}}\}$ for any sequences $(k_{1},\dots,k_{M})$ have negative weights or are not real-valued, no addition of $M$ nodes exists.

Hence we have proved the following lemma.

Lemma 8**.**

Let $X_{N}$, $W_{N}$ form a positive interpolatory quadrature rule, let $\rho$ (or a sequence of moments) be the density function, and let $M$ be given. Then the following statements are equivalent:

There exists a Patterson extension of $M$ nodes of the quadrature rule $X_{N}$, $W_{M}$ with solely non-negative weights; 2. 2.

There exist $M$ nodes $x_{N+1},\dots,x_{N+M}$ such that $X_{N}\cup\{x_{N+1},\dots,x_{N+M}\}$ forms the nodal set of a positive interpolatory quadrature rule.

As stated before, any algorithm that computes Patterson extensions can be used to verify whether $M$ nodes exist that can be added to the rule. If a Patterson extension with non-negative weights is found, say $x_{N+1},\dots,x_{N+M}$, Algorithm 3 can be used to explore all possible additions to the quadrature rule.

The algorithm based on the geometrical interpretation used in this article is outlined in Algorithm 4. By iterating over all possible sorted sequences $(k_{1},\dots,k_{M})$, this procedure can be used straightforwardly to verify whether there exist $M$ nodes that can be added to a given quadrature rule (though this is a costly procedure).

There are two special cases that are (for sake of simplicity) not incorporated in Algorithm 4. Firstly, if $w^{(N+m)}_{k}=0$ at the start of an iteration, the polynomial $\hat{\ell}_{M}(x)$ is not well-defined. This can be incorporated by selecting any non-zero $w^{(N+m)}_{k_{i}}$ at the start of the iteration. If no such $w^{(N+m)}_{k_{i}}$ exists, then all these weights are zero, which is the primary goal of the algorithm. Secondly, if $r_{k}\in X_{N}$ or ${\varepsilon}_{N+m}=0$, a quadrature rule is obtained that has higher degree than its number of nodes. This can be incorporated by combining all double nodes in $X_{N}$ and likewise adding the respective weights and by skipping any iteration that has ${\varepsilon}_{N+m}=0$.

4.2.3 Addition of multiple nodes

By combining the quadrature rule replacement of Section 3.2.3 (for $M=1$) and the replacement of the previous section (for $M>1$), we obtained a naive algorithm to firstly determine $M$ as small as possible such that there exists a positive interpolatory quadrature rule $X_{N+M}$ (Algorithm 4) and secondly to explore all such $M$ nodes (Algorithm 3, yielding the shaded areas of Figure 4).

Determining the number of nodes $M$ that can be added to an interpolatory quadrature rule can straightforwardly be done by solving (4.14) for each sequence of $k_{1},\dots,k_{M}$ with $k_{1}<\cdots<k_{M}$. This gives all locations where $M$ nodes have zero weight. If at any of these locations all nodes have non-negative weight, then $M$ nodes can be added to the rule. Otherwise, $M$ is increased and the process is repeated.

Often the value of $M$ is unknown a priori. Besides determining the $M$ nodes that can be added, the goal is also to determine $M$ as small as possible (this is also how we formulated the problem originally in Section 2.4). Algorithm 4 can be used to determine $M$, as results from previous iterations can be reused. To see this, suppose a quadrature rule is given and by applying Algorithm 4 it is known that no addition of at most $M-1$ nodes exist. Then during these calculations, all sequences of nodes have been determined that make $M-1$ weights zero. By initializing Algorithm 4 with these sequences, only the last iteration of the loop is necessary, which significantly reduces the computational expense.

It is required to repeatedly determine large numbers of polynomial roots in this algorithm. This is nearly impossible to do symbolically, except for some special cases (e.g. $M\leq 3$ or symmetric quadrature rules). Moreover determining the roots numerically can result in quick aggregation of numerical errors. We use variable precision arithmetic, i.e. determine the roots with a large number of significant digits.

For large $N$ this is a costly algorithm, as the number of sorted sequences of length $M$ equals

[TABLE]

which grows fast for large $N$. Therefore using this algorithm to compute all removals is slower than using existing techniques to compute the Patterson extension, albeit that it is able to reuse all additions of $M-1$ nodes to compute all additions of $M$ nodes.

If all sets of $M$ nodes have been determined that can be added to the quadrature rule, the techniques from Section 3.2.3 can be used to fully explore all nodes that can be added to the rule. This requires solving linear equalities, which can be done fast and accurately.

The possibility of adding $M$ nodes to the quadrature rule does not guarantee the possibility of adding $M+1$ nodes to the quadrature rule.

Example 7**.**

We revisit the quadrature rule example from Example 1, i.e.

[TABLE]

In Figure 5 regions are depicted where a single node can be added (similar to Figure 2) and regions where, upon adding a node from that region, another node can be added (this is the projection of Figure 4). The addition of the rightmost node with the latter property is depicted in Figure 5, demonstrating that there is a single node that can be added and that this is indeed a limiting case.

Notice that the intervals where a single node and where two nodes can be added are independent from each other. There exist pairs of nodes $x_{N+1},x_{N+2}$ firstly such that both $W_{N+1}$ and $W_{N+2}$ are all positive (in the right interval surrounded by squares), secondly such that $W_{N+1}$ is positive, but $W_{N+2}$ is not (the right interval surrounded by circles, outside the interval surrounded by squares), thirdly such that $W_{N+1}$ is not positive, but $W_{N+2}$ is (the left interval surrounded by squares), and finally such that both $W_{N+1}$ and $W_{N+2}$ are always negative (outside all intervals).

4.3 Constructing quadrature rules

Similar to the case of addition of a single node, the Patterson extension is obtained for specific choices of nodes that are added to the rule. In fact, the nodes determined with Algorithm 4 are a Patterson extension of a quadrature rule with a smaller number of nodes. As the Gaussian quadrature rule is a special case of the Patterson extension, this rule also follows from the framework discussed in this article. This is discussed in more detail in Section 4.3.1.

By repeatedly applying Algorithm 4, a sequence of nested quadrature rules can be determined. These rules and their properties are considered in Section 4.3.2.

4.3.1 Patterson extension

The boundary of the set that describes all possible additions is spanned by the Patterson extension (the open circles in Figure 2 and Figure 4). These nodes have the property that, upon adding them to the quadrature rule, a rule of degree $N+M$ is obtained with $M$ weights equal to zero. This is equivalent to the Patterson extension of the quadrature rule without those $M$ nodes with zero weight. For $M=1$, this was demonstrated in Section 3.3.1.

For general $M$, the Patterson extension can be deduced mathematically as follows. Let $X_{N}$, $W_{N}$ be a quadrature rule and, as before, let $x_{N+1},\dots,x_{N+M}$ be such that the following nodes form a quadrature rule of degree $N+M$:

[TABLE]

Furthermore, let $X_{N-M}$ be the nodes of an interpolatory quadrature rule of degree $N-M$ be as follows:

[TABLE]

Upon adding $\{x_{N+1},\dots,x_{N+M}\}$ to $X_{N-M}$, the nodes from (4.20) are obtained, that have degree $N+M$. Hence $M$ nodes are added to an interpolatory rule of degree $N-M$ and the obtained degree is $N+M$, which is by definition a Patterson extension. Notice that the obtained quadrature rule is interpolatory, but not necessarily positive.

The Gaussian quadrature rule can be deduced as special case from Algorithm 4. To see this, suppose $M=N+1$, which is the number of nodes of the rule under consideration. In that case, there is only a single sequence of $k_{1},\dots,k_{M}$, defined as follows up to a permutation:

[TABLE]

By applying Algorithm 4, the nodes from (4.20) are obtained with $M=N+1$, which are:

[TABLE]

Hence the $N+1$ nodes $x_{N+1},\dots,x_{2N+1}$ form a quadrature rule of degree $2N+1$, which is by definition the Gaussian quadrature rule. In other words, when adding a Gaussian quadrature rule to an existing quadrature rule and setting all existing weights to zero a valid addition is obtained.

Example 8**.**

To demonstrate where Patterson extensions occur in our work, reconsider the interpolatory quadrature rule with the nodes $X_{N}=\{-1,-1/6,1\}$. In Section 3.3.1 three different Patterson extensions related to this quadrature rule were discussed: $\{-1/3,1\}$, $\{-1/6,2\}$, or $\{-1,1/3\}$. All these rules are Patterson extensions (of smaller quadrature rules) with $M=1$. To obtain a Patterson extension with $M=2$ and subsequently a Gaussian quadrature rule, consider Algorithm 4 using $\{k_{1},k_{2},k_{3}\}=\{0,1,2\}$. The algorithm proceeds as follows:

In the first iteration, it follows that $\hat{\ell}_{1}(x)=x+5/3$ and therefore the following quadrature rule is obtained:

[TABLE]

Notice that the node $x_{N+1}=-5/3$ was obtained in Section 3.2.2, where we discussed that after adding this node one obtains $w_{0}^{(3)}=0$. 2. 2.

In the second iteration, it follows that $\hat{\ell}_{2}(x)=x^{2}+2/5x-1/5$. Here, the Patterson extension with $M=2$ of the quadrature rule with “nodes” $\{1\}$ is obtained. Hence the following rule is obtained (notice that the node $-5/3$ is removed):

[TABLE] 3. 3.

In the third iteration, it follows that $\hat{\ell}_{3}(x)=x^{3}-3/5x$, whose roots are the Gaussian quadrature rule or, equivalently, the Patterson extension with $M=3$ of the empty quadrature rule:

[TABLE]

In this specific example it is possible to determine all nodes symbolically, but for larger values of $M$ this is generally not possible.

Considering the nodes in a different order results into different intermediate Patterson extensions, but obviously the Gaussian quadrature rule is the rule that is finally obtained. These steps also demonstrate the possibility to store intermediate results: only the nodes of step 2 are necessary to deduce the nodes of step 3.

Specialized algorithms exist for specific distributions and specific values of $N$ and $M$ to construct Gaussian, Gauss–Kronrod, and Gauss–Patterson quadrature rules [14, 20], but it remains a challenging topic to determine the Patterson extension for general non-Gaussian quadrature rules. The algorithm presented in this article is not an alternative for these existing algorithms, but embeds the Patterson extension in the discussed framework and can be used to determine all $M$ nodes that can be added to a quadrature rule. If an efficient procedure to determine large numbers of Patterson extensions is available, it can be readily used to determine whether an extension for a specific $M$ exists. By consecutively replacing the new nodes (see Section 3.2.3) all $M$ nodes that can be added can be found.

4.3.2 Nested, positive, and interpolatory quadrature rule

Algorithm 4 provides a straightforward procedure to determine the minimal value of $M$ and the positive interpolatory quadrature rule nodes $X_{N+M}$ such that $X_{N}\subset X_{N+M}$. The replacement procedure for $M=1$ of Section 3.2.3 can be used to determine all possible nodes, given $M$. This is the original goal of the article as outlined in Section 2.4 and examples of such quadrature rules are depicted in Figure 6. Here, each quadrature rule is iteratively extended with a minimal number of nodes, and the nodes that are added are selected randomly from the set containing all $M$ nodes that can be added. There are two main differences with the quadrature rules obtained in Section 3.3.2, where an existing rule was used as basis for a larger quadrature rule: the rules obtained in this section are fully nested, but do add more than one node between two consecutive rules.

Both figures demonstrate that $M$ varies significantly and does not increase monotonically. This is in line with the conclusions drawn in the Section 4.2.3, as shown in Figure 5. Moreover for almost all $N$, the value of $M$ is significantly larger in case the Beta distribution is considered, which is related to the “bad” initial set of nodes for this distribution. A different initialization would lead to different values of $M$.

5 Numerical integration with positive quadrature rules

This article is concerned with the construction of quadrature rules with positive weights and two new quadrature rules have been introduced: one based on the consecutive replacement of single nodes (possibly resulting in a sequence of rules that is not nested) and one by randomly adding nodes ensuring positive weights. We briefly assess the numerical performance of these quadrature rules by means of the Genz test functions (see Table 1). The Genz test functions [12] are functions defined on $\Omega=[0,1]$ constructed specifically to test integration routines. Each function has a specific family attribute that is considered to be challenging for integration routines, that can be enlarged by a shape parameter $a$ and translated by a translation parameter $b$. We restrict ourselves to the uniform distribution, as in this case the exact value of the integral of the Genz functions is known analytically.

We consider the performance of the following four quadrature rules:

A quadrature rule that is determined by consecutively adding and replacing nodes originating from a Gaussian quadrature rule (see Figure 3). This rule was discussed in Section 3.3 and is a partially nested, positive, and interpolatory quadrature rule. The rule is initialized with the quadrature rule nodes $X_{N}=\{0,5/12,1\}$ (i.e. the nodes from the example as discussed before, translated to $[0,1]$). 2. 2.

A quadrature rule that is determined by consecutively randomly adding $M$ nodes to the rule such that the obtained rule is positive. Here $M$ is minimal, i.e. the smallest number of nodes is added for each $N$ (see Figure 6). This rule was discussed in Section 4.3 and is a nested, positive, and interpolatory quadrature rule. The rule is initialized in the same way as the quadrature rule of the previous point, i.e. using $X_{N}=\{0,5/12,1\}$. 3. 3.

The Clenshaw–Curtis quadrature rule [6], where the nodes $X_{N}$ are defined explicitly by (2.11). It is well known that these nodes have positive weights if the distribution under consideration is uniform, which is the case. This positive and interpolatory quadrature rule is nested for specific levels, i.e. $X_{N_{L}}\subset X_{N_{L+1}}$ with $N_{L}=2^{L}$ ($l=1,2,\dots$). 4. 4.

The Gaussian quadrature rule [14], where the nodes and weights are defined as the quadrature rule with $N+1$ nodes of degree $2N+1$. This quadrature rule is not nested, so refining the quadrature rule results in a significant number of new function evaluations.

The error measure $e_{N}$ is the absolute integration error, i.e.

[TABLE]

where $u=u_{g}$ with $g=1,\dots,6$, i.e. $u$ is one of the Genz test functions. To obtain meaningful results we select the parameters $a$ and $b$ randomly in the unit interval and repeat the experiment 100 times. This also affects the reduced quadrature rule: each experiment selects the node that is removed randomly and therefore 100 different sequences of nested quadrature rules are obtained. The errors reported here are averaged over the 100 experiments and are therefore denoted by $\overline{e}_{N}$.

It is instructive to compare the error with the upper bound that follows from the Lebesgue inequality (2.10):

[TABLE]

where we use that $\mu_{0}=1$ in our test cases. This error is determined using the algorithm of Remez [32, Chapter 3], with the implementation from chebfun [10]. Convergence results for the uniform distribution $\rho\equiv 1$ in $\Omega=[0,1]$ are gathered in Figure 7.

Notice that regardless of the function under consideration all quadrature rule errors remain far under the dashed line, that represents the right-hand side of (5.2). This shows that the bound from this inequality is far from sharp.

The first four Genz functions can be approximated well using polynomials, as they are analytic and have rapidly converging Chebyshev coefficients. The best approximation converges exponentially in these cases, which is also the case for the four quadrature rules under consideration. The quadrature rules determined using the framework of this article perform slightly worse than the Clenshaw–Curtis and the Gaussian quadrature rule. This is related to the fact that these rules exploit the structure of the underlying distribution to a large extent (e.g. symmetry and higher-order moments), whereas the rules in this work only optimize for the positivity of the weights. The Gaussian quadrature rule converges with the highest rate, which is related to its high polynomial degree (a rule of $N+1$ nodes has degree $2N+1$). However, the Gaussian rule is not nested, so to refine the estimate of the integral for increasing number of nodes the number of function evaluations increases significantly. If a computationally expensive function is considered, using a nested quadrature rule with fine granularity (such as the proposed rules) significantly reduces the cost of refining the quadrature rule estimate.

The fifth Genz test function is not differentiable and can therefore not be approximated well using a polynomial. This can be observed from the best approximation polynomial, that converges with order 1 (so we would expect that $\overline{e}_{N}\sim 1/N$). In this case the difference between the Gaussian rule and the other rules is significantly smaller, demonstrating that the high polynomial degree of Gaussian rules is less relevant if the integrand is not smooth.

The sixth Genz test function cannot be approximated accurately using a polynomial when considering the $\infty$-norm, as it is discontinuous. Hence the best approximation error remains constant. However, the approximation of the quadrature rules still converges with order $1/2$. In this case, there is a clear difference between the integration error (that is an averaged error) and the best approximation error (that is a uniform error).

6 Conclusion

In this article, a novel mathematical framework is presented for the construction of nested, positive, and interpolatory quadrature rules by using a geometrical interpretation. Given an existing quadrature rule, necessary and sufficient conditions have been derived for $M$ new nodes to form an interpolatory quadrature rule with positive weights. The conditions have been formulated as inequalities, which are explicit if $M=1$ and implicit if $M>1$.

The addition of a single node can be treated as a special case, which can be solved analytically. The analytical expression can be used to add nodes to and replace nodes within a quadrature rule. The addition of multiple nodes can be determined numerically and a naive algorithm is presented for this purpose. Based on the quadrature rules obtained by this algorithm, the set that encompasses all additions of $M$ nodes can be explored by iteratively replacing nodes.

The well-known Patterson extension of quadrature rules forms a special case of the framework, as it is obtained by constructing the quadrature rules with $M$ weights equal to zero. As such, our proposed framework and its geometrical interpretation are well embedded in existing theory on the addition of nodes to quadrature rules. The framework provides various possibilities to construct or adapt quadrature rules and two examples have been discussed: one based on consecutively adding and replacing one node and one based on consecutively adding multiple nodes.

Numerical integration using the two quadrature rules introduced in this work shows the key advantages of nested quadrature rules with positive weights: estimates computed using the quadrature rules are stable and nesting allows for computationally cheap refinements of the estimates. Existing quadrature rules, such as the Gaussian and the Clenshaw–Curtis quadrature rule, are not nested with the fine granularity as the rules in this work.

There are various options to further extend the framework set out in this article. The algorithm to determine whether multiple nodes exist that can be added to the quadrature rule depends on determining many polynomial roots and iterates over all possible sequences of nodes that can become zero. For a large number of nodes this is computationally very costly and therefore warrants the need to derive an efficient algorithm to determine these nodes. Moreover the framework set out in this article does not use the relations that exists between consecutive moments of a distribution [29], which can possibly be used to further extend the framework set out in this text.

Acknowledgments

This research is part of the Dutch EUROS program, which is supported by NWO domain Applied and Engineering Sciences and partly funded by the Dutch Ministry of Economic Affairs.