NA Simpson rule ($n=2$)

Section 8.5 Simpson rule ( $n=2$ )

Take

$x_0=a\text{,}$

$x_1=(a+b)/2$ and

$x_2=b\text{.}$ Then

$\begin{equation*} w_2 = \int_a^b\frac{x-\frac{a+b}{2}}{b-\frac{a+b}{2}}\frac{x-a}{b-a}dx = \frac{b-a}{6} = w_0 \end{equation*}$

and

$\begin{equation*} w_1 = \int_a^b\frac{x-a}{\frac{a+b}{2}-a}\frac{x-b}{\frac{a+b}{2}-b}dx = 2\frac{b-a}{3}. \end{equation*}$

Hence the method formula is

$\begin{equation*} Q_2(f,[a,b]) = \frac{b-a}{6}\left[f(a)+4f(\frac{a+b}{2})+f(b)\right]. \end{equation*}$

Simpson rule error. From the Newton-Cotes error formula (recall that

$n$ is even and so we need to use the remainder for

$n=3$ ), the error bound is given by

$\begin{equation*} \frac{(b-a)^5}{24}\displaystyle\max_{[a,b]}|f^{(4)}(x)|\int_0^2s(s-1)(s-2)(s-3)ds. \end{equation*}$

Hence,

$\begin{equation*} \left|\int_a^bf(x)dx-Q_2(f,[a,b])\right|\leq\frac{(b-a)^5}{90}\displaystyle\max_{[a,b]}|f^{(4)}(x)|. \end{equation*}$

After the usual decomposition procedure, we get finally that

$\begin{equation*} \left|\int_a^bf(x)dx-Q_{2,h}(f,[a,b])\right|\leq h^4\frac{b-a}{90}\max_{[a,b]}|f^{(4)}(x)|. \end{equation*}$

Convergence Speed. The formula above shows that the truncation error of the Riemann sum is

$O(h^4)\text{.}$ Below we visualize this scaling law by evaluating the error in the computation of

$\int_0^2\sin(x)dx$ with

$h=2^{-j}, 1\leq j\leq10\text{.}$ Due to the high convergence speed of this method, the round-off error shows already before

$h^{-3}$ but the error, before the round-off error kicks-off, goes down to about

$10^{-15}$ (for a much larger

$h$ than it was possible with the previous methods).

xxxxxxxxxx
 
from numpy import linspace,array,frompyfunc,sin,zeros                                                                                                                                         
from matplotlib import pyplot as plt 
n=15; a=0.; b=2.; N=1
s = zeros(n); err = zeros(n)
​
for i in range(n):
    h = (b-a)/N
    x = linspace(a,b,N+1)                                             # these are the x_i
    sim = ( sin(x[:-1]) + 4*sin((x[:-1]+x[1:])/2) + sin(x[1:]) ) / 6; # Simpson rule term
    s[i] = sum(h*sim)                                                 # this is the integral
    err[i] = abs(1-cos(2)-s[i])                                       # and this is the actual error
    N *= 2
​
I = range(n)
h = array([(b-a)/2**i for i in I])
plt.clf()
plt.loglog(h,err[I],label=r'$y=error(h)$')
plt.loglog(h,h**4,label=r'$y=h^4$')
plt.xlabel('log h') 
plt.ylabel('log y') 
plt.legend()
plt.show()

Section 8.5 Simpson rule (n=2n=2)

Section 8.5 Simpson rule ( $n=2$ )