Better Stirling Approximation Error

Question

Stirling's Approximation is given by

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

From a substantial improvement of the Stirling formula, we have an elegant approximation given by

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

We can write a simple function that has Stirling's, the improved version versus actual as

  nf[n_] := {Sqrt[2. Pi n] (n/E)^n, Sqrt[2. Pi n] (n/E + 1/(12 E n))^n, 1. n!}

My question, how can we find the percent error between Actual versus the two approximations and make it part of the function specification?

Maybe even plot the three variants.

Answer 1

n1[n_] = Sqrt[2 π n] (n/E)^n;
n2[n_] = Sqrt[2 π n] (n/E + 1/(12 E n))^n;

Series[n1[n] / n!, {n, ∞, 1}]
(*    1 - 1/(12 n) + O[1/n]^2    *)

Series[n2[n] / n!, {n, ∞, 4}]

(*    1 - 1/(1440 n^3) + O[1/n]^5    *)

We see that $n_2$ approximates $n!$ with a relative error that decreases as $n^{-3}$, whereas the error of $n_1$ decreases as $n^{-1}$ (the inverse of the well-known Stirling series $1+\frac{1}{12n}+\frac{1}{288n^2}+\ldots$).

Ramanujan's formula is better:

n3[n_] = Sqrt[π] (n/E)^n (8 n^3 + 4 n^2 + n + 1/30)^(1/6);
Series[n3[n] / n!, {n, ∞, 4}]
(*    1 + 11/(11520 n^4) + O[1/n]^5    *)

@GregHurst's suggestion is better still:

n4[n_] = Exp[-x] Sqrt[2π/x] (x Sqrt[1/(810 x^6) + x Sinh[1/x]])^x /. x -> n + 1;
Series[n4[n] / n!, {n, ∞, 7}]
(*    1 + 163/(340200 n^7) + O[1/n]^8    *)

Answer 2

I offer another way to compare two approximations of n!. We compare the Mortici formula (from the paper) to the two-term Stirling series (which the paper fails to do):

mortici = Sqrt[2 \[Pi] n] (n/E + 1/(12 E n))^n;
stirling2 = Sqrt[2 \[Pi] n] (n/E)^n (1 + 1/(12 n));
factorial = n!;
Series[( (* takes a few seconds *)
     mortici - factorial)/(stirling2 - factorial) // Simplify,
    {n, Infinity, 10}] // Normal // FullSimplify[#, n > 0] &;
Series[%, {n, Infinity, 2}]
(*  1/(5 n) + 77/(450 n^2) + O[1/n]^3  *)
{1/(5 n), (mortici - factorial)/(stirling2 - factorial)} /.
  {{n -> 100.}, {n -> 2500.`30}} // N
(*  1/(5n)     ratio
  {{0.002,   0.00201741},
   {0.00008, 0.0000800274}}
*)

Thus the Mortici formula has an absolute error approximately 1/(5n) times the error of the two-term Stirling formula.

---

Numerics note: It may be better to use logarithms and LogGamma, since the numbers can get very large. However, the arbitrary-precision numbers seem to work well enough up to n = 2500. In any case, here's the set-up for an logarithmic approach:

logmortici = mortici // Log // FullSimplify[#, n > 0] &;
logstirling2 = stirling2 // Log // FullSimplify[#, n > 0] &;
logfactorial = LogGamma[n + 1];
(* rel/abs errors *)
relmortici = Internal`Expm1[logmortici - logfactorial];
absmortici = relmortici factorial;
relstirling2 = Internal`Expm1[logstirling2 - logfactorial];
absstirling2 = relstirling2 factorial;

The function Internal`Expm1[x] computes Exp[x]-1 accurately when Exp[x] is near 1 (x near zero).

Note also that Exp[2500.] overflows machine-precision and is automatically promoted to arbitrary precision with $MachinePrecision digits (approx. 15.95). However, for the absolute error mortici - factorial for n = 2500., we need at least 19 digits to get a one-digit result.

Answer 3

The relative errors may be appended by:

nf[n_] := {t1 = Sqrt[2. Pi n] (n/E)^n, 
  t2 = Sqrt[2. Pi n] (n/E + 1/(12 E n))^n, 
  t3 = 1. n!, (t1 - t3)/t3, (t2 - t3)/t3
}

and plotted like:

DiscretePlot[nf[i][[-2 ;; -2]], {i, 1, 20}, PlotRange -> All]
DiscretePlot[nf[i][[-1 ;; -1]], {i, 1, 20}]