バーゼル問題の近似法

世界は2乗でできている自然にひそむ平方数の不思議 (ブルーバックス)

作者: 小島寛之
出版社/メーカー: 講談社
発売日: 2013/08/20
メディア: 新書
この商品を含むブログ (13件) を見る

この本に、オイラーのバーゼル問題に対する執念が書かれていました。バーゼル問題とは、

S = 1 + 1/2² + 1/3² + …

を求めよ、というもので、π²/6となります。これを求めるために、オイラーは小数点以下20桁まで求めたというのです。

これを何の気なしに読むと、二乗の逆数を延々と計算し続けたのかと思ってしまうかもしれませんが、そんなことはありません。なぜなら、途中の第n項で打ち切ると、真値に対する差は、

$\frac{1}{(n+1)^2} + \frac{1}{(n+2)^2} + \dots \approx \int_{n+1/2}^{\infty}{\frac{1}{x^2}dx}$

と近似できて、これはおよそ1/nです。すなわち、10²⁰程度まで計算を行わなければ小数点以下20桁までは求められないということになります。

ではEulerはどのように20桁まで求めたのでしょうか。それには上の近似を使うことが考えられます。

$S \approx 1 + \frac{1}{2^2} + \frac{1}{3^2} + \dots \frac{1}{n^2} + \frac{1}{n+1/2}$

これももちろん誤差があります。これを評価しましょう。誤差は、

$\sum_{k=n+1}^{\infty}{\int_{k-1/2}^{k+1/2}{(\frac{1}{x^2}}-\frac{1}{k^2})dx}$
$=\sum_{k=n+1}^{\infty}{(\frac{1}{k-1/2}-\frac{1}{k+1/2}-\frac{1}{k^2})}$
$\approx \sum_{k=n+1}^{\infty}{\frac{1}{4k^4}} \approx \frac{1}{12n^3}$

2000000くらいまで計算すれば20桁の精度が得られるようです。まだまだですね。
このように次々と近似をよくすることができます。この近似を一般化しましょう。

$S(f, n) \equiv \sum_{k=n+1}^{\infty}{f(k)}$

を近似します。

$S = \int_{n+1/2}^{\infty}{f(x)dx} + \sum_{k=n+1}^{\infty}{(f(k) - \int_{k-1/2}^{k+1/2}{f(x)dx})}$

また、

$g(k) \equiv f(k) - \int_{k-1/2}^{k+1/2}{f(x)dx}$
$I(f, n) \equiv \int_{n+1/2}^{\infty}{f(x)dx}$

と置くと、

$S(f, n) = I(f, n) + S(g, n)$

と近似を進めることができます。本当は一般的には近似が進むわけではないのですが、その辺には目をつむります。
これを繰り返し適用すると、

$f(x) = f_0(x)$
$f_{i+1}(x) = f_i(x) - \int_{x-1/2}^{x+1/2}{f_i(t)dt}$
$S(f, n) = I(f_0, n) + \dots + I(f_m, n) + S(f_{m+1}, n)$

となります。 $f_0(x) = x^{-2}$ でもう一度計算すると、

$f_1(x) = \frac{1}{x^2} - \int_{x-1/2}^{x+1/2}{\frac{1}{t^2}dt}$
$= \frac{1}{x^2} - (\frac{1}{x-1/2} - \frac{1}{x+1/2})$

ただ、f₂を求めるときにf₁を積分するとlogが出てきてややこしくなるので、Taylor展開をすることにします。

$= \frac{1}{x^2} - (\frac{1}{x-1/2} - \frac{1}{x+1/2}) = \frac{1}{x^2} - (\frac{1}{x} + \frac{1}{2x^2} + \frac{1}{4x^3} + \frac{1}{8x^4} + \dots) + (\frac{1}{x} - \frac{1}{2x^2} + \frac{1}{4x^3} - \frac{1}{8x^4} + \dots)$
$= \frac{1}{4x^4} + \dots$

この計算を次々に行うのはオイラーならやるかもしれませんが、我々には計算機代数という武器があるのでそれを使いましょう。

$\frac{a_1}{x} + \frac{a_2}{x^2} + \dots$

を、Pythonで[ 0, a₁, a₂, … ]と表します。Taylor展開は適当な項で打ち切ります。結果は、

$f_{0}(x) = \frac{1}{x^{2}}$
$f_{1}(x) = -\frac{1}{4x^{4}}-\frac{1}{16x^{6}}-\frac{1}{64x^{8}}-\frac{1}{256x^{10}}-\frac{1}{1024x^{12}}-\frac{1}{4096x^{14}}-\frac{1}{16384x^{16}}-\frac{1}{65536x^{18}}-\frac{1}{262144x^{20}}-\frac{1}{1048576x^{22}}$
$f_{2}(x) = \frac{5}{24x^{6}}+\frac{7}{32x^{8}}+\frac{123}{640x^{10}}+\frac{253}{1536x^{12}}+\frac{2041}{14336x^{14}}+\frac{1023}{8192x^{16}}+\frac{32759}{294912x^{18}}+\frac{65531}{655360x^{20}}+\frac{524277}{5767168x^{22}}$
$f_{3}(x) = -\frac{35}{96x^{8}}-\frac{63}{64x^{10}}-\frac{1023}{512x^{12}}-\frac{9581}{2560x^{14}}-\frac{56277}{8192x^{16}}-\frac{206941}{16384x^{18}}-\frac{107585809}{4587520x^{20}}-\frac{5643}{128x^{22}}$
$f_{4}(x) = \frac{35}{32x^{10}}+\frac{385}{64x^{12}}+\frac{11869}{512x^{14}}+\frac{121121}{1536x^{16}}+\frac{31384873}{122880x^{18}}+\frac{66637807}{81920x^{20}}+\frac{338972673}{131072x^{22}}$
$f_{5}(x) = -\frac{1925}{384x^{12}}-\frac{25025}{512x^{14}}-\frac{325325}{1024x^{16}}-\frac{64919855}{36864x^{18}}-\frac{297217817}{32768x^{20}}-\frac{29580588387}{655360x^{22}}$
$f_{6}(x) = \frac{25025}{768x^{14}}+\frac{525525}{1024x^{16}}+\frac{10635625}{2048x^{18}}+\frac{357271915}{8192x^{20}}+\frac{21852078885}{65536x^{22}}$
$f_{7}(x) = -\frac{875875}{3072x^{16}}-\frac{20845825}{3072x^{18}}-\frac{826090265}{8192x^{20}}-\frac{59806671925}{49152x^{22}}$
$f_{8}(x) = \frac{14889875}{4608x^{18}}+\frac{56581525}{512x^{20}}+\frac{9449114675}{4096x^{22}}$
$f_{9}(x) = -\frac{282907625}{6144x^{20}}-\frac{17823180375}{8192x^{22}}$
$f_{10}(x) = \frac{9901766875}{12288x^{22}}$

これを使ってSの近似式を書くと、

$S \approx \sum_{k=1}^n{\frac{1}{k^2}+\frac{1}{(n+1/2)}-\frac{1}{12(n+1/2)^{3}}+\frac{7}{240(n+1/2)^{5}}-\frac{31}{1344(n+1/2)^{7}}+\frac{127}{3840(n+1/2)^{9}}-\frac{2555}{33792(n+1/2)^{11}}+$
$\frac{1414477}{5591040(n+1/2)^{13}}-\frac{57337}{49152(n+1/2)^{15}}+\frac{118518239}{16711680(n+1/2)^{17}}-\frac{5749691557}{104595456(n+1/2)^{19}}+\frac{91546277357}{173015040(n+1/2)^{21}}$

n=100を代入してDecimalモジュールを使って高精度に計算すると、

3.1415926535 8979323846 2643383279 5028841971 699258655

と真値

3.1415926535 8979323846 2643383279 5028841971 693993751

に対して小数点以下42桁まで一致しました。
実際にオイラーがどのように計算したかわかりませんが、オイラーならこれくらい手で計算できそうな気がします。

from itertools import *
from fractions import Fraction
from decimal import Decimal, getcontext
import time

getcontext().prec = 50

def C(n, m):
    return reduce(lambda x, y: x * (n - y + 1) / y, range(1, m + 1), 1)

def integral(f):
    return [ 0 ] + [ -c / (n - 1) for n, c in enumerate(f[2:], 2) ] + [ 0 ]

def subtract(f, g):
    return [ a - b for a, b in izip(f, g) ]

def add(f, g):
    return [ a + b for a, b in izip(f, g) ]

# 1/(x-a)^n -> 1/x^n + na/x^{n+1} + n(n+1)/2a^2/x^{n+2} + ...
def expand(a, n):
    return [ 0 ] * n + [ C(k - 1, k - n) * a ** (k - n) for k in range(n, M + 1) ]

# f(x) -> f(x + a)
def shift(f, a):
    g = [ 0 ] * (M + 1)
    for k in range(1, M + 1):
        for l, c in enumerate(expand(-a, k)):
            g[l] += c * f[k]
    return g

def next_f(f):
    g = integral(f)
    half = Fraction(1, 2)
    return subtract(f, subtract(shift(g, half), shift(g, -half)))

def compute_fs(N):
    f0 = [ 0, 0, Fraction(1) ] + [ 0 ] * (M - 2)
    fs = [ f0 ]
    f = f0
    for _ in range(N):
        f = next_f(f)
        fs.append(f)
    return fs

def min_dim(f):
    for d, c in enumerate(f):
        if c != 0:
            return d

def TeX_term(c, dim, var):
    d = c.denominator
    n = c.numerator
    
    xd = (var + "^{%d}") % dim if dim > 1 else var
    
    if d == 1:
        return r"\frac{%d}{%s}" % (n, xd)
    elif n > 0:
        return r"\frac{%d}{%d%s}" % (n, d, xd)
    elif n < 0:
        return r"-\frac{%d}{%d%s}" % (-n, d, xd)

def TeX(f, var = "x"):
    min_d = min_dim(f)
    s = ""
    head = True
    for d, c in enumerate(f[min_d:], min_d):
        if c > 0 and not head:
            s += "+" + TeX_term(c, d, var)
        elif c != 0:
            s += TeX_term(c, d, var)
        head = False
    return s

def frac_to_dec(f):
    return Decimal(f.numerator) / f.denominator

def value(f, x):
    return reduce(lambda p, q: p / x + q, map(frac_to_dec, reversed(f)), 0)

def approx_S(N):
    L = 100
    fs = compute_fs(N)
    for k, f in enumerate(fs):
        print "[tex: f_{%d}(x) = %s]" % (k, TeX(f))
    g = reduce(lambda x, y: subtract(x, integral(y)), fs, [ 0 ] * (M + 1))
    print TeX(g, "(n+1/2)")
    res = value(g, L + 1 / Decimal(2))
    return sum(1 / Decimal(n * n) for n in range(1, L + 1)) + res

t0 = time.clock()
N = 10
M = N * 2 + 2
S = approx_S(N)
print (S * 6).sqrt()
print time.clock() - t0