Lecture 27: Error of quadrature

Last lecture we observed the following:

  1. The right-hand rule had error $O(n^{-1})$ for $n$ sample point
  2. The trapezium rule had error $O(n^{-2})$ for $n$ sample point
  3. The trapezium rule had error $O(n^{-\alpha})$ for all choices of $\alpha \geq 0$ for smooth, periodic functions

This lecture we set out to prove these results.

Right-hand rule convergence

Recall that the right-hand rule is defined by, for $x_k \triangleq kh$ and $h \triangleq 1/n$

$$\int_0^1 f(x) dx \approx h \sum_{k=1}^n f(x_k) = \sum_{k=1}^n \int_{x_{k-1}}^{x_k} (f(x) - f(x_k) ) dx$$

About the only tool available to analyse integrals is integration by parts:

$$ \begin{align*} \int_a^b u(x) v'(x) dx &= [u(x) v(x)]_a^b - \int_a^b u'(x) v(x) \cr & = u(b) v(b) - u(a)v(a) - \int_a^b u'(x) v(x) \end{align*} $$

Consider the error in the first panel. Taking $u(x) = f(x)$ and $v(x) = x$ in the integration by parts formula yields:

$$ \begin{align*} \int_0^h (f(x) - f(h)) dx = [(f(x) - f(h)) x]_0^h - \int_0^h f'(x) x dx \cr = - \int_0^h f'(x) x dx \end{align*} $$

Assume that $f'(x)$ is bounded in $[0,1]$, and define

$$M= \sup_{x \in [0,1]} |f'(x)|$$

We then get a bound on the error:

$$\begin{align*} \left|\int_0^h (f(x) - f(h)) dx\right| &= \left|\int_0^h f'(x) x dx\right| \leq \int_0^h |f'(x)| x dx \leq M \int_0^h x dx \cr &\leq {M h^2 \over 2} \end{align*} $$
    

This formula caries over to every other panel (by the same argument):

$$ \left|\int_{x_{k-1}}^{x_k} (f(x) - f(x_k)) dx \right| \leq {M h^2 \over 2} $$

Thus we have

$$ \begin{align*} \left\|\int_0^1 f(x) dx - h \sum_{k=1}^n f(x_k)\right\| &= \left\|\sum_{k=1}^n \int_{x_{k-1}}^{x_k} (f(x) - f(x_k) ) dx\right\| \cr &\leq \sum_{k=1}^n \left\|\int_{x_{k-1}}^{x_k} (f(x) - f(x_k) ) \right\| dx &\leq \sum_{k=1}^n {M h^2 \over 2} = {M h^2 n \over 2} = {M \over 2 n} = O(n^{-1}) \end{align*} $$

Thus we have proven that the right-hand rule converges like $O(n^{-1})$.

Trapezium rule observed convergence, revisited

Recall the trapezium rule, here implemented for general intervals $[a,b]$:

In [2]:
using PyPlot


function trap(f,a,b,n)
    h=(b-a)/n
    x=linspace(a,b,n+1)

    v=f(x)
    h/2*v[1]+sum(v[2:end-1]*h)+h/2*v[end]
end

trap(f,n)=trap(f,0.,1.,n);

Consider integration of the following simple function:

In [3]:
f=x->exp(x).*cos(x)

g=linspace(0.,1.,1000)
plot(g,f(g));

As in last lecture, we can compare the error of trap to the exact integral, approximated by quadgk.

In [5]:
ex=quadgk(f,0.,1.)[1]
exπ=quadgk(θ->cos(cos(θ-0.1)),0.,2π)[1]   # a periodic integral


ns=round(Int,logspace(1,5,100))   # integers spaced logarithmically apart
errT=zeros(length(ns))
errπ=zeros(length(ns))

for k=1:length(ns)
    n=ns[k]
    errT[k]=abs(trap(f,n-1)-ex    )   # error in non-periodic
    errπ[k]=abs(trap(θ->cos(cos(θ-0.1)),0.,2π,n-1)-exπ    )  # error in periodic
end

We see that for non-periodic functions we observe $O(n^{-2})$ convergence but for periodic functions we converge much faster:

In [6]:
loglog(ns,errT)   # blue curve
loglog(ns,errπ)   # green curve
loglog(ns,(1.0ns).^(-2))   # red curve, with same slope as blue curve
Out[6]:
1-element Array{Any,1}:
 PyObject <matplotlib.lines.Line2D object at 0x304c5ed90>

The defining property of smooth, periodic functions is that the function and all its derivatives match at 0 and 2π. Now consider a function where only some of he derivatives match:

In [8]:
h=x->10f(x).*x.^2.*(1-x).^2

plot(g,h(g));

We see here that the convergence rate has converged to $O(n^{-4})$:

In [11]:
exf=quadgk(f,0.,1.)[1]
exh=quadgk(h,0.,1.)[1]


ns=round(Int,logspace(1,4,100))
errf=zeros(length(ns))
errh=zeros(length(ns))

for k=1:length(ns)
    n=ns[k]
    errf[k]=abs(trap(f,n-1)-exf    )
    errh[k]=abs(trap(h,n-1)-exh    )
end

loglog(ns,errf)     # blue curve
loglog(ns,1./ns.^2) # green curve
loglog(ns,errh)     # red curve
loglog(ns,1./ns.^4) # aqua curve;

The conclusion is that derivatives at the endpoints dictate the convergence rate of the Trapezium rule.

Bernoulli Polynomial

To prove this, we need to be clever about our choice of functions when we integrate by parts. Define the first three Bernoulli Polynomials by

$$B_0(x) = 1, B_1(x) = x-{1 \over 2}, B_2(x) = x^2 -x + {1 \over 6}$$

as depticted here:

In [12]:
g=linspace(0.,1.,1000)
B0=x->ones(x)
B1=x->x-1/2
B2=x->x.^2-x+1/6
plot(g,B0(g))
plot(g,B1(g))
plot(g,B2(g));

These polynomials satisfy two important properties:

  1. $B_k'(x) = k B_{k-1}(x)$, just like monomials $x^k$
  2. $B_2(-1) = B_2(1) = {1 \over 6}

We Thus consider integration by parts with the Trapezium rule. Recall the Trapezium rule:

$$ \int_0^1 f(x) dx \approx h\left({f(x_0) \over 2} + \sum_{k=1}^{n-1} f(x_k) + {f(x_n) \over 2} \right) = \sum_{k=1}^n \int_0^1 \left(f(x_{k-1}) + (x-x_{k-1}) {f(x_k) - f(x_{k-1}) \over h}\right) dx $$

As in the right-hand rule, we begin with the error in the first panel:

$$\int_0^h \left[f(x) - \left(f(0) + x {f(h) - f(0) \over h}\right) \right] dx.$$

Taking

$$u(x) = f(x) - \left(f(0) + x {f(h) - f(0) \over h}\right)$$

and $v(x) = h B_1(x/h)$ (so that $v'(x) = 1$) in the integration by parts formula:

$$ \int_0^h \left[f(x) - \left(f(0) + x {f(h) - f(0) \over h}\right) \right] dx = [u(x) v(x)]_0^h - \int_0^h u(x) v(x) dx = - h \int_0^h \left[f'(x) - {f(h) - f(0) \over h}\right] B_1(x/h) dx. $$

Here we used that $u(0) = u(h) = 0$ to kill of the first terms.

Now take $u(x) = f'(x) - {f(h) - f(0) \over h}$ and $v(x) = h {B_2(x/h) \over 2}$ (so that $v'(x) = B_1(x/h)$):

$$\begin{align*} - h \int_0^h \left[f(x) - {f(h) - f(0) \over h}\right] B_1(x/h) dx &= -{h^2 \over 2}\left[\left(f'(h) - {f(h) - f(0) \over h}\right) B_1(1) - \left(f'(0) - {f(h) - f(0) \over h}\right) B_1(-1)\right] \cr & \qquad +{h^2 \over 2} \int_0^h f''(x) B_2(x/h) dx \cr & = - {f'(h) - f'(0) \over 12} h^2 +{h^2 \over 2} \int_0^h f''(x) B_2(x/h) dx \end{align*} $$

By the same logic, we have in each panel

$$ \int_{x_{k-1}}^{x_k} \left[f(x) - \left(f(x_{k-1}) + x {f(x_k) - f(x_{k-1}) \over h}\right) \right] dx = - {f'(x_k) - f'(x_{k-1}) \over 12} h^2 +{h^2 \over 2} \int_{x_{k-1}}^{x_k} f''(x) B_2((x-x_{k-1})/h) dx $$

Summing over every panel, and using the telescoping sum, gives us

$$ \begin{align*} \int_0^1 f(x) dx - h\left({f(x_0) \over 2} + \sum_{k=1}^{n-1} f(x_k) + {f(x_n) \over 2} \right) = -{h^2 \over 12} \sum_{k=1}^n \left(f'(x_k) - f'(x_{k-1}) - \int_{x_{k-1}}^{x_k} f''(x) B_2((x-x_{k-1})/h) dx\right) \cr -{h^2 \over 12} (f'(1) - f'(0)) + {h^2 \over 12} \sum_{k=1}^n \int_{x_{k-1}}^{x_k} f''(x) B_2((x-x_{k-1})/h) dx \end{align*} $$

But on each panel we have

$$\left|\int_{x_{k-1}}^{x_k} f''(x) B_2((x-x_{k-1})/h) dx\right| \leq \int_{x_{k-1}}^{x_k} |f''(x)| |B_2((x-x_{k-1})/h)| dx \leq M_2 {h \over 6}$$

where $M_2 = \sup_{x \in [0,1]} |f''(x)|$ and we used the fact that $|B_2(x)| \leq {1 \over 6}$ (Exercise). Thus we have

$$\left|\sum_{k=1}^n \int_{x_{k-1}}^{x_k} f''(x) B_2((x-x_{k-1})/h) dx\right| \leq {M_2 \over 6}$$

and hence we have

$$ \int_0^1 f(x) dx - h\left({f(x_0) \over 2} + \sum_{k=1}^{n-1} f(x_k) + {f(x_n) \over 2} \right) = -{1 \over 12 n^2} (f'(1) - f'(0)) + O(n^{-2}) = O(n^{-2}) $$

which proves our second observation.