ガンマ関数 - fortran66のブログ

Paul J. Nahin "Inside Interesting Integrals"

Inside Interesting Integrals: A Collection of Sneaky Tricks, Sly Substitutions, and Numerous Other Stupendously Clever, Awesomely Wicked, and Devilishly ... Lecture Notes in Physics) (English Edition)

作者:Nahin, Paul J.
Springer

Amazon

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 から組み込み関数になったガンマ関数と比較してみます。

ソース・プログラム

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

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
続行するには何かキーを押してください . . .