Skip to content

Instantly share code, notes, and snippets.

@ohno

ohno/GaussianIntegral.ipynb

Last active December 6, 2023 18:49
Show Gist options
  • Save ohno/0162699a248de54aa1d93906c265b1b7 to your computer and use it in GitHub Desktop.
Save ohno/0162699a248de54aa1d93906c265b1b7 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
GaussianIntegral.ipynb
{
"cells": [
{
"cell_type": "markdown",
"id": "5eb893df",
"metadata": {},
"source": [
"## はじめに\n",
"\n",
"時間に依存しない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}$ によって決定できる. 詳細は[こちら](https://zenn.dev/ohno/articles/c48920c9327b16)の記事を参照されたい. このとき用いられる基底関数 $\\phi$ は問題によって異なるが, 自乗可積分であることを要求すると, たいていはガウス関数か指数関数, あるいはそれらに少し手を加えたものが採用される. それらの基底関数を採用した場合, 行列要素 $H_{ij}$ と $S_{ij}$ の計算に[ガウス関数に関する積分公式](https://en.wikipedia.org/wiki/List_of_integrals_of_exponential_functions)や[指数関数に関する積分公式](https://en.wikipedia.org/wiki/List_of_integrals_of_Gaussian_functions)が利用される. 問題は, これらの公式が正しいことをどうやって確かめるのか？ということである. 不定積分の場合は数式処理システム（例えば[Wolfram Alpha](https://www.wolframalpha.com/input?i=d%2Fdx+1%2F2+x%5E2&lang=ja)）等を用いて微分してみれば簡単に検証できるが, 広義積分では同じ方法が使えない（適切な逆問題は存在しないのだろうか？）. そこで, このノートではガウス関数とその一般化に関する積分公式について数値積分を用いた検証を行う."
]
},
{
"cell_type": "markdown",
"id": "39fa3122",
"metadata": {},
"source": [
"## 環境"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "2c8a78c2",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Julia Version 1.8.3\n",
"Commit 0434deb161 (2022-11-14 20:14 UTC)\n",
"Platform Info:\n",
" OS: Windows (x86_64-w64-mingw32)\n",
" CPU: 8 × 11th Gen Intel(R) Core(TM) i7-1185G7 @ 3.00GHz\n",
" WORD_SIZE: 64\n",
" LIBM: libopenlibm\n",
" LLVM: libLLVM-13.0.1 (ORCJIT, tigerlake)\n",
" Threads: 1 on 8 virtual cores\n"
]
}
],
"source": [
"versioninfo()"
]
},
{
"cell_type": "markdown",
"id": "aca539aa",
"metadata": {},
"source": [
"## パッケージ\n",
"\n",
"ガウス積分の積分領域には無限大が含まれるため, 数値積分のパッケージでは100とか1000とか適当な大きな数を入れて有限区間を指定することになる. しかし, [QuadGK.jl](https://juliamath.github.io/QuadGK.jl/stable/)では積分領域に`Inf`や`-Inf`を使える（いい感じに処理してくれる）ことから, このノートでは全面的に[QuadGK.jl](https://juliamath.github.io/QuadGK.jl/stable/)を採用した. なお, 多変数の場合は[HCubature.jl](https://github.com/JuliaMath/HCubature.jl)や[MonteCarloIntegration.jl](https://github.com/ranjanan/MonteCarloIntegration.jl)の方が有利になるため, これらのテストも兼ねて3つのパッケージの比較も行う. また, 公式の実装に[SpecialFunctions.jl](https://specialfunctions.juliamath.org/stable/)と[GSL.jl](https://github.com/JuliaMath/GSL.jl)を使用する. 初めての場合は `using Pkg; Pkg.add(\"パッケージ名\")` を実行してインストールを行う."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "67e2ab55",
"metadata": {},
"outputs": [],
"source": [
"# using Pkg\n",
"# Pkg.add(\"QuadGK\")\n",
"# Pkg.add(\"HCubature\")\n",
"# Pkg.add(\"MonteCarloIntegration\")\n",
"# Pkg.add(\"SpecialFunctions\")\n",
"# Pkg.add(\"GSL\")"
]
},
{
"cell_type": "markdown",
"id": "9e4e78dd",
"metadata": {},
"source": [
"使用前は毎回 `using パッケージ名` を実行してパッケージを読み込む."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "4307ee85",
"metadata": {},
"outputs": [],
"source": [
"using Printf\n",
"using QuadGK\n",
"using HCubature\n",
"using MonteCarloIntegration\n",
"using SpecialFunctions\n",
"using GSL"
]
},
{
"cell_type": "markdown",
"id": "222015ce",
"metadata": {},
"source": [
"まずは簡単な積分を題材としてQuadGK.jl, HCubature.jl, MonteCarloIntegration.jlのそれぞれの使用例を比較する.\n",
"\n",
"$$\n",
"I\n",
"= \\int_0^2 x^3 ~\\mathrm{d}x\n",
"= \\left[ \\frac{1}{4}x^4 \\right]_0^2\n",
"= 4\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "2d7c61a4",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"I₁ = 4.0\n",
"I₂ = 3.999999999999999\n",
"I₃ = 3.998717327904603\n"
]
},
{
"data": {
"text/plain": [
"3.998717327904603"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"I₁ = quadgk(x -> x^3, 0, 2)[1] # QuadGK.jl\n",
"I₂ = hcubature(x -> x[1]^3, [0], [2])[1] # HCubature.jl\n",
"I₃ = vegas(x -> x[1]^3, [0], [2]).integral_estimate # MonteCarloIntegration.jl\n",
"\n",
"@show I₁\n",
"@show I₂\n",
"@show I₃"
]
},
{
"cell_type": "markdown",
"id": "4ac9b925",
"metadata": {},
"source": [
"いずれも概ね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変数の積分を考える.（たしかフビニの定理を使ってこんな感じに解けたはず.）\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"I\n",
"&= \\iint_D x^3 + y ~\\mathrm{d}x \\mathrm{d}y \\\\\n",
"&= \\int_0^1 \\left[ \\int_0^2 x^3 + y ~\\mathrm{d}x \\right] \\mathrm{d}y \\\\\n",
"&= \\int_0^1 \\left[ \\frac{1}{4}x^4 + xy \\right]_0^2 \\mathrm{d}y \\\\\n",
"&= \\int_0^1 4 + 2y \\mathrm{d}y \\\\\n",
"&= \\left[ 4y + 2\\cdot\\frac{1}{2}y^2 \\right]_0^1 \\\\\n",
"&= 5\n",
"\\end{aligned}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "1983269c",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"I₁ = 5.0\n",
"I₂ = 4.999999999999999\n",
"I₃ = 5.006635179715269\n"
]
},
{
"data": {
"text/plain": [
"5.006635179715269"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"I₁ = (\n",
" quadgk(y ->\n",
" quadgk(x ->\n",
" x^3 + y\n",
" , 0, 2, maxevals=10^3)[1]\n",
" , 0, 1, maxevals=10^3)[1]\n",
") # QuadGK.jl\n",
"I₂ = hcubature(x -> x[1]^3 + x[2], [0,0], [2,1])[1] # HCubature.jl\n",
"I₃ = vegas(x -> x[1]^3 + x[2], [0,0], [2,1]).integral_estimate # MonteCarloIntegration.jl\n",
"\n",
"@show I₁\n",
"@show I₂\n",
"@show I₃"
]
},
{
"cell_type": "markdown",
"id": "8992637f",
"metadata": {},
"source": [
"いずれも概ね5を返した. QuadGK.jlには1変数のルーチンしかないはずなので, 2重に使って2変数に対応する. 以下では積分公式を検証していく."
]
},
{
"cell_type": "markdown",
"id": "0c3c34bb",
"metadata": {},
"source": [
"## ガウス積分\n",
"\n",
"証明は非常に有名であり, [Wikipedia](https://ja.wikipedia.org/wiki/%E3%82%AC%E3%82%A6%E3%82%B9%E7%A9%8D%E5%88%86)にも載っているので省略する.\n",
"\n",
"$$\n",
"\\int_{-\\infty}^{\\infty} \\exp(-a x^2) ~\\mathrm{d}x = \\sqrt{\\frac{\\pi}{a}}\n",
"$$\n",
"\n",
"下記のように解析的な積分と数値積分との誤差率が0.01%以内であれば✔を表示する. $a$ の値を5通り変えて検証を行った."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "6ae6047c",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"a | analytical | QuadGK.jl \n",
"------ | ----------------- | -----------------\n",
" 0.01 | 17.724538509055 | 17.724538509055 ✔\n",
" 0.10 | 5.604991216398 | 5.604991216398 ✔\n",
" 1.00 | 1.772453850906 | 1.772453850906 ✔\n",
" 10.00 | 0.560499121640 | 0.560499121640 ✔\n",
"100.00 | 0.177245385091 | 0.177245385091 ✔\n"
]
}
],
"source": [
"println(\"a | analytical | QuadGK.jl \")\n",
"println(\"------ | ----------------- | -----------------\")\n",
"for a in [0.01, 0.1, 1.0, 10.0, 100.0]\n",
" analytical = sqrt(π/a)\n",
" numerical = quadgk(x -> exp(-a*x^2), -Inf, Inf, rtol=1e-12, maxevals=10^7)[1]\n",
" @printf(\"%6.2f | %17.12f | %17.12f %s\\n\", a, analytical, numerical, isapprox(analytical, numerical, rtol=1e-5) ? \"✔\" : \"✗\")\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "bf9a61fd",
"metadata": {},
"source": [
"## 一般化されたガウス積分\n",
"\n",
"1次元において, ガウス積分を下記のように一般化する.\n",
"\n",
"$$\n",
"\\mathrm{GGI}(n,a)\n",
"= \\int_0^{\\infty} x^{n} \\exp \\left(-a x^2\\right)\n",
"$$\n",
"\n",
"ここでは通常のガウス積分は$n=0$の特殊例であると考える. 森口繁一, 宇田川銈久, 一松信『岩波 数学公式Ⅰ』(岩波書店,1988,新装第3刷) p.233では$n$が偶数の場合と奇数の場合で分けた公式\n",
"\n",
"$$\n",
"\\int_0^{\\infty} x^{2n} \\exp \\left(-a x^2\\right)\n",
"= \\frac{(2n-1)!!}{2^{n+1}} \\sqrt{\\frac{\\pi}{a^{2n+1}}}\n",
"$$\n",
"\n",
"$$\n",
"\\int_0^{\\infty} x^{2n+1} \\exp \\left(-a x^2\\right)\n",
"= \\frac{n!}{2 a^{n+1}}\n",
"$$\n",
"\n",
"が与えられているので$n! = \\Gamma\\left( n+1 \\right)$および$\\Gamma\\left(\\frac{1}{2}+n\\right)=\\frac{(2 n-1) ! !}{2^n} \\sqrt{\\pi}$を用いてまとめると\n",
"\n",
"$$\n",
"\\int_0^{\\infty} x^{n} \\exp \\left(-a x^2\\right)\n",
"= \\frac{\\Gamma\\left( \\frac{n+1}{2} \\right)}{2 a^{\\frac{n+1}{2}}}\n",
"$$\n",
"\n",
"が得られる. [英語版のWikipedia](https://en.wikipedia.org/wiki/Gaussian_integral#Integrals_of_similar_form)に載っているのを見つけた. また $a$ が具体的な数値の場合なら[Wolfram Alpha](https://ja.wolframalpha.com/input?i=%E2%88%AB%5B0%2C%E2%88%9E%5D+x%5En+exp%28-7x%5E2%29+dx)でも計算できるので適当な素数を代入すれば公式として取り出せる."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "e759c1b0",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"a | n | analytical | QuadGK.jl \n",
"------ | --- | ----------------- | -----------------\n",
" 0.50 | 0 | 1.253314137316 | 1.253314137316 ✔\n",
" 0.50 | 1 | 1.000000000000 | 1.000000000000 ✔\n",
" 0.50 | 2 | 1.253314137316 | 1.253314137316 ✔\n",
" 0.50 | 3 | 2.000000000000 | 2.000000000000 ✔\n",
" 0.50 | 4 | 3.759942411947 | 3.759942411947 ✔\n",
" 0.50 | 5 | 8.000000000000 | 8.000000000000 ✔\n",
" 0.50 | 6 | 18.799712059733 | 18.799712059733 ✔\n",
" 0.50 | 7 | 48.000000000000 | 48.000000000002 ✔\n",
" 0.50 | 8 | 131.597984418128 | 131.597984418140 ✔\n",
" 0.50 | 9 | 384.000000000000 | 384.000000000000 ✔\n",
" 0.50 | 10 | 1184.381859763148 | 1184.381859763148 ✔\n",
" 1.00 | 0 | 0.886226925453 | 0.886226925453 ✔\n",
" 1.00 | 1 | 0.500000000000 | 0.500000000000 ✔\n",
" 1.00 | 2 | 0.443113462726 | 0.443113462726 ✔\n",
" 1.00 | 3 | 0.500000000000 | 0.500000000000 ✔\n",
" 1.00 | 4 | 0.664670194090 | 0.664670194090 ✔\n",
" 1.00 | 5 | 1.000000000000 | 1.000000000000 ✔\n",
" 1.00 | 6 | 1.661675485224 | 1.661675485224 ✔\n",
" 1.00 | 7 | 3.000000000000 | 3.000000000000 ✔\n",
" 1.00 | 8 | 5.815864198284 | 5.815864198284 ✔\n",
" 1.00 | 9 | 12.000000000000 | 12.000000000000 ✔\n",
" 1.00 | 10 | 26.171388892277 | 26.171388892277 ✔\n",
" 10.00 | 0 | 0.280249560820 | 0.280249560820 ✔\n",
" 10.00 | 1 | 0.050000000000 | 0.050000000000 ✔\n",
" 10.00 | 2 | 0.014012478041 | 0.014012478041 ✔\n",
" 10.00 | 3 | 0.005000000000 | 0.005000000000 ✔\n",
" 10.00 | 4 | 0.002101871706 | 0.002101871706 ✔\n",
" 10.00 | 5 | 0.001000000000 | 0.001000000000 ✔\n",
" 10.00 | 6 | 0.000525467927 | 0.000525467927 ✔\n",
" 10.00 | 7 | 0.000300000000 | 0.000300000000 ✔\n",
" 10.00 | 8 | 0.000183913774 | 0.000183913774 ✔\n",
" 10.00 | 9 | 0.000120000000 | 0.000120000000 ✔\n",
" 10.00 | 10 | 0.000082761198 | 0.000082761198 ✔\n",
"100.00 | 0 | 0.088622692545 | 0.088622692545 ✔\n",
"100.00 | 1 | 0.005000000000 | 0.005000000000 ✔\n",
"100.00 | 2 | 0.000443113463 | 0.000443113463 ✔\n",
"100.00 | 3 | 0.000050000000 | 0.000050000000 ✔\n",
"100.00 | 4 | 0.000006646702 | 0.000006646702 ✔\n",
"100.00 | 5 | 0.000001000000 | 0.000001000000 ✔\n",
"100.00 | 6 | 0.000000166168 | 0.000000166168 ✔\n",
"100.00 | 7 | 0.000000030000 | 0.000000030000 ✔\n",
"100.00 | 8 | 0.000000005816 | 0.000000005816 ✔\n",
"100.00 | 9 | 0.000000001200 | 0.000000001200 ✔\n",
"100.00 | 10 | 0.000000000262 | 0.000000000262 ✔\n"
]
}
],
"source": [
"# using SpecialFunctions # for gamma()\n",
"\n",
"GGI(n,a) = gamma((n+1)/2) / (2*a^((n+1)/2))\n",
"\n",
"println(\"a | n | analytical | QuadGK.jl \")\n",
"println(\"------ | --- | ----------------- | -----------------\")\n",
"for a in [0.5, 1.0, 10.0, 100.0]\n",
" for n in 0:10\n",
" analytical = GGI(n,a)\n",
" numerical = quadgk(x -> x^n * exp(-a*x^2), 0, Inf, maxevals=10^3)[1]\n",
" @printf(\"%6.2f | %3d | %17.12f | %17.12f %s\\n\", a, n, analytical, numerical, isapprox(analytical, numerical, rtol=1e-5) ? \"✔\" : \"✗\")\n",
" end\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "e5746911",
"metadata": {},
"source": [
"## もっと一般化されたガウス積分\n",
"\n",
"More Generalized Gaussian Integralを次のように定義する.\n",
"\n",
"\\begin{aligned}\n",
"\\mathrm{MGGI}(n,a,b)\n",
"&= \\int_0^{\\infty} x^{n} \\exp \\left(-a x -b x^2\\right) \\mathrm{d}x \\\\\n",
"&= \\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)\n",
"\\end{aligned}\n",
"\n",
"ここで $U$ は[合流型超幾何関数](https://en.wikipedia.org/wiki/Confluent_hypergeometric_function)である. これも $a$, $b$ が具体的な数値の場合なら[Wolfram Alpha](https://ja.wolframalpha.com/input?i=%E2%88%AB%5B0%2C%E2%88%9E%5D+x%5En+exp%28+-+7x+-+11x%5E2%29+dx)でも計算でき, 適当な素数を代入すれば公式として取り出せる."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "1af572d0",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"a | b | n | analytical | QuadGK.jl \n",
"------ | ------ | --- | ----------------- | -----------------\n",
" 0.10 | 0.10 | 0 | 2.365023857063 | 2.365023857063 ✔\n",
" 0.10 | 0.10 | 1 | 3.817488071469 | 3.817488071469 ✔\n",
" 0.10 | 0.10 | 2 | 9.916375249581 | 9.916375249581 ✔\n",
" 0.10 | 0.10 | 3 | 33.216693089895 | 33.216693089895 ✔\n",
" 0.10 | 0.10 | 4 | 132.137282198762 | 132.137282198762 ✔\n",
" 0.10 | 0.10 | 5 | 598.265220698515 | 598.265220698516 ✔\n",
" 0.10 | 0.10 | 6 | 3004.299444619790 | 3004.299444619781 ✔\n",
" 0.10 | 0.10 | 7 | 16445.806898645555 | 16445.806898645282 ✔\n",
" 0.10 | 0.10 | 8 | 96927.577112369851 | 96927.577112362909 ✔\n",
" 0.10 | 0.10 | 9 | 609368.487389638554 | 609368.487389637972 ✔\n",
" 0.10 | 0.10 | 10 | 4057056.726361821406 | 4057056.726361818612 ✔\n",
" 0.10 | 1.00 | 0 | 0.838361847809 | 0.838361847809 ✔\n",
" 0.10 | 1.00 | 1 | 0.458081907610 | 0.458081907610 ✔\n",
" 0.10 | 1.00 | 2 | 0.396276828524 | 0.396276828524 ✔\n",
" 0.10 | 1.00 | 3 | 0.438268066183 | 0.438268066183 ✔\n",
" 0.10 | 1.00 | 4 | 0.572501839477 | 0.572501839477 ✔\n",
" 0.10 | 1.00 | 5 | 0.847911040393 | 0.847911040393 ✔\n",
" 0.10 | 1.00 | 6 | 1.388859046672 | 1.388859046672 ✔\n",
" 0.10 | 1.00 | 7 | 2.474290168845 | 2.474290168845 ✔\n",
" 0.10 | 1.00 | 8 | 4.737292154909 | 4.737292154909 ✔\n",
" 0.10 | 1.00 | 9 | 9.660296067635 | 9.660296067635 ✔\n",
" 0.10 | 1.00 | 10 | 20.834799893709 | 20.834799893709 ✔\n",
" 0.10 | 10.00 | 0 | 0.275318798552 | 0.275318798552 ✔\n",
" 0.10 | 10.00 | 1 | 0.048623406007 | 0.048623406007 ✔\n",
" 0.10 | 10.00 | 2 | 0.013522822898 | 0.013522822898 ✔\n",
" 0.10 | 10.00 | 3 | 0.004794726486 | 0.004794726486 ✔\n",
" 0.10 | 10.00 | 4 | 0.002004449802 | 0.002004449802 ✔\n",
" 0.10 | 10.00 | 5 | 0.000948923048 | 0.000948923048 ✔\n",
" 0.10 | 10.00 | 6 | 0.000496367835 | 0.000496367835 ✔\n",
" 0.10 | 10.00 | 7 | 0.000282195075 | 0.000282195075 ✔\n",
" 0.10 | 10.00 | 8 | 0.000172317767 | 0.000172317767 ✔\n",
" 0.10 | 10.00 | 9 | 0.000112016441 | 0.000112016441 ✔\n",
" 0.10 | 10.00 | 10 | 0.000076982913 | 0.000076982913 ✔\n",
" 1.00 | 0.10 | 0 | 0.865392586515 | 0.865392586515 ✔\n",
" 1.00 | 0.10 | 1 | 0.673037067424 | 0.673037067424 ✔\n",
" 1.00 | 0.10 | 2 | 0.961777595453 | 0.961777595453 ✔\n",
" 1.00 | 0.10 | 3 | 1.921482696980 | 1.921482696980 ✔\n",
" 1.00 | 0.10 | 4 | 4.819250446898 | 4.819250446898 ✔\n",
" 1.00 | 0.10 | 5 | 14.333401705099 | 14.333401705099 ✔\n",
" 1.00 | 0.10 | 6 | 48.814252646962 | 48.814252646962 ✔\n",
" 1.00 | 0.10 | 7 | 185.930787918163 | 185.930787918163 ✔\n",
" 1.00 | 0.10 | 8 | 778.844903052866 | 778.844903052855 ✔\n",
" 1.00 | 0.10 | 9 | 3543.007001462193 | 3543.007001462197 ✔\n",
" 1.00 | 0.10 | 10 | 17332.985630067957 | 17332.985630067997 ✔\n",
" 1.00 | 1.00 | 0 | 0.545641360765 | 0.545641360765 ✔\n",
" 1.00 | 1.00 | 1 | 0.227179319617 | 0.227179319617 ✔\n",
" 1.00 | 1.00 | 2 | 0.159231020574 | 0.159231020574 ✔\n",
" 1.00 | 1.00 | 3 | 0.147563809331 | 0.147563809331 ✔\n",
" 1.00 | 1.00 | 4 | 0.165064626195 | 0.165064626195 ✔\n",
" 1.00 | 1.00 | 5 | 0.212595305563 | 0.212595305563 ✔\n",
" 1.00 | 1.00 | 6 | 0.306363912707 | 0.306363912707 ✔\n",
" 1.00 | 1.00 | 7 | 0.484603960337 | 0.484603960337 ✔\n",
" 1.00 | 1.00 | 8 | 0.829971714305 | 0.829971714305 ✔\n",
" 1.00 | 1.00 | 9 | 1.523429984196 | 1.523429984196 ✔\n",
" 1.00 | 1.00 | 10 | 2.973157722275 | 2.973157722275 ✔\n",
" 1.00 | 10.00 | 0 | 0.236502385706 | 0.236502385706 ✔\n",
" 1.00 | 10.00 | 1 | 0.038174880715 | 0.038174880715 ✔\n",
" 1.00 | 10.00 | 2 | 0.009916375250 | 0.009916375250 ✔\n",
" 1.00 | 10.00 | 3 | 0.003321669309 | 0.003321669309 ✔\n",
" 1.00 | 10.00 | 4 | 0.001321372822 | 0.001321372822 ✔\n",
" 1.00 | 10.00 | 5 | 0.000598265221 | 0.000598265221 ✔\n",
" 1.00 | 10.00 | 6 | 0.000300429944 | 0.000300429944 ✔\n",
" 1.00 | 10.00 | 7 | 0.000164458069 | 0.000164458069 ✔\n",
" 1.00 | 10.00 | 8 | 0.000096927577 | 0.000096927577 ✔\n",
" 1.00 | 10.00 | 9 | 0.000060936849 | 0.000060936849 ✔\n",
" 1.00 | 10.00 | 10 | 0.000040570567 | 0.000040570567 ✔\n",
" 10.00 | 0.10 | 0 | 0.099801188165 | 0.099801188165 ✔\n",
" 10.00 | 0.10 | 1 | 0.009940591748 | 0.009940591748 ✔\n",
" 10.00 | 0.10 | 2 | 0.001976353427 | 0.001976353427 ✔\n",
" 10.00 | 0.10 | 3 | 0.000588246113 | 0.000588246113 ✔\n",
" 10.00 | 0.10 | 4 | 0.000232995745 | 0.000232995745 ✔\n",
" 10.00 | 0.10 | 5 | 0.000115135036 | 0.000115135036 ✔\n",
" 10.00 | 0.10 | 6 | 0.000068141803 | 0.000068141803 ✔\n",
" 10.00 | 0.10 | 7 | 0.000046960922 | 0.000046960922 ✔\n",
" 10.00 | 0.10 | 8 | 0.000036916993 | 0.000036916993 ✔\n",
" 10.00 | 0.10 | 9 | 0.000032587264 | 0.000032587264 ✔\n",
" 10.00 | 0.10 | 10 | 0.000031901484 | 0.000031901484 ✔\n",
" 10.00 | 1.00 | 0 | 0.098109430732 | 0.098109430732 ✔\n",
" 10.00 | 1.00 | 1 | 0.009452846342 | 0.009452846342 ✔\n",
" 10.00 | 1.00 | 2 | 0.001790483654 | 0.001790483654 ✔\n",
" 10.00 | 1.00 | 3 | 0.000500428071 | 0.000500428071 ✔\n",
" 10.00 | 1.00 | 4 | 0.000183585126 | 0.000183585126 ✔\n",
" 10.00 | 1.00 | 5 | 0.000082930513 | 0.000082930513 ✔\n",
" 10.00 | 1.00 | 6 | 0.000044310249 | 0.000044310249 ✔\n",
" 10.00 | 1.00 | 7 | 0.000027240296 | 0.000027240296 ✔\n",
" 10.00 | 1.00 | 8 | 0.000018884388 | 0.000018884388 ✔\n",
" 10.00 | 1.00 | 9 | 0.000014539246 | 0.000014539246 ✔\n",
" 10.00 | 1.00 | 10 | 0.000012283518 | 0.000012283518 ✔\n",
" 10.00 | 10.00 | 0 | 0.086539258652 | 0.086539258652 ✔\n",
" 10.00 | 10.00 | 1 | 0.006730370674 | 0.006730370674 ✔\n",
" 10.00 | 10.00 | 2 | 0.000961777595 | 0.000961777595 ✔\n",
" 10.00 | 10.00 | 3 | 0.000192148270 | 0.000192148270 ✔\n",
" 10.00 | 10.00 | 4 | 0.000048192504 | 0.000048192504 ✔\n",
" 10.00 | 10.00 | 5 | 0.000014333402 | 0.000014333402 ✔\n",
" 10.00 | 10.00 | 6 | 0.000004881425 | 0.000004881425 ✔\n",
" 10.00 | 10.00 | 7 | 0.000001859308 | 0.000001859308 ✔\n",
" 10.00 | 10.00 | 8 | 0.000000778845 | 0.000000778845 ✔\n",
" 10.00 | 10.00 | 9 | 0.000000354301 | 0.000000354301 ✔\n",
" 10.00 | 10.00 | 10 | 0.000000173330 | 0.000000173330 ✔\n"
]
}
],
"source": [
"# using SpecialFunctions # for gamma()\n",
"# using GSL # for sf_hyperg_U()\n",
"\n",
"# 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)\n",
"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)\n",
"\n",
"println(\"a | b | n | analytical | QuadGK.jl \")\n",
"println(\"------ | ------ | --- | ----------------- | -----------------\")\n",
"for a in [0.1, 1.0, 10.0]\n",
" for b in [0.1, 1.0, 10.0]\n",
" for n in 0:10\n",
" analytical = MGGI(n,a,b)\n",
" numerical = quadgk(x -> x^n * exp(-a*x - b*x^2), 0, Inf, maxevals=10^3)[1]\n",
" @printf(\"%6.2f | %6.2f | %3d | %17.12f | %17.12f %s\\n\", a, b, n, analytical, numerical, isapprox(analytical, numerical, rtol=1e-5) ? \"✔\" : \"✗\")\n",
" end\n",
" end\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "0a3ecd4e",
"metadata": {},
"source": [
"## 2次元のガウス積分\n",
"\n",
"以下では, 上記の2つの一般化とは異なり, 次元を増やすという意味での一般化を考える. 2次元の積分は逐次的に$x,y$の順番で積分を計算していけば, 単に1次元のガウス積分の2乗になることがわかる. $x$による積分の部分は$y$による積分から見れば定数なので, くくり出せる. なお1変数のガウス積分の証明で2変数の積分を使っているのだから循環論法になっている気もするが, 他にも証明方法があったはずなので気にしない. \n",
"\n",
"\\begin{aligned}\n",
"\\iint_{\\mathbb{R}^2}\n",
"\\exp(-a r^2) ~\n",
"\\mathrm{d}\\pmb{r}\n",
"&=\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a x^2 -a y^2) ~\n",
"\\mathrm{d}x\n",
"\\mathrm{d}y \\\\\n",
"&=\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a x^2) \n",
"\\exp(-a y^2) ~\n",
"\\mathrm{d}x\n",
"\\mathrm{d}y \\\\\n",
"&=\n",
"\\left[\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a x^2) \n",
"\\mathrm{d}x\n",
"\\right]\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a y^2)\n",
"\\mathrm{d}y \\\\\n",
"&=\n",
"\\sqrt{\\frac{\\pi}{a}}\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a y^2)\n",
"\\mathrm{d}y \\\\\n",
"&=\n",
"\\sqrt{\\frac{\\pi}{a}}^2\n",
"\\end{aligned}"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "899bb743",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"a | analytical | QuadGK.jl \n",
"------ | ----------------- | -----------------\n",
" 0.01 | 314.159265358979 | 314.159265358887 ✔\n",
" 0.10 | 31.415926535898 | 31.415926535898 ✔\n",
" 1.00 | 3.141592653590 | 3.141592653590 ✔\n",
" 10.00 | 0.314159265359 | 0.314159265359 ✔\n",
"100.00 | 0.031415926536 | 0.031415926536 ✔\n"
]
}
],
"source": [
"println(\"a | analytical | QuadGK.jl \")\n",
"println(\"------ | ----------------- | -----------------\")\n",
"for a in [0.01, 0.1, 1.0, 10.0, 100.0]\n",
" analytical = sqrt(π/a)^2\n",
" numerical = (\n",
" quadgk(y ->\n",
" quadgk(x ->\n",
" exp(-a*x^2 -a*y^2)\n",
" , -Inf, Inf, maxevals=10^3)[1]\n",
" , -Inf, Inf, maxevals=10^3)[1]\n",
" )\n",
" @printf(\"%6.2f | %17.12f | %17.12f %s\\n\", a, analytical, numerical, isapprox(analytical, numerical, rtol=1e-5) ? \"✔\" : \"✗\")\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "7244dc31",
"metadata": {},
"source": [
"## 極座標系におけるガウス積分\n",
"\n",
"次式で与えられる極座標系への変換$(r,\\theta)\\mapsto(x,y)$の下でガウス積分を計算する.\n",
"\n",
"\\begin{aligned}\n",
" \\left\\{\n",
" \\begin{aligned}\n",
" x =& r \\cos \\theta \\\\\n",
" y =& r \\sin \\theta \\\\\n",
" \\end{aligned}\n",
" \\right.\n",
"\\end{aligned}\n",
"\n",
"ヤコビアンは下記のように計算できる.\n",
"\n",
"\\begin{aligned}\n",
" \\mathrm{d}x\n",
" \\mathrm{d}y\n",
" &=\n",
" \\left|\\frac{\\partial(x,y)}{\\partial(r,\\theta)}\\right|\n",
" \\mathrm{d}r\n",
" \\mathrm{d}\\theta \\\\\n",
" &=\n",
" \\left|\\det\\left(\\begin{array}{rr}\n",
" \\frac{\\partial x}{\\partial r} &\n",
" \\frac{\\partial x}{\\partial \\theta} \\\\\n",
" \\frac{\\partial y}{\\partial r} &\n",
" \\frac{\\partial y}{\\partial \\theta} \\\\\n",
" \\end{array}\\right)\\right|\n",
" \\mathrm{d}r\n",
" \\mathrm{d}\\theta \\\\\n",
" &=\n",
" \\left|\\det\\left(\\begin{array}{rr}\n",
" \\cos\\theta &\n",
" -r\\sin\\theta \\\\\n",
" \\sin\\theta &\n",
" r\\cos\\theta \\\\\n",
" \\end{array}\\right)\\right|\n",
" \\mathrm{d}r\n",
" \\mathrm{d}\\theta \\\\\n",
" &=\n",
" \\left|\n",
" \\cos\\theta \\times r\\cos\\theta - \\sin\\theta \\times (-r\\sin\\theta)\n",
" \\right|\n",
" \\mathrm{d}r\n",
" \\mathrm{d}\\theta \\\\\n",
" &=\n",
" \\left|\n",
" r\\cos^2\\theta + r\\sin^2\\theta\n",
" \\right|\n",
" \\mathrm{d}r\n",
" \\mathrm{d}\\theta \\\\\n",
" &=\n",
" \\left|\n",
" r(\\cos^2\\theta + \\sin^2\\theta)\n",
" \\right|\n",
" \\mathrm{d}r\n",
" \\mathrm{d}\\theta \\\\\n",
" &=\n",
" \\left|\n",
" r\n",
" \\right|\n",
" \\mathrm{d}r\n",
" \\mathrm{d}\\theta \\\\\n",
" &=\n",
" r\n",
" \\mathrm{d}r\n",
" \\mathrm{d}\\theta \\\\\n",
"\\end{aligned}\n",
"\n",
"なお, $0\\leq r<\\infty$であるから$|r|=r$としている. これを用いて, 先ほどと同様に1次元のガウス積分の2乗になることが示せる.\n",
"\n",
"\\begin{aligned}\n",
" \\iint_{\\mathbb{R}^2}\n",
" \\exp(-a r^2) ~\n",
" \\mathrm{d}\\pmb{r}\n",
" &=\n",
" \\int_{0}^{2\\pi}\n",
" \\int_{0}^{\\infty}\n",
" \\exp(-a r^2) ~\n",
" r ~\n",
" \\mathrm{d}r\n",
" \\mathrm{d}\\theta \\\\\n",
" &=\n",
" \\int_{0}^{2\\pi}\n",
" \\mathrm{d}\\theta\n",
" \\int_{0}^{\\infty}\n",
" r \\exp(-a r^2) ~\n",
" \\mathrm{d}r \\\\\n",
" &=\n",
" \\left[ \\theta \\right]_{0}^{2\\pi}\n",
" \\frac{0!}{2 a^{0+1}} \\\\\n",
" &=\n",
" 2\\pi \\frac{1}{2a} \\\\\n",
" &=\n",
" \\sqrt{\\frac{\\pi}{a}}^2 \\\\\n",
"\\end{aligned}"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "980b2347",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"a | analytical | QuadGK.jl \n",
"------ | ----------------- | -----------------\n",
" 0.01 | 314.159265358979 | 314.159265358979 ✔\n",
" 0.10 | 31.415926535898 | 31.415926535898 ✔\n",
" 1.00 | 3.141592653590 | 3.141592653590 ✔\n",
" 10.00 | 0.314159265359 | 0.314159265359 ✔\n",
"100.00 | 0.031415926536 | 0.031415926536 ✔\n"
]
}
],
"source": [
"println(\"a | analytical | QuadGK.jl \")\n",
"println(\"------ | ----------------- | -----------------\")\n",
"for a in [0.01, 0.1, 1.0, 10.0, 100.0]\n",
" analytical = sqrt(π/a)^2\n",
" numerical = (\n",
" quadgk(θ ->\n",
" quadgk(r ->\n",
" r * exp(-a*r^2)\n",
" , 0, Inf, maxevals=10^3)[1]\n",
" , 0, 2*π, maxevals=10^3)[1]\n",
" )\n",
" @printf(\"%6.2f | %17.12f | %17.12f %s\\n\", a, analytical, numerical, isapprox(analytical, numerical, rtol=1e-5) ? \"✔\" : \"✗\")\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "41d79464",
"metadata": {},
"source": [
"## 3次元のガウス積分\n",
"\n",
"3次元でも2次元と同様にして, 単に1次元のガウス積分の3乗となる.\n",
"\n",
"\\begin{aligned}\n",
"\\iiint_{\\mathbb{R}^3}\n",
"\\exp(-a r^2) ~\n",
"\\mathrm{d}\\pmb{r}\n",
"&=\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a x^2 -a y^2 -a z^2) ~\n",
"\\mathrm{d}x\n",
"\\mathrm{d}y\n",
"\\mathrm{d}z \\\\\n",
"&=\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a x^2) \n",
"\\exp(-a y^2) \n",
"\\exp(-a z^2) ~\n",
"\\mathrm{d}x\n",
"\\mathrm{d}y\n",
"\\mathrm{d}z \\\\\n",
"&=\n",
"\\left[\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a x^2) \n",
"\\mathrm{d}x\n",
"\\right] ~\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a y^2) \n",
"\\exp(-a z^2)\n",
"\\mathrm{d}y\n",
"\\mathrm{d}z \\\\\n",
"&=\n",
"\\sqrt{\\frac{\\pi}{a}}\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a y^2) \n",
"\\exp(-a z^2)\n",
"\\mathrm{d}y\n",
"\\mathrm{d}z \\\\\n",
"&=\n",
"\\sqrt{\\frac{\\pi}{a}}\n",
"\\left[\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a y^2) \n",
"\\mathrm{d}y\n",
"\\right]\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a z^2)\n",
"\\mathrm{d}z \\\\\n",
"&=\n",
"\\sqrt{\\frac{\\pi}{a}}^2\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a z^2)\n",
"\\mathrm{d}z \\\\\n",
"&=\n",
"\\sqrt{\\frac{\\pi}{a}}^3\n",
"\\end{aligned}"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "bea6f541",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"a | analytical | QuadGK.jl \n",
"------ | ----------------- | -----------------\n",
" 0.01 | 5568.327996831709 | 5568.327996829265 ✔\n",
" 0.10 | 176.085992288711 | 176.085992288711 ✔\n",
" 1.00 | 5.568327996832 | 5.568327996832 ✔\n",
" 10.00 | 0.176085992289 | 0.176085992289 ✔\n",
"100.00 | 0.005568327997 | 0.005568327997 ✔\n"
]
}
],
"source": [
"println(\"a | analytical | QuadGK.jl \")\n",
"println(\"------ | ----------------- | -----------------\")\n",
"for a in [0.01, 0.1, 1.0, 10.0, 100.0]\n",
" analytical = sqrt(π/a)^3\n",
" numerical = (\n",
" quadgk(z ->\n",
" quadgk(y ->\n",
" quadgk(x ->\n",
" exp(-a*x^2 -a*y^2 -a*z^2)\n",
" , -Inf, Inf, maxevals=10^3)[1]\n",
" , -Inf, Inf, maxevals=10^3)[1]\n",
" , -Inf, Inf, maxevals=10^3)[1]\n",
" )\n",
" @printf(\"%6.2f | %17.12f | %17.12f %s\\n\", a, analytical, numerical, isapprox(analytical, numerical, rtol=1e-5) ? \"✔\" : \"✗\")\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "a6323d0d",
"metadata": {},
"source": [
"## 球面座標系におけるガウス積分\n",
"\n",
"球面座標系への変換\n",
"\n",
"$$\n",
"\\begin{cases}\n",
" x = r\\sin\\theta\\cos\\varphi\\\\\n",
" y = r\\sin\\theta\\sin\\varphi\\\\\n",
" z = r\\cos\\theta\n",
"\\end{cases}\n",
"$$\n",
"\n",
"におけるヤコビアンは\n",
"\n",
"\\begin{aligned}\n",
" \\mathrm{d}x\n",
" \\mathrm{d}y\n",
" \\mathrm{d}z\n",
" &=\n",
" \\left|\\frac{\\partial(x,y,z)}{\\partial(r,\\theta,\\varphi)}\\right|\n",
" \\mathrm{d}r\n",
" \\mathrm{d}\\theta\n",
" \\mathrm{d}\\varphi \\\\\n",
" &=\n",
" \\left|\\det\\left(\\begin{array}{rr}\n",
" \\frac{\\partial x}{\\partial r} &\n",
" \\frac{\\partial x}{\\partial \\theta} &\n",
" \\frac{\\partial x}{\\partial \\varphi} \\\\\n",
" \\frac{\\partial y}{\\partial r} &\n",
" \\frac{\\partial y}{\\partial \\theta} &\n",
" \\frac{\\partial y}{\\partial \\varphi} \\\\\n",
" \\frac{\\partial z}{\\partial r} &\n",
" \\frac{\\partial z}{\\partial \\theta} &\n",
" \\frac{\\partial z}{\\partial \\varphi} \\\\\n",
" \\end{array}\\right)\\right|\n",
" \\mathrm{d}r\n",
" \\mathrm{d}\\theta\n",
" \\mathrm{d}\\varphi \\\\\n",
" &=\n",
" \\cdots \\\\\n",
" &=\n",
" r^2 \\sin(\\theta)\n",
" \\mathrm{d}r\n",
" \\mathrm{d}\\theta\n",
" \\mathrm{d}\\varphi \\\\\n",
"\\end{aligned}\n",
"\n",
"のように計算できるから, これを用いて積分領域$0\\leq r<\\infty, 0\\leq \\theta \\leq \\pi, \\leq \\varphi<2 \\pi$で積分すれば\n",
"\n",
"\\begin{aligned}\n",
" \\iiint_{\\mathbb{R}^3}\n",
" \\exp(-a r^2) ~\n",
" \\mathrm{d}\\pmb{r}\n",
" &=\n",
" \\int_{0}^{2\\pi}\n",
" \\int_{0}^{\\pi}\n",
" \\int_{0}^{\\infty}\n",
" \\exp(-a r^2) ~\n",
" r^2 \\sin\\theta ~\n",
" \\mathrm{d}r\n",
" \\mathrm{d}\\theta\n",
" \\mathrm{d}\\varphi \\\\\n",
" &=\n",
" 4\\pi\n",
" \\int_{0}^{\\infty}\n",
" r^2 \\exp(-a r^2) ~\n",
" \\mathrm{d}r \\\\\n",
" &=\n",
" 4\\pi \\frac{(2\\cdot1-1)!!}{2^{1+1}} \\sqrt{\\frac{\\pi}{a^{2\\cdot1+1}}} \\\\\n",
" &=\n",
" 4\\pi \\frac{1}{4} \\sqrt{\\frac{\\pi}{a^{3}}} \\\\\n",
" &=\n",
" \\sqrt{\\frac{\\pi}{a}}^3\n",
"\\end{aligned}\n",
"\n",
"と計算できる. なお, 計算には1次元のガウス積分を用いている."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "413fa0dd",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"a | analytical | QuadGK.jl \n",
"------ | ----------------- | -----------------\n",
" 0.01 | 5568.327996831709 | 5567.958330629872 ✗\n",
" 0.10 | 176.085992288711 | 176.085992503519 ✔\n",
" 1.00 | 5.568327996832 | 5.568327996978 ✔\n",
" 10.00 | 0.176085992289 | 0.176085992289 ✔\n",
"100.00 | 0.005568327997 | 0.005568327997 ✔\n"
]
}
],
"source": [
"println(\"a | analytical | QuadGK.jl \")\n",
"println(\"------ | ----------------- | -----------------\")\n",
"for a in [0.01, 0.1, 1.0, 10.0, 100.0]\n",
" analytical = sqrt(π/a)^3\n",
" numerical = (\n",
" quadgk(φ ->\n",
" quadgk(θ ->\n",
" quadgk(r ->\n",
" r^2 * sin(θ) * exp(-a*r^2)\n",
" , 0, Inf, maxevals=100)[1]\n",
" , 0, π, maxevals=10)[1]\n",
" , 0, 2*π, maxevals=20)[1]\n",
" )\n",
" @printf(\"%6.2f | %17.12f | %17.12f %s\\n\", a, analytical, numerical, isapprox(analytical, numerical, rtol=1e-5) ? \"✔\" : \"✗\")\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "87c13cdc",
"metadata": {},
"source": [
"最後にHCubature.jlとMonteCarloIntegration.jlも加えた結果を示す. モンテカルロ法では誤差率の閾値は甘めに1%とした."
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "61a76839",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"a | analytical | QuadGK.jl | HCubature.jl | MonteCarloIntegration.jl\n",
"------ | ----------------- | ----------------- | ----------------- | -----------------\n",
" 0.01 | 5568.327996831709 | 5568.327996831707 ✔ | 5568.327997006358 ✔ | 5567.410926491379 ✔ \n",
" 0.10 | 176.085992288711 | 176.085992288711 ✔ | 176.085992294186 ✔ | 176.025177198522 ✔ \n",
" 1.00 | 5.568327996832 | 5.568327996832 ✔ | 5.568327997105 ✔ | 5.564038220671 ✔ \n",
" 10.00 | 0.176085992289 | 0.176085992289 ✔ | 0.176085992297 ✔ | 0.175847087491 ✔ \n",
"100.00 | 0.005568327997 | 0.005568327997 ✔ | 0.005568327997 ✔ | 0.005570149391 ✔ \n"
]
}
],
"source": [
"println(\"a | analytical | QuadGK.jl | HCubature.jl | MonteCarloIntegration.jl\")\n",
"println(\"------ | ----------------- | ----------------- | ----------------- | -----------------\")\n",
"for a in [0.01, 0.1, 1.0, 10.0, 100.0]\n",
" analytical = sqrt(π/a)^3\n",
" I₁ = (\n",
" quadgk(φ ->\n",
" quadgk(θ ->\n",
" quadgk(r ->\n",
" r^2 * sin(θ) * exp(-a*r^2)\n",
" , 0, Inf, maxevals=200)[1]\n",
" , 0, π, maxevals=10)[1]\n",
" , 0, 2*π, maxevals=20)[1]\n",
" )\n",
" I₂ = hcubature(x -> x[1]^2 * sin(x[2]) * exp(-a*x[1]^2), [0,0,0], [100,π,2*π])[1]\n",
" I₃ = vegas(x -> x[1]^2 * sin(x[2]) * exp(-a*x[1]^2), [0,0,0], [100,π,2*π]).integral_estimate\n",
" @printf(\"%6.2f | %17.12f | %17.12f %s | %17.12f %s | %17.12f %s \\n\", a, analytical,\n",
" I₁, isapprox(analytical, I₁, rtol=1e-5) ? \"✔\" : \"✗\",\n",
" I₂, isapprox(analytical, I₂, rtol=1e-5) ? \"✔\" : \"✗\",\n",
" I₃, isapprox(analytical, I₃, rtol=1e-2) ? \"✔\" : \"✗\"\n",
" )\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "dd0472d1",
"metadata": {},
"source": [
"## 参考文献\n",
"\n",
"- [Wikipedia 日本語版](https://ja.wikipedia.org/wiki/%E3%82%AC%E3%82%A6%E3%82%B9%E7%A9%8D%E5%88%86)\n",
"- [Wikipedia 英語版](https://en.wikipedia.org/wiki/Gaussian_integral#Integrals_of_similar_form)\n",
"- 森口繁一, 宇田川銈久, 一松信『岩波 数学公式Ⅰ』(岩波書店,1988,新装第3刷) p.233\n",
"- [QuadGK.jlドキュメント](https://juliamath.github.io/QuadGK.jl/stable/)\n",
"- [HCubature.jlドキュメント](https://github.com/JuliaMath/HCubature.jl)\n",
"- [MonteCarloIntegration.jlドキュメント](https://github.com/ranjanan/MonteCarloIntegration.jl)\n",
"- Wolfram Alpha, [GGI](https://ja.wolframalpha.com/input?i=%E2%88%AB%5B0%2C%E2%88%9E%5D+x%5En+exp%28-7x%5E2%29+dx), [MGGI](https://ja.wolframalpha.com/input?i=%E2%88%AB%5B0%2C%E2%88%9E%5D+x%5En+exp%28+-+7x+-+11x%5E2%29+dx)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.8.3",
"language": "julia",
"name": "julia-1.8"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment