fortran66のブログ

fortran について書きます。

【メモ帳】exp(-x^4) の積分

(1/4)!

Youtube のおすすめで (1/4)! を計算する動画が出てきたのですが、 (1/4)! は \int_0^\infty \exp(-x^4)dx に等しいので見てみたところレムニスケート周率で表わされることが学べたので、メモっておくことにします。


www.youtube.com

exp(-x4) の積分

y=x^4 で置換すると dy = 4x^3dx なので、 \int_0^\infty \exp(-x^4)dx=\int_0^\infty \exp(-y)y^{-{3\over4}}dy/4={1\over4}\int_0^\infty \exp(-y)y^{{1\over4}-1}dy={1\over4}\Gamma({1\over4}) が示される。一般に ({1\over n})!=\int_0^\infty \exp(-x^n)dx となる。

fortran66.hatenablog.com

レムニスケート周率 \varpi

定義により、 \varpi=2 \int_0^1{dy\over\sqrt{1-y^4}}。その値は、

2.622057554292119810464839589891119413682754951431623162816821703...

oeis.org

Youtube での導出

  • Gamma 関数の相反公式より

\Gamma({1\over4}) \Gamma(1-{1\over4})=\Gamma({1\over4}) \Gamma({3\over4})={\pi\over\sin({\pi\over4})}

  • Beta 関数の性質から

B(1/4,1/2)={\Gamma({1\over4})\Gamma({1\over2})\over\Gamma({1\over4}+{1\over2})}={\Gamma({1\over4})\sqrt{\pi}\over\Gamma({3\over4})}

  • Beta 関数の定義から

ここでも y=x^4 で置換する。

B(1/4,1/2)=\int_0^1y^{{1\over4}-1}(1-y)^{{1\over2}-1}dy=\int_0^1x^{-3}(1-x^4)^{-{1\over2}}4x^3dx=4\int_0^1{1\over\sqrt{1-x^4}}dx=2\varpi

以上のことから、

B(1/4,1/2)={\Gamma({1\over4})^2\sqrt{\pi}\over\sqrt{2}\pi}={\Gamma({1\over4})^2\over\sqrt{2\pi}}

つまり

\Gamma({1\over4})^2=2\varpi\sqrt{2\pi}

したがって

({1\over4})!={1\over4}\Gamma({1\over4})={1\over4}\sqrt{2\varpi\sqrt{}2\pi}

結論

\int_0^\infty \exp(-x^4)dx=({1\over4})!={1\over4}\sqrt{2\varpi\sqrt{}2\pi}

寝言

\Gamma({1\over4})^2 が出てくるので、\int_0^\infty \exp(-x^2)dx の時のように、積分を二乗して球座標に行って解けないかとも思いましたが、今日は方角が悪くてやる気が出ません！Bata 関数の出方を見る限り難しい気もします。

Fortran で検算

[tex: \exp(-x4)] の積分は台形積分で、急速に減少する関数なので上限を 2 までにします。分割数は 32 。

    program Lemniscate
        implicit none
        real, parameter :: pi = 4 * atan(1.0)
        real, parameter :: aLemniscate = 2.622057554292119810
        integer, parameter :: n = 32
        integer :: i
        real :: x(0:n), y(0:n), s, xmax = 2.0, h 
        h = xmax / n
        
        print *, '1/4gamma(1/4)=              ', gamma(5.0/4.0), gamma(1.0/4.0)/4.0
        
        forall(i = 0:n) x(i) = h * i 
        y = exp(-x**4)
        
        s = (y(0) + y(n)) / 2 + sum(y(1:n-1)) ! Trapezoidal rule
        
        print *, '\int_0^\infty \exp(-x^4)dx= ', h * s

        print *, '1/4sqrt(2Lsqrt(2pi))=       ', sqrt(2 * aLemniscate * sqrt(2 * pi)) / 4
    end program Lemniscate

実行結果

 1/4gamma(1/4)=                0.9064025      0.9064025
 \int_0^\infty \exp(-x^4)dx=   0.9064025
 1/4sqrt(2Lsqrt(2pi))=         0.9064025

みんなのFortran―基礎から発展まで―

みんなのFortran―基礎から発展まで―

Amazon

C & Fortran 演習で学ぶ数値計算

C & Fortran 演習で学ぶ数値計算

Amazon

Modern Fortran: Building efficient parallel applications

Modern Fortran: Building efficient parallel applications

Amazon

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

Amazon