Ch5.2 Newton-Cotes Integration

Introduction

  • Newton-Cotes integration methods estimate the area under a curve \( f(x) \) by using panels of successively smaller width.
  • For panels, rectangles and trapezoids are often used.

title

Introduction

  • These integration methods correspond to finding area under interpolating polynomials for values of \( f(x) \) on nodes \( x_k \) for each panel, and then adding up (composite Newton-Cotes).
    • Midpoint rule uses constant interpolation on panels.
    • Trapezoid rule uses linear interpolation on panels.

titletitle

Midpoint Rule

title

  • This method begins by finding midpoint of rectangular panels.
  • To compute area of panel, \( f \) is evaluated at midpoint to obtain height, then multiplied by width.
  • Areas of panels are added up to obtain an estimate of area under \( f(x) \) over \( [a,b] \).

Midpoint Rule

  • In the formula below, \( f \) is evaluated at panel midpoints, then multiplied by width \( h \), and added together.

\[ \small{ \int_a^b f(x)dx \cong \sum_{k=1}^m f \left( a + h \left(k - \frac{1}{2} \right) \right)h, \,\,\, h = \frac{b-a}{m} } \]

title

Midpoint Rule

\[ \small{ \begin{aligned} \int_a^b f(x)dx &\cong \sum_{k=1}^m f \left( a + k \frac{(b-a)}{m} - \frac{(b-a)}{2m} \right) \frac{(b-a)}{m} \\ & = \sum_{k=1}^m f \left( a + \frac{(b-a)}{m} \left(k - \frac{1}{2} \right) \right) \frac{(b-a)}{m} \\ & = \sum_{k=1}^m f \left( a + h \left(k - \frac{1}{2} \right) \right)h, \,\,\, h = \frac{b-a}{m} \end{aligned} } \]

title

R Code: Midpoint Rule

midpt <- function (f, a, b, m = 100) {
nwidth = (b - a) / m
x = seq(a, b - nwidth,length.out = m) + nwidth / 2
y = f(x)
area = sum(y) * abs(b - a) / m
return ( area )
}

Example: Exact Value of Integral

\[ \small{ \int_0^1 x^2 dx = \frac{1}{3}x^3 |_0^1 = \frac{1}{3} } \]

title

Ex 1: Midpoint Rule

title

f <- function(x) { x^2 }
A1<-midpt(f,0,1,m = 10)
A2<-midpt(f,0,1,m = 100)
(A <- c(A1,A2))

[1] 0.332500 0.333325

(E <- abs(1/3 - A2))

[1] 8.333333e-06

Comments

  • As n increases, approximation of integral grows both more precise and more accurate.
  • This increasing precision and accuracy comes at cost of addition computing power.
  • More complex integrals will require more computing power to achieve less than desirable results.

title

Trapezoid Rule

  • Use formula for area of trapezoid for each panel, then add up.

\[ \small{ \begin{aligned} \int_a^b f(x)dx &\cong \sum_{k=1}^m \frac{f(c_k)+f(c_{k+1}) }{2}h \\ c_k & = a + h k, \,\,\, h = \frac{b-a}{m} \end{aligned} } \]

title

R Code: Trapezoid Rule

trap <- function (f, a, b, m = 100) {
x = seq(a, b, length.out = m + 1)
y = f(x)
p.area = sum((y[2:(m+1)] + y[1:m]))
p.area = p.area * abs(b - a) / (2 * m)
return (p.area )
}

Ex 2: Trapezoid Rule

title 

f <- function(x) { x^2 }
A1<-trap(f,0,1,m = 10)
A2<-trap(f,0,1,m = 100)
(A <- c(A1,A2))

[1] 0.33500 0.33335

(E <- abs(1/3 - A2))

[1] 1.666667e-05

Ex 3: Trapezoid vs Midpoint

title 

f <- function(x) { x^2 }
A1<-trap(f,0,1,m = 100)
A2<-midpt(f,0,1,m = 100)
c(A1,A2)

[1] 0.333350 0.333325

c(abs(1/3-A1),abs(1/3-A2))

[1] 1.666667e-05 8.333333e-06

Simpson's Rule

  • Simpson's Rule (1/3) uses a quadratic interpolating polynomial; requires odd number of nodes.
  • Simpson's 3/8 Rule uses a cubic interpolating polynomial.

title

Simpson's Rule: Formula

\[ \small{ \begin{aligned} \int_a^b f(x)dx &\cong \left[ f(a) + 4f(x_1) + 2f(x_2) + 4f(x_3) + 2f(x_4) + \cdots \right. \\ & \hspace{4in} \left. + 4f(x_{n-1}) + f(x_n) \right] h \\ h &= \frac{b-a}{3m} \end{aligned} } \]

title

R Code: Simpson's Rule

simp <- function (f, a, b, m = 100) {
x.ends = seq(a, b, length.out = m + 1)
y.ends = f(x.ends)
x.mids = (x.ends [2:( m + 1)] - x.ends [1:m]) / 2 + x.ends [1:m]
y.mids = f(x.mids )
p.area = sum(y.ends[2:(m+1)] + 4 * y.mids[1:m] + y.ends [1:m])
p.area = p.area * abs(b - a) / (6 * m)
return (p.area )
}

Ex 4: Simpson's Rule

title 

f <- function(x) { x^2 }
A1<-simp(f,0,1,m = 4)
A2<-simp(f,0,1,m = 8)
c(A1,A2)

[1] 0.3333333 0.3333333

c(abs(1/3-A1),abs(1/3-A2))

[1] 0 0

Ex 5: Simpson's Rule

title 

f <- function(x) { x^4 }
A1<-simp(f,0,1,m = 10)
A2<-simp(f,0,1,m = 100)
c(A1,A2)

[1] 0.2000008 0.2000000

c(abs(1/5-A1),abs(1/5-A2))

[1] 8.333333e-07 8.333334e-11

R Code: Simpson's 3/8 Rule

simp38 <- function (f, a, b, m = 100) {
x.ends = seq(a, b, length.out = m + 1)
y.ends = f(x.ends)
x.midh = (2 * x.ends[2:(m + 1)] + x.ends[1:m]) / 3
x.midl = (x.ends[2:(m + 1)] + 2 * x.ends[1:m]) / 3
y.midh = f(x.midh)
y.midl = f(x.midl)
p.area = sum(y.ends[2:(m+1)] + 3 * y.midh[1:m] 
             + 3*y.midl[1:m] + y.ends[1:m])
p.area = p.area * abs(b - a) / (8 * m)
return(p.area)
}

Ex 6: Simpson's 3/8 Rule

title 

f <- function(x) { x^4 }
A1<-simp38(f,0,1,m = 10)
A2<-simp38(f,0,1,m = 100)
c(A1,A2)

[1] 0.2000004 0.2000000

c(abs(1/5-A1),abs(1/5-A2))

[1] 3.703704e-07 3.703704e-11

Ex 7: Simpson 1/3 vs Simpson 3/8

title 

f <- function(x) { x^4 }
A1<-simp(f,0,1,m = 100)
A2<-simp38(f,0,1,m = 100)
c(A1,A2)

[1] 0.2 0.2

c(abs(1/5-A1),abs(1/5-A2))

[1] 8.333334e-11 3.703704e-11

5.2.1 Conclusion

  • Newton-Cotes integration methods estimate the area under a curve \( f(x) \) by using panels of successively smaller width.
    • Midpoint rule uses constant interpolation on panels.
    • Trapezoid rule uses linear interpolation on panels.
    • Simpson's 1/3 rule uses quadratic interpolation on panels.
    • Simpson's 3/8 rule uses cubic interpolation on panels.

title

Ch5.2.2: Error Formulas

  • The multipanel Newton-Cotes formulas use interpolating polynomials to fit data.
  • The integral of the polynomials provide the approximations to the actual integral.
  • The integral of error formula for LaGrange interpolating polynomials from Chapter 4 can be used to generate error formulas for numerical integrals.

\[ \small{ \begin{aligned} \int R_n(x) dx &= \int \frac{f^{n+1}(c)}{(n+1)!}(x-x_1)(x-x_2)\cdots(x-x_n) dx, \\ c &\in [x_1,x_n] \end{aligned} } \]

Error Formulas

  • The formulas below use \( h = b-a \).
  • Also, last formula is update to book formula.

\[ \small{ \begin{aligned} E_T & = \frac{h^3}{12m^2}f''(c) \\ E_M & = \frac{h^3}{24m^2}f''(c) \\ E_S & = \frac{h^5}{180m^4}f^{(4)}(c) \\ E_{S_{3/8}} & = \frac{h^5}{405m^4}f^{(4)}(c) \\ \end{aligned} } \]

Error Formulas: Discussion

  • All are exact for linear \( f(x) \).
  • Simpson formulas exact for cubic \( f(x) \).
  • Midpoint formula more accurate than trapezoid.
  • Cubic Simpson more accurate than quadratic Simpson.

\[ \small{ \begin{aligned} O(h^3): E_T &= \frac{h^3}{12m^2}f''(c), \,\,\, E_M = \frac{h^3}{24m^2}f''(c) \\ O(h^5): E_S & = \frac{h^5}{180m^4}f^{(4)}(c), \,\,\, E_{S_{3/8}} = \frac{h^5}{405m^4}f^{(4)}(c) \\ \end{aligned} } \]

Exact: Simpson's Rule

title 

f <- function(x) { x^2 }
A1<-simp(f,0,1,m = 4)
A2<-simp(f,0,1,m = 8)
c(A1,A2)

[1] 0.3333333 0.3333333

c(abs(1/3-A1),abs(1/3-A2))

[1] 0 0

Trapezoid vs Midpoint

title 

f <- function(x) { x^2 }
A1<-trap(f,0,1,m = 100)
A2<-midpt(f,0,1,m = 100)
c(A1,A2)

[1] 0.333350 0.333325

c(abs(1/3-A1),abs(1/3-A2))

[1] 1.666667e-05 8.333333e-06

Simpson 1/3 vs Simpson 3/8

title 

f <- function(x) { x^4 }
A1<-simp(f,0,1,m = 100)
A2<-simp38(f,0,1,m = 100)
c(A1,A2)

[1] 0.2 0.2

c(abs(1/5-A1),abs(1/5-A2))

[1] 8.333334e-11 3.703704e-11

5.2.2 Conclusion

  • Newton-Cotes methods exact for linear \( f(x) \).
  • Simpson formulas exact for cubic \( f(x) \).
  • Midpoint formula more accurate than trapezoid.
  • Cubic Simpson more accurate than quadratic Simpson.

\[ \small{ \begin{aligned} O(h^3): E_T &= \frac{h^3}{12m^2}f''(c), \,\,\, E_M = \frac{h^3}{24m^2}f''(c) \\ O(h^5): E_S & = \frac{h^5}{180m^4}f^{(4)}(c), \,\,\, E_{S_{3/8}} = \frac{h^5}{405m^4}f^{(4)}(c) \\ \end{aligned} } \]

Ch5.3: Gauss Quadrature

  • See Class Notes (pdf) for details of method.
  • Gauss quadrature achieves high accuracy for very few nodes.
    • Nodes and weights strategically computed.
    • Nodes and weights same for all \( f(x) \).
    • Linear and quadratic interpolation (\( n=2 \), \( n=3 \))
  • Greater accuracy and fewer flops than for Newton-Cotes.

title

R Code in Book for [a,b] = [0,1]

gaussint <- function(f, x, w) {
  y <- f(x)
  return (sum(y * w))
}

w <- c(1, 1)
x <- c(-1 / sqrt(3), 1 / sqrt(3))
f <- function(x) { x^2 }
gaussint(f, x, w)

[1] 0.6666667

Two-Point Gauss Quadrature on [a, b]

gaussquad2 <- function(f,a,b) {
  w <- c(1,1)
  x <- c(-1/sqrt(3),1/sqrt(3))
  d <- (b-a)/2
  s <- (b+a)/2
  u <- d*x + s
  y <- f(d*x + s)
  return(d*sum(y*w))
}

Example 8

title 

f <- function(x) { x^2 }
A1<-gaussquad2(f,0,1)
A2<-trap(f,0,1,m = 10)
c(A1,A2)

[1] 0.3333333 0.3350000

c(abs(1/3-A1),abs(1/3-A2))

[1] 0.000000000 0.001666667

Example 9

title 

f <- function(x) { x^4 }
A1<-gaussquad2(f,0,1)
A2<-trap(f,0,1,m = 10)
c(A1,A2)

[1] 0.1944444 0.2033300

c(abs(1/5-A1),abs(1/5-A2))

[1] 0.005555556 0.003330000

Three-Point Gauss Quadrature on [a, b]

gaussquad3 <- function(f,a,b) {
  w <- c(5/9,8/9,5/9)
  x <- c(-sqrt(3/5),0,sqrt(3/5))
  d <- (b-a)/2
  s <- (b+a)/2
  u <- d*x + s
  y <- f(d*x + s)
  return(d*sum(y*w))
}

Example 10

title 

f <- function(x) { x^2 }
A1<-gaussquad3(f,0,1)
A2<-simp(f,0,1,m = 10)
c(A1,A2)

[1] 0.3333333 0.3333333

c(abs(1/3-A1),abs(1/3-A2))

[1] 5.551115e-17 5.551115e-17

Example 11

title 

f <- function(x) { x^4 }
A1<-gaussquad3(f,0,1)
A2<-simp(f,0,1,m = 10)
c(A1,A2)

[1] 0.2000000 0.2000008

c(abs(1/5-A1),abs(1/5-A2))

[1] 0.000000e+00 8.333333e-07

Example 12

title 

f <- function(x){exp(-x^2)}
A1<-gaussquad3(f,1,1.5)
A2<-simp(f,1,1.5,m = 10)
c(A1,A2)

[1] 0.1093642 0.1093643

A <- 0.109364260812
c(abs(A-A1),abs(A-A2))

[1] 6.478000e-08 5.253018e-09

5.3 Conclusion

  • Gauss quadrature achieves high accuracy for very few nodes.
    • Nodes and weights strategically computed.
    • Nodes and weights same for all \( f(x) \).
    • Linear and quadratic interpolation (\( n=2 \), \( n=3 \))
  • Greater accuracy and fewer flops than for Newton-Cotes.

title