fortran66のブログ

fortran について書きます。

ガンマ関数

Paul J. Nahin "Inside Interesting Integrals"

amazon AA adblock だと見えませんw

Paul J. Nahin "Inside Interesting Integrals" という本を眺めていると、\int_0^\infty e^{-y^2}dy, \int_0^\infty e^{-y^3}dy, \int_0^\infty e^{-y^4}dy, ・・・, \int_0^\infty e^{-y^n}dy という積分がでてきて、\int_0^\infty e^{-y^2}dy=\sqrt{\pi}/2は Gaussian の積分で常識として、2 乗以外は、3 乗も 4 乗も使われるのを見たこともないし、どうでもいいですわ~と思ったのですが、これを y^n=x の変数変換をすると、おなじみのガンマ関数が出てきてぎゃふんとなりました。

ny^{n-1}dy=dx よりdy={1\over n}x^{{1\over n}-1}dxなので、\int_0^\infty x^{{1\over n}-1}e^{-x}{1\over n}dx={1\over n}\Gamma({1\over n})=\Gamma({1\over n}+1) となります。

これはガンマ関数を提示するのに、なかなか面白い導入部ではないかと思いました。さて、e^{-y^n}は、n が 1 より大きい正の整数の時は急激に減少する関数なので、積分範囲の上限が +∞ でも、実際には小さな値まで取れば十分になっています。

数値計算

ここでは数値積分を計算して、それを Fortran2008 から組み込み関数になったガンマ関数と比較してみます。

ソース・プログラム

台形積分をリチャードソン補外することで、シンプソン法相当の計算をしています。はじめからシンプソン法をやればいいかもしれませんが、もう一段やればロンバーグ法相当になるので、まぁいいでしょうw

amazon AA

    program p_Gamma
      implicit none
      real, parameter :: pi = 4 * atan(1.0)
      integer, parameter :: n = 2**8  
      integer :: i
      real :: a, b, x(n), y(n), s, t, h, f
      !
      a = 0.0
      b = 15.0
      h = (b - a) / (n - 1)
      forall (i = 1:n) x(i) = a + h * (i - 1)
      !  
      do i = 1, 10
        f = real(i)
        y = exp(-x**f)
        s = richard(h, y) ! ~ simpson's method
        t = gamma(1.0 / f + 1) !gamma(1.0/i) / i
        print *, i, t, s, abs(t - s) / t  
      end do    
    contains
      real function trap(y) ! Trapezoidal rule
        real, intent(in) :: y(:)
        integer :: n
        n = size(y)
        trap = 0.5 *(y(1) + y(n)) + sum(y(2:n-1))
      end function trap
  
      real function richard(h, y) ! Richardson extrapolation
        real, intent(in) :: h, y(:)
        real :: s1, s2
        s1 =     h * trap(y)
        s2 = 2 * h * trap(y(::2))
        richard = (4 * s1 - s2) / 3 
      end function richard
    end program p_Gamma

実行結果

Γ(1/n+1) を\int_0^\infty e^{-y^n}dy積分を数値的に計算して、組み込み関数の値と比較しています。積分範囲の∞を 15 にしたとき、相対誤差は単精度の限界のあたりに来ています。

           1   1.000000      0.9999995      4.7683716E-07
           2  0.8862270      0.8862269      6.7256636E-08
           3  0.8929795      0.8929799      4.0048832E-07
           4  0.9064025      0.9064024      6.5759579E-08
           5  0.9181687      0.9181687      0.0000000E+00
           6  0.9277194      0.9277192      1.9274572E-07
           7  0.9354376      0.9354377      1.2743693E-07
           8  0.9417427      0.9417428      1.2658371E-07
           9  0.9469653      0.9469649      4.4059956E-07
          10  0.9513507      0.9513475      3.3832430E-06
続行するには何かキーを押してください . . .