Riemann Zeta関数をお絵描き

これから暇を持て余している時を見つけて、今までpythonを使って書いたグラフとかをとりあえずは乗っけていこうかなと思います。第一回はRiemann Zeta関数の零点についてです。

１．Riemann Zeta関数とは。

いろんなところで定義や説明がかかれてるはずなので、例えばWikipediaや（リーマンゼータ関数 - Wikipedia）次の記事（リーマンのゼータ関数で遊び倒そう (Ruby編) - tsujimotterのノートブック）を参照された方がより多くの内容を見れると思います（実際僕はこの２つの記事で結構勉強させて頂きました、ありがとうございます）が、一応素人なりに書いておきます。

Riemann Zeta関数 を満たすに対して $$\zeta(s) = \sum_{n=0}^{\infty} \frac{1}{n^{s}}$$

これがRiemann Zeta関数と言われるものですが、今日のテーマはとなるようなを今はちょっと考えてみたいということです。本当は解析接続が何かということや解析接続を施す過程も書いておかなければ零点を考える上では十分ではないとも思ったのですが、とりあえずは上の定義域よりも広いところで（複素平面全体で）定義される同じような関数として以下の形もあげられます。

解析接続後のRiemann Zeta関数 任意のに対して $$\zeta(s) = \frac{1}{1-2^{1-s}} \sum_{n=0}^{\infty} \frac{1}{2^{n+1}} \sum_{k=0}^{n} (-1)^{k} {n \choose k} (k+1)^{-s}$$

詳しくはRiemann Zeta Function -- from Wolfram MathWorldをお読みください。すべてには目を通してないですがいろいろ書かれてて面白いです。一致の定理により解析接続した後のRiemann Zeta関数はユニークに決まる関数となっているので、この関数に対して以下零点をプロットしていきます。

またに負の偶数（に対して）を代入するとということが分かりますが、これは自明な零点と言われています。理由としてはRiemann Zeta関数がBernoulli数を使って表現できるところにあり（これはWikipediaの上の記事で確認できる）、Bernoulli数（）は が3以上の奇数の時にゼロになること（ベルヌーイ数 - Wikipedia）からある程度すぐ分かるということだと思います。

２．pythonを使って描く。

まず僕がどの環境でpythonを使っているかというと、jupyter notebookの中で使用しています。jupyterを使うとpythonの入出力とlatexを用いた文書作成の両方が一気にでき、出来上がったものが「データ＋文」ということで結構見やすくなります。僕は、参考として挙げる次の記事にもあるようにpythonとAnacondaをダウンロードし、そこからjupyter notebookを使っています。

データ分析で欠かせない！Jupyter Notebookの使い方【初心者向け】 | TechAcademyマガジン

これらは基本的にはすべて無料なので学生には非常に有難いですね。

さてようやくグラフを描くのですが、今回はRiemann Zeta関数についてはすでにあるものを使用しました。mpmathにはzetaという関数が用意されてますし、なんならcplotという僕が作ったものより性能がいいコードも準備されてるのでわざわざ自前でレベルカーブを描くほどでもないんですが、とりあえず今はRiemann Zeta関数についてはあるものを使ってグラフの書き方は自分で行うこととしました。またmpmathに慣れておらず、途中無理やり新しい配列に置き換えてるのでコードはかなり汚いかもしれません（汗）。

始めは絶対値付きのRiemann Zeta関数を3Dで見たときのグラフのコードです。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from mpmath import *
from mpl_toolkits.mplot3d import Axes3D

def three_D_plot(func,picname):
    
    x, y = np.meshgrid(np.linspace( -5, 5, 100), np.linspace( -11, 11, 220))
    z = x + y*1j
    w = np.abs(np.vectorize(func)(z))
    new = np.zeros((x.shape[0],x.shape[1]))
    for i in range(x.shape[0]):
        for j in range(x.shape[1]):
            new[i][j]=w[i][j]

    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot_surface(x, y, np.abs(new))
    ax.set_zlim3d(0, 2.)
    ax.autoscale(False)
    plt.axis([-6., 6., -12., 12.])
    plt.savefig(picname)
    plt.show()

three_D_plot(zeta,picname="zeta3D.png")

するとこんな絵が描けます。

はっきり言ってこの絵はあまり出来は良くなくて、例えばこの方（リーマンのゼータ関数の数値計算コード(複素平面) | シキノート）のように描ければベストですが今は描けていません。

さてこの上のグラフのレベルカーブを複素平面にプロットしたものが以下のコードです。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from mpmath import *

def contour_coloring(func, picname):
    
    x, y = np.meshgrid(np.linspace( -20, 10, 300), np.linspace( -40, 40, 400))
    z = x + y*1j
    w = np.abs(np.vectorize(func)(z))
    new = np.zeros((x.shape[0],x.shape[1]))
    for i in range(x.shape[0]):
        for j in range(x.shape[1]):
            new[i][j]=w[i][j]
    
    plt.figure(figsize=(6,8))
    cs = plt.pcolor(x, y, new, vmin=0, vmax=3)
    plt.colorbar(cs)
    plt.savefig(picname)
    plt.show()

contour_coloring(zeta,"RH.png")

すると次のような図が描けます。

なかなかきれいですよね。リーマン予想というのは「任意の非自明な零点（上のグラフだと濃い青の点）が、実部がの線上に乗っている」というもので、この範囲ではその主張は信じることができそうだなといえます。数学者によってもこの予想が正しいと信ずるかどうかはいろいろ議論があるようですが、少なくともこの狭い範囲ではそれが成り立ちそうだなと予想できますね。どちらのコードもグラフを描くようの関数を先に定義しているので、これをちょっと変えれば別の関数に対しても同じようなグラフが描けます。