2013-10-09

ゼータ関数

素数の分布に関するリーマン予想
オックスフォード白熱教室　第１回「素数の音楽を聴け」を録画したものを、家人が洗濯物を干す時間つぶしに観ていた(10分観ては中断、いつの間にか八重の桜の録画に切り替わったりする…)。僕は横で弁当箱をしまいつつ聴いていた。
ある自然数までにいくつの素数が存在するか
それを表す関数と関係するゼータ関数
「うまい具合にばらばら」に存在することと関係するんですよ、みたいなコメントも聞こえてきた
素数、自然数だけ、ばらばらになる、それがルール…
フィボナッチ数列が生物現象に現れるように、リーマン予想とかゼータ関数の解とかが生物の状態の背後にある、という風にストーリが作れるような気がしたけれど、すぐにはなんともわからない
ひとまず、リーマン予想、実部が0.5というのをRでなぞったりすると思うことがあるかも…と言うことで、描こうと思ったら…
Rにはゼータ関数がない？

x <- seq(0, 20, len=1001)
z <- 0.5 + x*1i
fr <- Re(zeta(z))
fi <- Im(zeta(z))
fa <- abs(zeta(z))
plot(x, fa, type="n", xlim = c(0, 20), ylim = c(-1.5, 2.5),
     xlab = "Imaginary part (on critical line)", ylab = "Function value",
     main = "Riemann's Zeta Function along the critical line")
lines(x, fr, col="blue")
lines(x, fi, col="darkgreen")
lines(x, fa, col = "red", lwd = 2)
points(14.1347, 0, col = "darkred")
legend(0, 2.4, c("real part", "imaginary part", "absolute value"),
       lty = 1, lwd = c(1, 1, 2), col = c("blue", "darkgreen", "red"))
grid()

Mathematicaにはあるらしい(こちらの図はMathematicaで描いたらしい)
RでもVGAMパッケージにzeta()関数があるが、実部が１以上に限定しているという…。それでは実部が0.5のゼータ関数のプロットができない
ということで関数を作るところから…
こちらで積分を数値計算すれば？とのコメントがあるので、これでやってみる
$\zeta(s) = \frac{2^{s-1}}{s-1} - 2^s\int_{0}^{\infty} \frac{\sin(s\times \arctan(x))}{(1+x^2)^{\frac{s}{2}}(e^{\pi x}+1)} dx$
Rにはゼータ関数がない？
- →pracmaパッケージにある

- 以下のソースとおおまかに合う(精度はpracmaパッケージの方が上のはず)

# これじゃ用が足りないようだ
my.zeta <- function(s,n){
	sum(1/(1:n)^s)
}
# こちらが積分式の数値積分版
my.zeta.2 <- function(s,n=1000,N=10000){
	term1 <- 2^(s-1)/(s-1)
	term2 <- 2^s
	x <- seq(from=0,to=n,length=N)
	term3 <- sum(sin(s*atan(x))/((1+x^2)^(s/2)*(exp(pi*x)+1))/N*n)
	term1-term2*term3
}


s <- seq(from=-30,to=30,length=1000)
ss <- complex(real=0.5,imaginary=s)
zs <- rep(0,length(ss))

for(i in 1:length(ss)){
	zs[i] <- my.zeta.2(ss[i])
}
matplot(s,cbind(Re(zs),Im(zs)),type="l",col=c(2,3))
abline(h=0)

x <- seq(0, 20, len=1001)
z <- 0.5 + x*1i
fr <- Re(zeta(z))
fi <- Im(zeta(z))
fa <- abs(zeta(z))
plot(x, fa, type="n", xlim = c(0, 20), ylim = c(-1.5, 2.5),
     xlab = "Imaginary part (on critical line)", ylab = "Function value",
     main = "Riemann's Zeta Function along the critical line")
lines(x, fr, col="blue")
lines(x, fi, col="darkgreen")
lines(x, fa, col = "red", lwd = 2)
points(14.1347, 0, col = "darkred")
legend(0, 2.4, c("real part", "imaginary part", "absolute value"),
       lty = 1, lwd = c(1, 1, 2), col = c("blue", "darkgreen", "red"))
grid()