Rainbow Zyop/文

双纽线积分

$$ I_{0}:=\int_{0}^{1}\frac{1}{\sqrt{1-x^4}}\,\mathrm{d}x $$

以及

$$ I_2:=\int_{0}^{1}\frac{x^2}{\sqrt{1-x^4}}\,\mathrm{d}x $$

来自 Jacob Bernoulli 1694年的研究。这是两个无法用初等函数表示的积分。自这个积分被Bernoulli提出后，有很多数学家都在研究这两个积分。这里我们提其中的两个研究：其一是数值计算，其二是分析学的研究。

[注：题图是 James Stirling 1730年的著作，有Springer 2003年的英译本。这位是数值计算的大前辈。阶乘的Stirling公式正是以他的名字命名。]

---

Stirling，1692年生于苏格兰。1724年，Stirling前往伦敦，那时科学巨人Newton已经到了生命的最后几年。Stirling受Newton青睐，在1735年离开伦敦之前在学术方面是非常活跃的(特别是数学)。他最重要的著作Methodus Differentialis就是在伦敦期间写成。

Stirling在这本书中考虑了双纽线积分以及其他若干级数及积分的数值计算问题。Stirling的目标是对这些积分以及级数进行精度尽可能高的数值计算。我们从一个简单的例子出发，来描述Stirling的思路。

学过一点微积分的人一定十分清楚：

$$ \int_{0}^{1}\frac{1}{\sqrt{1-x^2}}\,\mathrm{d}x=\frac{\pi}{2} $$

假设我们不知道圆周率的值，那么，如何估计左边积分的值呢？

计算这类积分最直接的方法并不是利用各种梯形法则和Simpson法则，而是利用幂级数。Newton早就知道，对被积函数作幂级数展开，可以用来计算积分。二项级数

$$ \frac{1}{\sqrt{1-x^2}}=1+\frac{1}{2}x^2+\frac{1}{2}\cdot\frac{3}{4}x^4+\cdots $$

在区间$[0,1)$上是收敛的。对幂级数作逐项积分，我们可以得到积分的值为

$$ 1+\frac{1}{2}\cdot\frac{1}{3}+\frac{1}{2}\cdot\frac{3}{4}\cdot\frac{1}{5}+\cdots $$

能不能用这个级数来进行求和呢？理论上是可以的，但是通过简单的分析，可以知道，级数第k项满足渐进公式 $a_{k}\sim ck^{-3/2}$ 。也就是说，如果我们取前N项求和，那么和式与积分的误差在 $N^{-1/2}$ 的量级。那么，为了计算这个积分到小数点后7位，我们就需要对 $10^{14}$ 项求和。这显然不合适。为了计算此类级数，我们就需要一些数值加速的方法。

我们设上面级数各项为 $a_i(i=0,1,2,\cdots,n,\cdots)$ 。那么 $a_0=1,a_{i+1}=\frac{(i+1/2)^2}{(i+1)(i+3/2)}a_i$. Stirling指出，如果从第$N$项开始求和，那么级数和可以用和式首项乘以一个线性因子来近似，也就是

$$ \sum_{i=N}^{\infty}a_i=(p+qN)a_N+\sum_{i=N}^{\infty}b_i $$

$p,q,b_n$ 待定。

从上面的关系式我们可以得到

$$b_{N}=\left(1-p-qN+\frac{(N+1/2)^2}{(N+1)(N+3/2)}(p+q(N+1))\right)a_N.$$

为了使待定的 $b_{n}$ 尽可能地小，我们就需要调整 $p，q$，使得上式括号里的式子通分之后分子中含N的项的次数尽可能地小。取$p=4/3，q=2$就可以消去分子中所有关于$N$的项。

化简后我们可以得到 $b_{N}=\frac{1}{3(N+1)(N+3/2)}a_N$ 。新的级数各项趋于0的速度明显比原级数的速度快，从 $N^{-3/2}$ 加速到 $N^{-7/2}$ 。注意到 $b_{N+1}=\frac{(N+1/2)^2}{(N+2)(N+5/2)}b_N $,所以我们可以用同样的方法加速新级数的收敛。

回到双纽线积分。通过幂级数展开我们容易得到

$$ I_{0}=1+\frac{1}{2}\cdot\frac{1}{5}+\frac{1}{2}\cdot\frac{3}{4}\cdot\frac{1}{9}+\cdots $$ $$ I_{2}=\frac{1}{3}+\frac{1}{2}\cdot\frac{1}{7}+\frac{1}{2}\cdot\frac{3}{4}\cdot\frac{1}{11}+\cdots $$

这样的级数第N项仍然是以 $N^{-3/2}$ 的速度趋于0的。这样的收敛速度对于级数的数值计算极为不利，因此Stirling在计算这两个积分时就采用了上面提到的数值加速方法。

Stirling首先取 $I_0$ 对应级数的前9项的和。和为 $\frac{19148291449}{15750758400}\approx 1.21570599730613606$。然后Stirling运用上面的算法对从第十项开始的和式进行反复加速。

原级数各项满足递推关系 $a_0=1,a_{i+1}=\frac{(i+1/2)(i+1/4)}{(i+1)(i+5/4)}a_i$。

用Stirling的方法，第一次加速后的数列 $b_{i}$ 满足递推式 $b_{i+1}=\frac{(i+1/2)(i+1/4)}{(i+2)(i+5/4+1)}b_i$ ，与 $a_{i}$ 的关系为 $b_{N}=\frac{1}{2\times1-1/2}\times\frac{1-1/2}{1-1/2}\times\frac{1-1/4}{N+1}\times\frac{1/4-1/2+1}{N+1/4+1}a_N$

$a_i$ 求和近似值为 $(p+qN)a_N=\frac{(N+1)(2-1/2)-1(1-1/4)}{(2-1/2)(1-1/2)}a_N$

第二次加速后的数列满足递推式 $c_{i+1}=\frac{(i+1/2)(i+1/4)}{(i+3)(i+5/4+2)}c_i$，与 $b_{i}$ 的关系为 $c_{N}=\frac{2}{2\times2-1/2}\times\frac{2-1/2}{3-1/2}\times\frac{2-1/4}{N+2}\times\frac{1/4-1/2+2}{N+1/4+2}b_N$，$b_i$ 求和近似值为 $(p^\prime+q^\prime N)b_N=\frac{(N+3)(4-1/2)-2(2-1/4)}{(3-1/2)(4-1/2)}b_N$

以此类推。

Stirling对数列进行了9次加速，各次加速所得的和的近似值列表如下：

$$\begin{equation}\begin{split}0.09524164973078547\\0.00008069254020837\\0.00000043223735955\\0.00000000522484720\\0.00000000010371186\\0.00000000000290174\\0.00000000000010478\\0.00000000000000461\\0.00000000000000023\\ \end{split}\end{equation}$$

把这些近似值与原级数前9项的和加在一起，就得到Stirling的结果：

$$I_{0}\approx1.31102877714605987$$

尽管Stirling计算出的结果准确到了小数点后15位，但是他书中没有数值计算的误差分析。这是那个时代数值计算的一大缺陷。我们把它交给有兴趣的读者来完成。用相同的算法，Stirling也同时得到了 $I_2 \approx0.59907011736779611$ 。

---

Stirling在他的书中还计算了很多有趣的级数。他计算出了 $\sum_{n=1}^{\infty}\frac{1}{n^2}\approx1.64493406684822643.$ 但数值计算是代替不了理论分析的。Euler自己在1735年就这样写道[原文与英译文来自Andre Weil的书An Approach Through History from Hammurapi to Legendre，下面中译文就来自Weil的英译]：

Tantopere iam pertractatae et investigatae sunt series reciprocae potestatum numerorum naturalium, ut vix probabile videatur de iis novi quicquam inveniri posse... Ego etiam iam saepius...has series diligenter sum persecutus neque tamen quicquam aliud sum assecutus, nisi ut earum summam......proxime veram definiverim...

关于自然数的平方的倒数和已经有大量的研究，以致于人们不太可能找到什么新发现……我也一样，尽管尝试了很多次，除了得到级数的近似值外并没有新发现……但是峰回路转，

... Deductus sum autem nuper omnino inopinato ad elegantem summae huius seriei $1+\frac{1}{4}+\frac{1}{9}+\frac{1}{16}+\mathrm{etc.}$ expressionem, quae a circuli quadratura pendet...

但是相当出人意料的是，我找到了这个级数基于圆面积的一个漂亮的公式……

根据Weil的记载，Euler马上把自己的发现通知了住在Basel的Daniel Bernoulli以及当时住在Edinburgh的Stirling。Stirling在1738年10月26号给Maclaurin的信中这样写道：

I have lately had a letter from Mr Euler at Petersburgh,...his last one is full of great many ingenious things, but it is long and I am not quite master of all the particulars.

另一个Euler的发现也没有被Stirling捕捉到。在1738年12月20日写给Johann Bernoulli的信件中，Euler对自己的老师提到[译文仍然来自Weil的同一本书]：

Observavi nuper insignem elasticae rectangulae proprietatem, in qua si abscissa ponatur $x$, est applicata $=\int \frac{xxdx}{\sqrt{a^4-x^4}}$ et longitudo curvae $=\int \frac{aadx}{\sqrt{a^4-x^4}}$ quae expressiones ita sunt comparatae, ut inter se comparari nequeant. At si abscissa sumatur $= a$, inveni rectangulum sub applicata et arcu comprehensum aequale esse areae circuli cuius diameter sit abscissa = a; quae observatio mihi quidem notatu maxime digna videtur.

我新近发现了关于弹性曲线的一个引人注目的性质，该曲线由[$f(x)=\int_{0}^{x}\frac{u^2\,\mathrm{d}u}{\sqrt{a^4-u^4}}$ 定义]，其弧长为[$ s(x)=\int_{0}^{x}\frac{a^2\,\mathrm{d}u}{\sqrt{a^4-u^4}} $]。两个积分互相之间[一般]没有关联，但是如果取x=a，那么有[$ f(a)s(a)=\frac{1}{4}\pi a^2 $]，这对我而言是最值得注意的发现。

如果我们用Stirling的计算检验一下(令$a=1$)，就可以发现， $I_{0}\times I_{2}$ 与 $\pi/4$ 的差绝对值在 $10^{-16}$ 量级。这是巧合吗？Euler在1775年发表了一个简洁的证明：[Euler应当早有证明，只是很晚才发表出来]

证明：记 $I_n=\int_{0}^{1}\frac{x^n}{\sqrt{1-x^4}}\,\mathrm{d}x$ 。根据分部积分，我们有

$$ \begin{equation}\begin{split}I_{n+3}&=\int_{0}^{1}\frac{x^{n+3}}{\sqrt{1-x^4}}\,\mathrm{d}x\\&=\int_{0}^{1}-\frac{x^{n}}{2}\,\mathrm{d}\sqrt{1-x^4}\\&=\frac{n}{2}\int_{0}^{1}x^{n-1}\sqrt{1-x^4}\,\mathrm{d}x\\&=\frac{n}{2}\left(I_{n-1}-I_{n+3}\right)\end{split}\end{equation} $$

所以有递推式 $I_{n+4}=\frac{n+1}{n+3}I_{n},n\in\mathbb{N}$。

我们注意到， $I_{1}=\frac{\pi}{4},I_{3}=\frac{1}{2}$. 根据递推式，我们有

$$ I_{4n+4}I_{4n+6}=\frac{I_0I_2}{4n+5} $$ $$ I_{4n+5}I_{4n+7}=\frac{1}{4n+6}\frac{\pi}{4} $$

根据积分 $I_{n}$ 的单调性以及递推公式，我们立刻可以得到上面两个等式左边在n趋于无穷时比值为1。因此我们已经证明了 $I_{0}I_{2}=\pi/4$。

[2018. 1. 10 补注：Euler公布自己的这个公式比一般数学史中提到的1786年要早得多。他在1744年的变分法专著中有一个关于弹性曲线的附录(Additamentum 1)。此附录第27节中就已经给出了这个等式(它与Gauss还产生了奇妙的关联)。根据Weil的数论史，更早一些，也就是1739年1月，Euler就已经在圣彼得堡宣读过这个结果。]

---

Euler的发现仅仅是冰山一角。这也是本文起名为“预篇”的原因。这冰山一角在18世纪中叶仅仅崩落了一点点。更多的发现，要等到19世纪的英雄们出场。