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.

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]:

Which is roughly equivalent to Wolfram Alpha's answer.