Simpson’s Rule for Approximating Definite Integrals in R

Simpson’s rule is another closed Newton-Cotes formula for approximating integrals over an interval with equally spaced nodes. Unlike the trapezoidal rule, which employs straight lines to approximate a definite integral, Simpson’s rule uses the third Lagrange polynomial, \(P_3(x)\) to approximate the definite integral and as such can give exact results when approximating integrals of up to third-degree polynomials. Simpson’s rule is also generally more performant and often converges faster compared to the trapezoid rule.

Simpson’s rule is defined as:

\( \int^{x_2}_{x_0} f(x) \space dx = \frac{h}{3}[f(x_0) + 4f(x_1) + f(x_2)] – \frac{h^2}{90}f^{(4)}(\epsilon) \)

Where \(x_1 = x_0 + h\), \(h = (x_2 – x_0)/2\) and \(\frac{h^2}{90}f^{(4)}(\epsilon)\) is the error term. The fourth order error term indicating Simpson’s rule will return an exact result when the polynomial in question has a degree of three or less.

Simpson’s Rule Example

Let’s see if Simpson’s rule is more accurate in approximating the following function compared to the result obtained from the Trapezoidal rule.

\( \int_0^{\frac{\pi}{4}} x \space sin \space x \space dx \)

Approximating the integral with Simpson’s rule:

\( h = \frac{\pi/4 – 0}{2} = \frac{\pi}{8}, \qquad x_1 = 0 + h = \frac{\pi}{8} \) \( \int_0^{\frac{\pi}{4}} x \space sin \space x \space dx \approx \frac{\pi/8}{3} [0 + 0.6011177 + 0.5553604] = .1513826 \)

Integrating the function over the interval \([0, \frac{\pi}{4}]\) yields \(\frac{\pi – 4}{4\sqrt{2}} = 0.15175\), which gives an approximation error of \(-0.0003674\). Compare this to the approximation found by the trapezoidal rule, \(0.2180895\), which results in an error value of \(0.0663395\).

Visualizing how the two methods approximate the integral over the interval gives a good example of why Simpson’s rule is typically more accurate compared to the trapezoidal rule.

# Replicate the polynomial in question in R as a function.
f <- function(x) {
  return(x * sin(x))
}

# Construct the points using the intervals for plotting
df <- data.frame(cbind(c(0, pi/4, pi/4, 0), c(0, f(pi/4), 0, 0)))
curvedf <- data.frame(cbind(c(0, pi/8, pi/4), c(0, f(pi/8), f(pi/4))))


# As Simpson's rule uses the third Lagrange polynomial for approximation, we first 
# interpolate the polynomial using Lagrange's method. The `lagrange.poly()` function 
# is from a previous post on interpolating polynomials with the Lagrange approach, found
# here: http://wp.me/p4aZEo-5RJ

lagrange.poly(curvedf$X1, curvedf$X2)
## [1] "0.058260083543634*x + 0.826137273909778*x**2"


f2 <- function(x) {
  return(0.058260083543634*x + 0.826137273909778*x^2)
}

# Plot the function and the trapezoidal and Simpson's rule approximations.
ggplot(data = data.frame(x = 0), mapping = aes(x = x)) + 
  stat_function(fun = f, color = 'blue') + xlim(c(0,pi/4)) + 
  stat_function(fun = f2) + xlim(c(0,pi/4)) + 
  geom_segment(aes(x = 0, y = 0, xend = pi/4, yend = f(pi/4))) + 
  geom_segment(aes(x = pi/4, y = 0, xend = pi/4, yend = f(pi/4))) + 
  geom_polygon(data = df, aes(x=X1, y=X2), fill = 'blue', alpha = 0.2) + 
  geom_area(stat = 'function', fun = f2, fill = 'black', alpha = 0.4, xlim = c(0, pi/4))

The blue shaded area represents the approximation of the integral over the interval by the trapezoidal rule while the blue-grayish shaded area is the approximation of Simpson’s rule. The blue curve is the polynomial \(x \space sin(x)\). We can see the approximation given by Simpson’s rule is much more close to the actual definite integral compared to the approximation given by the trapezoidal rule.

We can also write a function to approximate a definite integral using Simpson’s rule:

simpsons.rule <- function(f, a, b) {
  if (is.function(f) == FALSE) {
    stop('f must be a function with one parameter (variable)')
  }
  
  h <- (b - a) / 2
  x0 <- a
  x1 <- a + h
  x2 <- b
  
  s <- (h / 3) * (f(x0) + 4 * f(x1) + f(x2))
  
  return(s)
}

simpsons.rule(f, 0, pi/4)
## [1] 0.1513826

Simpson’s Rule Second Example

Approximate the following definite integral using Simpson’s rule and compare again to the result of the trapezoid rule, which we found in a previous post.

\( \int^1_0 x^2 \space e^{-x} \space dx \)
\( h = \frac{(x_2 – x_0)}{2} = \frac{1}{2}, \qquad x_1 = x_0 + h = 0 + \frac{1}{2} = \frac{1}{2} \)

Applying Simpson’s rule:

\( \int^1_0 x^2 \space e^{-x} \space dx \approx \frac{1/2}{3}[0 + 0.6065308 + 0.3678794] = 0.1624017 \)

The integral over the interval \([0, 1]\) is \(2 – \frac{5}{e} = 0.1606028\), thus the error is equal to 0.0017989. The approximation resulting from the trapezoidal rule is \(0.1839397\), giving an error of 0.0233369.

Using the simpsons.rule() function we created earlier to arrive at the same result:

f2 <- function(x) {
  return(x^2 * exp(-x))
}

simpsons.rule(f2, 0, 1)
## [1] 0.1624017

Composite Simpson’s Rule

As with other Newton-Cotes formulas for approximating definite integrals, Simpson’s rule becomes inaccurate as the integration interval increases and if the nodes are unequally spaced. Similar to the composite trapezoidal rule, the composite approach of Simpson’s rule involves separating the integration interval into an even integer \(n\) (the number of subintervals, \(n\) in the composite Simpson’s rule approach is required to be even, unlike the composite trapezoidal rule), and then summing each subinterval. More formally, let \(h = (b – a) / n\) and \(x_j = a + jh\) for \(j = 0, 1, 2, \cdots, n\). The composite Simpson’s rule is then defined as:

\( \int_a^b f(x) \space dx = \frac{h}{3} \Bigg[f(a) + 2 \sum^{(n/2)-1}_{j=1} f(x_{2j}) + 4 \sum^{n/2}_{j=1} f(x_{2j – 1}) + f(b) \Bigg] – \frac{b – a}{180} h^f f^{(4)} (\mu) \)

Expanding the summations in Simpson’s rule:

\( \int_a^b f(x) \space dx = \frac{h}{3} \Bigg[f(a) + 4f(x_1) + 2f(x_2) + 4f(x_3) + 2f(x_4) + \cdots + 4f(x_{n-1}) + f(b) \Bigg] \)

Where \(a \leq \mu \leq b\). Let’s say we are interested in approximating the following integral using the Composite Simpson’s rule with \(n = 8\) subintervals:

\( \int_0^2 e^{2x} \space sin \space 3x \space dx \)

Approximating the integral with the Composite Simpson’s rule proceeds as follows:

\( \int_0^2 e^{2x} \space sin \space 3x \space dx = \frac{1}{12} \Bigg[e^{2(2)} \space sin \space 3(2) + 2 \sum^{3}_{j=1} e^{2x_{2j}} \space sin \space 3(x_{2j}) + 4 \sum^4_{j=1} e^{2x_{2j-1}} \space sin \space 3(x_{2j-1}) + e^{2(0)} \space sin \space 3(0) \Bigg] = -14.18334 \)

The following is an implementation of the Composite Simpson’s rule.

composite.simpson <- function(f, a, b, n) {
  if (is.function(f) == FALSE) {
    stop('f must be a function with one parameter (variable)')
  }
  
  h <- (b - a) / n
  
  xj <- seq.int(a, b, length.out = n + 1)
  xj <- xj[-1]
  xj <- xj[-length(xj)]

  approx <- (h / 3) * (f(a) + 2 * sum(f(xj[seq.int(2, length(xj), 2)])) + 4 * sum(f(xj[seq.int(1, length(xj), 2)])) + f(b))
  
  return(approx)
  
}

Replicate the function in R.

f3 <- function(x) {
  return(exp(2 * x) * sin(3 * x))
}

Approximate the definite integral with the composite.simpson() function.

composite.simpson(f3, 0, 2, 8)
## [1] -14.18334

References

Burden, R. L., & Faires, J. D. (2011). Numerical analysis (9th ed.). Boston, MA: Brooks/Cole, Cengage Learning.

Weisstein, Eric W. “Simpson’s Rule.” From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/SimpsonsRule.html

The post Simpson’s Rule for Approximating Definite Integrals in R appeared first on Aaron Schlegel.