外挿法で円周率を求める

外挿法は，粗い計算結果から極限値を推定する方法です．この記事では，円周率近似の例を出発点にして，多角形近似，誤差展開，Richardson extrapolation の考え方を整理します．題材は円周率ですが，見ているものは有限要素法にも現れる「誤差の主項を消す」という考え方です．

数値計算では，ある量 \( a_0 \) を直接求めることが難しいため，パラメータ \( h>0\) に依存する近似値 \( a(h)\) を計算し，
\[ a_0 = \lim_{h\to 0} a(h) \] として目的の値を得ようとすることが多くあります．ただし，\( h\) を小さくするだけでは計算量が増えます．また，丸め誤差も無視できなくなります．外挿法の基本的な考え方は，いくつかの \( h\) に対する値 \[ a(h_0),\ a(h_1),\ a(h_2),\ldots \] を使い，\( h=0 \) での値を推定することです．重要なのは，単に点を結ぶことではありません．背後に \[ a(h)=a_0+a_1h^q+a_2h^{2q}+a_3h^{3q}+\cdots \] のような漸近展開があるとき，外挿は特に強力になります．つまり，外挿法は「誤差の構造を利用する方法」です．

半径1の円を考えます．円周の長さは \(2\pi\) です．この円に内接する正 \(n\) 角形を考えると，一辺の長さは \[ 2\sin\frac{\pi}{n} \] です．したがって，正 \(n\) 角形の周長の半分は

\[ T_n := n\sin\frac{\pi}{n} \] となります．これは円周の半分 \(\pi\) より小さいので， \[ T_n < \pi \] です．一方，外接正 \(n\) 角形を用いれば \[ \pi < U_n := n\tan\frac{\pi}{n} \] も得られます．したがって， \[ n\sin\frac{\pi}{n} < \pi < n\tan\frac{\pi}{n} \] という上下評価が得られます．

ここでは実行用の計算では \(\sin(\pi/n)\) を直接使いません．円周率を求めたいのに，計算の中で \( \pi\) を使うと，ちょっと手品のタネが見えてしまうからです．

単位円に内接する正 \(n\) 角形を考え，その一辺の長さを \(s_n\) とします．中心角は \(2\pi/n\) なので， \[ s_n=2\sin\frac{\pi}{n} \] です．特に，正六角形では \[ s_6=2\sin\frac{\pi}{6}=1 \] となります．

正 \(n\) 角形から正 \(2n\) 角形へ移るには， \[ s_{2n}=2\sin\frac{\pi}{2n} \] を求めればよいです．半角公式を使うと， \[ s_{2n} = \sqrt{2-\sqrt{4-s_n^2}} \] が得られます．しかし，この形では \(n\) が大きいときに桁落ちが起こりやすいです． そこで有理化して， \[ s_{2n} = \frac{s_n}{\sqrt{2+\sqrt{4-s_n^2}}} \] と計算します．これが数値的に安定な形です．

以上により，正六角形の一辺 \(s_6=1\) から出発して，平方根だけを用いる更新式により \(s_{12},s_{24},s_{48},\ldots\) を順に計算できます．そして \[ \pi_n:=\frac{n s_n}{2} \] とおけば，これは単位円の周長と直径の比，すなわち \(\pi\) の近似値になります．実際，理論的には \(s_n=2\sin(\pi/n)\) なので， \[ \pi_n=n\sin\frac{\pi}{n} \] です．これは内接正 \(n\) 角形による近似であり， \[ \pi_n<\pi,\quad \pi_n\to\pi\quad(n\to\infty) \] が成り立ちます．

重要なのは，実際の計算では \(\pi\) や \(\sin(\pi/n)\) を使っていない点です．\(\pi\) を含む式はあくまで理論的な説明であり，計算自体は正六角形の値 \(s_6=1\) と平方根の反復だけで行われます．したがって，ここで得られる \(\pi_n\) は，\(\pi\) を前もって知っているから出てくる値ではなく，多角形近似から作られる純粋な近似列です．

ここまでで，\(\pi\) に近づく計算可能な列 \[ \pi_6,\pi_{12},\pi_{24},\ldots \] が得られました．しかし，この列はそのままでは比較的ゆっくり収束します．外挿法の考え方は，この近似列の誤差構造を利用して，より高精度な近似値を作ることにあります．

理論を調べるために，ここでは一時的に \[ T_n = n\sin\frac{\pi}{n} \] という表示を使います．\(h=1/n\) とおくと， \[ T(h) = \frac{1}{h}\sin(\pi h) \] です．Taylor 展開より， \[ T(h) = \pi - \frac{\pi^3}{3!}h^2 + \frac{\pi^5}{5!}h^4 - \frac{\pi^7}{7!}h^6 + \cdots \] となります．つまり，\(T(h)\) は \[ T(h)=\pi+c_2h^2+c_4h^4+c_6h^6+\cdots \] という偶数冪の誤差展開を持っています．この「偶数冪で展開される」という事実が，外挿で大きく効きます．

いま， \[ T_n = \pi + c_2h^2 + c_4h^4 + O(h^6), \quad h=\frac{1}{n}, \] とします．辺数を倍にすると \(h\) は \(h/2\) になるので， \[ T_{2n} = \pi + c_2\left(\frac{h}{2}\right)^2 + c_4\left(\frac{h}{2}\right)^4 + O(h^6) \] です．ここで \[ S_n := \frac{4T_{2n}-T_n}{3} \] とおくと，\(c_2h^2\) の項が消えます．実際， \[ S_n = \pi + O(h^4) \] となります．もとの多角形近似は \(O(h^2)\) でしたが，この線形結合により \(O(h^4)\) の近似に改善されます．

外挿法の要点
\(h\) を小さくするだけでなく，\(h\) と \(h/2\) の計算結果を組み合わせて，誤差の主要項を消します．これが Richardson extrapolation の最も基本的な形です．

さらに高次の誤差項も順に消したい場合は，外挿表を作ります．まず， \[ R_{k,0}=T_{6\cdot 2^k} \] とおきます．そして \(j\geq 1\) に対して， \[ R_{k,j} =\frac{4^jR_{k,j-1}-R_{k-1,j-1}}{4^j-1} \] と定めます．ここで \(4^j\) が出てくるのは，誤差展開が \[ h^2,\ h^4,\ h^6,\ldots \] という偶数冪で進むからです．

以下のコードでは，正六角形から始めて，辺数を倍々に増やします．計算本体では \(\pi\) を使いません．最後の誤差表示だけ，確認用に math.pi を使っています．

import math


def inscribed_polygon_pi_sequence(levels=8):
    """内接正多角形による円周率近似を返す．"""
    n = 6
    s = 1.0  # 正六角形の一辺の長さ
    values = []

    for k in range(levels):
        pi_n = n * s / 2.0
        values.append((k, n, pi_n))

        # 半角公式による安定な更新式．
        # s_n = 2 sin(theta) とすると，
        # s_{2n} = 2 sin(theta/2)
        #        = s_n / sqrt(2 + 2 cos(theta))．
        cos_theta = math.sqrt(1.0 - (s / 2.0) ** 2)
        s = s / math.sqrt(2.0 + 2.0 * cos_theta)
        n *= 2

    return values


def richardson_table(values):
    """偶数冪の誤差展開を仮定した Richardson 外挿表を作る．"""
    table = []

    for k, (_, _, value) in enumerate(values):
        row = [value]

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

        table.append(row)

    return table


def main():
    levels = 8
    values = inscribed_polygon_pi_sequence(levels)
    table = richardson_table(values)

    print(f"{'k':>2} {'n':>7} {'polygon_pi':>20} {'raw_error':>12} "
          f"{'extrapolated_pi':>20} {'extrap_error':>12}")
    print("-" * 78)

    for k, n, pi_n in values:
        extrapolated = table[k][-1]
        raw_error = abs(math.pi - pi_n)
        extrap_error = abs(math.pi - extrapolated)

        print(f"{k:2d} {n:7d} {pi_n:20.15f} {raw_error:12.3e} "
              f"{extrapolated:20.15f} {extrap_error:12.3e}")


if __name__ == "__main__":
    main()

 k       n           polygon_pi    raw_error      extrapolated_pi extrap_error
------------------------------------------------------------------------------
 0       6    3.000000000000000    1.416e-01    3.000000000000000    1.416e-01
 1      12    3.105828541230249    3.576e-02    3.141104721640332    4.879e-04
 2      24    3.132628613281238    8.964e-03    3.141592453897650    1.997e-07
 3      48    3.139350203046867    2.242e-03    3.141592653577892    1.190e-11
 4      96    3.141031950890509    5.607e-04    3.141592653589793    0.000e+00
 5     192    3.141452472285462    1.402e-04    3.141592653589793    4.441e-16
 6     384    3.141557607911857    3.505e-05    3.141592653589793    0.000e+00
 7     768    3.141583892148318    8.761e-06    3.141592653589792    8.882e-16

正 \(96\) 角形による素朴な近似では，誤差はまだ約 \(5.6\times 10^{-4}\) です．しかし，外挿表の対角成分を使うと，同じ程度のデータから倍精度計算の限界に近い値まで到達します． ここで誤差が 0.000e+00 と表示される部分がありますが，数学的に誤差がゼロという意味ではありません．倍精度浮動小数点数で見える範囲では，math.pi との差が表示されなくなっているという意味です．

この例で重要なのは，円周率そのものではありません．重要なのは，近似値の誤差が \[ c_2h^2+c_4h^4+c_6h^6+\cdots \] という形を持つと分かれば，その構造を利用して精度を上げられるという点です．

これは有限要素法でも同じです．たとえば，ある点 \(P\) における Ritz 射影の誤差が \[ R_hu(P)=u(P)+h^2 e(P)+O(h^4) \] のような展開を持てば，\(h\) と \(h/2\) の近似解を組み合わせることで，\(h^2 e(P)\) の主誤差を消すことができます．ただし PDE の場合には，解の正則性，領域の角，境界近傍のメッシュ構造，安定性評価などが絡みます．そこが円周率の例との大きな違いです．

• R. Rannacher, Special Topics in Numerics III: Extrapolation and Defect Correction, Lecture Notes, Heidelberg University, SS 2017, Section 6.1.1 and 6.2.
• L. F. Richardson and J. A. Gaunt, The deferred approach to the limit, Philosophical Transactions of the Royal Society of London. Series A, 1927.