[Math] Higher Order Terms in Stirling’s Approximation

Question

Some websites and books give stirling approximation as

$$n! = \sqrt{2 \pi n} \left( \frac{n}{e} \right)^n \left( 1 + O \left(\frac{1}{n} \right)\right)$$

However when I check their derivations most do it analytically up to

$$n! \approx \sqrt{2 \pi n} \left(\frac{n}{e}\right)^n$$

So I tried proving the $$O\left(\frac {1}{n}\right) $$ term but to no avail.

Here are my steps, hopefully someone can point me in the right direction

First I began by defining
$$n!=\Gamma(n+1) = \int{e^{-x}x^{n} }dx $$
$$=\int{dxe^{nlnx-x}}$$

The maximum of the integrand is at $$x =n$$

Next, I will expand the powers by a taylor approximation about x=n

$$f(x) = -x+nlnx$$
$$f'(x) = -1 +nx^{-1}$$
$$f''(x) = -nx^{-2}$$
$$f^{3}(x)=2nx^{-3}$$

By letting $$ x = n+\epsilon$$

I get

$$f(n+\epsilon) = -n + nlnn + 0\epsilon + -\frac{1}{n^2}\epsilon^2 + \frac{2}{n^{2}\epsilon^3} $$

Plugging back into the integration

$$\Gamma(n+1) = n! = n^{n}e^{-n}\int{e^{-\frac{\epsilon^2}{2n}+\frac{\epsilon^3}{3n^2}+…}}d\epsilon +…$$

Normally if it were up to the second derivative term of the taylor's expansion, it would be a gaussian integral. But if there are more terms, how would I do it?
Or is there another approach?

Answer 1

The Euler-MacLaurin summation formula can be used to approximate $n!$ conveniently.

The following is taken from Concrete Mathematics by D.E. Knuth ad R.L. Graham. We find in section 9.5 among other representations this version of the Euler-MacLaurin summation formula for integer values $a\leq b$

\begin{align*} \sum_{a\leq k < b}f(k)&=\int_a^bf(x)dx-\frac{1}{2}\left.f(x)\right|_a^{b}+\sum_{k=1}^{m}\frac{B_{2k}}{(2k)!}\left.f^{(2k-1)}(x)\right|_a^b\tag{1}\\ &\qquad+ \mathcal{O}\left((2\pi)^{-2m}\right)\int_a^b\left|f^{(2m)}(x)\right|dx \end{align*}

with $B_k$ the Bernoulli numbers starting with \begin{align*} B_0=1,B_1=-\frac{1}{2},B_2=\frac{1}{6},B_4=-\frac{1}{30},B_6=\frac{1}{42},\ldots \end{align*} and $B_k=0$ for odd $k$ greater $1$.

It is also shown, that if the function $f$ has the nice property $$f^{2m}(x)\geq 0$$ for $a\leq x \leq b$ then it turns out that the remainder in (1) can be replaced with \begin{align*} \Theta_m\frac{B_{2m+2}}{(2m+2)!}\left.f^{(2m+1)}(x)\right|_a^b\qquad\qquad \text{for some }0<\Theta<1 \end{align*}

and we obtain \begin{align*} \sum_{a\leq k < b}f(k)&=\int_a^bf(x)dx-\frac{1}{2}\left.f(x)\right|_a^{b}+\sum_{k=1}^{m}\frac{B_{2k}}{(2k)!}\left. f^{(2k-1)}(x)\right|_a^b\\ &\qquad+\Theta_m\frac{B_{2m+2}}{(2m+2)!}\left.f^{(2m+1)}(x)\right|_a^b\tag{2} \end{align*}

---

Example: $f(x)=\ln x$

In order to find the Stirling formula for $n!$ we consider $f(x)=\ln x$ and derive following approximation for $(n-1)!$

\begin{align*} \sum_{1\leq k < n}\ln k&=n\ln n - n +\sigma-\frac{1}{2}\ln n+\sum_{k=1}^m\frac{B_{2k}}{2k(2k-1)n^{2k-1}}\\ &\qquad+\Theta\frac{B_{2m+1}}{(2m+2)(2m+1)n^{2m+1}}\qquad\qquad 0<\Theta<1 \end{align*}

Another instructive example in section 9.6, namely the asymptotic estimation for $\sum_k\binom{2n}{k}$ provides a formula for $\sigma$ \begin{align*} e^\sigma=\sqrt{2\pi} \end{align*}

Now we have everything we need to estimate $n!$ as precise as we like:

Let's assume we want to approximate $n!$ up to terms of order $n^{-3}$. We put $m=2$ and observe

\begin{align*} \sum_{1\leq k < n}\ln k=n\ln n-n+\sqrt{2\pi}-\frac{1}{2}\ln n+\frac{1}{12n}-\frac{1}{360n^3} +\mathcal{O}\left(\frac{1}{n^5}\right)\tag{3} \end{align*}

Using the generating function of $\exp(x)=\sum_{n\geq 0}\frac{x^n}{n!}$ we obtain from (3) \begin{align*} n!&=\exp\left(n\ln n-n+\sqrt{2\pi}+\frac{1}{2}\ln n+\frac{1}{12n}-\frac{1}{360n^3} +\mathcal{O}\left(\frac{1}{n^5}\right)\right)\tag{4}\\ &=\sqrt{2\pi n}\left(\frac{n}{e}\right)^n \exp\left(\frac{1}{12n}-\frac{1}{360n^3}+\mathcal{O}\left(\frac{1}{n^5}\right)\right)\\ &=\sqrt{2\pi n}\left(\frac{n}{e}\right)^n \left(1+\left(\frac{1}{12n}-\frac{1}{360n^3}\right)+\frac{1}{2}\left(\frac{1}{12n}-\frac{1}{360n^3}\right)^2\right.\\ &\qquad\qquad\qquad\qquad\left.+\frac{1}{6}\left(\frac{1}{12n}-\frac{1}{360n^3}\right)^3+\mathcal{O}\left(\frac{1}{n^4}\right)\right)\\ &=\sqrt{2\pi n}\left(\frac{n}{e}\right)^n \left(1+\left(\frac{1}{12n}-\frac{1}{360n^3}\right)+\frac{1}{2}\left(\frac{1}{144n^2}\right)\right.\\ &\qquad\qquad\qquad\qquad\left.+\frac{1}{6}\left(\frac{1}{1728n^3}\right)+\mathcal{O}\left(\frac{1}{n^4}\right)\right) \end{align*}

Observe the change from $-\frac{1}{2}\ln n$ to $+\frac{1}{2}\ln n$ from (3) to (4) due to the change from $(n-1)!$ to $n!$.

We finally conclude \begin{align*} n!=\sqrt{2\pi n}\left(\frac{n}{e}\right)^n \left(1+\frac{1}{12n}+\frac{1}{288n^2}-\frac{139}{51840n^3}+\mathcal{O}\left(\frac{1}{n^4}\right)\right) \end{align*}