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.

Atom Builder

This app lets you drag and drop electrons, protons, and neutrons to create atoms with different charges, elements, and atomic masses. You can also enter the element symbol, charge and atomic mass and it will build the atom for you.

Note, however, it only does the first 20 elements.

Learning the Periodic Table: A Prototype

My middle school class is about to cover some very basic chemistry so I’ve asked them to memorize the first 20 elements in their correct order on the periodic table. To help, I’ve put together this interactive exercise where they drag an icon of the element to its correct place on the table. It says the name of the element whenever you start dragging a tile with the symbol. It’s also timed so students can quantify and compare how good they are.

Prototype exercise for learning the first 20 elements of the periodic table.

In this first prototype the elements are presented in order, but I figure that additional levels could have:

  • The elements come up at random (done).
  • Have the elements come up by vertical column (group) (done).
  • Instead of tiles with the elements, have diagrams with their electron configuration.
  • The program say the name of the element and the student has to click on the right cell in the table.

This was put together using HTML5 and Javascript. KineticJS was particularly useful. It should, in theory, work in any browser (but I have only tested it in Firefox and Google Chrome) and on touch-screen tablets as well.

Programming Agents

An artful spaceship.

For my programming elective, I thought I’d try to get students programming autonomous agents that we could all compile together into some type of game. Since it was, for most of them a first class in programming, this has proven a little too ambitious, and I’ve had to adjust a bit to build up to it.

Step 1: Build a spaceship/vehicle/agent.

A simple ship made of a sphere, ring, and pyramid, with a clear front end (pointed to the left: vector=(1,0,0).
  • We’re using VPython so we can build 3d ships. Students get to exercise their creativity a bit and learn how to place things in 3d co-ordinate space.
  • The ship must have a front end and all the parts must be put into a single frame so the whole thing can be moved as a unit.
  • Although it’s not required, I’ve added in the beginnings of a trail so we can track the motion of the spaceship
from visual import *

f = frame()
r = ring(frame=f, thickness=.25, axis=(0,1,0))
cabin = sphere(frame=f, radius=0.6, color=color.red)
front = pyramid(frame=f, height=1,width=1,length=2, color=color.blue)
trail = curve(pos=[f.pos,])

Step 2: Make the spaceship a class.

  • Putting the spaceship in a separate class, away from the rest of the code, makes it easier to move from one program to another. That means we can import everyone’s ships into a single program at the end.
from visual import *

class spaceship1:

    def __init__(self):

        self.f = frame()
        self.r = ring(frame=self.f, thickness=.25, axis=(0,1,0))
        self.cabin = sphere(frame=self.f, radius=0.6, color=color.red)
        self.front = pyramid(frame=self.f, height=1,width=1,length=2, color=color.blue)
        self.trail = curve(pos=[self.f.pos,])

ss = spaceship1()

Make the spaceship move

  • We’re going to control the movement of the spaceship using the arrow keys. Left and right to turn; up and down to accelerate and decelerate.
  • To allow movement, we’ll put the keyboard controls into an infinite loop, but to allow maximum flexibility (and to show how to use create methods in a class) we’re putting the actual movement in as methods in the class: the turning method is named turn_xz and movement is move_xz since all the motion is going to be restricted in the xz plane for now.
from visual import *

class spaceship1:

    def __init__(self):

        self.f = frame()
        self.r = ring(frame=self.f, thickness=.25, axis=(0,1,0))
        self.cabin = sphere(frame=self.f, radius=0.6, color=color.red)
        self.front = pyramid(frame=self.f, height=1,width=1,length=2, color=color.blue)
        self.trail = curve(pos=[self.f.pos,])

    def turn_xz(self, turn_rate):
        self.f.axis = rotate(self.f.axis, turn_rate, (0,1,0))

    def move_xz(self, speed):
        self.f.pos = self.f.pos + self.f.axis*speed
        self.trail.append(pos=self.f.pos)
                          
ss = spaceship1()
turn = 0.0
speed=0.0

while 1:
    rate(20)  #slows down execution of the loop

    ss.turn_xz(turn)
    ss.move_xz(speed)

    if scene.kb.keys: # event waiting to be processed?
        s = scene.kb.getkey() # get keyboard info

        if s == 'left':
            turn = turn + 0.01
        if s == 'right':
            turn = turn - 0.01

        if s == 'up':
            speed = speed + 0.01
        if s == 'down':
            speed = speed - 0.01

        if s == ' ':
            if turn <> 0:
                turn = 0
            elif speed <> 0:
                speed = 0

Fire missiles

  • We add in missiles by creating a class very similar to the spaceship. For now our missile is just a ball but I’m putting it into a frame anyway in case later on I want it to be a composite object.
  • Since we can fire multiple missiles, in the code I create a list to hold all the missiles.
  • Missiles are fired using the “a” key.
from visual import *

class spaceship1:

    def __init__(self):

        self.f = frame()
        self.r = ring(frame=self.f, thickness=.25, axis=(0,1,0))
        self.cabin = sphere(frame=self.f, radius=0.6, color=color.red)
        self.front = pyramid(frame=self.f, height=1,width=1,length=2, color=color.blue)
        self.trail = curve(pos=[self.f.pos,])

    def turn_xz(self, turn_rate):
        self.f.axis = rotate(self.f.axis, turn_rate, (0,1,0))

    def move_xz(self, speed):
        self.f.pos = self.f.pos + self.f.axis*speed
        self.trail.append(pos=self.f.pos)
          
class missile:
    def __init__(self, axis, pos):

        self.f = frame(axis=axis, pos=pos)
        self.r = sphere(frame=self.f, radius=0.2, color=color.green)
        self.speed = 0.5
        

    def move_xz(self):
        self.f.pos = self.f.pos + self.f.axis*self.speed

                
ss = spaceship1()
turn = 0.0
speed=0.0

missiles = []  #list for missiles

while 1:
    rate(20)  #slows down execution of the loop

    ss.turn_xz(turn)
    ss.move_xz(speed)

    #move missiles
    for i, m in enumerate(missiles):
        m.move_xz()

    if scene.kb.keys: # event waiting to be processed?
        s = scene.kb.getkey() # get keyboard info

        if s == 'left':
            turn = turn + 0.01
        if s == 'right':
            turn = turn - 0.01

        if s == 'up':
            speed = speed + 0.01
        if s == 'down':
            speed = speed - 0.01

        if s == ' ':
            if turn <> 0:
                turn = 0
            elif speed <> 0:
                speed = 0

        #fire missile
        if s == 'a':
            missiles.append(missile(ss.f.axis, ss.f.pos))

Final adjustments

  • Finally, I’m going to put in a boundary, and set it up to delete the missiles if they go out of bounds.
  • I’m also going to set an option so that the scene follows the ship, which makes it easier to keep track of things.
from visual import *

class spaceship1:

    def __init__(self):

        self.f = frame()
        self.r = ring(frame=self.f, thickness=.25, axis=(0,1,0))
        self.cabin = sphere(frame=self.f, radius=0.6, color=color.red)
        self.front = pyramid(frame=self.f, height=1,width=1,length=2, color=color.blue)
        self.trail = curve(pos=[self.f.pos,])

    def turn_xz(self, turn_rate):
        self.f.axis = rotate(self.f.axis, turn_rate, (0,1,0))

    def move_xz(self, speed):
        self.f.pos = self.f.pos + self.f.axis*speed
        self.trail.append(pos=self.f.pos)
                          

class missile:
    def __init__(self, axis, pos):

        self.f = frame(axis=axis, pos=pos)
        self.r = sphere(frame=self.f, radius=0.2, color=color.green)
        self.speed = 0.5
        

    def move_xz(self):
        self.f.pos = self.f.pos + self.f.axis*self.speed

class bounds:
    def __init__(self, boundary_distance):
        self.boundary_distance = boundary_distance
        self.border = curve(pos=[(-boundary_distance,0,-boundary_distance),
                                 (boundary_distance,0,-boundary_distance),
                                 (boundary_distance,0,boundary_distance),
                                 (-boundary_distance,0,boundary_distance),
                                 (-boundary_distance,0,-boundary_distance)])
    def in_bounds(self, pos):
        if ( pos.x < -self.boundary_distance or
             pos.x > self.boundary_distance or
             pos.z < -self.boundary_distance or
             pos.z > self.boundary_distance):
            return false
        else:
            return true
        

ss = spaceship1()
turn = 0.0
speed=0.0
l_track_ship = 1
scene.range=10
scene.forward = (0,-1,0)
missiles = []

boundary_distance = 100
boundary = bounds(boundary_distance)

while 1:
    rate(20)

    if l_track_ship > 0:
        scene.center = ss.f.pos
    
    ss.turn_xz(turn)
    ss.move_xz(speed)

    #move missiles
    del_list = []
    for i, m in enumerate(missiles):
        m.move_xz()

        #delete missile if out of bounds
        if boundary.in_bounds(m.f.pos) == False:
            del_list.append(i)
    for i in del_list:
        missiles[i].f.visible = False
        del missiles[i]
            
    
    #if m <> None:
    #    m.move_xz()

    if scene.kb.keys: # event waiting to be processed?
        s = scene.kb.getkey() # get keyboard info

        if s == 'left':
            turn = turn + 0.01
        if s == 'right':
            turn = turn - 0.01

        if s == 'up':
            speed = speed + 0.01
        if s == 'down':
            speed = speed - 0.01

        if s == ' ':
            if turn <> 0:
                turn = 0
            elif speed <> 0:
                speed = 0

        #fire missile
        if s == 'a':
            missiles.append(missile(ss.f.axis, ss.f.pos))
    

Next steps

Now that we’re cooking with charcoal — we have a functioning program — the next steps will be setting up a target to shoot at, and, if the students are ready, letting them program their ships to move autonomously. If they’re not ready then, as an example, I’ll program some opposition.

Practicing Plotting Points on the Co-ordinate Plain

Pre-Algebra class starts next week, so in preparation for one of the early lessons on how to plot x,y co-ordinates, I put together an interactive plotter that lets students drag points onto the co-ordinate plain.

Students practice plotting points by dragging the red dot to the coordinates given.

Usage

The program generates random coordinate pairs within the area of the chart (or you can enter values of the coordinates yourself):

  • Clicking the “Show Point” button will place a yellow dot at the point.
  • When you’re confident you understand how the coordinate pairs work, you can practice by dragging the red dot to where you think the point is and the program will tell you if you’re right or not.

About the Program

This interactive application uses the jQuery and KineticJS javascript libraries. The latter library in particular is useful for making the HTML5 canvases interactive, so you can click on points on the graph and drag them.

When I have some time, after classes settle down, I’ll see if I can figure out how to embed this type of app into this (WordPress) blog. KineticJS is based off HTML5 canvases, which is what I use for the other interactive graphs I’ve posted, so it shouldn’t be terribly hard (at least in principle).

Switch: Converting Audio Files for Sound Effects

Adding sound to webpages is pretty easy with HTML5’s Audio tag (as Jean-Baptiste Jung demonstrates), but since different browsers can only handle certain, different types of sound files, you’ll often need to convert your sound into .wav, .ogg, and .mp3 formats so they can work for everyone.

Switch Sound File Converter does just what it says, and is free for non-commercial use.

Learning Matricies by Programming a Matrix Solver

One of my pre-Calculus students convinced me that the best way for him to learn how to work with matrices was for him to program a matrix solver. I helped him create a Gaussian Elimination solver in Python (which we’ve been using since last year).

Gaussian Elimination Matrix Solver by Alex Shine (comments by me).

from visual import *

'''The coefficient matrix (m) can be any sized square matrix 
   with an additional column for the solution'''
m = [ [1,-2,1,-4],[0,1,2,4],[2,3,-2,2]]

'''Convert the input matrix, m (which is a list) into an array'''
a = array(m,float)
print "Input matrix:"
print a
print

'''Get the shape of the matrix (number of rows and columns)'''
(r,c) = a.shape 

rs = (1-r)

'''Solve'''
for j in range(r):

    print "Column #:", j
    for i in range(r):
        if a[i,j] <> 0:
            a[i,:] = a[i,:] / a[i,j]
    print a
    print

    for i in range(rs+j,j):
        if a[i,j] <> 0:
            a[i,:] = a[j,:] - a[i,:]
    print a
    print

print "Solution"
for i in range (r):
    a[i,:] = a[i,:] / a[i,i]
    print "Variable", i, "=", a[i,-1]

print
print "Solution Matrix:"
print a

The code above solves the following system of equations:


  x - 2y +  z = -4 
       y + 2z =  4 
 2x + 3y - 2z =  2 

Which can be written in matrix form as such:

 \left[ \begin{array}{ccc} 1 & -2 & 1 \\ 0 & 1 & 2 \\ 2 & 3 & -2 \end{array} \right]  \left[ \begin{array}{c} x \\ y \\ z \end{array} \right] =  \left[ \begin{array}{c} -4 \\ 4 \\ 2 \end{array} \right]

You use the solver by taking the square matrix on the left hand side of the equation and combining it with the column on the hand side as an additional column:

 \left[ \begin{array}{cccc} 1 & -2 & 1 & -4 \\ 0 & 1 & 2 & 4\\ 2 & 3 & -2 & 2\end{array} \right]

This is entered into the program as the line:

m = [ [1,-2,1,-4],[0,1,2,4],[2,3,-2,2]]

When you run the above program you should get the results:

>>> ================================ RESTART ================================
>>> 
Input matrix:
[[ 1. -2.  1. -4.]
 [ 0.  1.  2.  4.]
 [ 2.  3. -2.  2.]]

Column #: 0
[[ 1.  -2.   1.  -4. ]
 [ 0.   1.   2.   4. ]
 [ 1.   1.5 -1.   1. ]]

[[ 1.  -2.   1.  -4. ]
 [ 0.   1.   2.   4. ]
 [ 0.  -3.5  2.  -5. ]]

Column #: 1
[[-0.5         1.         -0.5         2.        ]
 [ 0.          1.          2.          4.        ]
 [-0.          1.         -0.57142857  1.42857143]]

[[ 0.5         0.          2.5         2.        ]
 [ 0.          1.          2.          4.        ]
 [ 0.          0.          2.57142857  2.57142857]]

Column #: 2
[[ 0.2  0.   1.   0.8]
 [ 0.   0.5  1.   2. ]
 [ 0.   0.   1.   1. ]]

[[-0.2  0.   0.   0.2]
 [ 0.  -0.5  0.  -1. ]
 [ 0.   0.   1.   1. ]]

Solution
Variable 0 = -1.0
Variable 1 = 2.0
Variable 2 = 1.0

Solution Matrix:
[[ 1. -0. -0. -1.]
 [-0.  1. -0.  2.]
 [ 0.  0.  1.  1.]]
>>>

Be aware that:

  • The code is designed to take any size of matrix.
  • The matrix you put in can not have any zeros on its diagonal, so some manipulation is often necessary before you can use the code.

Other notes:

  • The negative zeros (-0) that show up especially in the solution matrix may not look pretty but do not affect the solution.
  • The code imports the vpython module in the first line but what it really needs is the numpy module, which vpython imports, for the arrays.

The next step is to turn this into a function or a class that can be used in other codes, but it’s already proved useful. My calculus students compared their solutions for the coefficients of a quadratic equation that they had to solve for their carpet friction experiment, which was great because their first answers were wrong.

A calculus student uses the matrix solver. Mr. Shine is now trying to convert the solver into an iPhone app.