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 [35]:
# Lower bound
a = 0
# Upper bound
b = 4
# Number of intervals
n = 4
# Difference in x between intervals
deltaX = (b - a) / n
In [36]:
# Function being integrated
def f(x):
    return 1.0 / (x + 1.0)
In [37]:
# Determines what value to plug into f(x) for each i
def currentX(i):
    return a + (i * deltaX)
In [38]:
# 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 [39]:
# 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[39]:
1.622222222222222

Which is roughly equivalent to Wolfram Alpha's answer.