# 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