Stirling's Formula and Asymptotic Series

By 苏剑林 | April 15, 2016

Stirling's approximation, or Stirling's formula, was initially proposed as an approximation for factorials: $$n!\sim \sqrt{2\pi n}\left(\frac{n}{e}\right)^n$$ The symbol $\sim$ means: $$\lim_{n\to\infty}\frac{\sqrt{2\pi n}\left(\frac{n}{e}\right)^n}{n!}=1$$ Further increasing the precision of Stirling's formula leads to the so-called Stirling series: $$n!=\sqrt{2\pi n}\left(\frac{n}{e}\right)^n\left(1+\frac{1}{12n}+\frac{1}{288n^2}\dots\right)$$ Unfortunately, this is an asymptotic series.

This article will discuss an improved derivation of Stirling's formula and its asymptotic series, and explain why asymptotic series are "asymptotic."

Stirling's Formula

There are several ways to derive the Stirling series. Wikipedia provides a derivation based on the Euler-Maclaurin summation formula. However, it has some shortcomings. Because it calculates the series sum $n!=\sum_{k=1}^n \ln k$, it cannot directly yield the constant value $\sqrt{2\pi}$, and this line of reasoning cannot be directly extended to the series corresponding to the Gamma function ($\Gamma(x+1)$, as a generalization of the factorial, shares the same approximation formula, but the summation-based approach cannot derive this point).

The English Wikipedia also provides a derivation using Laplace's method. I have performed further calculations based on this, which can directly yield the asymptotic series for the Gamma function, and I would like to share it with everyone. The idea of this method is to consider: $$\Gamma(x+1)=\int_0^{\infty}e^{-t}t^x dt$$ For any positive integer $n$, $\Gamma(n+1)=n!$, so we only need to focus on estimating the above integral. Let $t=xs$, substituting this in gives: $$\Gamma(x+1)=x^{x+1}\int_0^{\infty}e^{-x(s-\ln s)} ds$$ At this point, Wikipedia simply states to use Laplace's method. However, general readers may not be familiar with Laplace's method, and the method itself typically only yields a finite approximation rather than a full series solution. My improvement lies here. Let $s=e^u$, substituting this in gives: $$\Gamma(x+1)=x^{x+1}\int_{-\infty}^{\infty}e^{-x(e^u-u)+u} du$$ Further rewriting this as: $$\Gamma(x+1)=x\left(\frac{x}{e}\right)^{x}\int_{-\infty}^{\infty}e^{-xu^2/2-x(e^u-1-u-u^2/2)+u} du$$ In fact, when $x$ is sufficiently large, the main part of the integral is contributed by $\int_{-\infty}^{\infty}e^{-xu^2/2}du=\sqrt{2\pi/x}$, which directly yields Stirling's formula: $$\Gamma(x+1)\approx \sqrt{2\pi x}\left(\frac{x}{e}\right)^{x}$$

Stirling Series

Expand the factor $e^{-x(e^u-1-u-u^2/2)+u}$ of the integrand into a power series of $u$: $$e^{-x(e^u-1-u-u^2/2)+u}=\sum_{k=0}^{\infty}a_k u^k$$ Then multiply by $e^{-xu^2/2}$ and integrate to obtain the Stirling series. To discern the order of each term and facilitate computer derivation, some handling is required. The processing method refers to "Implicit Function Solution for a Nonlinear Difference Equation," which primarily involves manually introducing parameters.

First, it is easy to obtain: $$\int_{-\infty}^{\infty}e^{-xu^2/2}u^kdu=c_k u^{-(k+1)/2}$$ where $c_k$ is independent of $x$. It can be seen that the series should be in orders of $x^{-1/2}$, and one $u$ contributes an $x^{-1/2}$, so we introduce parameters $u\to qu$ and $x\to x/q^2$ (effectively $x^{-1/2}\to qx^{-1/2}$, because as stated, the order is based on $x^{-1/2}$), i.e., consider: $$\int_{-\infty}^{\infty}e^{-xu^2/2-x/q^2\cdot(e^{qu}-1-qu-q^2 u^2/2)+q u} du$$ When $q=1$, this is exactly the case we need. Expanding the integrand as a power series of $q$ gives: $$e^{-u^2 x}+q e^{-u^2 x} \left(u-\frac{u^3 x}{6}\right)+\frac{1}{72} q^2 u^2 e^{-\frac{1}{2} u^2 x} \left(u^4 x^2-15 u^2 x+36\right)\dots$$ Substituting $q=1$ and integrating term by term, we get: $$\sqrt{\frac{2 \pi }{x}}\left(1+\frac{1}{12x}+\frac{1}{288x^2}\dots\right)$$ This yields the Stirling series: $$\Gamma(x+1)\approx \sqrt{2\pi x}\left(\frac{x}{e}\right)^{x}\left(1+\frac{1}{12x}+\frac{1}{288x^2}\dots\right)$$

The Mathematica code is:

Integrate[ Normal[Series[ Exp[-x*u^2/2 - x/q^2*(Exp[q*u] - 1 - q*u - q^2*u^2/2) + q*u], {q, 0, 5}]] /. q -> 1, {u, -Infinity, Infinity}, Assumptions -> x > 0]/Sqrt[2*Pi]/Sqrt[x] // Expand

To speed up the calculation, you can perform an Expand before integrating:

Integrate[ Normal[Series[ Exp[-x*u^2/2 - x/q^2*(Exp[q*u] - 1 - q*u - q^2*u^2/2) + q*u], {q, 0, 16}]] /. q -> 1 // Expand, {u, -Infinity, Infinity}, Assumptions -> x > 0]/Sqrt[2*Pi]/Sqrt[x] // Expand

If further acceleration is needed, you could write in the integration results manually instead of using the built-in integration function. But for our learning exploration, this step is unnecessary.

Why are Asymptotic Series Asymptotic?

If you pick a fixed $x$ and substitute it into the series to calculate infinitely many terms, the result will inevitably be infinite! Because this is an asymptotic series, its radius of convergence is 0. Taking the first few terms might give good results, but if you actually sum an infinite number of terms, the result will diverge.

Why do asymptotic series occur? Let's start with a simple example we've discussed before, considering the integral: $$I(\varepsilon)=\int_{-\infty}^{\infty}e^{-x^2-\varepsilon x^4}dx$$ If we expand and calculate it: \begin{aligned} &\int_{-\infty}^{\infty}e^{-x^2}\left(1-\varepsilon x^4+\frac{1}{2}\varepsilon^2 x^8\dots\right)dx\\ =&\sqrt{\pi}\left(1-\frac{3}{4}\varepsilon+\frac{105}{32}\varepsilon^2\dots\right) \end{aligned} This is also an asymptotic series. Why? This is somewhat like the "short-board effect"—the radius of convergence of a composite function's Taylor series depends on the function with the smallest radius of convergence. We expand $e^{-\varepsilon x^4}$ as a power series of $x$; although theoretically, the radius of convergence of this series is infinite, this only holds in the limit form. For a series where we take only finite terms, this expansion is not valid over the entire range of real numbers; as $x$ becomes larger, the allowable range for $\varepsilon$ becomes smaller. Our integration is performed between negative and positive infinity, using infinite $x$, so the value of $\varepsilon$ can only be 0. Therefore, the series is asymptotic.

In contrast, if we consider the series expansion of: $$\int_{-M}^{M}e^{-x^2-\varepsilon x^4}dx$$ No matter how large $M$ is, the resulting series will always have an infinite radius of convergence. This is the difference between finite and infinite.

So why do asymptotic series still possess some accuracy? That is because the primary contribution to the integral $\int_{-\infty}^{\infty}e^{-x^2-\varepsilon x^4}dx$ comes from values near $x=0$. The contribution at infinity is greatly weakened by the factor $e^{-x^2}$. Therefore, in the first few terms, the series still has decent accuracy. However, although weakened, it still exists; when infinitely many weak terms are summed together, they become infinite~

Using the same reasoning, we can explain the asymptotic nature of the Stirling series—in fact, it only converges at $x\to\infty$.