Geraldine - 1 year ago 167
Python Question

# Simpson's rule in Python

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.

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.

``````# 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)
``````

You probably forget to initialize `x` before the second loop, also, starting conditions and number of iterations are off. Here is the correct way:

``````def simpson(f, a, b, n):
h=(b-a)/n
k=0.0
x=a + h
for i in range(1,n/2 + 1):
k += 4*f(x)
x += 2*h

x = a + 2*h
for i in range(1,n/2):
k += 2*f(x)
x += 2*h
return (h/3)*(f(a)+f(b)+k)
``````

Your mistakes are connected with the notion of a loop invariant. Not to get into details too much, it's generally easier to understand and debug cycles which advance at the end of a cycle, not at the beginning, here I moved the `x += 2 * h` line to the end, which made it easy to verify where the summation starts. In your implementation it would be necessary to assign a weird `x = a - h` for the first loop only to add `2 * h` to it as the first line in the loop.

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download