Using Ipython Notebook, Basic Numerical Integration: the Trapezoid Rule

Using Ipython Notebook,

Basic Numerical Integration: the Trapezoid Rule

In [1]:

1+1

Out[1]:

2

In [2]:

%matplotlib inline

import numpy as np

import matplotlib.pyplot as plt

In [3]:

def f(x):

    return (x-4)*(x-15)*(x-17)+85

x = np.linspace(0, 10, 200)

y = f(x)

In [4]:

a, b = 1, 14

xint = x[np.logical_and(x>=a, x<=b)][::30]

yint = y[np.logical_and(x>=a, x<=b)][::30]

In [5]:

plt.plot(x, y, lw=2)

plt.axis([0, 10, 0, 140])

plt.fill_between(xint, 0, yint, facecolor='gray', alpha=0.4)

plt.text(0.5 * (a + b), 30,r"$\int_a^b f(x)dx$", horizontalalignment='center', fontsize=20);

In [6]:

from __future__ import print_function

from scipy.integrate import quad, trapz

integral, error = quad(f, 1, 14)

print("The integral is:", integral, "+/-", error)

print("The trapezoid approximation with", len(xint), "points is:", trapz(yint, xint))

The integral is: 1875.25 +/- 3.5310617657e-11
The trapezoid approximation with 6 points is: 534.431021088

In [ ]: