##// END OF EJS Templates
some fixes
some fixes

File last commit:

r6035:3077781f
r7449:2e4aae94
Show More
trapezoid_rule.ipynb
117 lines | 21.0 KiB | text/plain | TextLexer

Basic numerical integration: the trapezoid rule

A simple illustration of the trapezoid rule for definite integration:

$$ \int_{a}^{b} f(x)\, dx \approx \frac{1}{2} \sum_{k=1}^{N} \left( x_{k} - x_{k-1} \right) \left( f(x_{k}) + f(x_{k-1}) \right). $$
First, we define a simple function and sample it between 0 and 10 at 200 points
In [1]:
def f(x):
    return (x-3)*(x-5)*(x-7)+85

x = linspace(0, 10, 200)
y = f(x)

Choose a region to integrate over and take only a few points in that region

In [2]:
a, b = 1, 9
xint = x[logical_and(x>=a, x<=b)][::30]
yint = y[logical_and(x>=a, x<=b)][::30]

Plot both the function and the area below it in the trapezoid approximation

In [3]:
plot(x, y, lw=2)
axis([0, 10, 0, 140])
fill_between(xint, 0, yint, facecolor='gray', alpha=0.4)
text(0.5 * (a + b), 30,r"$\int_a^b f(x)dx$", horizontalalignment='center', fontsize=20);
No description has been provided for this image

Compute the integral both at high accuracy and with the trapezoid approximation

In [4]:
from scipy.integrate import quad, trapz
integral, error = quad(f, 1, 9)
print "The integral is:", integral, "+/-", error
print "The trapezoid approximation with", len(xint), "points is:", trapz(yint, xint)
The integral is: 680.0 +/- 7.54951656745e-12
The trapezoid approximation with 6 points is: 621.286411141
In [ ]: