32) Numerical Intergration

32) Numerical Intergration#

Last time:

Splines conditioning
Boundary Value Problems
Higher dimensions

Today:

Numerical Integration
1.1 Midpoint rule
1.2 Trapezoid rule
Extrapolation
Polynomial interpolation for integration

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

CosRange (generic function with 1 method)

1. Numerical Integration#

We’re interested in computing definite integrals

\int_{a}^{b} f (x) d x

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

Observations:

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$ ?

Some test functions#

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)

zip(Function[f_expx, f_dtanh], Function[F_expx, F_dtanh])

plot(integrands, xlims=(-1, 1))

Fundamental Theorem of Calculus#

Let $f (x)$ be a continuous function and define $F (x)$ by

F (x) = \int_{a}^{x} f (s) d s .

Then

F (x)

is uniformly continuous, differentiable, and

F^{'} (x) = f (x) .

We say that

F

is an antiderivative of

f

. This implies that

\int_{a}^{b} f (x) d x = F (b) - F (a) .

We will test the accuracy of our integration schemes using an antiderivative provided in our tests.

Method of Manufactured Solutions#

Analytically integrating an arbitrary function is hard
- tends to require trickery (e.g., u substitution)
- not always possible to express in closed form (e.g., elliptic integrals)
- sometimes needs special functions $\erf x = \frac{2}{\sqrt{π}} \int_{0}^{x} e^{- t^{2}} d t$
- don’t know when to give up
On the other hand, analytic differentation
- involves straightforward application of the product rule and chain rule.

So if we just choose an arbitrary function $F$ (the antiderivative), we can

compute $f = F^{'}$
numerically integrate $\int_{a}^{b} f$ and compare to $F (b) - F (a)$

Newton-Cotes methods#

Main idea: Approximate $f (x)$ using piecewise polynomials (an interpolation problem) on equispaced nodes and integrate the polynomials.

1.1 Midpoint Riemann Sum method #

function fint_midpoint(f, a, b; n=20)
    dx = (b - a) / n
    x = LinRange(a + dx/2, b - dx/2, n)
    sum(f.(x)) * dx
end

for (f, F) in tests
    a, b = -2, 2
    I_num = fint_midpoint(f, a, b, n=20)
    I_analytic = F(b) - F(a)
    println("$f: $I_num error=$(I_num - I_analytic)")
end

f_expx: 10.885522849146847 error=-0.03044402970425253
f_dtanh: 1.9285075531458646 error=0.00045239299423083246

How does the accuracy change as we use more points?

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
plot_accuracy(fint_midpoint, tests, 2 .^ (0:10))

1.2 Trapezoid Rule#

The trapezoid rule uses piecewise linear functions on each interval.

\begin{array}{r} \begin{aligned} \int_{a}^{b} [f (a) + \frac{f (b) - f (a)}{b - a} (x - a)] & = f (a) (x - a) + \frac{f (b) - f (a)}{2 (b - a)} (x - a)^{2} |_{x = a}^{b} \\ = f (a) (b - a) + \frac{f (b) - f (a)}{2 (b - a)} (b - a)^{2} \\ = \frac{b - a}{2} (f (a) + f (b)) . \end{aligned} \end{array}

Can you get to the same result using a geometric argument?
What happens when we sum over a bunch of adjacent intervals?

Code for the Trapezoid Rule:

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

plot_accuracy(fint_trapezoid, tests, 2 .^ (1:10))

2. Extrapolation#

Let’s switch our plot around to use $h = Δ x$ instead of number of points $n$ .

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, .1 * hs.^k, label="\$h^{$k}\$")
    end
    p
end

plot_accuracy_h (generic function with 1 method)

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

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_{2 h}$ , then we can extrapolate to $h = 0$ .

For sufficiently small $h$ , we can neglect $O (h^{3})$ and write

\begin{array}{r} \begin{aligned} I_{h} - I_{0} & = c h^{2} \\ I_{2 h} - I_{0} & = c (2 h)^{2} . \end{aligned} \end{array}

Subtracting these two lines, we have

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

which can be solved for $c$ as

c = \frac{I_{h} - I_{2 h}}{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_{2 h}}{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)

Observations:

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

F_sin(x) = sin(30x)
f_sin(x) = 30*cos(30x)
plot_accuracy_h(fint_richardson, zip([f_sin], [F_sin]), 2 .^ (2:10), ref=1:5)

3. Polynomial interpolation for integration#

x = LinRange(-1, 1, 100)
P = vander_legendre(x, 10)
plot(x, P, legend=:none)

Main ideas:

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

Questions:

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, \forall m \neq n .

Integration using Legendre polynomials#

function plot_accuracy_n(fint, tests, ns; ref=[1,2])
    a, b = -2, 2
    p = plot(xscale=:log10, yscale=:log10, xlabel="n", ylabel="error",
        legend=:bottomright)
    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
    p
end

plot_accuracy_n (generic function with 1 method)

function fint_legendre(f, a, b; n=20)
    x = CosRange(-1, 1, n)
    P = vander_legendre(x)
    x_ab = (a+b)/2 .+ (b-a)/2*x
    c = P \ f.(x_ab)
    (b - a) * c[1]
end

fint_legendre(x -> 1 + x, -1, 1, n=2)

2.0

p = plot_accuracy_h(fint_legendre, tests, 4:20, ref=1:5)