Simpson revival

[Albert Chan]

In Python, Romberg's method corrections can be simplified with generator.

Code:

from __future__ import division

def simpson(f, a, b, n=2):          # Simpson generator
    h = (b-a) / 2
    t0 = (f(a)+f(b)) * h
    while 1:
        s = sum(f(a + i*h) for i in xrange(1, n, 2))
        t1 = t0/2 + s*h             # simpson's partial sum
        yield t1 + (t1-t0)/3        # simpson's formula
        t0, n, h = t1, 2*n, h/2

Code:

def _romberg(gen, y, k):            # Romberg generator
    for z in gen:
        yield z + (z-y)/k
        y = z

Code:

def romberg(f, a, b, verbal=False, k=15):
    "Romberg's Method integration"
    gen = simpson(f, a, b, n=2)
    x = gen.next()
    while 1:
        y = gen.next()
        x = y + (y-x) / k
        if verbal: print x
        if x == y: return x         # stop if converged
        gen = _romberg(gen, y, k)   # Romberg corrections
        k = 4 * k + 3

Redo previous post example:

>>> romberg(lambda x: 1/x, 1, 2, verbal=True)
0.693174603175
0.693147477645
0.693147181917
0.693147180562
0.69314718056
0.69314718056
0.69314718055994506