Gaussian quadrature

import numpy as np
import sympy as sp
import plotly.graph_objects as go

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

Optimal Evaluation Nodes: The Gauss Points¶

Up until this point, we have been using techniques for which the data points are fixed, (hopefully) evenly spaced, and at the limits of our subdomains.

Gaussian quadrature assumes you are free to choose the evaluation points $x_i$ . The goal is then to find special sets of points $x_i$ (called Gauss points) and weights $c_i$ for which the integral is computed with maximal polynomial degree precision.

The 2-Point Gauss-Legendre Formula¶

Let’'s propose a formula for the integral of a function $f(x)$ in the range $[-1, 1]$ :

I = \int_{-1}^1 f(x) \,\mathrm{d}x \approx c_0 f(x_0) + c_1 f(x_1)

(1)

where the unknowns are $c_0, c_1, x_0, x_1$ (4 degrees of freedom). We require this formula to integrate monomials $1, x, x^2, x^3$ exactly:

\begin{aligned} \int_{-1}^1 1 \,\mathrm{d}x &= 2 = c_0 + c_1 \\ \int_{-1}^1 x \,\mathrm{d}x &= 0 = c_0 x_0 + c_1 x_1 \\ \int_{-1}^1 x^2 \,\mathrm{d}x &= \frac{2}{3} = c_0 x_0^2 + c_1 x_1^2 \\ \int_{-1}^1 x^3 \,\mathrm{d}x &= 0 = c_0 x_0^3 + c_1 x_1^3 \end{aligned}

(2)

Solving this system yields:

c_0 = c_1 = 1, \quad x_0 = -\frac{1}{\sqrt{3}}, \quad x_1 = \frac{1}{\sqrt{3}}

(3)

That is a remarkable result! Any cubic polynomial can be integrated exactly merely through the sum of the function evaluated at two points!

def two_point_gaussian_quadrature(f):
    return f(-1.0 / np.sqrt(3)) + f(1.0 / np.sqrt(3))

x = sp.var("x")

def quad_tester(expr):
    fcn = sp.lambdify(x, expr, modules=["numpy"])
    fI = float(sp.integrate(expr, (x, -1, 1)))
    num_val = float(two_point_gaussian_quadrature(fcn))
    print(f"Function: {expr} | Exact: {fI:.6f} | Error: {num_val - fI:.2e}")

quad_tester(x**2)
quad_tester(x**3 - x**2 + sp.pi * x + 17)
quad_tester(-10**9 * x**3 + sp.exp(1) * x**2 - sp.sqrt(sp.pi) * x + 42)

Transformation to an Arbitrary Domain $[a, b]$ ¶

The preceding case specified integration limits on $[-1, 1]$ . To apply Gaussian quadrature to arbitrary integration limits $[a, b]$ , we perform a coordinate change:

x = \frac{b + a}{2} + \frac{b - a}{2} \xi, \quad \mathrm{d}x = \frac{b - a}{2} \,\mathrm{d}\xi

(4)

where $\xi \in [-1, 1]$ .

def flex_gq(f, a, b):
    x0 = -1.0 / np.sqrt(3)
    x1 = 1.0 / np.sqrt(3)
    x_d0 = ((b + a) + (b - a) * x0) / 2.0
    x_d1 = ((b + a) + (b - a) * x1) / 2.0
    dx = (b - a) / 2.0
    return (f(x_d0) + f(x_d1)) * dx

def flex_gq_tester(expr, a, b):
    fcn = sp.lambdify(x, expr, modules=["numpy"])
    exact = float(sp.integrate(expr, (x, a, b)))
    approx = float(flex_gq(fcn, a, b))
    print(f"Integral of {expr} on [{a}, {b}]: Exact = {exact:.6f}, Approx = {approx:.6f}, Error = {approx - exact:.2e}")

flex_gq_tester(x**2, 0, 1)
flex_gq_tester(x**3 - x**2 + sp.pi * x + 17, 0, 1)

The $n$ -point Gauss-Legendre Formula¶

Jinkies!

The method may be generalized to $n$ -point quadrature that integrates polynomials up to degree $2n - 1$ exactly:

I = \sum_{i=1}^n c_i f(x_i)

(5)

with $c_i$ given by Legendre polynomial weights:

Integration order	Coefficient ( $c_i$ )	Gauss Point ( $x_i$ )	Error Term
2	$1, 1$	$\pm \frac{1}{\sqrt{3}}$	$\sim f^{(4)}$
3	$\frac{8}{9}, \frac{5}{9}, \frac{5}{9}$	$0, \pm \sqrt{\frac{3}{5}}$	$\sim f^{(6)}$
4	$0.347855, 0.652145$	$\pm 0.861136, \pm 0.339981$	$\sim f^{(8)}$
5	$0.236927, 0.478629, 0.568889$	$\pm 0.906180, \pm 0.538469, 0$	$\sim f^{(10)}$

Error Analysis for Gauss-Legendre Quadrature¶

Note that the error is given in terms of higher derivatives of $f$ , since the step size is no longer an issue. In general, the error is:

Error = \frac{2^{2n+1} (n!)^4}{(2n+1)[(2n)!]^3} f^{(2n)}(\xi)

(6)

for some $\xi \in [-1, 1]$ .

Multidimensional Gauss-Legendre Quadrature¶

Due to the separable nature of multidimensional integration, extending this concept to $N$ dimensions is an exercise in repeated application:

\int_{-1}^1 \int_{-1}^1 f(\alpha, \beta) \,\mathrm{d}\alpha \,\mathrm{d}\beta = \int_{-1}^1 \left[\int_{-1}^1 f(\alpha, \beta) \,\mathrm{d}\beta\right] \mathrm{d}\alpha \approx \sum_i c_i \sum_j c_j f(\alpha_i, \beta_j)

(7)

Quadrilateral Elements¶

An example of a quadrilateral with a 3rd integration order is shown below, illustrating the mapping between computational parent coordinates (left) and real physical coordinates (right):