ロンバーグ積分

ロンバーグ積分とは，まず合成台形公式で近似し、次に Richardson の補外法で近似をより正確にする。

Romberg's method

に，C で書かれたプログラムがある。

これを，Python に書き換えると以下のようになる。

for をベクトル計算にするなどにより，ずいぶん短くなる？

import scipy as sp

def Romberg(f, a, b, max_steps=1000, acc=1e-14):
    Rp = sp.empty(max_steps)  # previous row
    Rc = sp.empty(max_steps)  # current row
    h = b - a
    Rp[0] = (f(a) + f(b)) * h / 2
    for i in range(1, max_steps):
        h /= 2
        Rc[0] = h * sum(f(sp.arange(1, 2 ** i, 2) * h)) + Rp[0] / 2  # R(i,0)
        for j in range(1, i + 1):
            n_k = 4 ** j
            Rc[j] = (n_k * Rc[j - 1] - Rp[j - 1]) / (n_k - 1)  # R(i,j)
        if i > 1 and abs(Rp[i - 1] - Rc[i]) < acc:
            return Rc[i - 1]
        Rp, Rc = Rc, Rp # swap Rn and Rc
    return Rp[max_steps - 1]  # best guess

def Pi(x):
    return 4 / (1 + x**2)

a = Romberg(Pi, 0, 1, max_steps=9, acc=1e-14)
a　# 3.141592653589794
sp.pi - a  # -8.881784197001252e-16