Next Step from Simpson's Rule: Third Degree Approximation

Simpson's Rule, which is a extension of Trapezoidal Rule, which itself is a extension of Riemann Sum, is a approximation of true integral using second degree polynomial (rather than first degree polynomial of Trapezoidal and zeroth degree of Riemann Sum).
Simpson uses the fact that
$$p(x) = Ax^2 + Bx + C, \\
\int_{a}^b p(x) dx = \frac{b-a}{6}(p(a)+4p(\frac{a+b}{2})+p(b))$$ and extends to approximate that
$$\int_{a}^bf(x) dx \approx \frac{b-a}{3n}[f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+\dots +4f(x_{n-1})+f(x_n)]$$ Now then, surely if second order polynomial yielded a closer approximation than first, then third order polynomial would be better approximation than second as well.

What does it mean to use Third-degree Polynomial?

For zeroth degree polynomial, we pick a single point in the interval to represent the curve. We cannot know about the slope or any concavity with a single point, and so a zeroth degree polynomial (a constant) is needed.
This is what Riemann Sum uses to approximate an integral.

Zeroth Degree Polynomial approximating the curve in an interval using one point.

For first degree polynomial, we simply picked two points and connected it with a straight line. Between two points, we cannot deduce any concavity and can only find the slope, meaning first-degree polynomial would be needed.
This is what Trapezoidal Sum uses to approximate an integral.

First Degree Polynomial approximating the curve in an interval using two points.

For second degree polynomial, we picked three points and fitted a curve. Now for three points, the point can show concavity, indicated that we need second degree polynomial to fit all three points.
This is what Simpson's Rule uses to approximate an integral.

Second Degree Polynomial approximating the curve in an interval using three points.

Now then, Third-degree polynomial would mean picking 4 points, then fitting a curve between those points. Since 4 points can indicate at most 1 inflection, a third-degree polynomial would be needed to accurately connect them.
This is what we will use to approximate an integral.

Third Degree Polynomial approximating the curve in an interval using four points.

Now we simply must show how four points picked on a third-degree polynomial can indicate the area under the interval.
First let's see how Thomas Simpson was able to do this.

Simpson's Rule

We first know that
$$\begin{align*}
p(x) &= Ax^2 + Bx + C\\
\int_{a}^bp(x)dx &= \left[ A\frac{x^3}{3} + B\frac{x^2}{2} + Cx \right]_a^b \\
&= A\frac{b^3-a^3}{3} + B\frac{b^2-a^2}{2} + C(b-a) \\
& = \frac{b-a}{6}[ 2A(b^2+ba+a^2) + 3(b+a)  + 6C ] \\
&\text{and for this step my Textbook states}\\
&\text{''by expansion and collection of terms, the expression inside the bracket becomes''} \\
&= \frac{b-a}{6}[ (Aa^2+Ba+C) + 4(A(\frac{b+a}{2})^2 + B(\frac{b+a}{2}) + C) + (Ab^2+Bb+C) ] \\
\therefore&= \frac{b-a}{6}[p(a)+4p(\frac{b+a}{2})+p(b)]
\end{align*}$$ So, though textbook was not very clear about how to get to last step, we have a general idea on how to begin.

Area Under Third Degree Polynomial

Let's first begin with integrating
$$\begin{align*}
q(x) &= Ax^3+Bx^2+Cx+D \\
\int_a^bq(x)dx &= \left[ A\frac{x^4}{4} + B\frac{x^3}{3} + C\frac{x^2}{2} +Dx \right]_a^b \\
&= A\frac{b^4-a^4}{4} + B\frac{b^3-a^3}{3} + C\frac{b^2-a^2}{2} + D(b-a) \\
&\text{now to turn this expression in terms of $q(x)$,}\\
&\text{we first need to reduce the fourth order terms to third orders} \\
&= \frac{b-a}{4!}\left(6A(b^3+b^2a+ba^2+a^3) + 8B(b^2+ba+a^2)+12C(b+a) + 24D\right)
\end{align*}$$ Now we arrive at the crucial step of expressing this in terms of $q(x)$. Beware: if we start mindlessly "expanding and collecting the terms," we will simply arrive at Simpson's Rule again.
Let's use the fact that we are trying to pick four distinct points and use those to express the integral. What four points should we pick between a and b?
Ideally, we want them to be evenly spaced throughout the interval, so we can choose
$$a,\frac{b+2a}{3},\frac{2b+a}{3},b$$ And so, we want to express the expression above as
$$\frac{b-a}{4!}\left( Kq(a) + Pq(\frac{b+2a}{3}) + Pq(\frac{2b+a}{3}) + Kq(b) \right)$$ where $K$ and $P$ are some real number 'weights,' just like 1 and 4 in Simpson's Rule.
The terms $q(a)$ and $q(b)$ has the same weights and the $q(\frac{b+2a}{3})$ and $q(\frac{2b+a}{3})$ also has the same weights. This is because if the curve was to be flipped horizontally, we still expect the formula to work and produce the same results, and since point at $b$ will become $a$ and so on and so on, we can deduce that these terms will have the same weights.
Then, what will the value for $K$ and $P$ be?

Let's first expand what the values for each points would be. We will do this in a tabular manner:
$$\begin{align*}
&q(b) &&q(\frac{b+2a}{3}) &&q(\frac{2b+a}{3}) &&q(a) \\
&Ab^3 &&\frac{A}{27}(b^3+6b^2a + 12ba^2 + 8a^3) &&\frac{A}{27}(8b^3 + 12b^2a + 8ba^2 + a^3) &&Aa^3\\
&Bb^2 &&\frac{B}{9}(b^2+4ba + 4a^2) &&\frac{B}{9}(4b^2+4ba+a^2) && Ba^2 \\
&Cb &&\frac{C}{3}(b+2a) && \frac{C}{3}(2b+a) && Ca \\
&D&& D &&D && D
\end{align*}$$
Let's add the terms for $q(\frac{b+2a}{3}) $ and $q(\frac{2b+a}{3})$ and find that
$$\begin{align*}
& \frac{A}{27}(9b^3 + 18b^2a + 18ba^2 + 9a^3) &=& \frac{A}{3}(b^3 + 2b^2a + 2ba^2 + a^3) \\
&\frac{B}{9}(5b^2+8ba + 5a^2) &=& \frac{B}{9}(5b^2+8ba + 5a^2) \\
&\frac{C}{3}(3b+3a) &=& C(b+a)\\
&2D &=& 2D \end{align*}$$ We notice that the LCM of denominators is 9, so let us set
$$P = 9L$$ So we find that
$$\begin{align*}
K(q(b)+q(a)) &= AK(b^3+a^3) + BK(b^2+a^2) + CK(b+a) + DK\\
P(q(\frac{b+2a}{3})+q(\frac{2b+a}{3})) &= 3AL(b^3 + 2b^2a + 2ba^2 + a^3) + BL(5b^2+8ba + 5a^2) + 9CL(b+a) + 18DL
\end{align*}$$
We know that
$$K(q(b)+q(a)) + P(q(\frac{b+2a}{3})+q(\frac{2b+a}{3})) = 6A(b^3+b^2a+ba^2+a^3) + 8B(b^2+ba+a^2)+12C(b+a) + 24D$$
Let's solve for $K$ and $L$ by using the all the terms that carry $A$.
$$\begin{align*}
\therefore 6A(b^3+b^2a+ba^2+a^3) &= AK(b^3+a^3) + 3AL(b^3 + 2b^2a + 2ba^2 + a^3)\\
6b^3 + 6b^2a + 6ba^2 + 6a^3 &= (K+3L)b^3 + 6Lb^2a + 6Lba^2 + (K+3L)a^3\\
\end{align*}$$ From this, it is easy to see that $K=3$ and $L = 1$.

We can solve for terms carrying $B$ or $C$ or even $D$ and reach the same result:
$$\begin{align*}
\therefore 8B(b^2+ba+a^2) &= BK(b^2+a^2) + BL(5b^2+8ba + 5a^2) \\
8b^2 + 8ba + 8a^2 &= (K+5L)b^2 + 8Lba + (K+5L)a^2
\end{align*}$$ Again we find $K=3$ and $L=1$.

$$\begin{align*}
\therefore 12C(b+a) &= CK(b+a) + 9CL(b+a) \\
12b + 12a &= (K+9L)b + (K+9L)a
\end{align*}$$ Same result here.

And finally:
$$\begin{align*}
\therefore 24D &= DK + 18DL \\
24 &= K+18L \end{align*}$$
And because $P = 9L$, we find that $P = 9$.
Therefore we discover that
$$\begin{align*}
q(x) &= Ax^3+Bx^2+Cx+D\\
\int_{a}^bq(x)dx &= \frac{b-a}{24}\left( 3q(a) +9q(\frac{2a+b}{3}) + 9q(\frac{a+2b}{3}) + 3q(b) \right) \\
&= \frac{b-a}{8}\left(q(a) + 3q(\frac{2a+b}{3}) + 3q(\frac{a+2b}{3}) + q(b) \right) \\
&\blacksquare
\end{align*}$$

Using this fact, we can build our approximation for integral of fuction $f$.

Third Degree Approximation

Notice that each interval has 3 increments. This means that if you wish to use $n$ number of intervals between range $[a,b]$, then
$$\Delta x = \frac{b-a}{3n}$$ And points you pick will be:
$$a=x_0<x_1<x_2<x_3<x_4\dots<x_{3n-2}<x_{3n-1}<x_{3n}=b$$ The points $x_{3i},x_{3i+1},x_{3i+2},x_{3i+3}$ for $i \in \mathbb{N}_0$ can be used to find single Third-order Approximation in interval of $[x_{3i},x_{3i+3}]$.
Using the formula above, we find that
$$\begin{align*}
\int_{x_{3i}}^{x_{3i+3}}f(x)dx &\approx \frac{x_{3i+3}-x_{3i}}{8} \left( f(x_{3i}) + 3f(x_{3i+1}) + 3f(x_{3i+2}) + f(x_{3i+3}) \right) \\
&\approx \frac{3 \Delta x}{8} \left( f(x_{3i}) + 3f(x_{3i+1}) + 3f(x_{3i+2}) + f(x_{3i+3})\right) \\
\therefore & \approx \frac{b-a}{8n}\left( f(x_{3i}) + 3f(x_{3i+1}) + 3f(x_{3i+2}) + f(x_{3i+3}) \right)
\end{align*}$$ Now adding up all the interval, we will have found the approximation that
$$\int_{a}^bf(x)dx \approx \frac{b-a}{8n}\left[ f(x_0) + 3f(x_1)+3f(x_2) + 2f(x_3) + 3f(x_4)+\dots + 3f(x_{3n -1}) + f(x_{3n}) \right] \\ \blacksquare$$

Notice the coefficient pattern of $$1\,3\,3\,2\,3\,3\,2\dots2\,3\,3\,1$$ A couple of 3s separated by a single 2, with 1 on start and end.

Visualizing Approximation

You can visualize the process as such:

Visualization of Simpson's Rule. In the three intervals, picks three points to draw Second Degree Polynomial (in red). The true area (int green) is approximated with area of polynomial in each intervals (in red).

Visualization of Third Degree Approximation. In the three intervals, picks four points to draw a Third Degree Polynomial (in red). The true area (in green) is approximated with area of polynomial in each intervals (in red).

Notice that for larger intervals, Third Degree Approximation is much closer to the actual curve itself.

Downside to Higher Degree Polynomial Approximations

So we have seen how for larger intervals, higher order polynomials produce closer approximation to actual curves, but these do come with their drawbacks.
Firstly, more points are needed to compute the approximations. Though this may not be much of a problem for computers, when approximating integral by hand, using Third Degree Approximation may not be practical compared to using Simpson's Rule or even just Riemann Sum. And for some curves, Second Degree Polynomial is enough to closely approximate the curve. For example, in the images above, approximation of Third Degree Polynomial was not that much better than that of Simpson's Rule in the third interval.
This leads to our second point: as we use more and more intervals and their range approach 0, the curve will become more and more linear within that interval, making usage of higher degree polynomial meaningless. Again, in the images above, though Third Degree Polynomial was much closely approximate the first interval than Simpson's Rule, if we simply had divided that interval into two, then Simpson's Rule would have produced a result not that different than that of Third Degree Approximation.
For these reasons, Third Degree Approximations will be most practical and useful when computing integral of large range with least number of intervals.

The Next Steps

The next steps from Third Degree Approximation would of course be Fourth Degree Approximation, but will these really be better than Third Degree Approximations or even Simpson's Rule? As we use finer and finer intervals, maybe Trapezoidal Sum was all we ever needed. 

But of course, for broader intervals and quicker and closer approximations, maybe higher degree approximations do have their benefits. And so, I will leave finding formula for Fourth Degree Approximations and maybe even finding patterns for n-th Degree Approximations as a fun exercise for you readers.