Ch 5.4.2 Romberg's Method

Joseph Edwards
12/3/20

Introduction

  • Romberg's Method is an extension of Newton-Cotes Integration
  • Not adaptive (in most implementations), uses the trapezoidal rule to provide an integral estimate
  • Uses Richardson Extrapolation to significantly reduce the error of our estimate
    alt text

Richardson Extrapolation

  • Richardson extrapolation converts O(h) results into ones that are O(h2 ), O(h3 ), O(h4 ), etc
    alt text
  • Romberg Integration uses this idea to create better integral approximations using previous results

Recursive Implementation

alt text alt text

  • S is Simpson's Rule, T is the trapezoid function, f is the function being integrated, and m is the number of panels.\[ I_{0,0} \]is a one paneled trapezoid rule over the integral and\[ I_{j,0} \]is the trapezoid rule with 2j panels over the integral.
  • In R, the results are stored in a matrix

Implementation in R

romberg <- function(f,a,b,m,tab = FALSE) {
  R <- matrix(NA, nrow = m, ncol = m)

  R[1, 1] <- trap(f,a,b,m = 1)
  for (j in 2:m) {
    R[j,1] <- trap(f,a,b,m = 2^(j-1))
    for (k in 2:j) {
      k4 <- 4^(k-1)
      R[j,k] <- k4 * R[j,k-1] - R[j-1,k-1]
      R[j,k] <- R[j,k] / (k4-1)
    }
  }
  if (tab == TRUE) 
    return (R)
  return (R[m,m])
}

R Code

romberg(sin, 0, pi, m = 5, tab = TRUE)

             [,1]     [,2]     [,3]     [,4] [,5]
[1,] 1.923607e-16       NA       NA       NA   NA
[2,] 1.570796e+00 2.094395       NA       NA   NA
[3,] 1.896119e+00 2.004560 1.998571       NA   NA
[4,] 1.974232e+00 2.000269 1.999983 2.000006   NA
[5,] 1.993570e+00 2.000017 2.000000 2.000000    2
  • Each value not in the leftmost column is calculated using the value to the left, and the one above that. As each approximation is calculated, the error becomes smaller and smaller

R Code Continued

f <- function(x) {sin(x)^2 + log(x)}
romberg(f, 1, 10, m = 3)

[1] 18.48497

romberg(f, 1, 10, m = 5)

[1] 18.52473

romberg(f, 1, 10, m = 10)

[1] 18.52494

Conclusion/Comments

  • Romberg's Method is a very powerful application of Richardson Extrapolation for approximating integrals with finite domains
  • It does have limitations, such as the function needing to be evaluable at equally spaced points with defined endpoints
  • Other methods, such as the midpoint rule, can be used in place of the trapezoidal rule to avoid this issue.
    alt text