Computing wrong value of integral using simpsons rule

Question

I wrote the following piece of code for Simpson's rule integration to approximate sin. The equation is in this attachment. I wrote separate loops for the even and odd terms, as they are shown grouped in the attachment.

import math 

x1=0
x2=math.pi
N=6
delta_x=(x2-x1)/N
f=math.sin
sum=f(0)

#odd summation
for i in range (1,N+1):
    sum=sum+f(x1+(2*i+1)*delta_x)

sum=4*sum
print(sum)

#even summation

for i in range(2,N+1):
    sumeven=0
    sumeven=sumeven+f(x1+(2*i)*delta_x)
sumeven=2*sumeven

sumeven=sumeven+f(N)
print(sumeven)

integral=(delta_x/3)*(sum+sumeven) 
print(integral)   

But when I print the values, it gives me very small negative numbers.

Can anyone see what is wrong with my code?

Answer 1

First of all, learn some basic debugging. You've done only two simple print lines for this code, neither of which traces the more detailed problems. See this lovely debug blog for help. I added some instrumentation, cleaned up some style, and ran the code again, tracing the logic and data flow somewhat better. Diagnosis below.

sum_odd = f(0)

# odd summation
for i in range (1, N+1):
    x_val = x1 + (2*i+1)*delta_x
    sum_odd = sum_odd + f(x_val)
    print (x_val/math.pi, "pi", f(x_val), sum_odd)

sum_odd = 4*sum_odd

print(sum_odd)

# even summation
for i in range(2, N+1):
    sum_even = 0
    x_val = x1 + (2*i)*delta_x
    sum_even = sum_even + f(x_val)
    print (x_val/math.pi, "pi", f(x_val), sum_even)

Output:

0.5 pi 1.0 1.0
0.8333333333333333 pi 0.5000000000000003 1.5000000000000004
1.1666666666666665 pi -0.4999999999999997 1.0000000000000007
1.5 pi -1.0 6.661338147750939e-16
1.8333333333333333 pi -0.5000000000000004 -0.4999999999999998
2.1666666666666665 pi 0.4999999999999993 -4.996003610813204e-16
-1.9984014443252818e-15
0.6666666666666666 pi 0.8660254037844387 0.8660254037844387
1.0 pi 1.2246467991473532e-16 1.2246467991473532e-16
1.3333333333333333 pi -0.8660254037844384 -0.8660254037844384
1.6666666666666665 pi -0.866025403784439 -0.866025403784439
2.0 pi -2.4492935982947064e-16 -2.4492935982947064e-16
-0.27941549819892636
-0.04876720424671586

DIAGNOSIS

You have two immediate problems:

1. You initialize sum_even inside the loop, which discards prior computations for that half of the series. Move it before the loop.
2. You integrated over twice the specified range; this, alone, will get you a tiny value (ideally, this would be 0, the integral for a full cycle of the function).

Fix the initialization, fix the limits, and you should see good results, more like this:

0.16666666666666666 pi 0.49999999999999994 0.49999999999999994
0.5 pi 1.0 1.5
0.8333333333333333 pi 0.5000000000000003 2.0000000000000004
8.000000000000002
0.3333333333333333 pi 0.8660254037844386 0.8660254037844386
0.6666666666666666 pi 0.8660254037844387 1.7320508075688772
1.0 pi 1.2246467991473532e-16 1.7320508075688774
3.184686116938829
1.9520959854268212

Also, I strongly recommend that you work through tutorials on debugging and coding style; these will help you in future work. I know from experience. :-)