哪里有数,哪里就有美.
- Proclus
知其然,更知其所以然.
-中国先哲数学实验上海交大数学系怎样 计算
π 的值?
你也许能写出 π = 3.1415926535
实际问题
π―圆周率,我们十分熟悉的常数.
用 Matlab 容易求出 π 到几百位
>> digits(100)
>> vpa(pi)
但你会计算 π 的值吗？你又能用几种方法计算？
ans =
3.1415926535897932384626433832795028841971693
99375105820974944592307816406286208998628034
825342117068
刘徽割园法
? 从正六边形开始，逐步求边长与面积
? 递推法
,26
n
n
a的正多边形边长为设边数为?
oA
B
C
D
22222
)( ODOCADDCADAC?+=+=如图
2
2
22
1
42)
2
(11)
2
(
n
nn
n
a
aa
a=
+=
+
相应正多边形面积
422
1
2
1
1
nn
n
aa
ADOCS =?=×=
+
? π 的值
n
n
n
n
aS 2326
1
1
=≈
+
+
π
(刘徽计算到96边形面积，得到 π≈3.141)
用Matlab计算
function y=calpi(n)
syms a;
for i=1:n
a=sqrt(2-sqrt(4-a^2));
end
a=subs(a,'a','1');
y=3*2^n*vpa(a,n+5);
n=5 3.14145240288
（ 刘徽 3位 ）
n=11 3.1415926198105522176
（ 祖冲之 7位 ）
n=16 3.1415926535797449417818112
（ 韦达 10位 ）
n=61 3.1415926535897932384626433832754
（ 鲁道夫 29 (34)位 ）
LL+
+?+?=
12
1
)1(
5
1
3
1
1
4
1
n
n
π
LL+
+?+?=
12
)1(
53
arctan
12
1
53
n
xxx
xx
n
n
LL+?+?+?=
+
22142
2
)1(1
1
1
nn
xxx
x
积分导出取 x=1
利用幂级数计算用Matlab计算
function y=calpi1(k)
for n=1:k
a(n)=(-1).^(n-1)./(2*n-1);
end;
4*sum(a)
创建 m文件 calpi1.m,内容如下,
>> calpi1(1000)
ans =
3.14059265383979
>> calpi1(10000)
ans =
3.14149265359003
在命令窗口中输入如下命令：
>> calpi1(15000)
ans =
3.14152598692319
>> calpi1(20000)
ans =
3.14154265358982
问题,能不能算得更快一点、更精确一点？
43
1
arctan
2
1
arctan
π
=+
LL+
+?+?=
12
1
53
)
2
1
(
12
)1(
)
2
1
(
5
1
)
2
1
(
3
1
2
1
[4
n
n
n
π
Machin公式
4239
1
arctan
5
1
arctan4
π
=?
简单公式
])
3
1
(
12
)1(
)
3
1
(
5
1
)
3
1
(
3
1
3
1
12
1
53
LL+
+?+?+
n
n
n
用Matlab计算
function y=calpi2(k)
for n=1:k
a(n)=(-1).^(n-1)*(1/2).^(2*n-1)./(2*n-1)+(-1).^(n-
1)*(1/3).^(2*n-1)./(2*n-1);
end;
vpa(4*sum(a))
创建m文件 calpi2.m,内容如下,
>> digits(30) %保留小数点后 30位
>> calpi2(10)
ans =
3.14159257960635063255949717131
在命令窗口中输入如下命令：
>> calpi2(20)
ans =
3.14159265358975625659354591335
>> calpi2(50)
ans =
3.14159265358979323846264338328
π==
∫
dx
x
A
1
0
2
1
1
4
＋
2
1
1
)(
x
xy
+
=设 将[0,1]区间 n等分，取 x
k
=k/n,
])(2[
2
0121 nn
yyyyy
n
A +++++=?
L梯形法
)(2)[(
6
1
224220?
+++++?
mm
yyyyy
m
SimpsonL法
y
k
= 1/ (1+x
k
2
)
)](4
1231?
利用数值积分方法
++++
m
yyyL
用Matlab计算
function y=fun(x)
y=4./(1+x.^2);
function y=calpi3(k)
for n=1:k-1
a(n)=2*fun(n/k);
end;
vpa(1/(2*k)*(sum(a)+fun(0)+fun(1)))
创建 m 文件 fun.m，内容如下,
创建 m 文件 calpi3.m,内容如下,
>> digits(30) %保留小数点后 30位
>> calpi3(100)
ans =
3.14157598692312900467982217378
在命令窗口中输入如下命令：
>> calpi3(500)
ans =
3.14159198692313035294887413329
>> calpi3(10000)
ans =
3.14159265192314007819618382200
针与平行线相交的次数为 n
Monte Carlo 法从Buffon落针实验谈起,
纸上一组平行线距离为1，
将长度为1的针多次地扔到纸上。若扔针次数为 m,而其中
Buffon指出,π 的数值与 m/n 有关，他由此求出 π 的近似值为3.142
设计方案
π ＝4 m/n
计算机模拟：产生区间 [0,1]上数目为 n 的一组在正方形 0< x <1,0< y<1
上随机的投大量的点，那么落在四分之一园内的点数数 m与在正方形内的点数 n
之比 m/n 应为这两部分图形面积之比 ＝ π /4,故随机数 (x,y)，计算满足 x
2
+ y
2
<1 的点数 m
用Matlab计算
function y=calpi4(k)
m=0;
for n=1:k
if rand(1)^2+rand(1)^2<=1
m=m+1;
end;
end;
4*m/k
创建 m文件calpi4.m，内容如下,
>>calpi4(10000)
ans =
3.13240000000000
>>calpi4(50000)
ans =
3.14728000000000
>>calpi4(100000)
ans =
3.14608000000000
在命令窗口中输入如下命令：
其他方法
? 1/π 的展开式
Ramanujan 公式
∑
∞
=
+
=
1
42
396)!(
)263901103()!4(
9801
221
n
n
n
nn
π
? 算术几何平均值迭代法
n
n
n
n
nnn
nn
n
baMbab
ba
aba
∞→∞→
++
===
+
=== limlim,,
2
,
2
1
,1
1100
2
0
22
2
)(21
1
M
ba
n
nn
n
∑
∞
=

=
π
计算 π 的意义
? 反映数学和计算技术发展的一个侧面
“历史上一个国家所算得的圆周率的准确程度，可以作为衡量这个国家当时数学发展水平的指标。,
3.1415926 ＜ π ＜ 3.1415927
113
355
7
22
<<π
(领先世界900余年 )
? 引发新的概念、方法和思想,产生新的问题位数 2035 100万 10亿 2061亿 12411万亿年代 1949 1973 1989 1999 2002
人工计算：实验法 → 几何法 → 分析法计算机方法,
最高记录：808位(1948 )
? 测试或检验超级计算机的各项性能( Super PI)
作业
8
1
arctan
5
1
arctan
2
1
arctan
4
++=
π
1.验证公式
2.利用积分
dxx
n
∫
2
0
sin
π
为奇数n
n
n
,
2!!
!)!1( π
=
推导公式
LL
12
2
12
2
5
4
3
4
3
2
1
2
2 +
=
n
n
n
nπ
3,Buffon落针实验中，若扔针次数为 m,
而其中针与平行线相交的次数为 n,试导出
π 与 m/n 的关系实验任务
1,用反正切函数的幂级数展开式结合有关公式
2,用数值积分计算 π，分别用梯形法和 Simpson
简单公式和 Machin公式所用的项数,
求 π，若要精确到40位、50位数字，试比较法精确到 10位数字，用 Simpson法精确到15
位数字,
看能否求得 5位精确数字？
4,设计方案用计算机模拟 Buffon实验
3,用 Monte Carlo 法计算 π，除了加大随机数，
在随机数一定时可重复算若干次后求平均值，
5,利用学习过的知识 (或查阅资料 )，提出其他计算 π的方法 (先用你学过的知识证明 )，然后
6.对你在实验中应用的计算 π的方法进行比较讨论实践这方法,
谢谢各位！