So I was doing calculus homework tonight.  Numerical estimation of integrals, using left endpoint, right endpoint, midpoint and trapezoidal estimation.  I thought that it was really tedious to write all that out by hand, so I thought I'd have the computer do it for me.  Voila.  A python script that will numerically estimate integrals.  Although you could just go to Wolfram Alpha.

Edit: added simpson's rule

69 lines of beauty:

#Numerical estimation of integrals
from math import *

def main():
    f = lambda x: x**3
   
    print left(f)
    print right(f)
    print midpoint(f)
    print trapezoid(f)
    f = lambda x: abs(8-x)
    print midpoint(f, 4, 7.0, 13.0)
    print trapezoid(f, 4, 7.0, 13.0)
    print simpson(f, 4, 7.0, 13.0)
    f = lambda x: e**(-5*x**2)
    print trapezoid(f, 4, 0.0, 1.0)
    print midpoint(f, 4, 0.0, 1.0)

    print simpson(f, 4, 0.0, 1.0)

def left(function, n = 20, start = 0.0, end = 10.0):
    s = 0
    dist = end-start
    deltaX = dist/n
    for x in range(n):
        s += function(start + deltaX*x)
    return s*deltaX

def right(function, n = 20, start = 0.0, end = 10.0):
    s = 0
    dist = end-start
    deltaX = dist/n
    for x in range(n):
        s += function(start + deltaX*(x+1))
    return s*deltaX

def midpoint(function, n = 20, start = 0.0, end = 10.0):
    s = 0
    dist = end-start
    deltaX = dist/n
    for x in range(n):
        s += function(start + 1/2.0*(deltaX*x+deltaX*(x+1)))
    return s*deltaX

def trapezoid(function, n = 20, start = 0.0, end = 10.0):
    s = 0
    dist = end-start
    deltaX = dist/n
    for x in range(n+1):
        if x in (0, n):
            s += function(start + deltaX*(x))
        else:
            s += 2*function(start + deltaX*(x))
    return s*deltaX/2

def simpson(function, n = 20, start = 0.0, end = 10.0):
    s = 0
    dist = end-start
    deltaX = dist/n
    for x in range(n+1):
        if x in (0, n):
            s += function(start + deltaX*x)
        elif x%2 == 1:
            s += 4*function(start + deltaX*x)
        else:
            s += 2*function(start + deltaX*x)

    return s*deltaX/3

if __name__ == "__main__":
    main()

Sooo, I've been playing around today with Bézier curves, and trying to implement those in Python and Pygame.  I finally got a piece of code that works, with much help from wikipedia's article on Bézier curves.  It still only has very basic functionality, but I had fun while doing it.  The curve follows the cursor.

Here's a screenshot of the script running:

Here's my code:


import pygame, sys
from pygame.locals import *
##Globals
BLACK = (  0,   0,   0)
WHITE = (255, 255, 255)
BLUE =  (  0,   0, 255)
GREEN = (  0, 255,   0)
RED =   (255,   0,   0)

FPS = 30
lx = 0
ly = 1


def main():
    global FPSCLOCK, SCREEN, circleList, mousex, mousey, screen_dims
    screen_dims = (900, 600)
    FPSCLOCK = pygame.time.Clock()
    pygame.display.init()
    SCREEN = pygame.display.set_mode(screen_dims)
    pygame.display.set_caption("BezierCurves")
    mousex, mousey = (0,0)
    SCREEN.fill(BLACK)
    drawBezier([(0,0), (100,200),(200,0)])
    x = 0
    while True:
        SCREEN.fill(BLACK)
        #print mousex, mousey
        drawBezier([(300,400), (mousex,mousey),(400, 100)])
        draw()
        checkForQuit()
        updateGameState()        
        pygame.display.update()
        FPSCLOCK.tick(FPS)

def drawBezier(pointList, color = WHITE):
    global SCREEN
    if len(pointList) < 2:
        raise ValueError

    d = listDistance(pointList)
    if d == 0:
        d=1
    #print d

    curve = pygame.PixelArray(SCREEN)
    for t in range(int(d + 1)):
        p = bezierPoint(pointList, t/d)
        #print p
        if inside((0,0), screen_dims, p):
            curve[int(p[lx])][int(p[ly])] = color
    del curve




def bezierPoint(pointList, t):
    """
    Recursive function to draw high order bezier curves.
    """
    if len(pointList) == 2:
        return tBetween(pointList[0], pointList[1], t )
    else:
        qPoints = [tBetween(pointList[p], pointList[p+1], t)
                   for p in range(len(pointList)-1)]
        return bezierPoint(qPoints, t)
    

def draw():
    pass

def tBetween(p1, p2, t):
    """
    Returns a point at a certain percentage t on a line between p1 and p2
    Acceptable values for t can be from 0 to 1.
    """
    if t > 1 or t < 0:
        raise ValueError
    return (p1[lx]+(p2[lx] - p1[lx])*t, p1[ly]+(p2[ly]-p1[ly])*t)

def listDistance(pointList):
    """
    this returns the distance of straight lines drawn
    between all points in pointList. For the purpose of
    estimating how many pixels need to be drawn to form
    a bezier curve without spaces
    """
    s = sum([distance(pointList[p], pointList[p+1])
                   for p in range(len(pointList)-1)])
    return s

def distance(p1, p2):
    dx = p1[lx] - p2[lx]
    dy = p1[ly] - p2[ly]
    return abs(dx**2 + dy**2)**.5

def inside(dimsMin, dimsMax, point):
    if point[lx] < dimsMin[lx] or point[lx] > dimsMax[lx]:
        return False
    if point[ly] < dimsMin[ly] or point[ly] > dimsMax[ly]:
        return False
    return True

def updateGameState():
    global mousex, mousey
    for event in pygame.event.get(MOUSEMOTION):
        mousex, mousey = event.pos
        #print mousex, mousey

def terminate():
    pygame.quit()
    sys.exit()

def checkForQuit():
    for event in pygame.event.get(QUIT):
        terminate()
    for event in pygame.event.get(KEYUP):
        if event.key == K_ESCAPE:
            terminate()
        pygame.event.post(event)


    

if __name__ == "__main__":
    main()