Is the Féjer-Clenshaw-Curtis quadrature actually nested?

Methodology
1

Hi all,

I was playing with FejerAlgorithm to illustrate classical quadrature rules and I noticed an unexpected behavior. According to different sources (e.g., Chapter 9 from Sullivan, T.J. Introduction to uncertainty quantification. Vol. 63. Springer, 2015.), these quadratures should be nested, unlike Gaussian quadratures. Meaning, for example, that a rule of size 2n should contain the nodes from a rule of size n.
My understanding is that different versions of this rule were proposed, affecting the formulas of their weights. However the nodes should follow Chebyshev nodes (i.e., Chebyshev points of the first kind). The expression of these nodes on the interval [-1, 1] vary depending on the sources:

  • Expression 1: on the OT documentation and in (Trefethen, L.N. “Is Gauss quadrature better than Clenshaw–Curtis?.” SIAM review 50.1 (2008)): x_{k}=\cos \left({\frac {k\pi}{n}} \right)
  • Expression 2: on their Wikipedia page: x_{k}=\cos \left({\frac {2k-1}{2n}}\pi \right)

When using the FejerAlgorithm OpenTURNS class, I do not manage to recover the nested property.
Here is a MRE using the expression 1 and the FejerAlgorithm class, maybe I misunderstood the constructor:

import numpy as np
import matplotlib.pyplot as plt
import openturns as ot

# Expression 1
n1 = 16
n2 = 32
cheb_nodes_n1 = np.cos((np.arange(n1)) * np.pi / (n1))
cheb_nodes_n2 = np.cos((np.arange(n2)) * np.pi / (n2))
# FejerAlgorithm
ot_fejer_nodes_n1 = ot.FejerAlgorithm([n1]).getNodes()
ot_fejer_nodes_n2 = ot.FejerAlgorithm([n2]).getNodes()
#
plt.figure(figsize=(7, 2))
plt.scatter(cheb_nodes_n1, [0] * n1, marker=".", color="C2")
plt.scatter(cheb_nodes_n2, [0.02] * n2, marker=".", color="C2")
plt.scatter(ot_fejer_nodes_n1, [0.1] * n1, marker=".", color="C3")
plt.scatter(ot_fejer_nodes_n2, [0.12] * n2, marker=".", color="C3")
plt.xlim((-1.1, 1.05))
plt.ylim((-0.02, 0.14))
plt.text(1.06, 0.0, "Expression 1")
plt.text(1.06, 0.1, "FejerAlgorithm")
plt.grid(True);

Thanks for your feedback!
Elias

2

Hi Elias,

The FejerAlgorithm class provides 3 different quadrature rules: Clenshaw-Curtis, Fejer1 and Fejer2, as described here:

By default, FejerAlgorithm generates the Clenshaw Curtis rule, with nodes equal to x_k=\cos\left(\frac{k\pi}{N-1}\right) for k=0,\dots,N-1. To have nested quadratures you have to choose N and 2N-1 nodes. It is illustrated by the following code:

import numpy as np
import matplotlib.pyplot as plt
import openturns as ot

# Expression 1
n1 = 16
n2 = 2*n1-1
cheb_nodes_n1 = np.cos((np.arange(n1)) * np.pi / (n1-1))
cheb_nodes_n2 = np.cos((np.arange(n2)) * np.pi / (n2-1))
# FejerAlgorithm
ot_fejer_nodes_n1 = ot.FejerAlgorithm([n1]).getNodes()
ot_fejer_nodes_n2 = ot.FejerAlgorithm([n2]).getNodes()
#
plt.figure(figsize=(7, 2))
plt.scatter(cheb_nodes_n1, [0] * n1, marker=".", color="C2")
plt.scatter(cheb_nodes_n2, [0.02] * n2, marker=".", color="C2")
plt.scatter(ot_fejer_nodes_n1, [0.1] * n1, marker=".", color="C3")
plt.scatter(ot_fejer_nodes_n2, [0.12] * n2, marker=".", color="C3")
plt.xlim((-1.1, 1.05))
plt.ylim((-0.02, 0.14))
plt.text(1.06, 0.0, "Expression 1")
plt.text(1.06, 0.1, "FejerAlgorithm")
plt.grid(True);
plt.show()

ClenshawCurtis
ClenshawCurtis700×200 9.26 KB

The two other rules provided by the FejerAlgorithm class are ot.FejerAlgorithm.FEJERTYPE1: x_k=\cos\left(\frac{(2k+1)\pi}{2N}\right) for k=0,\dots,N-1 and ot.FejerAlgorithm.FEJERTYPE2: x_k=\cos\left(\frac{(k+1)\pi}{N+1}\right) for k=0,\dots,N-1. The ot.FejerAlgorithm.FEJERTYPE1 rules are not nested, but the ot.FejerAlgorithm.FEJERTYPE2 rules are. You have to choose N and 2N+1 nodes.
It is illustrated by the following code:

import numpy as np
import matplotlib.pyplot as plt
import openturns as ot

# Expression 1
n1 = 16
n2 = 2*n1+1
cheb_nodes_n1 = np.cos((np.arange(n1)+1) * np.pi / (n1+1))
cheb_nodes_n2 = np.cos((np.arange(n2)+1) * np.pi / (n2+1))
# FejerAlgorithm
ot_fejer_nodes_n1 = ot.FejerAlgorithm([n1], ot.FejerAlgorithm.FEJERTYPE2).getNodes()
ot_fejer_nodes_n2 = ot.FejerAlgorithm([n2], ot.FejerAlgorithm.FEJERTYPE2).getNodes()
#
plt.figure(figsize=(7, 2))
plt.scatter(cheb_nodes_n1, [0] * n1, marker=".", color="C2")
plt.scatter(cheb_nodes_n2, [0.02] * n2, marker=".", color="C2")
plt.scatter(ot_fejer_nodes_n1, [0.1] * n1, marker=".", color="C3")
plt.scatter(ot_fejer_nodes_n2, [0.12] * n2, marker=".", color="C3")
plt.xlim((-1.1, 1.05))
plt.ylim((-0.02, 0.14))
plt.text(1.06, 0.0, "Expression 1")
plt.text(1.06, 0.1, "FejerAlgorithm")
plt.grid(True);
plt.show()

Fejer2
Fejer2700×200 9.29 KB

Is it clearer now ;-)?

Cheers

Régis

3

Hi Régis,

It’s cristal clear as usual!

Thanks for taking the time.
Elias