fortran66のブログ

fortran について書きます。

【寝言】夏至

至点 


www.youtube.com

Euler の γ

 \displaystyle \gamma := \lim_{n\rightarrow\infty}\left(\sum_{i=1}^n{1\over i} - \log(n)\right)

収束が遅い(~O(1/n) ) Euler の γ ですが、引くlog(n) を定義からずらして log(n+1/2) にすると 定義より早く O(1/24n2) で収束するのが草。

計算による比較

小さい項からの和

    program euler_gamma
        implicit none
        integer, parameter :: kd = kind(0.0d0)
        real(kd) :: s
        integer :: i, n = 10**7
        s = 0.0_kd
        do i = n, 1, -1
            s = s + 1.0_kd / i
        end do
        print *, s - log(real(n)), s - log(n + 0.5_kd), 1.0_kd / real(n)**2 / 24
   end  program euler_gamma    

カスを集めるのに小さい方から足しました。

定義通り、log(n+1/2)、誤差見積もり

  0.577215967910643       0.577215664901544       4.166666650965334E-016

組み手和(Kahan の和)

   block      
        real(kd) :: s, t, u, v 
        integer :: i, n = 10**7
        s = 0.0_kd
        v = 0.0_kd
        do i = 1, n
            u = 1.0_kd / i - v
            t = s + u
            v = (t - s) - u
            s = t
        end do
        print *, s - log(real(n)), s - log(n + 0.5_kd), 1.0_kd / real(n)**2 / 24
    end block    

カス項が効いているようなので、もう少し真面目な総和計算するとより良い値が得られました。

  0.577215967910632       0.577215664901534       4.166666650965334E-016

Wolfram language

N[EulerGamma, 20]
0.57721566490153286061

もう少し収束の速いもの

 \displaystyle \gamma = {1\over2}+\sum_{n=1}^{\infty}\sum_{m=2^{n-1}}^{2^n-1}{1\over2m(2m+1)(2m+2)}

    program euler_gamma
        implicit none
        integer, parameter :: kd = kind(0.0d0)
        integer :: i, j, n, m, m0, m1, nmax
        real(kd) :: s, t, x2
        nmax = 27
        s = 0.5_kd
        m1 = 1
        do n = 1, nmax
            t = 0.0_kd
            m0 = m1
            m1 = 2**n
            do m = m0, m1 - 1 
                x2 = 2 * m
                t = t + n / (x2 * (x2 + 1) * (x2 + 2))   
            end do 
            s = s + t
            print *, n, m1, s, t
            if (t < 1.0e-16_kd) exit
        end do    
        print *, 'gamma =', s
    end program euler_gamma
           1           2  0.541666666666667       4.166666666666666E-002
           2           4  0.564285714285714       2.261904761904762E-002
           3           8  0.572991591741592       8.705877455877456E-003
           4          16  0.575914184398886       2.922592657293993E-003
           5          32  0.576829154090275       9.149696913897526E-004
           6          64  0.577103770405349       2.746163150734348E-004
           7         128  0.577183875992416       8.010558706751096E-005
           8         256  0.577206763957738       2.288796532201553E-005
           9         512  0.577213201244027       6.437286288142040E-006
          10        1024  0.577214989382304       1.788138277449749E-006
          11        2048  0.577215481120550       4.917382461222898E-007
          12        4096  0.577215615230996       1.341104457486260E-007
          13        8192  0.577215651552576       3.632158007173262E-008
          14       16384  0.577215661331463       9.778887010694266E-009
          15       32768  0.577215663950808       2.619344739581377E-009
          16       65536  0.577215664649300       6.984919308599662E-010
          17      131072  0.577215664834837       1.855369191549276E-010
          18      262144  0.577215664883949       4.911271389529140E-011
          19      524288  0.577215664896910       1.296029950023478E-011
          20     1048576  0.577215664900320       3.410605131646578E-012
          21     2097152  0.577215664901215       8.952838470575934E-013
          22     4194304  0.577215664901450       2.344791028008255E-013
          23     8388608  0.577215664901511       6.128431095930439E-014
          24    16777216  0.577215664901527       1.598721155460254E-014
          25    33554432  0.577215664901531       4.163336342344230E-015
          26    67108864  0.577215664901532       1.082467449009531E-015
          27   134217728  0.577215664901533       2.810252031082286E-016
 gamma =  0.577215664901533

思ったより収束が早くないです。