0 Replies - 1140 Views - Last Post: 14 October 2012 - 08:28 AM Rate Topic: -----

#1 Nemeenem  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 1
  • Joined: 14-October 12

Numeric Integration(Simpson's Rule)

Posted 14 October 2012 - 08:28 AM

So I have some concern about my code and its precision of calculations,

Relevant information:
-Language: Python 3.2.3
-Simpson's rule:Posted Image
-Error:
Posted Image

Here is my code, where limits is a list of two floats with the limits of integration, dx is the step(b-a)/n as a float, intervals is the number of subinvtervals used in simpsons rule, and func is the parsed equation
def simpsonsApprox(dx, limits, intervals, func):
    
    x = limits[0]
    summation = float(eval(func))

    for i in range(1,intervals,2):
        x = limits[0] + (dx*i)
        summation += 4*float(eval(func))

    for i in range(2,intervals-1,2):
        x = limits[0] + (dx*i)
        summation += 2*float(eval(func))
        
    x = limits[1]
    summation += float(eval(func))
    summation = dx * summation / 3
    return summation

My concern is as follows:
According to the error formula above for the equation 'sqrt(1-x^2)-x' to reach precision of 14 digits(E_s = 10^-14) you must use around 737 subintervals. The integral of sqrt(1-x^2)-x from 0 to 1/sqrt(2) is exactly pi/8.

In my tests using ~737 intervals I can only get 12 digits of pi/8 when according to my calculations I should be getting 14.

I did research and found that a python float is the same as a C double, with precision of 53 bits which should be around 16 digits.

I'm not trying to get 16 digits or beyond so I'm confused as to why I'm not getting my anticipated results.

I found that using more subintervals, something around 1400, I could get 14 digits. But, this does not follow the formula above.

The questions I pose:
Is the rounding error at the calculation of each subinterval the reason for this result?
Is this due to a logical flaw in my code?(Although it seems to work fine)
Could the use of the eval() function be at fault?
I was told that a float should work fine for this, and I wouldn't need to use anything more precise e.g. Decimal

Any input on this matter would be greatly appreciated, I hope that I explained the problem well enough that it's not necessary to understand the whole calculus part.

Regards, NEM

Is This A Good Question/Topic? 0
  • +

Page 1 of 1