Stirling's Formula and Asymptotic Series

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 is derived by summing the series n! = \sum_{k=1}^n \ln k, it cannot directly provide the constant value \sqrt{2\pi}, and the results of this approach cannot be directly generalized 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).

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, we have \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 it typically only yields a finite approximation rather than a 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 This 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 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

By expanding the factor e^{-x(e^u-1-u-u^2/2)+u} of the integrand into a series in u: e^{-x(e^u-1-u-u^2/2)+u}=\sum_{k=0}^{\infty}a_k u^k Then, multiplying by e^{-xu^2/2} and integrating, we can obtain the Stirling series. To identify the order of each term and facilitate computer derivation, some processing is required. The processing method refers to "An Implicit Function Solution to a Nonlinear Difference Equation", mainly involving the manual introduction of parameters.

First, it is easy to obtain: \int_{-\infty}^{\infty}e^{-xu^2/2}u^kdu=c_k x^{-(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 each u contributes an x^{-1/2}. Therefore, we introduce parameters u \to qu and x \to x/q^2 (effectively x^{-1/2} \to qx^{-1/2} because, as mentioned, the order is 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 The case q=1 is exactly what we need. Expanding the integrand as a series in q, we get: 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 can pre-write the integration results yourself instead of using the built-in integration function. However, for our exploratory study, this is unnecessary.

Why are Asymptotic Series Asymptotic?

If you fix a value for x and substitute it into the series, calculating infinitely many terms will inevitably result in infinity! This is because it is an asymptotic series with a radius of convergence of 0. Taking the first few terms might yield good results, but if you truly calculate the sum of infinitely many terms, the result will diverge.

Why do asymptotic series occur? Let’s start with a simple example we have 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" (Liebig’s Law of the Minimum)—the radius of convergence of the Taylor series of a composite function depends on the function with the smallest radius of convergence. We expand e^{-\varepsilon x^4} as a series in x; although theoretically the radius of convergence of this series is infinite, this only holds in the limit. For a series with a finite number of terms, this expansion is not valid over the entire range of real numbers. As x becomes larger, the valid range for \varepsilon becomes smaller. Since our integration is performed between positive and negative infinity, involving infinite x, 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 result is always a series with an infinite radius of convergence. This is the difference between the finite and the infinite.

So why do asymptotic series still possess a certain degree of accuracy? This is because the main contribution to the integral \int_{-\infty}^{\infty}e^{-x^2-\varepsilon x^4}dx comes from values near x=0. The contribution from 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 superimposed, they become infinite.

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

When reposting, please include the original address of this article: https://kexue.fm/archives/3731

For more detailed reposting matters, please refer to: Scientific Space FAQ