『岩波数学公式I』の例の数値積分をJuliaのQuadGKパッケージを使って精度保証する

00Iwanami.md

『岩波数学公式I』の例の数値積分をJuliaのQuadGKパッケージを使って精度保証する

2023年12月25日にJulia v1.10.0がリリースされました。おめでとうございます。

ちょうどJulia Advent Calendar 2023 https://qiita.com/advent-calendar/2023/julia に 積分のパッケージの紹介 https://zenn.dev/ohno/articles/440234fbb2adec がありました。

そのうちの一つ、QuadGKパッケージ https://juliamath.github.io/QuadGK.jl/stable/ を使って 『岩波数学公式I』のp.240の例の数値積分を精度保証しようとおもいます。 永らくゼロだとされ、1988年に訂正され、2018年にTwitterで少し話題になっていた積分です。 2016年に作ったRubyのRangeを積分範囲としてブロックで与えられた関数をシンプソンの公式で数値積分するライブラリを使って2022年に挑戦したのですが、 精度保証がなく、JITを使っても35秒もかかっていました。 https://gist.github.com/t-nissie/b6ef8d39229a2534498b

JuliaやQuadGKのインストールなどについてのメモ

Windowsなら適当なところ（たとえばDドライブ直下）にインストールして、PATHを通す。 Cygwinから使っている。

export PATH=/cygdrive/d/Julia-1.10.0/bin:$PATH

Macなら/Applications/Julia-1.10.app/Contents/Resources/julia/binにPATHを通す。 iTerm2から使っている。

PATHは ~/.bashrc などに設定する。juliaコマンドでREPLが起動する。REPLはexit()で抜ける。 ]でパッケージモード。会社の中などでプロキシを利用するならhttp_proxy, https_proxyを適宜設定しておく。

add QuadGK

でパッケージをインストール。バックスペースキーまたはdeleteキーでパッケージモードを抜ける。

ネイピア数＝自然対数の底の入力は\euler[TAB]。

QuadGKの精度保証はrtolとatol。

数値積分の戦略

戦略を練る

$$ \int_{0}^{1} \frac{1}{\log|\log x|} dx = -0.154479641320 $$

という定積分。曲者。被積分関数を $f(x) = \frac{1}{\log|\log x|}$ としてプロットすると、

となって、 $1/e$ で発散しているので直接数値積分はできない。 座標 $(1/e,0)$ の点のまわりで回転対称に近いので、 $x=1/e$ で $f(x)$ をひっくり返したもの $f(2/e-x)$ を足して、

$$ \int_{0}^{1/e} [f(x)+f(2/e-x)] dx + \int_{2/e}^{1} f(x) dx $$

がおそらく求める積分。 ただし $x=1/e$ 近傍では「大きな数−大きな数」になっているので、 区間 $(\frac{1}{e}-10^{-2}, \frac{1}{e}+10^{-2})$ では級数展開

$$f(x)+f(2/e-x) = -\frac{e^2}{12}(x-1/e)^2 - \frac{2e^4}{45}(x-1/e)^4 -\frac{859e^6}{30240}(x-1/e)^6 + O((x-1/e)^8)$$

のほうを数値積分する。級数展開は https://wolframalpha.com が教えてくれた。自分では導けていない。

Rubyでは気を遣った他のところはQuadGKなら自動でやってくれる。

数値積分を実行

いざ

julia> using QuadGK

julia> f = x -> 1.0/log(-log(x))
#1 (generic function with 1 method)

julia> g = x -> (abs(x-1/ℯ)<1e-2) ? -(x-1/ℯ)^2*ℯ^2/12 -2*(x-1/ℯ)^4*ℯ^4/45 -859*(x-1/ℯ)^6*ℯ^6/30240 : f(x)+f(2/ℯ-x)
#3 (generic function with 1 method)

julia> intf = quadgk(f, 2/ℯ, 1, atol=1e-15)
(-0.1363271202013918, 7.877473262818042e-16)

julia> intg = quadgk(g, 0, 1/ℯ, atol=1e-15)
(-0.018152521118649536, 7.668599698560629e-16)

julia> intf[1]+intg[1]
-0.15447964132004133

julia> intf[2]+intg[2]
1.554607296137867e-15

julia> exit()

まとめ

しくみはよくわからないのですがJuliaのQuadGKパッケージを使って例の数値積分を簡単に精度保証することができました。 Rubyでも同じようなことがしたいです。