斯特灵(stirling)公式与渐近级数
By 苏剑林 | 2016-04-15 | 58794位读者 |斯特灵近似,或者称斯特灵公式,最开始是作为阶乘的近似提出
$$n!\sim \sqrt{2\pi n}\left(\frac{n}{e}\right)^n$$
符号$\sim$意味着
$$\lim_{n\to\infty}\frac{\sqrt{2\pi n}\left(\frac{n}{e}\right)^n}{n!}=1$$
将斯特灵公式进一步提高精度,就得到所谓的斯特灵级数
$$n!=\sqrt{2\pi n}\left(\frac{n}{e}\right)^n\left(1+\frac{1}{12n}+\frac{1}{288n^2}\dots\right)$$
很遗憾,这个是渐近级数。
相关资料有:
https://zh.wikipedia.org/zh-cn/斯特灵公式
https://en.wikipedia.org/wiki/Stirling%27s_approximation
本文将会谈到斯特灵公式及其渐近级数的一个改进的推导,并解释渐近级数为什么渐近。
斯特灵公式 #
斯特灵级数有几种推导方法,维基上给出了基于欧拉-麦克劳林求和公式的推导。然而,它有一些不足的地方。因为它是通过求级数和$n!=\sum_{k=1}^n \ln k$的方法来求的,它不能直接给出常数值$\sqrt{2\pi}$,而且这种思路的结果不能直接推广到伽马函数对应的级数($\Gamma(x+1)$作为阶乘的推广,具有同样的近似公式,但是基于级数求和的思路不能得出这一点。)。
英文维基还给出了利用拉普拉斯方法的一种推导。我在此基础上做了进一步的推算,能够直接得到伽马函数的渐近级数,在这里跟大家分享一下。这种方法的思路是考虑
$$\Gamma(x+1)=\int_0^{\infty}e^{-t}t^x dt$$
对于任意正整数$n$,都有$\Gamma(n+1)=n!$,因此我们只需要致力于上述积分的估计就行了。设$t=xs$,代进入得到
$$\Gamma(x+1)=x^{x+1}\int_0^{\infty}e^{-x(s-\ln s)} ds$$
到这里,维基上就直接说利用拉普拉斯方法了。但是,一般的读者并不熟悉拉普拉斯方法,而且通过拉普拉斯方法只能得到一个有限的近似,没有办法得到级数解。笔者的改进正是在这里,设$s=e^u$,代入得到
$$\Gamma(x+1)=x^{x+1}\int_{-\infty}^{\infty}e^{-x(e^u-u)+u} du$$
进一步改写成
$$\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$$
事实上,当$x$充分大的时候,积分的主要部分由$\int_{-\infty}^{\infty}e^{-xu^2/2}du=\sqrt{2\pi/x}$贡献,从而直接得到斯特灵公式
$$\Gamma(x+1)\approx \sqrt{2\pi x}\left(\frac{x}{e}\right)^{x}$$
斯特灵级数 #
将被积函数的因子$e^{-x(e^u-1-u-u^2/2)+u}$展开为$u$的级数
$$e^{-x(e^u-1-u-u^2/2)+u}=\sum_{k=0}^{\infty}a_k u^k$$
然后乘上$e^{-xu^2/2}$再进行积分,就可以得到斯特灵级数。为了辩明各项的阶,并且方便计算机推导,需要做一些处理。处理方法参考《一个非线性差分方程的隐函数解》,主要还是人工引入参数。
首先,可以很容易得到
$$\int_{-\infty}^{\infty}e^{-xu^2/2}u^kdu=c_k u^{-(k+1)/2}$$
其中$c_k$与$x$无关。可以看成,级数应该以$x^{-1/2}$为阶的,而一个$u$则贡献了一个$x^{-1/2}$,故引入参数$u\to qu$,$x\to x/q^2$(事实上是$x^{-1/2}\to qx^{-1/2}$,因为已经说了,以$x^{-1/2}$为阶),即考虑
$$\int_{-\infty}^{\infty}e^{-xu^2/2-x/q^2\cdot(e^{qu}-1-qu-q^2 u^2/2)+q u} du$$
当$q=1$时正是我们需要的情况。将被积函数展开为$q$的级数,得到
$$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$$
代入$q=1$,然后逐项积分,得到
$$\sqrt{\frac{2 \pi }{x}}\left(1+\frac{1}{12x}+\frac{1}{288x^2}\dots\right)$$
这就得到了斯特灵级数
$$\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)$$
Mathematica代码为
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
为了加快计算速度,可以在积分之前做一次Expand:
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
如果需要进一步加速计算,可以自己预先写入积分的计算结果,而不是用内置的积分函数。但是就我们的学习探究而言,此举就没有必要了。
渐近级数为何渐近? #
如果取定一个$x$,代入到级数中,算无穷多项,结果必然是无穷大!因为这个是渐近级数,收敛半径是0,取前几项可能效果还不错,但是如果真的计算无穷多项的和,结果就会出现发散情况。
为什么会出现渐近级数呢?我们先从一个我们以前讨论过的简单的例子,考虑积分
$$I(\varepsilon)=\int_{-\infty}^{\infty}e^{-x^2-\varepsilon x^4}dx$$
如果将它展开计算:
$$\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}$$
这也是渐近级数。为什么呢?这有点像短板效应——复合函数的泰勒级数的收敛半径取决于收敛半径最小的那个函数。我们展开$e^{-\varepsilon x^4}$为$x$的级数,虽然在理论上,这个级数的收敛半径为无穷,但这只是在极限的形式下成立。而对于取有限项的级数,这个展开式并非在整个实数范围内有效,当$x$取得越大的时候,$\varepsilon$的取值范围就越小。而我们的积分是在正负无穷之间进行的,用到了无穷的$x$,这时候$\varepsilon$的取值只能为0。因此级数的渐近的。
对比之下,如果考虑
$$\int_{-M}^{M}e^{-x^2-\varepsilon x^4}dx$$
的级数展开式,不管$M$多大,得到的总是一个收敛半径为无穷的级数。这就是有限与无限的差别。
那为什么渐近级数又具有一定的精确度?那是因为,对积分$\int_{-\infty}^{\infty}e^{-x^2-\varepsilon x^4}dx$的主要贡献是在$x=0$附近的值,无穷远处的贡献被因子$e^{-x^2}$大大削弱了,因此,在前几项的时候,级数还是有不错的精确度。然而,虽然被削弱,但毕竟还是存在的,当无穷多个弱项叠加起来时,却成为了无穷大~
用同样的理由,可以解释斯特灵级数的渐近性——事实上它只在$x\to\infty$处收敛。
转载到请包括本文地址:https://kexue.fm/archives/3731
更详细的转载事宜请参考:《科学空间FAQ》
如果您还有什么疑惑或建议,欢迎在下方评论区继续讨论。
如果您觉得本文还不错,欢迎分享/打赏本文。打赏并非要从中获得收益,而是希望知道科学空间获得了多少读者的真心关注。当然,如果你无视它,也不会影响你的阅读。再次表示欢迎和感谢!
如果您需要引用本文,请参考:
苏剑林. (Apr. 15, 2016). 《斯特灵(stirling)公式与渐近级数 》[Blog post]. Retrieved from https://kexue.fm/archives/3731
@online{kexuefm-3731,
title={斯特灵(stirling)公式与渐近级数},
author={苏剑林},
year={2016},
month={Apr},
url={\url{https://kexue.fm/archives/3731}},
}
May 4th, 2016
苏剑林,你个傻逼。
请文明用语,谢谢!
我惊呆了
March 12th, 2022
不同的方法,风景依旧迷人。
鞍点近似
对于积分 $I(\lambda)=\int_C e^{-\lambda g(z)}dz$ , 要让被积函数下降最快,就要找到一条最速下降路径 $C^\ast$
将$g(z)$在其驻点$z_0$展开
$$g(z)\sim g(z_0)+\frac12g''(z_0)(z-z_0)^2$$
设$z=u+iv$,由于$g(z)$的海森矩阵为
$$H[g(z)]=\begin{pmatrix}
\frac{\partial^2 g}{\partial u^2}&\frac{\partial^2 g}{\partial u\partial v}\\\\
\frac{\partial^2 g}{\partial v\partial u}&\frac{\partial^2 g}{\partial v^2}
\end{pmatrix}=
\begin{pmatrix}
2&0\\\\
0&-2
\end{pmatrix}$$
其为不定矩阵,故$z_0$是$g(z)$的鞍点,于是
$$I(\lambda)=\int_C e^{-\lambda g(z)}dz\sim e^{-\lambda g(z_0)}\int_{C^\ast}e^{-\frac{\lambda}{2}g''(z_0)(z-z_0)^2}dz$$
接下来找最陡路径
设$g(z)=\phi(z)+i\psi(z)$, $g(z)$的解析性意味着$\nabla\phi\cdot\nabla\psi=0$,即两者正交。
则在沿着$\psi=\eta\in\mathrm{const}$的路径$C^\ast$上,$\frac{d\phi}{ds}=|\nabla\phi|$,即此时$\phi$变化最快。
设$z-z_0=x+iy$,则$(z-z_0)^2=x^2-y^2+i2xy$,所以$\begin{cases}\phi=x^2-y^2\\\\\psi=2xy\end{cases}$
通过鞍点的最陡路径只有一条,$z=z_0$意味着$z-z_0=x+iy=0$,即$\begin{cases}x=0\\\\y=0\end{cases}$
在实轴路径上$\phi=x^2$上升最快,在虚轴路径上$\phi=-y^2$下降最快,我们取实轴路径,这样积分收敛最快
此时$z-z_0=x$,因此
$$I(\lambda)=\int_C e^{-\lambda g(z)}dz\sim e^{-\lambda g(z_0)}\int_{-\infty}^{+\infty}e^{-\frac{\lambda}{2}g''(z_0)x^2}dx$$
最终可得
$$\int_C e^{-\lambda g(z)}dz\sim e^{-\lambda g(z_0)}\sqrt\frac{2\pi}{\lambda g''(z_0)}$$
Γ函数的渐近行为
$$\Gamma(n+1)=\int_0^\infty x^ne^{-x}dx=\int_0^\infty e^{-(x-n\ln x)}dx$$
这里$g(x)=x-n\ln x$,其中鞍点$x_0$满足$g'(x_0)=1-\frac nx=0$,即$x_0=n$
所以$g(x_0)=n-n\ln n,\ g''(x_0)=\frac1n$,于是有
$$\Gamma(n+1)\sim\sqrt{2n\pi}e^{-n+n\ln n}=\sqrt{2n\pi}\left(\frac{n}{e}\right)^n$$
强!复变玩得好也让人赏心悦目。
March 13th, 2022
很早就关注博主了,你对数学的思考与方法,也让我对数学分析有了不同的理解,加油