Stirling's Formula and Asymptotic Series

斯特灵(stirling)公式与渐近级数 苏剑林 · 2016-04-15 · Mathematics series integral analysis

Stirling's approximation, or Stirling's formula, was originally proposed as an approximation of the factorial:

$$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$$

Pushing Stirling's formula to higher precision gives what is called the 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.

Related references:

https://zh.wikipedia.org/zh-cn/斯特灵公式

https://en.wikipedia.org/wiki/Stirling%27s_approximation

This post will discuss Stirling's formula together with an improved derivation of its asymptotic series, and will explain why the asymptotic series is asymptotic. more

Stirling's Formula

There are several ways to derive the Stirling series. Wikipedia gives a derivation based on the Euler–Maclaurin summation formula. However, it has some shortcomings. Because it works by summing a series $n!=\sum_{k=1}^n \ln k$, it cannot directly give the value of the constant $\sqrt{2\pi}$, and this line of reasoning cannot be directly generalized to the corresponding series for the gamma function ($\Gamma(x+1)$, as the generalization of the factorial, satisfies the same kind of approximation formula, but a derivation based on summing series cannot establish this).

The English Wikipedia also gives a derivation using Laplace's method. Building on that, I carried out a further derivation that directly yields the asymptotic series of the gamma function, which I'd like to share here. The idea of this method is to consider

$$\Gamma(x+1)=\int_0^{\infty}e^{-t}t^x dt$$

For any positive integer $n$, we have $\Gamma(n+1)=n!$, so we only need to focus on estimating the integral above. Setting $t=xs$ and substituting in, we get

$$\Gamma(x+1)=x^{x+1}\int_0^{\infty}e^{-x(s-\ln s)} ds$$

At this point, Wikipedia simply invokes Laplace's method. But ordinary readers are not familiar with Laplace's method, and moreover Laplace's method only gives a limited approximation—it cannot produce a full series solution. My improvement comes in precisely at this step: setting $s=e^u$ and substituting in, we obtain

$$\Gamma(x+1)=x^{x+1}\int_{-\infty}^{\infty}e^{-x(e^u-u)+u} du$$

which can be further rewritten 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 contribution to the integral comes from $\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}$$

The Stirling Series

Expanding the factor $e^{-x(e^u-1-u-u^2/2)+u}$ in the integrand as a series in $u$,

$$e^{-x(e^u-1-u-u^2/2)+u}=\sum_{k=0}^{\infty}a_k u^k$$

and then multiplying by $e^{-xu^2/2}$ and integrating gives the Stirling series. In order to keep track of the order of each term, and to make the derivation more amenable to computer algebra, some further processing is needed. The method follows An Implicit-Function Solution to a Nonlinear Difference Equation; the main trick is still to manually introduce an auxiliary parameter.

First, it is easy to see that

$$\int_{-\infty}^{\infty}e^{-xu^2/2}u^kdu=c_k u^{-(k+1)/2}$$

where $c_k$ is independent of $x$. We can regard the series as being organized by order in $x^{-1/2}$, with each factor of $u$ contributing one power of $x^{-1/2}$. So we introduce a parameter $u\to qu$, $x\to x/q^2$ (more precisely $x^{-1/2}\to qx^{-1/2}$, since as already noted, the ordering is by powers of $x^{-1/2}$), i.e., we consider

$$\int_{-\infty}^{\infty}e^{-xu^2/2-x/q^2\cdot(e^{qu}-1-qu-q^2 u^2/2)+q u} du$$

Setting $q=1$ recovers exactly the case we need. Expanding the integrand as a series in $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 gives

$$\sqrt{\frac{2 \pi }{x}}\left(1+\frac{1}{12x}+\frac{1}{288x^2}\dots\right)$$

which 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\cdot u^2/2 - x/q^2\cdot (Exp[q\cdot u] - 1 - q\cdot u - q^2\cdot u^2/2) +
q\cdot u], {q, 0, 5}]] /. q -> 1, {u, -Infinity, Infinity},
Assumptions -> x > 0]/Sqrt[2\cdot Pi]/Sqrt[x] // Expand

To speed up the computation, one can perform an Expand before the integration:

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

If even faster computation is needed, one could pre-compute the integrals analytically instead of relying on the built-in integration function. But for our purposes here—exploration and learning—that isn't necessary.

Why Is an Asymptotic Series Asymptotic?

If we fix a particular $x$ and substitute it into the series, summing infinitely many terms, the result is bound to diverge to infinity! This is because the series is asymptotic, with radius of convergence zero: taking the first few terms may work quite well, but if we actually sum infinitely many terms, the result diverges.

Why do asymptotic series arise in the first place? Let's start with a simple example we've discussed before, the integral

$$I(\varepsilon)=\int_{-\infty}^{\infty}e^{-x^2-\varepsilon x^4}dx$$

If we expand it term by term:

$$\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 again an asymptotic series. Why? It's a bit like a "weakest-link" effect: the radius of convergence of the Taylor series of a composite function is limited by whichever inner function has the smallest radius of convergence. We expand $e^{-\varepsilon x^4}$ as a series in $x$; although in theory this series has an infinite radius of convergence, that statement only holds in a limiting sense. For a truncated (finite-term) series, the expansion is not valid over the whole real line—the larger $x$ is, the smaller the range of validity for $\varepsilon$. But our integral runs from minus infinity to infinity, so it makes use of infinitely large values of $x$, for which $\varepsilon$ can only be treated as zero. That's why the series is asymptotic.

By contrast, if we consider the series expansion of

$$\int_{-M}^{M}e^{-x^2-\varepsilon x^4}dx$$

then no matter how large $M$ is, we always get a series with infinite radius of convergence. This is the essential difference between the finite and the infinite cases.

So why does an asymptotic series still have decent accuracy? That's because the main contribution to the integral $\int_{-\infty}^{\infty}e^{-x^2-\varepsilon x^4}dx$ comes from values near $x=0$, and the contribution from far away is greatly suppressed by the factor $e^{-x^2}$. So, for the first several terms, the series still gives good accuracy. However, although suppressed, these contributions don't vanish entirely, and when infinitely many weak terms are added together, the sum blows up to infinity.

The same reasoning explains the asymptotic nature of the Stirling series—in fact, it only converges at $x\to\infty$.