2023-04-14 Quadrature#

Last time#

  • Midpoint and trapezoid rules

  • Extrapolation

Today#

  • Polynomial interpolation for integration

  • Gauss quadrature

using LinearAlgebra
using Plots
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

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

function plot_accuracy(fint, tests, ns; ref=[1,2])
    a, b = -2, 2
    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

function fint_trapezoid(f, a, b; n=20)
    dx = (b - a) / (n - 1)
    x = LinRange(a, b, n)
    fx = f.(x)
    fx[1] /= 2
    fx[end] /= 2
    sum(fx) * dx
end

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

Integration#

We’re interested in computing definite integrals

\[ \int_a^b f(x) dx \]

and will usually consider finite domains \(-\infty < a <b < \infty\).

  • Cost: (usually) how many times we need to evaluate the function \(f(x)\)

  • Accuracy

    • compare to a reference value

    • compare to the same method using more evaluations

  • Consideration: how smooth is \(f\)?

Extrapolation#

Let’s switch our plot around to use \(h = \Delta x\) instead of number of points \(n\).

plot_accuracy_h(fint_trapezoid, tests, 2 .^ (2:10))
../_images/2023-04-14-quadrature_4_0.svg

Extrapolation math#

The trapezoid rule with \(n\) points has an interval spacing of \(h = 1/(n-1)\). Let \(I_h\) be the value of the integral approximated using an interval \(h\). We have numerical evidence that the leading error term is \(O(h^2)\), i.e.,

\[ I_h - I_0 = c h^2 + O(h^3) \]
for some as-yet unknown constant \(c\) that will depend on the function being integrated and the domain of integration. If we can determine \(c\) from two approximations, say \(I_h\) and \(I_{2h}\), then we can extrapolate to \(h=0\). For sufficiently small \(h\), we can neglect \(O(h^3)\) and write
\[\begin{split}\begin{split} I_h - I_0 &= c h^2 \\ I_{2h} - I_0 &= c (2h)^2 . \end{split}\end{split}\]
Subtracting these two lines, we have
\[ I_{h} - I_{2h} = c (h^2 - 4 h^2) \]
which can be solved for \(c\) as
\[ c = \frac{I_{h} - I_{2h}}{h^2 - 4 h^2} . \]
Substituting back into the first equation, we solve for \(I_0\) as
\[ I_0 = I_h - c h^2 = I_h + \frac{I_{h} - I_{2h}}{4 - 1} .\]
This is called Richardson extrapolation.

Extrapolation code#

function fint_richardson(f, a, b; n=20)
    n = div(n, 2) * 2 + 1
    h = (b - a) / (n - 1)
    x = LinRange(a, b, n)
    fx = f.(x)
    fx[[1, end]] /= 2
    I_h = sum(fx) * h
    I_2h = sum(fx[1:2:end]) * 2h
    I_h + (I_h - I_2h) / 3
end
plot_accuracy_h(fint_richardson, tests, 2 .^ (2:10), ref=1:5)
../_images/2023-04-14-quadrature_7_0.svg
  • we now have a sequence of accurate approximations

  • it’s possible to apply extrapolation recursively

  • works great if you have a power of 2 number of points

    • and your function is nice enough

Quadrature form#

At the end of the day, we’re taking a weighted sum of function values. We call \(w_i\) the quadrature weights and \(x_i\) the quadrature points or abscissa.

\[\int_a^b f(x) \approx \sum_{i=1}^n w_i f(x_i) = \mathbf w^T f(\mathbf x)\]
function quad_trapezoid(a, b; n=20)
    dx = (b - a) / (n - 1)
    x = LinRange(a, b, n)
    w = fill(dx, n)
    w[[1, end]] /= 2
    x, w
end
quad_trapezoid (generic function with 1 method)
x, w = quad_trapezoid(-1, 1)

w' * cos.(x) - fint_trapezoid(cos, -1, 1)
2.220446049250313e-16

Polynomial interpolation for integration#

x = LinRange(-1, 1, 100)
P = vander_legendre(x, 10)
plot(x, P[:,end], legend=:none)
../_images/2023-04-14-quadrature_13_0.svg

Idea#

  • Sample the function \(f(x)\) at some points \(x \in [-1, 1]\)

  • Fit a polynomial through those points

  • Return the integral of that interpolating polynomial

Question#

  • What points do we sample on?

  • How do we integrate the interpolating polynomial?

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.\]

Integration using Legendre polynomials#

function quad_legendre(a, b; n=20)
    x = CosRange(-1, 1, n)
    P = vander_legendre(x)
    x_ab = (a+b)/2 .+ (b-a)/2*x
    w = (b - a) * inv(P)[1,:]
    x_ab, w
end

function fint_legendre(f, a, b; n=20)
    x, w = quad_legendre(a, b, n=n)
    w' * f.(x)
end

fint_legendre(x -> 1 + x, -1, 1, n=4)
2.0
p = plot_accuracy(fint_legendre, tests, 4:50, ref=1:5)
plot!(p, xscale=:linear)
../_images/2023-04-14-quadrature_17_0.svg

Doing better#

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

\[ 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\).

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

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

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

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

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\)?

Gauss-Legendre in code#

  1. Solve for the points, compute the weights

  • Use a Newton solver to find the roots. You can use the recurrence to write a recurrence for the derivatives.

  • Create a Vandermonde matrix and extract the first row of the inverse or (using more structure) the derivatives at the quadrature points.

  1. Use duality of polynomial roots and matrix eigenvalues.

  • A fascinating mathematical voyage; something you might see more in a graduate linear algebra class.

function fint_gauss(f, a, b; n=4)
    """Gauss-Legendre integration using Golub-Welsch algorithm"""
    beta = @. .5 / sqrt(1 - (2 * (1:n-1))^(-2))
    T = diagm(-1 => beta, 1 => beta)
    D, V = eigen(T)
    w = V[1,:].^2 * (b-a)
    x = (a+b)/2 .+ (b-a)/2 * D
    w' * f.(x)
end
fint_gauss(sin, -2, 3, n=4)
0.5733948071694299
plot_accuracy(fint_gauss, tests, 3:30, ref=1:4)
plot!(xscale=:linear)
../_images/2023-04-14-quadrature_21_0.svg

\(n\)-point Gauss exactly integrates polynomials of degree \(2n-1\)#

plot_accuracy(fint_gauss, [(x -> 12(x-.2)^11, x->(x-.2)^12)], 3:12)
plot!(xscale=:linear)
../_images/2023-04-14-quadrature_23_0.svg
plot_accuracy(fint_legendre, [(x -> 11x^10, x->x^11)], 3:12)
plot!(xscale=:linear)
../_images/2023-04-14-quadrature_24_0.svg

FastGaussQuadrature.jl#

using FastGaussQuadrature

n = 100
x, q = gausslegendre(n)
scatter(x, q, label="Gauss-Legendre", ylabel="weight", xlims=(-1, 1))
scatter!(gausslobatto(n)..., label="Gauss-Lobatto")
../_images/2023-04-14-quadrature_26_0.svg