Skip to content

Instantly share code, notes, and snippets.

@ohno
Created July 12, 2023 03:22
Show Gist options
  • Save ohno/18b7014012bef5ad7c21acbb474e29a1 to your computer and use it in GitHub Desktop.
Save ohno/18b7014012bef5ad7c21acbb474e29a1 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "c5092c41",
"metadata": {},
"outputs": [],
"source": [
"using QuadGK"
]
},
{
"cell_type": "markdown",
"id": "2770d7f6",
"metadata": {},
"source": [
"\\begin{align}\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",
"\\int_{-\\infty}^{\\infty}\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a y^2) \n",
"\\exp(-a z^2)\n",
"\\left[\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a x^2) \n",
"\\mathrm{d}x\n",
"\\right] ~\n",
"\\mathrm{d}y\n",
"\\mathrm{d}z \\\\\n",
"&=\n",
"\\sqrt{\\frac{\\pi}{a}}\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a z^2)\n",
"\\left[\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a y^2) \n",
"\\mathrm{d}y\n",
"\\right]\n",
"\\mathrm{d}z \\\\\n",
"&=\n",
"\\sqrt{\\frac{\\pi}{a}}\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a z^2)\n",
"\\left[\n",
"\\int_{-\\infty}^{\\infty}\n",
"\\exp(-a y^2) \n",
"\\mathrm{d}y\n",
"\\right]\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{align}"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "dc350e2c",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"evals \tnumerical\n",
"10^3\t1.001530395550591\n",
"20^3\t0.9999875777006089\n",
"50^3\t0.9999938466553702\n",
"100^3\t1.0000001156363314\n",
"200^3\t1.0000000000003697\n",
"500^3\t1.0000000000003697\n",
"1000^3\t1.0000000000003697\n",
"∞\t1.0\n"
]
}
],
"source": [
"a = π\n",
"\n",
"println(\"evals \\tnumerical\")\n",
"for n in [10, 20, 50, 100, 200, 500, 1000]\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=n)[1]\n",
" , -Inf, Inf, maxevals=n)[1]\n",
" , -Inf, Inf, maxevals=n)[1]\n",
" )\n",
" println(\"$n^3\\t\", numerical)\n",
"end\n",
"println(\"∞\\t\", sqrt(π/a)^3)"
]
}
],
"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