数値微分（3）

前回までは、数学的な補足をして、それをPythonで実装したが、従来の書き方、Cとなんら変わらない書き方になっている。なぜ関数プログラミングは重要かは、関数型なら徹底的に部品化ができて、数値微分をその例に出している。これをPythonで書けるかどうか試してみたい。

まず、

[ h, h/2, h/4, ... ]

という無限列を作る。無限列を作るにはジェネレータを使う。

def gen_half(x):

    yield x

    while True:

        x /= 2

        yield x

さらに、1/2にする演算を関数にすると、次のようになる。

def Repeat(f, a):

    yield a

    while True:

        a = f(a)

        yield a
def halve(x):

    return x / 2
h_list = Repeat(halve, h0)

Repeatは、aとfから

[ a, f(a), f(f(a)), f(f(f(a))), ... ]

という無限列を作る。これで、当初の目的の無限列を生成できる。
次に、

[ easydiff(f, x, h), easydiff(f, x, h/2), easydiff(f, x, h/4), ... ]

という無限列を作る。ここではたぶんmapを使う意味はない。

def easydiff(f, x, h):

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

    h_list = Repeat(halve, h0)

    while True:

        yield easydiff(f, x, h_list.next())

微分の近似の列が得られたので、最後に収束判定を部品化する。無限列の前2つの差がeps以下なら2番目を返す、そうでなければ無限列の最初をポップした無限列に同じ操作をする。
ただ、このような操作はPythonにはたぶんできない。しかたがないので、次のようなコードでごかました。

def within(eps, l):

    b = l.next()

    while True:

        a = b

        b = l.next()

        if abs(a - b) <= eps:

            return b

ともかく、これで部品化できた。

from math import *
def Repeat(f, a):

    yield a

    while True:

        a = f(a)

        yield a
def halve(x):

    return x / 2
def gen_half(x):

    yield x

    while True:

        x /= 2

        yield x
def easydiff(f, x, h):

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

    h_list = Repeat(halve, h0)

    while True:

        yield easydiff(f, x, h_list.next())
def within(eps, l):

    b = l.next()

    while True:

        a = b

        b = l.next()

        if abs(a - b) <= eps:

            return b
def diff(eps, f, x):

    return within(eps, differentiate(1.0, f, x))
def f(x):

    return sqrt(x * x + x + 1)
x = 1.0

eps = 1e-7

print diff(eps, f, x)

print sqrt(3) / 2       # reference