jiyang1971的个人博客分享 http://blog.sciencenet.cn/u/jiyang1971

博文

计算方法之祖冲之的精度 精选

已有 6896 次阅读 2018-10-31 09:27 |个人分类:大众物理学|系统分类:科普集锦

 

  

山巅一寺一壶酒,尔乐苦煞吾。

 

$\pi \approx 3.1415926535 \approx \frac{22}{7} \approx \frac{355}{113}$

 


祖冲之(429年-500年)是我国古代伟大的数学家。他改进了刘徽(约225年—约295年)的割圆术,把圆周率计算到小数点后第7位,即$3.1415926 < \pi <3.1415927$,还给出了约率($22/7 \approx 3.143$) 和密率($355/113 \approx 3.1415929$)。《算经十书》收录过祖冲之的《缀术》五卷,然而“学官莫能究其深奥,故废而不理”(《隋书》)。祖冲之到底是怎么计算$\pi$的,我们现在并不清楚。


往事越千年。随着科学知识的普及,每个大学一年级的学生,都可以轻松地把$\pi$计算到祖冲之的精度,但是,他们中的一些人并没有意识到自己的能力。

中学生的方法

对于半径为1的圆,其周长$2\pi$和面积$\pi$都可以用来计算圆周率。这个圆的内切正$n$边形和外切正$n$边形的周长和面积,在$n$趋于无穷大的时候,都趋近于圆的周长和面积。

内切正$n$边形的周长是$2n\sin\frac{\pi}{n}$,外切正$n$边形的周长是$2n\tan\frac{\pi}{n}$,对于所有的$n$,都有

$n\sin\frac{\pi}{n} < \pi <n\tan\frac{\pi}{n}$

当然,我们不能用这个公式来计算$\pi$值,否则就循环论证了。但是,我们可以选择适当的正$n$边形,并逐渐增大$n$。最常用的方法是从正6边形开始,然后逐步加倍,12、24、48、96,等等。这是因为内切正6边形的边长是1,外切正6边形的边长是$\sqrt{3}/3$,都是很好算的数值。

内切正$6\cdot 2^k$边形的边长是$a_k=2 \sin \theta _k$,其中,$\theta _k =\frac{\pi}{6\cdot2^k}$。如果知道了$a_k$,就很容易知道

$a_{k+1}=2\sin \theta _{k+1}= 2\sin\theta _k /2 $

$=2\sqrt{\frac{1 - cos \theta_k}{2}}=\sqrt{2}\frac{\sin \theta _k}{\sqrt{1 + \cos \theta_k}}$

$=\sqrt{2}\frac{\sin \theta _k}{\sqrt{1 + \sqrt{1-\sin ^2\theta_k}}}=\sqrt{2}\frac{a_k/2}{\sqrt{1+\sqrt{1-a^2_k/4}}}$

再根据$a_0=2\sin \frac{\pi}{6}=1$,就可以递推地得到所有的$a_{k}$,而$\pi =6\cdot 2^k a_k$。

同样的方法,可以得到外切正$6\cdot 2^k$边形的边长$b_k$的递推公式是

$b_{k+1}=2\tan \theta _{k+1}= 2\tan \theta _k /2$

$=\frac{1-\cos \theta _k }{1+\cos _k }=\frac{b_k}{\sqrt{1+b^2_k/4}+1}$

再根据$b_0=2\tan \frac{\pi}{6}=2\sqrt{3}/3$,就可以递推地得到所有的$b_{k}$,而$\pi =6\cdot 2^k b_k$。

利用这些公式,就可以得到$\pi$值如下 

$k$

正多边形

内切$\pi$值

外切$\pi$值

0

6边形

3.00

3.4641016

1

12边形

3.1058285

3.2153903

2

24边形

3.1326286

3.1596599

3

48边形

3.1393502

3.1460862

4

96边形

3.1410320

3.1427146

5

192边形

3.1414525

3.1418730

6

384边形

3.1415576

3.1416627

7

768边形

3.1415839

3.1416102

8

1536边形

3.1415905

3.1415970

7

3072边形

3.1415921

3.1415937

8

6144边形

3.1415925

3.1415929

9

12288边形

3.1415926

3.1415927

10

24576边形

3.1415926

3.1415927

 

表面看来,刚才这些计算涉及了三角函数,其实不用三角函数也可以得到相应的递推公式,因为只涉及到勾股定理。比如说,对于内切正多边形,AB的长度为$a_k$,AC的长度为$a_{k+1}$,则

$a_{k+1}=AC=\sqrt{AD^2+CD^2}=\sqrt{(a_k/2)^2+(1-\sqrt{1-(a_k/2)^2})^2}=\sqrt{2}\frac{a_k/2}{\sqrt{1+\sqrt{1-a^2_k/4}}}$

这里只涉及开平方和加减乘除,原则上谁都可以算的,但是计算太繁琐了,如果没有计算器,是很难算的。中国古代没有我们现在的阿拉伯数字系统,计算起来很不方便,虽然说有人能够“照位运筹如飞,人眼不能逐”( 沈括《梦溪笔谈·技艺》),恐怕也不能采用蛮力来计算$\pi$值。


 

 

大学生的方法

在单位圆里,弧边梯形AOCB的面积减去三角形BOC的面积,就是扇形角AOB(张角为$\phi$)的面积($\phi /2$),由此可以得到圆周率$\pi$。

AOCB的面积是$\int ^{\sin \phi}_0 \sqrt{1-x^2}dx$,而BOC的面积是$\sin \phi \cos \phi /2$,扇形角AOB的面积就是

$\phi =\int ^{\sin \phi}_0 \sqrt{1-x^2}dx - \sin \phi \cos \phi /2$

把积分里的$\sqrt{1-x^2}$用二项式定理展开,就可以得到

$\phi\approx \int ^{\sin \phi}_0 dx(x-\frac{1}{2}x^2 -\frac{1}{8}x^4-\frac{1}{16}x^6-\frac{5}{128}x^8)-\frac{1}{4}\sin 2\phi$

$=x(1-\frac{1}{6}x^2-\frac{1}{32}x^4-\frac{1}{112}x^6-\frac{5}{1152}x^8)^{\sin \phi}_0 -\frac{1}{4}\sin 2\phi$

选择一些简单的$\phi$值,就可以得到$\pi$值

$\phi$角度

弧度

$\pi$值

90

$\pi/2$

3.1552579

60

$\pi/3$

3.1295483

30

$\pi/6$

3.1392663

15

$\pi/12$

3.1414185

7.5

$\pi/24$

3.1415813

3.75

$\pi/48$

3.1415919

1.875

$\pi/96$

3.1415926

0.9375

$\pi/192$

3.1415927

0.46875

$\pi/384$

3.1415927

 


需要注意的是,上面这两种方法都是野蛮处理,计算只涉及加减乘除和开方(用于计算$\sin$ 值),原则上来说,可以用手算出来的(开方可以用迭代法得到),尽管我这里的结果来自于计算机(我手算过15度,确认是可以算的,可是后来发现结果有点错儿)。这些方法甚至都没有用到刘徽不等式,也没有给出$\pi$值的上下限,从证明技巧上来说,还远不如古人,更不用说符合现代数学的严格要求了。然而,科学知识的普及,使得我们这个时代的普通人也能轻松达到古代杰出人物的计算水平,却是无可置疑的。

 

 

 



http://blog.sciencenet.cn/blog-1319915-1143713.html

上一篇:计算方法之锦鲤吉祥
下一篇:计算方法之圆周率计算的补充说明

17 信忠保 张磊 杨正瓴 黄永义 史晓雷 赵克勤 夏力钢 武夷山 韩枫 郑永军 杨荣佳 张江敏 徐传胜 刘浔江 薛亮 吕秀齐 李颖业

该博文允许注册用户评论 请点击登录 评论 (29 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备14006957 )

GMT+8, 2018-11-20 07:51

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部