历史上的谜案——刘徽有没有使用外推法？

话说当年我国古代数学家刘徽创立“割圆术”计算圆周率的事迹，在今天已被不少学生知晓；虽不能说家喻户晓，但是也为各教科书以及老师津津乐道。和古希腊的“数学之神”阿基米德同出一辙，刘徽也是使用圆的内接、外切正多边形来逼近圆形的；不一样的是，刘徽使用的方法是计算半径为1的圆的内接、外切正多边形的面积，而阿基米德计算的则是直径为1的圆的内接、外切正多边形的周长。两者的计算效果有什么区别呢？其实阿基米德的方法应该更快一点，阿基米德算到正n边形所得到的值，相当于刘徽算到正2n边形了。

在此我们不再对两者的计算方法进行区分，因为两者的本质都是一样的。按照现代数学的写法，“割圆术”的理论依据是
$$lim_{n\to \infty} n \sin(\frac{\pi}{n})=\pi\tag{1}$$
当然，刘徽不可能有现代计算正弦函数值的公式（现在计算正弦函数值一般用泰勒级数展开，而泰勒级数展开需要用到$\pi$的值），甚至在他那个时代就连笔墨也没有，据我所知即使是后来的祖冲之推算圆周率时，唯一的计算工具也只是现在称为“算筹”的小棍。不过刘徽还是凭借着超强的毅力，利用递推的方法逐步求圆周率。

设$\pi(n)=n sin(\frac{\pi}{n})=2n sin(\frac{\pi}{2n})cos(\frac{\pi}{2n})$，而$\pi(2n)=2n sin(\frac{\pi}{2n})$，结合$sin^2 \theta+cos^2 \theta=1$，不难解得
$$\pi(2n)=n\sqrt{2-2\sqrt{1-(\frac{\pi(n)}{n})^2}}\tag{2}$$这就是割圆术需要的递推公式。

刘徽用(2)式一直算，从n=6开始，算到了192边形的面积，即相当于$\pi(96)$的值，即是3.14103...，于是按照正常逻辑，刘徽最大取值3.14或3.141，不能再精确了，再取多几位也没有什么意义。但是，刘徽却创造了一种“捷法”，进一步求得了圆周率的近似值。首先他求出了半径为10的内接正96边形的面积为$313\frac{584}{625}$，内接正192边形的面积为$314\frac{64}{625}$（分别相当于$\pi(48),\pi(92)$），求出两者差

$$314\frac{64}{625}-313\frac{584}{625}=\frac{105}{625}$$

然后“以十二觚之幂为率消息”，“取此分寸之三十六”，即$\frac{36}{625}$，加上到正192边形的面积，得到
$$314\frac{64}{625}+\frac{36}{625}=314.16$$
于是他求出了精度达5位小数的圆周率！引号中的话是刘徽原话，没人知道刘徽的捷法是什么。只知道这种方法可以根据数列的有限个值，推导出该数列的极限。这就有点像现在被称为“外推法”的思想了，于是有人认为：刘徽已经使用了外推法！而且后来祖冲之也是改进了外推法（而不是改进了割圆术）才推导得到圆周率为3.1415926，否则凭借那时候的计算能力，不大可能有这样的结果。真想是否这样？我们不得而知。

外推加速圆周率

外推法是怎样起作用的呢？我们知道：
$$\sin x=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+...$$
令$h=\frac{1}{n}$，就有
$$\begin{aligned}\pi(n)=\pi-\frac{h^2 \pi^3}{3!}+\frac{h^4 \pi^5}{5!}-\frac{h^6 \pi^7}{7!}+... \\ \pi(2n)=\pi-(\frac{1}{2})^2 \frac{h^2 \pi^3}{3!}+(\frac{1}{2})^4 \frac{h^4 \pi^5}{5!}-(\frac{1}{2})^6 \frac{h^6 \pi^7}{7!}+...\end{aligned}$$

于是有$4\pi(2n)-\pi(n)=3\pi-(\frac{3}{4})\frac{h^4 \pi^5}{7!}...$，读者可以看到，我们消去了$h^2$项，使精度达到了$O(h^3)$，于是可以使用
$$\pi^{\cdot \cdot }=\frac{1}{3}[4\pi(2n)-\pi(n)]$$
来对已经求得的近似值就行修正。对于刘徽的数据，我们立马可以求得$\pi^{**}=3.141584$，接近3.1416，怪不得有人会认为刘徽使用了外推法！

我们还可以继续写出
$$\pi(4n)=\pi-(\frac{1}{4})^2 \frac{h^2 \pi^3}{3!}+(\frac{1}{4})^4 \frac{h^4 \pi^5}{5!}-(\frac{1}{4})^6 \frac{h^6 \pi^7}{7!}+...$$
用$\pi(n),\pi(2n),\pi(4n)$三者的线性组合来消去$h^2,h^4$项，得到
$$\pi^{\cdot \cdot }=\frac{1}{45}[64\pi(4n)-20\pi(2n)+\pi(n)]$$

这个公式精度为$O(h^5)$。如果你还嫌不够，我还可以写出
$$\pi^{\cdot \cdot }=\frac{1}{2835}[4096\pi(8n)-1344\pi(4n)+84\pi(2n)-\pi(n)]$$

对于最后一个公式，精度为$O(h^7)$，我们只需要令n=3，就可以得到$\pi^{**}=3.141592651...$。换句话说，计算到正24边形就可以达到祖冲之的精度！

外推法综述

对于电子计算机时代，即使使用外推法，用割圆术求圆周率仍然还是显得“过时”了，本文旨在通过这个实例，来总结一下这个星期对外推法的一点粗浅研究。数学家Romberg被公认为是外推法的创始人，真正的外推法并不是像上述那样通过增加线性组合的项来实现提高精度的，而是通过对修正值的多次反复“再修正”来实现的。具体情形可以参考《外推法及其应用》。

至于外推法的应用，可以说非常非常广泛。因为外推法要提高精度，仅仅涉及到线性组合这一基本运算，因此对于数值计算有着重要的意义。《计算方法导论》中的很多公式，如Simpson积分公式、解常微分方程的RK方法等等，都是外推法的实例，具体大家可以参考相关书籍，如《外推法及其应用》等。

转载到请包括本文地址：https://kexue.fm/archives/1292

更详细的转载事宜请参考：《科学空间FAQ》