Juliaによる積分公式の検算

JuliaLang Advent Calendar 2023 3日目です🎄

はじめに

時間に依存しないSchrödinger方程式に変分法を適用する, すなわちエネルギーを波動関数の汎関数 $E[\psi] = \frac{\langle\psi|\hat{H}|\psi\rangle}{\langle\psi|\psi\rangle}$ とみなして最小化問題を考えると, 試行波動関数 $\psi=\sum c_i\phi_i$ の線形パラメータ $c_i$ は一般化固有値問題 $\boldsymbol{H}\boldsymbol{C} = E\boldsymbol{S}\boldsymbol{C}$ によって決定できる. 詳細はこちらの記事を参照されたい. このとき用いられる基底関数 $\phi$ は問題によって異なるが, 自乗可積分であることを要求すると, たいていはガウス関数か指数関数, あるいはそれらに少し手を加えたものが採用される. それらの基底関数を採用した場合, 行列要素 $H_{ij}$ と $S_{ij}$ の計算にガウス関数に関する積分公式や指数関数に関する積分公式が利用される. 問題は, これらの公式が正しいことをどうやって確かめるのか？ということである. 不定積分の場合は数式処理システム（例えばWolfram Alpha）等を用いて微分してみれば簡単に検証できるが, 広義積分では同じ方法が使えない（適切な逆問題は存在しないのだろうか？）. そこで, このノートではガウス関数とその一般化に関する積分公式について数値積分を用いた検証を行う.

環境

versioninfo()

Julia Version 1.8.3
Commit 0434deb161 (2022-11-14 20:14 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 8 × 11th Gen Intel(R) Core(TM) i7-1185G7 @ 3.00GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, tigerlake)
  Threads: 1 on 8 virtual cores

パッケージ

ガウス積分の積分領域には無限大が含まれるため, 数値積分のパッケージでは100とか1000とか適当な大きな数を入れて有限区間を指定することになる. しかし, QuadGK.jlでは積分領域にInfや-Infを使える（いい感じに処理してくれる）ことから, このノートでは全面的にQuadGK.jlを採用した. なお, 多変数の場合はHCubature.jlやMonteCarloIntegration.jlの方が有利になるため, これらのテストも兼ねて3つのパッケージの比較も行う. また, 公式の実装にSpecialFunctions.jlとGSL.jlを使用する. 初めての場合は using Pkg; Pkg.add("パッケージ名") を実行してインストールを行う.

using Pkg
Pkg.add("QuadGK")
Pkg.add("HCubature")
Pkg.add("MonteCarloIntegration")
Pkg.add("SpecialFunctions")
Pkg.add("GSL")

使用前は毎回 using パッケージ名 を実行してパッケージを読み込む.

using Printf
using QuadGK
using HCubature
using MonteCarloIntegration
using SpecialFunctions
using GSL

まずは簡単な積分を題材としてQuadGK.jl, HCubature.jl, MonteCarloIntegration.jlのそれぞれの使用例を比較する.

$$I = \int_0^2 x^3 ~\mathrm{d}x = \left[ \frac{1}{4}x^4 \right]_0^2 = 4$$

I₁ = quadgk(x -> x^3, 0, 2)[1]                      # QuadGK.jl
I₂ = hcubature(x -> x[1]^3, [0], [2])[1]            # HCubature.jl
I₃ = vegas(x -> x[1]^3, [0], [2]).integral_estimate # MonteCarloIntegration.jl

@show I₁
@show I₂
@show I₃

# output
I₁ = 4.0
I₂ = 3.999999999999999
I₃ = 4.000495416149978

いずれも概ね4を返した. どのパッケージも積分値だけでなく誤差やその他の情報を一緒に返してくるため, 最後に[1]や.integral_estimateを加えて積分値を抜き出している. また, HCubature.jl, MonteCarloIntegration.jlは多変数を前提として開発されているので, 1変数の場合でもxは配列であると考えて $x^3$ は x[1]^3 と書く. 次に $D=\left\{(x,y) | 0\leq x \leq 2, 0\leq y \leq 1 \right\}$ として2変数の積分を考える.（たしかフビニの定理を使ってこんな感じに解けたはず.）

$$\begin{aligned} I &= \iint_D x^3 + y ~\mathrm{d}x \mathrm{d}y \\ &= \int_0^1 \left[ \int_0^2 x^3 + y ~\mathrm{d}x \right] \mathrm{d}y \\ &= \int_0^1 \left[ \frac{1}{4}x^4 + xy \right]_0^2 \mathrm{d}y \\ &= \int_0^1 4 + 2y \mathrm{d}y \\ &= \left[ 4y + 2\cdot\frac{1}{2}y^2 \right]_0^1 \\ &= 5 \end{aligned}$$

I₁ = (
    quadgk(y ->
    quadgk(x ->
    x^3 + y
    , 0, 2, maxevals=10^3)[1]
    , 0, 1, maxevals=10^3)[1]
)                                                              # QuadGK.jl
I₂ = hcubature(x -> x[1]^3 + x[2], [0,0], [2,1])[1]            # HCubature.jl
I₃ = vegas(x -> x[1]^3 + x[2], [0,0], [2,1]).integral_estimate # MonteCarloIntegration.jl

@show I₁
@show I₂
@show I₃

# output
I₁ = 5.0
I₂ = 4.999999999999999
I₃ = 5.00545791767893

いずれも概ね5を返した. QuadGK.jlには1変数のルーチンしかないはずなので, 2重に使って2変数に対応する. 以下では積分公式を検証していく.

ガウス積分

証明は非常に有名であり, Wikipediaにも載っているので省略する.

$$\int_{-\infty}^{\infty} \exp(-a x^2) ~\mathrm{d}x = \sqrt{\frac{\pi}{a}}$$

下記のように解析的な積分と数値積分との誤差率が0.01%以内であれば✔を表示する. $a$ の値を5通り変えて検証を行った.

println("a      | analytical        | QuadGK.jl        ")
println("------ | ----------------- | -----------------")
for a in [0.01, 0.1, 1.0, 10.0, 100.0]
    analytical = sqrt(π/a)
    numerical  = quadgk(x -> exp(-a*x^2), -Inf, Inf, rtol=1e-12, maxevals=10^7)[1]
    @printf("%6.2f | %17.12f | %17.12f %s\n", a, analytical, numerical, isapprox(analytical, numerical, rtol=1e-5) ? "✔" :  "✗")
end

a      | analytical        | QuadGK.jl        
------ | ----------------- | -----------------
  0.01 |   17.724538509055 |   17.724538509055 ✔
  0.10 |    5.604991216398 |    5.604991216398 ✔
  1.00 |    1.772453850906 |    1.772453850906 ✔
 10.00 |    0.560499121640 |    0.560499121640 ✔
100.00 |    0.177245385091 |    0.177245385091 ✔

一般化されたガウス積分

1次元において, ガウス積分を下記のように一般化する.

$$\mathrm{GGI}(n,a) = \int_0^{\infty} x^{n} \exp \left(-a x^2\right)$$

ここでは通常のガウス積分は$n=0$の特殊例であると考える. 森口繁一, 宇田川銈久, 一松信『岩波 数学公式Ⅰ』(岩波書店,1988,新装第3刷) p.233では$n$が偶数の場合と奇数の場合で分けた公式

$$\int_0^{\infty} x^{2n} \exp \left(-a x^2\right) = \frac{(2n-1)!!}{2^{n+1}} \sqrt{\frac{\pi}{a^{2n+1}}}$$

$$\int_0^{\infty} x^{2n+1} \exp \left(-a x^2\right) = \frac{n!}{2 a^{n+1}}$$

が与えられているので$n! = \Gamma\left( n+1 \right)$および$\Gamma\left(\frac{1}{2}+n\right)=\frac{(2 n-1) ! !}{2^n} \sqrt{\pi}$を用いてまとめると

$$\int_0^{\infty} x^{n} \exp \left(-a x^2\right) = \frac{\Gamma\left( \frac{n+1}{2} \right)}{2 a^{\frac{n+1}{2}}}$$

が得られる. 英語版のWikipediaに載っているのを見つけた. また $a$ が具体的な数値の場合ならWolfram Alphaでも計算できるので適当な素数を代入すれば公式として取り出せる.

# using SpecialFunctions # for gamma()

GGI(n,a) = gamma((n+1)/2) / (2*a^((n+1)/2))

println("a      |   n | analytical        | QuadGK.jl        ")
println("------ | --- | ----------------- | -----------------")
for a in [0.5, 1.0, 10.0, 100.0]
    for n in 0:10
        analytical = GGI(n,a)
        numerical  = quadgk(x -> x^n * exp(-a*x^2), 0, Inf, maxevals=10^3)[1]
        @printf("%6.2f | %3d | %17.12f | %17.12f %s\n", a, n, analytical, numerical, isapprox(analytical, numerical, rtol=1e-5) ? "✔" :  "✗")
    end
end

a      |   n | analytical        | QuadGK.jl        
------ | --- | ----------------- | -----------------
  0.50 |   0 |    1.253314137316 |    1.253314137316 ✔
  0.50 |   1 |    1.000000000000 |    1.000000000000 ✔
  0.50 |   2 |    1.253314137316 |    1.253314137316 ✔
  0.50 |   3 |    2.000000000000 |    2.000000000000 ✔
  0.50 |   4 |    3.759942411947 |    3.759942411947 ✔
  0.50 |   5 |    8.000000000000 |    8.000000000000 ✔
  0.50 |   6 |   18.799712059733 |   18.799712059733 ✔
  0.50 |   7 |   48.000000000000 |   48.000000000002 ✔
  0.50 |   8 |  131.597984418128 |  131.597984418140 ✔
  0.50 |   9 |  384.000000000000 |  384.000000000000 ✔
  0.50 |  10 | 1184.381859763148 | 1184.381859763148 ✔
  1.00 |   0 |    0.886226925453 |    0.886226925453 ✔
  1.00 |   1 |    0.500000000000 |    0.500000000000 ✔
  1.00 |   2 |    0.443113462726 |    0.443113462726 ✔
  1.00 |   3 |    0.500000000000 |    0.500000000000 ✔
  1.00 |   4 |    0.664670194090 |    0.664670194090 ✔
  1.00 |   5 |    1.000000000000 |    1.000000000000 ✔
  1.00 |   6 |    1.661675485224 |    1.661675485224 ✔
  1.00 |   7 |    3.000000000000 |    3.000000000000 ✔
  1.00 |   8 |    5.815864198284 |    5.815864198284 ✔
  1.00 |   9 |   12.000000000000 |   12.000000000000 ✔
  1.00 |  10 |   26.171388892277 |   26.171388892277 ✔
 10.00 |   0 |    0.280249560820 |    0.280249560820 ✔
 10.00 |   1 |    0.050000000000 |    0.050000000000 ✔
 10.00 |   2 |    0.014012478041 |    0.014012478041 ✔
 10.00 |   3 |    0.005000000000 |    0.005000000000 ✔
 10.00 |   4 |    0.002101871706 |    0.002101871706 ✔
 10.00 |   5 |    0.001000000000 |    0.001000000000 ✔
 10.00 |   6 |    0.000525467927 |    0.000525467927 ✔
 10.00 |   7 |    0.000300000000 |    0.000300000000 ✔
 10.00 |   8 |    0.000183913774 |    0.000183913774 ✔
 10.00 |   9 |    0.000120000000 |    0.000120000000 ✔
 10.00 |  10 |    0.000082761198 |    0.000082761198 ✔
100.00 |   0 |    0.088622692545 |    0.088622692545 ✔
100.00 |   1 |    0.005000000000 |    0.005000000000 ✔
100.00 |   2 |    0.000443113463 |    0.000443113463 ✔
100.00 |   3 |    0.000050000000 |    0.000050000000 ✔
100.00 |   4 |    0.000006646702 |    0.000006646702 ✔
100.00 |   5 |    0.000001000000 |    0.000001000000 ✔
100.00 |   6 |    0.000000166168 |    0.000000166168 ✔
100.00 |   7 |    0.000000030000 |    0.000000030000 ✔
100.00 |   8 |    0.000000005816 |    0.000000005816 ✔
100.00 |   9 |    0.000000001200 |    0.000000001200 ✔
100.00 |  10 |    0.000000000262 |    0.000000000262 ✔

もっと一般化されたガウス積分

More Generalized Gaussian Integralを次のように定義する.

$$\begin{aligned} \mathrm{MGGI}(n,a,b) &= \int_0^{\infty} x^{n} \exp \left(-a x -b x^2\right) \mathrm{d}x \\ &= \frac{a}{2^{n+2} b^{n/2+1}} U\left(\frac{n}{2}+1, \frac{3}{2}, \frac{a^2}{4 b}\right) \Gamma(n+1) \end{aligned}$$

ここで $U$ は合流型超幾何関数である. これも $a$, $b$ が具体的な数値の場合ならWolfram Alphaでも計算でき, 適当な素数を代入すれば公式として取り出せる.

# using SpecialFunctions # for gamma()
# using GSL # for sf_hyperg_U()

# MGGI(n,a,b) = a*gamma(1+n)/(2^(2+n)*b^(1+n/2)) * HypergeometricFunctions.U(1+n/2,3/2,a^2/4/b)
MGGI(n,a,b) = a/(2^(n+2)*b^(n/2+1)) * GSL.sf_hyperg_U(n/2+1,3/2,a^2/4/b) * gamma(n+1)

println("a      | b      |   n | analytical        | QuadGK.jl        ")
println("------ | ------ | --- | ----------------- | -----------------")
for a in [0.1, 1.0, 10.0]
    for b in [0.1, 1.0, 10.0]
        for n in 0:10
            analytical = MGGI(n,a,b)
            numerical  = quadgk(x -> x^n * exp(-a*x - b*x^2), 0, Inf, maxevals=10^3)[1]
            @printf("%6.2f | %6.2f | %3d | %17.12f | %17.12f %s\n", a, b, n, analytical, numerical, isapprox(analytical, numerical, rtol=1e-5) ? "✔" :  "✗")
        end
    end
end

a      | b      |   n | analytical        | QuadGK.jl        
------ | ------ | --- | ----------------- | -----------------
  0.10 |   0.10 |   0 |    2.365023857063 |    2.365023857063 ✔
  0.10 |   0.10 |   1 |    3.817488071469 |    3.817488071469 ✔
  0.10 |   0.10 |   2 |    9.916375249581 |    9.916375249581 ✔
  0.10 |   0.10 |   3 |   33.216693089895 |   33.216693089895 ✔
  0.10 |   0.10 |   4 |  132.137282198762 |  132.137282198762 ✔
  0.10 |   0.10 |   5 |  598.265220698515 |  598.265220698516 ✔
  0.10 |   0.10 |   6 | 3004.299444619790 | 3004.299444619781 ✔
  0.10 |   0.10 |   7 | 16445.806898645555 | 16445.806898645282 ✔
  0.10 |   0.10 |   8 | 96927.577112369851 | 96927.577112362909 ✔
  0.10 |   0.10 |   9 | 609368.487389638554 | 609368.487389637972 ✔
  0.10 |   0.10 |  10 | 4057056.726361821406 | 4057056.726361818612 ✔
  0.10 |   1.00 |   0 |    0.838361847809 |    0.838361847809 ✔
  0.10 |   1.00 |   1 |    0.458081907610 |    0.458081907610 ✔
  0.10 |   1.00 |   2 |    0.396276828524 |    0.396276828524 ✔
  0.10 |   1.00 |   3 |    0.438268066183 |    0.438268066183 ✔
  0.10 |   1.00 |   4 |    0.572501839477 |    0.572501839477 ✔
  0.10 |   1.00 |   5 |    0.847911040393 |    0.847911040393 ✔
  0.10 |   1.00 |   6 |    1.388859046672 |    1.388859046672 ✔
  0.10 |   1.00 |   7 |    2.474290168845 |    2.474290168845 ✔
  0.10 |   1.00 |   8 |    4.737292154909 |    4.737292154909 ✔
  0.10 |   1.00 |   9 |    9.660296067635 |    9.660296067635 ✔
  0.10 |   1.00 |  10 |   20.834799893709 |   20.834799893709 ✔
  0.10 |  10.00 |   0 |    0.275318798552 |    0.275318798552 ✔
  0.10 |  10.00 |   1 |    0.048623406007 |    0.048623406007 ✔
  0.10 |  10.00 |   2 |    0.013522822898 |    0.013522822898 ✔
  0.10 |  10.00 |   3 |    0.004794726486 |    0.004794726486 ✔
  0.10 |  10.00 |   4 |    0.002004449802 |    0.002004449802 ✔
  0.10 |  10.00 |   5 |    0.000948923048 |    0.000948923048 ✔
  0.10 |  10.00 |   6 |    0.000496367835 |    0.000496367835 ✔
  0.10 |  10.00 |   7 |    0.000282195075 |    0.000282195075 ✔
  0.10 |  10.00 |   8 |    0.000172317767 |    0.000172317767 ✔
  0.10 |  10.00 |   9 |    0.000112016441 |    0.000112016441 ✔
  0.10 |  10.00 |  10 |    0.000076982913 |    0.000076982913 ✔
  1.00 |   0.10 |   0 |    0.865392586515 |    0.865392586515 ✔
  1.00 |   0.10 |   1 |    0.673037067424 |    0.673037067424 ✔
  1.00 |   0.10 |   2 |    0.961777595453 |    0.961777595453 ✔
  1.00 |   0.10 |   3 |    1.921482696980 |    1.921482696980 ✔
  1.00 |   0.10 |   4 |    4.819250446898 |    4.819250446898 ✔
  1.00 |   0.10 |   5 |   14.333401705099 |   14.333401705099 ✔
  1.00 |   0.10 |   6 |   48.814252646962 |   48.814252646962 ✔
  1.00 |   0.10 |   7 |  185.930787918163 |  185.930787918163 ✔
  1.00 |   0.10 |   8 |  778.844903052866 |  778.844903052855 ✔
  1.00 |   0.10 |   9 | 3543.007001462193 | 3543.007001462197 ✔
  1.00 |   0.10 |  10 | 17332.985630067957 | 17332.985630067997 ✔
  1.00 |   1.00 |   0 |    0.545641360765 |    0.545641360765 ✔
  1.00 |   1.00 |   1 |    0.227179319617 |    0.227179319617 ✔
  1.00 |   1.00 |   2 |    0.159231020574 |    0.159231020574 ✔
  1.00 |   1.00 |   3 |    0.147563809331 |    0.147563809331 ✔
  1.00 |   1.00 |   4 |    0.165064626195 |    0.165064626195 ✔
  1.00 |   1.00 |   5 |    0.212595305563 |    0.212595305563 ✔
  1.00 |   1.00 |   6 |    0.306363912707 |    0.306363912707 ✔
  1.00 |   1.00 |   7 |    0.484603960337 |    0.484603960337 ✔
  1.00 |   1.00 |   8 |    0.829971714305 |    0.829971714305 ✔
  1.00 |   1.00 |   9 |    1.523429984196 |    1.523429984196 ✔
  1.00 |   1.00 |  10 |    2.973157722275 |    2.973157722275 ✔
  1.00 |  10.00 |   0 |    0.236502385706 |    0.236502385706 ✔
  1.00 |  10.00 |   1 |    0.038174880715 |    0.038174880715 ✔
  1.00 |  10.00 |   2 |    0.009916375250 |    0.009916375250 ✔
  1.00 |  10.00 |   3 |    0.003321669309 |    0.003321669309 ✔
  1.00 |  10.00 |   4 |    0.001321372822 |    0.001321372822 ✔
  1.00 |  10.00 |   5 |    0.000598265221 |    0.000598265221 ✔
  1.00 |  10.00 |   6 |    0.000300429944 |    0.000300429944 ✔
  1.00 |  10.00 |   7 |    0.000164458069 |    0.000164458069 ✔
  1.00 |  10.00 |   8 |    0.000096927577 |    0.000096927577 ✔
  1.00 |  10.00 |   9 |    0.000060936849 |    0.000060936849 ✔
  1.00 |  10.00 |  10 |    0.000040570567 |    0.000040570567 ✔
 10.00 |   0.10 |   0 |    0.099801188165 |    0.099801188165 ✔
 10.00 |   0.10 |   1 |    0.009940591748 |    0.009940591748 ✔
 10.00 |   0.10 |   2 |    0.001976353427 |    0.001976353427 ✔
 10.00 |   0.10 |   3 |    0.000588246113 |    0.000588246113 ✔
 10.00 |   0.10 |   4 |    0.000232995745 |    0.000232995745 ✔
 10.00 |   0.10 |   5 |    0.000115135036 |    0.000115135036 ✔
 10.00 |   0.10 |   6 |    0.000068141803 |    0.000068141803 ✔
 10.00 |   0.10 |   7 |    0.000046960922 |    0.000046960922 ✔
 10.00 |   0.10 |   8 |    0.000036916993 |    0.000036916993 ✔
 10.00 |   0.10 |   9 |    0.000032587264 |    0.000032587264 ✔
 10.00 |   0.10 |  10 |    0.000031901484 |    0.000031901484 ✔
 10.00 |   1.00 |   0 |    0.098109430732 |    0.098109430732 ✔
 10.00 |   1.00 |   1 |    0.009452846342 |    0.009452846342 ✔
 10.00 |   1.00 |   2 |    0.001790483654 |    0.001790483654 ✔
 10.00 |   1.00 |   3 |    0.000500428071 |    0.000500428071 ✔
 10.00 |   1.00 |   4 |    0.000183585126 |    0.000183585126 ✔
 10.00 |   1.00 |   5 |    0.000082930513 |    0.000082930513 ✔
 10.00 |   1.00 |   6 |    0.000044310249 |    0.000044310249 ✔
 10.00 |   1.00 |   7 |    0.000027240296 |    0.000027240296 ✔
 10.00 |   1.00 |   8 |    0.000018884388 |    0.000018884388 ✔
 10.00 |   1.00 |   9 |    0.000014539246 |    0.000014539246 ✔
 10.00 |   1.00 |  10 |    0.000012283518 |    0.000012283518 ✔
 10.00 |  10.00 |   0 |    0.086539258652 |    0.086539258652 ✔
 10.00 |  10.00 |   1 |    0.006730370674 |    0.006730370674 ✔
 10.00 |  10.00 |   2 |    0.000961777595 |    0.000961777595 ✔
 10.00 |  10.00 |   3 |    0.000192148270 |    0.000192148270 ✔
 10.00 |  10.00 |   4 |    0.000048192504 |    0.000048192504 ✔
 10.00 |  10.00 |   5 |    0.000014333402 |    0.000014333402 ✔
 10.00 |  10.00 |   6 |    0.000004881425 |    0.000004881425 ✔
 10.00 |  10.00 |   7 |    0.000001859308 |    0.000001859308 ✔
 10.00 |  10.00 |   8 |    0.000000778845 |    0.000000778845 ✔
 10.00 |  10.00 |   9 |    0.000000354301 |    0.000000354301 ✔
 10.00 |  10.00 |  10 |    0.000000173330 |    0.000000173330 ✔

2次元のガウス積分

以下では, 上記の2の一般化とは異なり, 次元を増やすという意味での一般化を考える. 2次元の積分は逐次的に$x,y$の順番で積分を計算していけば, 単に1次元のガウス積分の2乗になることがわかる. $x$による積分の部分は$y$による積分から見れば定数なので, くくり出せる. なお1変数のガウス積分の証明で2変数の積分を使っているのだから循環論法になっている気もするが, 他にも証明方法があったはずなので気にしない.

$$\begin{aligned} \iint_{\mathbb{R}^2} \exp(-a r^2) ~ \mathrm{d}\pmb{r} &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \exp(-a x^2 -a y^2) ~ \mathrm{d}x \mathrm{d}y \\ &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \exp(-a x^2) \exp(-a y^2) ~ \mathrm{d}x \mathrm{d}y \\ &= \left[ \int_{-\infty}^{\infty} \exp(-a x^2) \mathrm{d}x \right] \int_{-\infty}^{\infty} \exp(-a y^2) \mathrm{d}y \\ &= \sqrt{\frac{\pi}{a}} \int_{-\infty}^{\infty} \exp(-a y^2) \mathrm{d}y \\ &= \sqrt{\frac{\pi}{a}}^2 \end{aligned}$$

println("a      | analytical        | QuadGK.jl        ")
println("------ | ----------------- | -----------------")
for a in [0.01, 0.1, 1.0, 10.0, 100.0]
    analytical = sqrt(π/a)^2
    numerical  = (
        quadgk(y ->
        quadgk(x ->
        exp(-a*x^2 -a*y^2)
        , -Inf, Inf, maxevals=10^3)[1]
        , -Inf, Inf, maxevals=10^3)[1]
    )
    @printf("%6.2f | %17.12f | %17.12f %s\n", a, analytical, numerical, isapprox(analytical, numerical, rtol=1e-5) ? "✔" :  "✗")
end

a      | analytical        | QuadGK.jl        
------ | ----------------- | -----------------
  0.01 |  314.159265358979 |  314.159265358887 ✔
  0.10 |   31.415926535898 |   31.415926535898 ✔
  1.00 |    3.141592653590 |    3.141592653590 ✔
 10.00 |    0.314159265359 |    0.314159265359 ✔
100.00 |    0.031415926536 |    0.031415926536 ✔

極座標系におけるガウス積分

次式で与えられる極座標系への変換$(r,\theta)\mapsto(x,y)$の下でガウス積分を計算する.

$$\begin{aligned} \left\{ \begin{aligned} x =& r \cos \theta \\ y =& r \sin \theta \\ \end{aligned} \right. \end{aligned}$$

ヤコビアンは下記のように計算できる.

$$\begin{aligned} \mathrm{d}x \mathrm{d}y &= \left|\frac{\partial(x,y)}{\partial(r,\theta)}\right| \mathrm{d}r \mathrm{d}\theta \\ &= \left|\det\left(\begin{array}{rr} \frac{\partial x}{\partial r} & \frac{\partial x}{\partial \theta} \\ \frac{\partial y}{\partial r} & \frac{\partial y}{\partial \theta} \\ \end{array}\right)\right| \mathrm{d}r \mathrm{d}\theta \\ &= \left|\det\left(\begin{array}{rr} \cos\theta & -r\sin\theta \\ \sin\theta & r\cos\theta \\ \end{array}\right)\right| \mathrm{d}r \mathrm{d}\theta \\ &= \left| \cos\theta \times r\cos\theta - \sin\theta \times (-r\sin\theta) \right| \mathrm{d}r \mathrm{d}\theta \\ &= \left| r\cos^2\theta + r\sin^2\theta \right| \mathrm{d}r \mathrm{d}\theta \\ &= \left| r(\cos^2\theta + \sin^2\theta) \right| \mathrm{d}r \mathrm{d}\theta \\ &= \left| r \right| \mathrm{d}r \mathrm{d}\theta \\ &= r \mathrm{d}r \mathrm{d}\theta \\ \end{aligned}$$

なお, $0\leq r<\infty$であるから$|r|=r$としている. これを用いて, 先ほどと同様に1次元のガウス積分の2乗になることが示せる.

$$\begin{aligned} \iint_{\mathbb{R}^2} \exp(-a r^2) ~ \mathrm{d}\pmb{r} &= \int_{0}^{2\pi} \int_{0}^{\infty} \exp(-a r^2) ~ r ~ \mathrm{d}r \mathrm{d}\theta \\ &= \int_{0}^{2\pi} \mathrm{d}\theta \int_{0}^{\infty} r \exp(-a r^2) ~ \mathrm{d}r \\ &= \left[ \theta \right]_{0}^{2\pi} \frac{0!}{2 a^{0+1}} \\ &= 2\pi \frac{1}{2a} \\ &= \sqrt{\frac{\pi}{a}}^2 \\ \end{aligned}$$

println("a      | analytical        | QuadGK.jl        ")
println("------ | ----------------- | -----------------")
for a in [0.01, 0.1, 1.0, 10.0, 100.0]
    analytical = sqrt(π/a)^2
    numerical  = (
        quadgk(θ ->
        quadgk(r ->
        r * exp(-a*r^2)
        , 0, Inf, maxevals=10^3)[1]
        , 0, 2*π, maxevals=10^3)[1]
    )
    @printf("%6.2f | %17.12f | %17.12f %s\n", a, analytical, numerical, isapprox(analytical, numerical, rtol=1e-5) ? "✔" :  "✗")
end

a      | analytical        | QuadGK.jl        
------ | ----------------- | -----------------
  0.01 |  314.159265358979 |  314.159265358979 ✔
  0.10 |   31.415926535898 |   31.415926535898 ✔
  1.00 |    3.141592653590 |    3.141592653590 ✔
 10.00 |    0.314159265359 |    0.314159265359 ✔
100.00 |    0.031415926536 |    0.031415926536 ✔

3次元のガウス積分

3次元でも2次元と同様にして, 単に1次元のガウス積分の3乗となる.

$$\begin{aligned} \iiint_{\mathbb{R}^3} \exp(-a r^2) ~ \mathrm{d}\pmb{r} &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \exp(-a x^2 -a y^2 -a z^2) ~ \mathrm{d}x \mathrm{d}y \mathrm{d}z \\ &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \exp(-a x^2) \exp(-a y^2) \exp(-a z^2) ~ \mathrm{d}x \mathrm{d}y \mathrm{d}z \\ &= \left[ \int_{-\infty}^{\infty} \exp(-a x^2) \mathrm{d}x \right] ~ \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \exp(-a y^2) \exp(-a z^2) \mathrm{d}y \mathrm{d}z \\ &= \sqrt{\frac{\pi}{a}} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \exp(-a y^2) \exp(-a z^2) \mathrm{d}y \mathrm{d}z \\ &= \sqrt{\frac{\pi}{a}} \left[ \int_{-\infty}^{\infty} \exp(-a y^2) \mathrm{d}y \right] \int_{-\infty}^{\infty} \exp(-a z^2) \mathrm{d}z \\ &= \sqrt{\frac{\pi}{a}}^2 \int_{-\infty}^{\infty} \exp(-a z^2) \mathrm{d}z \\ &= \sqrt{\frac{\pi}{a}}^3 \end{aligned}$$

println("a      | analytical        | QuadGK.jl        ")
println("------ | ----------------- | -----------------")
for a in [0.01, 0.1, 1.0, 10.0, 100.0]
    analytical = sqrt(π/a)^3
    numerical  = (
        quadgk(z ->
        quadgk(y ->
        quadgk(x ->
        exp(-a*x^2 -a*y^2 -a*z^2)
        , -Inf, Inf, maxevals=10^3)[1]
        , -Inf, Inf, maxevals=10^3)[1]
        , -Inf, Inf, maxevals=10^3)[1]
    )
    @printf("%6.2f | %17.12f | %17.12f %s\n", a, analytical, numerical, isapprox(analytical, numerical, rtol=1e-5) ? "✔" :  "✗")
end

a      | analytical        | QuadGK.jl        
------ | ----------------- | -----------------
  0.01 | 5568.327996831709 | 5568.327996829265 ✔
  0.10 |  176.085992288711 |  176.085992288711 ✔
  1.00 |    5.568327996832 |    5.568327996832 ✔
 10.00 |    0.176085992289 |    0.176085992289 ✔
100.00 |    0.005568327997 |    0.005568327997 ✔

球面座標系におけるガウス積分

球面座標系への変換

$$\begin{cases} x = r\sin\theta\cos\varphi\\ y = r\sin\theta\sin\varphi\\ z = r\cos\theta \end{cases}$$

におけるヤコビアンは

$$\begin{aligned} \mathrm{d}x \mathrm{d}y \mathrm{d}z &= \left|\frac{\partial(x,y,z)}{\partial(r,\theta,\varphi)}\right| \mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= \left|\det\left(\begin{array}{rr} \frac{\partial x}{\partial r} & \frac{\partial x}{\partial \theta} & \frac{\partial x}{\partial \varphi} \\ \frac{\partial y}{\partial r} & \frac{\partial y}{\partial \theta} & \frac{\partial y}{\partial \varphi} \\ \frac{\partial z}{\partial r} & \frac{\partial z}{\partial \theta} & \frac{\partial z}{\partial \varphi} \\ \end{array}\right)\right| \mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= \cdots \\ &= r^2 \sin(\theta) \mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ \end{aligned}$$

のように計算できるから, これを用いて積分領域$0\leq r<\infty, 0\leq \theta \leq \pi, \leq \varphi<2 \pi$で積分すれば

$$\begin{aligned} \iiint_{\mathbb{R}^3} \exp(-a r^2) ~ \mathrm{d}\pmb{r} &= \int_{0}^{2\pi} \int_{0}^{\pi} \int_{0}^{\infty} \exp(-a r^2) ~ r^2 \sin\theta ~ \mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= 4\pi \int_{0}^{\infty} r^2 \exp(-a r^2) ~ \mathrm{d}r \\ &= 4\pi \frac{(2\cdot1-1)!!}{2^{1+1}} \sqrt{\frac{\pi}{a^{2\cdot1+1}}} \\ &= 4\pi \frac{1}{4} \sqrt{\frac{\pi}{a^{3}}} \\ &= \sqrt{\frac{\pi}{a}}^3 \end{aligned}$$

と計算できる. なお, 計算には1次元のガウス積分を用いている.

println("a      | analytical        | QuadGK.jl        ")
println("------ | ----------------- | -----------------")
for a in [0.01, 0.1, 1.0, 10.0, 100.0]
    analytical = sqrt(π/a)^3
    numerical  = (
        quadgk(φ ->
        quadgk(θ ->
        quadgk(r ->
        r^2 * sin(θ) * exp(-a*r^2)
        , 0,  Inf, maxevals=100)[1]
        , 0,   π, maxevals=10)[1]
        , 0, 2*π, maxevals=20)[1]
    )
    @printf("%6.2f | %17.12f | %17.12f %s\n", a, analytical, numerical, isapprox(analytical, numerical, rtol=1e-5) ? "✔" :  "✗")
end

a      | analytical        | QuadGK.jl        
------ | ----------------- | -----------------
  0.01 | 5568.327996831709 | 5567.958330629872 ✗
  0.10 |  176.085992288711 |  176.085992503519 ✔
  1.00 |    5.568327996832 |    5.568327996978 ✔
 10.00 |    0.176085992289 |    0.176085992289 ✔
100.00 |    0.005568327997 |    0.005568327997 ✔

最後にHCubature.jlとMonteCarloIntegration.jlも加えた結果を示す. モンテカルロ法では誤差率の閾値は甘めに1%とした.

println("a      | analytical        | QuadGK.jl            | HCubature.jl         | MonteCarloIntegration.jl")
println("------ | ----------------- | -----------------    | -----------------    | -----------------")
for a in [0.01, 0.1, 1.0, 10.0, 100.0]
    analytical = sqrt(π/a)^3
    I₁  = (
        quadgk(φ ->
        quadgk(θ ->
        quadgk(r ->
        r^2 * sin(θ) * exp(-a*r^2)
        , 0,  Inf, maxevals=200)[1]
        , 0,   π, maxevals=10)[1]
        , 0, 2*π, maxevals=20)[1]
    )
    I₂ = hcubature(x -> x[1]^2 * sin(x[2]) * exp(-a*x[1]^2), [0,0,0], [100,π,2*π])[1]
    I₃ = vegas(x -> x[1]^2 * sin(x[2]) * exp(-a*x[1]^2), [0,0,0], [100,π,2*π]).integral_estimate
    @printf("%6.2f | %17.12f | %17.12f %s | %17.12f %s | %17.12f %s \n", a, analytical,
        I₁, isapprox(analytical, I₁, rtol=1e-5) ? "✔" :  "✗",
        I₂, isapprox(analytical, I₂, rtol=1e-5) ? "✔" :  "✗",
        I₃, isapprox(analytical, I₃, rtol=1e-2) ? "✔" :  "✗"
    )
end

a      | analytical        | QuadGK.jl            | HCubature.jl         | MonteCarloIntegration.jl
------ | ----------------- | -----------------    | -----------------    | -----------------
  0.01 | 5568.327996831709 | 5568.327996831707 ✔ | 5568.327997006358 ✔ | 5568.478496127373 ✔ 
  0.10 |  176.085992288711 |  176.085992288711 ✔ |  176.085992294186 ✔ |  176.120007666057 ✔ 
  1.00 |    5.568327996832 |    5.568327996832 ✔ |    5.568327997105 ✔ |    5.561248789741 ✔ 
 10.00 |    0.176085992289 |    0.176085992289 ✔ |    0.176085992297 ✔ |    0.176020030094 ✔ 
100.00 |    0.005568327997 |    0.005568327997 ✔ |    0.005568327997 ✔ |    0.005567615795 ✔ 

参考文献

• Wikipedia 日本語版
• Wikipedia 英語版
• 森口繁一, 宇田川銈久, 一松信『岩波 数学公式Ⅰ』(岩波書店,1988,新装第3刷) p.233
• QuadGK.jlドキュメント
• HCubature.jlドキュメント
• MonteCarloIntegration.jlドキュメント
• Wolfram Alpha, GGI, MGGI

おわりに

いかがでしたか？ 元々は開発したパッケージの紹介を予定していましたが, 公式パッケージに登録する際になんと名前が似ているパッケージが既にあるという理由で弾かれてしまい,　使用方法やドキュメントに変更が生じるため紹介を延期しました. (間に合えば後日のアドベントカレンダーで紹介する予定です.)

数値積分については, 数式とそのまま対応したコードが書けるQuadGK.jlを気に入って使用しています. 積分の評価回数を減らしたり収束閾値を甘くするなどの工夫をしていますが, やはり3次元の積分が時間を食います. パッケージのテストに使っているためGitHub Actionsの無料時間の枠内で不安を感じない程度には高速化したい思っており, モンテカルロ積分の導入を検討し始めたところです. とはいえ, モンテカルロ積分でも定義域を広くとりすぎるとサンプリングがうまくいかないようなエラーもあるので, うまく積分領域のカットオフを決められないかと考えています.

最後に, 広義積分の公式について数値積分の他に何かよい検算の方法をご存じの方はご教授頂ければ幸いです.

https://gist.github.com/ohno/0162699a248de54aa1d93906c265b1b7