Integration Exercise 1

Imports



In [38]:

    
%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 [39]:

    
def trapz(f, a, b, N):
    """Integrate the function f(x) over the range [a,b] with N points."""
    # YOUR CODE HERE
    area=0
    sub_area=0
    h=((b-a)/N)
    for k in range (N+1):
        sub_area=0.5*h*((f(a+(k-1)*h)+f(a+k*h)))
        area = area+sub_area
    return area
f= lambda x: x**2

trapz(f,0,1,1000)









    Out[39]:





0.3333335004999998



In [40]:

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



In [42]:

    
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)



In [ ]:

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 [43]:

    
resultf,errorf = integrate.quad(f,0,1)
resultg,errorg = integrate.quad(g,0,np.pi)

print(resultf, errorf, resultg, errorg)









    



0.33333333333333337 3.700743415417189e-15 2.0 2.220446049250313e-14



In [ ]:

    
assert True # leave this cell to grade the previous one