计算pai(07)

在命令窗口中输入如下命令： >> digits(30) %保留小数点后30位 >> calpi3(100) ans = 3.14157598692312900467982217378 >> calpi3(500) ans = 3.14159198692313035294887413329 >> calpi3(10000) ans = 3.14159265192314007819618382200
在命令窗口中输入如下命令： >> digits(30) %保留小数点后30位 >> calpi2(10) ans = 3.14159257960635063255949717131 >> calpi2(20) ans = 3.14159265358975625659354591335 >> calpi2(50) ans = 3.14159265358979323846264338328
π
2 2 4 4 2n 2n = ⋅ ⋅ ⋅ ⋅L ⋅ ⋅ L 2 1 3 3 5 2n − 1 2n + 1
3. Buffon落针实验中，若扔针次数为m, 而其中针与平行线相交的次数为n，试导出
π 与m/n 的关系
实验任务
1. 用反正切函数的幂级数展开式结合有关公式
求π，若要精确到40位、50位数字，试比较 简单公式和Machin公式所用的项数.
>> digits(100) >> vpa(pi) ans = 3.1415926535897932384626433832795028841971693 99375105820974944592307816406286208998628034 825342117068
但你会计算π 的值吗？你又能用几种方法计算？
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))
问题: 能不能算得更快一点、更精确一点？ 简单公式
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) n −1 1 2 n −1 1 1 1 3 1 1 5 + − ( ) + ( ) −L + ( ) + L] 3 3 3 5 3 2n − 1 3
哪里有数,哪里就 有美.
- Proclus
知其然,更知其所以 然.
-中国先哲
数学实验
怎样计算 π 的值 ?
上海交大数学系
说明： 本实验内容参考了中国科技大学 李尚志等编写的数学实验教材
Fra Baidu bibliotek际问题
π―圆周率, 我们十分熟悉的常数. 你也许能写出 π = 3.1415926535 用Matlab 容易求出π 到几百位
22 355 <π < 7 113
3.1415926 ＜ π ＜ 3.1415927 (领先世界900余年)
人工计算：实验法→ 几何法→ 分析法 最高记录：808位(1948 ) 计算机方法:
位数 年代 2035 1949 100万 1973 10亿 1989 2061亿 1999 12411万亿 2002
用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)))
计算π的方法(先用你学过的知识证明)，然后 实践这方法.
6.对你在实验中应用的计算π的方法 进行比较讨论
谢谢各位！
方法2
利用数值积分
1 设 y ( x) = 1+ x2
1 A = 4∫ dx = π 2 01 x ＋
1
将[0,1]区间n等分，取xk=k/n, yk= 1/ (1+xk2)
2 梯形法 ⇒ A = [2( y1 + y2 + L + yn −1 ) + y0 + yn ] n 1 Simpson法 ⇒ [( y0 + y2 m ) + 2( y2 + y4 + L + y2 m − 2 ) 6m + 4( y1 + y3 + L + y2 m −1 )]
用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
方法1
利用幂级数表达式
1 2 4 n −1 2 n − 2 = 1 − x + x − L + ( −1) x +L 2 1+ x
x x n −1 x + − L + ( − 1) +L arctan x = x − 3 5 2n − 1
3 5 2 n −1
π
1 1 1 n −1 = 1 − + − L + ( − 1) +L 4 3 5 2n − 1
方法3
Monte Carlo 法
从Buffon落针实验谈起： 纸上一组平行线距离为1， 将长度为1的针多次地扔到 纸上。若扔针次数为m,而其中 针与平行线相交的次数为n Buffon指出：π 的数值与 m/n 有关，他由此 求出π 的近似值为3.142
设计方案
在正方形 0< x <1, 0< y<1 上随机的投大量的点，那么 落在四分之一园内的点数 数m与在正方形内的点数n 之比 m/n 应为这两部分图形 面积之比＝π /4, 故
在命令窗口中输入如下命令： >>calpi4(10000) ans = 3.13240000000000 >>calpi4(50000) ans = 3.14728000000000 >>calpi4(100000) ans = 3.14608000000000
计算 π 的意义
反映数学和计算技术发展的一个侧面 “历史上一个国家所算得的圆周率的准确程 度，可以作为衡量这个国家当时数学发展水 平的指标。”
测试或检验超级计算机的各项性能(Super PI) 引发新的概念、方法和思想 ,产生新的问题
作业
1.验证公式 2.利用积分
π
1 1 1 = arctan + arctan + arctan 5 8 4 2
∫
推导公式
π
2 0
(n − 1)!! π sin x dx = ⋅ , n为奇数 n!! 2
n
2. 用数值积分计算π，分别用梯形法和Simpson
法精确到10位数字，用Simpson法精确到15 位数字.
3. 用Monte Carlo 法计算π，除了加大随机数，
在随机数一定时可重复算若干次后求平均值， 看能否求得5位精确数字？
4. 设计方案用计算机模拟Buffon实验
5. 利用学习过的知识(或查阅资料)，提出其他
π ＝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