To illustrate the power of series expansions and to continue exploring the centrality of Gaussian integrals, we will develop an algebraic approximation to the factorial function \(n!\) for large \(n\) that is good to parts per million already by \(n = 10\).
You are quite familiar with factorials, \[ n! = n \times (n-1) \times \cdot \times 2 \times 1 \] with the understanding that \(0! = 1\). Euler’s \(\Gamma\) function generalizes the factorial function to nonintegral values of the argument. It is defined by \begin{equation}\label{eq:Gamma} \boxed{ \Gamma(n+1) = \int_0^\infty x^n e^{-x} \dd{x} } \end{equation} where the extra 1 in the argument is courtesy of Legendre.
To see that \(\Gamma(n+1)\) corresponds to the factorial function, integrate by parts: \[ \Gamma(n+1) = \underbrace{\left. n x^{n-1} e^{-x} \right|_0^\infty}_{\text{vanishes}} + \int_0^\infty n x^{n-1} e^{-x}\dd{x} = n \Gamma(n) \] or \begin{equation}\label{eq:recurrence} \Gamma(n+1) = n \Gamma(n) \end{equation}
When \(n = 1\), we have a straightforward integral to do \[ \Gamma(1) = \int_0^\infty x^0 e^{-x} \dd{x} = \left. -e^{-x} \right|_0^\infty = 1 \] Hence, the recursion relation \(\Gamma(n+1) = n \Gamma(n)\) along with the termination condition \(\Gamma(1) = 1\) proves that \(\Gamma(n+1) = n!\) for nonnegative integer \(n\).
What about when \(n\) is non-integral? For instance, what about \(n = -\frac12\)? That is, can we evaluate \[ I = \int_0^\infty \frac{1}{\sqrt{t}} e^{-t} \dd{t} \] Let \(x = \sqrt{t}\) or \(x^2 = t\), so that \(\dd{t} = 2 x \dd{x}\), giving \[ I = \int_0^\infty \frac1x e^{-x^2} 2x \dd{x} = 2\int_0^\infty e^{-x^2}\dd{x} = \int_{-\infty}^{\infty} e^{-x^2}\dd{x} \] While it may not look like we are any closer to an evaluation, it is undeniably true that squaring this expression yields \[ I^2 = \int_{-\infty}^\infty \dd{x} \int_{-\infty}^\infty \dd{y} \, e^{-x^2} e^{-y^2} = \int_{-\infty}^\infty \dd{x} \int_{-\infty}^\infty \dd{y} \; e^{-(x^2+y^2)} \] since both \(x\) and \(y\) are dummy variables of integration. On the other hand, we can read this expression as the integral over the \(xy\) plane of the integrand \(e^{-x^2-y^2} = e^{-r^2}\), where \(r\) is the distance from the origin. Rather than integrating in cartesians, we can use polar coordinates: \[ I^2 = \int_0^{2\pi} \int_0^\infty e^{-r^2} \, r \dd{r} \dd{\theta} \] Making the \(u\) substitution \(u = r^2\), so that \(\dd{u} = 2 r \dd{r}\), we can rewrite this double integral as \[ I^2 = \int_0^{2\pi} \int_0^\infty e^{-u} \frac{\dd{u}}{2} \dd{\theta} = \pi \] Therefore, \[ I = \Gamma\qty(\frac12) = \int_{-\infty}^\infty e^{-x^2}\dd{x} = \sqrt{\pi} \] We can then use the recurrence relation of Eq. (\ref{eq:recurrence}) to deduce that \begin{align} \Gamma\qty(\frac32) &= \frac12 \Gamma\qty(\frac12) = \frac{\sqrt{\pi}}{2} \notag \\ \Gamma\qty(\frac52) &= \frac32 \Gamma\qty(\frac32) = \frac{3\sqrt{\pi}}{4} \notag \end{align} and in general \[ \Gamma\qty(n+\frac12) = \frac{2n-1}{2} \times \frac{2n-3}{2} \times \cdots \frac{1}{2} \Gamma\qty(\frac12) = \frac{(2n-1)!!}{2^n} \sqrt{\pi} \] where the double factorial is defined by \[ n!! = n \times (n-2) \times \cdots \times \begin{cases} 1 & n\text{ odd} \\ 2 & n\text{ even} \end{cases} \]
The gamma function \(\Gamma(x)\) for small arguments \(x\). For large values of \(x\), \(\Gamma(x) = (x-1)!\).
The figure was generated with the following Python code:
import matplotlib.pyplot as plt
from scipy.special import gamma
f,a = plt.subplots()
x = np.linspace(0.5, 4, 100)
a.plot(x, gamma(x))
a.set_xlabel("$x$")
a.set_ylabel(r"$\Gamma(x)$");
Factorials grow extremely rapidly with their argument; it could be handy to have an approximate expression. It might be easier to estimate \(\ln n!\), though, since the logarithm converts all the multiplications to addition: \[ \ln n! = \ln 1 + \ln 2 + \cdots + \ln n \]
Figure 1 — A way to estimate \(\ln x!\) would be to compute the area under the red curve, although that sure looks like an overestimate, since there is area outside the blue bars. By contrast, the area under the green curve would be an underestimate, since it is entirely contained in the bars.
The area under the red curve would be \[ \ln n!_{\text{red}} = \int_{1}^{n+1} \ln x \dd{x} = x\ln x - x \bigg|_1^{n+1} = (n+1) \ln (n+1) - n - 1 \] The upper limit on the integral is \(n+1\) because the bar for \(\ln n\) starts at \(n\) and goes to \(n+1\). The area under the green curve would be \[ \ln n!_{\text{green}} = \int_{2}^{n+1} \ln (x-1) \dd{x} = \int_1^{n} \ln x’ \dd{x’} = x’\ln x’ - x’ \bigg|_1^{n} = n \ln n - n \] The true value lies above green but below red. So, crudely, \(\ln n! \approx n \ln n - n\) and \(n! \approx n^n e^{-n}\), although that’s going to be a modest underestimate. As a spot check, we can evaluate \(6! = 720\): \[ 6^6 e^{-6} = 115.65 \] Wow, that’s not very good, is it? Well, the red estimate would be \(7^7 e^{-7} = 751\), so at least the upper bound seems pretty reasonable.
We know from above that \(\Gamma(n+1) = n!\) for nonnegative integers \(n\). That is, \[ n! = \int_0^{\infty} x^n e^{-x} \dd{x} \] The integrand is the product of a rapidly increasing function of \(x\) (at large \(n\)) and a rapidly decreasing function of \(x\). There’s a peak somewhere. Maybe we could develop an approximation by expanding about that peak.
Figure 2 — The integrand of \(\Gamma(x)\), shown here for \(n = 6\). It is peaked in the vicinity of \(x = n\).
Actually, it will be smoother if we expand the logarithm of the integrand around the peak, and then exponentiate our series expansion. So, the steps are:
\begin{equation} y = n \ln x - x \qqtext{so} y’ = \frac{n}{x} - 1 \qqtext{and so} y’ = 0 \qqtext{when} x = n \end{equation} Evidently, the peak at \(y' = 0\) occurs when \(x = n\), so we will look to expand about that point.
\begin{align} i &= 2 & \qquad y’’ &= -\frac{n}{x^2} &\qquad y’’ (n) &= -\frac{1}{n} \notag \\ i &= 3 & \qquad y’’’ &= \frac{2! \; n}{x^3} &\qquad y’’’ (n) &= \frac{2}{n^2} \notag \\ i &= 4 & \qquad y^{\text{iv}} &= -\frac{3! \; n}{x^4} &\qquad y^{\text{iv}} (n) &= \frac{3!}{n^3} \notag \\ i &= 5 & \qquad y^{\text{v}} &= \frac{4! \; n}{x^5} &\qquad y^{\text{v}} (n) &= \frac{4!}{n^4} \notag \end{align}
I think I see the pattern here. So \begin{equation}\label{eq:Taylor} y = n \ln n - n - \frac{(x-n)^2}{2!\;n} + \frac{2!(x-n)^3}{3!\;n^2} - \frac{3!(x-n)^4}{4!\; n^3}+ \cdots \end{equation} We can simplify the notation by letting \(\xi = x-n\): \begin{equation}\label{eq:ybetter} y = n \ln n - n - \frac{\xi^2}{2n} + \frac{\xi^3}{3n^2} - \frac{\xi^4}{4n^3} + \cdots \end{equation}
To get our feet wet, let’s just include the first term in this series with \(\xi\) dependence. Then we have \begin{align} n! &\approx \int_{-n}^\infty \exp\left[ n\ln n - n - \frac{\xi^2}{2n} \right] \dd{\xi} \notag &= n^n e^{-n} \int_{-n}^\infty e^{-\alpha \xi^2} \dd{\xi} \notag \end{align} where I have defined \(\alpha \equiv 1/2n\). If we could let the lower limit go off to \(-\infty\), we would know that the integral was \(\sqrt{\pi/\alpha} = \sqrt{2 \pi n}\). Looking at Fig. 2, it might seem like a reasonable approximation; the integrand seems really small there. Our new approximation to \(n!\) is \begin{equation}\label{eq:nfact2} n! \approx \sqrt{2 \pi n} \; n^n e^{-n} \end{equation} When \(n = 6\), this approximation gives 710.1, which is a heck of a lot closer to the true value of 720 than our first approximation.
Figure 3 — The integrand of \(\Gamma(7) = 6!\) (black curve) and the approximate integrand from Eq. (\ref{eq:ybetter}) when we include only the term quadratic in \(\xi = x-n\).
Looking at the approximate integrand for \(n = 6\) in Fig. 3, it seems clear that we are overcounting area to the left of the peak and undercounting to its right. It seems likely that for larger values of \(n\) we are probably going to be doing better. Let’s make a quick check:
| $$n$$ | $$n!$$ | $$n^n e^{-n} / n!$$ | $$\sqrt{2\pi n} \, n^n \, e^{-n} /n!$$ |
|---|---|---|---|
| 1 | 1 | 0.367879 | 0.922137 |
| 2 | 2 | 0.270671 | 0.959502 |
| 3 | 6 | 0.224042 | 0.972702 |
| 4 | 24 | 0.195367 | 0.979424 |
| 5 | 120 | 0.175467 | 0.983493 |
| 6 | 720 | 0.160623 | 0.986220 |
| 7 | 5040 | 0.149003 | 0.988174 |
| 8 | 40320 | 0.139587 | 0.989643 |
| 9 | 362880 | 0.131756 | 0.990787 |
| 10 | 3628800 | 0.125110 | 0.991704 |
So, that hunch seems to be right; by \(n = 10\), our approximation that includes the Gaussian integral is closer than 1% to the true value.
We now need to return to Eq. (\ref{eq:ybetter}). Let \begin{equation}\label{eq:q} q = \frac{\xi^3}{3n^2} - \frac{\xi^4}{4n^3} + \frac{\xi^5}{5n^4} - \frac{\xi^6}{6n^5} + \cdots \end{equation} Then, \begin{equation}\label{eq:nf2} n! \approx n^n e^{-n} \int_{-\infty}^{\infty} e^{-\alpha n^2} e^q \dd{\xi} \end{equation} If we now expand \(e^q\) in a Taylor series, \[ e^q = 1 + q + \frac{q^2}{2!} + \frac{q^3}{3!} + \cdots \] then we will have terms of the form \begin{equation}\label{eq:blah} I_{\mu} = \int_{-\infty}^{\infty} e^{-\alpha \xi^2} \xi^{\mu} \dd{\xi} \end{equation}
Before sweating the details, note that if \(\mu\) is odd, then the integrand is odd and so the integral evaluates to zero. So, we only need to worry about even values of \(\mu = 2m\). We can evaluate these by differentiating under the integral sign:
\begin{align}
I_{2m} &= \int_{-\infty}^{\infty} \xi^{2m} e^{-\alpha \xi^2} \dd{\xi}
= \int_{-\infty}^{\infty} \left( -\pdv{}{\alpha}\right)^m \;e^{-\alpha \xi^2} \dd{\xi} \notag \\
&= \left( -\pdv{}{\alpha}\right)^m \int_{-\infty}^{\infty} e^{\alpha \xi^2} \dd{\xi} = \left( -\pdv{}{\alpha}\right)^m \sqrt{\frac{\pi}{\alpha}}
= \sqrt{\frac{\pi}{\alpha}} \frac{1}{2} \times \frac{3}{2} \times \cdots \times \frac{2m-1}{2} \alpha^{-m} \notag \\
&= \sqrt{2 \pi n} \frac{(2m-1)!!}{2^m} (2n)^m = \sqrt{2 \pi n} \; n^m \; (2m-1)!!
\end{align}
So, how do we keep track of which terms are important as we expand the exponential of \(q\)? Keeping in mind that the odd powers of \(\xi\) integrate to zero, the terms we need from \(q\) are all negative: \begin{align} \label{eq:qterms} q: &\quad -\sqrt{2\pi n} \left[ \frac{1}{4n^{3}} n^{2} (3!!) + \frac{1}{6n^{5}} n^{3} (5!!) + \frac{7!! \; n^{4}}{8n^{7}} + \cdots \right]\\ &\quad -\sqrt{2 \pi n} \left[ \frac{3}{4n} + \frac{5}{2n^{2}} + \frac{105}{8n^{3}} + \cdots \right] \end{align} What about the next term, \(q^{2}/2\). What are the terms with even powers of \(\xi\)? \begin{align} q^{2} & = \frac{\xi^{6}}{(3n^2)^{2}} + \frac{\xi^{8}}{(4n^{3})^{2}} + 2\frac{\xi^{3}}{3n^{2}} \frac{\xi^{5}}{5n^{4}} + \cdots \\ &= \frac{\xi^{6}}{9n^{4}} + \frac{\xi^{8}}{n^{6}} \left( \frac{1}{16} + \frac{2}{15} \right) \end{align} On integrating, the first term will yield \begin{equation}\label{eq:blah2} \frac{1}{2} \sqrt{2 \pi n} \frac{1}{9n^4} n^3 (5!!) = \sqrt{2\pi n} \frac{5}{6n} \end{equation} So, noting that \(\frac{5}{6} - \frac{3}{4} = \frac{1}{12}\), through the terms proportional to \(1/n\), we get \begin{equation}\label{eq:Stirling} n! \approx \sqrt{2 \pi n} \, n^n e^{-n} \left(1 + \frac{1}{12n} + \cdots \right) \end{equation} which is called Stirling’s approximation. As you can see from the table below, the correction makes the approximation converge extremely well, being off by a mere 32 parts per million at the rather modest value of \(n = 10\).
| $$n$$ | $$n!$$ | $$\frac{n^n e^{-n}}{n!}$$ | $$\frac{\sqrt{2\pi n} \, n^n \, e^{-n} }{n!}$$ | $$\frac{\sqrt{2\pi n} \, n^n \, e^{-n}}{n!} \left(1 + \frac{1}{12n} \right) $$ |
|---|---|---|---|---|
| 1 | 1 | 0.367879 | 0.922137 | 0.998982 |
| 2 | 2 | 0.270671 | 0.959502 | 0.999481 |
| 3 | 6 | 0.224042 | 0.972702 | 0.999721 |
| 4 | 24 | 0.195367 | 0.979424 | 0.999829 |
| 5 | 120 | 0.175467 | 0.983493 | 0.999885 |
| 6 | 720 | 0.160623 | 0.986220 | 0.999917 |
| 7 | 5040 | 0.149003 | 0.988174 | 0.999938 |
| 8 | 40320 | 0.139587 | 0.989643 | 0.999952 |
| 9 | 362880 | 0.131756 | 0.990787 | 0.999961 |
| 10 | 3628800 | 0.125110 | 0.991704 | 0.999968 |
For those up to a challenge, show that the next term in the series is \[ n! \approx \sqrt{2 \pi n} n^n e^{-n} \left(1 + \frac{1}{12 n} + \frac{1}{288 n^2} + \cdots \right) \]
Perhaps a more insightful way to show the error arising from the various versions of the Stirling series is to take the absolute difference between the true value and the approximate value, and divide that by the true value.
Figure 4 — Relative error in the various approximations to the factorial function. The dots represent the true values; the dashed curves are power laws proportional to \(1/n\), \(1/n^2\), and \(1/n^3\) chosen to pass through the points at \(n = 30\).