gistfile1.py

#!/usr/bin/python
# -*- coding: utf-8 -*-

'''
・ 台形公式(Trapezoidal rule)
・ シンプソンの公式(Simpson's rule)
'''


import math


def simpsons(src, dst, func=None, split=50):
    if func == None:
        func = lambda x: x

    area = 0
    if src == dst:
        return area
    elif src > dst:
        src, dst = dst, src

    h = float(dst - src) / (split * 2)
    xo, xe = 0, 0
    for n in range(0, split - 1):
        xo += func(src + h * (2 * n + 1))
        xe += func(src + h * (2 * n + 2))
    xo += func(dst - h)
    area = (func(src) + func(dst) + 4 * xo + 2 * xe) * h / 3
    return area


def trapezoidal(src, dst, func=None, split=50):
    if func == None:
        func = lambda x: x

    area = 0
    if src == dst:
        return area
    elif src > dst:
        src, dst = dst, src

    h = float(dst - src) / split
    for n in xrange(split):
        '''
            台形の公式: (上底 + 下底) * 高さ / 2
        '''
        upper = func(src + (h * n))
        lowwer = func(src + (h * (n + 1)))
        area += ((upper + lowwer) * h / 2)
    return area


def main():
    fx = lambda x: math.sqrt(abs(x))

    print 'Trapezoidal rule'
    print 'x(0->10): %f' % trapezoidal(0, 10)
    print 'x(0->10): %f' % trapezoidal(0, 10, func=lambda x: -x)
    for split in [1, 10, 50, 100] + range(500, 2000, 500):
        print '√|x|(0->20|%4d): %f' % (split,
            trapezoidal(0, 20, func=fx, split=split))
    print "Simpson's rule"
    print 'x(0->10): %f' % simpsons(0, 10)
    print 'x(0->10): %f' % simpsons(0, 10, func=lambda x: -x)
    for split in [1, 10, 50, 100] + range(500, 2000, 500):
        print '√|x|(0->20|%4d): %f' % (split,
            simpsons(0, 20, func=fx, split=split))


if __name__ == '__main__':
    main()