Simpson's Rule

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

Definition

Simpson's rule uses a quadratic polynomial on each subinterval of a partition to approximate the function $f(x)$ and to compute the definite integral. This is an improvement over the trapezoid rule which approximates $f(x)$ by a straight line on each subinterval of a partition.

The formula for Simpson's rule is

$$ S_N(f) = \frac{\Delta x}{3} \sum_{i=1}^{N/2} \left( f(x_{2i-2}) + 4 f(x_{2i-1}) + f(x_{2i}) \right) $$

where $N$ is an even number of subintervals of $[a,b]$, $\Delta x = (b - a)/N$ and $x_i = a + i \Delta x$.

Error Formula

Theorem Let $S_N(f)$ denote Simpson's rule

$$ S_N(f) = \frac{\Delta x}{3} \sum_{i=1}^{N/2} \left( f(x_{2i-2}) + 4 f(x_{2i-1}) + f(x_{2i}) \right) $$

where $N$ is an even number of subintervals of $[a,b]$, $\Delta x = (b - a)/N$ and $x_i = a + i \Delta x$. The error bound is

$$ E_N^S(f) = \left| \ \int_a^b f(x) \, dx - S_N(f) \ \right| \leq \frac{(b-a)^5}{180N^4} K_4 $$

where $\left| \ f^{(4)}(x) \ \right| \leq K_4$ for all $x \in [a,b]$.

Implementation

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

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

    Simpson's rule approximates the integral \int_a^b f(x) dx by the sum:
    (dx/3) \sum_{k=1}^{N/2} (f(x_{2i-2} + 4f(x_{2i-1}) + f(x_{2i}))
    where x_i = a + i*dx and dx = (b - a)/N.

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

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

    Examples
    --------
    >>> simps(lambda x : 3*x**2,0,1,10)
    1.0
    '''
    if N % 2 == 1:
        raise ValueError("N must be an even integer.")
    dx = (b-a)/N
    x = np.linspace(a,b,N+1)
    y = f(x)
    S = dx/3 * np.sum(y[0:-1:2] + 4*y[1::2] + y[2::2])
    return S

Let's test our function on integrals for which we know the exact value. For example, we know

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

simps(lambda x : 3*x**2,0,1,10)
1.0

Test our function again with the integral

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

simps(np.sin,0,np.pi/2,100)
1.0000000003382361

scipy.integrate.simps

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.simps computes the approximation of a definite by Simpson's 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.simps returns the approximation of the integral using Simpson's rule.

Examples

Approximate ln(2)

Find a value $N$ which guarantees that Simpson's rule approximation $S_N(f)$ of the integral

$$ \int_1^2 \frac{1}{x} dx $$

satisfies $E_N^S(f) \leq 0.0001$.

Compute

$$ f^{(4)}(x) = \frac{24}{x^5} $$

therefore $\left| \, f^{(4)}(x) \, \right| \leq 24$ for all $x \in [1,2]$ and so

$$ \frac{1}{180N^4} 24 \leq 0.0001 \ \Rightarrow \ \frac{20000}{15N^4} \leq 1 \ \Rightarrow \ \left( \frac{20000}{15} \right)^{1/4} \leq N $$

Compute

(20000/15)**0.25
6.042750794713537

Compute Simpson's rule with $N=8$ (the smallest even integer greater than 6.04)

approximation = simps(lambda x : 1/x,1,2,8)
print(approximation)
0.693154530655

We could also use the function scipy.integrate.simps to compute the exact same result

N = 8; a = 1; b = 2;
x = np.linspace(a,b,N+1)
y = 1/x
approximation = spi.simps(y,x)
print(approximation)
0.693154530655

Verify that $E_N^S(f) \leq 0.0001$

np.abs(np.log(2) - approximation) <= 0.0001
True

Exercises

Under construction