.ipynb .pdf

repository open issue

Binder Live Code

Contents

2022-04-13 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

fint_trapezoid (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 $-\mathrm{\infty }<a<b<\mathrm{\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=\mathrm{\Delta }x$ instead of number of points $n$.

plot_accuracy_h(fint_trapezoid, tests, 2 .^ (2:10))

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{array}{r}\begin{array}{rl}{I}_h-I_0& =c{h}^2\\ I_{2h}-I_0& =c(2h{)}^2.\end{array}\end{array}$$

Subtracting these two lines, we have

$$I_h-I_{2h}=c(h^2-4h^2)$$

which can be solved for $c$ as

$$c=\frac{I_h-I_{2h}}{h^2-4h^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)

• 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, 4)
plot(x, P)
plot!(x -> 0, color=:black, label=:none)

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,\mathrm{\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:20, ref=1:5)

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, and 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:20, ref=1:4)

$n$-point Gauss exactly integrates polynomials of degree $2n-1$

plot_accuracy(fint_gauss, [(x -> 11x^10, x->x^11)], 3:12)
plot!(xscale=:linear)

plot_accuracy(fint_legendre, [(x -> 11x^10, x->x^11)], 3:12)
plot!(xscale=:linear)

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")

Trefethen, Six Myths of Polynomial Interpolation and Quadrature

@time gausslegendre(1000000);

  0.023108 seconds (10 allocations: 22.888 MiB)

Transforming integrals

Suppose we have a strictly monotone differentiable function $\varphi :(-\mathrm{\infty },\mathrm{\infty })\to (-1,1)$. Then with $x=\varphi (s)$, our integral transforms as

$${\int }_{-1}^{1}f(x)\mathrm{d}x={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}f(\varphi (s)){\varphi }^{\prime }(s)\mathrm{d}s.$$

previous

2022-04-11 Integration

next

2022-04-15 Transformed Quadrature