数値微分（2）

前回のグラフを見ると、両対数グラフで直線的に真値に近づいていく。だから、よりよい微分の近似値が得られる。Taylor展開して、

f(x+h) = f(x) + hf'(x) + h2/2!f''(x) + ...

f'(x) - (f(x+h) - f(x)) / h = h/2!f''(x) + ...

だから、近似と真値との差は、近似でhに比例する。しかし、f''(x)が0ならより高次になる。近似であることは無視して、

Ak ≡ (f(x+h/2k) - f(x)) / (h/2k)

として、

f'(x) - A0 = hnB

f'(x) - A1 = (h/2)nB

f'(x) - A2 = (h/4)nB

2n(f'(x) - A1) = f'(x) - A0

2n(f'(x) - A2) = f'(x) - A1

(f'(x) - A1)2 = (f'(x) - A0)(f'(x) - A2)

f'(x) = (A0 + A2 - 2A1) / (A12 - A0A2)

2n = (A0 - A1) / (A1 - A2)

これで、テキストと合っている。
コードにすると、

from math import *
def f(x):

    return sqrt(x * x + x + 1)
def diff(f, x, h):

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

    a0 = diff(f, x, h)

    a1 = diff(f, x, h * 0.5)

    a2 = diff(f, x, h * 0.25)

    return (a0 * a2 - a1 * a1) / (a0 + a2 - a1 * 2)
h = 1

while h > 1e-8:

    print "%e %.16f" % (h, diff2(f, 1, h))

    h /= 2.0

print sqrt(3) / 2

新しい数値微分を赤でプロットすると、

新しいほうは、2乗で収束していることがわかる。