Quadrature Rules

In this tutorial we investigate how recurrence coefficients of orthogonal polynomials lead to quadrature rules.

We want to solve the integral

\[I = \int_{-1}^{1} f(t) w(t) \mathrm{d} t,\]

with the weight function

\[w(t) = (1-t)^a (1+t)^b\]

for all $t \in [-1, 1]$ and $a, b > -1$. For the function $f$ we choose

\[f(t) = t^2.\]

To solve the integral we do the following:

  1. Choose number of nodes $N$;
  2. Generate recurrence coefficients;
  3. Generate quadrature rule from those recurrence coefficients.

We will compare Gauss quadrature to Gauss-Radau quadrature and Gauss-Lobatto quadrature.

Make sure to check out this tutorial too.

Let's begin:

using PolyChaos, LinearAlgebra
my_f(t) = t^2
a, b = 1.23, 3.45 # shape parameters of Jacobi weight
int_exact = 0.353897; # reference value
0.353897

Now we compute $N$ recurrence coefficients.

N = 4
α, β = rm_jacobi(N+1, a, b)
([0.33233532934131743, 0.1791854079858716, 0.1120747682907886, 0.0767199517952717, 0.055815332777486494], [1.5640214854839802, 0.11582724334265602, 0.1679536499742999, 0.19497630805780736, 0.210634847295037])

Gauss

The first quadrature rule is Gauss quadrature. This method goes back to Golub and Welsch.

n_gauss, w_gauss = gauss(N, α, β)
int_gauss = dot(w_gauss, my_f.(n_gauss))
print("first point:\t $(n_gauss[1])\n")
print("end point:\t $(n_gauss[end])\n")
print("error Gauss:\t $(int_gauss - int_exact)\n")
first point:	 -0.5166972439999838
end point:	 0.8101563565074797
error Gauss:	 4.202394473518112e-7

Since Gauss quadrature has a degree of exactness of $2N-1$, the value of the integral is exact.

Gauss-Radau

Gauss-Radau quadrature is a variant of Gauss quadrature that allows to specify a value of a node that has to be included. We choose to include the right end point $t = 1.0$.

n_radau, w_radau = radau(N-1, α, β, 1.)
int_radau = dot(w_radau, my_f.(n_radau))
print("first point:\t $(n_radau[1])\n")
print("end point:\t $(n_radau[end])\n")
print("error Radau:\t $(int_radau - int_exact)")
first point:	 -0.42969928403284124
end point:	 1.0
error Radau:	 4.2023944762936694e-7

Gauss-Lobatto

Next, we look at Gauss-Lobatto quadrature, which allows to include two points. We choose to include the left and end point of the interval, which are $t \in [-1.0, 1.0]$.

n_lob, w_lob = lobatto(N-2, α, β, -1., 1.)
int_lob = dot(w_lob, my_f.(n_lob))
print("first point:\t $(n_lob[1])\n")
print("end point:\t $(n_lob[end])\n")
print("error Lobatto:\t $(int_lob - int_exact)")
first point:	 -1.0
end point:	 1.0
error Lobatto:	 4.2023944774038924e-7

There are other quadratures that we subsume as all-purpose quadrature rules. These include Fejér's first and second rule, and Clenshaw-Curtis quadrature.

Fejér's First Rule

Fejér's first rule does not include the end points of the interval.

n_fej, w_fej = fejer(N)
int_fej = dot(w_fej, my_f.(n_fej).*(1 .- n_fej).^a.*(1 .+ n_fej).^b)
print("first point:\t $(n_fej[1])\n")
print("end point:\t $(n_fej[end])\n")
print("error Fejer:\t $(int_fej-int_exact)")
first point:	 0.9238795325112867
end point:	 -0.9238795325112867
error Fejer:	 -0.05060511879836699

Fejér's Second Rule

Fejér's second rule does include the end points of the interval.

n_fej2, w_fej2 = fejer2(N)
int_fej2 = dot(w_fej2, my_f.(n_fej2).*(1 .- n_fej2).^a.*(1 .+ n_fej2).^b)
print("first point:\t $(n_fej2[1])\n")
print("end point:\t $(n_fej2[end])\n")
print("error Fejer2:\t $(int_fej2 - int_exact)")
first point:	 1.0
end point:	 -1.0
error Fejer2:	 0.12124113856343288

Clenshaw-Curtis

Clenshaw-Curtis quadrature is similar to Féjer's second rule, as in it includes the end points of the integration interval. For the same number of nodes it is also more accurate than Féjer's rules, generally speaking.

n_cc, w_cc = clenshaw_curtis(N)
int_cc = dot(w_cc, my_f.(n_cc).*(1 .- n_cc).^a.*(1 .+ n_cc).^b)
print("first point:\t\t $(n_cc[1])\n")
print("end point:\t\t $(n_cc[end])\n")
print("error Clenshaw-Curtis:\t $(int_cc - int_exact)")
first point:		 1.0
end point:		 -1.0
error Clenshaw-Curtis:	 0.026213510850746302

As we can see, for the same number of nodes $N$, the quadrature rules based on the recurrence coefficients can greatly outperform the all-purpose quadratures. So, whenever possible, use quadrature rules based on recurrence coefficients of the orthogonal polynomials relative to the underlying measure. Make sure to check out this tutorial too.