34) Transformed Quadrature

34) Transformed Quadrature#

Last time#

  • Polynomial interpolation for integration

  • Gauss quadrature

Today#

  1. Recap on Gauss Quadrature

  2. Activity: six myths of polynomial interpolation and quadrature

  3. Singular integrals and tanh-sinh quadrature

  4. Adaptive integration

  5. Integration in multiple dimensions

using LinearAlgebra
using Plots
using LaTeXStrings
default(linewidth=4, legendfontsize=12)

function vander_legendre(x, k=nothing)
    if isnothing(k)
        k = length(x) # Square by default
    end
    m = length(x)
    Q = ones(m, k)
    Q[:, 2] = x
    for n in 1:k-2
        Q[:, n+2] = ((2*n + 1) * x .* Q[:, n+1] - n * Q[:, n]) / (n + 1)
    end
    Q
end

CosRange(a, b, n) = (a + b)/2 .+ (b - a)/2 * cos.(LinRange(-pi, 0, n))

F_expx(x) = exp(2x) / (1 + x^2)
f_expx(x) = 2*exp(2x) / (1 + x^2) - 2x*exp(2x)/(1 + x^2)^2

F_dtanh(x) = tanh(x)
f_dtanh(x) = cosh(x)^-2

F_rsqrt(x) = 2 * sqrt(x + 1)
f_rsqrt(x) = 1 / sqrt(x + 1)

integrands = [f_expx, f_dtanh, f_rsqrt]
antiderivatives = [F_expx, F_dtanh, F_rsqrt]
tests = zip(integrands, antiderivatives)

function plot_accuracy(fint, tests, ns; ref=[1,2])
    a, b = -1, 1
    p = plot(xscale=:log10, yscale=:log10, xlabel="n", ylabel="error")
    for (f, F) in tests
        Is = [fint(f, a, b, n=n) for n in ns]
        Errors = abs.(Is .- (F(b) - F(a)))
        scatter!(ns, Errors, label=f)
    end
    for k in ref
        plot!(ns, ns.^(-1. * k), label="\$n^{-$k}\$")
    end
    p
end

plot_accuracy (generic function with 1 method)

1. Recap on Gauss Quadrature#

Newton–Cotes quadrature corresponds to equispaced points, Clenshaw–Curtis quadrature to Chebyshev points,and Gauss quadrature to Legendre points.#

We want to approximate

\[\int_{-1}^1 f(x) dx \approx \sum_{i=1}^n w_i f(x_i) = \mathbf w^T f(\mathbf x)\]

We call \(w_i\) the quadrature weights and \(x_i\) the quadrature points or abscissa.

How do we choose the points?

Recall that the Legendre polynomials \(P_0(x) = 1\), \(P_1(x) = x\), …, are pairwise orthogonal

\[\int_{-1}^1 P_m(x) P_n(x) = 0, \quad \forall m\ne n.\]

How does this work?#

Suppose a polynomial \(P(x)\) on the interval \([-1,1]\) can be written as

\[ P(x) = P_n(x) q(x) + r(x) \]

where \(P_n(x)\) is the \(n\) th Legendre polnomials and both \(q(x)\) and \(r(x)\) are polynomials of maximum degree \(n-1\). We know we can write this form by using long division of polynomials.

  • Why is \(\int_{-1}^1 P_n(x) q(x) = 0\)?

    • Since Gaussian quadrature is exact on the polynomial \(r(x)\), since it is just integration of the interpolating polynomial of degree \(n-1\), which is identical to \(r(x)\).

  • Can every polynomials of degree \(2n-1\) be written in the above form?

    • Yes, by long division of polynomials.

  • How many roots does \(P_n(x)\) have on the interval?

    • \(n\) complex-valued, considering their multiplicity (Fundamental Theorem of Algebra: every non-zero, single-variable, degree \(n\) polynomial with complex coefficients has, counted with multiplicity, exactly \(n\) complex roots.).

  • Can we choose points \(\{x_i\}\) such that the first term is 0?

    • At the roots \(x_i\) of the \(n\)-th Legendre polynomial, \(P(x_i) = r(x_i)\), since \(P_n (x_i ) = 0\) for all \(i\).

If \(P_n(x_i) = 0\) for each \(x_i\), then we need only integrate \(r(x)\), which is done exactly by integrating its interpolating polynomial. How do we find these roots \(x_i\)?

We approximate

\[ \int_{-1}^{-1} f(x) dx \approx \sum_{i=1}^{n} w_i f(x_i) \]

where \(w_i\) are the weights defined as

\[ w_i = \int_{-1}^{1}L_i(x) dx \quad i=1, \dots, n \]

the Lagrange interpolating polynomial with \(n\) interpolating nodes, and the quadrature nodes \(x_i\) are the roots of the \(n\)-th Legendre polynomial on \([-1,1]\).

Note:

To approximate integrals on a general interval \([a,b]\), the problem needs to be translated to \([-1,1]\) with a simple change of variables. Using the substitution \(t = (2x - a - b) / (b-a)\), we find

\[ \int_a^b f(x) dx = \int_{-1}^{1} f \left( \frac{(b-a)t + b + a}{2} \right) \frac{(b-a)}{2} dt \]

Julia packages#

Recall the Julia package FastGaussQuadrature.jl that offers different quadrature formulae.

using FastGaussQuadrature

n = 10
x, q = gausslegendre(n)
scatter(x, q, label="Gauss-Legendre", ylabel="weight", xlims=(-1, 1))
scatter!(gausslobatto(n)..., label="Gauss-Lobatto")

Spectral Element Method (SEM)#

Gauss-Lobatto nodes are preferred in the quadrature formulae for the Spectral Element Method (SEM), a high-order (“spectral” or “exponential”) method for the solution of Partial Differential Equations.

In finite element formulations, the weak form of a Partial Differential Equation (PDE) is evaluated on a subdomain \(\Omega_e\) (element) and the local results are composed into a larger system of equations that models the entire problem on the global domain \(\Omega\). The weak-form involves the \(L^2\) inner product (integral) of a variational formulation with suitable basis functions and trial/test functions.

Gauss-Lobatto nodes are preferred because they include nodes at the endpoints \(-1\) and \(1\), hence ensuring continuity of the basis functions across element boundaries.

The more general formulation of SEM, that instead of special points, uses equispaced points is called the Finite Element Method (FEM).

Finite Element Methods (FEM)#

Main ideas:

  • Need to integrate the product of basis functions over every element

  • Choose a quadrature (e.g., Gauss) to minimize number of points for sufficient accuracy

  • Bad effects if insufficient points

2. Activity 34.1: Read Trefethen, Six Myths of Polynomial Interpolation and Quadrature and discuss.#

# how fast/slow if the FastGaussQuadrature.jl package implementation?

@time gausslegendre(10000000);

  0.305080 seconds (26 allocations: 228.883 MiB, 12.02% gc time, 1.12% compilation time)

3. Singular integrands#

plot([sqrt log x->.5*x^(-.5)], xlim=(0, 2), ylim=(-4, 4), label = [L"\sqrt{x}" L"\log{x}" L"\frac{1}{2x^5}"])

function fint_gauss(f, a, b, n)
    x, w = gausslegendre(n)
    x = (a+b)/2 .+ (b-a)/2*x # change of variables
    w *= (b - a)/2 # change of variables
    w' * f.(x)
end

fint_gauss (generic function with 1 method)

plot(3:4:100,
    n -> abs(fint_gauss(x -> .5*x^(-.5), 0, 1, n) - 1),
    marker=:auto, yscale=:log10, xscale=:log10, label = L"\int_{0}^{1}\frac{1}{2x^5} = 1")
plot!(n -> 1/n,  label = L"\frac{1}{n}", title = "Error")

Tanh-Sinh quadrature#

Main idea: make everything smooth.

When functions have singularities near the endpoints, it is usually more efficient to integrate via a change of variables. Suppose we have a strictly monotone differentiable function \(\phi: (-\infty, \infty) \to (-1, 1)\).

Then with \(x = \phi(s)\), our integral transforms as

\[ \int_{-1}^1 f(x) \mathrm dx = \int_{-\infty}^\infty f(\phi(s)) \phi'(s) \mathrm d s . \]

The tanh-sinh method uses a transformation such that \(\phi'(s) \to 0\) faster than the singularity \(f(\phi(s))\) grows, such that the integrand goes to 0 at finite \(s\).

tanhsinh(s) = tanh(pi/2*sinh(s))

function dtanhsinh(s)
    ds = 1
    t = pi/2 * sinh(s)
    dt = pi/2 * cosh(s) * ds
    (1 - tanh(t)^2) * dt
end

p = plot([tanhsinh], color=:black, label="tanhsinh(s)",
    xlims=(-3, 3),
    xlabel="s", title="tanh-sinh function and integrands")
for f in integrands
    plot!(s -> f(tanhsinh(s))*dtanhsinh(s), label="$f ∘ tanhsinh")
end
p

Implementation#

The function below implements tanh-sinh quadrature on the interval \((-1,1)\). Given the number of points, we need to choose both the limits of integration (we can’t afford to integrate all the way to infinity) and the spacing. Here we make an arbitrary choice to integrate on the interval \((-L, L)\) where \(L = \log n\). The grid spacing thus scales as \(h \approx 2 \log n / n\).

Modify the quadrature so it can be used to integrate on an arbitrary interval \((a,b)\).

function fint_tanhsinh(f, a, b; n=9)
    L = log(n)
    h = 2 * L / (n - 1)
    s = LinRange(-L, L, n)
    x = tanhsinh.(s)
    w = h * dtanhsinh.(s)

    ## Challenge: modify the weights w and points x to integrated on (a,b), not (-1, 1)

    w' * f.(x)
end

fint_tanhsinh (generic function with 1 method)

plot_accuracy(fint_tanhsinh, tests, 9:4:60, ref=[2,3,4])
plot!(xscale=:identity,)

# If you complete the challenge above uncomment the following two lines

# @assert fint_tanhsinh(log, 0, 1, n=20) ≈ -1
# println("Tests pass")

Functions that are rough or discontinuous in the interior#

If a function has rough or singular behavior in the interior of the domain, we’ll get low accuracy results.

using Pkg
Pkg.add("SpecialFunctions")

using SpecialFunctions

f_heaviside(x) = 1.0 * (x > 0)
F_heaviside(x) = max(x, 0) # ramp function is the antiderivative of the Heaviside function

f_mollify(x, sigma=.02) = 1/(sigma * sqrt(2*π)) * exp(-.5*(x/sigma)^2)
F_mollify(x, sigma=.02) = .5*erf(x/(sigma*sqrt(2))) # erf is the Gauss error function

rough_tests = [(f_heaviside, F_heaviside), (f_mollify, F_mollify)]
plot([f_heaviside, f_mollify], xlim=(-1, 1), label = ["Heaviside" "Mollify"])

   Resolving package versions...

  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`

plot_accuracy(fint_tanhsinh, rough_tests, 5:4:800, ref=[2,3,4]); plot!(legend=:bottomleft)

4. Adaptive integration#

H-adaptive integration computes the integral, adaptively subdividing the integration area/volume into smaller and smaller pieces until convergence is achieved to the desired tolerance.

A Julia package: HCubature.jl#

using Pkg
Pkg.add("HCubature")

using HCubature

@show F_heaviside(1 - .3) - F_heaviside(-1 - .3)
hquadrature(x -> f_heaviside(x - 0.3), -1, 1, maxevals=1000) #  The return value of hcubature is a tuple (I, E) of the estimated integral I and an estimated error E.

   Resolving package versions...

  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`

F_heaviside(1 - 0.3) - F_heaviside(-1 - 0.3) = 0.7

(0.700000008169331, 9.264907000048497e-9)

hquadrature(f_mollify, -1, 1, maxevals=1000)

(1.0000000000000002, 6.773756033519694e-10)

F_mollify(1) - F_mollify(-1)

1.0

Quadrature/cubature for multi-dimensional integration in Python: quadpy#

Notes on integration#

  • Transforms can make the integrand smooth

  • Transforms can make the domain shape more convenient

  • Adaptive integration

  • Curse of dimensionality

    • Sparse grids (Smolyak quadrature)

    • Adaptive randomized methods (Markov Chain Monte Carlo)

Contents