Romberg積分入門

第1回では円周率の多角形近似，第2回では数値微分を扱いました．どちらでも重要だったのは，近似値そのものではなく，その背後にある誤差展開でした．第3回では，定積分を題材にします．台形公式の誤差展開を利用して，低次の誤差項を順に消す方法が Romberg 積分です．

今回の目的は，次の定積分 \[ I=\int_a^b f(x)\,dx \] を高精度に近似することです．もっとも基本的な数値積分公式の一つに，台形公式があります．台形公式は単純で実装しやすい一方で，そのままでは2次精度です．しかし，台形公式の誤差には \[ h^2,\ h^4,\ h^6,\ldots \] という偶数冪の展開が現れます．そのため，Richardson extrapolation と非常に相性がよいです．重要なのは，台形公式を単に細かくするだけでなく，Euler--Maclaurin 展開に基づいて \(h=0\) へ外挿する点です．

区間 \([a,b]\) を \(N\) 等分します． \[ h=\frac{b-a}{N},\quad x_j=a+jh \quad (j=0,1,\ldots,N) \] とおきます．このとき，合成台形公式は \[ T(h) = h\left\{ \frac12 f(a) +\sum_{j=1}^{N-1} f(x_j) +\frac12 f(b) \right\} \] で定義されます．これは，各小区間 \([x_j,x_{j+1}]\) 上でグラフを直線で近似し，その下の台形の面積を足し合わせる方法です．したがって，考え方そのものは非常に素朴です．しかし，この素朴な公式が外挿法と組み合わさることで，高精度な方法になります．

外挿法を使うには，近似値 \(T(h)\) がどのような誤差展開をもつかを知る必要があります．ここで現れるのが Euler--Maclaurin 公式です．\(f\) が十分滑らかならば，台形公式には次の形の展開があります． \[ T(h) = \int_a^b f(x)\,dx + \sum_{k=1}^m h^{2k} \frac{B_{2k}}{(2k)!} \left\{ f^{(2k-1)}(b)-f^{(2k-1)}(a) \right\} + O(h^{2m+2}). \] ここで \(B_{2k}\) は Bernoulli 数です．特に最初の項だけを見ると， \[ T(h) = I + \frac{h^2}{12} \left\{ f'(b)-f'(a) \right\} +O(h^4) \] です． したがって，台形公式の誤差は基本的に \(O(h^2)\) です．

誤差展開の最初の部分だけを取り出して， \[ T(h)=I+c_2h^2+c_4h^4+O(h^6) \] と書きます．同じ公式を \(h/2\) で計算すると， \[ T(h/2) = I+c_2\frac{h^2}4+c_4\frac{h^4}{16}+O(h^6) \] です．この二つを \[ R_1(h):=\frac{4T(h/2)-T(h)}3 \] と組み合わせると，\(c_2h^2\) の項が消えます．実際， \[ R_1(h)=I+O(h^4) \] です．つまり，台形公式を二つのメッシュ幅で計算するだけで，2次精度の近似から4次精度の近似を作れます．これは第1回の円周率近似や，第2回の中心差分で見た外挿と同じ考え方です．

1段の外挿で \(O(h^2)\) の誤差を消せるなら，さらに次の \(O(h^4)\) の誤差も消したくなります．そのために，台形公式を \[ h,\quad \frac h2,\quad \frac h4,\quad \frac h8,\ldots \] で順に計算し，外挿を表の形に整理します．これが Romberg 表です．

\(R_{i,0}\) を台形公式の値とします．たとえば \(h_i=(b-a)/2^i\) として， \[ R_{i,0}=T(h_i) \] です．第1列目は，単なる台形公式の値です．次に，外挿した値を \[ R_{i,k} \] と書きます．更新式は \[ R_{i,k} = R_{i,k-1} + \frac{R_{i,k-1}-R_{i-1,k-1}}{4^k-1}, \quad k=1,\ldots,i. \] です．分母に \(4^k-1\) が出る理由は，台形公式の誤差展開が \(h^2\) の冪で進むからです．\(k\) 段目の外挿では，主に \(h^{2k}\) の項を消しています．\(h\) を半分にすると \[ h^{2k} \quad\text{は}\quad \left(\frac h2\right)^{2k} = \frac{h^{2k}}{4^k} \] になります．そのため，補正係数に \(4^k-1\) が現れます．Romberg 表は，次のような三角形の表になります．

一般には，右下の対角成分 \[ R_{0,0},\ R_{1,1},\ R_{2,2},\ldots \] が高精度な近似値になります．

ここでは \[ I=\int_0^1 e^x\,dx \] を Romberg 積分で計算します．この積分は厳密に \[ I=e-1 \] と分かるので，誤差を確認しやすい例です．次のコードは，合成台形公式を計算し，その値から Romberg 表を作ります．HTMLの中にコードをそのまま掲載しているので，この部分だけをコピーしてPython ファイルとして実行することもできます．

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
Romberg integration.

This code accompanies the article:
外挿法と誤差展開（第3回）：台形公式から Romberg 積分へ

The goal is to show how the composite trapezoidal rule can be improved
by Richardson extrapolation.

We use the example

    I = ∫_0^1 exp(x) dx = e - 1.

The composite trapezoidal rule has an even-power error expansion:

    T(h) = I + c_2 h^2 + c_4 h^4 + c_6 h^6 + ...

Romberg integration removes these error terms successively.
"""

import math


def composite_trapezoidal(f, a, b, n):
    """
    Composite trapezoidal rule.

    Let h=(b-a)/n and x_j=a+jh. Then

        T_n = h { 1/2 f(a) + sum_{j=1}^{n-1} f(x_j) + 1/2 f(b) }.
    """
    h = (b - a) / n

    total = 0.5 * f(a) + 0.5 * f(b)

    for j in range(1, n):
        x_j = a + j * h
        total += f(x_j)

    return h * total


def romberg_table(f, a, b, max_level):
    """
    Build a Romberg table.

    We use n=1,2,4,8,... subintervals.
    Thus h is halved at each new row.

    The first column is the composite trapezoidal rule:
        R[i][0] = T(h_i).

    The remaining columns are computed by Richardson extrapolation:
        R[i][k] = R[i][k-1]
                  + (R[i][k-1] - R[i-1][k-1]) / (4^k - 1).

    The factor 4^k appears because the trapezoidal rule has an expansion
    in powers of h^2.
    """
    table = []

    for i in range(max_level + 1):
        n = 2 ** i

        row = [composite_trapezoidal(f, a, b, n)]

        for k in range(1, i + 1):
            correction = (row[k - 1] - table[i - 1][k - 1]) / (4 ** k - 1)
            row.append(row[k - 1] + correction)

        table.append(row)

    return table


def print_romberg_table(table, exact):
    """
    Print the Romberg table and the errors of the diagonal entries.
    """
    print("=" * 88)
    print("Romberg integration for I = integral_0^1 exp(x) dx")
    print("Exact value: e - 1")
    print("=" * 88)

    max_level = len(table) - 1

    header = " i   N       T(h)=R[i,0]"
    for k in range(1, max_level + 1):
        header += f"        R[i,{k}]"
    print(header)
    print("-" * 88)

    for i, row in enumerate(table):
        n = 2 ** i
        line = f"{i:2d} {n:4d} "
        for value in row:
            line += f" {value:15.12f}"
        print(line)

    print()
    print("Errors of the diagonal entries R[i,i]")
    print("-" * 88)
    print(f"{'i':>2} {'N':>4} {'R[i,i]':>18} {'error':>14} {'rate':>8}")

    previous_error = None
    previous_h = None

    for i, row in enumerate(table):
        n = 2 ** i
        h = 1.0 / n
        diagonal = row[-1]
        error = abs(diagonal - exact)

        if previous_error is None:
            rate_text = "   ---"
        else:
            rate = math.log(previous_error / error) / math.log(previous_h / h)
            rate_text = f"{rate:8.3f}"

        print(f"{i:2d} {n:4d} {diagonal:18.14f} {error:14.6e} {rate_text}")

        previous_error = error
        previous_h = h


def main():
    f = math.exp
    a = 0.0
    b = 1.0
    exact = math.e - 1.0

    table = romberg_table(f, a, b, max_level=5)

    print_romberg_table(table, exact)


if __name__ == "__main__":
    main()

実行結果の例

========================================================================================
Romberg integration for I = integral_0^1 exp(x) dx
Exact value: e - 1
========================================================================================
 i   N       T(h)=R[i,0]        R[i,1]        R[i,2]        R[i,3]        R[i,4]        R[i,5]
----------------------------------------------------------------------------------------
 0    1   1.859140914230
 1    2   1.753931092465  1.718861151877
 2    4   1.727221904558  1.718318841922  1.718282687925
 3    8   1.720518592164  1.718284154700  1.718281842218  1.718281828795
 4   16   1.718841128580  1.718281974052  1.718281828675  1.718281828460  1.718281828459
 5   32   1.718421660316  1.718281837562  1.718281828462  1.718281828459  1.718281828459  1.718281828459

Errors of the diagonal entries R[i,i]
----------------------------------------------------------------------------------------
 i    N             R[i,i]          error     rate
 0    1   1.85914091422952   1.408591e-01    ---
 1    2   1.71886115187659   5.793234e-04    7.926
 2    4   1.71828268792476   8.594657e-07    9.397
 3    8   1.71828182879453   3.354852e-10   11.323
 4   16   1.71828182845908   3.375078e-14   13.279
 5   32   1.71828182845904   6.661338e-16    5.663

最初の列 \(R[i,0]\) は，単なる台形公式です．行が下がるにつれて分割数 \(N\) が増え，\(h\) が半分になっています．一方，右側の列は外挿によって作られた値です．右下の対角成分 \(R[i,i]\) を見ると，急速に厳密値 \(e-1\) に近づくことが分かります．

Romberg 積分で起きていることは，次の一言で説明できます．

したがって，Romberg 積分は単なる「台形公式の改良版」ではありません． むしろ，

• 低次の近似公式を使う，
• 誤差展開を調べる，
• 主要誤差を線形結合で消す，
• 高次精度の近似を作る，

という外挿法の典型例です．

Romberg 積分がうまく働くためには，\(f\) が十分滑らかで，Euler--Maclaurin 型の誤差展開が有効であることが重要です．不連続な関数，端点で特異性をもつ関数，急激に振動する関数では，表の右下に進めば必ず高精度になるとは限りません．これは有限要素法でも同じです．外挿法は魔法ではなく，誤差展開があって初めて強力になります．

• 合成台形公式は，定積分を計算する基本的な公式です．
• 台形公式の誤差は，Euler--Maclaurin 展開により偶数冪で表されます．
• そのため，台形公式には \(h^2\) に関する Richardson extrapolation が自然に使えます．
• Romberg 積分は，台形公式と外挿法を組み合わせた高精度積分法です．
• ただし，十分な滑らかさと誤差展開の存在が重要です．

第4回では，この考え方を常微分方程式に移します．非stiffな問題に対する Gragg 型の外挿法が扱われています．数値積分より少し難しくなりますが，「低次の方法を基礎にして，誤差展開を利用して精度を上げる」という主役は同じです．

• Rolf Rannacher, Special Topics in Numerics III: Extrapolation and Defect Correction, Lecture Notes, Heidelberg University, 2017.