数値微分(1)

なぜ関数プログラミングは重要か
http://www.sampou.org/haskell/article/whyfp.html
の 4.2. 数値微分 に面白いことが書いてある。
微分の近似値を計算するのに、f'(x) ≒ (f(x + h) - f(x)) / hを使うが、当然hが小さくなるほど真値に近くなる。しかし、あまり小さすぎると数値誤差が大きくなってしまう。どれくらいのhが最適なのだろう。
感覚的にはわかっていたが、そういえば試したことはない。そこで次のようなコードで試してみた。


from math import *

def f(x):
return sqrt(x * x + x + 1)

def diff(f, x, h):
return (f(x + h) - f(x)) / h

h = 1
while h > 1e-15:
print "%e %.16f" % (h, diff(f, 1, h))
h /= 2.0
print sqrt(3) / 2

関数は、√(x2+x+1)で微分すると、(2x+1)/2√(x2+x+1)、x=1で√3/2となる。
hと真値との差(の絶対値)をプロットした。

これは、関数によって大きく異なるのだろうが、グラフを見る限り、1e-7くらいまでは順調に真値に近づいている。1.49e-8で最も真値に近づいているが、これはたまたまだろう。