Ch19.5 Part 1: Newton-Cotes Numerical Integration

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

• 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.

You know, people say they pick their nose.

But I feel like I was just born with mine.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

#Define the function to integrate
#The second function below is taken from book. 
def f(x):                          
    #return 1.0*x**2 + 1.0         
    return np.exp(-x**2)          

#Left sum for n rectangles over [a,b]
#n+1 nodes gives n intervals and n rectangles
def leftsum(a,b,n):
    h = (b-a)/n            
    nodes=np.zeros(n+1)     
    yvalues=np.zeros(n+1)
    #Compute nodes and yvalues
    for k in range(0,n+1): #from 0 to n
        nodes[k] = a + k*h     
        yvalues[k] = f(nodes[k])   
    #Initialize S and add up yvalues
    S = 0                    
    for k in range(0,n):    
        S = S + yvalues[k]     
    #Left sum value
    S = h*S        

#These lines also indented 4 spaces
#as part of midsum definition.

#Compute Percent Relative Error (PRE)
Q = quad(f,a,b) #[scipy value, error estimate]
Q = Q[0]  #scipy value only
PRE = (Q - S)/Q*100  #PRE
print('Q = %.6f, S = %.6f, PRE = %.6f' % (Q, S,PRE))  

#Plot commands
plt.plot(nodes,yvalues,'-o')   
plt.xlabel('x')
plt.legend(('f'),loc = 0)
plt.show() 

import numpy as np
import matplotlib.pyplot as plt

def leftsum(a,b,n):
    #Left sum for n rectangles over [a,b]
    #n+1 nodes gives n intervals and n rectangles
    f = lambda x: np.exp(-x**2) #This is f(x)
    nodes=np.linspace(a,b,n+1) #n+1 nodes
    xn = nodes[0:n]    #n left nodes 
    fn = f(xn)         #left sum y-values
    h = (b-a)/n        #interval width
    S = sum(fn)*h      #vectorized sum
    print(S)
    plt.plot(xn,fn,'-o')   
    plt.xlabel('x')
    plt.legend(('f'),loc = 0)
    plt.show() 

#Midpoint sum for n rectangles over [a,b]
#n+1 nodes gives n intervals and n rectangles
def midsum(a,b,n):
    h = (b-a)/n         #interval width        
    nodes=np.zeros(n+1) #initalize nodes  
    yvalues=np.zeros(n+1)
    #Compute nodes and yvalues
    for k in range(0,n+1):  #from 0 to n
        nodes[k] = a + (k+1/2)*h     
        yvalues[k] = f(nodes[k])   
    #Initialize S and add up yvalues
    S = 0                    
    for k in range(0,n):    
        S = S + yvalues[k]     
    #Midpoint sum value
    S = h*S   

import numpy as np
import matplotlib.pyplot as plt

def midsum(a,b,n):
    #Midpoint sum for n rectangles over [a,b]
    f = lambda x: np.exp(-x**2) #This is f(x)
    N = 2*n+1           #Need this for midpoints
    nodes=np.linspace(a,b,N) #2n+1 nodes in [a,b]
    xn = nodes[1:N:2]   #n midpts (step size = 2) 
    fn = f(xn)          #midpoint y-values
    h = (b-a)/n         #interval width
    S = sum(fn)*h       #vectorized sum
    print(S)
    plt.plot(xn,fn,'-o')   
    plt.xlabel('x')
    plt.legend(('f'),loc = 0)
    plt.show() 

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

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

\[ \small{ T_n = \frac{h}{2}\left[ f(x_0) + 2f(x_1) + 2f(x_2) + \ldots + 2f(x_{n-1}) + f(x_n) \right]} \]

#Simpson's Rule for n intervals over [a,b]
#n+1 nodes gives n intervals; n must be even

#Compute nodes and yvalues
for k in range(0,n+1):  #from 0 to n
    nodes[k] = a + k*h     
    yvalues[k] = f(nodes[k])   
#Prepare summations
S4 = 0                    
for k in range(1,n,2):    
    S4 = S4 + 4*yvalues[k]   
S2 = 0                    
for k in range(2,n-1,2):    
    S2 = S2 + 2*yvalues[k]   
#Simpson Rule value
S = S4 + S2 + yvalues[0] + yvalues[n]
S = S*(h/3)    

import numpy as np
import matplotlib.pyplot as plt 

def Simp(a,b,n):
    #Simpson's Rule for n intervals over [a,b]
    #n+1 nodes gives n intervals; n must be even
    f = lambda x: np.exp(-x**2) #This is f(x)
    xn = np.linspace(a,b,n+1) #n+1 nodes
    fn = f(xn)         #y-values at nodes
    h = (b-a)/n        #interval width
    S =  4*sum(fn[1:n:2]) + 2*sum(fn[2:n-1:2]) 
    S = ( S + fn[0] + fn[n])*(h/3) 
    print(S)
    plt.plot(xn,fn,'-o')   
    plt.xlabel('x')
    plt.legend(('f'),loc = 0)
    plt.show() 

• 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}{12n^2}f''(c), \,\,\, E_M = \frac{h^3}{24n^2}f''(c) \\ O(h^5): E_{S_{1/3}} & = \frac{h^5}{180n^4}f^{(4)}(c), \,\,\, E_{S_{3/8}} = \frac{h^5}{405n^4}f^{(4)}(c) \\ \end{aligned} } \]