# ゼータ関数

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のゼータ関数のプロットができない
• ということで関数を作るところから…
• こちら積分数値計算すれば？とのコメントがあるので、これでやってみる
• 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()