-Language: Python 3.2.3
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 summation = float(eval(func)) for i in range(1,intervals,2): x = limits + (dx*i) summation += 4*float(eval(func)) for i in range(2,intervals-1,2): x = limits + (dx*i) summation += 2*float(eval(func)) x = limits 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.