Programming Numerical Integration with Python (and Javascript)

Numerically integrating the area under the curve using four trapezoids.

I gave a quick introduction to programming for my calculus class, which has been working on numerical integration.

Numerical integration is usually used for functions that can’t be integrated (or not easily integrated) but for this example we’ll use a simple parabolic function so we can compare the numerical results to the analytical solution (as seen here).

With the equation:

f(x) = -\frac{1}{4} x^2 + x + 4

To find the area under the curve between x = 1 and x = 5 we’d find the definite integral:

Area = \int_{_1}^{^5} \left(-\frac{x^2}{4} + x + 4 \right) \,dx

which gives the result:

Area = 17 \frac{2}{3} = 17.6\bar{6}

For numerical integration, we break the area of concern into a number of trapezoids, find the areas of all the trapezoids and add them up.

We’ll define the left and right boundaries of the area as a and b, and we can write the integral as:

Area = \int_{_a}^{^b} f(x) \,dx

The left and right boundaries of the area we’re interested in are defined as a and b respectively. The area of each trapezoid is defined as An.

We also have to choose a number of trapezoids (n) or the width of each trapezoid (dx). Here we choose four trapezoids (n = 4), which gives a trapezoid width of one (dx = 1).

The area of the first trapezoid can be calculated from its width (dx) and the height of the two upper ends of the trapezoid (f(x0) and f(x1).

So if we define the x values of the left and right sides of the first trapezoids as x0 and x1, the area of the first trapezoid is:

A_1 = \frac{f(x_{_0})+f(x_{_1})}{2} dx

For this program, we’ll set the trapezoid width (dx) and then calculate the number of trapezoids (n) based on the width and the locations of the end boundaries a and b. So:

n = \frac{b-a}{dx}

and the sum of all the areas will be:

\displaystyle\sum\limits_{i=1}^{n} \frac{f(x_{i-1})+f(x_{i})}{2} dx

We can also figure out that since the x values change by the same value (dx) for every trapezoid, it’s an arithmetic progression, so:

x_{i-1} = a + (i-1) dx

and,

x_{i} = a + i \cdot dx

so our summation becomes:

\displaystyle\sum\limits_{i=1}^{n} \frac{f(a+(i-1)dx)+f(a + i \cdot dx)}{2} dx

Which we can program with:

numerical_integration.py

# the function to be integrated
def func(x):
    return -0.25*x**2 + x + 4

# define variables
a = 1.          # left boundary of area
b = 5.          # right boundary of area
dx = 1          # width of the trapezoids

# calculate the number of trapezoids
n = int((b - a) / dx)

# define the variable for area
Area = 0

# loop to calculate the area of each trapezoid and sum.
for i in range(1, n+1):
    #the x locations of the left and right side of each trapezpoid
    x0 = a+(i-1)*dx
    x1 = a+i*dx

    #the area of each trapezoid
    Ai = dx * (func(x0) + func(x1))/ 2.

    # cumulatively sum the areas
    Area = Area + Ai

#print out the result.
print "Area = ", Area

And the output looks like

>>> 
Area =  17.5
>>>

While the programming is pretty straightforward, it was a bit of a pain getting Python to work for one of my students who is running Windows 8. I still have not figured out a way to get it to work properly, so I’m considering trying to do it using Javascript.

Update

The javascript functions for numerical integration:

function numerically_integrate(a, b, dx, f) {
	
	// calculate the number of trapezoids
	n = (b - a) / dx;
	
	// define the variable for area
	Area = 0;
	
	//loop to calculate the area of each trapezoid and sum.
	for (i = 1; i <= n; i++) {
		//the x locations of the left and right side of each trapezpoid
		x0 = a + (i-1)*dx;
		x1 = a + i*dx;
		
		// the area of each trapezoid
		Ai = dx * (f(x0) + f(x1))/ 2.;
		
		// cumulatively sum the areas
		Area = Area + Ai	
		
	} 
	return Area;
}

//define function to be integrated
function f(x){
	return -0.25*Math.pow(x,2) + x + 4;
}

// define variables
a = 1;		// left boundary of area
b = 5;		// right boundary of area
dx = 1;		// width of the trapezoids
	
// print out output
alert("Area = "+ numerically_integrate(a, b, dx, f));

This is a demonstration of a full html file that uses the function, and should work in any modern browser (download files: numerical-integration.zip).

Update 2

I’ve added the above javascript code to the embeddable graphs to allow it to calculate and display numerical integrals: you can change the values in the interactive graph below.

Leave a Reply