Editing Gaussian integral (section)

==Computation==

===By polar coordinates===
A standard way to compute the Gaussian integral, the idea of which goes back to Poisson,<ref name="york.ac.uk">{{cite web |title=The Probability Integral | last=Lee | first=Peter M. |url=https://www.york.ac.uk/depts/maths/histstat/normal_history.pdf }}</ref> is to make use of the property that:

<math display="block">\left(\int_{-\infty}^{\infty} e^{-x^2}\,dx\right)^2 = \int_{-\infty}^{\infty} e^{-x^2}\,dx \int_{-\infty}^{\infty} e^{-y^2}\,dy = \int_{-\infty}^{\infty}  \int_{-\infty}^{\infty} e^{-\left(x^2+y^2\right)}\, dx\,dy. </math>

Consider the function <math>e^{-\left(x^2 + y^2\right)} = e^{-r^{2}}</math>on the plane <math>\mathbb{R}^2</math>, and compute its integral two ways:
# on the one hand, by [[double integration]] in the [[Cartesian coordinate system]], its integral is a square: <math display="block">\left(\int e^{-x^2}\,dx\right)^2;</math>
# on the other hand, by [[shell integration]] (a case of double integration in [[polar coordinates]]), its integral is computed to be <math>\pi</math>

Comparing these two computations yields the integral, though one should take care about the [[improper integral]]s involved.

<math display="block">\begin{align}
   \iint_{\R^2} e^{-\left(x^2 + y^2\right)}dx\,dy
   &= \int_0^{2\pi} \int_0^{\infty} e^{-r^2}r\,dr\,d\theta\\[6pt]
   &= 2\pi \int_0^\infty re^{-r^2}\,dr\\[6pt]
   &= 2\pi \int_{-\infty}^0 \tfrac{1}{2} e^s\,ds && s = -r^2\\[6pt]
   &= \pi \int_{-\infty}^0 e^s\,ds \\[6pt]
   &= \pi \, \left[ e^s\right]_{-\infty}^{0} \\[6pt]
   &= \pi \,\left(e^0 - e^{-\infty}\right) \\[6pt]
   &= \pi \,\left(1 - 0\right) \\[6pt]
   &=\pi,
 \end{align}</math>
where the factor of {{mvar|r}} is the [[Jacobian determinant]] which appears because of the [[list of canonical coordinate transformations|transform to polar coordinates]] ({{math|''r'' ''dr'' ''dθ''}} is the standard measure on the plane, expressed in polar coordinates [[Wikibooks:Calculus/Polar Integration#Generalization]]), and the substitution involves taking {{math|1=''s'' = −''r''<sup>2</sup>}}, so {{math|1=''ds'' = −2''r'' ''dr''}}.

Combining these yields
<math display="block">\left ( \int_{-\infty}^\infty e^{-x^2}\,dx \right )^2=\pi,</math>
so
<math display="block">\int_{-\infty}^\infty e^{-x^2} \, dx = \sqrt{\pi}.</math>

====Complete proof====
To justify the improper double integrals and equating the two expressions, we begin with an approximating function:
<math display="block">I(a) = \int_{-a}^a e^{-x^2}dx.</math>

If the integral
<math display="block">\int_{-\infty}^\infty e^{-x^2} \, dx</math>
were [[absolutely convergent]] we would have that its [[Cauchy principal value]], that is, the limit
<math display="block">\lim_{a\to\infty} I(a) </math>
would coincide with
<math display="block">\int_{-\infty}^\infty e^{-x^2}\,dx.</math>
To see that this is the case, consider that

<math display="block">\int_{-\infty}^\infty \left|e^{-x^2}\right| dx < \int_{-\infty}^{-1} -x e^{-x^2}\, dx + \int_{-1}^1 e^{-x^2}\, dx+ \int_{1}^{\infty} x e^{-x^2}\, dx < \infty .</math>

So we can compute
<math display="block">\int_{-\infty}^\infty e^{-x^2} \, dx</math>
by just taking the limit
<math display="block">\lim_{a\to\infty} I(a).</math>

Taking the square of <math>I(a)</math> yields

<math display="block">\begin{align}
I(a)^2 & = \left ( \int_{-a}^a e^{-x^2}\, dx \right ) \left ( \int_{-a}^a e^{-y^2}\, dy \right ) \\[6pt]
& = \int_{-a}^a \left ( \int_{-a}^a e^{-y^2}\, dy \right )\,e^{-x^2}\, dx \\[6pt]
& = \int_{-a}^a \int_{-a}^a e^{-\left(x^2+y^2\right)}\,dy\,dx.
\end{align}</math>

Using [[Fubini's theorem]], the above double integral can be seen as an area integral
<math display="block">\iint_{[-a, a] \times [-a, a]} e^{-\left(x^2+y^2\right)}\,d(x,y),</math>
taken over a square with vertices {{math|{(−''a'', ''a''), (''a'', ''a''), (''a'', −''a''), (−''a'', −''a'')}<nowiki/>}} on the ''xy''-[[Cartesian plane|plane]].

Since the exponential function is greater than 0 for all real numbers, it then follows that the integral taken over the square's [[incircle]] must be less than <math>I(a)^2</math>, and similarly the integral taken over the square's [[circumcircle]] must be greater than <math>I(a)^2</math>. The integrals over the two disks can easily be computed by switching from Cartesian coordinates to [[list of canonical coordinate transformations|polar coordinates]]:

<math display="block">\begin{align}
x &= r \cos \theta, &
y &= r \sin\theta
\end{align}</math>
<math display="block">
\mathbf J(r, \theta) = 
\begin{bmatrix}
  \dfrac{\partial x}{\partial r} & \dfrac{\partial x}{\partial\theta}\\[1em]
  \dfrac{\partial y}{\partial r} & \dfrac{\partial y}{\partial\theta} \end{bmatrix}
= \begin{bmatrix}
  \cos\theta &           -  r\sin \theta \\
  \sin\theta & \hphantom{-} r\cos \theta
\end{bmatrix}
</math>
<math display="block">d(x,y) = \left|J(r, \theta)\right| d(r,\theta) = r\, d(r,\theta).</math>
<math display="block">\int_0^{2\pi} \int_0^a re^{-r^2} \, dr \, d\theta < I^2(a) < \int_0^{2\pi} \int_0^{a\sqrt{2}} re^{-r^2} \, dr\, d\theta.</math>

(See [[list of canonical coordinate transformations|to polar coordinates from Cartesian coordinates]] for help with polar transformation.)

Integrating,
<math display="block">\pi \left(1-e^{-a^2}\right) <  I^2(a) < \pi \left(1 - e^{-2a^2}\right). </math>

By the [[squeeze theorem]], this gives the Gaussian integral
<math display="block">\int_{-\infty}^\infty e^{-x^2}\, dx = \sqrt{\pi}.</math>

===By Cartesian coordinates===
A different technique, which goes back to Laplace (1812),<ref name="york.ac.uk" /> is the following. Let
<math display="block">\begin{align}
y & = xs \\
dy & = x\,ds.
\end{align}</math>

Since the limits on {{mvar|s}} as {{math|''y'' → ±∞}} depend on the sign of {{mvar|x}}, it simplifies the calculation to use the fact that {{math|''e''<sup>−''x''<sup>2</sup></sup>}} is an [[even function]], and, therefore, the integral over all real numbers is just twice the integral from zero to infinity. That is,

<math display="block">\int_{-\infty}^{\infty} e^{-x^2} \, dx = 2\int_{0}^{\infty} e^{-x^2}\,dx.</math>

Thus, over the range of integration, {{math|''x'' ≥ 0}}, and the variables {{mvar|y}} and {{mvar|s}} have the same limits. This yields:
<math display="block">\begin{align}
I^2 &= 4 \int_0^\infty \int_0^\infty e^{-\left(x^2 + y^2\right)} dy\,dx \\[6pt]
&= 4 \int_0^\infty \left( \int_0^\infty e^{-\left(x^2 + y^2\right)} \, dy \right) \, dx \\[6pt]
&= 4 \int_0^\infty \left( \int_0^\infty e^{-x^2\left(1+s^2\right)} x\,ds \right) \, dx \\[6pt]
\end{align}</math>
Then, using [[Fubini's theorem]] to switch the [[order of integration (calculus)|order of integration]]:
<math display="block">\begin{align}
I^2 &= 4 \int_0^\infty \left( \int_0^\infty e^{-x^2\left(1 + s^2\right)} x \, dx \right) \, ds  \\[6pt]
&= 4 \int_0^\infty \left[ \frac{e^{-x^2\left(1+s^2\right)} }{-2 \left(1+s^2\right)} \right]_{x=0}^{x=\infty} \, ds \\[6pt]
&= 4 \left (\frac{1}{2} \int_0^\infty \frac{ds}{1+s^2} \right) \\[6pt]
&= 2 \arctan(s)\Big |_0^\infty \\[6pt]
&= \pi.
\end{align}</math>

Therefore, <math>I = \sqrt{\pi}</math>, as expected.

=== By [[Laplace's method]] ===

In Laplace approximation, we deal only with up to second-order terms in Taylor expansion, so we consider <math>e^{-x^2}\approx 1-x^2 \approx (1+x^2)^{-1}</math>.  

In fact, since <math>(1+t)e^{-t} \leq 1</math> for all <math>t</math>, we have the exact bounds:<math display="block">1-x^2 \leq e^{-x^2} \leq (1+x^2)^{-1}</math>Then we can do the bound at Laplace approximation limit:<math display="block">\int_{[-1, 1]}(1-x^2)^n dx \leq \int_{[-1, 1]}e^{-nx^2} dx \leq \int_{[-1, 1]}(1+x^2)^{-n} dx</math>  

That is,
<math display="block">2\sqrt n\int_{[0, 1]}(1-x^2)^n dx \leq \int_{[-\sqrt n, \sqrt n]}e^{-x^2} dx \leq 2\sqrt n\int_{[0, 1]}(1+x^2)^{-n} dx</math>  

By trigonometric substitution, we exactly compute those two bounds: <math>2\sqrt n(2n)!!/(2n+1)!!</math> and <math>2\sqrt n (\pi/2)(2n-3)!!/(2n-2)!!</math>  

By taking the square root of the [[Wallis formula]], <math display="block">\frac \pi 2 = \prod_{n=1} \frac{(2n)^2}{(2n-1)(2n+1)}</math>we have <math>\sqrt \pi = 2 \lim_{n\to \infty} \sqrt{n} \frac{(2n)!!}{(2n+1)!!}</math>, the desired lower bound limit. Similarly we can get the desired upper bound limit.
Conversely, if we first compute the integral with one of the other methods above, we would obtain a proof of the Wallis formula.