0

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?

Prune
  • 76,765
  • 14
  • 60
  • 81
  • I haven't really dug into it, but the first thing I'd check is whether you're getting integer division somewhere. At least if you're on Python 2. – Peter DeGlopper Dec 06 '17 at 18:38
  • Just glancing at your code, you use `sumeven=0` *inside* your loop, i.e. reseting the `sumeven` value to `0` on each iteration, making you get only the last value, not the sum of the sequence. – juanpa.arrivillaga Dec 06 '17 at 18:41

1 Answers1

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. :-)

Prune
  • 76,765
  • 14
  • 60
  • 81
  • I'm sorry, I fail to see how I integrated over twice the specified range. –  Dec 06 '17 at 20:41
  • @Soffie: That's one reason I put in the extra `print` statements, to show the actual values you're using. Look at how you calculate the `x` value in each loop: your `2*i` term drives it to twice the specified upper limit. – Prune Dec 06 '17 at 21:36
  • ah ok- but I need 2*i+1 and the 2*i for odd and even respectively so I can sum up the delta_x's multiplied by(1,3,5)...and (2,4,6..) for the odd and even terms. Do I need to go about it a different way? –  Dec 06 '17 at 21:53
  • Yes. Check your algebra, and compare to the `x` values you actually need. – Prune Dec 06 '17 at 22:54