Calculus: Integrals II
Simpson's Rule
So far we've used rectangles and trapezoids to approximate integrals, but what about more advanced shapes? What if the top of the top of the slice of area we're computing could more closely follow the curve of the function? This is getting quite close to asking whether we can estimate the area under the curve by solving for the area under the curve, but of course knowing the exact answer precludes the need for an estimate in the first place. Taking a second look at the methods we've already covered will give us a hint as to what to do next.
The left Riemann sum, right Riemann sum, and midpoint method all estimate the value of the integral using the value of known integrals. Namely, they used the integrals of constant functions, which can be calculated directly as a rectangle. The trapezoid rule calculates the integral of a linear function, which can be calculated directly as a trapezoid. What if we could approximate the area under the curve using a quadratic function?
The answer is yes, we can do that! In fact, the method of approximating an integral by summing the integrals of quadratic functions is called Simpson's rule, and it will be explained in just a moment. But first, to answer a burning question - isn't this cheating? Where does the area underneath a quadratic curve come from if not from taking its integral? The short answer is that figuring out the exact solution to an integral is indeed required for deriving this method. However, this method is grouped in this integrals section because it follows the theme of approximating integrals with sums.
Approximating an integral using Simpson's rule is only a bit more complicated than using the trapezoid rule and is given by the following formula:
$$\int\limits_{a}^{b} f(x) \, dx = \dfrac{\Delta x}{3}\left(f(x_0)+ 4f(x_1) + 2(f_2) + 4f(x_3) + 2f(x_4) + \ldots + 4f(x_{n-1}) + f(x_{n}) \right)$$
The first and last terms have coefficients of 1, and the coefficients of the interior terms alternate between 4 and 2, always both beginning and ending with 4. The derivation is left as an exercise - with all the steps done out, of course!
An important note: Simpson's rule fits a 2nd degree polynomial through the sample points, so it can exactly calculate the area under any polynomial up to and including a quadratic. The rule operates just as the trapezoid rule can exactly solve for the area under any linear function, including constant functions.
Problems
Estimate $\displaystyle\int\limits_{0}^{2}4x+5 \, dx$ where $N = 4$. Do you expect Simpson's rule to provide an exact answer?
Since the integrand is a polynomial of degree 2 or lower, we should expect Simpson's rule to produce the exact result. Let's do out the sum and compare it to the actual result:
Step 1: Calculate $\Delta x$
$\Delta x = \dfrac{b-a}{N} \\ \Delta x = \dfrac{2-0}{4} \\ \Delta x = \dfrac{1}{2} \\ $
Step 2: Calculate sum:
$ \displaystyle\int\limits_{0}^{2}4x+5 \, dx \approx \dfrac{\Delta x}{3}\left( f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_{3}) + f(x_4) \right) \\ \displaystyle\int\limits_{0}^{2}4x+5 \, dx \approx \dfrac{\Delta x}{3}\left( (4(0)+5) + 4(4(0.5)+5) + 2(4(1)+5) + 4(4(1.5)+5) + (4(2)+5) \right) \\ \displaystyle\int\limits_{0}^{2}4x+5 \, dx \approx \dfrac{1}{6}\left( (4(0)+5) + 4(4(0.5)+5) + 2(4(1)+5) + 4(4(1.5)+5) + (4(2)+5) \right) \\ \displaystyle\int\limits_{0}^{2}4x+5 \, dx \approx \dfrac{1}{6}\left( 5 + 4(7) + 2(9) + 4(11) + 13 \right) \\ \displaystyle\int\limits_{0}^{2}4x+5 \, dx \approx \dfrac{1}{6}\left( 5 + 28 + 18 + 44 + 13 \right) \\ \displaystyle\int\limits_{0}^{2}4x+5 \, dx \approx \dfrac{1}{6}\left( 108 \right) \\ \displaystyle\int\limits_{0}^{2}4x+5 \, dx \approx 18 \\$
Let's now solve the integral analytically and compare the results:
$\displaystyle\int\limits_{0}^{2}4x+5 \, dx = \left. \dfrac{1}{2}(4)x^2 + 5x \right|_{0}^{2} \\ \displaystyle\int\limits_{0}^{2}4x+5 \, dx = \left. 2x^2 + 5x \right|_{0}^{2} \\ \displaystyle\int\limits_{0}^{2}4x+5 \, dx = 2(2)^2 + 5(2) - 2(0)^2 - 5(0) \\ \displaystyle\int\limits_{0}^{2}4x+5 \, dx = 8 + 10 - 0 - 0 \\ \displaystyle\int\limits_{0}^{2}4x+5 \, dx = 18 \\ $
As expected, Simpson's rule produces the exact result. How lovely.
Estimate $\displaystyle\int\limits_{-1}^{1}x^2 \, dx$ where N=4. Do you expect Simpson's rule to provide an exact answer?
Simpson's rule estimates the integral by fitting 2nd degree polynomials to the sample points and taking the integral of the resulting function. Since the integrand here is a quadratic, we should expect Simpson's rule to give us an exact answer. Let's take the integral and compare it against the exact solution to check our understanding:
Step 1: Calculate $\Delta x$:
$ \Delta x = \dfrac{b-a}{N} \\ \Delta x = \dfrac{1-(-1)}{4} \\ \Delta x = \dfrac{1}{2}\\ $
Step 2: Calculate sum:
$ \displaystyle\int\limits_{-1}^{1}x^2 \, dx \approx \dfrac{\Delta x}{3}\left(f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + f(x_4)\right) \\ \displaystyle\int\limits_{-1}^{1}x^2 \, dx \approx \dfrac{1}{6}\left((-1)^2 + 4(-0.5)^2 + 2(0)^2 + 4(0.5)^2 + 1^2\right) \\ \displaystyle\int\limits_{-1}^{1}x^2 \, dx \approx \dfrac{1}{6}\left(1 + 4(.25) + 0 + 4(0.25) + 1\right) \\ \displaystyle\int\limits_{-1}^{1}x^2 \, dx \approx \dfrac{1}{6}\left(1 + 1 + 0 + 1 + 1\right) \\ \displaystyle\int\limits_{-1}^{1}x^2 \, dx \approx \dfrac{1}{6}(4) \\ \displaystyle\int\limits_{-1}^{1}x^2 \, dx \approx \dfrac{2}{3} \\ $
Now let's calculate the area analytically and compare it to the answer above:
$ \displaystyle\int\limits_{-1}^{1} x^2 \, dx = \left. \dfrac{1}{3}x^3 \right|_{-1}^{1} \\ \displaystyle\int\limits_{-1}^{1} x^2 \, dx = \dfrac{1}{3}1^3 - \dfrac{1}{3}(-1)^3 \\ \displaystyle\int\limits_{-1}^{1} x^2 \, dx = \dfrac{1}{3} + \dfrac{1}{3}\\ \displaystyle\int\limits_{-1}^{1} x^2 \, dx = \dfrac{2}{3} \\ $
Look at that, a perfect match. It's almost like we know what we're doing.
Estimate $\displaystyle\int\limits_{0}^{\pi} \sin(x) \, dx$ where $N=4$. Will Simpson's rule give an exact answer?
Simpson's rule evaluates the area under second degree polynomials, but the sine function is not such a polynomial. Therefore, the result won't be exact. Let's solve for the estimate:
Step 1: Calculate $\Delta x$
$ \Delta x = \dfrac{b-a}{N} \\ \Delta x = \dfrac{\pi-0}{4} \\ \Delta x = \dfrac{\pi}{4} \\ $
Step 2: Calculate sum:
$ \displaystyle\int\limits_{0}^{\pi} \sin(x) \, dx \approx \dfrac{\Delta x}{3}\left( f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + f(x_4)\right) \\ \displaystyle\int\limits_{0}^{\pi} \sin(x) \, dx \approx \dfrac{\pi}{12}\left( f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + f(x_4)\right) \\ \displaystyle\int\limits_{0}^{\pi} \sin(x) \, dx \approx \dfrac{\pi}{12}\left( \sin(0) + 4\sin\left(\dfrac{\pi}{4}\right) + 2\sin\left(\dfrac{\pi}{2}\right) + 4\sin\left(\dfrac{3\pi}{4}\right) + \sin(\pi)\right) \\ \displaystyle\int\limits_{0}^{\pi} \sin(x) \, dx \approx \dfrac{\pi}{12}\left( 0 + 4\left(\dfrac{1}{\sqrt{2}}\right) + 2(1) + 4\left(\dfrac{1}{\sqrt{2}}\right) + 0\right) \\ \displaystyle\int\limits_{0}^{\pi} \sin(x) \, dx \approx \dfrac{\pi}{12}\left( 2 + \dfrac{8}{\sqrt{2}} \right) \\ \displaystyle\int\limits_{0}^{\pi} \sin(x) \, dx \approx \dfrac{\pi}{12}\left( 2 + 4\sqrt{2} \right) \\ \displaystyle\int\limits_{0}^{\pi} \sin(x) \, dx \approx \dfrac{\pi(1 + 2\sqrt{2})}{6} \\ $
Estimate $\displaystyle\int\limits_{0}^{2\pi} \cos(2\pi) \, d\theta$ where $N=8$
Step 1: Calculate $\Delta \theta$
$ \Delta \theta = \dfrac{b-a}{N} \\ \Delta \theta = \dfrac{2\pi - 0}{8} \\ \Delta \theta = \dfrac{\pi}{4} \\ $
Step 2: Calculate sum:
$ \displaystyle\int\limits_{0}^{2\pi} \cos(2\theta) \, d\theta \approx \dfrac{\Delta \theta}{3}\left( f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + 2f(x_4) + \ldots + 4f(x_7) + f(x_8) \right) \\ \displaystyle\int\limits_{0}^{2\pi} \cos(2\theta) \, d\theta \approx \dfrac{\pi}{12}\left( \cos(2(0)) + 4\cos\left(2\left(\dfrac{\pi}{4}\right)\right) + 2\cos\left(2\left(\dfrac{\pi}{2}\right)\right) + \\ 4\cos\left(2\left(\dfrac{3\pi}{4}\right)\right) + 2\cos(2(\pi)) + 4\cos\left(2\left(\dfrac{5\pi}{4}\right)\right) + 2\cos\left(2\left(\dfrac{3\pi}{2}\right)\right) + 4\cos\left(2\left(\dfrac{7\pi}{4}\right)\right) + \cos(2(2\pi))\right) \\ \phantom{0} \\ \displaystyle\int\limits_{0}^{2\pi} \cos(2\theta) \, d\theta \approx \dfrac{\pi}{12}\left( \cos(0) + 4\cos\left(\dfrac{\pi}{2}\right) + 2\cos\left(\pi\right) + \\ 4\cos\left(\dfrac{3\pi}{2}\right) + 2\cos(2\pi) + 4\cos\left(\dfrac{5\pi}{2}\right) + 2\cos(3\pi) + 4\cos\left(\dfrac{7\pi}{2}\right) + \cos(4\pi)\right) \\ \phantom{0} \\ \displaystyle\int\limits_{0}^{2\pi} \cos(2\theta) \, d\theta \approx \dfrac{\pi}{12} ( 1 + 4(0) + (2)(-1) + 4(0) + 2(1) + 4(0) - 2(1) + 4(0) + 1 ) \\ \displaystyle\int\limits_{0}^{2\pi} \cos(2\theta) \, d\theta \approx \dfrac{\pi}{12} ( 1 + 0 - 2 + 0 + 2 + 0 - 2 + 0 + 1 ) \\ \displaystyle\int\limits_{0}^{2\pi} \cos(2\theta) \, d\theta \approx \dfrac{\pi}{12} ( 0 ) \\ \displaystyle\int\limits_{0}^{2\pi} \cos(2\theta) \, d\theta \approx 0 \\ $
Estimate $\displaystyle\int\limits_{-2}^{2}x^3 + 2x \, dx$ where $N=6$
Step 1: Calculate $\Delta x$
$ \Delta x = \dfrac{2-(-2)}{6} \\ \Delta x = \dfrac{4}{6} \\ \Delta x = \dfrac{2}{3} \\ $
Step 2: Calculate sum:
$ \displaystyle\int\limits_{-2}^{2} x^3 + 2x \, dx \approx \dfrac{\Delta x}{3} \left( f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + 2f(x_4) + 4f(x_5) + f(x_6) \right) \\ \displaystyle\int\limits_{-2}^{2} x^3 + 2x \, dx \approx \dfrac{2}{9} \left( \left((-2)^3+2(-2)\right) + 4\left(\left(\dfrac{-4}{3}\right)^3+2\left(\dfrac{-4}{3}\right)\right) + 2\left(\left(\dfrac{-2}{3}\right)^3+2\left(\dfrac{-2}{3}\right)\right) + \\ 4((0)^3+2(0)) + 2\left(\left(\dfrac{2}{3}\right)^3+2\left(\dfrac{2}{3}\right)\right) + 4\left(\left(\dfrac{4}{3}\right)^3+2\left(\dfrac{4}{3}\right)\right) + \left(2^3+2(2)\right) \right) \\ \displaystyle\int\limits_{-2}^{2} x^3 + 2x \, dx \approx \dfrac{2}{9} \left( -\left((2)^3+2(2)\right) - 4\left(\left(\dfrac{4}{3}\right)^3+2\left(\dfrac{4}{3}\right)\right) - 2\left(\left(\dfrac{2}{3}\right)^3+2\left(\dfrac{2}{3}\right)\right) + \\ 4(0) + 2\left(\left(\dfrac{2}{3}\right)^3+2\left(\dfrac{2}{3}\right)\right) + 4\left(\left(\dfrac{4}{3}\right)^3+2\left(\dfrac{4}{3}\right)\right) + \left(2^3+2(2)\right) \right) \\ \displaystyle\int\limits_{-2}^{2} x^3 + 2x \, dx \approx \dfrac{2}{9} (0) \\ \displaystyle\int\limits_{-2}^{2} x^3 + 2x \, dx \approx 0 \\ $
Estimate $\displaystyle\int\limits_{1}^{2} \dfrac{2}{x} \, dx$ where $N=4$
A calculator may be of service for this problem.
Step 1: Calculate $\Delta x$
$ \Delta x = \dfrac{2-1}{4} \\ \Delta x = \dfrac{1}{4} \\ $
Step 2: Calculate sum:
$ \displaystyle\int\limits_{1}^{2} \dfrac{2}{x} \, dx \approx \dfrac{\Delta x}{3} \left( f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + f(x_4) \right) \\ \displaystyle\int\limits_{1}^{2} \dfrac{2}{x} \, dx \approx \dfrac{1}{12} \left( \dfrac{2}{1} + 4\left(\dfrac{2}{1+1/4}\right) + 2\left(\dfrac{2}{1+1/2}\right) + 4\left(\dfrac{2}{1+3/4}\right) + \dfrac{2}{2} \right) \\ \displaystyle\int\limits_{1}^{2} \dfrac{2}{x} \, dx \approx \dfrac{1}{12} \left( 2 + 4\left(\dfrac{8}{5}\right) + 2\left(\dfrac{4}{3}\right) + 4\left(\dfrac{8}{7}\right) + 1 \right) \\ \displaystyle\int\limits_{1}^{2} \dfrac{2}{x} \, dx \approx \dfrac{1}{12} \left( 2 + \dfrac{24}{5} + \dfrac{8}{3} + \dfrac{24}{7} + 1 \right) \\ \displaystyle\int\limits_{1}^{2} \dfrac{2}{x} \, dx \approx \dfrac{1}{12} \left( \dfrac{1459}{105} \right) \\ \displaystyle\int\limits_{1}^{2} \dfrac{2}{x} \, dx \approx \dfrac{1459}{1260} \\ $
Estimate $\displaystyle\int\limits_{-4}^{1} 3x^4 + 7x + 2x \, dx$ where $N = 6$
A calculator may be of service for this problem.
Step 1: Calculate $\Delta x$
$ \Delta x = \dfrac{-1-(-4)}{6} \\ \Delta x = \dfrac{3}{6} \\ \Delta x = \dfrac{1}{2} \\ $
Step 2: Calculate sum:
$ \displaystyle\int\limits_{-4}^{1} 3x^4 + 7x + 2x \, dx \approx \dfrac{\Delta x}{3} \left( f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + 2f(x_4) + 4f(x_5) + f(x_6) \right) \\ \displaystyle\int\limits_{-4}^{1} 3x^4 + 7x + 2x \, dx \approx \dfrac{1}{6} \left( \left( 3(-4)^4 + 7(-4) + 2 \right) + 4\left( 3\left(\dfrac{-7}{2}\right)^4 + 7\left(\dfrac{-7}{2}\right) + 2 \right) + \\ 2\left( 3(-3)^4 + 7(-3) + 2 \right) + 4\left( 3\left(\dfrac{-5}{2}\right)^4 + 7\left(\dfrac{-5}{2}\right) + 2 \right) + \\ 2\left( 3(-2)^4 + 7(-2) + 2 \right) + 4\left( 3\left(\dfrac{-3}{2}\right)^4 + 7\left(\dfrac{-3}{2}\right) + 2 \right) + \left(3(-1)^4 + 7(-1) + 2 \right) \right) \\ \displaystyle\int\limits_{-4}^{-1} 3x^4 + 7x + 2 \, dx \approx \dfrac{1}{6} \left( \left( 768 - 28 + 2 \right) + 4\left( \dfrac{7203}{16}-\dfrac{49}{2} + 2 \right) + 2\left( 243 - 21 + 2 \right) + \\ 4\left( \dfrac{1875}{16} - \dfrac{35}{2} + 2 \right) + 2\left( 48 - 14 + 2 \right) + 4\left( \dfrac{243}{16} - \dfrac{21}{2} + 2 \right) + \left(3 - 7 + 2 \right) \right) \\ \displaystyle\int\limits_{-4}^{-1} 3x^4 + 7x + 2 \, dx \approx \dfrac{1}{6} \left( (742) + 4\left( \dfrac{6843}{16} \right) + 2(224) + 4\left(\dfrac{1627}{16}\right) + 2( 36) + 4\left( \dfrac{107}{16} \right) + \left(-2\right) \right) \\ \displaystyle\int\limits_{-4}^{-1} 3x^4 + 7x + 2 \, dx \approx \dfrac{1}{6} \left( 742 + \dfrac{6843}{4} + 448 + \dfrac{1627}{4} + 72 + \dfrac{107}{4} - 2 \right) \\ \displaystyle\int\limits_{-4}^{-1} 3x^4 + 7x + 2 \, dx \approx \dfrac{1}{6} \left( \dfrac{13617}{4} \right) \\ \displaystyle\int\limits_{-4}^{-1} 3x^4 + 7x + 2 \, dx \approx \dfrac{4539}{8} \\ $
Estimate $\displaystyle\int\limits_{0}^{8} 2^x \, dx$ where $N=4$
Step 1: Calculate $\Delta x$
$ \Delta x = \dfrac{b-a}{N} \\ \Delta x = \dfrac{8-0}{4} \\ \Delta x = 2 \\ $
Step 2: Calculate sum:
$ \displaystyle\int\limits_{0}^{8} 2^x \, dx \approx \dfrac{\Delta x}{3} \left( f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + f(x_4) \right) \\ \displaystyle\int\limits_{0}^{8} 2^x \, dx \approx \dfrac{2}{3} \left( 2^0 + 4(2^2) + 2(2^4) + 4(2^6) + 2^8 \right) \\ \displaystyle\int\limits_{0}^{8} 2^x \, dx \approx \dfrac{2}{3} \left( 1 + 16 + 32 + 256 + 256 \right) \\ \displaystyle\int\limits_{0}^{9} 2^x \, dx \approx \dfrac{2}{3} (561) \\ \displaystyle\int\limits_{0}^{9} 2^x \, dx \approx 374 \\ $
Estimate $\displaystyle\int\limits_{0}^{8} 2^x \, dx$ where $N=8$.
Step 1: Calculate $\Delta x$
$ \Delta x = \dfrac{b-a}{N} \\ \Delta x = \dfrac{8-0}{8} \\ \Delta x = 1 \\ $
Step 2: Calculate sum:
$ \displaystyle\int\limits_{0}^{8} 2^x \, dx \approx \dfrac{\Delta x}{3} \left( f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + 2f(x_4) + \ldots + 4f(x_7) + f(x_8) \right) \\ \displaystyle\int\limits_{0}^{8} 2^x \, dx \approx \dfrac{1}{3} \left( 2^0 + 4(2^1) + 2(2^2) + 4(2^3) + 2(2^4) + 4(2^5) + 2(2^6) + 4(2^7) + 2^8 \right) \\ \displaystyle\int\limits_{0}^{8} 2^x \, dx \approx \dfrac{1}{3} \left( 1 + 8 + 8 + 32 + 32 + 128 + 128 + 512 + 256 \right) \\ \displaystyle\int\limits_{0}^{9} 2^x \, dx \approx \dfrac{1}{3} (1105) \\ \displaystyle\int\limits_{0}^{9} 2^x \, dx \approx \dfrac{1105}{3} \\ $
Show that evaluating $\displaystyle\int\limits_{a}^{b} \alpha x^2 + \beta x + \gamma \, dx$ produces the same result as evaluating the integral analytically. Note that this requires knowledge of how to evaluate integrals analytically, which is covered in the next section.
Hint: Solve analytically first.
Because the analytically solved integral will produce an easily simplifiable result, it will be easier to solve for it first. Solving using Simpson's rule will produce a big slew of summation terms, and it will be easier to rearrange these terms to match something simple than to try to do the reverse.
Evaluate analytically:
$ \displaystyle\int\limits_{a}^{b} \alpha x^2 + \beta x + \gamma \, dx = \left. \dfrac{1}{3}\alpha x^3 + \dfrac{1}{2}\beta x^2 + \gamma x \right|_{a}^{b} \\ \displaystyle\int\limits_{a}^{b} \alpha x^2 + \beta x + \gamma \, dx = \left( \dfrac{1}{3}\alpha b^3 + \dfrac{1}{2}\beta b^2 + \gamma b \right) - \left( \dfrac{1}{3}\alpha a^3 + \dfrac{1}{2}\beta a^2 + \gamma a \right) \\ \displaystyle\int\limits_{a}^{b} \alpha x^2 + \beta x + \gamma \, dx = \dfrac{1}{3}\alpha (b^3-a^3) + \dfrac{1}{2}\beta (b^2-a^2) + \gamma (b-a) \\ $
Evaluate using Simpson's rule:
$ \displaystyle\int\limits_{a}^{b} \alpha x^2 + \beta x + \gamma \, dx \approx \dfrac{\Delta x}{3}\left( f(x_0) + 4f(x_1) + f(x_2) \right) \\ \displaystyle\int\limits_{a}^{b} \alpha x^2 + \beta x + \gamma \, dx \approx \dfrac{b-a}{6}\left( \left(\alpha a^2 + \beta a + \gamma \right) + 4\left( \alpha\left(\dfrac{a+b}{2}\right)^2 + \beta\left(\dfrac{a+b}{2}\right) + \gamma \right) + \left(\alpha b^2 + \beta b + \gamma \right) \right) \\ \displaystyle\int\limits_{a}^{b} \alpha x^2 + \beta x + \gamma \, dx \approx \dfrac{b-a}{6}\left( \alpha\left(a^2 + 4\left(\dfrac{a+b}{2}\right)^2 + b^2\right) + \beta\left(a+4\left(\dfrac{a+b}{2}\right)+b\right) + 6\gamma \right) \\ \displaystyle\int\limits_{a}^{b} \alpha x^2 + \beta x + \gamma \, dx \approx \dfrac{b-a}{6}\left( \alpha\left(a^2 + 4\dfrac{(a+b)^2}{2^2} + b^2\right) + \beta(a+2(a+b)+b) + 6\gamma \right) \\ \displaystyle\int\limits_{a}^{b} \alpha x^2 + \beta x + \gamma \, dx \approx \dfrac{b-a}{6}\left( \alpha\left(a^2 + (a+b)^2 + b^2\right) + \beta(3a+3b) + 6\gamma \right) \\ \displaystyle\int\limits_{a}^{b} \alpha x^2 + \beta x + \gamma \, dx \approx \dfrac{b-a}{6}\left( \alpha\left(a^2 + a^2 + 2ab + b^2 + b^2\right) + \beta(3a+3b) + 6\gamma \right) \\ \displaystyle\int\limits_{a}^{b} \alpha x^2 + \beta x + \gamma \, dx \approx \dfrac{b-a}{6}\left( \alpha\left(2a^2 + 2ab + 2b^2\right) + \beta(3a+3b) + 6\gamma \right) \\ \displaystyle\int\limits_{a}^{b} \alpha x^2 + \beta x + \gamma \, dx \approx \dfrac{b-a}{6}\left( 2\alpha\left(a^2 + ab + b^2\right) + 3\beta(a+b) + 6\gamma \right) \\ \displaystyle\int\limits_{a}^{b} \alpha x^2 + \beta x + \gamma \, dx \approx \dfrac{1}{6}\left( 2\alpha\left(\left(a^2 + ab + b^2\right)(b-a)\right) + 3\beta(a+b)(b-a) + 6\gamma(b-a) \right) \\ \displaystyle\int\limits_{a}^{b} \alpha x^2 + \beta x + \gamma \, dx \approx \dfrac{1}{6}\left( 2\alpha\left(a^2b + ab^2 + b^3 - a^3 - a^2b - ab^2\right) + 3\beta(b^2-a^2) + 6\gamma(b-a) \right) \\ \displaystyle\int\limits_{a}^{b} \alpha x^2 + \beta x + \gamma \, dx \approx \dfrac{1}{6}\left( 2\alpha\left(b^3 - a^3\right) + 3\beta\left(b^2-a^2\right) + 6\gamma(b-a) \right) \\ \displaystyle\int\limits_{a}^{b} \alpha x^2 + \beta x + \gamma \, dx \approx \dfrac{1}{3}\alpha\left(b^3 - a^3\right) + \dfrac{1}{2}\beta\left(b^2-a^2\right) + \gamma(b-a) \\ $
Just as anticipated! Isn't math just the darned coolest thing?
Derive Simpson's rule. First, fit a parabola between 3 sample points, $x_{-1}, x_{i},$ and $x_{i+1}$. Then, integrate the result. Then sum up the terms to get the final rule.
Hint: Use Lagrange interpolation.
Step 1: Fit a parabola through the sample points:
Using Lagrange's method for fitting an $n$th degree polynomial, call it $g(x)$, through $n+1$ points, such as those on a function $f(x)$, gives us the following:
$ g_i(x) = f(x_{i-1})\dfrac{(x-x_{i})(x-x_{i+1})}{(x_{i-1}-x_{i})(x_{i-1}-x_{i+1})} + f(x_{i})\dfrac{(x-x_{i-1})(x-x_{i+1})}{(x_{i}-x_{i-1})(x_{i}-x_{i+1})} + f(x_{i+1})\dfrac{(x-x_{i-1})(x-x_{i})}{(x_{i+1}-x_{i-1})(x_{i+1}-x_{i})} \\ $
Step 2: Integrate:
$ \displaystyle\int\limits_{x_{i-1}}^{x_{x+1}} g_i(x) \, dx = \displaystyle\int\limits_{x_{i-1}}^{x_{x+1}} f(x_{i-1})\dfrac{(x-x_{i})(x-x_{i+1})}{(x_{i-1}-x_{i})(x_{i-1}-x_{i+1})} + f(x_{i})\dfrac{(x-x_{i-1})(x-x_{i+1})}{(x_{i}-x_{i-1})(x_{i}-x_{i+1})} + \\ f(x_{i+1})\dfrac{(x-x_{i-1})(x-x_{i})}{(x_{i+1}-x_{i-1})(x_{i+1}-x_{i})} \, dx \\ $
Before we get going, it will be helpful to split this integral up into the sum of three smaller integrals. We can then evaluate each integral independently to keep the page neater.
$ \displaystyle\int\limits_{x_{i-1}}^{x_{x+1}} g_i(x) \, dx = \displaystyle\int\limits_{x_{i-1}}^{x_{x+1}} f(x_{i-1})\dfrac{(x-x_{i})(x-x_{i+1})}{(x_{i-1}-x_{i})(x_{i-1}-x_{i+1})} \, dx + \displaystyle\int\limits_{x_{i-1}}^{x_{x+1}} f(x_{i})\dfrac{(x-x_{i-1})(x-x_{i+1})}{(x_{i}-x_{i-1})(x_{i}-x_{i+1})} \, dx + \\ \displaystyle\int\limits_{x_{i-1}}^{x_{x+1}} f(x_{i+1})\dfrac{(x-x_{i-1})(x-x_{i})}{(x_{i+1}-x_{i-1})(x_{i+1}-x_{i})} \, dx \\ \phantom{0} \\ \displaystyle\int\limits_{x_{i-1}}^{x_{x+1}} g_i(x) \, dx = A + B + C \\ $
Step 2 A:
$ A = \displaystyle\int\limits_{x_{i-1}}^{x_{x+1}} f(x_{i-1})\dfrac{(x-x_{i})(x-x_{i+1})}{(x_{i-1}-x_{i})(x_{i-1}-x_{i+1})} \, dx \\ A = \displaystyle\int\limits_{x_{i-1}}^{x_{x+1}} f(x_{i-1})\dfrac{ x^2 - (x_{i}+x_{i+1})x + x_{i}x_{i+1} } {(x_{i-1}-x_{i})(x_{i-1}-x_{i+1})}\, dx \\ A = \left. f(x_{i-1}) \dfrac{ \frac{1}{3}x^3 - \frac{1}{2}(x_{i}+x_{i+1})x^2 + x_{i}x_{i+1}x } {(x_{i-1}-x_{i})(x_{i-1}-x_{i+1})} \right|_{x_{i-1}}^{x_{i+1}} \\ A = f(x_{i-1}) \dfrac{ \frac{1}{3}\left(x_{i+1}^3-x_{i-1}^3\right) - \frac{1}{2}(x_{i}+x_{i+1})\left(x_{i+1}^2-x_{i-1}^2\right) + x_{i}x_{i+1}(x_{i+1}-x_{i-1}) } {(x_{i-1}-x_{i})(x_{i-1}-x_{i+1})} \\ A = \dfrac{1}{6}f(x_{i-1}) \dfrac{ 2\left(x_{i+1}^3-x_{i-1}^3\right) - 3(x_{i}+x_{i+1})\left(x_{i+1}^2-x_{i-1}^2\right) + 6x_{i}x_{i+1}(x_{i+1}-x_{i-1}) } {(x_{i-1}-x_{i})(x_{i-1}-x_{i+1})} \\ A = \dfrac{1}{6}f(x_{i-1}) \dfrac{ 2\left(x_{i+1}^2+x_{i+1}x_{i-1}+x_{i-1}^2\right)\left(x_{i+1}-x_{i-1}\right) - 3\left(x_{i}+x_{i+1}\right)\left(x_{i+1}+x_{i-1}\right)\left(x_{i+1}-x_{i-1}\right) + 6x_{i}x_{i+1}(x_{i+1}-x_{i-1}) } {(x_{i-1}-x_{i})(x_{i-1}-x_{i+1})} \\ A = \dfrac{1}{6}f(x_{i-1}) \dfrac{ -2\left(x_{i+1}^2+x_{i+1}x_{i-1}+x_{i-1}^2\right) + 3\left(x_{i}+x_{i+1}\right)\left(x_{i+1}+x_{i-1}\right) - 6x_{i}x_{i+1} } {x_{i-1}-x_{i}} \\ A = \dfrac{1}{6}f(x_{i-1}) \dfrac{ -2\left(x_{i+1}^2+x_{i+1}x_{i-1}+x_{i-1}^2\right) + 3\left(x_{i}x_{i+1} + x_{i}x_{i-1} + x_{i+1}^2 + x_{i+1}x_{i-1}\right) - 6x_{i}x_{i+1} } {x_{i-1}-x_{i}} \\ A = \dfrac{1}{6}f(x_{i-1}) \dfrac{ \left(-2x_{i+1}^2-2x_{i+1}x_{i-1}-2x_{i-1}^2\right) + \left(3x_{i}x_{i+1} + 3x_{i}x_{i-1} + 3x_{i+1}^2 + 3x_{i+1}x_{i-1}\right) - 6x_{i}x_{i+1} } {x_{i-1}-x_{i}} \\ A = \dfrac{1}{6}f(x_{i-1}) \dfrac{ \left(x_{i+1}^2 + x_{i+1}x_{i-1}-2x_{i-1}^2\right) + \left(-3x_{i}x_{i+1} + 3x_{i}x_{i-1} \right) } {x_{i-1}-x_{i}} \\ $
It seems we might be stuck at this point. However, remember that because the intervals are evenly spaced, $x_i = \dfrac{x_{i-1}+x_{i+1}}{2}$. Let's make a substitution:
$ A = \dfrac{1}{6}f(x_{i-1}) \dfrac{ \left(x_{i+1}^2 + x_{i+1}x_{i-1}-2x_{i-1}^2\right) + \left(-3\left(\frac{x_{i-1}+x_{i+1}}{2}\right)x_{i+1} + 3\left(\frac{x_{i-1}+x_{i+1}}{2}\right)x_{i-1} \right) } {x_{i-1}-\frac{x_{i-1}+x_{i+1}}{2}} \\ A = \dfrac{1}{6}f(x_{i-1}) \dfrac{ \left(2x_{i+1}^2 + 2x_{i+1}x_{i-1}-4x_{i-1}^2\right) + \left(-3\left(x_{i-1}+x_{i+1}\right)x_{i+1} + 3\left(x_{i-1}+x_{i+1}\right)x_{i-1} \right) } {2x_{i-1}-(x_{i-1}+x_{i+1})} \\ A = \dfrac{1}{6}f(x_{i-1}) \dfrac{ \left(2x_{i+1}^2 + 2x_{i+1}x_{i-1}-4x_{i-1}^2\right) + \left(-3x_{i-1}x_{i+1}-3x_{i+1}^2 + 3x_{i-1}^2 + 3x_{i-1}x_{i+1} \right) } {x_{i-1}-x_{i+1}} \\ A = \dfrac{1}{6}f(x_{i-1}) \dfrac{ \left(-x_{i+1}^2 + 2x_{i+1}x_{i-1}-x_{i-1}^2\right) } {x_{i-1}-x_{i+1}} \\ A = \dfrac{1}{6}f(x_{i-1}) \dfrac{ -(x_{i+1}-x_{i-1})(x_{i+1}-x_{i-1}) } {x_{i-1}-x_{i+1}} \\ A = \dfrac{1}{6}f(x_{i-1}) \dfrac{ (x_{i-1}-x_{i+1})(x_{i+1}-x_{i-1}) } {x_{i-1}-x_{i+1}} \\ A = \dfrac{1}{6}f(x_{i-1}) (x_{i+1}-x_{i-1}) \\ $
Remember that $\Delta x = \dfrac{b-a}{N}$. Since $N=2$, we can make the following simplification:
$ A = \dfrac{\Delta x}{3}f(x_{i-1}) \\ $
Wow, that was painful. But if Simpson could do it, so can we.
Step 2 B:
$ B = \displaystyle\int\limits_{x_{i-1}}^{x_{x+1}} f(x_{i})\dfrac{(x-x_{i-1})(x-x_{i+1})}{(x_{i}-x_{i-1})(x_{i}-x_{i+1})} \, dx \\ B = \displaystyle\int\limits_{x_{i-1}}^{x_{x+1}} f(x_{i})\dfrac{x^2 -(x_{i-1}+x_{i+1})x + x_{i-1}x_{i+1}}{(x_{i}-x_{i-1})(x_{i}-x_{i+1})} \, dx \\ B = \left. f(x_{i})\dfrac{\frac{1}{3}x^3 -\frac{1}{2}(x_{i-1}+x_{i+1})x^2 + (x_{i-1}x_{i+1})x}{(x_{i}-x_{i-1})(x_{i}-x_{i+1})} \right|_{x_{i-1}}^{x_{i+1}} \\ B = f(x_{i})\dfrac{\frac{1}{3}(x_{i+1}^3-x_{i-1}^3) -\frac{1}{2}(x_{i-1}+x_{i+1})(x_{i+1}^2-x_{i-1}^2) + (x_{i-1}x_{i+1})(x_{i+1}-x_{i-1})}{(x_{i}-x_{i-1})(x_{i}-x_{i+1})} \\ B = f(x_{i})\dfrac{\frac{1}{3}(x_{i+1}^3-x_{i-1}^3) -\frac{1}{2}(x_{i-1}+x_{i+1})(x_{i+1}^2-x_{i-1}^2) + (x_{i-1}x_{i+1})(x_{i+1}-x_{i-1})}{(x_{i}-x_{i-1})(x_{i}-x_{i+1})} \\ $
Once again, we need to use $x_i = \frac{x_{i-1}+x_{i+1}}{2}$ to make a substitution:
$ B = f(x_{i})\dfrac{\frac{1}{3}(x_{i+1}^3-x_{i-1}^3) -\frac{1}{2}(x_{i-1}+x_{i+1})(x_{i+1}^2-x_{i-1}^2) + (x_{i-1}x_{i+1})(x_{i+1}-x_{i-1})}{\left(\frac{x_{i-1}+x_{i+1}}{2}-x_{i-1}\right)\left(\frac{x_{i-1}+x_{i+1}}{2}-x_{i+1}\right)} \\ B = 4f(x_{i})\dfrac{\frac{1}{3}(x_{i+1}^3-x_{i-1}^3) -\frac{1}{2}(x_{i-1}+x_{i+1})(x_{i+1}^2-x_{i-1}^2) + (x_{i-1}x_{i+1})(x_{i+1}-x_{i-1})}{(x_{i-1}+x_{i+1}-2x_{i-1})(x_{i-1}+x_{i+1}-2x_{i+1})} \\ B = 4f(x_{i})\dfrac{\frac{1}{3}(x_{i+1}^3-x_{i-1}^3) -\frac{1}{2}(x_{i-1}+x_{i+1})(x_{i+1}^2-x_{i-1}^2) + (x_{i-1}x_{i+1})(x_{i+1}-x_{i-1})}{(x_{i+1}-x_{i-1})(x_{i-1}-x_{i+1})} \\ B = 4f(x_{i})\dfrac{\frac{1}{3}(x_{i+1}^2+x_{i-1}x_{i+1}+x_{i-1}^2)(x_{i+1}-x_{i-1}) -\frac{1}{2}(x_{i-1}+x_{i+1})(x_{i+1}+x_{i-1})(x_{i+1}-x_{i-1}) + (x_{i-1}x_{i+1})(x_{i+1}-x_{i-1})}{(x_{i+1}-x_{i-1})(x_{i-1}-x_{i+1})} \\ B = 4f(x_{i})\dfrac{\frac{1}{3}(x_{i+1}^2+x_{i-1}x_{i+1}+x_{i-1}^2) -\frac{1}{2}(x_{i-1}+x_{i+1})(x_{i+1}+x_{i-1}) + (x_{i-1}x_{i+1})}{x_{i-1}-x_{i+1}} \\ B = \dfrac{2}{3}f(x_{i})\dfrac{2(x_{i+1}^2+x_{i-1}x_{i+1}+x_{i-1}^2) -3(x_{i-1}+x_{i+1})(x_{i+1}+x_{i-1}) + 6(x_{i-1}x_{i+1})}{x_{i-1}-x_{i+1}} \\ B = \dfrac{2}{3}f(x_{i})\dfrac{2x_{i+1}^2+2x_{i-1}x_{i+1}+x_{i-1}^2 - 3x_{i-1}^2 -6x_{i-1}x_{+1} - 3x_{i+1}^2 + 6x_{i-1}x_{i+1}}{x_{i-1}-x_{i+1}} \\ B = \dfrac{2}{3}f(x_{i})\dfrac{ 2x_{i-1}x_{i+1} - x_{i-1}^2 - x_{i+1}^2 }{x_{i-1}-x_{i+1}} \\ B = \dfrac{2}{3}f(x_{i})\dfrac{ -(x_{i-1}^2 - 2x_{i-1}x_{i+1} + x_{i+1}^2) }{x_{i-1}-x_{i+1}} \\ B = \dfrac{2}{3}f(x_{i})\dfrac{ -(x_{i-1}-x_{i+1})^2 }{x_{i-1}-x_{i+1}} \\ B = \dfrac{2}{3}f(x_{i})\dfrac{ (x_{i-1}-x_{i+1})^2 }{x_{i+1}-x_{i-1}} \\ B = \dfrac{2}{3}f(x_{i})\left(x_{i-1}-x_{i+1} \right) \\ $
Remembering again that $\Delta x = \dfrac{b-a}{N}$. Since $N=2$, we can make the following simplification:
$ B = \dfrac{4\Delta x}{3}f(x_{i-1}) \\ $
Okay, just one more nasty integral and we'll be out of the Fire Swamp.
Step 2 C:
$ C = \displaystyle\int\limits_{x_{i-1}}^{x_{x+1}} f(x_{i+1})\dfrac{(x-x_{i-1})(x-x_{i})}{(x_{i+1}-x_{i-1})(x_{i+1}-x_{i})} \, dx \\ C = \displaystyle\int\limits_{x_{i-1}}^{x_{x+1}} f(x_{i+1})\dfrac{x^2-(x_{i}+x_{i-1})x + x_{i-1}x_{i}}{(x_{i+1}-x_{i-1})(x_{i+1}-x_{i})} \, dx \\ C = \left. f(x_{i+1})\dfrac{\frac{1}{3}x^3-\frac{1}{2}(x_{i}+x_{i-1})x^2 + x_{i-1}x_{i}x}{(x_{i+1}-x_{i-1})(x_{i+1}-x_{i})} \right|_{x_{i-1}}^{x_{i}} \\ C = f(x_{i+1})\dfrac{\frac{1}{3}(x_{i+1}^3-x_{i-1}^3)-\frac{1}{2}(x_{i}+x_{i-1})(x_{i+1}^2-x_{i-1}^2) + x_{i-1}x_{i}(x_{i+1}-x_{i-1})}{(x_{i+1}-x_{i-1})(x_{i+1}-x_{i})} \\ C = \dfrac{1}{6}f(x_{i+1})\dfrac{2(x_{i+1}^3-x_{i-1}^3)-3(x_{i}+x_{i-1})(x_{i+1}^2-x_{i-1}^2) + 6x_{i-1}x_{i}(x_{i+1}-x_{i-1})}{(x_{i+1}-x_{i-1})(x_{i+1}-x_{i})} \\ C = \dfrac{1}{6}f(x_{i+1})\dfrac{2(x_{i-1}^2+x_{i-1}x_{i+1}+x_{i+1}^2)(x_{i+1}-x_{i-1})-3(x_{i}+x_{i-1})(x_{i+1}+x_{i-1})(x_{i+1}-x_{i-1}) + 6x_{i-1}x_{i}(x_{i+1}-x_{i-1})}{(x_{i+1}-x_{i-1})(x_{i+1}-x_{i})} \\ C = \dfrac{1}{6}f(x_{i+1})\dfrac{2(x_{i-1}^2+x_{i-1}x_{i+1}+x_{i+1}^2)-3(x_{i}+x_{i-1})(x_{i+1}+x_{i-1}) + 6x_{i-1}x_{i}}{x_{i+1}-x_{i}} \\ C = \dfrac{1}{6}f(x_{i+1})\dfrac{2\left(x_{i-1}^2+x_{i-1}x_{i+1}+x_{i+1}^2\right) - 3\left(x_{i}x_{i+1} + x_{i-1}x_{i} + x_{i-1}x_{i+1} + x_{i-1}^2\right) + 6x_{i-1}x_{i}}{x_{i+1}-x_{i}} \\ C = \dfrac{1}{6}f(x_{i+1})\dfrac{\left(2x_{i-1}^2+2x_{i-1}x_{i+1}+2x_{i+1}^2\right) + \left(-3x_{i}x_{i+1} - 3x_{i-1}x_{i} - 3x_{i-1}x_{i+1} - 3x_{i-1}^2\right) + 6x_{i-1}x_{i}}{x_{i+1}-x_{i}} \\ C = \dfrac{1}{6}f(x_{i+1})\dfrac{\left(-x_{i-1}^2-x_{i-1}x_{i+1}+2x_{i+1}^2\right) + \left(-3x_{i}x_{i+1} + 3x_{i-1}x_{i} \right) }{x_{i+1}-x_{i}} \\ $
For the third and final time, we need to use $x_i = \frac{x_{i-1}+x_{i+1}}{2}$ to make a substitution:
$ C = \dfrac{1}{6}f(x_{i+1})\dfrac{\left(-x_{i-1}^2-x_{i-1}x_{i+1}+2x_{i+1}^2\right) + \left(-3\left(\frac{x_{i-1}+x_{i+1}}{2}\right)x_{i+1} + 3x_{i-1}\left(\frac{x_{i-1}+x_{i+1}}{2}\right)\right) }{x_{i+1}-\left(\frac{x_{i-1}+x_{i+1}}{2}\right)} \\ C = \dfrac{1}{6}f(x_{i+1})\dfrac{\left(-2x_{i-1}^2-2x_{i-1}x_{i+1}+4x_{i+1}^2\right) + \left(-3\left(x_{i-1}+x_{i+1}\right)x_{i+1} + 3x_{i-1}\left(x_{i-1}+x_{i+1}\right)\right) }{2x_{i+1}-x_{i-1}-x_{i+1}} \\ C = \dfrac{1}{6}f(x_{i+1})\dfrac{\left(-2x_{i-1}^2-2x_{i-1}x_{i+1}+4x_{i+1}^2\right) + \left(-3x_{i-1}x_{i+1}-3x_{i+1}^2 + 3x_{i-1}^2 + 3x_{i-1}x_{i+1} \right) }{x_{i+1}-x_{i-1}} \\ C = \dfrac{1}{6}f(x_{i+1})\dfrac{x_{i+1}^2 -2 x_{i-1}x_{i+1} + x_{i-1}^2 }{x_{i+1}-x_{i-1}} \\ C = \dfrac{1}{6}f(x_{i+1})\dfrac{(x_{i+1} - x_{i-1})^2}{x_{i+1}-x_{i-1}} \\ C = \dfrac{1}{6}f(x_{i+1})(x_{i+1} - x_{i-1}) \\ $
Remembering once last time that $\Delta x = \dfrac{b-a}{N}$. Since $N=2$, we can make the following simplification:
$ C = \dfrac{\Delta x}{3}f(x_{i+1}) \\ $
Step 2 D: Reassemble the integral:
$ \displaystyle\int\limits_{x_{i-1}}^{x_{i+1}} g_i(x) \, dx = A + B + C \\ \displaystyle\int\limits_{x_{i-1}}^{x_{i+1}} g_i(x) \, dx = \dfrac{\Delta x}{3}f(x_{i-1}) + \dfrac{4\Delta x}{3}f(x_{i}) + \dfrac{\Delta x}{3}f(x_{i+1}) \\ \displaystyle\int\limits_{x_{i-1}}^{x_{i+1}} g_i(x) \, dx = \dfrac{\Delta x}{3}\left(f(x_{i-1}) + 4f(x_{i}) + f(x_{i+1}) \right) \\ $
Step 3: Assemble the sum:
Since the integral spans 2 sampling widths, the next integral should start at $i+2$ so as not to estimate a portion of the same region twice. Writing this out gives the following:
$ \displaystyle\int\limits_{x_{i-1}}^{x_{i+3}} g_i(x) \, dx + \displaystyle\int\limits_{x_{i-1}}^{x_{i+3}} g_{i+2}(x) \, dx = \dfrac{\Delta x}{3}\left(f(x_{i-1}) + 4f(x_{i}) + f(x_{i+1}) \right) + \dfrac{\Delta x}{3}\left(f(x_{i+2-1}) + 4f(x_{i+2}) + f(x_{i+2+1}) \right) \\ \displaystyle\int\limits_{x_{i-1}}^{x_{i+3}} g_i(x) \, dx + \displaystyle\int\limits_{x_{i-1}}^{x_{i+3}} g_{i+2}(x) \, dx = \dfrac{\Delta x}{3}\left(f(x_{i-1}) + 4f(x_{i}) + f(x_{i+1}) \right) + \dfrac{\Delta x}{3}\left(f(x_{i+1}) + 4f(x_{i+2}) + f(x_{i+3}) \right) \\ \displaystyle\int\limits_{x_{i-1}}^{x_{i+3}} g_i(x) \, dx + \displaystyle\int\limits_{x_{i-1}}^{x_{i+3}} g_{i+2}(x) \, dx = \dfrac{\Delta x}{3}\left(f(x_{i-1}) + 4f(x_{i}) + 2f(x_{i+1}) + 4f(x_{i+2}) + f(x_{i+3}) \right) \\ $
Summing the areas under successive $2$nd degree polynomials in this fashion gives the desired result.