Integration Exercise 1

Imports



In [27]:

    
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy import integrate

Trapezoidal rule

The trapezoidal rule generates a numerical approximation to the 1d integral:

$$ I(a,b) = \int_a^b f(x) dx $$

by dividing the interval $[a,b]$ into $N$ subdivisions of length $h$:

$$ h = (b-a)/N $$

Note that this means the function will be evaluated at $N+1$ points on $[a,b]$. The main idea of the trapezoidal rule is that the function is approximated by a straight line between each of these points.

Write a function trapz(f, a, b, N) that performs trapezoidal rule on the function f over the interval $[a,b]$ with N subdivisions (N+1 points).



In [28]:

    
def trapz(f, a, b, N):
    """Integrate the function f(x) over the range [a,b] with N points."""
    pts = np.linspace(a, b, N + 1)
    vals = f(pts)
    h = (b - a) / (1.0 * N)
    
    area = .5 * h * sum(vals[0:N] + vals[1:(N+1)])
    return area
    #raise NotImplementedError()



In [29]:

    
f = lambda x: x**2
g = lambda x: np.sin(x)



In [30]:

    
I = trapz(f, 0, 1, 1000)
assert np.allclose(I, 0.33333349999999995)
J = trapz(g, 0, np.pi, 1000)
assert np.allclose(J, 1.9999983550656628)

Now use scipy.integrate.quad to integrate the f and g functions and see how the result compares with your trapz function. Print the results and errors.



In [31]:

    
def compare(f, a, b, N):
    trapint = trapz(f, a, b, N)
    quadint = integrate.quad(f, a, b)[0]
    
    print("Trapezoid Rule: %f" % trapint)
    print("Scipy: %f" % quadint)
    print("Error: %f" % (quadint - trapint))

compare(f, 0, 1, 1000)
compare(g, 0, 1, 1000)
#raise NotImplementedError()









    



Trapezoid Rule: 0.333334
Scipy: 0.333333
Error: -0.000000
Trapezoid Rule: 0.459698
Scipy: 0.459698
Error: 0.000000



In [32]:

    
assert True # leave this cell to grade the previous one