努力で数論

僕が勉強した内容をまとめておきたくて作ってみました。

素数定理

今日は自分が趣味で作ったpdfを以下に紹介したいと思います。内容は素数分布の近似を与える素数定理についてで、その証明は解析的整数論の一つの金字塔と言えると思います。まだ作成途中であったり、ミスはあるのですがとりあえずというところまで来たので先に公開してしまおうと思いましたｗｗ。英語ではありますが読みたい方はぜひリンクからダウンロードしていただければと思います。

 

https://drive.google.com/file/d/1Lxx0YvlhBJ1NXbhEnZpEp7ajdlQnpEnK/view?usp=sharing

Euler-Mascheroni 定数

本日は級数の分野よりハーモニックシリーズとオイラーの定数について、一つご紹介です。初めに参考文献としては

Euler–Mascheroni constant - Wikipedia

Digamma function - Wikipedia

です。

オイラーの定数 以下の \lim_{n \to \infty}は有限な値に収束し、それはオイラーの定数と呼ばれる。 $$\gamma = \lim _{n \to \infty} \left( - \log n + \sum _{k=1}^{n} \frac {1}{k} \right)$$


定理の主張としては、ハーモニックシリーズは \log nのオーダーで収束しその差はオイラーの定数分しかないということになっております。Wikipediaによると、このオイラーの定数は 0.57721566490153286060651209008240243104215933593992... であり、この記事ではこの極限の値を数値的に計算をする方法を見ていきたいというものです。誤差分を初めに 10^{-6}などとして置いて、極限値がこの誤差分分より小さくなる大きな nを探していくということでもいいのですが（ \epsilon - N論法的なこと）、オイラーの定数は積分計算でも得られるので、まずはその準備としてディガンマ関数というものを用意しましょう。ここではガンマ関数はWeierstrassの無限積表示をされているものと仮定しましょう（https://okenta.hatenablog.jp/entry/2019/07/31/160242）。すると、

$$\psi(z) = \frac{d}{dz} \log (\Gamma(z)) = \frac{\Gamma'(z)}{\Gamma(z)} = \lim _{n \to \infty} (\log n - \sum_{k =0}^ {n} \frac{1}{z+k})$$ $$\psi(1) = \lim _{n \to \infty} (\log n - \sum_{k = 0}^ {n-1} \frac{1}{k+1} - \frac{1}{n+1}) = - \gamma$$

のようにオイラーの定数が得られます。他方、ここでガンマ関数の微分を積分表示の微分と捉えてみると

$$\Gamma'(z) = \frac{d}{dz} \int_{0}^{\infty} t^{z-1} e^{-t} dt = \int_{0}^{\infty} \frac{d}{dz} t^{z-1} e^{-t} dt = \int_{0}^{\infty} t^{z-1} \log t \cdot e^{-t} dt.$$

であり、 \Gamma(1) = 1なのでこれらからオイラーの定数は次のように表されます。

$$\gamma = - \frac{\Gamma'(1)}{\Gamma(1)} = - \int_{0}^{\infty} \log t \cdot e^{-t} dt$$

変数変換すると

$$\gamma = - \int_{0}^{1} \log \left( \log \frac{1}{n} \right) dn$$

なので、この積分を数値計算してみましょう。数値計算はガウスの求積法により精度の良い計算が出来ます。他にも有名な定積分の数値計算としては trapezoidal rule や simpson' s rule、Euler-Maclaurinの公式などがあります。ガウスの求積法はscipyに内蔵される積分計算にも使われているのでかなり良い精度を保証してくれるようです。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def Legendre(N,x):
    P = np.zeros((N+1,len(x)))
    Q = np.zeros((N+1,len(x)))
    P[0], P[1] = 1 + 0*x, x
    Q[1] = 1 + 0*x
    for n in range(1,N):
        P[n+1] = ((2*n+1)*x*P[n] - n*P[n-1])/(n+1)
        Q[n+1] =    (2*n+1)*P[n] + Q[n-1] 
    return P, Q

def quadrature(N,imax=10):
    i = np.arange(N)
    x = np.cos(np.pi*(4*i+3)/(4*N+2))
    for i in range(imax):
        P, Q = Legendre(N,x)
        x -= P[N]/Q[N]
    w = 2/((1-x**2)*(Q[N]**2))
    return x, w 

def integrate(y,w):
    return np.sum(y*w)

def scale(a,b,x,w):
    t  = x*(b-a)/2+(b+a)/2
    dt = (b-a)*w/2
    return t, dt

def f(t):
    return np.log(np.log(1/t))

def DYI_I(a,b,N): 
    x,w = quadrature(N)
    t,dt = scale(a,b,x,w)
    return integrate(f(t),dt)

-DYI_I(0,1,10000)

するとDYI_I は次の値を返してくれます。

$$0.5772156588982373$$

これは文字通りDIYですのでargumentの3個目の値を大きなものにすると精度が高くなります。僕のパソコンでは時間がかかってそれは面倒なので、scipyにお願いしましょう。

import scipy.integrate as integrate
import scipy.special as special

def integrand(n):
    return np.log(np.log(1/n))

eulerconst = - integrate.quad(integrand, 0, 1)[0]
print("gamma = ",eulerconst)

するとこれは

$$gamma = 0.5772156649132887$$

を返します。こちらの方が初めのウィキの値と近くDYI_Iよりは精度がいいですね。本日はオイラーの定数に注目してガウスの求積法を実装してみましたがいかがだったでしょうか。上でfの部分を別の関数に置き換えても定積分の数値計算は簡単にできると思うのでそちらもやってみたい方はぜひ。ガウスの求積法についての参考文献も以下に挙げておきます。

ガウス求積 - Wikipedia

Fibonacci数列についての考察

複素関数と解析接続を学ぶ中で以下の動画を見て、pythonを使ってもタマキさんが動画内で描いていた図と似たようなものを描けそうだなと思ったので、今回はそれをトライしてみました。まずは以下の動画を見て頂けると内容が分かるかと思います。

www.youtube.com

Fibonacci数列とは。

フィボナッチ数列 a_{0} = a_{1} = 1を満たす数列 \{ a_{n} \}について以下の漸化式 a_{n+2} = a_{n+1} + a_{n} を満たす。このとき、フィボナッチ数列は 1,1,2,3,5,8,...となる。


このとき上の漸化式は \phi = \frac{1 + \sqrt{5}}{2}, \psi = \frac{1 - \sqrt{5}}{2}に対して、以下の関係式を満たす。

a_{n+2} - \phi a_{n+1} = \psi (a_{n+1} - \phi a_{n})\\ a_{n+2} - \psi a_{n+1} = \phi (a_{n+1} - \psi a_{n})

\Longleftrightarrow

a_{n+2} - \phi a_{n+1} = \psi^{n} (a_{1} - \phi a_{0}) = \psi^{n} (1 - \phi) = \psi^{n+1} a_{n+2} - \psi a_{n+1} = \phi^{n} (a_{1} - \psi a_{0}) = \phi^{n} (1 - \psi) = \phi^{n+1}

\Longleftrightarrow

(\psi - \phi)a_{n+1} = \psi^{n+1} - \phi^{n+1}

\Longleftrightarrow

\sqrt{5} a_{n+1} = \phi^{n} - \psi^{n}

\Longleftrightarrow

a_{n} = \frac{1}{\sqrt{5}} \bigg \{ \bigg( \frac{1 + \sqrt{5}}{2} \bigg) ^{n} - \bigg( \frac{1 - \sqrt{5}}{2} \bigg)^{n} \bigg \}

Fibonacci数列を係数に持つ級数と解析接続

フィボナッチ数列を係数に持つ級数の関数 \{ a_{n} \}をフィボナッチ数列にするとき、 F(x)を次のように定義する F(x) = a_{0} + a_{1}x + a_{2}x^2 + ... = \sum_{n=0}^{\infty}a_{n}x^{n}


この級数は次の収束半径 Rを持つ。

R = \lim_{n \to \infty} \bigg| \frac{a_{n}}{a_{n+1}} \bigg| = \lim_{n \to \infty} \frac{ \bigg| \frac{1}{\sqrt{5}} \bigg \{ \bigg(\frac{1 + \sqrt{5}}{2} \bigg)^{n} - \bigg(\frac{1 - \sqrt{5}}{2} \bigg)^{n} \bigg \} \bigg| }{ \bigg| \frac{1}{\sqrt{5}} \bigg \{ \bigg(\frac{1 + \sqrt{5}}{2} \bigg)^{n+1} - \bigg(\frac{1 - \sqrt{5}}{2} \bigg)^{n+1} \bigg \} \bigg| } = \lim_{n \to \infty} 2 \frac{ \bigg| \big(1 + \sqrt{5} \big)^{n} - \big(1 - \sqrt{5} \big)^{n} \bigg| }{ \bigg| \big(1 + \sqrt{5} \big)^{n+1} - \big(1 - \sqrt{5} \big)^{n+1} \bigg| } = \lim_{n \to \infty} 2 \frac{ \bigg| \frac{1}{1+\sqrt{5}} - \frac{1}{1+\sqrt{5}} \bigg( \frac{1- \sqrt{5}}{1+ \sqrt{5}} \bigg)^{n} \bigg| }{ \bigg| 1 - \bigg(\frac{1 - \sqrt{5}}{1 + \sqrt{5}} \bigg)^{n+1} \bigg| } = \lim_{n \to \infty} 2 \frac{ \bigg| \frac{1}{1+\sqrt{5}}\bigg| }{|1|} = \frac{2}{1+\sqrt{5}} < 1

この時、この収束半径内で級数は収束することが言えます。さて、他方

xF(x) = a_{0}x + a_{1}x^2 + a_{2}x^3 + ...

なので

F(x) + xF(x) = a_{0} + (a_{0} + a_{1})x + (a_{1} + a_{2})x^2 + ... = a_{1} + a_{2}x + a_{3}x^2 + ... = \frac{1}{x} (a_{1}x + a_{2}x^2 + a_{3}x^3 + ... ) = \frac{1}{x} (F(x) - a_{0}) = \frac{1}{x} (F(x) - 1)

\Longrightarrow

F(x) = \frac{1}{1-x-x^2}

を満たす。 F(x)は x=\frac{1 \pm \sqrt{5}}{2}を除く任意の点で正則で、先の級数表示の解析接続したものとなっています。

pythonでグラフを作る

さて、この関数の実部と虚部を複素平面上からの写像で3Dのグラフに表す。これはミックスされた形でタマキさんの動画でも流れているものかと思います。

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

def three_D_plots(func,picname):
    
    x, y = np.meshgrid(np.linspace( -5, 5, 500), np.linspace( -5, 5, 500))
    z = x + y*1j
    w = func(z)
    fig1 = plt.figure()
    ax = Axes3D(fig1)
    ax.plot_surface(x, y, w.real, color="blue")
    ax.set_zlim3d(0, 15.)
    ax.autoscale(False)
    plt.axis([-6., 6., -6., 6.])
    plt.title("real_part")
    
    fig2 = plt.figure()
    ax = Axes3D(fig2)
    ax.plot_surface(x, y, w.imag, color="orange")
    ax.set_zlim3d(0, 15.)
    ax.autoscale(False)
    plt.axis([-6., 6., -6., 6.])
    plt.title("imag_part")
    
    plt.savefig(picname)
    plt.show()

def fibonacci(s):
    return 1/(1-s-s**2)

three_D_plots(fibonacci,"fibonacci_real.png","fibonacci_img.png")

f:id:okenta:20190801021155p:plain

f:id:okenta:20190801021223p:plain

さて、今までやってきた通りこれを今度は絶対値付きのグラフにしたり、2Dの図を作ったりすると

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from mpl_toolkits.mplot3d import Axes3D
from scipy.special import gamma

def three_D_plot(func,picname):
    x, y = np.meshgrid(np.linspace( -5, 5, 500), np.linspace( -5, 5, 500))
    z = x + y*1j
    w = func(z)
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot_surface(x, y, np.abs(w))
    ax.set_zlim3d(0, 15.)
    ax.autoscale(False)
    plt.axis([-6., 6., -6., 6.])
    plt.savefig(picname)
    plt.show()
    
def plot_complex(complex_func,picname1,picname2):

    x, y = np.meshgrid(np.linspace(-4, 4, 400), np.linspace(-4, 4, 400))
    z = x + y*1j

    angles = (np.angle(complex_func(z)) + 2*np.pi) % (2*np.pi)
    plt.figure(figsize=(6,5))
    im_1 = plt.pcolor(x, y, angles, cmap="hsv", vmin=0, vmax=6)
    plt.colorbar(im_1)
    plt.title('argument')
    plt.savefig(picname1)

    plt.figure(figsize=(6,5))
    im_2 = plt.pcolor(x, y, np.abs(complex_func(z)), vmin=0, vmax=4)
    plt.title('absolute value')
    plt.colorbar(im_2)
    plt.savefig(picname2)

    plt.show()
    

three_D_plot(fibonacci,"fibonacci3D.png")
plot_complex(fibonacci,"fibonacci_arg.png","fibonacci_abs.png")

f:id:okenta:20190801021939p:plain

f:id:okenta:20190801021955p:plain

f:id:okenta:20190801022011p:plain

のようになります。タマキさんの今回の動画の面白いところはサムネにもある式がどういう意味をもつのかということだと思うので解析接続についてや無限級数の収束半径などは飛ばしてもとても楽しめます！その他Riemann Zeta関数についても似た内容の動画がアップされているかと思うので、ぜひ興味のある方はそちらもチェックなさってみてはいかがでしょう。

複素関数のグラフについて、様々な例

今まで作ってきた複素関数のプロット例をいくつか紹介したいと思います。一応、グラフのプロット関数をもう一度書いておきます。

複素関数のグラフを描く関数。

まず与えられた複素関数の絶対値付きの値を返す、3Dのグラフを描くプロット関数を下に記します。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from mpl_toolkits.mplot3d import Axes3D
from scipy.special import gamma

def three_D_plot(func,picname):
    x, y = np.meshgrid(np.linspace( -5, 5, 500), np.linspace( -5, 5, 500))
    z = x + y*1j
    w = func(z)
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot_surface(x, y, np.abs(w))
    ax.set_zlim3d(0, 15.)
    ax.autoscale(False)
    plt.axis([-6., 6., -6., 6.])
    plt.savefig(picname)
    plt.show()

次にこの複素関数を偏角ごとに２Dに描く関数と、絶対値付きの値をレベルカーブごとに２Dにプロットする関数を下に書いておきます。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.special import gamma

def plot_complex(complex_func,picname1,picname2):

    x, y = np.meshgrid(np.linspace(-4, 4, 400), np.linspace(-4, 4, 400))
    z = x + y*1j

    angles = (np.angle(complex_func(z)) + 2*np.pi) % (2*np.pi)
    plt.figure(figsize=(6,5))
    im_1 = plt.pcolor(x, y, angles, cmap="hsv", vmin=0, vmax=6)
    plt.colorbar(im_1)
    plt.title('argument')
    plt.savefig(picname1)

    plt.figure(figsize=(6,5))
    im_2 = plt.pcolor(x, y, np.abs(complex_func(z)), vmin=0, vmax=4)
    plt.title('absolute value')
    plt.colorbar(im_2)
    plt.savefig(picname2)

    plt.show()

いくつかの具体例。

例①

このプロット関数達を用いてまずガンマ関数を前回のようにプロットしてみる。ガンマ関数自体については一つ前のの記事を読んでみて下さい。（Gamma関数をお絵描き - 努力で数論）

three_D_plot(gamma,picname="gamma3D.png")

するとグラフは以下のようになる。

f:id:okenta:20190731152732p:plain

plot_complex(gamma,"gamma_arg.png","gamma_abs.png")

するとグラフは以下のようになる。

f:id:okenta:20190731155324p:plain

f:id:okenta:20190731155340p:plain

例②

さて別の複素関数でトライしてみよう。

例2 任意の z \in \mathbb{C}に対して $$f(z) = z$$ 


def identity(z):
    return z

three_D_plot(identity,picname="identity3D.png")

するとグラフは以下のようになる。

f:id:okenta:20190731223106p:plain

plot_complex(identity,"identity_arg.png","identity_abs.png")

するとグラフは以下のようになる。

f:id:okenta:20190731223524p:plain

f:id:okenta:20190731223541p:plain

例③

次の例は逆数を返す関数。ちなみにこの関数は上の例2とは違って z=0で極を持つのでそこを除けば正則な関数となっている。こういったビジュアルの要素であるグラフを通じ、極の理解が深まるはずである。

例3 0以外の任意の z \in \mathbb{C}に対して $$f(z) = \frac{1}{z}$$


まず絶対値付きの場合は

def inverse(z):
    return 1/z

three_D_plot(inverse,picname="inverse3D.png")

f:id:okenta:20190731224055p:plain

となり、次に2Dの図は

plot_complex(inverse,"inverse_arg.png","inverse_abs.png")

f:id:okenta:20190731224227p:plain

f:id:okenta:20190731224245p:plain

となる。

例④

複素指数関数は次のように定義でき、任意の z \in \mathbb{C}に対して正則である。

例4 任意の z \in \mathbb{C}に対して $$f(z) = e^{z}$$ 


three_D_plot(np.exp,picname="exp3D.png")

f:id:okenta:20190731224617p:plain

plot_complex(np.exp,"exp_arg.png","exp_abs.png")

f:id:okenta:20190731224906p:plain

f:id:okenta:20190731224930p:plain

グラフは上のようになる。

例⑤

複素三角関数の一例として sin(z)のグラフを描く。これも任意の z \in \mathbb{C}に対して正則である。

例5 任意の z \in \mathbb{C}に対して $$f(z) = sin(z)$$ 


three_D_plot(np.sin,picname="sin3D.png")

するとグラフは以下のようになる。

f:id:okenta:20190731225240p:plain

plot_complex(np.sin,"sin_arg.png","sin_abs.png")

するとグラフは以下のようになる。

f:id:okenta:20190731225334p:plain

f:id:okenta:20190731225355p:plain

Gamma関数をお絵描き

本日もpythonで描いたグラフを乗っけていこうと思います。第二回はガンマ関数のグラフについてです。

Gamma関数とは。

まずガンマ関数の定義ですが、これもWikipedia（ガンマ関数 - Wikipedia）を参照されるのがベターでしょうが、一応一般的によくある定義と解析接続したものであるWeierstrassの表示を紹介したいと思います。

Gamma関数 Re(z) > 0を満たす z \in \mathbb{C}に対して $$\Gamma (z) = \int_0^{\infty} t^{z-1} e^{-t}dt $$


これも定義域が複素平面全てではないのでグラフを描くにあたっては都合がよくありません。そこで次の解析接続後の、Weierstrassの無限積表示の定義も書いておきます。

Gamma関数（Weierstrassの無限積表示） 原点と負の整数における極を除く任意の z \in \mathbb{C}に対して $$\Gamma (z) = \bigg[ z e^{\gamma z} \prod _{r=1}^{\infty} \bigg( 1 + \frac{z}{r} \bigg) e^{-z/r} \bigg]^{-1}$$ ただし \gammaはオイラーの定数である。


この表示の逆数、即ち \frac{1}{\Gamma(x)}は任意の z \in \mathbb{C}において正則（複素微分可能で導関数が連続）です。他にもガンマ関数に関する定理、命題は次のページによくまとまっています。

Gamma Function -- from Wolfram MathWorld

pythonを使って描く。

まず zに任意の実数を代入してガンマ関数を返す時のグラフが次のコードである。

from scipy.special import gamma, zeta
from matplotlib import pyplot as plt
%matplotlib inline
import numpy as np

x = np.linspace( -5., 5., 500)
y = gamma(x)
plt.axis([-5., 5., -5., 5.])
plt.plot(x, y,"o")
plt.savefig("gamma_real.png")
plt.show()

すると得られる画像が次の通り。

f:id:okenta:20190731152247p:plain

さて次に絶対値付きガンマ関数の3Dのグラフのコードが次の通りです。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from mpl_toolkits.mplot3d import Axes3D
from scipy.special import gamma

def three_D_plot(func,picname):
    x, y = np.meshgrid(np.linspace( -5, 5, 500), np.linspace( -5, 5, 500))
    z = x + y*1j
    w = func(z)
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot_surface(x, y, np.abs(w))
    ax.set_zlim3d(0, 15.)
    ax.autoscale(False)
    plt.axis([-6., 6., -6., 6.])
    plt.savefig(picname)
    plt.show()

three_D_plot(gamma,picname="gamma3D.png")

するとグラフは次のようになります。

f:id:okenta:20190731152732p:plain

この図からも分かるように虚部がゼロの部分の断面図が先の図で実軸対象に折り返したものに一致することが確認できます。次にこのグラフの偏角の様子とレベルカーブを複素平面にプロットした図が次のコードが与える2つのグラフです。1つ目のカラフルな図が偏角によって色が変化していて、2つ目がレベルカーブです。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.special import gamma

def plot_complex(complex_func,picname1,picname2):

    x, y = np.meshgrid(np.linspace(-4, 4, 400), np.linspace(-4, 4, 400))
    z = x + y*1j

    angles = (np.angle(complex_func(z)) + 2*np.pi) % (2*np.pi)
    plt.figure(figsize=(6,5))
    im_1 = plt.pcolor(x, y, angles, cmap="hsv", vmin=0, vmax=6)
    plt.colorbar(im_1)
    plt.title('argument')
    plt.savefig(picname1)

    plt.figure(figsize=(6,5))
    im_2 = plt.pcolor(x, y, np.abs(complex_func(z)), vmin=0, vmax=4)
    plt.title('absolute value')
    plt.colorbar(im_2)
    plt.savefig(picname2)

    plt.show()

plot_complex(gamma,"gamma_arg.png","gamma_abs.png")

f:id:okenta:20190731155324p:plain

f:id:okenta:20190731155340p:plain

上の図から実部が負の整数で虚部がゼロの時、黄色っぽくなっていることからここで極を取っていることが分かりますね。さて今まで書いてきたガンマ関数ですが今回はpythonのscipyというライブラリにもともとあるガンマ関数を使っています。もちろんmpmathにもガンマ関数はあって前回も書いた通り、cplotという便利なプロット用の関数も用意されてるのでそちらを使っても同じものが描けます。最後にそれを一応乗せておきましょう。

from mpmath import *
cplot(gamma, points=50000, file="gamma2D.png")

グラフとしては上の自前のものをミックスしたような形になっているのが得られます。ちなみに黒の濃淡が絶対値の大小を表しています。

f:id:okenta:20190731160152p:plain

Riemann Zeta関数をお絵描き

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

１．Riemann Zeta関数とは。

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

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


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

解析接続後のRiemann Zeta関数 任意の s \in \mathbb{C}に対して $$\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関数はユニークに決まる関数となっているので、この関数に対して以下零点をプロットしていきます。

また sに負の偶数（ n \in \mathbb{N}に対して s = -2n）を代入すると \zeta(-2n) = 0ということが分かりますが、これは自明な零点と言われています。理由としてはRiemann Zeta関数がBernoulli数を使って表現できるところにあり（これはWikipediaの上の記事で確認できる）、Bernoulli数（ B_{n}）は nが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")

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

f:id:okenta:20190731020832p:plain

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

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

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")

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

f:id:okenta:20190731021951p:plain

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