simpson

gistfile1.py

import math

def eval_func(function, lo, hi, n):
    # Takes a function, and evaluates it in n points evenly distributed from lo to hi
    step = abs(lo-hi)/float(n)
    ret = []
    for i in xrange(n+1):
        ret.append(function(lo + i*step))
    return ret

def simpson(function, lo, hi, n):
    # Takes four arguments, the function to approximate,
    # upper and lower integration limit, and the n=2m number of points to evaluate at.
    # Returns an approximated floating point number
    function_val = eval_func(function, lo, hi, n) # List of function values
    s_0 = function_val[0] + function_val[n]
    s_1 = 0
    s_2 = 0
    for i in xrange(1, len(function_val)-1):
        if i % 2 == 0:
            s_2 += function_val[i]
        else:
            s_1 += function_val[i]
    h = (hi-lo)/n
    return (h/3)*(s_0 + 4*s_1 + 2*s_2)

def main():
    f = lambda x: math.exp(-(x**2))
    print simpson(f, -1.0, 1.0, 100)
main()