I am trying to program some simple math concepts in order to familiarize myself with the Python language, and one of the first things I tried was a derivative calculator, but I ran into a problem as values got smaller and smaller.
from math import *
def f(x):
return x**2;
def lim(x,h):
return ((f(x + h)-f(x))/h);
a = float(input("At what point? "))
b = 1
for i in range(20):
print(lim(a,b))
b = b/10
Of course, the derivative is defined as the limit as h approaches zero of the function denoted in my code as "lim(x,h)".
The loop in this code is less than ideal, but initially I tried to devise some ways to have h approach zero as close as possible, and I found that when I did that I started getting terrible results. As such, I made this code to try to see where the problem began.
I noticed that after h = 1e-7, the value that my code gives starts to increasingly deviate from the expected result of 2x, and furthermore it appears to deviate seemingly devoid of any pattern. For example, when I input 7 as my x value, the output increasingly approaches 14, as expected, up until h=1e-7. After that point it departs from the trend of approaching closer and closer to the expected value, and around h = 1e-14 there are some wild deviations that were, at least for me, incredibly unexpected.
For reference, this is the output of the code with input 7:
At what point? 7
15.0
14.099999999999966
14.009999999999678
14.00100000000748
14.000099999975646
14.000009999648453
14.00000100204579
14.00000016360536
13.99999973727972
14.00000115836519
13.999965631228402
13.99982352268125
13.997691894473972
14.068746168049982
13.500311979441902
14.210854715202002
0.0
0.0
0.0
0.0
I don't think there is anything wrong with this code, it should give the value of (f(x+h) - f(x))/h as h approaches zero. Is there some limitation on computer mathematics to deal with numbers when they get this small? I noticed the same pattern happens when I try inputting the numbers into Wolfram alpha's engine.
Let me give an example. For the value x = 7 and h = 1e-14, doing the math by hand I would expect a return of 14 + 1e-14.
However when I plug that problem into Wolfram alpha, I get the exact same unexpected result that I get from my code of 13.500311979441902 (pictures of math worked out by hand and wolfram attached).
What causes this problem? How would I go about creating a derivative calculator with very high degrees of accuracy?
Python float use a fixed size data, like IEEE 754 double precision. Thus about 264 different values are encodable.
Calculations like x + h, x**2, a-b, a/b often result, not in the exact math result, but a nearby float. This calculation error becomes more important for OP's f(x + h)-f(x) as |h| becomes much smaller than |x|. @jasonharper
Eventually f(x + h)-f(x) becomes 0.0, as OP has seen.
This can be somewhat mitigated by using a power of 2 as that will induce less error. (float is encoded using powers-of-2) Yet even that will eventually result in f(x + h)-f(x) after a few dozen iterations.
Once |f(x + h)-f(x)| < |f(x)|*about_1e-15, h is too small to certainly form a proper lim(x,h) return value - low 50s iterations.
As for me, I'd consider lim(x, x*pow(2, -52/2)) for quick and dirty derivative.
(Details omitted for brevity. Much depends of the complexity of f() and value of x)
# b = 1
b = a
# for i in range(20):
for i in range(20 * 3):
print(lim(a,b))
# b = b/10
b = b/2
Output
At what point?
7
15.0
14.5
14.25
14.125
14.0625
14.03125
14.015625
14.0078125
14.00390625
14.001953125
14.0009765625
14.00048828125
14.000244140625
14.0001220703125
14.00006103515625
14.000030517578125
14.000015258789062
14.000007629394531
14.000003814697266
14.000001907348633
14.000000953674316
14.000000476837158
14.000000238418579
14.00000011920929
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
16.0
16.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0