Error bounds for numerical integrals

Suppose we wish to estimate

to an accuracy of 0.00001. In this problem, well work through that process for several choices of summation procedure. For the the computations, feel free to use code in our Basic Integration notebook.

Left sum

• Find a value n large enough so that the corresponding left hand sum is within our desired accuracy level.
• Type out the resulting sum in mathematical notation.
• Evaluate the resulting sum in Python.

Trapezoidal sum

• Find a value n large enough so that the corresponding Trapezoidal sum is within our desired accuracy level.
• Type out the resulting sum in mathematical notation.
• Evaluate the resulting sum in Python.

Simpson’s rule

• Find a value n large enough so that the corresponding sum provided by Simpson’s rule is within our desired accuracy level.
• Type out the resulting sum in mathematical notation.
• Evaluate the resulting sum in Python.

Gaussian Quadrature

Approximate the integral using 3-point Gaussian quadrature.

SciPy

For crying out loud, just compute the integral already using scipy.integrate.quad.

For the left hand sum we know the max error would be m_1\dfrac{(b-a)^2}{n} where m_1 is the largest value our function can take over the interval [a,b], for f(x) = e^{-x^2} the largest value would be 1, thus for this problem m_1 = 1. Now if we want an accuracy of 0.00001 we just need to solve \dfrac{4}{n} = 0.0001. Thus n = 400000.

Now we just need to sum \sum_{i=0}^{n} f(x_i)\Delta x, in python we will use

def left_sum(f,a,b,n):
   xs = np.linspace(a,b,n+1)
   ys = f(-xs**2)
   return sum([ys[i] for i in range(n)])*(b-a)/(n)


 left_sum(np.exp, 0, 2, 400001)

and we get 0.88208384496707015

Trapezoidal Sum
We actually need n = 50

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import quad, dblquad, odeint
from time import time

import warnings
warnings.filterwarnings(‘ignore’)

def trapezoidal_sum(f,a,b,n):
xs = np.linspace(a,b,n+1)
ys = f(xs)
return sum([(ys[i+1]+ys[i]) for i in range(n)])(b-a)/(2n)
def f(x):
return np.exp(-1*(x**2))

truevalue = quad(f, 0,2)
print(truevalue)
bestN = 0
foundValue = False
for i in range(1000):
i = i+1
y = trapezoidal_sum(f,0,2,i)
error = np.abs(y-truevalue[0])
if error<10**-5 and foundValue == False:
bestN = i
print(“bestN is”,i)
foundValue = True
print(bestN)
y = trapezoidal_sum(f,0,2,50)
print(y)
print(np.abs(y-truevalue))

Executing this code we get the following printed
TrueValue =(0.8820813907624215, 9.793070696178202e-15)
bestN is 50
50
Estimate = 0.882071625027
Error [ 9.76573581e-06 8.82071625e-01]

Here the first line is the true value of the integral
the second line shows n=50 is the first value where the error of trapezoidal rule estimation drops below 10^-5, the third line shows the computed value using the trapezoidal rule with this n. The last line shows the error, which is just under 10^-5. We have parts a) and c) answered. To write this in mathematical notation

\frac{2}{50} \sum_{k=1}^{50} \frac{e^{-x_{k-1}^2} +e^{-x_{k}^2}}{2}

Where x_k = k\frac{2}{50}

We can simplify this to \frac{1}{50}\sum_{k=1}^{50} e^{-x_{k-1}^2} +e^{-x_{k}^2}

For accuracy of 0.0001, n = 400,000 if 4/n=0.00001

Left sum

first I defined f(x)

  def f(x): 
      return np.exp(-x**2)

Then I used the following code for the left sum

 def left_sum(f,a,b,n):
     xs = np.linspace(a,b,n+1)
     ys = f(xs)
     return sum([ys[i] for i in range(n)])*(b-a)/(n)
  left_sum(f, 0, 2, 400000)

to get 0.88208384497316439

Trapezoidal sum:

 def trapezoidal_sum(f,a,b,n):
     xs = np.linspace(a,b,n+1)
     ys = f(xs)
     return sum([(ys[i+1]+ys[i]) for i in range(n)])*(b-a)/(2*n)
 trapezoidal_sum(f, 0, 2, 400000)

to get 0.88208139076227221

Simpson’s rule:

 def simps(f,a,b,n):
     xs = np.linspace(a,b,n+1)
     ys = f(xs)
      return sum([(ys[i]+4*ys[i+1]+ys[i+2]) for i in range(0,n,2)])*(b-a)/(3*n)
 simps(f, 0, 2, 400000)

to get 0.88208139076242131

funmanbobyjo, i am confused. You solved for n to be 400,000 but then your python code only has it make 100 rectangles still. Or am I reading something wrong?