# Numerical Integration with Simpson's Rule¶

## Formula¶

Unlike the Riemann sum methods, Stewart does not give a sigma formula for Simpson's Rule. Instead, the textbook defines it as a series:

$$\int_{a}^b f(x) dx \approx S_n = \frac{\Delta x}{3} [f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + ... + 2f(x_{n-2}) + 4f(x_{n-1}) + f(x_n)]$$

Where $n$ is even and $\Delta x = \frac{b - a}{n}$.

The authors say "note the pattern of the coefficients" (pg. 520), which feels like a cop-out. What does $x_0$ mean? How do I find the actual values I'm supposed to plug into the function being integrated? Attempting to convert the formula into sigma form made these blindspots obvious.

$$\frac{\Delta x}{3} \sum_{i = 0}^{n} \Delta x \ S(i)$$

Where $n$ is even and $\Delta x = \frac{b - a}{n}$.

While it is possible to start the summation at 1 (to match the formulas for the left and right Riemann sums), I think it is simpler to start at 0.

What is $S(i)$? This is a special piecewise function that determines the coefficient (if any) to multiply $f(x_{i})$ by.

$$S(i) = \begin{cases} f(x_{i}), & \text{if i is 0 or n} \\ 2 \ f(x_{i}), & \text{if i is even} \\ 4 \ f(x_{i}), & \text{if i is odd} \end{cases}$$

Note that precedence is important here. I don't know how to represent this in math notation, but programatically the logic is:

In [ ]:
# Will not actually run because n, x_i, and f(x) are not defined
def S(i):
if (i == 0) or (i == n):
return f(x_i)
elif (i % 2 == 0):
return 2 * f(x_i)
else:
# odd
return 4 * f(x_i)


Finally, $x_{i}$ is a function that returns the value to plug into $f(x)$ given the current index $i$.

$$x_{i} = a + i \ \Delta x$$

The notation, $x_{i}$, is somewhat confusing because it looks like a variable but is actually a function.

## Worked Example¶

This is from an old uOttawa Calculus I exam.

Approximate the integral $\int_{0}^4 \frac{1}{x+1} dx$ numerically to 4 decimal places with n = 4 subdivisions with Simpson's Rule.

In :
# Lower bound
a = 0
# Upper bound
b = 4
# Number of intervals
n = 4
# Difference in x between intervals
deltaX = (b - a) / n

In :
# Function being integrated
def f(x):
return 1.0 / (x + 1.0)

In :
# Determines what value to plug into f(x) for each i
def currentX(i):
return a + (i * deltaX)

In :
# Function that applies Simpson's coefficients (based on i)
def S(i):
# if zeroth and last values
if (i == 0) or (i == n):
# no coefficient
return f(currentX(i))

# if even
elif (i % 2 == 0):
# coefficient of 2
return 2 * f(currentX(i))

# if odd
else:
# coefficient of 4
return 4 * f(currentX(i))


Note that the return statements are different from the previously defined $S(i)$ function, here $currentX()$ is used to determine the value to submit to f(x).

In :
# Computation
simpsonSum = 0.0

# Loop from 0 to n (inclusive)
for i in range(0, n+1):
simpsonSum += S(i)

# Don't forget to complete the formula by multiplying the sum by (deltaX / 3.0)!
(deltaX / 3.0)*simpsonSum

Out:
1.622222222222222

Which is roughly equivalent to Wolfram Alpha's answer.