实验2π的计算

. 22||||2

12188 .

224 4242))21

(11()21(|||||| .22 . . 11121

111221212

2

21122N N ∈?=?=????===?=∈--=--=--+=+====++-+n a S a OE AB S S n a a a a a EF AF AE a S a n S a n n n OAE n n n n n n ，

，故又注意到，相应的边数，故，

注意到，边长和面积，则步单位圆内接多边形的第分别是

，设值也就越接近圆周率近单位圆的面积，其数边形的面积就越来越接那么内接正多

正多边形的边数开始，逐步成倍地增加从单位圆的内接正方形Δ实验2：π的计算

学院：机械与动力工程学院 姓名：唐子彦 学号：515020910137

一、实验目的

通过求π的近似值，了解历史上计算π值的一些方法，包括刘徽割圆术、级数展开、数值积分和Monte Carlo 法等. 复习微积分中相关知识，比较它们的差异，了解计算方法对提高计算效率的意义.

二、实验内容

（一）利用鲁道夫割圆术计算π值

【问题】德国人鲁道夫一生计算圆周率，他同样是用圆的内接多边形逼近圆周. 不过，他是从正方形开始成倍增加边数. 试推导出他计算所采用的递推公式，然后求π的近似值到10位和20位有效数字. 【解】

图1.1

.

1 ,21 ,21 421 ,

2 1112???>∈?==?????>∈--==---n n a n S n n a n a n n n n n 且，且，

的递推公式：

由此可得鲁道夫割圆术N N 为通过“鲁道夫法”计算π的近似值，现编写M 文件如下：

为求出指定有效数字位数的π值，还需编写M 文件如下：

运行结果如下：

图1.2

【答】由图1.2不难发现，“鲁道夫法”计算π值效率较高，n=16时即可得到π的10位有效数字，n=32时即可得到π的20位有效数字. 然而，由于中途过程涉及开方等运算，因此计算较为复杂繁琐

.

（二）利用幂级数展开式计算π值

【问题】简单公式

31arctan 21arctan

4+=π

，Machin 公式239

1arctan 51arctan 44-=π，

以及公式8

1

arctan 51arctan 21arctan 4++=π. 试验证上述三个公式（分别记为公式

1、2、3），并利用反正切函数的幂级数展开式求π值，比较上述三公式的计算效率. 此外，再找出一种利用幂级数展开式求π的方法并验证之. 【解】

.

14

2

4

4

04

tan 13121131

21tan tan 1tan tan )tan(.

3

1tan ,21tan 31arctan ,21arctan )1(得证，公式，故且，因为，则记π

βαπ

π

π

βαπβαβαβαβαβα=

+=

+

<

+<==?-+

=

-+=+====.

24

42

03

404

tan 1239111912012391

119120tan 4tan 1tan 4tan )4tan(119

1204tan 125)5

(1512tan 1tan 22tan .

239

1

tan ,51tan 2391arctan ,51arctan )2(22得证，公式，故且，因为，则再记ππππ=

-<

-<

-<==?+-

=

+-=-=?=-?

=

-=?====y x y x y x y x y x x x x x y x y x .

34

.

431tan 1tan 1)4(tan 3

1

8

15118151tan tan 1tan tan )tan(.

8

1

tan ,51tan ,21tan 81arctan )3(得证，公式即

，故且，

因为，则再记z x z x z x z x z x z x z ++=+=-=+-=-=?-+

=-+=+====απ

απαααπα

为利用上述三个公式求出π值，现编写以下三个M文件：

第一个M文件，对应于公式1：

第二个M文件，对应于公式2：

第三个M文件，对应于公式3：

为比较它们计算指定有效数字位数的π值的效率，还需编写M文件如下（详见第五页）：

运行结果如下（每五个为一组）：

图2.1 图2.2

图2.3 图2.4 从上述四幅图中可以看出，公式1、3的计算效率基本相同，而公式2的计算效率高于其他两个.

然而，从有效数字位数m=15开始，所得结果中公式1、2、3所需项数不再增加，这可能是MATLAB本身的原因. 个人猜测：当MATLAB计算到一个与π极为

相近的数时，可能将其自动补全为π，而没有继续计算. 为试图解决该问题，可

改用C++进行编程，所需的CPP文件如下：

运行结果（每五个为一组）见第9页.

从运行结果来看，有效数字位数m=1~16均可得到正确的项数n1、n2、n3. 然而当m=17时，或许是由于C++语言中long double类型的计算精度有限，无法进行高精度的浮点运算，导致循环条件恒为真，即程序进入死循环，无法得出正确结果. 对于m>17的情形更是如此（参见图2.8）.

图2.5 图2.6

左图：图2.7 下图：图2.8

.

4)2

1(!)!2)(12(!)!12(2121arcsin 621.

!)!2)(12(!

)!12(d !)!2(!)!12(11d arcsin :)1,1(.

11

!)!2(!)!12(1],[)(:)1,1(],[!)!2(!)!12()( .)1,1(!)!

2(!)!12(1).

1,1(11 1).1,1(11

!)!2(!)!12(lim .

11

21lim 21lim 121!)!2(!)!12(21:.

!)!2(!)!12(lim lim 2,!)!2(!)!12(!)21

(12,

0!

)2

1(43!21211)]

(1[11

.

)2

1(!)!2)(12(!)!12(2121arcsin 6 4.)4(11

20

1112202

2],[122122

22221021

0422

1

22

112得证，公式，即得到取可求积定理”知：再由“函数项级数逐项且，则记内内闭一致收敛在第二定理，级数据无意义，故收敛域为时，又当，即收敛区间为，则收敛半径故，且注意到，则记证明如下：：公式利用以下公式计算π值除此之外，我们还可以∑?∑∑?∑∑

∏∏∑∞=+∞

=∞

=+∞

=∞

=∞

→∞→∞→∞→∞→-=-=-∞=++-+===+-+=-+=-=-∈?-?-+∈-??-=--+--±=-==-==+=+<-≤

∈?-==???????=-=+-==?

+++?+?++=-+=-+-+==n n x

n n n n x

b a n n n n n n n

k

k n n n n k k n n ____

n k i n n n i n n n n n x x n n n x t t n n t t

x x x

x n n b a C x u b a x n n x u x n n Abel x

x r k k ρn n n n n n

n k k a ρk

n k k k i k n a x n i x x x x n n n πρ

π

N

为比较公式4与前三个公式的计算效率，现编写M 文件如下：

输入图2.9所示命令行，运行结果如下：

图2.10

【答】对比以上四公式的计算结果不难发现，公式1、3计算效率大致相同且较低，公式2的计算效率最高，公式4的计算效率介于它们之间比公式1、3略高。

(三)数值积分计算π值

【问题】利用数值积分计算π，分别用“梯形法”和“Simpson 法”精确到10位有效数字，再用“Simpson 法”精确到15位有效数字. 【解】在本题中，“梯形法”和“Simpson 法”的计算原理如下：

为利用数值积分计算π值，现编写如下两个M 文件：

].

42)[(32

].

42)[(61.11

2,,2,1,0 ,22]1,0[)2(23

(2.

23(21)(21]1)(21[.11,,2,1,0 ,[0,1])1(].

1,0[,11

.41d 4111

12220111122202

01

0000122

1

0def 2∑∑∑∑∑∑∑∑?-==--==-=-===++++≈+++≈+=?==

-≈-=--=?+≈?+=?==

∈+===+=m k m

k k k m m k m

k k k m k k k n k k n k n k k n k n k k k k

k k y y y y m y y y y m A x y m k m k x m Sim pson y n y n y y y n n y y A x y n k n k x n x x y A A x x ππππ

故计算，可以得到：曲线用抛物线代替进行将相邻两个区间的函数，则等分，令将区间法

故算，可以得到：曲线用直线段代替来计将每个小区间上的函数，则等分，令将区间梯形法

记，则由于

为计算指定有效数字位数的π值，还需编写以下两个M文件：(1)梯形法：

(2)Simpson法：

运行结果见第14页图3.1.

图3.1

【答】由图3.1知，利用“梯形法”和“Simpson 法”计算10位有效数字的π值所需的最少划分份数分别为12910和20；利用“Simpson 法”计算15位有效数字的π值所需的最少划分份数为92.

此外不难发现，当指定相同有效数字位数时，“Simpson 法”所需的划分份数远少于“梯形法”. 换言之，当划分份数一定时，“Simpson 法”的计算精度远高于“梯形法”.

综上所述，“Simpson 法”的计算效率远高于“梯形法”.

（四）利用Monte Carlo 法计算π值

【问题】(1)试用Monte Carlo 法计算π的5位有效数字； (2)设计方案试用计算机程序模拟Buffon 实验.

【解】首先将Monte Carlo 法和Buffon 实验的原理进行一些简要说明.

文件如下：

为此，需编写的点数个，并计算其中满足共上的随机数组生区间现用计算机程序模拟产从而的面积之比，即的面积与之比应趋于与投点总数内的点数则落在区域内随机地投大量的点，在区域，单位圆内，记正方形在坐标平面法：区域总面积

指定区域面积

率几何概型”，即待求概上述两种方法均基于“M .1),(]1,0[.

4,4 }.1|),{(}10,10|),{( )1(.

2222m y x n y x n

m

n m A B n m B A y x y x B y x y x A xOy Carlo Monte ≤+==≤+=≤≤≤≤==ππ

由于取n=10000运行该程序1次和取n=1000运行该程序10次再求平均值的本质没有区别，故每次进行投点实验时，仅运行该程序一次. 运行结果如下：

图4.1

从图4.1中可以看出，随着投点次数不断增加，所得的π值与真实值的误差总体上不断减小. 本次实验中，当投点次数为10亿，可得π的5位有效数字.

然而，投点10万次，即可求出π的4位有效数字；而投点1000万次，仅可求出π的3位有效数字. 这说明增加投点次数，不能说计算π的精度必定提高

.

上述公式的推导如下：

近似值：够大时，可以求出π的足，则当数为，针和平行线的相交次设投针次数为多次扔到该平面上的针随机地将长度为的平行线多条间距为在一个平面内画出足够实验：

.2 .)( .)2(n

m

d l n n m d l l d Buffon ?≈≤

π

图4.2

图4.3

.

2

2sin 2.

sin 2||).

,()2

2( )0( .2.42

d l d

l l l

h h d

h d h d

≤≤≤≤≤≤-<≤，故要求在该带形区域内，即又由于针的中点必须落相交的约束条件为：不难求出，针与该直线，设针的坐标为为中点到该直线的距离记，针的

针与该直线的夹角记为）见图所形成的带形区域（详两侧

可以只考虑一条直线及，故各直线地位相同，由于平行线数量足够多θθθπθθ

.22d )]sin 2(sin 2[3.40n

m d l d l d l l S S P m n ?≈?=--==≈?πππθθθπ

总内可知：由图

为利用计算机程序进行随机模拟，现编写M 文件如下：

运行结果如下：

图4.4

【答】尽管Monte Carlo 法和Buffon

实验的原理不复杂，操作也简单，但计算

.

]!

)!12(!)!2( !)!12(!)!2([2.

lim !)!12(!

)!2( !)!12(!)!2(lim !)!12(!)!22( !)!12(!)!2(lim 12212lim

.

lim }{.2

0]1)2()2([!)!12(!)!2( !)!12(!)!2(.2

212!)!12(!)!2( !)!12(!)!2(!)!12(!)!22( !)!12(!)!2(2!)!12(!)!2( !)!12(!)!2( .

!

)!12(!

)!22(2!)!2(!)!12(!)!12(!)!2( .

d sin d sin

sin

.

sin sin sin ,1sin 0:]2

,0[11def 12

220

122022

01

212212公式得证，故由“夹逼性”知：，故又存在单调有界定理”知：单调递增有上界，由“即且，显然又代入公式得：据“积分不等式”知：，故由于Wallis n n n n a n n n n n n n n n n a a a a a a k k n n n n n n n n n n n n n n n n n n n n n n n n x x x x x x x x n x x n n n n n n n n n n n n n n

k n n

n n n n ∏∏???∞

=∞

→∞→∞→∞→∞

→+=-+-++?-==+?-=--?-=-+>≤<=-=+?--+?+?-=--?-≤≤+?-?--≤?-≤+≤≤≤≤∈?≤≤∈?π

ππππ

π

π

πN π效率很低，必须通过多次（往往达成百上千亿次）实验方可得到较精确的结果. （五）利用“Wallis 公式”计算π值

【问题】试利用积分?=????

???+=+=?-=?,2,1,01

2 , !)!12(!)!2(2 ,2!)!2(!)!12(d sin 20k k n k k k

n k k x x n ，ππ

推导公

式：

??+?-??????=

122122543442122n n

n n π

并利用该公式计算π值，观察计算效率.

【解】

为利用Wallis 公式计算π值，并观察其计算效率，现编写M 文件如下：

运行结果如下：

图5.1

【答】图5.1表明，利用“Wallis 公式”计算π值，每增加一位有效数字，所需的项数n 就要提高一个数量级. 因此，“Wallis 公式”计算π的效率较低，大致和反正切函数公式

1arctan 4

=π

的幂级数展开式的计算效率处于同一数量级.

（六）e 的计算

【问题】试利用多种方法计算自然对数的底数e 的近似值. 【解】以下四种方法均可用于计算e 值.

(1)利用公式1：n n n

e )1

1(lim +=∞→.

证：这是e 的定义式，故证明从略. 现编写M 文件如下：