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

**URL:** <https://openturns.discourse.group/t/is-the-fejer-clenshaw-curtis-quadrature-actually-nested/299>\
**Category:** Methodology\
**Created:** [June 20, 2023, 11:57am UTC](https://openturns.discourse.group/t/is-the-fejer-clenshaw-curtis-quadrature-actually-nested/299 "2023-06-20T11:57:41Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![efekhari27](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/efekhari27/32/92_2.png) [@efekhari27](https://openturns.discourse.group/u/efekhari27)\
**Post date:** [June 20, 2023, 11:57am UTC](https://openturns.discourse.group/t/is-the-fejer-clenshaw-curtis-quadrature-actually-nested/299/1 "2023-06-20T11:57:41Z")

</div>

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:

```auto
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

---

<div class="post-metadata">

**Author:** ![regislebrun](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/regislebrun/32/154_2.png) [@regislebrun](https://openturns.discourse.group/u/regislebrun)\
**Post date:** [June 21, 2023, 5:47pm UTC](https://openturns.discourse.group/t/is-the-fejer-clenshaw-curtis-quadrature-actually-nested/299/2 "2023-06-21T17:47:27Z")

</div>

Hi Elias,

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

> **[Clenshaw–Curtis quadrature](https://en.wikipedia.org/wiki/Clenshaw%E2%80%93Curtis_quadrature)**
>
> Clenshaw–Curtis quadrature and Fejér quadrature are methods for numerical integration, or "quadrature", that are based on an expansion of the integrand in terms of Chebyshev polynomials. Equivalently, they employ a change of variables 
>   
>     
>       
> x
> =
> cos
> ⁡
> θ
>       
>     
> {\\displaystyle x=\\cos \\theta }
>   
> and use a discrete cosine transform (DCT) approximation for the cosine series. Besides having fast-converging accuracy comparable to Gaussian qu...

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:

```auto
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](https://global.discourse-cdn.com/free1/uploads/openturns/original/1X/ddd9a631cf9defa60b6c8d98bad33601caa0daa4.png)

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:

```auto
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](https://global.discourse-cdn.com/free1/uploads/openturns/original/1X/26a9bb9d2a354b25eba31369a5ab1e5e13c7c1e1.png)

Is it clearer now ;-)?

Cheers

Régis

---

<div class="post-metadata">

**Author:** ![efekhari27](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/efekhari27/32/92_2.png) [@efekhari27](https://openturns.discourse.group/u/efekhari27)\
**Post date:** [June 21, 2023, 6:33pm UTC](https://openturns.discourse.group/t/is-the-fejer-clenshaw-curtis-quadrature-actually-nested/299/3 "2023-06-21T18:33:56Z")

</div>

Hi Régis,

It’s cristal clear as usual!

Thanks for taking the time.  
Elias