title: Simpson's rule in Python
body: <p>For a numerical methods class, I need to write a program to evaluate a definite integral with Simpson's composite rule. I already got this far (see below), but my answer is not correct. I am testing the program with f(x)=x, integrated over 0 to 1, for which the outcome should be 0.5. I get 0.78746... etc.
  I know there is a Simpson's rule available in Scipy, but I really need to write it myself.</p>
  
  <p>I suspect there is something wrong with the two loops. I tried "for i in range(1, n, 2)" and "for i in range(2, n-1, 2)" before, and this gave me a result of 0.41668333... etc.
  I also tried "x += h" and I tried "x += i*h". The first gave me 0.3954, and the second option 7.9218.</p>
  
  <pre><code># Write a program to evaluate a definite integral using Simpson's rule with
  # n subdivisions
  
  from math import *
  from pylab import *
  
  def simpson(f, a, b, n):
      h=(b-a)/n
      k=0.0
      x=a
      for i in range(1,n/2):
          x += 2*h
          k += 4*f(x)
      for i in range(2,(n/2)-1):
          x += 2*h
          k += 2*f(x)
      return (h/3)*(f(a)+f(b)+k)
  
  def function(x): return x
  
  print simpson(function, 0.0, 1.0, 100)
  </code></pre>