fortran66のブログ

fortran について書きます。

【メモ帳】オイラーの定数 γ

オイラーの定数 γ 速めの収束

irresistible integrals をパラパラ眺めていたら、§9.4.1 に \sum_{k=1}^n{1\over k}-\log(n+0.5) がやはりオイラー定数 γ に収束して、しかも収束が速いとありました。

Irresistible Integrals: Symbolics, Analysis and Experiments in the Evaluation of Integrals

Irresistible Integrals: Symbolics, Analysis and Experiments in the Evaluation of Integrals

本来の定義の \sum_{k=1}^n{1\over k}-\log(n) は収束が遅いので有名で、収束誤差が \mathcal{O}(1/2n) なのですが、\log に 0.5 足したものは \mathcal{O}(1/24n^2) の誤差で収束するとのことです。証明もそれほど面倒ではない模様です。

こんな、積分の中点公式みたいなもので加速されて出てくるなんて。お前はまだガンマを知らない、という気持ちです。

お前はまだグンマを知らない DVD

お前はまだグンマを知らない DVD

以下で、計算して確かめてみます。

実行結果

1000 まで足して、小数点以下 7 桁位まで求まっています。比較として本来の定義で計算した値も求めています。こちらは小数点以下 3 桁ほどまで求まっていて、理論通りの結果となっています。

文献値 0.57721 56649 (小数点以下 10 桁)

         100  0.577219790140490       0.582207331651529
         200  0.577216701374822       0.579713581573410
         300  0.577216126324240       0.578881405643301
         400  0.577215924668091       0.578465144068524
         500  0.577215831235245       0.578215331568328
         600  0.577215780449559       0.578048766753451
         700  0.577215749814171       0.577929780547829
         800  0.577215729924379       0.577840534693221
         900  0.577215716284744       0.577771117576444
        1000  0.577215706526555       0.577715581568206
続行するには何かキーを押してください . . .

プログラム

    program euler_const
        implicit none
        integer, parameter :: kd = kind(1.0d0)
        real(kd) :: g, s
        integer :: i
        s = 0.0_kd
        do i = 1, 10**3
            s = s + 1.0_kd / real(i, kind = kd)
            if (mod(i, 10**2) == 0) print *, i, s - log(i + 0.5_kd), s - log(i + 0.0_kd)
        end do
    end program euler_const