In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy import integrate
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 [12]:
def trapz(f, a, b, N):
"""Integrate the function f(x) over the range [a,b] with N points."""
h = (b-a)/N
I = 0
for i in range(N):
I += ((f(a+(b-a)*i/N)+f(a+(b-a)*(i+1)/N))/2)*h
print(I)
return I
In [13]:
f = lambda x: x**2
g = lambda x: np.sin(x)
In [14]:
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 [17]:
L = integrate.quad(f, 0,1)
K = integrate.quad(g, 0, np.pi)
print(L,K)
In [18]:
assert True # leave this cell to grade the previous one