Notebook

Euler変換によるζ(1.5)の計算¶

黒木玄
2021-03-25

In [1]:

using Plots
using SpecialFunctions
SpecialFunctions.lgamma(x::Real) = logabsgamma(x)[1]
logrelerr(x, y) = log10(abs(x/y - 1))

Out[1]:

logrelerr (generic function with 1 method)

Euler変換¶

交代級数のEuler変換については

を参照. さらに

ディリクレ級数の滑らかなカットオフ

も参考になるはず. 交代ゼータ函数は極を持たないので, 滑らかなカットオフが正確な値に収束してくれる.

交代ゼータ函数 $\eta(s)$ は次のように定義される:

$\eta(s) = \sum_{n=1}^\infty \frac{(-1)^{n-1}}{n^s}$

正確にはこれの解析接続を交代ゼータ函数と呼ぶことにする. $\eta(s)$ は複素平面上で極を持たない. Riemannのゼータ函数

$\zeta(s) = \sum_{n=1}^\infty \frac{1}{n^s}$

との関係は,

$\eta(s) = \sum_{n=1}^\infty \frac{1}{n^s} - 2 \sum_{n=1}^\infty \frac{1}{(2n)^s} = \zeta(s) - 2\frac{1}{2^s}\zeta(s) = (1 - 2^{1-s})\zeta(s)$

なので, $\eta(s)$ が計算できれば $\zeta(s)$ は

$\zeta(s) = \frac{\eta(s)}{1 - 2^{1-s}}$

で計算できる.

$\displaystyle \sum_{k=0}^\infty (-1)^k a_k$ のEuler変換の結果は

$\lim_{L\to\infty} \sum_{k=0}^{L-1} w^{(L)}_k (-1)^k a_k, \quad w^{(L)}_k = \sum_{n=k}^{L-1} 2^{-n-1}\binom{n}{k}.$

$w^{(L)}_k$ を次の $k$ の滑らかな函数で近似できる:

$w_{\mathrm{smooth}}^{(L)}(k) = \frac{1}{2}\mathrm{erfc}\left(\frac{k-L/2}{\sqrt{L/2}}\right), \quad \mathrm{erfc}(x) = \int_x^\infty\frac{e^{-t^2}}{\sqrt{\pi}}\,dt$

http://www.kurims.kyoto-u.ac.jp/~ooura/papers/CET-R99.pdf のp.2を見よ.

In [2]:

binomover2np1(n, k) = exp(lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1) - (n+1)*log(2))
eulerweight(L, k) = sum(binomover2np1(n, k) for n in k:L-1; init=0.0)

struct EulerWeight{T} W::T end

function EulerWeight(; Lmax=2^6)
    W = eulerweight.(1:Lmax, (0:Lmax-1)')
    EulerWeight(W)
end

Lmax = 2^6
(w::EulerWeight)(L, k) = w.W[L, k+1]
w = EulerWeight(; Lmax)

Out[2]:

EulerWeight{Matrix{Float64}}([0.5 0.0 … 0.0 0.0; 0.75 0.25 … 0.0 0.0; … ; 0.9999999999999999 0.9999999999999999 … 1.084202172485508e-19 0.0; 0.9999999999999999 0.9999999999999999 … 3.523657060577715e-18 5.42101086242753e-20])

In [3]:

w_smooth(L, k) = (1/2)*erfc((k-L/2)/sqrt(L/2))

k = 0:Lmax-1
x = range(extrema(k)...; length=1001)
scatter(k, w.(Lmax, k); label="Euler's weights")
plot!(x, w_smooth.(Lmax, x); label="smooth version")

Out[3]:

In [4]:

eta_naive(s, L=2^6) = sum((-1)^(n-1)/n^s for n in 1:L )
eta_euler(s, L=2^6; w=w) = sum(w(L, n-1)*(-1)^(n-1)/n^s for n in 1:L)
etasmooth(s, L=2^6; w=w_smooth) = eta_euler(s, L; w)
etacutoff(s, L=2^6; N=30) = sum(exp(-(n/N)^2)*(-1)^(n-1)/n^2 for n in 1:L)
@doc eta

Out[4]:

eta(x)

Dirichlet eta function $\eta(s) = \sum^\infty_{n=1}(-1)^{n-1}/n^{s}$ .

In [5]:

L = Lmax
eta_naive_val = eta_naive(1.5, L)
eta_euler_val = eta_euler(1.5, L)
eta_smoothval = etasmooth(1.5, L)
eta_exact_val = eta(1.5)
@show eta_naive_val
@show eta_euler_val
@show eta_smoothval
@show eta_exact_val;

eta_naive_val = 0.7641819041811663
eta_euler_val = 0.7651470246254075
eta_smoothval = 0.7651470246254077
eta_exact_val = 0.7651470246254081

In [6]:

zeta_naive_val = eta_naive_val/(1 - 2^(1-1.5))
zeta_euler_val = eta_euler_val/(1 - 2^(1-1.5))
zeta_smoothval = eta_smoothval/(1 - 2^(1-1.5))
zeta_exact_val = eta_exact_val/(1 - 2^(1-1.5))
@show zeta_naive_val
@show zeta_euler_val
@show zeta_smoothval
@show zeta_exact_val;

zeta_naive_val = 2.6090802213754345
zeta_euler_val = 2.6123753486854864
zeta_smoothval = 2.6123753486854873
zeta_exact_val = 2.612375348685488

In [7]:

ENV["LINES"] = 100
L = 1:Lmax
eta_naive_val = eta_naive.(1.5, L)
eta_euler_val = eta_euler.(1.5, L)
eta_smoothval = etasmooth.(1.5, L)
eta_exact_val = eta(big"1.5")
plot(; legend=:bottomright, xlabel="L")
plot!(L, eta_euler_val; label="Euler transform")
plot!(L, eta_smoothval; label="smooth version", ls=:dash)
plot!(L, eta_naive_val; label="naive sum", ls=:dashdot)
hline!([eta_exact_val]; label="exact value", c=:black, ls=:dot)
title!("eta(1.5)"; titlefontsize=10)

Out[7]:

In [8]:

naive_logrelerr = logrelerr.(eta_naive_val, eta_exact_val)
euler_logrelerr = logrelerr.(eta_euler_val, eta_exact_val)
smoothlogrelerr = logrelerr.(eta_smoothval, eta_exact_val)
plot(; legend=:right, xlabel="L", ytick=-100:100)
plot!(L, euler_logrelerr; label="Euler transform")
plot!(L, smoothlogrelerr; label="smooth version", ls=:dash)
plot!(L, naive_logrelerr; label="naive sum", ls=:dashdot)
title!("log10 relative error"; titlefontsize=10)

Out[8]:

Euler-Maclaurinの公式¶

近似公式

$\zeta(s) \approx \sum_{n=1}^{a-1}\frac{1}{n^s} - \frac{a^{1-s}}{1-s} + \frac{1}{2a^s} + \sum_{k=2}^n \frac{B_k}{ka^{s+k-1}}\binom{-s}{k-1}$

の詳細については

13 Euler-Maclaurinの和公式

を参照せよ. 実際の計算時には $k\geqq 2$ かつ $k$ が奇数のとき $B_k = 0$ になることに注意せよ.

$\zeta(s) \approx \sum_{n=1}^{a-1}\frac{1}{n^s} - \frac{a^{1-s}}{1-s} + \frac{1}{2a^s} + \sum_{k=1}^m \frac{B_{2k}}{2k a^{s+2k-1}}\binom{-s}{2k-1}$

In [9]:

mydiv(a, b) = a/b
mydiv(a::Rational, b::Rational) = a//b
mydiv(a::Integer, b::Integer) = a ÷ b

function mybinom(n, k)
    k < 0 && return zero(n)
    iszero(k) && return one(n)
    b = one(n)
    for j in 1:k
        b = mydiv(b*(n - k + j), j)
    end
    b
end

struct Bernoulli{T} B::T end
function Bernoulli(nmax::Integer = 30, T::Type = Int64)
    B = zeros(Rational{T}, nmax + 1)
    B[begin]   = 1     # B_0
    B[begin+1] = -1//2 # B_1
    for n in T(2):2:nmax
        B[begin+n] = -1//(n+1) * sum(mybinom(n+1, j)*B[begin+j] for j in 0:n-1)
    end
    Bernoulli(B)
end
(B::Bernoulli)(n) = B.B[begin+n]
B = Bernoulli()
@show B.(0:30)

function zeta_em(s::AbstractFloat, a=10, m=9; B=B)
    n = 2m
    z = sum(1/n^s for n in 1:a-1; init=zero(s))
    z += -a^(1 - s)/(1 - s)
    n == 0 && return z
    z += 1/(2*a^s)
    n == 1 && return z
    z += -sum(typeof(s)(B(k))/(k*a^(s+k-1))*mybinom(-s, k-1) for k in 2:2:n)
    z
end

zeta_em(s::Real, a=10, m=9; B=B) = zeta_em(float(s), a, m; B)

B.(0:30) = Rational{Int64}[1//1, -1//2, 1//6, 0//1, -1//30, 0//1, 1//42, 0//1, -1//30, 0//1, 5//66, 0//1, -691//2730, 0//1, 7//6, 0//1, -3617//510, 0//1, 43867//798, 0//1, -174611//330, 0//1, 854513//138, 0//1, -236364091//2730, 0//1, 8553103//6, 0//1, -23749461029//870, 0//1, 8615841276005//14322]

Out[9]:

zeta_em (generic function with 6 methods)

In [10]:

zeta_em____2 = zeta_em(2, 10, 9; B)
zeta_exact_2 = zeta(2)
@show zeta_em____2
@show zeta_exact_2;

zeta_em____2 = 1.6449340668482264
zeta_exact_2 = 1.6449340668482273

In [11]:

zeta_em____1_5 = zeta_em(1.5, 9, 8; B)
zeta_exact_1_5 = zeta(1.5)
@show zeta_em____1_5
@show zeta_exact_1_5;

zeta_em____1_5 = 2.6123753486854886
zeta_exact_1_5 = 2.6123753486854886

In [12]:

function plot_em_relerrs(s)
    a = 1:10
    zeta_em____val = zeta_em.(s, a, a .- 1; B)
    zeta_exact_val = zeta(BigFloat(s))
    em_logrelerr = logrelerr.(zeta_em____val, zeta_exact_val)
    plot(2a, em_logrelerr; label="", xlabel="number of terms", xtick=2a, ytick=-100:100)
    title!("log10 relative errors of Euler-Maclaurin approximaton of zeta($s)"; titlefontsize=12)
end

Out[12]:

plot_em_relerrs (generic function with 1 method)

In [13]:

plot_em_relerrs(4)

Out[13]:

In [14]:

plot_em_relerrs(3)

Out[14]:

In [15]:

plot_em_relerrs(2)

Out[15]:

In [16]:

plot_em_relerrs(1.5)

Out[16]:

In [17]:

plot_em_relerrs(0.5)

Out[17]:

In [18]:

plot_em_relerrs(0.001)

Out[18]:

In [ ]: