怎样计算圆周率

用Matlab计算
创建m文件 calpi1.m，内容如下： function y=calpi1(k) for n=1:k a(n)=(-1).^(n-1)./(2*n-1); end; 4*sum(a)
在命令窗口中输入如下命令： >> calpi1(1000) ans = 3.14059265383979 >> calpi1(10000) ans = 3.14149265359003 >> calpi1(15000) ans = 3.14152598692319 >> calpi1(20000) ans = 3.14154265358982
在命令窗口中输入如下命令： >>calpi4(10000) ans = 3.13240000000000 >>calpi4(50000) ans = 3.14728000000000 >>calpi4(100000) ans = 3.14608000000000
其他方法
1/π 的展开式 Ramanujan 公式
2
相应正多边形面积
a 1 1 a S n +1 = OC × AD = ⋅ n = n 2 2 2 4
π 的值
π ≈ 6 ⋅ 2 n +1 ⋅ S n +1 = 3 ⋅ 2 n an
(刘徽计算到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);
在命令窗口中输入如下命令： >> digits(30) %保留小数点后30位 >> calpi2(10) ans = 3.14159257960635063255949717131 >> calpi2(20) ans = 3.14159265358975625659354591335 >> calpi2(50) ans = 3.14159265358979323846264338328
2 2 ∞ (4n)!(1103 + 26390n) = ∑ π 9801 n =1 (n!) 2 396 4 n 1
算术几何平均值迭代法
an + bn 1 a0 = 1, b0 = , an +1 = , bn +1 = anbn , M = lim an = lim bn n →∞ n →∞ 2 2
ቤተ መጻሕፍቲ ባይዱ
但你会计算π 的值吗？你又能用几种方法计算？
刘徽割园法
从正六边形开始，逐步求 边长与面积 递推法
设边数为6 ⋅ 2 n 的正多边形边长为an ,
o B C D A
如图 AC 2 = AD 2 + DC 2 = AD 2 + (OC − OD) 2
an 2 ⎛ a ⎞ 2 ⎜1 − 1 − ( n ) 2 ⎟ = 2 − 4 − a n an +1 = ( ) + ⎜ 2 2 ⎟ ⎝ ⎠
1.验证公式 2.利用积分
π
1 1 1 = arctan + arctan + arctan 5 8 4 2
∫
推导公式
π
2 0
(n − 1)!! π sin x dx = ⋅ , n为奇数 n!! 2
n
π
2 2 4 4 2n 2n = ⋅ ⋅ ⋅ ⋅L ⋅ ⋅ L 2 1 3 3 5 2n − 1 2n + 1
1 =
2 2 1 − ∑ 2 n (an − bn ) n =0 ∞
π
2M 2
计算 π 的意义
反映数学和计算技术发展的一个侧面 “历史上一个国家所算得的圆周率的准确程 度，可以作为衡量这个国家当时数学发展水 平的指标。”
22 355 <π < 7 113
3.1415926 ＜ π ＜ 3.1415927 (领先世界900余年)
问题: 能不能算得更快一点、更精确一点？
简单公式
1 1 π arctan + arctan = 2 3 4
1 1 1 3 1 1 5 ( − 1) n − 1 1 2 n − 1 π = 4[ − ( ) + ( ) − L + ( ) +L 2 3 2 5 2 2n − 1 2
1 1 1 3 1 1 5 ( − 1) n −1 1 2 n −1 + − ( ) + ( ) −L + ( ) + L] 3 3 3 5 3 2n − 1 3
Machin公式
π 1 1 = 4 arctan − arctan 5 239 4
用Matlab计算
创建m文件 calpi2.m，内容如下： function y=calpi2(k) for n=1:k a(n)=(-1).^(n-1)*(1/2).^(2*n-1)./(2*n-1)+(-1).^(n1)*(1/3).^(2*n-1)./(2*n-1); end; vpa(4*sum(a))
在命令窗口中输入如下命令： >> digits(30) %保留小数点后30位 >> calpi3(100) ans = 3.14157598692312900467982217378 >> calpi3(500) ans = 3.14159198692313035294887413329 >> calpi3(10000) ans = 3.14159265192314007819618382200
3. 用Monte Carlo 法计算π，除了加大随机数， 在随机数一定时可重复算若干次后求平均值， 看能否求得5位精确数字？ 4. 设计方案用计算机模拟Buffon实验
5. 利用学习过的知识(或查阅资料)，提出其他
计算π的方法(先用你学过的知识证明)，然后 实践这方法. 6.对你在实验中应用的计算π的方法 进行比较讨论
人工计算：实验法→ 几何法→ 分析法 最高记录：808位(1948 ) 计算机方法:
位数 年代 2035 1949 100万 1973 10亿 1989 2061亿 1999 12411万亿 2002
测试或检验超级计算机的各项性能(Super PI) 引发新的概念、方法和思想 ,产生新的问题
作业
Monte Carlo 法
从Buffon落针实验谈起： 纸上一组平行线距离为1， 将长度为1的针多次地扔到 纸上。若扔针次数为m,而其中 针与平行线相交的次数为n Buffon指出：π 的数值与 m/n 有关，他由此 求出π 的近似值为3.142
设计方案
在正方形 0< x <1, 0< y<1 上随机的投大量的点，那么 落在四分之一园内的点数 数m与在正方形内的点数n 之比 m/n 应为这两部分图形 面积之比＝π /4, 故
π ＝4 m/n
计算机模拟：产生区间[0,1]上数目为n 的一组 随机数(x, y)，计算满足 x2 + y2 <1 的点数m
用Matlab计算
创建m文件calpi4.m，内容如下： 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
1 = 1 − x 2 + x 4 − L + ( − 1) n −1 x 2 n − 2 + L 1+ x2
积分导出
x3 x5 x 2 n −1 + − L + ( − 1) n −1 +L arctan x = x − 3 5 2n − 1
取 x=1
π
1 1 1 n −1 = 1 − + − L + ( − 1) +L 4 3 5 2n − 1
用Matlab计算
创建 m 文件 fun.m，内容如下： function y=fun(x) y=4./(1+x.^2); 创建 m 文件calpi3.m，内容如下： 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)))
3. Buffon落针实验中，若扔针次数为m, 而其中针与平行线相交的次数为n，试导出
π 与m/n 的关系
实验任务
1. 用反正切函数的幂级数展开式结合有关公式 求π，若要精确到40位、50位数字，试比较 简单公式和Machin公式所用的项数. 2. 用数值积分计算π，分别用梯形法和Simpson 法精确到10位数字，用Simpson法精确到15 位数字.
哪里有数,哪里就 有美.
- Proclus
知其然,更知其所以 然.
-中国先哲
数学实验
怎样计算 π 的值 ?
上海交大数学系
实际问题
π―圆周率, 我们十分熟悉的常数. 你也许能写出 π = 3.1415926535 用Matlab 容易求出π 到几百位
>> digits(100) >> vpa(pi) ans = 3.1415926535897932384626433832795028841971693 99375105820974944592307816406286208998628034 825342117068
利用数值积分方法
1 A = 4∫ dx = π 2 01 x ＋
1
1 设 y ( x) = 1+ x2
将[0,1]区间n等分，取xk=k/n, yk= 1/ (1+xk2)