03 π と e とEuler定数 γ
黒木玄2018-04-18～2019-04-03, 2023-06-22
Copyright 2018,2023 Gen Kuroki
License: MIT https://opensource.org/licenses/MIT (https://opensource.org/licenses/MIT)
Repository: https://github.com/genkuroki/Calculus (https://github.com/genkuroki/Calculus)
このファイルは次の場所できれいに閲覧できる:
http://nbviewer.jupyter.org/github/genkuroki/Calculus/blob/master/03%20pi%2C%20e%2C%20and%20Euler's%20γ.ipynb(http://nbviewer.jupyter.org/github/genkuroki/Calculus/blob/master/03%20pi%2C%20e%2C%20and%20Euler's%20%CE%B3.
https://genkuroki.github.io/documents/Calculus/03%20pi%2C%20e%2C%20and%20Euler's%20%CE%B3.pdf(https://genkuroki.github.io/documents/Calculus/03%20pi%2C%20e%2C%20and%20Euler's%20%CE%B3.pdf)
このファイルは Julia Box (https://juliabox.com) で利用できる.
自分のパソコンにJulia言語 (https://julialang.org/)をインストールしたい場合にはWindowsへのJulia言語のインストール (http://nbviewer.jupyter.org/gist/genkuroki/81de23edcae631a995e19a2ecf946a4f)
Julia v1.1.0 の Windows 8.1 へのインストール(https://nbviewer.jupyter.org/github/genkuroki/msfd28/blob/master/install.ipynb)
を参照せよ. 前者は古く, 後者の方が新しい.
論理的に完璧な説明をするつもりはない. 細部のいい加減な部分は自分で訂正・修正せよ.

目次
1 準備: 曲線の長さ
 2 円周率2.1 円周率の定義1
2.2 正弦函数の微積分
2.3 円周率の定義2
2.4 正接函数の微積分
2.5 Machinの公式による円周率の計算
 3 Napier数
3.1 対数函数と指数函数
3.2 Napier数の定義
 4 Euler定数
4.1 Euler定数の定義
4.2 条件収束交代級数への応用
4.3 Euler定数がガンマ函数の無限積表示に出て来ること
4.4 Euler定数が高次元単体上の一様分布に関連して自然に出て来ること
4.5 ζ(s) - 1/(s-1) の s→0 での極限がEuler定数に等しいこと

In [1]:



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

using Base.MathConstants
using Base64
using Printf
using Statistics
const e = ℯ
endof(a) = lastindex(a)linspace(start, stop, length) = range(start, stop, length=length)​
using Plots#gr(); ENV["PLOTS_TEST"] = "true"#clibrary(:colorcet)#clibrary(:misc)
default(fmt=:png)​function pngplot(P...; kwargs...)sleep(0.1)
pngfile = tempname() * ".png"
savefig(plot(P...; kwargs...), pngfile)showimg("image/png", pngfile)end
pngplot(; kwargs...) = pngplot(plot!(; kwargs...))​
showimg(mime, fn) = open(fn) do fbase64 = base64encode(f)display("text/html", """<img src="data:$mime;base64,$base64">""")
end​
using SymPy#sympy.init_printing(order="lex") # default#sympy.init_printing(order="rev-lex")​
using SpecialFunctions
using QuadGKIn [2]:



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

# Override the Base.show definition of SymPy.jl:# https://github.com/JuliaPy/SymPy.jl/blob/29c5bfd1d10ac53014fa7fef468bc8deccadc2fc/src/​@eval SymPy function Base.show(io::IO, ::MIME"text/latex", x::SymbolicObject)
print(io, as_markdown("\\displaystyle " * sympy.latex(x, mode="plain", fold_short_fra
end@eval SymPy function Base.show(io::IO, ::MIME"text/latex", x::AbstractArray{Sym})
function toeqnarray(x::Vector{Sym})a = join(["\\displaystyle " * sympy.latex(x[i]) for i in 1:length(x)], "\\\\")"""\\left[ \\begin{array}{r}$a\\end{array} \\right]"""end
function toeqnarray(x::AbstractArray{Sym,2})sz = size(x)a = join([join("\\displaystyle " .* map(sympy.latex, x[i,:]), "&") for i in 1:sz"\\left[ \\begin{array}{" * repeat("r",sz[2]) * "}" * a * "\\end{array}\\right]"
end
print(io, as_markdown(toeqnarray(x)))
end

1 準備: 曲線の長さ(𝑥(𝑡),𝑦(𝑡)), 𝑎 < 𝑡 < 𝑏 の長さ ℓ は定積分𝑏
′ (𝑡)2 + 𝑦′ (𝑡)2⎯ 𝑑𝑡ℓ = ∫𝑎 √⎯𝑥⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯
で定義される. 被積分函数 √ 𝑥′ (𝑡) 2 + 𝑦′ (𝑡) 2 は点 (𝑥(𝑡),𝑦(𝑡)) の運動の速さを意味する.
高校数学IIIの復習: 曲線の長さは速さの時間積分で定義される. すなわち, 曲線

2 円周率
2.1 円周率の定義1
円周率𝜋 の定義は単位半円周の長さである. 単位半円周 𝐶 は
𝐶 = {(𝑥,𝑦) ∈ ℝ2 ∣ 𝑥2 + 𝑦2 = 1, 𝑥 > 0}と定義される. これは

𝑥(𝑡) = √⎯1⎯⎯⎯⎯⎯⎯⎯− 𝑡2⎯, 𝑦(𝑡) = 𝑡
によって

𝐶 = {(𝑥(𝑡),𝑦(𝑡)) ∣ −1 < 𝑡 < 1}と表わされ, ちょっとした計算で⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯ 1√𝑥′ (𝑡)2 + 𝑦′ (𝑡)2 = √⎯1⎯⎯⎯⎯⎯⎯⎯− 𝑡2⎯となることがわかる(自分で計算してみよ). ゆえに, 円周率は

1 𝑑𝑡
𝜋 = ∫−1 ⎯⎯⎯⎯⎯⎯⎯⎯√1 − 𝑡2⎯と表わされる. これを円周率の定義だと思ってもよい.
In [3]:



1
2
3
4
5
6
7

xx(t) = sqrt(1-t^2)
yy(t) = tt = -1:0.01:1
plot(size=(500, 400), aspect_ratio=1, xlim=(-1.1,1.1), ylim=(-1.1,1.1))
plot!(xx.(t), yy.(t), label="C", legend=:topleft)
hline!([0], color="gray", ls=:dot, label="")vline!([0], color="gray", ls=:dot, label="")

1
2
3
4

t = symbols("t", real=true)x = sqrt(1-t^2)
y = t
simplify(diff(x,t)^2+diff(y,t)^2)

Out[3]:

In [4]:

Out[4]:

In [5]:

Out[5]:− 𝑡2 1− 1𝜋

1 integrate(1/sqrt(1-t^2), (t,-1,1))

2.2 正弦函数の微積分
前節の計算より,−1 ≦ 𝑦 ≦ 1 に対して,

𝑦 𝑑𝑡
𝐹(𝑦) = ∫0 √⎯1⎯⎯⎯⎯⎯⎯⎯− 𝑡2⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯
は 𝐶 上の点 (1,0) から (√ 1 − 𝑦2 ,𝑦) までの弧の長さの ±1 倍になる. 符号は 𝑦 > 0 なら + になり, 𝑦 < 0 なら − になる. すなわち, 𝜃(𝑦) は弧度法の意味での角度である. 単位円上の 𝑦 座標を角度 𝜃 の函数とみなしたものを 𝑦 = sin𝜃 と定義するのであった. これは, 𝑦 = sin𝜃 は 𝜃 = 𝐹(𝑦) の逆函数であることを意味する. 𝐹(𝑦) は arcsin𝑦 と書かれることが多い.

𝑑𝜃 = 𝐹 ′(𝑦) = 1 なので, その逆函数の導函数は
𝑑𝑦− 𝑦2⎯√⎯1⎯⎯⎯⎯⎯⎯⎯⎯𝑑sin𝜃 = 𝑑𝑦 = √⎯1⎯⎯⎯⎯⎯⎯⎯⎯
2⎯ = √⎯1⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯−− sin2 𝜃⎯.
𝑦𝑑𝜃 𝑑𝜃
これで sin𝜃 の導函数の計算の仕方もわかった.
2.3 円周率の定義2𝑡

傾き の直線

𝑦 = 𝑡𝑥 と単位半円周 𝐶 の交点を (𝑋,𝑌) = (𝑋(𝑡),𝑌(𝑡)) と書くと,
1 ⎯ , 𝑌 = 𝑌(𝑡) = ⎯⎯⎯⎯⎯⎯⎯⎯𝑡 ⎯.
𝑋 = 𝑋(𝑡) = √⎯1⎯⎯⎯⎯⎯⎯⎯
2+𝑡√1 + 𝑡2

このとき,⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯√𝑋 ′(𝑡)2 + 𝑌 ′ (𝑡)2 = 1 +1 𝑡2となる(自分で工夫して楽な計算の仕方を見付けるように努力してみよ). ゆえに円周率は𝜋 = ∫−∞ 1 +𝑑𝑡𝑡2∞と表わされる. これを円周率の定義だと思ってもよい.

問題: 積分∞ 𝑑𝑡
𝑦 と置換すると,∫−∞ 1 + 𝑡2 で 𝑡 = √⎯1⎯⎯⎯⎯⎯⎯⎯⎯
− 𝑦2⎯
∞ 𝑑𝑡
1 𝑑𝑦
=
∫−∞ 1 + 𝑡2 ∫−1 √⎯1⎯⎯⎯⎯⎯⎯⎯⎯− 𝑦2⎯□⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯
注意: 𝑥 = 1/√ 1 + 𝑡2 , 𝑦 = 𝑡/√ 1 + 𝑡2 のとき, 𝑥2 + 𝑦2 = 1 より 𝑥 = √ 1 − 𝑦2 となるので, 𝑡 = 𝑦/𝑥 = 𝑦/√ 1 − 𝑦2 . □.
𝑡 とおくと, 𝑑𝑦 = 𝑑𝑡 となることも確認せよ. □問題: 𝑦 =⎯⎯⎯⎯⎯⎯⎯⎯∫ √⎯1⎯⎯⎯⎯⎯⎯⎯⎯− 𝑦2⎯ ∫ 1 + 𝑡2
√1 + 𝑡2⎯
注意: 𝑚 > 0 であるとする. 曲線 𝑥𝑚 + 𝑦𝑚 = 1 と直線 𝑦 = 𝑡𝑥 の 𝑥 > 0 における交点は(𝑥,𝑦) = ( (1 + 1𝑡𝑚 )1/𝑚 , (1 + 𝑡𝑡𝑚 )1/𝑚 ) と表わされる. さらに, 𝑥𝑚 (1 + 𝑡𝑚 ) = 1 より 𝑚𝑥𝑚−1 (1 + 𝑡𝑚 )𝑑𝑥 + 𝑚𝑥𝑚 𝑡𝑚−1 𝑑𝑡 = 0と 𝑦 = 𝑡𝑥 より,
𝑑𝑥 = 𝑑𝑥 = 𝑑𝑥 = − 𝑥𝑚 𝑡𝑚−1 𝑑𝑡 = − 𝑡𝑚−2 𝑑𝑡 .
(1 − 𝑥𝑚 )1/𝑚 𝑦 𝑡𝑥 𝑡𝑥𝑚 (1 + 𝑡𝑚 ) 1 + 𝑡𝑚
上の円の場合は本質的に 𝑚 = 2 の特別な場合である. □が得られることを確認せよ.

In [6]:

Out[6]:

In [7]:

Out[7]:

1
2
3
4

1

t = symbols("t", real=true)
X = 1/sqrt(1+t^2)
Y = t*X
simplify(sqrt(diff(X,t)^2+diff(Y,t)^2))

𝑡2 + 1
1 y = symbols("y", real=true)
2 t = y/sqrt(1-y^2)
3 simplify(diff(t,y)/(1+t^2))

1 ⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯√1 − 𝑦2In [8]:Out[8]:

In [9]:

1 t = symbols("t", positive=true)
2 y = t/√(1+t^2)
3 simplify(diff(y,t)/√(1-y^2))1𝑡2 + 1
1
2
3
4
5



Out[9]:

t = symbols("t", positive=true)m = symbols("m", positive=true)x = 1/(1+t^m)^(1/Sym(m))
y = t*x
simplify(diff(x,t)/y)− 𝑡𝑚𝑡𝑚−2+ 1

以下の問題は https://twitter.com/cc226158/status/1043059177889783810(https://twitter.com/cc226158/status/1043059177889783810) より.
問題:

𝑑𝑥∫1 (9 − 𝑥3 )1/3 を求めよ.
2

1/3
𝑥3 + 𝑦3 = 9 と 𝑦 = 𝑡𝑥 の 𝑥 > 0 における交点の 𝑥 座標は 𝑥 = ( 1 +9 𝑡3 ) になる. これを逆に解くと
3 1/3
𝑡 = (9 −𝑥𝑥 ) なので 𝑥 = 1 のとき 𝑡 = 2 となり, 𝑥 = 2 のとき 𝑡 = 1/2 になる. 上の注意と同様の計算によって𝑑𝑥
𝑑𝑥 𝑡𝑑𝑡(9 − 𝑥3 )1/3 = 𝑦 = − 1 + 𝑡3 .
略解:ゆえに, 置換積分によって𝑑𝑥 = 2 𝑡𝑑𝑡 = 𝜋 ⎯⎯ ≈ 0.6046∫1 (9 − 𝑥3 )1/3 ∫1/2 1 + 𝑡3 3√3となることが若干の計算によって確かめられる. □In [10]:Out[10]:In [11]:Out[11]:

In [12]:

Out[12]:

In [13]:

Out[13]:2

1 x = symbols("x", positive=true)
2 simplify(integrate(1/(9-x^3)^(1/Sym(3)), (x, 1, 2)))

1 1|
1 1|
√3 ⎯⎯3 (2 𝐹1 ( 3 ,4 3 ||| 19 ) − 22 𝐹1 ( 3 ,4 3 ||| 89 ))
3 |
3 |
−
3
1 t = symbols("t", positive=true)
2 x = (9/(1+t^3))^(1/Sym(3))
3 diff(x, t)/(t*x)− 𝑡3 +𝑡 1

1 t = symbols("t", positive=true)
2 simplify(integrate(t/(1+t^3), t))⎯⎯ √3(2𝑡−1)
2 − 𝑡 + 1) √3 atan ( 3 )
log
𝑡
log(𝑡
+
1)
(
− 3 +
+
6
3
1 t = symbols("t", positive=true)
2 simplify(integrate(t/(1+t^3), (t, 1/Sym(2), 2)))

√⎯⎯3𝜋
9

2.4 正接函数の微積分
前節の計算より,𝑎𝜃 = 𝐺(𝑎) = ∫0 1 +𝑑𝑡𝑡2

(1,0)𝑎𝑦 = 𝑎𝑥は単位半円周上の点から傾き の直線と単位半円周等しくなる. すなわち, 弧度法の意味での角度に等しい.

𝐶 の交点までの弧の長さの ±1 倍になり, 符号 𝑎 の符号に

𝑎 = tan𝜃 が 𝜃 = 𝐺(𝑎) の逆函数になっていることを意味して𝐺(𝑎) arctan𝑎𝑑𝜃 = 𝐺′ (𝑎) = 1 なので, 𝑑tan𝜃 = 𝑑𝑎 = 1 + 𝑎2 = 1 + tan2 𝜃.
𝑑𝑎𝑑𝜃 𝑑𝜃
1 + 𝑎2正接函数の定義は傾きを角度の函数とみたものであった. これはいる.
はと書かれることが多い.

以下の公式もよく使われる:In [14]:𝑑𝑡 𝜋∫0 1 + 𝑡2 = 4 ,1

𝑑𝑡 𝜋∫0 1 + 𝑡2 = 2 .
∞1 t = symbols("t", real=true)
2 integrate(1/(1+t^2), (t,0,1)) |> display
3 integrate(1/(1+t^2), (t,0,oo)) |> display

𝜋
4
𝜋
2

2.5 Machinの公式による円周率の計算|𝑥| < 𝑥1 のとき 𝑥 ∞
∞
∞
𝑥
2𝑘+1
arctan𝑥 = ∫0 1 +𝑑𝑡𝑡2 = ∫0 ∑(−1)𝑘𝑡2𝑘 𝑑𝑡 = ∑(−1)𝑘 ∫0 𝑡2𝑘 𝑑𝑡 = ∑(−1)𝑘 2𝑘𝑥 + 1 .
𝑘=0
𝑘=0
𝑘=0
詳しくはノート「06 Taylorの定理」を参照せよ. □arctanのMaclaurin展開:

4arctan 15 − arctan 2391 = 𝜋4 .

Machinの公式:

証明:

tan の加法定理

cos𝛽 + cos𝛼 sin𝛽 tan𝛼 + tan𝛽tan(𝛼 + 𝛽) = sin𝛼
cos𝛼 cos𝛽 − sin𝛼 sin𝛽 = 1 − tan𝛼 tan𝛽より,

− tan𝛽2tan𝛼tan(𝛼 − 𝛽) = 1tan𝛼+ tan𝛼 tan𝛽 , tan(2𝛼) = 1 − tan2 𝛼 .

𝑎 = tan𝛼, 𝑏 = tan𝛽 とおくと, これらは次のように書き直される:
arctan𝑎 − arctan𝑏 = arctan 1𝑎+−𝑎𝑏𝑏 , 2arctan𝑎 = arctan 1 −2𝑎𝑎2 .
特に, 𝑎 = 𝑚/𝑛, 𝑏 = (𝑚 − 𝑛)/(𝑚 + 𝑛) のとき
𝑎 − 𝑏 = 𝑚(𝑚 + 𝑛) − 𝑛(𝑚 − 𝑛) = 𝑚2 + 𝑛2 = 1, 2𝑎 = 2𝑚𝑛
1 + 𝑎𝑏 𝑛(𝑚 + 𝑛) + 𝑚(𝑚 − 𝑛) 𝑚2 + 𝑛2
1 − 𝑎2 𝑛2 − 𝑚2
なので

ゆえに,

問題: 円周率arctan 𝑚𝑛 − arctan 𝑚𝑚 −+ 𝑛𝑛 = arctan1 = 𝜋4 , 2arctan 𝑚𝑛 = arctan 𝑛22𝑚𝑛− 𝑚2
4arctan 15 = 2arctan 125 = arctan 120
119 ,1 = arctan1 = 𝜋 .
4arctan 15 − arctan 2391 = arctan 120−
arctan119
239
4

𝜋 を小数点以下第3桁まで求めよ.

解答例: Machinの公式より,

𝜋 = 16arctan 15 − 4arctan 2391 .

□arctan𝑥 は |𝑥| < 1 で次のようにTaylor展開される:
∞
𝑠
2𝑘+1
3
5
7
arctan𝑥 = ∫0 1 +𝑑𝑡𝑡2 = ∑(−1)𝑘 2𝑘𝑥 + 1 = 𝑥 − 𝑥3 + 𝑥5 − 𝑥7 + ⋯𝑘=0

0 < 𝑥 < 1 のとき, 𝑥, 𝑥3 /3, 𝑥5 /5,… は正の実数の単調減少数列で 0 に収束するものなので,
𝑛
𝑛+1
2𝑘+1
2𝑘+1 𝑛−1
2𝑘+1
𝑛 is odd ⟹ ∑(−1)𝑘 2𝑘𝑥 + 1 < arctan𝑥 < ∑(−1)𝑘 2𝑘𝑥 + 1 < ∑(−1)𝑘 2𝑘𝑥 + 1 ,
𝑘=0
𝑘=0
𝑘=0
𝑛−1
𝑛+1
𝑛
2𝑘+1
2𝑘+1
2𝑘+1
𝑛 is even ⟹ ∑(−1)𝑘 2𝑘𝑥 + 1 < ∑(−1)𝑘 2𝑘𝑥 + 1 < arctan𝑥 < ∑(−1)𝑘 2𝑘𝑥 + 1 .
𝑘=0
𝑘=0
𝑘=0
ゆえに, 0 < 𝑥 < 1 のとき,
𝑛−1
2𝑛+1
2𝑘+1
𝑛 is odd ⟹ − 2𝑛𝑥 + 1 < arctan𝑥 − ∑(−1)𝑘 2𝑘𝑥 + 1 < 0,
𝑘=0
𝑛−1
2𝑘+1
2𝑛+1
𝑛 is even ⟹ 0 < arctan𝑥 − ∑(−1)𝑘 2𝑘𝑥 + 1 < 2𝑛𝑥 + 1 .
𝑘=0
𝑅(𝑥,𝑛) = (−1)𝑛 2𝑛𝑥+𝑛 1 とおくと,
16𝑅(1/2,3) = −0.000029⋯, 4𝑅(1/239,1) = −0.0000000097⋯
これらは小数点以下第2桁まで求めたい場合の誤差として十分に小さい. 𝑎 = 1/5 , 𝑏 = 1/239 のとき
3
5
16 (𝑎 − 𝑎3 + 𝑎5 ) = 3.158357333⋯,
4𝑏 = 0.01673640⋯⋯,3
5
16 (𝑎 − 𝑎3 + 𝑎5 ) − 4𝑏 = 3.14162⋯
したがって, Machinの公式と上の誤差の評価より, 𝜋 = 3.141⋯ であることがわかる. □In [15]:

1 # https://twitter.com/genkuroki/status/1024160463729225733
2 showimg("image/jpg", "machin.jpg")

In [16]:



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

R(x,n) = (-1)^n*x^(2n+1)/(2n+1)​@show a, b = 1/5, 1/239@show n, m = 3, 1
println()@show 16R(a,n)@show 4R(b,m)@show 16R(a,n)-4R(b,n)
println()@show 16a@show 16(a-a^3/3)@show 16(a-a^3/3+a^5/5)@show 4b@show 16(a-a^3/3+a^5/5)-4b
println()@show π - (16(a-a^3/3+a^5/5)-4b)@show float(π)(a, b) = (1 / 5, 1 / 239) = (0.2, 0.0041841004184100415)(n, m) = (3, 1) = (3, 1)
16 * R(a, n) = -2.925714285714287e-5
4 * R(b, m) = -9.766636714833518e-8
16 * R(a, n) - 4 * R(b, n) = -2.9257142857130044e-5
16a = 3.2
16 * (a - a ^ 3 / 3) = 3.1573333333333333
16 * ((a - a ^ 3 / 3) + a ^ 5 / 5) = 3.1583573333333335
4b = 0.016736401673640166
16 * ((a - a ^ 3 / 3) + a ^ 5 / 5) - 4b = 3.1416209316596935
π - (16 * ((a - a ^ 3 / 3) + a ^ 5 / 5) - 4b) = -2.8278069900355973e-5
float(π) = 3.141592653589793
Out[16]: 3.141592653589793
In [17]:



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

# 小数点以下第12桁まで正確に計算​
R(x,n) = (-1)^n*x^(2n+1)/(2n+1)
f(x,n) = sum(k->(-1)^k*x^(2k+1)/(2k+1), 0:n-1)​@show a, b = 1/5, 1/239@show n, m = 8, 2
println()@show 16R(a, n)@show 4R(b, m)@show 16R(a, n) - 4R(b, m)
println()@show 16f(a, n) - 4f(b, m)
println()@show π - (16f(a, n) - 4f(b, m))@show float(π);(a, b) = (1 / 5, 1 / 239) = (0.2, 0.0041841004184100415)(n, m) = (8, 2) = (8, 2)
16 * R(a, n) = 1.2336188235294129e-12
4 * R(b, m) = 1.025889257698589e-12
16 * R(a, n) - 4 * R(b, m) = 2.077295658308238e-13
16 * f(a, n) - 4 * f(b, m) = 3.1415926535896284
π - (16 * f(a, n) - 4 * f(b, m)) = 1.6475709685437323e-13
float(π) = 3.141592653589793
In [18]:

1 BigFloat(π)

Out[18]: 3.141592653589793238462643383279502884197169399375105820974944592307816406286198

3 Napier数
3.1 対数函数と指数函数
このノートでは𝑥 > 0 の(自然)対数函数 𝑦 = log𝑥 を

𝑥
𝑦 = log𝑥 = ∫1 𝑑𝑡𝑡で定義し, これの逆函数で指数函数

𝑥 = exp𝑦 を定義しておく.

log1 = 0, exp0 = 1 となる.
𝑑𝑦 = 𝑑log𝑥 = 1 が成立し, その逆函数については 𝑑𝑥 = 𝑑exp𝑦 = 𝑥 = exp𝑦.
さらに,𝑑𝑥 𝑑𝑥 𝑥𝑑𝑦 𝑑𝑦
これらの条件によって log𝑥 と exp𝑦 は 意的に定まることを示せるが ここではその証明には触れない
問題: 上の対数函数の定義だけを用いて log(𝑥1 𝑥2 ) = log 𝑥1 + log 𝑥2 (𝑥1 , 𝑥2 > 0) となることを示せ.
このとき,

証明:

2つ目の項の積分変数をlog(𝑥1 𝑥2 ) = ∫1𝑥1𝑥2 𝑑𝑡𝑥1𝑥2 𝑑𝑡𝑥2

𝑡 = 𝑥1 𝑠 で置換すると,

𝑑𝑡
𝑡 = ∫1 𝑡 + ∫𝑥𝑥1

𝑥1𝑥2 𝑑𝑡

1

𝑑𝑠
𝑥1 𝑑𝑠
𝑡 = ∫1 𝑥1 𝑠 = ∫1 𝑠 .
𝑥2

𝑡.

∫𝑥 1
𝑠 を改めて 𝑡 とおけば欲しい結果が得られる. □問題: 上の問題の結果を用いて exp(𝑦1 + 𝑦2 ) = exp(𝑦1 )exp(𝑦2 ) となることを示せ.
証明: 𝑥𝑖 = exp 𝑦𝑖 とおくと, 𝑦𝑖 = log 𝑥𝑖 であり, 上の問題の結果より, 𝑦1 + 𝑦2 = log(𝑥1 𝑥2 ) なのでexp(𝑦1 + 𝑦2 ) = explog(𝑥1 𝑥2 ) = 𝑥1 𝑥2 = exp(𝑦1 )exp(𝑦2 ). □𝑎 > 0 に対して, 𝑎𝑦 を
𝑎𝑦 = exp(𝑦log𝑎)によって定義する. このとき, 整数 𝑛 について
𝑎1 = exp(log𝑎)) = 𝑎, 𝑎𝑛 = exp(𝑛log𝑎)) = (𝑛 個の exp(log𝑎)) の積) = (𝑛 個の 𝑎 の積).
後者で上の問題の結果を使った. これより, 整数の 𝑛 について, 上の定義に基いた 𝑎𝑛 と 𝑛 個の 𝑎 の積が一致することがわかった. 上のように定義しても混乱は生じない.

3.2 Napier数の定義
このノートでは Napier数𝑒をによって定義する. すなわち𝑒は𝑒 = exp(1)𝑒 𝑑𝑡∫1 𝑡 = 1という条件で定義されているものとする.
定理:

証明:

lim (1 + 𝑦𝑛 )𝑛 = exp(𝑦).
𝑛→∞

𝑦 𝑛
𝑦
(1 + 𝑛 ) = exp[𝑛log(1 + 𝑛 )]である. 𝑑(log𝑥)/𝑑𝑥 = 1/𝑥 の 𝑥 = 1 の場合より, ℎ → 0 のとき,log(1 + ℎ) = log(1 + ℎ) − log(1) → 1.
ℎℎ
これを ℎ = 𝑦/𝑛 , 𝑛 → ∞ の場合に適用すると,

log(1 + 𝑦/𝑛) = 𝑦 log(1 + 𝑦/𝑛) → 𝑦.
1/𝑛
𝑦/𝑛

exp の連続性を用いると, 𝑛 → ∞ のとき,
𝑦 𝑛𝑦(1 + ) = exp[𝑛log(1 + )] → exp(𝑦).
上の定理の 𝑦 = 1 の特別な場合として次が得られる.
ゆえに,

系:

1 𝑛 = 𝑒.
1+lim𝑛→∞ ( 𝑛 )以上によって, 対数函数と指数函数とNapier数
In [19]:□□𝑒 の理論が微分積分学を用いてすべて再構成された.

1 BigFloat(e)Out[19]: 2.718281828459045235360287471352662497757247093699959574966967627724076630353555

4 Euler定数以上によって, 数学における基本的な定数

𝜋 と 𝑒 の定義が明瞭になった.

しかし, それだけで終わってしまうと, 高校の数学で習うことを単に再構成したに過ぎなくなってしまう.
そこでこのノートではその次に知っておくべきだと思われるEuler定数𝛾 について説明しよう.

Euler定数とその歴史については次の文献が詳しい:
https://arxiv.org/abs/1303.1856 (https://arxiv.org/abs/1303.1856)

4.1 Euler定数の定義
数列

𝑎𝑛 , 𝑏𝑛 を次のように定める:

𝑎𝑛 = 11 + 12 + ⋯ + 1𝑛 − log(𝑛 + 1)
𝑏𝑛 = 11 + 12 + ⋯ + 1𝑛 − log𝑛.

𝑦 = 1/𝑥 は狭義単調減少函数なので
𝑛+1𝑛log(𝑛 + 1) − log𝑛 = ∫𝑛 𝑑𝑥𝑥 < 1𝑛 < ∫𝑛−1 𝑑𝑥𝑥 = log𝑛 − log(𝑛 − 1),𝑛+1
log(𝑛 + 1) = ∫1 𝑑𝑥𝑥 < 11 + 12 + ⋯ + 1𝑛となる(図を描いて確認せよ). 後者より,

0 < 𝑎𝑛 < 𝑏𝑛が成立していることがわかり, 前者より,

𝑎𝑛 − 𝑎𝑛−1 = 1𝑛 + log𝑛 − log(𝑛 + 1) > 0,𝑏𝑛 − 𝑏𝑛−1 = 1𝑛 + log(𝑛 − 1) − log𝑛 < 0.となり, 𝑎𝑛 は単調増加し, 𝑏𝑛 は単調減少することがわかる. (図を描けば, 𝑎𝑛 が単調増加し, 𝑏𝑛 が単調減少することは見た目的に明らかになる.)𝑎𝑛 も 𝑏𝑛 も収束することがわかる.
さらに, 0 < log(𝑛 + 1) − log𝑛 < 1/𝑛 より 𝑎𝑛 と 𝑏𝑛 の収束先が同じであることもわかる.
それらの収束先を 𝛾 と書き, Euler定数と呼ぶ.
𝛾 = 𝑛→∞lim ( 11 + 12 + ⋯ + 1𝑛 − log(𝑛 + 1)) = 𝑛→∞lim ( 11 + 12 + ⋯ + 1𝑛 − log𝑛) .
したがって,

この定義は次のように書き換えられる:1 + 1 + ⋯ + 1 = log𝑛 + 𝛾 + 𝑜(1).
1 2
𝑛

ここで

𝑜(1) は 𝑛 → ∞ で 0 に収束する項を表す.

考え方: 解析学では(調べたい量) = (相対的にわかりやすい量) + (誤差)の形式で考える. Euler定数の定義は, 調べたいもの1 + 1 +⋯+ 1
1 2𝑛を相対的にわかりやすい量𝑛log𝑛 + 𝛾

log𝑛 と比較して第2項の定数項 𝛾 は小さいとみ
1/1 + 1/2 + ⋯ + 1/𝑛 を log𝑛 で大雑把に近似するより, log𝑛 + 𝛾

プラス誤差の形式で認識することであるとみなせる. が大きなとき, 第1項のなせ, 第2項 と比較して誤差項は小さいとみなせる.
で近似する方がより精密な近似になる.

𝛾

𝑜(1)

このように新たに項を付け加えることによって近似の精度を逐次的に高めて行くことは, 解析学における典型的なパターンであり, あらゆる場所で出会うことになる.

□

注意: digamma函数

𝜓(𝑠) (𝑠 > 0) を

𝜓(𝑠) = − 𝑛→∞lim [ 1𝑠 + 𝑠 +1 1 + ⋯ + 𝑠 +1 𝑛 − log𝑛]と定めることができる. これはガンマ函数がと書けることを使えば

𝑛𝑠 𝑛!
Γ(𝑠) = 𝑛→∞lim 𝑠(𝑠 + 1)⋯(𝑠+ 𝑛)𝜓(𝑠) = 𝑑logΓ(𝑠)/𝑑𝑠 とも表わされる. ガンマ函数に関するこの公式についてはガンマ分布の中心極限定理とStirlingの公式 (https://genkuroki.github.io/documents/20160501StirlingFormula.pdf)

𝜓(1) = −𝛾 であることがわかり, digamma函数の定義より 𝜓(𝑠 + 1) = 1/𝑠 + 𝜓(𝑠)𝑛 について𝜓(𝑛 + 1) = 11 + 12 + ⋯ + 1𝑛 − 𝛾であることがわかる. この公式はdigamma函数 𝜓(𝑠 + 1) は調和級数からEuler定数を引いたものの連続的な補間になっていることを意味している. したがって, 𝜓(𝑠 + 1) の 𝑠 → ∞ での漸近挙動がよくわかれば, 調和級数の漸近挙動もよくわかることになの第8.2節を参照せよ. Euler定数の定義よりであることがわかるので, 帰納的に非負の整数る.

𝑠 → ∞ においてlogΓ(𝑠 + 1) = 𝑠log𝑠 − 𝑠 + 12 log𝑠 + log √⎯2𝜋⎯⎯⎯ + 12𝑠1 + 𝑂 ( 𝑠12 ) .

ガンマ函数については次のStirlingの近似公式がよく知られている:

これを𝑠 で微分することによって,𝜓(𝑠 + 1) = log𝑠 + 2𝑠1 − 12𝑠1 2 + 𝑂 ( 𝑠13 )以上をまとめると次が得られる:(𝑠 → ∞).

1 + 1 + ⋯ + 1 = log𝑛 + 𝛾 + 1 − 1 + 𝑂 1
2𝑛 12𝑛2 ( 𝑛3 ) (𝑠 → ∞).
1 2
𝑛
1 + 1 + ⋯ + 1 の 𝑛 → ∞ における漸近挙動の詳細はガンマ函数に関するStirlingの公式の対数微分にこのように調和級数
1 2𝑛よって得られる. □In [20]:



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

function plot_eulergamma(N, n)f(x) = 1/xp = plot(xlim=(0,N+2), ylim=(0, 1.2))
x = 1.0:0.05:n
plot!(p, x, f.(x), color=:blue, label="log($n)", fill=(0, 1.0, :pink))
x = n:0.05:n+1
plot!(p, x, f.(x), color=:blue, label="log($(n+1))", fill=(0, 0.5, :pink))
plot!(p, [1,1], [0,1], color=:red, label="1/1+1/2+...+1/$n")
for k in 1:N
x = k:0.05:k+1
g(x) = 1/kc = ifelse(k > n, :pink, :red)
plot!(p, x, g.(x), color=c, label="")plot!(p, [k+1,k+1], [0,1/k], color=c, label="")end
x = 0.1:0.05:N+2
plot!(p, x, f.(x), color=:blue, label="y=1/x")
plot!(title="n = $n", titlefontsize=10)
p
end
​
N = 6
p1 = plot_eulergamma(N, N-1)
p2 = plot_eulergamma(N, N)
plot(p1, p2, size=(700, 300))

Out[20]:In [21]:

1 aa(n) = sum(k->1/k, 1:n) - log(n+1)
2 bb(n) = sum(k->1/k, 1:n) - log(n)3 [(n, aa(n), bb(n)) for n in [1:9; @.(10^(1:7))]]

Out[21]: 16-element Vector{Tuple{Int64, Float64, Float64}}:(1, 0.3068528194400547, 1.0)
(2, 0.4013877113318902, 0.8068528194400547)(3, 0.4470389722134427, 0.7347210446652235)(4, 0.47389542089923276, 0.6970389722134425)(5, 0.49157386410527826, 0.6738954208992329)(6, 0.5040898509446865, 0.6582405307719448)(7, 0.5134156011773068, 0.6469469938018293)(8, 0.520632565520923, 0.6384156011773068)(9, 0.5263831609742078, 0.6317436766320341)(10, 0.5310729811698831, 0.6263831609742079)(100, 0.5722570007983601, 0.5822073316515279)(1000, 0.5767160812351229, 0.5777155815682065)(10000, 0.5771656690678633, 0.5772656640681966)(100000, 0.5772106649431983, 0.5772206648931988)(1000000, 0.5772151649019452, 0.5772161649014453)(10000000, 0.5772156149015366, 0.5772157149015307)
In [22]:

1 eulergamma

Out[22]: γ = 0.5772156649015...

4.2 条件収束交代級数への応用
次の2つの数列を考える:

𝐴𝑛 = 11 + 13 + ⋯ + 2𝑛 1− 1 , 𝐵𝑛 = 12 + 14 + ⋯ + 2𝑛1 .

このとき

2𝑛
𝑛
𝐴𝑛 + 𝐵𝑛 = ∑ 𝑘1 = log(2𝑛) + 𝛾 + 𝑜(1), 2𝐵𝑛 = ∑ 𝑘1 = log𝑛 + 𝛾 + 𝑜(1).
𝑘=1

𝑘=1

𝑜(1) は 𝑛 → ∞ で 0 に収束する項を意味する. ゆえに,
𝐴𝑛 = log2 + 12 log𝑛 + 𝛾2 + 𝑜(1), 𝐵𝑛 = 12 log𝑛 + 𝛾2 + 𝑜(1).
これより, 𝑛 → ∞ で𝐴𝑛 − 𝐵𝑛 → log2.
ここで

すなわち,

1 − 1 + 1 − 1 + ⋯ = log2.
1 2 3 4𝑝,𝑞 について𝐴𝑝𝑛 = log2 + 12 (log𝑝 + log𝑛) + 𝛾2 + 𝑜(1), 𝐵𝑞𝑛 = 12 (log𝑞 + log𝑛) + 𝛾2 + 𝑜(1)より一般に正の整数

なので,

例:

𝑛 → ∞ のとき

𝐴𝑝𝑛 − 𝐵𝑞𝑛 → log2 + 12 log 𝑝𝑞 .

𝑝 = 2, 𝑞 = 1 のとき

1 + 1 − 1 + 1 + 1 − 1 + ⋯ = 3 log2.
1 3 2 5 7 42

𝑝 = 1, 𝑞 = 2 のとき

1 − 1 − 1 + 1 − 1 − 1 + ⋯ = 1 log2.
1 2 4 3 6 8
2

考え方: 以上のような計算が可能になったのは,が大雑把に

1 + 1 +⋯+ 1
1 2
𝑛

log𝑛 で近似できることだけではなく, より精密に
1 + 1 + ⋯ + 1 = log𝑛 + 𝛾 + 𝑜(1)1 2𝑛の形で漸近挙動がわかっているからである.
In [23]:



1
2
3
4
5
6
7

A(n) = sum(k->1/(2k-1), 1:n)
B(n) = sum(k->1/(2k), 1:n)​n = 10^7@show S = A(2n) - B(n)@show T = 3/2*log(2)@show S - T;

S = A(2n) - B(n) = 1.0397207458399187
T = (3 / 2) * log(2) = 1.0397207708399179
S - T = -2.499999918192941e-8
In [24]:

1 sum(k->1/(4k-3)+1/(4k-1)-1/(2k), 1:10^7) - T

Out[24]: -2.500000007010783e-8

4.3 Euler定数がガンマ函数の無限積表示に出て来ること𝑠 > 0 のとき, ガンマ函数 Γ(𝑠) がによって定義され, digamma函数

𝜓(𝑠) が∞Γ(𝑠) = ∫0 𝑒−𝑥 𝑥𝑠−1 𝑑𝑥□′
1 ∞ 𝑒−𝑥 𝑥𝑠−1 log𝑥𝑑𝑥
𝜓(𝑠) = 𝑑𝑠𝑑 logΓ(𝑠) = ΓΓ(𝑠)(𝑠) = Γ(𝑠)∫0によって定義される.
ガンマ函数について次が成立してることが知られている(Weierstrassの公式):

1 = 𝑒𝛾𝑠 𝑠 ∞ 1 + 𝑠 𝑒−𝑠/𝑛 .
[( 𝑛 ) ]
∏
Γ(𝑠)
𝑛=1

この公式の証明についてはガンマ分布の中心極限定理とStirlingの公式 (https://genkuroki.github.io/documents/20160501StirlingFormula.pdf)の第8.2節を参照せよ.
Weierstrassの公式より∞−𝜓(𝑠) = 𝛾 + 1𝑠 − ∑ [ 1𝑛 − 𝑛 +1 𝑠 ] .
𝑛=1

特に

∞

−𝜓(1) = − ∫0 𝑒−𝑥 log𝑥𝑑𝑥 = 𝛾 = 0.5772156649015⋯.

Euler定数𝛾 は次の積分表示を持つことも示せる:∞
𝛾 = ∫0 ( 1 −𝑒−𝑡𝑒−𝑡 − 𝑒𝑡−𝑡 ) 𝑑𝑡.

この公式の証明については上の文献の第8.7節を参照.
In [25]:

1 -digamma(big(1))

Out[25]: 0.5772156649015328606065120900824024310421593359399235988057672348848677267776685

4.4 Euler定数が高次元単体上の一様分布に関連して自然に出て来ること

𝑟 > 0 に対して, 𝑛 − 1 次元単体 Δ𝑛−1 (𝑟) がΔ𝑛−1 (𝑟) = {(𝑥1 ,…, 𝑥𝑛 ) ∈ ℝ𝑛≧0 ∣ 𝑥1 + ⋯ + 𝑥𝑛 = 𝑟}と定義される. ここで ℝ ≧0 は非負の実数全体の集合を表す. 例えば, 𝑛 = 3 のとき
Δ2 (𝑟) = {(𝑥,𝑦,𝑧) ∈ ℝ3 ∣ 𝑥,𝑦,𝑧 ≧ 0, 𝑥 + 𝑦 + 𝑧 = 𝑟}.
は点 (𝑟,0,0), (0,𝑟,0), (0,0,𝑟) を頂点とする正三角形になる.
Δ𝑛−1 (𝑛) 上の一様分布(例えば 𝑛 = 3 の場合には正三角形上に一様に分布する確率分布)を考える.
そのとき, Δ𝑛−1 (𝑛) の点の第 𝑖 座標 𝑥𝑖 の分布(例えば 𝑥1 の分布)は 𝑛 → ∞ のとき確率密度函数が exp(−𝑥) の指数分布に近付くことを証明できる.

𝑛 → ∞ で log𝑥𝑖 の期待値は ∫0∞ 𝑒−𝑥 log𝑥𝑑𝑥 = 𝜓(1) = −𝛾 = −0.5772⋯ に近付く.
このような形でも自然にEuler定数 𝛾 が自然に出て来る.
ゆえに,

このような現象は「統計力学」ではよく知られている. 詳しくは次のノートを参照せよ:
Kullback-Leibler情報量とSanovの定理 (https://genkuroki.github.io/documents/20160616KullbackLeibler.pdf)

問題: 次のセルの rand_Uniform_on_Simplex 函数がどうしてなっているかを説明せよ.

(𝑥1 ,…, 𝑥𝑛 ) (𝑥1 , 𝑥1 + 𝑥𝑛−1
2 ,…, 𝑥1 + ⋯ + 𝑥𝑛−1 )
𝑛−1
[0,1]
(𝑛 − 1)!

Δ𝑛−1 (1) 上の一様分布に従う擬似乱数を発生させる函数に
Δ𝑛−1 (1)
𝑛−1
□

ヒント:
を
に対応させる写像を考えよ. その写像を通して, 単体
上の
一様分布と
次元立方体
上の一様分布をどのように関係付けられるかについて考えよ.
次元立方体はその成分の並び方の順序で分類することによって個の部分集合に分割される(境界の重なりは無視してよい).

In [26]:



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

# r = 1 の場合の n-1 次元単体上の一様分布に従う擬似乱数を発生させる函数
function rand_Uniform_on_Simplex(n)Z = sort(rand(n-1))
return vcat(Z, 1.0) - vcat(0.0, Z)
end​
function rand_Uniform_on_Simplex(n, N)
X = zeros(n, N)
for j in 1:N
X[:,j] = rand_Uniform_on_Simplex(n)endreturn X
end​# x_1 の分布と log x_i 達の平均の分布をプロットする函数
function plot_Uniform_on_Simplex(n; N=10^4)
X = n*rand_Uniform_on_Simplex(n, N)
ElogX = vec(mean(log.(X), dims=1))​
p1 = plot(xlabel="x", xlim=(0, 7.5), title="n = $n")
histogram!(p1, X[1,:], normed=true, alpha=0.5, label="x_1")
plot!(p1, x->exp(-x), 0, 7.5, lw=2, label="exp(-x)")​p2 = plot(xlabel="x", legend=:topleft, legendfontsize=7)
xlims!(-eulergamma-0.3, -eulergamma+0.3)
histogram!(p2, ElogX, normed=true, alpha=0.5, linealpha=0.1, label="mean of log x_i"
vline!([-eulergamma], lw=2, label="x = $(round(Float64(-eulergamma),digits=5))")​plot(p1, p2, size=(700, 300), layout=grid(1,2,widths=[0.35, 0.65]))end

Out[26]: plot_Uniform_on_Simplex (generic function with 1 method)
In [27]:



1 plot_Uniform_on_Simplex(100)Out[27]:

In [28]:

Out[28]:

1 plot_Uniform_on_Simplex(1000)

In [29]:

1 @time plot_Uniform_on_Simplex(10000)
4.683882 seconds (450.77 k allocations: 6.794 GiB, 6.02% gc time)Out[29]:𝑛−1

Δ𝑛−1 (𝑛)

𝑥1

𝑛→∞

問題:
次元単体
上の一様分布に関するの分布がで確率密度函数がとを示せ. 直観的に明瞭な説明をしていれば必ずしも論理的に厳密な議論でなくてもよい.

exp(−𝑥) の指数分布に近付くこ

Δ𝑛−1 (𝑟) を ℝ𝑛−1 = {(𝑥1 ,…, 𝑥𝑛−1 )} への射影の像を 𝐷𝑛−1 (𝑟) と書くことにする:
𝐷𝑛−1 (𝑟) = {(𝑥1 ,…, 𝑥𝑛−1 ) ∈ ℝ𝑛−1≧0 ∣ 𝑥1 + ⋯ + 𝑥𝑛−1 ≦ 𝑟}.
Δ𝑛−1 (𝑟) 上の一様分布を 𝐷𝑛−1 (𝑟) に射影して得られる分布は 𝐷𝑛−1 (𝑟) 上の一様分布になる.
(𝑟) の体積は 𝑟𝑛−1 /(𝑛 − 1)! に等しい. そして,0 ≦ 𝑎 ≦ 𝑟 のとき, 𝐷𝑛−1 (𝑟) を 𝑥1 ≧ 𝑎 に制限して得られる部分集合は
𝐷𝑛−1
𝑛−1
𝑛−1
𝐷 (𝑟 − 𝑎) に合同なのでその体積は (𝑟 − 𝑎) /(𝑛 − 1)! に等しい. ゆえに, 𝑥1 ≧ 𝑎 となる確率は 𝑃(𝑥1 ≧ 𝑎) は/(𝑛 − 1)! = 1 − 𝑎 𝑛−1
𝑃(𝑥1 ≧ 𝑎) = (𝑟 −𝑟𝑛−1𝑎)𝑛−1/(𝑛 − 1)! ( 𝑟 )となるので, 𝑟 = 𝑛 のとき, 𝑛 → ∞ で𝑃(𝑥1 ≧ 𝑎) = (1 − 𝑎𝑛 )𝑛−1 → exp(−𝑎).
一方, 確率密度函数 exp(−𝑥) において 𝑥 ≧ 𝑎 となる確率は∞∫𝑎 exp(−𝑥)𝑑𝑥 = [−exp(−𝑥)]∞𝑎 = exp(−𝑎).
以上を比較することによって, 𝐷𝑛−1 (𝑛) 上の一様分布に関する 𝑥1 の分布は 𝑛 → ∞ で確率密度函数 exp(−𝑥) の指数分布に近付くことがわかる. □注意: 以上の議論において 𝐷𝑛−1 (𝑟) に射影したことは本質的ではなく, Δ𝑛−1 (𝑟) の体積(面積)が 𝑟𝑛−1 に比例することだけを使っても同様の議論を遂行できる. 𝐷𝑛−1 (𝑟) に射影するメリットは比例定数が 1/(𝑛 − 1)! になることが易しく求まることに過ぎない. その比例定数は確率を求めるときに分子分母でキャンセルして消えるので議論の本質には何も影響しない. □問題: 𝐷1 (𝑟), 𝐷2 (𝑟), 𝐷3 (𝑟) を図示し, それぞれの長さ, 面積, 体積がそれぞれ 𝑟 , 𝑟2 /2 , 𝑟3 /6 であることを確認せよ. □解答例: 簡単のため,

In [30]:



1 plot(size=(300,300), aspect_ratio=1)
2 plot!([0,1], [0,0], lw=3, color=:blue, label="D^1(1)")
3 plot!(xlims=(-0.1,1.1), ylims=(-0.5,0.5))

Out[30]:In [31]:

Out[31]:

1
2
3
4
5

plot(size=(300,300), aspect_ratio=1)
plot!([0,1], [1,0], color=:blue, fill=(0, 0.5), label="D^2(1)")
plot!([0,0], [0,1], color=:blue, label="")plot!([0,1], [0,0], color=:blue, label="")plot!(xlims=(-0.1,1.1), ylims=(-0.1,1.1))

In [32]:

1
2
3
4
5
6
7
8
9



x = 0:0.01:1
y = 0:0.01:1
f(x,y) = ifelse(x+y>1, NaN, 1-x-y)​
pyplot()p = plot(size=(350, 300), zlims=(0,1))
surface!(x, y, f, color=:rainbow, camera=(120, 45))
display(p)gr();

sys:1: UserWarning: Z contains NaN values. This may result in rendering artifacts.

4.5 ζ(s) - 1/(s-1) の s→0 での極限がEuler定数に等しいこと
以下ではこの先に進まないと知らないはずのことを勝手に用いる.
問題:ℎ > 0 のときとなることを示せ.

1 − 1 = 1.
∑(𝑛=1 𝑛ℎ (𝑛 + 1)ℎ )
∞

□

解答略.
問題:

ℎ > 0 ならば, 𝑛 → ∞ で

1 1 1
1
1
𝑛1+ℎ − ℎ ( 𝑛ℎ + (𝑛 + 1)ℎ ) = 𝑂 ( 𝑛2+ℎ )となることを示せ.
解答例: 二項展開より

1 −ℎ = 1 − ℎ + ℎ(ℎ + 1) + ⋯ = 1 − ℎ + 𝑂 1
1
+
( 𝑛)
𝑛 2𝑛2
𝑛 ( 𝑛2 )なので−ℎ(左辺) = ℎ𝑛1ℎ ( ℎ𝑛 − 1 + (1 + 1𝑛 ) )1 .
= ℎ𝑛1ℎ ( ℎ𝑛 − 1 + 1 − ℎ𝑛 + 𝑂 ( 𝑛12 )) = 𝑂 ( 𝑛2+ℎ)
定義: 𝑠 > 1 に対して, 𝜁(𝑠) を∞𝜁(𝑠) = ∑ 𝑛1𝑠□𝑛=1と定める. この級数は絶対収束している.
問題:ℎ > 0 ならば∞∞ (log𝑛)𝑘 (log𝑛)𝑘+1 − (log(𝑛 + 1))𝑘+1𝑘𝜁(1 + ℎ) − 1ℎ = ∑ (−ℎ)
𝛾
,
𝛾
=
𝑘
𝑘
∑( 𝑛 +)𝑘!
𝑘+1となることを示せ.

□

𝑘=0

𝑛=1

解答略.

𝑛 に関する総和に等しい.
1/𝑛ℎ = exp(−ℎlog𝑛) と∞ 𝑥𝑘
1/(𝑛 + 1)ℎ = exp(−ℎlog(𝑛 + 1)) に exp(𝑥) = ∑ 𝑘! を適用し, −ℎ に関する同次の項をまとめて整理せよ. そのとき 1/𝑛1+ℎ
𝑘=1
は exp(−ℎlog𝑛)/𝑛 だとみなせ. □問題: 上の問題の 𝛾0 がEulerの定数 𝛾 に等しいことを示せ.
ヒント: この問題の左辺は上の問題の左辺の解答例:∞
𝛾0 = ∑ ( 1𝑛 + log𝑛 − log(𝑛 + 1))
𝑛=1= ( 11 − log2) + ( 12 + log2 − log3) + ( 13 + log3 − log4) + ⋯であり, この級数の 𝑛 = 𝑁 までの和は
1 + 1 + ⋯ + 1 − log(𝑁 + 1)1 2𝑁に等しいので, その級数はEuler定数 𝛾 に収束する. □以上によって𝜁(1 + ℎ) = 1ℎ + 𝛾 + 𝑂(ℎ) (ℎ → 0)であることが示された.
問題:

𝜁(1 + ℎ) = 1ℎ + 𝛾 + 𝑂(ℎ) (ℎ → 0) を以上とは別の方法を使って示してみよ. □例えばEuler-Maclaurinの和公式を使っても証明可能である.
問題:ℎ ↘ 0 のとき, 𝜁(1 + ℎ) − 1ℎ → 𝛾 となることを数値的に確認せよ. □In [33]:



1
2
3
4
5
6
7

N = 25@show h = big"0.1"^N@show c = zeta(1+h)-1/h@show g = big(eulergamma)@show round(c, digits=N)@show round(g, digits=N)@show c - g;

h = big"0.1" ^ N = 1.0000000000000000000000000000000000000000000000000000000000000000000000
00000053e-25
c = zeta(1 + h) - 1 / h = 0.577215664901532860606512097712841345076300974995098998824212132
3324325714408041
g = big(eulergamma) = 0.5772156649015328606065120900824024310421593359399235988057672348848
677267776685
round(c, digits = N) = 0.577215664901532860606512100000000000000000000000000000000000000000
0000000000015
round(g, digits = N) = 0.577215664901532860606512100000000000000000000000000000000000000000
0000000000015
c - g = 7.630438914034141639055175400018444897447564844663135574900855315781988976401014e-2
7
In [34]:



1 BigFloat(eulergamma)

Out[34]: 0.5772156649015328606065120900824024310421593359399235988057672348848677267776685