1数学建模实验-圆周率的计算

Ä 1630年,最后一位用古典方法求π的人 年 最后一位用古典方法求π 格林伯格也只求到了π的第 位小数 格林伯格也只求到了 的第39位小数 的第
8-23
分析方法
从十七世纪中叶起， 从十七世纪中叶起，人们开始用更先进的 分析方法来求π的近似值，其中应用的主 分析方法来求π的近似值， 要工具是收敛的无穷乘积和无穷级数， 要工具是收敛的无穷乘积和无穷级数，在 本节中我们将介绍一些用此类方法求π 本节中我们将介绍一些用此类方法求π近 似值的实例。 似值的实例。
19-23
设计方案
在正方形 0< x <1, 0< y<1 上随机的投大量的点，那么 落在四分之一园内的点数 数m与在正方形内的点数n 之比m/n应为这两部分图形 面积之比＝π/4,故 π＝4 m/n 计算机模拟：产生区间[0,1]上数目为n的一组 随机数(x,y)，计算满足x2+y2<1的点数m
>>ans=
3.1415925796063512110
>> vpa(4*symsum((-1)^(k-1)/(2*k-1)*(1/2)^(2*k1),1,20)+4*symsum((-1)^(k-1)/(2*k-1)*(1/3)^(2*k-1),1,20),20)
>>ans=
3.1415926535897574098
11-23
麦琴( Ä 麦琴(Machin)给出 )
π
1 1 = 4 arctan − arctan (Machin公式) 公式) 公式 4 5 239
此式求得了π的第 位小数且全部正确 此式求得了 的第100位小数且全部正确 的第
12-23
方法1 利用幂级数表达式
1 2 4 n−1 2n−2 =1− x + x −L+ (−1) x +L 2 1+ x
6边形 边形
12边形 边形
24边形 边形
圆
6-23
96边形和圆外切 Ä 阿基米德曾用圆内接 96边形和圆外切 96边形夹逼的方法证明了 96边形夹逼的方法证明了
223 22 <π < 71 7
由 sin θ < θ < tan θ 和 θ = π 96 导出
公元5世纪 世纪， Ä 公元 世纪，祖冲之指出
数学建模实验（ 数学建模实验（一）
怎样计算 π 的值 ?
贵州师范大学 周晓军
1-23
哪里有数,哪里 就有美.
- Proclus
Байду номын сангаас
知其然,更知其所 以然.
-中国先哲
2-23
实际问题
π―圆周率, 我们十分熟悉的常数. 你也许能写出 π = 3.1415926535 用Matlab容易求出π到几百位
>>vpa(pi,100) >>ans=
10-23
Ä 在中学数学中证明过下面的等式
π
1 1 = arctan1 = arctan + arctan 4 2 3
C A B D
1 1 arctan 和 arctan 3 2 的展开式的收敛速度
都比 arctan 1 快得多
左边三个正方形 组成的矩形中， 组成的矩形中， 由∠A = ∠B + ∠C 和∠C = ∠D 可得
16-23
方法2 利用数值积分
1 设 y(x) = 1+ x2
1 A = 4∫ dx = π 2 01 x ＋
1
将[0,1]区间 n 等分，取 xk=k/n, yk= 1/ (1+ k2) (1+x
2 梯形法⇒ A = [2( y1 + y2 +L+ yn−1) + y0 + yn ] n
2 Simpson法⇒ [( y0 + y2m ) + 2( y2 + y4 +L+ y2m−2 ) 3m + 4( y1 + y3 +L+ y2m−1)]17-23
x3 x5 x2n−1 n−1 arctan x = x − + −L+ (−1) +L 3 5 2n −1
π
1 1 1 n−1 =1− + −L+ (−1) +L 4 3 5 2n −1
13-23
用Matlab计算 Matlab计算 lab
>>syms k >>vpa(4*symsum((-1)^(k-1)/(2*k-1),1,1000),20) >>ans=
14-23
问题: 能不能算得更快一点、更精确 一点？ 简单公式
1 1 π arctan + arctan = 2 3 4
1 1 1 3 1 1 5 (−1 n−1 1 2n−1 ) π = 4[ − ( ) + ( ) −L+ ( ) −1 +L 2 3 2 5 2 2n −1 2
1 1 13 1 15 (−1 n−1 1 2n−1 ) + − ( ) + ( ) −L+ ( ) +L ] 3 3 3 5 3 2n −1 3
5.利用学习过的知识(或查阅资料)，提出其他
计算π的方法
23-23
9-23
Ä 1656年，沃里斯(Wallis)证明 年 沃里斯( ) 2 2 4 4 6 6 π = 2 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ L 1 3 3 5 5 7 ∞ 2k 2k = 2 ⋅∏ ⋅ 2k + 1 k =1 2 k − 1 请编程序实现上述求法
3.1405926538397929260
>>vpa(4*symsum((-1)^(k-1)/(2*k-1),1,10000),20) >>ans=
3.1414926535900432385
>>vpa(4*symsum((-1)^(k-1)/(2*k-1),1,100000),20) >>ans=
3.1414926535900432385
4-23
古 典 方 法 分 析 方 法 其 它 方 法 Ä 概率方法 Ä 数值积分方法
5-23
古典方法
的近似值呢？显然， 用什么方法来计 算π的近似值呢？显然，不可能仅根 据圆周率的定义，用圆的周长去除以直径。起先， 据圆周率的定义，用圆的周长去除以直径。起先，人们 采用的都是用圆内接正多边形和圆外切正多边形来逼近 的古典方法。 的古典方法。
Machin公式
1 1 π 4arctan − arctan = 5 239 4
15-23
再用Matlab
>>syms k >> vpa(4*symsum((-1)^(k-1)/(2*k-1)*(1/2)^(2*k1),1,10)+4*symsum((-1)^(k-1)/(2*k-1)*(1/3)^(2*k-1),1,10),20)
Matlab计算：
>>syms t >>f=inline(1/(1+t^2)); >>n=1000;x=0+(0:n)/n; >>vpa(2*(2*sum(f(x(2:n)))+f(x(1))+f(x(n+1)))/n) ans = 3.1415924869231237615 >>m=10;X=0+(0:2*m)/(2*m); >>vpa(2*(2*sum(f(X(3:2:2*m1)))+4*sum(f(X(2:2:2*m)))+f(X(1))+f(X(2*m+1)))/(3*m ),20) ans = 18-23 3.1415926529697850711
1 1 π arctan + arctan = 2 3 4
1 1 π 4arctan − arctan = 5 239 4
2. 用数值积分计算π，分别用梯形法和Simpson 法精确到10位数字，用Simpson法 精确到15位数字.
要求：写出实验报告（包括：算法、程 序、结果）
22-23
3. 用Monte Carlo 法计算π，除了加大随机数， 在随机数一定时可重复算若干次后求平均值， 看能否求得5位精确数字？ 4. 设计方案用计算机模拟Buffon实验
方法3 Monte Carlo 法
从Buffon落针实验谈起： 平行线距离为1，针长度为1; 设针中心到较近平行线距离 为y ，针与平行线夹角θ , 针 与平行线相交充要条件为
1 π y ≤ sin θ， ≤θ ≤ 0 2 2
1 2
θ
D o
π
2
(左图正弦曲 线下方面积D)
次数很大,落针应均匀分布，落针中心在D与 总数之比为D面积与总面积之比为2/π
3.1415926 < π < 3.1415927
比西方得到同样结 果几乎早了1000年 果几乎早了 年
7-23
十五世纪中叶，阿尔·卡西给出 卡西给出π Ä 十五世纪中叶，阿尔 卡西给出π的16 位小数， 位小数，打破了祖冲之的纪录 Ä 1579年,韦达证明 年
3.1415926535 < π < 3.1415926537
3.1415926535897932384626433832795028841971 6939937510582097494459230781640628620899 8628034825342117068 但是你会计算π的值吗？你又能用几种 方法计算π的值？
3-23
π的计算 的计算
Ä 圆周率是人类获得的最古老的数学概念 之一，早在大约3700年前（即公元前1700 之一，早在大约3700年前（即公元前1700 3700年前 年左右） 256/81（ 年左右）的古埃及人就已经在 用256/81（ 约3.1605）作为π的近似值了。几千年来 3.1605）作为π的近似值了。 ，人们一直没有停止过求π的努力。 人们一直没有停止过求π的努力。
20-23
Matlab模拟：
m=0;n=100000; x=rand(1,n);y=rand(1,n); for i=1:n if x(i)^2+y(i)^2<1 m=m+1; end; end; 4*m/n
21-23
实验任务
1. 用反正切函数的幂级数展开式结合有关公式 求π，若要精确到以40位、50位数字，试比较 简单公式和Machin公式所用的项数.