Trapezoid Rule

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Definition

The definite integral of $f(x)$ is equal to the (net) area under the curve $y=f(x)$ over the interval $[a,b]$. Riemann sums approximate definite integrals by using sums of rectangles to approximate the area.

The trapezoid rule gives a better approximation of a definite integral by summing the areas of the trapezoids connecting the points

$$ (x_{i-1},0), (x_i,0), (x_{i-1},f(x_{i-1})), (x_i,f(x_i)) $$

for each subinterval $[x_{i-1},x_i]$ of a partition. Note that the area of each trapezoid is

$$ (x_i - x_{i-1}) f(x_{i-1}) + \frac{1}{2}(x_i - x_{i-1}) (f(x_i) - f(x_{i-1})) = \frac{1}{2}(f(x_i) + f(x_{i-1}))(x_i - x_{i-1}) $$

For example, we can use a single trapezoid to approximate:

$$ \int_0^1 e^{-x^2} \, dx $$

f = lambda x: np.exp(-x**2)
x = np.linspace(-0.5,1.5,100)
y = f(x)
plt.plot(x,y)

x0 = 0; x1 = 1;
y0 = f(x0); y1 = f(x1);
plt.fill_between([x0,x1],[y0,y1])

plt.xlim([-0.5,1.5]); plt.ylim([0,1.5]);
plt.show()

png

A = 0.5*(x1 - x0)*(y1 + y0)
print("Trapezoid area:",A)
Trapezoid area: 0.683939720586

For a partitition $x_0 = a < x_1 < \cdots < x_n = b$, the trapezoid rule formula simplifies to

$$ \int_a^b f(x) \, dx \approx \frac{1}{2} \sum_{i = 1}^n \left( f(x_i) + f(x_{i-1}) \right) ( x_i - x_{i-1} ) $$

In particular, given an integer $N$, the trapezoid rule for $N$ subintervals of $[a,b]$ of equal length is

$$ T_N(f) = \frac{\Delta x}{2} \sum_{i=1}^N (f(x_i) + f(x_{i-1})) $$

where $\Delta x = (b - a)/N$ is the length of the subintervals and $x_i = a + i \Delta x$. Notice that the trapezoid is the average of the left and right Riemann sums

$$ T_N(f) = \frac{\Delta x}{2} \sum_{i=1}^N (f(x_i) + f(x_{i-1})) = \frac{1}{2} \left( \sum_{i=1}^N f(x_i) \Delta x + \sum_{i=1}^N f(x_{i-1}) \Delta x \right) $$

For example, let's plot the trapezoids for $f(x)=\frac{1}{1 + x^2}$ on $[0,5]$ with $N=10$.

f = lambda x : 1/(1+x**2)
a = 0; b = 5; N = 10
n = 10 # Use n*N+1 points to plot the function smoothly

x = np.linspace(a,b,N+1)
y = f(x)

X = np.linspace(a,b,n*N+1)
Y = f(X)

plt.plot(X,Y)

for i in range(N):
    xs = [x[i],x[i],x[i+1],x[i+1]]
    ys = [0,f(x[i]),f(x[i+1]),0]
    plt.fill(xs,ys,'b',edgecolor='b',alpha=0.2)

plt.title('Trapezoid Rule, N = {}'.format(N))
plt.show()

png

Let's compute the sum of areas of the trapezoids:

y = f(x)
y_right = y[1:] # Right endpoints
y_left = y[:-1] # Left endpoints
dx = (b - a)/N
T = (dx/2) * np.sum(y_right + y_left)
print(T)
1.37310408123

We know the exact value

$$ \int_0^5 \frac{1}{1 + x^2} dx = \arctan(5) $$

and we can compare the trapezoid rule to the value

I = np.arctan(5)
print(I)
1.37340076695
print("Trapezoid Rule Error:",np.abs(I - T))
Trapezoid Rule Error: 0.000296685714906

Error Formula

Theorem. Let $T_N(f)$ denote the trapezoid rule

$$ T_N(f) = \frac{\Delta x}{2} \sum_{i=1}^N (f(x_i) + f(x_{i-1})) $$

where $\Delta x = (b-a)/N$ and $x_i = a + i \Delta x$. The error bound is

$$ E_N^T(f) = \left| \ \int_a^b f(x) \ dx - T_N(f) \ \right| \leq \frac{(b-a)^3}{12 N^2} K_2 $$

where $\left| \ f''(x) \, \right| \leq K_2$ for all $x \in [a,b]$.

Implementation

Let's write a function called trapz which takes input parameters $f$, $a$, $b$ and $N$ and returns the approximation $T_N(f)$. Furthermore, let's assign default value $N=50$.

def trapz(f,a,b,N=50):
    '''Approximate the integral of f(x) from a to b by the trapezoid rule.

    The trapezoid rule approximates the integral \int_a^b f(x) dx by the sum:
    (dx/2) \sum_{k=1}^N (f(x_k) + f(x_{k-1}))
    where x_k = a + k*dx and dx = (b - a)/N.

    Parameters
    ----------
    f : function
        Vectorized function of a single variable
    a , b : numbers
        Interval of integration [a,b]
    N : integer
        Number of subintervals of [a,b]

    Returns
    -------
    float
        Approximation of the integral of f(x) from a to b using the
        trapezoid rule with N subintervals of equal length.

    Examples
    --------
    >>> trapz(np.sin,a=0,b=np.pi/2,N=1000)
    0.99999979438323316
    '''
    x = np.linspace(a,b,N+1)
    y = f(x)
    y_right = y[1:] # Right endpoints
    y_left = y[:-1] # Left endpoints
    dx = (b - a)/N
    T = (dx/2) * np.sum(y_right + y_left)
    return T

Let's test our function on an integral where we know the answer

$$ \int_0^{\pi/2} \sin x \ dx = 1 $$

trapz(np.sin,0,np.pi/2,1000)
0.99999979438323316

Let's test our function again:

$$ \int_0^1 3 x^2 \ dx = 1 $$

trapz(lambda x : 3*x**2,0,1,10000)
1.0000000050000002

And once more:

$$ \int_0^1 x \ dx = \frac{1}{2} $$

trapz(lambda x : x,0,1,1)
0.5

scipy.integrate.trapz

The SciPy subpackage scipy.integrate contains several functions for approximating definite integrals and numerically solving differential equations. Let's import the subpackage under the name spi.

import scipy.integrate as spi

The function scipy.integrate.trapz computes the approximation of a definite by the trapezoid rule. Consulting the documentation, we see that all we need to do it supply arrays of $x$ and $y$ values for the integrand and scipy.integrate.trapz returns the approximation of the integral using the trapezoid rule. The number of points we give to scipy.integrate.trapz is up to us but we have to remember that more points gives a better approximation but it takes more time to compute!

Examples

Approximate ln(2)

Find a value $N$ which guarantees that the trapezoid rule approximation $T_N(f)$ of the integral

$$ \int_1^2 \frac{1}{x} \, dx = \ln(2) $$

satisfies $E_N^T(f) \leq 10^{-8}$.

For $f(x) = \frac{1}{x}$, we compute $f''(x) = \frac{2}{x^3} \leq 2$ for all $x \in [1,2]$ therefore the error formula implies

$$ \left| \int_1^2 \frac{1}{x} \, dx - T_N(f) \right| \leq \frac{2}{12N^2} $$

Then $E_N^T \leq 10^{-8}$ is guaranteed if $\frac{1}{6N^2} \leq 10^{-8}$ which implies

$$ \frac{10^4}{\sqrt{6}} \leq N $$

10**4/np.sqrt(6)
4082.4829046386303

We need 4083 subintervals to guarantee $E_N^T(f) \leq 10^{-8}$. Compute the approximation using our own implementation of the trapezoid rule:

approximation = trapz(lambda x : 1/x,1,2,4083)
print(approximation)
0.693147184309

We could also use scipy.integrate.trapz to get the exact same result:

N = 4083
x = np.linspace(1,2,N+1)
y = 1/x
approximation = spi.trapz(y,x)
print(approximation)
0.693147184309

Let's verify that this is within $10^{-6}$

np.abs(approximation - np.log(2)) < 10**(-8)
True

Success! However, a natural question arises: what is the actual smallest $N$ such that the trapezoid rule gives the estimate of $\ln (2)$ to within $10^{-8}$?

for n in range(1,4083):
    approx = trapz(lambda x : 1/x,1,2,n)
    if np.abs(approx - np.log(2)) < 10e-8:
        print("Accuracy achieved at N =",n)
        break
Accuracy achieved at N = 791

Fresnel Integral

Fresnel integrals are examples of nonelementary integrals: antiderivatives which cannot be written in terms of elementary functions. There are two types of Fresnel integrals:

$$ S(t) = \int_0^t \sin(x^2) dx \ \ \text{and} \ \ C(t) = \int_0^t \cos(x^2) dx $$

Use the trapezoid rule to approximate the Fresnel integral

$$ S(1) = \int_0^1 \sin(x^2) dx $$

such that the error is less than $10^{-5}$.

Compute the derivatives of the integrand

$$ f(x) = \sin(x^2) \ \ , \ \ f'(x) = 2x\cos(x^2) $$

$$ f''(x) = 2\cos(x^2) - 4x^2\sin(x^2) \ \ , \ \ f'''(x) = -12x\sin(x^2) - 8x^3\cos(x^2) $$

Since $f'''(x) \leq 0$ for $x \in [0,1]$, we see that $f''(x)$ is decreasing on $[0,1]$. Values of $f''(x)$ at the endpoints of the interval are

x = 0
2*np.cos(x**2) - 4*x**2*np.sin(x**2)
2.0
x = 1
2*np.cos(x**2) - 4*x**2*np.sin(x**2)
-2.2852793274953065

Therefore $\left| \, f''(x) \, \right| \leq 2.2852793274953065$ for $x \in [0,1]$. Use the error bound formula to find a good choice for $N$

$$ \frac{(b-a)^3}{12 N^2} K_2 \leq 10^{-5} \Rightarrow \sqrt{\frac{10^5(2.2852793274953065)}{12}} \leq N $$

np.sqrt(10**5 * 2.2852793274953065 / 12)
137.99997969490511

Let's compute the integral using the trapezoid rule with $N=138$ subintervals

x = np.linspace(0,1,139)
y = np.sin(x**2)
I = spi.trapz(y,x)
print(I)
0.310273030322

Therefore the Fresnel integral $S(1)$ is approximately

$$ S(1) = \int_0^1 \sin(x^2) \, dx \approx 0.310273030322 $$

with error less than $10^{-5}$.

Logarithmic Integral

The Eulerian logarithmic integral is another nonelementary integral

$$ \mathrm{Li}(t) = \int_2^t \frac{1}{\ln x} dx $$

Let's compute Li(10) such that the error is less than 0.01. Compute derivatives of the integrand

$$ f(x) = \frac{1}{\ln x} \ \ , \ \ f'(x) = -\frac{1}{x(\ln x)^2} \ \ , \ \ f''(x) = \frac{\ln x + 2 }{x^2\ln x} $$

Plot $f''(x)$ on the interval $[2,10]$

x = np.linspace(2,10,100)
y = (np.log(x) + 2) / (x**2 * np.log(x))
plt.plot(x,y)
plt.show()

png

Clearly $f''(x)$ is decreasing and is approaching 0 therefore

$$ \left| \, f''(x) \, \right| \leq \frac{\ln (2) + 2}{4 \ln (2)} $$

for $x \in [2,\infty)$.

(np.log(2) + 2)/(4*np.log(2))
0.9713475204444818

Use the error formula:

$$ \frac{(b-a)^3}{12 N^2} K_2 \leq 0.01 \Rightarrow \frac{8^3}{12 N^2} 0.9713475204444818 \leq 0.01 \Rightarrow \sqrt{ \frac{8^3}{12(0.01)} 0.9713475204444818} \leq N $$

np.sqrt(8**3 * 0.9713475204444818 / 0.12)
64.377139476912063

Compute the trapzoid rule with $N=65$

x = np.linspace(2,10,66)
y = 1/np.log(x)
I = spi.trapz(y,x)
print(I)
5.12172369518

Therefore the Eulerian logarithmic integral is

$$ \mathrm{Li}(10) = \int_2^{10} \frac{1}{\ln x} dx \approx 5.12172369518 $$

such that the error is less than 0.01.

Exercises

Under construction