hplgit/scipro-primer

Switch branches/tags
Nothing to show
Fetching contributors…
Cannot retrieve contributors at this time
62 lines (51 sloc) 1.56 KB
 def Simpson(f, a, b, n=500): """ Return the approximation of the integral of f from a to b using Simpson's rule with n intervals. """ if a > b: print 'Error: a=%g > b=%g' % (a, b) return None # Check that n is even if n % 2 != 0: print 'Error: n=%d is not an even integer!' % n n = n+1 # make n even h = (b - a)/float(n) sum1 = 0 for i in range(1, n/2 + 1): sum1 += f(a + (2*i-1)*h) sum2 = 0 for i in range(1, n/2): sum2 += f(a + 2*i*h) integral = (b-a)/(3*n)*(f(a) + f(b) + 4*sum1 + 2*sum2) return integral def h(x): return (3./2)*sin(x)**3 from math import sin, pi def application(): print 'Integral of 1.5*sin^3 from 0 to pi:' for n in 2, 6, 12, 100, 500: approx = Simpson(h, 0, pi, n) print 'n=%3d, approx=%18.15f, error=%9.2E' % \ (n, approx, 2-approx) def test_Simpson(): """Test exact integration of quadratic polynomials.""" a = 1.5 b = 2.0 n = 8 g = lambda x: 3*x**2 - 7*x + 2.5 # test integrand G = lambda x: x**3 - 3.5*x**2 + 2.5*x # integral of g exact = G(b) - G(a) approx = Simpson(g, a, b, n) success = abs(exact - approx) < 1E-14 msg = 'Simpson: %g, exact: %g' % (approx, exact) assert success, msg def experiments(): """Some other tests with known answers.""" from math import sin, pi, exp, log n = 200 print 'sin from 0 to pi:', Simpson(sin, 0, pi, n) print 'exp from 0 to log(3):', Simpson(exp, 0, log(3), n) test_Simpson() application() experiments()