Python/maths/simpson_rule.py

'''
Numerical integration or quadrature for a smooth function f with known values at x_i

This method is the classical approch of suming 'Equally Spaced Abscissas' 

method 2: 
"Simpson Rule"

'''
from __future__ import print_function


def method_2(boundary, steps):
# "Simpson Rule"
# int(f) = delta_x/2 * (b-a)/3*(f1 + 4f2 + 2f_3 + ... + fn)
	h = (boundary[1] - boundary[0]) / steps
	a = boundary[0]
	b = boundary[1]
	x_i = makePoints(a,b,h)
	y = 0.0
	y += (h/3.0)*f(a)
	cnt = 2
	for i in x_i:
		y += (h/3)*(4-2*(cnt%2))*f(i)	
		cnt += 1
	y += (h/3.0)*f(b)
	return y

def makePoints(a,b,h):
	x = a + h
	while x < (b-h):
		yield x
		x = x + h

def f(x): #enter your function here
	y = (x-0)*(x-0)
	return y

def main():
	a = 0.0 #Lower bound of integration
	b = 1.0	#Upper bound of integration
	steps = 10.0		#define number of steps or resolution
	boundary = [a, b]	#define boundary of integration
	y = method_2(boundary, steps)
	print('y = {0}'.format(y))

if __name__ == '__main__':
        main()