数値微分（5）

ここまでやると、あとはテキストに従って書くだけ。まず、次のようにimproveという関数を定義して、新たな近似を作る。

def order(s):

    a = s.next()

    b = s.next()

    c = s.next()

    return int(log((a - c) / (b - c) - 1, 2.0) + 0.5)
def elimerror(n, l):

    b = l.next()

    while True:

        a = b

        b = l.next()

        yield (b * (2 ** n) - a) / (2 ** n - 1)
def improve(s):

    return elimerror(order(s), s)
def diff(eps, f, x):

    return within(eps, improve(differentiate(1.0, f, x)))

と思ったらそうはいかなかった。どこがまずいかというと、elimerrorの引数でsが2回使われているところである。orderのためにsを3回使って、そのあとにelimerrorを使ったときは4つ目からになってしまう。前にもこんなことがあった。
そこで、先頭に戻れるようにするクラスを作る。

class inf_list:

    def __init__(self, g):

        self.estimated = [ ]

        self.pos = 0

        self.g = g

    

    def __iter__(self):

        return self

    

    def next(self):

        if self.pos < len(self.estimated):

            ret = self.estimated[self.pos]

        else:

            ret = self.g.next()

            self.estimated.append(ret)

        self.pos += 1

        return ret

    

    def rewind(self):

        self.pos = 0

一度評価したらestimatedに入れる。posを持っておき、これがestimatedの長さ以上なら評価する。

def counter():

    n = -1

    while True:

        n += 1

        print "counter :", n

        yield n
g = inf_list(counter())

for i in range(2):

    print g.next()

g.rewind()

for i in range(3):

    print g.next()

これを実行すると、

counter : 0

0

counter : 1

1

0

1

counter : 2

2

これでようやく正しい加速バージョンが書ける。

def order(s):

    a = s.next()

    b = s.next()

    c = s.next()

    s.rewind()      // add

    return int(log((a - c) / (b - c) - 1, 2.0) + 0.5)
def improve(s):

    g = inf_list(s)

    return elimerror(order(g), g)

これで収束が速くなる。
さらに収束を加速するために、improveを繰り返し適用し、それぞれの2番目の値を取り出す。

def second(l):

    l.next()

    return l.next()
def super(s):

    return Map(second, Repeat(improve, s))
def diff(eps, f, x):

    return within(eps, super(differentiate(1.0, f, x)))

いくつかにinf_listをつけなければならない。
これであってるかなあ。
ソースは下で、結果はこう。

0.869295562859

0.866208623819

0.866030200973

0.866025459008

0.866025404029

0.866025404029

0.866025403784

from math import *

from functools import partial
def Repeat(f, a):

    yield a

    while True:

        a = f(a)

        yield a
def Map(f, l):

    while True:

        yield f(l.next())
def halve(x):

    return x / 2
def easydiff(f, x, h):

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

    return inf_list(Map(partial(easydiff, f, x), Repeat(halve, h0)))
def order(s):

    a = s.next()

    b = s.next()

    c = s.next()

    s.rewind()

    return int(log((a - c) / (b - c) - 1, 2.0) + 0.5)
def elimerror(n, l):

    b = l.next()

    while True:

        a = b

        b = l.next()

        yield (b * (2 ** n) - a) / (2 ** n - 1)
class inf_list:

    def __init__(self, g):

        self.estimated = [ ]

        self.pos = 0

        self.g = g

    

    def __iter__(self):

        return self

    

    def next(self):

        if self.pos < len(self.estimated):

            ret = self.estimated[self.pos]

        else:

            ret = self.g.next()

            self.estimated.append(ret)

        self.pos += 1

        return ret

    

    def rewind(self):

        self.pos = 0
def improve(s):

    s.rewind()

    g = inf_list(s)

    return inf_list(elimerror(order(g), g))
def second(l):

    l.next()

    return l.next()
def super(s):

    return Map(second, Repeat(improve, s))
def within(eps, l):

    b = l.next()

    while True:

        a = b

        b = l.next()

        print b

        if abs(a - b) <= eps:

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

    return within(eps, super(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