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.

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.
- Simpson's 1/3 rule uses quadratic interpolation on panels.
- Simpson's 3/8 rule uses cubic interpolation on panels.

Humor



What did the finger say to the thumb?

I'm in glove with you.

Midpoint Rule

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} } \]

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} } \]

Midpoint Rule with R

R program from book:

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

Midpoint Rule: Code Check

a <- 0; b <- 3; m <- 3
h = (b-a)/m 
(x <- seq(a,b-h,length=m))

[1] 0 1 2

(x <- x + h/2 )

[1] 0.5 1.5 2.5

f<-function(x){x^2}
(y <- f(x))

[1] 0.25 2.25 6.25

sum(y)*abs(b-a)/m

[1] 8.75

Example 1: Exact Value of Integral

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

Example 1: Midpoint Rule

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.

Trapezoid Rule

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

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

\[ \small{ T_m = \frac{h}{2}\left[ f(x_0) + 2f(x_1) + 2f(x_2) + \ldots + 2f(x_{m-1}) + f(x_m) \right]} \]

Trapezoid Rule with R

trap <- function (f, a, b, m = 100) {
x <- seq(a,b,length = 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 )}

\[ \small{ T_m = \frac{b-a}{2n}\left[ f(x_0) + 2f(x_1) + 2f(x_2) + \ldots + 2f(x_{m-1}) + f(x_m) \right]} \]

Trapezoid Rule: Check Code Chunk

a <- 0
b <- 6
m <- 6
h = (b-a)/m 
(x<-seq(a,b,length=m+1))

[1] 0 1 2 3 4 5 6

f <- function(x){x^2}

(y1 <- f(x[2:(m+1)]))

[1]  1  4  9 16 25 36

(y2 <- f(x[1:m]))

[1]  0  1  4  9 16 25

(y3 <- y1+y2 )

[1]  1  5 13 25 41 61

Example 2: Exact Value of Integral

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

Example 2: Trapezoid Rule

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

Example 3: Trapezoid vs Midpoint

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.

Simpson's 1/3 Rule: Formula

\[ \small{ \int_a^b f(x)dx \cong S_m } \]

\[ \small{ \begin{aligned} S_m & = \left[ f(x_1) + 4f(x_2) + 2f(x_3) + 4f(x_4) + 2f(x_5) + \cdots + 4f(x_{m}) + f(x_{m+1}) \right] h \\ h &= \frac{b-a}{3m} \end{aligned} } \]

Simpson's Rule with R

Program in book:

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 )}

\[ \small{ S_m = \left[ f(x_1) + 4f(x_2) + 2f(x_3) + 4f(x_4) + 2f(x_5) + \cdots + 4f(x_{m}) + f(x_{m+1}) \right] h } \]

Simpson's Rule: Code Check

a <- 0; b <- 8; m <- 8; f <- function(x){2*x+1}

\[ \small{ \left[ f(0) + 4f(1) + 2f(2) + 4f(3) + 2f(4) + 4f(5) + 2f(6) + 4f(7) + f(8) \right] h } \]

Simpson's Rule: Check Code Chunk 1

a <- 0; b <- 8; m <- 8; f <- function(x){2*x+1}

(x.ends <- seq(a,b,length.out = m + 1))

[1] 0 1 2 3 4 5 6 7 8

(y.ends <- f(x.ends))

[1]  1  3  5  7  9 11 13 15 17

\[ \small{ f(0) + 4f(1) + 2f(2) + 4f(3) + 2f(4) + 4f(5) + 2f(6) + 4f(7) + f(8) } \]

Simpson's Rule: Check Code Chunk 2

(x.ends <- seq(a,b,length.out = m + 1))

[1] 0 1 2 3 4 5 6 7 8

(x.mids <- (x.ends[2:(m+1)] - x.ends[1:m])/2)

[1] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5

(x.mids <- x.mids + x.ends[1:m])

[1] 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5

\[ \small{ f(0) + 4f(1) + 2f(2) + 4f(3) + 2f(4) + 4f(5) + 2f(6) + 4f(7) + f(8) } \]

Simpson's Rule: Check Code Chunk 3

(y.mids <- f(x.mids))

[1]  2  4  6  8 10 12 14 16

(y.ends <- f(x.ends))

[1]  1  3  5  7  9 11 13 15 17

(A <- sum(y.ends[2:(m+1)] + 4*y.mids[1:m] 
             + y.ends [1:m])*abs(b-a)/(6*m))

[1] 72

\[ \small{ f(0) + 4f(1) + 2f(2) + 4f(3) + 2f(4) + 4f(5) + 2f(6) + 4f(7) + f(8) } \]

Simpson's Rule: Code Check

(A <- sum(y.ends[2:(m+1)] + 4*y.mids[1:m] 
             + y.ends [1:m])*abs(b-a)/(6*m))

[1] 72

\[ \small{ \left[ f(0) + 4f(1) + 2f(2) + 4f(3) + 2f(4) + 4f(5) + 2f(6) + 4f(7) + f(8) \right] h } \]

Example 4: Exact Value of Integral

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

Example 4: Simpson's Rule

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

Example 5: Simpson's Rule

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

Simpson's 3/8 Rule with R

Simpson's Rule uses a quadratic interpolating polynomial.
Simpson's 3/8 Rule uses a cubic interpolating polynomial.

simp38 <- function (f, a, b, m = 100) {
x.ends = seq(a,b,length = 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) }

Example 6: Simpson's 3/8 Rule

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

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

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

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

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} } \]