Simpson's rule can be derived by approximating the function $f(x)$ on the interval $[a,b]$ by connecting three consecutive points of the partition with a parabola, then treating the next set of three points in the partition in the same way, and continuing until you reach the end of the interval $[a,b]$, for a partition of $n+1$ equidistant points on the interval with mesh size $h=\frac{b-a}{n}$. So we assume that $n$ is even, say $n=2m$. The figure below visualizes the situation for $n=2$.

The parabola $P_2(x)$ can be found via Lagrange interpolation and the formula is

$$\begin{aligned}P_2(x) &= \;\;\;\,\,\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}f(x_0)\\ \\ &\phantom{=}{}+\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}f(x_1) \\ \\ &\phantom{=} {}+\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}f(x_2)\end{aligned}$$ with $x_0=a, x_1=a+h, x_2=a+2h$. In other words: $$\begin{aligned}P_2(x)&=\frac{(x-a-h)(x-a-2h)}{2h^2}f(x_0)\\ \\ &\phantom{=}-\frac{(x-a)(x-a-2h)}{h^2}f(x_1)\\ \\ &\phantom{=}+\frac{(x-a)(x-a-h)}{2h^2}f(x_2)\end{aligned}$$ Then we get with the substitution rule for $y=\frac{x-a}{h}$ that $$\begin{aligned}\int_a^bP_2(x)\,\dd x &= f(x_0)\int_0^2\frac{(y-1)(y-2)}{2}dy-f(x_1)\int_0^2y(y-2)\,dy\\ &\phantom{=}+f(x_2)\int_0^2y(y-1)\,dy\\ \\ &= \frac{h}{3}\bigl(f(x_0)+4f(x_1)+f(x_2)\bigr)\end{aligned}$$

We can now look at a refined partition. The points of the partition are $x_k=a+k h$ with $k=0,1,\ldots 2m$ and we can divide them into groups of 3: $\{x_{2k-2},x_{2k-1},x_{2k}\}$ for $k=1,\ldots m$. At each sub-interval $[x_{2k-2},x_{2k-1}, x_{2k}]$ we approximate the integral

$$\int_{x_{2k-2}}^{x_{2k}}f(x)\,\dd x$$ with the area $I_k$ under the parabola through the points $\bigl((x_{2k-2},f(x_{2k-2}\bigr))$, $\bigl(x_{2k-1},f(x_{2k-1})\bigr)$ and $\bigl(x_{2k},f(x_{2k})\bigr)$. See the figure below to get an impression of how this works.

Above we have already derived the approximation formula for each sub-interval:

$$I_k=\frac{h}{3}\bigl(f(x_{2k-2})+4f(x_{2k-1})+f(x_{2k})\bigr)$$ Summing over all sub-intervals we get $$\begin{aligned}I&=\sum_{k=1}^mI_k=\sum_{k=1}^m \frac{h}{3}\bigl(f(x_{2k-2})+4f(x_{2k-1})+f(x_{2k})\bigr) \\ \\ &= \frac{h}{3}\bigl(f(a)+f(b) +2\sum_{k=1}^{m-1} f(a+2kh)+4\sum_{k=1}^{m} f(a+(2k-1)h)\bigr)\end{aligned}$$ We have thus found the following result.

In case we allow points halfway through a sub-interval, we can also rewrite Simpson's rule, as

$$\int_a^bf(x)\,\dd x \approx \frac{h}{6}\Bigl(f(a)+f(b)+2\sum_{k=1}^{n-1} f(a+kh)+4\sum_{k=1}^{n} f\bigl(a+(k-\tfrac{1}{2})h\bigr)\Bigr)$$ Then you can better see that

Simpson's formula ${}=\tfrac{1}{3}\cdot{}$ trapezoid rule ${}+\tfrac{2}{3}\cdot{}$ midpoint-Riemann sum

We already know that the truncation error of the midpoint Riemann sum and the trunctation error of the trapezoidal rule are quadratic in mesh size $h$; so Simpson's rule is at least quadratic in $h$, but in reality Simpson's rule is an even better approximation.

---

Programming task

Write a function SimpsonRule(f,a,b,n) that implements Simpson's rule.

Apply your function with $n=100$ in the following two cases:

• $\displaystyle \int_0^1\frac{4}{x^2+1}\,\dd x=\pi$
  
  
• $\displaystyle \int_0^{\pi}\sin(x)\,\dd x=2$