Accuracy of Simpson's Rule

In [1]:
import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as pt

In [2]:
simpson_weights = np.array([1/6, 4/6, 1/6])
simpson_weights

Out[2]:
array([ 0.16666667,  0.66666667,  0.16666667])


Here is a function and its definite integral:

$$\text{int_f}(x)=\int_0^x f(\xi)d\xi$$

In [1]:
alpha = 4

def f(x):
return np.cos(alpha*x)
def int_f(x):
return 1/alpha*(np.sin(alpha*x)-np.sin(alpha*0))


Plotted:

In [53]:
plot_x = np.linspace(0, 1, 200)

pt.plot(plot_x, f(plot_x), label="f")
pt.fill_between(plot_x, 0*plot_x, f(plot_x),alpha=0.3)
pt.plot(plot_x, int_f(plot_x), label="$\int f$")
pt.grid()
pt.legend(loc="best")

Out[53]:
<matplotlib.legend.Legend at 0x7fc5410f2128>


This here plots the function, the interpolant, and the area under the interpolant:

In [61]:
# fix nodes
h = 1
x = np.linspace(0, h, 3)

# find interpolant
V = np.array([
1+0*x,
x,
x**2
]).T
a, b, c = la.solve(V, f(x))
interpolant = a + b*plot_x + c*plot_x**2

# plot
pt.plot(plot_x, f(plot_x), label="f")
pt.plot(plot_x, interpolant, label="Interpolant")
pt.fill_between(plot_x, 0*plot_x, interpolant, alpha=0.3, color="green")
pt.plot(x, f(x), "og")
pt.grid()
pt.legend(loc="best")

Out[61]:
<matplotlib.legend.Legend at 0x7fc540e509b0>

In [55]:
true_val = int_f(h)
error = (f(x).dot(simpson_weights) * h - true_val)/true_val
print(error)

0.161228524153


• Compare the error for $h=1,0.5,0.25$. What order of accuracy do you observe?