数学实验之8.pi的近似计算
东华理工大学《数学实验》报告
东华理工大学《数学实验》报告学号姓名成绩实验名称定积分的近似计算一、实验目的及要求掌握Excel产生均匀随机数的方法，了解Excel的一些常用数学函数；了解通过蒙特卡罗（Monte-Carlo）方法模拟随机事件的方法。
二、实验环境多媒体计算机WindowsXP操作系统Excel软件三、问题背景：Monte-Carlo方法的基本想法是首先建立与所解决问题有相似性的概率模型，利用这种相似性把这种概率模型的某些特征与数学分析问题的解答联系起来，然后对模型进行随机模拟或统计抽样，再利用所得结果求出这些特征的统计值做为原来的分析问题的近似解，它在解决数学问题中有广泛的应用。
下面我们利用Monte-Carlo方法来计算定积分。
考虑定积分，不妨设当这时等于由曲线与直线x=a，x=b所围成的区域G的面积，为求它的值，我们在G外作一个矩形（四边所在直线为x=a,x=b,y=c,y=d）。
若在矩形区域内随机投掷一点，落点的两个坐标是独立的且在区间[a,b],[c,d]上是均匀分布，那么每个点落在曲线y=f(x)以下的区域G内的概率p显然等于这个区域的面积与矩形面积之比。
设想在矩形内随机投掷N个点，若这N个点中有n个点落入G中，则G的面积的估计值为，此即为定积为的近似值I。
四、实验内容利用Monte-Carlo方法计算：的近似值五、实验过程：（1）选定矩形区域（2）确定实验次数，即N值，不妨取N=100（3）选定单元格A1,输入x,选定A2,输入函数(4)选定A2,按”ctrl+C”,选定A3～A101,按”ctrl+V”。
此时A2～A101为区间[0，2]的均匀随机数（5）用与（3）（4）步同样的方法在B2～B101产生[0，1]的均匀随机数作为y(6)在C2输入”=if(b2<sqrt(1-a2^2/4),0,1)”, 再选定C2,按”ctrl+C”,选定C3～C101,按”ctrl+V”(7)选定D1,输入函数“=frequency(C2:C101,0),些时对C3～C101列中等于0的数作频数统计，记此数值为n,也即计算有n个点落在区域G中（8）的近似值为六、实验记录N n 的近似值100100010000七、思考与深入：（1）异常情况分析：（2）根据所得结果猜测椭圆的面积是。
实验8 无穷级数与数值计算
西安理工大学应用数学系
（2）相应的Matlab代码为 syms x taylor((1+x)^0.5, x, 7, 1)
运行结果为 ans =2^(1/2) + (2^(1/2)*(x - 1))/4 - (2^(1/2)*(x - 1)^2)/32 +
(2^(1/2)*(x - 1)^3)/128 - (5*2^(1/2)*(x - 1)^4)/2048 + (7*2^(1/2)*(x - 1)^5)/8192 - (21*2^(1/2)*(x 1)^6)/65536
(n
3
1)!
105,取
n8
即可。
西安理工大学应用数学系
同理取 x 1
2
,令误差
2n
1
3 (n
1)!
105
,取 n 6 即可。
计算e的Matlab代码为
syms x y y=taylor(exp(x),x,9)
a=subs(y,x,'1') %将符号表达式y中的变量x换成变量1 e=vpa(a,10) % 将符号表达式a转化成10位精度的近似解
1 x 1 1 x 1 x2 1 3 x3 (1)n1 1 35 (2n 3) xn （6.6）
2 24 246 取
2 46 2n
于是有
, 代入(6.6), 有
1 x 1 1 x 1 x2 13 x3 (1)n1 1 3 5 (2n 3) xn
西安理工大学应用数学系
hold on for i=1:length(x)
y(i)=x(i)-x(i)^3/6+x(i)^5/120; plot(x(i),y(i),'.'),pause(.005) end hold on for i=1:length(x) y(i)=x(i)-x(i)^3/6+x(i)^5/120-
实验二 素数问题
练习九 在二维坐标面上标出点列 ( n, π ( n)), n = 1,2,L, N ,(取 不同的 N ,如1000,10000等).也可以用折线将点连 起来.观察 π (n) 趋于无穷的趋势,将它同 y = x , y = x 比较,你会有什们结论?类似地观察点列 ( n, π ( n) / n) 和 ( n, π ( n) / n ) 以及 ( n, π ( n) /( n / Log( n))) .你能据此猜 测趋于无穷的极限的阶吗?
Mersenne素数是极其稀少的.借助大型计算机, 截止96年11月,数学家仅发现了34个Mersenne素数. 它们对应的 n是：
2,3,5,7,13,17,19,31,61,89,107,127,521,607,1279,2203, 2281,3217,4253,4423,9689,9941,11213,19937,21701, 23209,44497,86243,123049,216091,756839,859433, 1257787,1398269.
早在十七十八世纪，数学家Fermat和Ruler等 就研究过这类公式．1640年Fermat在给Mersenne 的信中指出，对所有的整数 n, Fn = 2 2 + 1 永远是素数. 的确F0 = 3, F1 = 5, F2 = 7, F3 = 257, F4 = 65537 ,都是素数. 然而,1732年,大数学家Ruler指出,F5 = 4294967297 不是素数,他并且找到了F5 的因子分解.此后,人们分 别证明了 F6 与 F7都是合数,并得到了它们的素因子分 解．实际上，有人猜测 Fn 当 n > 4时都是合数.
进一步的问题 关于素数,存在许许多多富有挑战性的问题,吸 引众多的数学家及业余爱好者．下面我们介绍几个 供有兴趣的同学参阅. Goldbach猜想 Bertrand猜想 大整数的素因子分解 完全数 孪生素数 青一色数的素性
数学实验之素数
素数探秘一、实验目的素数(Prime)是构造所有数的“基本材料”，犹如化学上的化学元素和物理学中的基本粒子，有关素数的许多看似简单却极富刺激性的奇妙问题，向一代代数学家提出了挑战，始终吸引着他们的目光.本实验通过对若干素数问题的基本认识，激发学生对数论中研究课题的兴趣，并体会探索数学奥妙的艰巨性.二、实验内容1素数的判别欧几里得给出素数的定义：如果一个大于1的自然数只能被1及它本身整除，则称该数为素数. 否则称为合数(1即不是素数也不是合数).(1) Mathematica的素数函数与筛法Mathematica系统提供了两个常用的与素数有关的函数：Prime[n]返回从第一个素数2数起的第n个素数.PrimeQ[n] 判断自然数n是否为素数，是则返回True，否则返回False.例如，要输出第25个素数并判断67013是否为素数，输入{Prime[25],PrimeQ[67021]}输出结果为{97,True}.使用系统函数输出某个指定范围内的所有素数，只要定义如下的函数即可.bprime[m_Integer,n_Integer]:=Select[Table[k,{k,m,n}],PrimeQ]例如，要输出100—115之间的素数，这些上面的命令后，输入bprime[15,30]输出结果为{101,103,107,109,113}最早人们是如何寻找素数呢？2000多年前，希腊学者埃拉托色尼(Eratosthenes 公元前约284-192年)给出了一个寻找素数的简便方法—筛法：写下从2、3、…、N，注意到2是一个素数，划去后面所有2的倍数，越过2，第一个没有被划去的数是3，它是第二个素数，接下来再划掉所有3的倍数，3之后没有被划去的数是5，然后再划掉除5外所有5的倍数，以此类推. 显然，划掉的都是较小整数的倍数，它们都不是素数，都被筛掉了. 而素数永远不会被筛掉，它们就是要寻找的不超过N的所有素数.问题1 编写一个利用筛法寻找素数的程序，并输出2—100内的所有素数.解输入Sieve[n_Integer]:=Module[{t,i,temp},t=Range[2,n];For[i=2,i<=Floor[Sqrt[n]],i++,t=Select[t,Mod[#,i]!=0&];If[PrimeQ[i],t=Prepend[t,i]]];Sort[t]]Sieve[100]执行后输出100以内的25个素数：{2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97}(2) 素数判断判断一个正整数n 是否是素数，只要将n 分别除以2n 都不能被整除，则n 为素数. 这种判断一个整数是否为素数的方法，称为试除法，我们作下面的实验：问题2 自定义一个用试除法判断素数的函数，判断67021是否是素数.解 将n 分别除以2n 都不能被整除，则2mod(, )0i n i =≠，其中，为. 下面的素数判断函数，当n 是素数，返回“True ”；若不是素数，返回“False ”. 输入 Prnum[n_Integer]:=If[Product[Mod[n ,i],{i,2,Floor[Sqrt[n ]]}]!=0,True,False]Prnum[67021]执行后得到True.问题3 输出某个指定范围内的所有素数.解 可以定义如下的函数 BetweenPrime[m_Integer,n_Integer]:=Module[{t ={}},For[k =m ,k <=n ,k ++,If[Product[Mod[k ,i ],{i ,2,Floor[Sqrt[k]]}]==0,Continue[],t =Append[t ,k ]]];t ]BetweenPrime[50,100]执行后得到{53,59,61,67,71,73,79,83,89,97}(3) 素数的无穷性关于素数首先提出的问题是有没有最大的素数，或者素数是否有穷尽？欧几里得在《几何原本》中提出如下命题：“素数的数目大于任何指定的素数集合中的素数的数目”. 即任何有限的素数集合都不可能包含全部的素数，素数有无穷多.欧几里得使用非常优美的数学推理方法证明了这个命题. 首先，设所有的素数按照自小到大的次序列成：12n p p p <<<， (1) 证明的诀窍在于研究数：121n N p p p =+. (2)显然，n N p >. 若N 是一个素数，则说明在n p 之后确实还有一个素数，它不属于原有的素数集合. 另一方面，若N 是一个合数，则它必定能够被某个素数整除，不妨设这个素数为p ，由于N 除以 (1,2,,)i p i n =所得的余数都是1，所以，素数p 必为不同于各个i p 的另一个素数. 因此，不论N 是否是素数，我们总能找到一个新素数. 例如对3，7，13，371312742137N =⋅⋅+==⋅，137是一个新素数.问题4 计算121n n N p p p =+，3,4,n =，判断它们是否是素数，如果不是素数，求出其素因子. 解 为使程序简单起见，假定 (1,2,,)i p i n =为前n 个素数，输入Pnum[n_Integer]:=Module[{s =2,t ={}},For[i =2,i<=n,i ++,s *=Prime[i];t =Append[t,Apply[Times,s]+1]];GridBox[Transpose[{t ,PrimeQ/@t ,FactorInteger/@t }],RowLines →True,ColumnLines →True]//DisplayForm ]为了使得输出结果便于观察，程序中使用GridBox[]函数. 取一个较小的n ，如n =8，调用上述程序得到Pnum[8] //DisplayForm 7True7,131True31,1211True211,12311True2311,130031False59,1,509,1510511False19,1,97,1,277,9699691False 347,1,27953, 从输出结果看到前四个结果2317⋅+=，…，23571112311⋅⋅⋅⋅+=都是相应素数集合之外的新素数，而接下来的23571113130031⋅⋅⋅⋅⋅+=等三个数都是个合数，30031的两个素因子59和509都是新素数.2 素数公式与梅森素数(1) 多项式形的素数生成公式利用上述方法判断一个整数是否是素数以及寻找新的素数，对于不太大的n ，看起来并不难. 但是，要判断一个非常大的数是否是素数，可就不简单了. 人们一直试图找到一个能生成所有素数的公式，十七世纪 法国数学家费马(Fermat )曾断言，对任意正整数n ，221n n F =+永远是一个素数. 不难计算对4n ≤，费马的推测确实是正确的. 1732年欧拉提出反例，指出325216416700417F =+=⨯是一个合数，费马的断言是错误的. 此后，有人甚至推断当5n ≥，n F 都是合数.1772年欧拉给出了一个公式，对于整数n ，多项式的值是素数. 然而这个推断也是错误的，事实上，当40n =时，2240404141++=不是素数.计算发现多项式241n n -+、217n n ++(勒让德Legendre )、272491n n ++等，对于某些整数，其值为素数，对于另一些整数，其值为合数.问题5 计算241n n ++，0,1,50n =，判断对于哪些n (n 可取负整数)它们是素数，如果不是素数，求出其素数因子.解 输入f[x_Integer]:=x ^2+x +41Table[{k ,f[k ],PrimeQ[f[k ]],FactorInteger[f[k ]]},{k,0,50}]//ColumnForm执行后得到(仅列出部分输出结果){0,41,True ,{{41,1}}}...........................{39,1601,True ,{{1601,1}}}{40,1681,False,{{41,2}}}{41,1763,False,{{41,1},{43,1}}}...........................{49,2491,False ,{{47,1},{53,1}}}{50,2591,True ,{{2591,1}}}(2) 梅森素数历史上，关于梅森素数的研究也是用公式寻找素数的一个实例：1644年法国神父梅森(M .Mersenne )指 出，对于素数p =2，3，5，7，13，17，19，31，67，127和257，21p -都是素数，而对于其它小于257的素数p ，21p -都是合数. 后来证明梅森的这个推断并不正确，但在17世纪时能得出这一结论，也实属不易. 如今我们把形如21n -的数称为梅森数，若p 为素数且21p -也是素数，称为梅森素数，记为21p p M =-.对任何素数p ，如何判断p M 是不是素数？1930年数学家卢卡斯(E .Lucas )和莱默(D .Lehmer )提出一个有效的检验法：设p 3≥为一个素数，检验法的伪代码为：赋值:4s =；循环变量i 从3到p ，计算2: 2 mod 21p s s =--；若0s = 则21p -是素数，否则是合数. 问题6 编写一个判断p M 是否是梅森素数的程序，检验梅森的结论.解 由于23M =是显然的素数，不用检验. 定义如下检验程序： Mersenne[p_]:=Module[{s =4,t },m =2^p -1;For[i =3,i <=p ,i ++,s =Mod[s ^2-2,m ]];If[s ==0,t =True,t =False];t] (*函数的自变量p 必须是素数*)S ={3,5,7,13,17,19,31,67,127,257}Table[Mersenne[S[[k ]]],{k,Length[S]}] (*调用程序计算*)执行后得到{True ,True ,True ,True ,True ,True ,True ,False ,True ,False}由此结果知道67M 和257M 不是素数. 此外，不超过257的素数共有55个，其中p =61，89和107对应的61M 、89M 和107M 也都是素数，可惜梅森并未检索到.梅森素数是很稀少的，因为随着p 增大，21p -增长速度非常快，现在人们借助于计算机寻找21p -型的梅森素数，旨在发现最大的素数，迄今为止，仅仅找到了43个梅森素数，2005年12月美国数学家柯蒂斯·库珀发现的最大的梅森素数为3040245721-，它有9152052位数.数学家一直努力寻找生成素数的公式，但截至目前，并没有一个函数或是方程式可以有效地产生所有的素数.3 素数分布素数在自然数中的分布是很不规则的，但是人们发现，随着数的增大，素数变得越来越稀疏，十八世纪以来，素数的统计分布问题成为数学家们研究的一个重要课题，为此引进素数个数函数()x π，它是不超过实数x 的素数的个数，例如，(10)4π=，(100)25π=，()x π是一个单调增函数，由欧几里得定理x →+∞时，()x π→+∞.而欧拉又发现()lim 0x x x π→∞=.我们通过如下的实验观察这一事实.问题7 编写程序定义一个统计某个范围内素数个数的函数，计算100以内、100—200、...、900—1000之间素数的个数，再统计10000至100000内，每间隔10000的范围内的素数个数.解 定义Pinum[k_,m_]:=Module[{t =0},For[j =k ,j <=m ,j ++,If[PrimeQ[j ],t +=1]];t ];Table[Pinum[k ,k +100],{k ,0,1000,100}] (*以100为间隔的范围内素数个数*)执行后得到{25,21,16,16,17,14,16,14,15,14,16}再输入如下的命令，输出10000到100000内，每间隔10000的范围内素数的个数.Table[{k ,”-”,k +10000,r =Pinum[k ,k +10000]},{k,10000,100000,10000}]//TableForm执行后得到如下结果://TableForm10000 — 20000 103320000 — 30000 98330000 — 40000 95840000 — 50000 93050000 — 60000 92460000 — 70000 87870000 — 80000 90280000 — 90000 87690000 — 100000 879100000 — 110000 861Mathematica 系统提供了计算不超过x 的素数个数的函数：PrimePi[x] 返回不超过x 的素数个数.问题5.3.7的第二段程序若改为如下形式，输出结果完全相同. Table[{k ,r =PrimePi[k ],r /k //N },{k,10000,100000,10000}]//TableForm为了解()x π的性质，18世纪的许多数学家都试图找出一个简单的解析函数近似地表示()x π. 1792年年仅15岁的高斯，曾猜测函数()x π的渐进表达式为()x π～ln x x. (3) 后来，他又给出更精确的对数积分式()x π～2()ln x dt Li x t=⎰. (4) 法国数学家勒让德又猜测对于较大的x ， ()x π～ln x x B+， (5) 其中 1.03866B =-(称为勒让德常数).关于()x π的这些猜测的理论证明，最终由两位年龄相仿的法国数学家阿达玛(J.Hadamard 1865—1963)和比利时数学家普桑(A.Poussin 1866—1962)几乎同时(于1896年)各自独立完成，由于该问题在数论中非常重要，数学家们称为“素数定理”.作为实验习题，建议读者用图示方法表达以上三个渐进表达式的逼近情况.4 有关素数的其他问题素数理论中不乏富于挑战性的、至今仍未解决的难题，我们仅列出几个供读者参考.(1) 完全数早在公元前希腊数学家发现数6有一个特性，它等于它自身因子的和：6=1+2+3. 又如28=1+2+4+7+ 14，这种数称为完全数.除了6，28之外，下一个完全数是496，如何找出其他的完全数呢？欧几里得证明了：若21n -是一个素 数，则数12(21)n n --是完全数. 这种形式的完全数显然都是偶数，称为偶完全数. 后来欧拉又证明了凡偶完全数必呈12(21)n n --的形式，其中21n -是一个素数. 据此就在完全数与梅森素数之间建立起联系，显然偶完全数也是非常稀少的.偶完全数还有一些奇特的性质：(ⅰ) 任何偶完全数的个位数字必为6或8；(ⅱ) 除6以外的偶完全数都可以表为几个连续奇数的奇数次方之和，例如332813=+，33334961357=+++；(ⅲ) 偶完全数n 的所有因子(含1和n 本身)的倒数和恒等于2.除了偶完全数外，有没有奇完全数呢？至今没有答案.(2) 孪生素数孪生素数指的是差为2的素数对，例如{3，5}、{5，7}、{11，13}、{41，43}等. 1000多年前，有人就提出了孪生素数猜想：“存在无穷多对孪生素数p和2p ”，可是至今也没人能证明它.问题8 编写程序查找一个指定范围内的孪生素数.解将指定范围内的素数顺次两两交叉分成素数组，如10以内的素数2，3，5，7，分成{2，3}，{3，5}，{5，7}，从每对素数中找出两数之差等于-2的数对，就是问题的解.输入twinp[m_,n_]:=Module[{a,b},a=Select[Prime[Range[PrimePi[m],PrimePi[n]]],PrimeQ];b=Select[Partition[a,2,1],Subtract@@#==-2&]]twinp[300,500]执行后得到{{311,313},{347,349},{419,421},{431,433},{461,463}}(3) 哥德巴赫猜想1742年6月7日哥德巴赫写信给数学家欧拉，提出了以下猜想:任何一个不小于6的偶数，都可以表示成两个奇素数之和.例如：6=3+3，8=3+5，14=3+11，…. 这就是所谓的哥德巴赫猜想(Goldbach's Conjecture)，哥德巴赫猜想是小学生都能明白的问题，其理论证明却也是数论中乃至数学中最棘手的难题，即所谓的“1+1”问题. 为证明猜想的正确性，200多年来，数学家们付出了无数的心血和劳动，至今仍未解决. 我国数学家陈景润，1966年证明了“1+2”，即“任何一个大偶数，都可以表为一个素数及二个素数的乘积之和”，是目前最接近于哥德巴赫猜想的一个结论.问题9 编写程序验证哥德巴赫猜想的正确性.解给出如下程序，程序要求输入的自变量必须取偶数值：Goldbach[n_]:=Module[{a,t={}},If[Mod[n,2]!=0,Print["n must be an even number!"]];a=Select[Prime[Range[PrimePi[n]]],PrimeQ];For[i=1,i<=Length[a]-1,i++,For[j=i+1,j<=Length[a],j++,If[a[[i]]+a[[j]]==n,t=Append[t,{a[[i]],a[[j]]}],Continue[]]]];t]对24—30的4个偶数，应用该程序验证，为输出结果便于观察，使用GridBox[]函数：p=Table[2k,{k,12,15}];GridBox[Table[{p[[i]],Goldbach[p[[i]]]},{i,Length[p]}],RowLines->True,ColumnLines->True]//DisplayForm执行后得到(仅列出部分结果)//DisplayForm245,19,7,17,11,263,23,7,19285,23,11,17307,23,11,19,13,素数理论虽然是纯数学研究的课题，看起来似乎没有什么实用价值，但是，近年来在大批量信息处理以及密码学等问题中显示其广阔的应用前景，感兴趣的可以参考有关著作.三、实验习题习题1 设12{,,,}n p p p 为一个任意的素数集合，通常称形如121n n N p p p =+的数为欧几里得数. 改写问题5.3.3的程序，试之能适应于任意n 个素数的情况，验证欧几里得定理.习题2 做实验回答下列问题：设n 为整数，问(1) 是否存在无穷多形如21n +的素数？(2) 是否存在无穷多形如22n +的素数？(3) 若1n >，在n 与2n 之间是否一定能找到一个素数？习题3 据说用欧拉给出的多项式241n n ++能够产生的素数是最多的，你是否能够找到类似的多项式，由它产生尽可能多的素数. 并对同一个n ，统计它们产生的素数的个数.习题4 从斐波那契数列中找出是素数的项，推测斐波那契数列中素数项是否有无穷多?习题5 证明欧几里得定理：若21n -是一个素数，则数12(21)n n --是完全数.习题 6 除了{3，5}这对孪生素数外，所有的孪生素数都是61n ±的形式，从该种形式的数对中选出孪生素数.提示：用Select[]函数.习题7 根据完全数的定义，编写一个搜索完全数的程序.参考程序：为了缩短搜索时间，仅在个位数为6和8的整数中寻找. 对于过大的n 用这个程序寻找完全数仍然需要较长的时间，例如取n 为10000，找到第四个完全数时大约需要2分多钟. PrefectNum[n_Integer]:=Module[{l ,s0,t ,s ={},S ={}},For[i =1,i <=n ,i ++,If[Last[IntegerDigits[i]]==6||Last[IntegerDigits[i]]==8,s =Append[s ,i ]]]; (*产生不超过n 的所有个位数为6和8的整数*)l =Length[s ];For[m =1,m <l ,m ++,s0=0;t ={};For[i =1,i <s[[m]],i ++,If[Mod[s[[m]],i ]==0,s0+=i ;t =Append[t ,i ]]];If[s0==s[[m ]],S =Append[S ,{s[[m]],t }]]];Return[S ]]习题8 根据完全数与梅森素数的关系，确定一个梅森素数p M ，验证12p p M -确实是一个偶完全数.。
计算pi
实验四你会用几种方法计算PI（圆周率）的值一、问题分析若想计算π的值，就要将跟π有关联的联系在一起，找到与π近似等价的式子，利用计算其值来得到π的值，还有对于含有π的面积、体积等关系式，可以尽量使用较规则的图形来代替进行面积、体积的求解。
二、模型建立2.1数值积分法找一个积分值等于π的定积分,则只要利用定积分计算出的值，就可以得到π的近似值。
2.2幂级数法利用arctanx的泰勒展开式，计算π的近似值。
当x=1时，arctan1=2.3迭代法1976年的迭代算法：2.4 随机模拟法（蒙特卡罗方法）用随机模拟求单位圆面积向单位正方形随机投n块小石头，n很大时小石头大致均匀第分布在正方形中，如果有k块落在单位圆内，单位圆面积的近似值三、解决问题所需的基本理论和方法(1)对于定积分,则只要计算出的值，就可以得到π的近似值,也就是计算出与直线y=0,x=0,x=1所围成的曲边梯形，而对于此类计算往往采用数值积分梯形公式计算。
梯形公式：将积分区间n等分将所有梯形面积加起来得到Trapz(x)：输出数组x，输出按梯形公式x的积分(单位步长)Trapz(x,y):计算y对x的梯形积分，其中x、y定义函数关系y=f(x)(2)利用arctanx的泰勒展开式，计算π的近似值。
函数taylor用于实现Taylor级数r=taylor(f,n,v),指定自变量v和阶数nr= taylor(f,n,v,a),指定自变量v、阶数n，计算f在a的级数(3)级数法由于利用arctanx的幂级数展开法的收敛较慢，可采用公式的计算来求pi值。
(4)特殊公式（BBP）四、设计算法、编程求解4.1数值积分法梯形公式Matlab代码：format longx=0:0.1:1; % x=0:0.01:1; x=0:0.02:1; x=0:0.001:1; x=0:0.0001:1;y=sqrt(1-x.^2);pi=4*trapz(x,y)4.2幂级数法Matlab代码：(1)format longsyms xf=atan(x);t= taylor(f,10,x,0); % t= taylor(f,100,x,0); t= taylor(f,500,x,0);t= taylor(f,1000,x,0); t= taylor(f,10000,x,0); x=1;pi=4*eval(t)(2)format longsyms xf=atan(x);t= taylor(f,10,x,0);x=1/5;s1=eval(t);x=1/239;s2=eval(t);pi=16*s1-4*s2当n=10时，pi =3.141592682404399format longa=1;b=1/sqrt(2);s=1/2;for n=1:1:10n,y=a;a=(a+b)/2;b=sqrt(b*y);c=a^2-b^2;s=s-2^n*c;pi=2*a^2/send4.4蒙特卡罗方法Matlab代码：format longs=0;n=10; % n=100; n=1000; n=10000; n=100000; n=1000000 for i=1:na=rand(1,2);if a(1)^2+a(2)^2<=1s=s+1;endendpi=4*s/n4.5 BBP公式format longsyms xy=1/16^x*(4/(8*x+1)-2/(8*x+4)-1/(8*x+5)-1/(8*x+6));s=0;for x=0:1:10s1=eval(y);x,s=s+s1end五、分析求解结果由上表可知，蒙特卡罗方法计算出的pi值与真实值的误差相差较大并且收敛速度很慢；对于级数法，但由于所选择的的级数方法、公式不同，得到的结果也就不同，收敛速度较慢，而的收敛速度就较快；数值积分法和迭代法准确度较高，但数值积分法的收敛速度没有迭代法快、精度高，所以一般情况下采用迭代法求近似值较准确。
pi的计算 实验报告
Ramanujan公式1914年，印度数学家Srinivasa Ramanujan在他的论文里发表了一系列共14条圆周率的计算公式，这是其中之一。
这个公式每计算一项可以得到8位的十进制精度。
1985年Gosper用这个公式计算到了圆周率的17,500,000位。
3、AGM(Arithmetic-Geometric Mean)算法Gauss-Legendre公式：初值：重复计算：最后计算：这个公式每迭代一次将得到双倍的十进制精度，比如要计算100万位，迭代20次就够了。
1999年9月Takahashi和Kanada用这个算法计算到了圆周率的206,158,430,000位，创出新的世界纪录。
4、Borwein四次迭代式：初值：重复计算：最后计算：这个公式由Jonathan Borwein 和Peter Borwein 于1985年发表，它四次收敛于圆周率。
5、Bailey-Borwein-Plouffe 算法014211()1681848586n n n n n n π∞==---++++∑这个公式简称BBP 公式，由David Bailey, Peter Borwein 和Simon Plouffe 于1995年共同发表。
它打破了传统的圆周率的算法，可以计算圆周率的任意第n 位，而不用计算前面的n-1位。
这为圆周率的分布式计算提供了可行性。
1997年，Fabrice Bellard 找到了一个比BBP 快40％的公式：第三部分：对于π的几种计算的研究和讨论： 1、数值积分法（I ）利用积分公式⎰-=10214dx x π计算πn=10 ans =; n=20 ans =; n=50 ans =; n=100 ans =; n=200 ans =; n=500 ans =; n=1000 ans =; n=2000 ans =;半径为1的圆称为单位圆，它的面积等于π。
只要计算出单位圆的面积，就算出了π。
多种方法计算Pi并且精确度比较
多种⽅法计算Pi并且精确度⽐较多种⽅法计算圆周率并⽐较精确度【摘要】本⽂介绍了多种⽅法求圆周率的近似值并对各种⽅法进⾏精确度的⽐较得出具体情况选择的⽅法，且通过mathematica 编程模拟实验过程，得出各种⽅法的特点。
【关键字】圆周率数值积分法泰勒级数法蒙特卡罗法拉马努⾦公式法0.引⾔平⾯上圆的周长与直径之⽐是⼀个常数，称为圆周率，记作π。
在很长的⼀段时期，计算π的值是数学上的⼀件重要的事情。
有数学家甚⾄说：“历史上⼀个国家所得的圆周率的准确程度，可以作为衡量⼀个国家当时数学发展的⼀⾯旗帜。
”⾜以见圆周率扮演的是⾓⾊是如此举⾜轻重。
π作为经常使⽤的数学常数，它的计算已经持续了2500多年了，到今天都依然在进⾏着，中间涌现出许多的计算⽅法，它们都各有千秋，在此，我们选择⼏种较典型的⽅法，包括数值积分法，泰勒级数法，蒙特卡罗法，韦达公式法，拉马努⾦公式法以及迭代法来和⼤家⼀起体验π的计算历程，同时利⽤mathematica 通过对各种⽅法作精确度的⽐较得出选择的优先顺序，为相关的理论研究提供⼀定的参考价值。
1.数值积分法以单位圆的圆⼼为原点建⽴直⾓坐标系，则单位圆在第⼀象限内的部分G 是⼀个扇形，由曲线y=211x+(x ∈[0,1])及两条坐标轴围成，它的⾯积S=4π。
算出了S 的近似值，它的4倍就是π的近似值。
（41112π=+?dx x ）⽤⼀组平⾏于y 轴的直线x=i x (1≤i ≤n-1，a=0x <1x <2x <...地看作抛物线，就得到⾟普森公式。
梯形公式：S ≈]2...[10121y y y y y n ab n +++++-- ⾟普森公式：S ≈)]...(4)...(2)[(62123211210--+++++++++-n n n y y y y y y y y n abMathematica 程序如下：n=1000;y[x_]:=4/(1+x^2);s1=(Sum[y[k/n],{k,1,n-1}]+(y[0]+y[1])/2)/n;s2=(y[0]+y[1]+2*Sum[y[k/n],{k,1,n-1}]+4*Sum[y[(k-1/2)/n],{k,1,n}])/(6*n); Print[{N[s1,20],N[s2,20],N[Pi,30]}]注：以上s1,s2分别是⽤梯形共识和⾟普森公式计算出的π。
圆周率的数学实验算法设计
3. 级数方法计算圆周率
算法原理：
由于 1 1+x2
=1- x2+x4- …+(- 1)n- x1 2n-2+…，
两边关于[0,x]上积分可得：
arctanx=x- x3 + x5 - …+(- 1)n-1 x2n-1 …
35
2n- 1
令 x=1 即：π =1- 1 + 1 - …+(- 1)n-1 1 …从而：
与直径之比。π 的计算伴随着人类的进步而发展，许多数学家在其计算
上花费了巨大精力。在中国有刘徽、祖冲之等，在国外有阿基米德、卡西
等。近现代科学家如华罗庚、闵嗣鹤、严士健等在其数论论文中也对圆
周率问题进行了探讨。有些数学家甚至说：“历史上一个国家所算得的
圆周率的准确程度，可以作为衡量一个国家当时数学发展的一面旗
基金项目：本文受华东交通大学校研基金资助。 作者简介：王广超（1978- ），男，汉族，山东微山人，硕士，讲师，华东交通大学基础科学学院。
—8—
帜。”
在信息技术发展日新月异的今天，数学教育面临着巨大的挑战。信
息时代以计算机作为工具，许多无法求解的问题有可能得到解决 ，计算
机技术应用领域不断涌现出新的要求 , 使得大量新兴的数学正在被有
效地采用。面对这样的形势，数学实验是将数学知识、数学建模思想与
计算机应用三者紧密结合的一体的教学模式。本文设计了 π 的多种计
for i=1:n
if sin(a(k))*sin((a(k)+b(k))/2)<0
a(k+1)=a(k);
b(k+1)=(a(k)+b(k))/2;
k=k+1;
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2k 1
2n 1
计算一下要精确到Pi的200位小数需要取 级数的多少项？
2015/10/27
上一页
下一页
主 页
He Renbin
利用级数计算Pi
2、欧拉的两个级数（1748年发现） 2 1 2 6 k 1 k 2 1 8 k 0 (2k 1) 2 这两个级数收敛也非常缓慢，计算时实 用价值不大。
2015/10/27
上一页
下一页
主 页
He Renbin
“割圆术”中学问多
我国2000多年前的《周髀算经》称“周三径 一”，这是π的第一个近似值，叫做“古率”。 据说，汉代大科学家、文学家张衡，有“圆 周率一十之面”的推算。清代李潢考证这句话意 思为π≈sqrt(10)。 魏晋间刘徽由圆内接正六边形依次倍增到正 192边形，计算周长与值径之比，得 3.141024< π<3.142704 实际应用时取3.14，或分数值157/50。
将积积分区 n 等分，即 x i i/n, i 0,1,, n 。 将所有这所有这些梯形起来就得到 1 n 1 f(x i ) f(x i 1 ) S n i 0 2
2015/10/27
上一页
下一页
主 页
He Renbin
数值积分法计算Pi
2、辛普森（Simpson）公式
下一页
主 页
He Renbin
拉马努金（Ramanujan)公式
另一个经过改进的计算公式为：
1 12 640320
3 2
(1) (6n)! 13591409 545140134n 3 3n ( n ! ) ( 3 n )! 640320 n 0
n
级数每增加一项，可提高14位小数的 精度。
2015/10/27
2
上一页
下一页
主 页
He Renbin
“割之弥细，失之弥少，割之又 割，则与圆合体而无所失矣。”
面积与边长有如下关系：
S6( n1) 2 4 a6 n
2
圆面积S与多边形的面积Sn之间有如下关系：
S 2 n S 2S 2 n S n
2015/10/27
上一页
He Renbin
韦达（VieTa）公式
1、从sint开始
t t sin t 2 cos( ) sin( ) 2 2 t t t 4 cos( ) cos( ) sin( ) 2 4 4 t t t t 8 cos( ) cos( ) cos( ) sin( ) 2 4 8 8
2015/10/27
另一种用蒙特卡罗法来计算Pi的方法是 1777年法国数学家蒲丰（Buffon）提出的随 机掷针实验。其步骤如下： （1）取一张白纸，在上面画出许多间 距为d的等距平行线。 （2）取一根长度为l(l<d)的均匀直针， 随机地向画有平行线的纸上掷去，一共掷n 次。观察针和直线相交的次数m。
2015/10/27
2015/10/27
上一页
下一页
主 页
He Renbin
利用级数计算Pi
加速效果非常明显！
2015/10/27
上一页
下一页
主 页
He Renbin
蒙特卡罗（Monte Carlo）法
单位圆的面积等于Pi，使用蒙特卡罗法， 即用随机投点的方法来求出这个面积Pi的近 似值。具体方法如下： 在平面直角坐标系中，以O(0,0)， A(1,0)，C(1,1)，B(0,1)为四个顶点作一个正 方形，其面积S＝1。以原点O为圆心的单位 圆在这个正方形内的部分是圆心角为直角 的扇形，面积为S1＝Pi/4。
1、莱布尼茨级数（1674年发现）
(1) 4 k 0 2k 1
k
2015/10/27
上一页
下一页
主 页
He Renbin
利用级数计算Pi
1844年，数学家达什在没有计算机 的情况下利用此式算出了Pi的前200位小 数。使用误差估计式 n (1) k 1
r ( n) 4
2015/10/27
上一页
下一页
主 页
He Renbin
韦达（VieTa）公式
1593年，韦达首次给出了计算Pi的 精确表达式：
2 2 2 2 2 2 2 2 2 韦达公式看起来有些神秘，其实它 的导出过程所用的都是朴实简洁的数学 方法。 2
2015/10/27
上一页
下一页
主 页
2nl md
特别取针的长度L为d/2时，π=n/m。
2015/10/27
上一页
下一页
主 页
He Renbin
拉马努金（Ramanujan)公式
目前，计算pi的一个极其有效的公式为
2 2 (4n)! 1103 26390n 4 4n 9801 n 0 (n!) 396 1
这个级数收敛得非常快，级数每增加 一项，可提高大约8位小数的精度。
下一页
主 页
He Renbin
刘徽不等式
借助于计算机来完成刘徽的工作： a(1)=sqrt(3);b(1)=3*sqrt(3)/2; for i=2:6 a(i)=sqrt(2-sqrt(4-a(i-1)^2)); b(i)=3*2^(i-2)*a(i); c(i)=2*b(i)-b(i-1); end n=[3,6,12,24,48,96]; size(b) result=[n;a;b]
2015/10/27
上一页
下一页
主 页
He Renbin
蒙特卡罗（Monte Carlo）法
在这个正方形内随机地投入n个点，设 其中有m个点落在单位扇形内。则
2015/10/27
m S1 4m , n S 4 n
随机投点如何来实现？
上一页 下一页 主 页
He Renbin
蒲丰（Buffon）掷针实验
2015/10/27
上一页
下一页
主 页
He Renbin
利用级数计算Pi
观察级数可知，x的值越接近于0，级 数收敛越快。由此可以考虑令
1 1 x tan α , arctan 5 5 2 tan 2x 5 tan 2 2 2 1 tan 1 x 12 2 tan 2 120 tan 4 1 2 1 tan 2 119
上一页
下一页
主 页
He Renbin
韦达（VieTa）公式
所以，对任意N，总有
sin t 2 N t N t sin( N ) cos( n ) t t 2 n 1 2 sin t t 令N , 有 ＝ cos( n ） n 1 t 2 2 取t , 得到 ＝ cos( n 1） 2 n 1 2
令yi f ( xi ) 1 Si ( yi 1 4 y 1 yi ) i n 2
n 1 n 1 S [( y0 yn ) 2 yi 4 y 1 ] i 6n i 1 i 1 2
2015/10/27
上一页
下一页
主 页
He Renbin
利用级数计算Pi
下一页
0 3.4019 3.2117 3.1594 3.1461 3.1427
He Renbin
主 页
割圆术的意义
刘徽创立的割圆术，其意义不仅在 于计算出了Pi的近似值，而且还在于提 供了一种研究数学的方法。这种方法相 当于今天的“求积分”，后者经16世纪 英国的牛顿和德国的莱布尼茨作系统总 结而得名。鉴于刘徽的巨大贡献，所以 不少书上把他称做“中国数学史上的牛 顿”，并把他所创造的割圆术称为“徽 术”。
思考： 如何利用韦达公式构造 出一种迭代算法？
2015/10/27
上一页
下一页
主 页
He Renbin
数值积分法计算Pi
定积分
4 dx 2 1 x 0
计算出这个积分的数值，也就得到了Pi 的值。
1
2015/10/27
上一页
下一页
主 页
He Renbin
数值积分法计算Pi
1、梯形公式
2015/10/27
上一页
下一页
主 页
He Renbin
实验指导
世界上数学家们一致公认： “历史上一个国家计算圆周率的准 确度，可以作为衡量这个国家当时 数学水平的一个标志。”
2015/10/27
上一页
下一页
主 页
He Renbin
π值——算法美的追求
π作为圆周率的符号，是由著名数学家欧勒于 公元1737年首先使用的。古代的希伯来人，在描 述所罗门庙宇中的“熔池”时曾经这样写道： “池为圆形，对径为十腕尺，池高为五腕尺，其 周长为三十腕尺。”可见，古希伯来人认为圆周 率等于3。不过，那时的建筑师们，似乎没有人不 明白，圆周长与直径的比要比3大一些。 公元前3世纪古希腊大数学家阿基米德求出了 223/71<π <22/7。
2015/10/27
上一页
下一页
主 页
He Renbin
韦达（VieTa）公式
2、从cos(pi/4)开始
2 cos( ) 4 2 cos(
8
)
cos
4 2
1
2 1 2 2
2 2 2
2015/10/27
上一页
下一页
主 页
He Renbin
韦达（VieTa）公式
3、使用VieTa公式计算Pi的近似值
上一页
下一页
主 页
He Renbin
蒲丰（Buffon）掷针实验