Trapezoidal Rule - Numerical Analysis

#%config InlineBackend.figure_format = 'svg'
from pylab import *

8.3.1Single interval¶

Consider a function $f : [a,b] \to \re$ and let us approximate this by a polynomial of degree one, say $p_1(x)$ . This can be obtained by interpolation

p_1(a) = f(a), \qquad p_1(b) = f(b)

(1)

Then the approximate integral is

I_1(f) = \int_a^b p_1(x) \ud x = \frac{b-a}{2}[f(a) + f(b)]

(2)

a, b = 0.0, 1.0
f = lambda x: 1.0 + x + 0.5* sin(pi*x)
x = linspace(a,b,100)
plot([a,a],[0,f(a)],'k--')
plot([b,b],[0,f(b)],'k--')
plot([a,b],[0,0],'bo-')
plot(x, f(x), 'r-', label='f(x)')
plot([a,b],[f(a),f(b)],'sk-',label='$p_1(x)$')
d = 0.05
text(a,-d,"a",ha='center',va='top')
text(b,-d,"b",ha='center',va='top')
text(a,f(a)+0.06,"f(a)",ha='center',va='bottom')
text(b,f(b)+0.06,"f(b)",ha='center',va='bottom')
axis([-0.1, 1.1, -0.3, 2.2])
legend(), xlabel('x');

Note that this is the area under the straight line curve $p_1(x)$ which has the shape of a (rotated) trapezoid. To study the error in this approximation, we start with the error in interpolation (30)

f(x) - p_1(x) = (x-a)(x-b) f[a,b,x]

(3)

Then, the error in the integral is

\begin{aligned} E_1(f) &= \int_a^b [f(x) - p_1(x)] \ud x \\ &= \int_a^b (x-a) (x-b) f[a,b,x] \ud x \end{aligned}

(4)

Using integral mean value theorem

E_1(f) = f[a,b,\xi] \int_a^b (x-a)(x-b) \ud x, \qquad \textrm{for some $\xi \in [a,b]$}

(6)

Finally, using the properties of divided differences (65)

\begin{aligned} E_1(f) &= \left[ \half f''(\eta) \right] \left[ -\frac{1}{6} (b-a)^3 \right], \qquad \textrm{for some $\eta \in [a,b]$} \\ &= -\frac{1}{12} (b-a)^3 f''(\eta) \end{aligned}

(7)

The error will be small only if $b-a$ is small.

8.3.2Composite trapezoidal rule¶

Divide $[a,b]$ into $n$ equal intervals by the partition

a = x_0 < x_1 < \ldots < x_n = b \qquad \textrm{with spacing} \qquad h = \frac{b-a}{n}

(8)

and

x_j = a + j h, \qquad j=0,1,2,\ldots, n

(9)

figure(figsize=(8,5))
xp = linspace(0.0,1.0,10)
x  = linspace(0,1,100)
f  = lambda x: 2 + sin(2*pi*x)*cos(2*pi*x)
plot([0,1],[0,0],'k+-')
plot(xp,0*xp,'o')
plot(x,f(x),'r-',label='f(x)')
plot(xp,f(xp),'sk-')
for xx in xp:
    plot([xx,xx],[0,f(xx)],'b--')
axis([-0.1,1.1,-0.3,2.7])
d = 0.1
text(xp[0],-d,'$x_0$',va='top',ha='center')
text(xp[1],-d,'$x_1$',va='top',ha='center')
text(xp[-2],-d,'$x_{n-1}$',va='top',ha='center')
text(xp[-1],-d,'$x_n$',va='top',ha='center')
text(0,d,'$a$',ha='center',backgroundcolor='white')
text(1,d,'$b$',ha='center',backgroundcolor='white')
legend(), xlabel('x');

Let us use the Trapezoidal rule in each interval $[x_{j-1},x_j]$

I(f) = \int_a^b f(x) \ud x = \sum_{j=1}^n \int_{x_{j-1}}^{x_j} f(x) \ud x

(10)

From the previous section, we know that

\int_{x_{j-1}}^{x_j} f(x) \ud x = \frac{h}{2}[f(x_{j-1}) + f(x_j)] - \frac{h^3}{12}f''(\eta_j), \qquad \eta_j \in [x_{j-1}, x_j]

(11)

Hence

\begin{aligned} I(f) &= \clr{red}{ \sum_{j=1}^n \frac{h}{2}[f(x_{j-1}) + f(x_j)] } - \sum_{j=1}^n \frac{h^3}{12}f''(\eta_j) \\ &= \clr{red}{ I_n(f) } - \sum_{j=1}^n \frac{h^3}{12}f''(\eta_j) \end{aligned}

(12)

where

I_n(f) = h [ \shalf f_0 + f_1 + f_2 + \ldots + f_{n-1} + \shalf f_n]

(13)

is the composite Trapezoidal rule. The error is

E_n(f) = I(f) - I_n(f) = - \frac{h^3}{12} \sum_{j=1}^n f''(\eta_j) = - \frac{n h^3}{12} \underbrace{\frac{1}{n} \sum_{j=1}^n f''(\eta_j)}_{M_n}

(14)

where

\min_{x \in [a,b]} f''(x) \le \min_{1 \le j \le n} f''(\eta_j) \le M_n \le \max_{1 \le j \le n} f''(\eta_j) \le \max_{x \in [a,b]} f''(x)

(15)

If $f''(x)$ is continuous, then it attains all values between the minimum and maximum.

\implies \exists \eta \in [a,b] \qquad \textrm{such that} \qquad M_n = f''(\eta)

(16)

Hence, after using $nh=b-a$ , the error is given by

E_n(f) = - \frac{1}{12}(b-a) h^2 f''(\eta), \qquad \textrm{for some $\eta \in [a,b]$}

(17)

The error goes to zero as $n \to \infty$ since $h \to 0$ .

Example 1

The next functions compute integral with composite trapezoidal rule.

# n intervals, n+1 points
def trapz1(a,b,n,f):
    h = (b-a)/n
    x = linspace(a,b,n+1)
    y = f(x)
    res = h * (sum(y[1:n]) + 0.5*(y[0] + y[n]))
    return res

\begin{align} f(x) &= x, \qquad x \in [0,1] \\ \int_0^1 f(x) \ud x &= \half \end{align}

(20)

f = lambda x: x
print("Integral = ", trapz1(0.0,1.0,10,f))

Integral =  0.5

\begin{align} f(x) &= x^2, \qquad x \in [0,1] \\ \int_0^1 f(x) \ud x &= \frac{1}{3} \end{align}

(21)

f = lambda x: x**2
print("Integral = ", trapz1(0.0,1.0,10,f))

Integral =  0.3350000000000001

\begin{align} f(x) &= \exp(x)\cos(x), \qquad x \in [0,\pi] \\ \int_0^1 f(x) \ud x &= -\frac{1}{2}(1+\exp(\pi)) \end{align}

(22)

f = lambda x: exp(x)*cos(x)
qe = -0.5*(1.0 + exp(pi)) # Exact integral

n,N = 4,10
e = zeros(N)
for i in range(N):
    e[i] = trapz1(0.0,pi,n,f) - qe
    if i > 0:
        print('%6d %24.14e %14.5e'%(n,e[i],e[i-1]/e[i]))
    else:
        print('%6d %24.14e'%(n,e[i]))
    n = 2*n

     4    -1.26567653098185e+00
     8    -3.11816113365945e-01    4.05905e+00
    16    -7.76577835071954e-02    4.01526e+00
    32    -1.93958006245669e-02    4.00385e+00
    64    -4.84778281250620e-03    4.00096e+00
   128    -1.21187271271594e-03    4.00024e+00
   256    -3.02963615787633e-04    4.00006e+00
   512    -7.57406187865683e-05    4.00002e+00
  1024    -1.89351368735657e-05    4.00000e+00
  2048    -4.73378310594796e-06    4.00000e+00

The last column shows $E_n(f)/E_{2n}(f)$ , which is close to 4, showing that the error is $\order{h^2}$ .

Remark 4

The scipy.integrate module provides functions that perform the Trapezoidal integration.

from scipy.integrate import trapezoid
f = lambda x: x**2
a,b,n = 0.0,1.0,11
x = linspace(a,b,n)
y = f(x)
print("Integral = ", trapezoid(y,x))

Integral =  0.33499999999999996

For uniformly space points we can just pass the spacing instead of x.

h = (b-a)/(n-1)
print("Integral = ", trapezoid(y,dx=h))

Integral =  0.3350000000000001

8.3.3Asymptotic error estimate¶

We can estimate the error as

|E_n(f)| \le \frac{1}{12}(b-a) h^2 \max_{x \in [a,b]}|f''(x)|

(23)

provided we can estimate the second derivative. A more useful error estimate can be derived as follows.

Using the error of the composite rule, we can compute the limit

\begin{aligned} \lim_{n \to \infty} \frac{E_n(f)}{h^2} &= - \frac{1}{12} \lim_{n \to \infty} h \sum_{j=1}^n f''(\eta_j), \qquad \eta_j \in [x_{j-1}, x_j] \\ &= -\frac{1}{12} \int_a^b f''(x) \ud x \\ &= -\frac{1}{12}[f'(b) - f'(a)] \end{aligned}

(24)

Hence for large $n$

E_n(f) \approx -\frac{h^2}{12}[f'(b) - f'(a)] =: \tilde E_n(f)

(25)

8.3.4Corrected trapezoidal rule¶

Using the asymptotic error estimate

\begin{gather} I(f) - I_n(f) \approx \tilde E_n(f) = -\frac{h^2}{12}[f'(b) - f'(a)] \\ \implies I(f) \approx I_n(f) -\frac{h^2}{12}[f'(b) - f'(a)] \end{gather}

(28)

we can define a corrected Trapezoidal rule as

\begin{align} C_n(f) &= I_n(f) -\frac{h^2}{12}[f'(b) - f'(a)] \\ &= h \left[ \shalf f_0 + f_1 + f_2 + \ldots + f_{n-1} + \shalf f_n \right] - \frac{h^2}{12}[f'(b) - f'(a)] \end{align}

(29)

Example 2

Consider the integral

I = \int_0^\pi \ee^x \cos x \ud x

(30)

The exact value is

I = -\half (1 + \ee^\pi) \approx -12.0703463164

(31)

The next function implements both trapezoid and the corrected rule.

# n = number of intervals
# There are n+1 points
def trapz2(a,b,n,f,df):
    h = (b-a)/n
    x = linspace(a,b,n+1)
    y = f(x)
    res = sum(y[1:n]) + 0.5*(y[0] + y[n])
    return h*res, h*res - (h**2/12)*(df(b) - df(a))

f  = lambda x: exp(x)*cos(x)
df = lambda x: exp(x)*(cos(x) - sin(x))
qe = -0.5*(1.0 + exp(pi)) # Exact integral

n,N = 4,10
e1,e2 = zeros(N),zeros(N)
for i in range(N):
    e1[i],e2[i] = trapz2(0.0,pi,n,f,df) - qe
    if i > 0:
        print('%6d %24.14e %10.5f %24.14e %10.5f' %
              (n,e1[i],e1[i-1]/e1[i],e2[i],e2[i-1]/e2[i]))
    else:
        print('%6d %24.14e %10.5f %24.14e %10.5f' %
              (n,e1[i],0,e2[i],0))
    n = 2*n

     4    -1.26567653098185e+00    0.00000    -2.47437900765224e-02    0.00000
     8    -3.11816113365945e-01    4.05905    -1.58292813961225e-03   15.63166
    16    -7.76577835071954e-02    4.01526    -9.94872006128134e-05   15.91087
    32    -1.93958006245669e-02    4.00385    -6.22654792081789e-06   15.97791
    64    -4.84778281250620e-03    4.00096    -3.89293344227326e-07   15.99449
   128    -1.21187271271594e-03    4.00024    -2.43329250082525e-08   15.99863
   256    -3.02963615787633e-04    4.00006    -1.52084034255040e-09   15.99966
   512    -7.57406187865683e-05    4.00002    -9.50493017626286e-11   16.00054
  1024    -1.89351368735657e-05    4.00000    -5.94013727095444e-12   16.00120
  2048    -4.73378310594796e-06    4.00000    -3.73034936274053e-13   15.92381

As seen in last column, the error in the corrected rule converges at a rate of 16, i.e., as $\order{h^4}$ .