結果は1億個の素数リストを作成し、それを用いてオイラー積の計算を行い、ζ(2)=1.6449340668…と11桁の正しい値を得ることができました。
オイラーが行わなかった計算にあえてチャレンジしてみましたが、今回はオイラーが行ったゼータの数値計算の再現を紹介しようと思います。
目次
コード1「python_math15.ipynb」(Jupyter Notebook用）
まず準備として、ζ(2)からζ(10)までの数値計算をmpmathライブラリによって実行しておきます。
ζ(s)はmpmath.zeta(s)により得られます。
関数sympy.N()は任意の桁で値を評価します。
sympy.N(mpmath.zeta(s), 30)
とすれば、ζ(s)を30桁の精度で計算してくれます。
オイラー・マクローリンの公式のコーディング
今から300年前、オイラーが手計算によりゼータの数値計算をいかに行ったか。
その1つの答えがオイラー・マクローリンの公式です。
左辺のシグマの値を右辺の式により数値計算します。まず自然数nを引数とする数列f(n)を実数xを引数とする関数f(x)とします。
右辺第1項はf(x)の定積分です。
ざっくり解説するならば、左辺シグマは数列の和なのでグラフに表せば短冊を集めたガタガタ形をしています。
それを右辺ではまずnを実数xとして、nの区間と同じ区間で曲線が囲む面積でガタガタ形を近似してしまいます。
右辺の残りの2項でその補正を行っているということです。最後の項は誤差項です。補正はk階微分と関・ベルヌーイ数を用いて計算します。
使用した関数
- f(x)積分区間[a, b]における定積分 sympy.integrate(f,(x,a,b))
- 関・ベルヌーイ数B(n) sympy.bernoulli(n)
- 階乗n! sympy.factorial(n)
- 関数fをxでk階微分 sympy.diff(f, x, k)
- 微分係数の計算 数値代入subs(x, n)メソッド
これらを用いてオイラー・マクローリンの公式をコーディングします。
オイラー・マクローリンの公式を適用するにあたり、右辺の最初の数項だけ除いたところから公式を適用するのが誤差を防ぐポイントです。
そこで最初の10項だけは数列をそのまま数値計算行います。その部分が変数Aです。変数I2が右辺第2項部分、そして変数Sが第3項部分のシグマの値です。
定義する関数emzeta(s, n)のsはζ(s)、第2引数nは第３項部分シグマの項数です。
このJupyter Notebook用のプログラムファイルは次からダウンロードできます
https://drive.google.com/file/d/14GBXskRXn_o0XamiW7zto2uFIwGWDCsk/view?usp=sharing
結果
以上の結果を次にまとめました。
第3項部分シグマの項数を10万とした場合の結果です。16桁ほどの精度で正しいことがわかります。
ところで、Pythonによる計算からはわかりづらいのですが、なぜオイラー・マクローリンの公式が手計算を可能にするのか。
それは、nが100000としたところで、右辺の計算は初項と末項（n=100000）の2つだけの計算に帰着していることから、さほど計算量が爆発的に大きくならないということです。
それにしても、実際にオイラー・マクローリンの公式を手計算で実行しようとすると、多くの準備計算が必要なのは事実です。
しかし、それでも電子計算機がなかった時代にゼータの数値計算を可能にするこの公式には絶大な偉力があります。
Pythonでコーディングしてみると、オイラー・マクローリンの公式も簡単に見えてくるのが不思議です。
ぜひこのコードの関数部分に自由に関数を入れて、公式を試してみてください。
桜井進デザイン 数式Tシャツ
私はオイラー・マクローリンの公式にとても感動したので、Tシャツにデザインして公式を着ています。
πをはじめゼータなど10種類デザインしました。
シグマヴァージョン オイラー・マクローリンの公式
ゼータ関数ヴァージョン リーマンの関数等式
