Exercise: trapezoidal integration

In this exercise, you are tasked with implementing the simple trapezoid rule formula for numerical integration. If we want to compute the definite integral

$$ \int_{a}^{b}f(x)dx $$

we can partition the integration interval $[a,b]$ into smaller subintervals. We then approximate the area under the curve for each subinterval by calculating the area of the trapezoid created by linearly interpolating between the two function values at each end of the subinterval:

For a pre-computed $y$ array (where $y = f(x)$ at discrete samples) the trapezoidal rule equation is:

$$ \int_{a}^{b}f(x)dx\approx\frac{1}{2}\sum_{i=1}^{n}\left(x_{i}-x_{i-1}\right)\left(y_{i}+y_{i-1}\right). $$

In pure python, this can be written as:

def trapz_slow(x, y):
    area = 0.
    for i in range(1, len(x)):
        area += (x[i] - x[i-1]) * (y[i] + y[i-1])
    return area / 2

Exercise 2

Part 1

Create two arrays $x$ and $y$, where $x$ is a linearly spaced array in the interval $[0, 3]$ of length 10, and $y$ represents the function $f(x) = x^2$ sampled at $x$.



In [1]:

    
import numpy as np

x = np.linspace(0, 3, 10)
y = x ** 2

print x
print y









    



[ 0.          0.33333333  0.66666667  1.          1.33333333  1.66666667
  2.          2.33333333  2.66666667  3.        ]
[ 0.          0.11111111  0.44444444  1.          1.77777778  2.77777778
  4.          5.44444444  7.11111111  9.        ]

Part 2

Use indexing (not a for loop) to find the 9 values representing $y_{i}+y_{i-1}$ for $i$ between 1 and 10.

Hint: What indexing would be needed to get all but the last element of the 1d array y. Similarly what indexing would be needed to get all but the first element of a 1d array.



In [2]:

    
y_roll_sum = y[:-1] + y[1:]
print y_roll_sum









    



[  0.11111111   0.55555556   1.44444444   2.77777778   4.55555556
   6.77777778   9.44444444  12.55555556  16.11111111]

Part 3

Write a function trapz(x, y), that applies the trapezoid formula to pre-computed values, where x and y are 1-d arrays. The function should not use a for loop.



In [3]:

    
def trapz(x, y):
    return 0.5 * np.sum((x[1:] - x[:-1]) * (y[:-1] + y[1:]))

Part 4

Verify that your function is correct by using the arrays created in #1 as input to trapz. Your answer should be a close approximation of $\int_0^3 x^2$ which is $9$.



In [4]:

    
trapz(x, y)









    Out[4]:





9.0555555555555554

Part 5 (extension)

numpy and scipy.integrate provides many common integration schemes. Find the documentation for NumPy's own version of the trapezoidal integration scheme and check its result with your own:



In [5]:

    
print np.trapz(y, x)









    



9.05555555556

Part 6 (extension)

Write a function trapzf(f, a, b, npts=100) that accepts a function f, the endpoints a and b and the number of samples to take npts. Sample the function uniformly at these points and return the value of the integral.

Use the trapzf function to identify the minimum number of sampling points needed to approximate the integral $\int_0^3 x^2$ with an absolute error of $<=0.0001$. (A loop is necessary here)



In [6]:

    
def trapzf(f, a, b, npts=100):
    x = np.linspace(a, b, npts)
    y = f(x)
    return trapz(x, y)

def x_squared(x):
    return x ** 2

abs_err = 1.0
n_samples = 0
expected = 9
while abs_err > 0.0001:
    n_samples += 1
    integral = trapzf(x_squared, 0, 3, npts=n_samples)
    abs_err = np.abs(integral - 9)

print 'Minimum samples for absolute error less than or equal to 0.0001:', n_samples









    



Minimum samples for absolute error less than or equal to 0.0001: 214