数学实验之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）根据所得结果猜测椭圆的面积是。

素数探秘一、实验目的素数(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（圆周率）的值一、问题分析若想计算π的值，就要将跟π有关联的联系在一起，找到与π近似等价的式子，利用计算其值来得到π的值，还有对于含有π的面积、体积等关系式，可以尽量使用较规则的图形来代替进行面积、体积的求解。

二、模型建立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值与真实值的误差相差较大并且收敛速度很慢；对于级数法，但由于所选择的的级数方法、公式不同，得到的结果也就不同，收敛速度较慢，而的收敛速度就较快；数值积分法和迭代法准确度较高，但数值积分法的收敛速度没有迭代法快、精度高，所以一般情况下采用迭代法求近似值较准确。

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并且精确度⽐较多种⽅法计算圆周率并⽐较精确度【摘要】本⽂介绍了多种⽅法求圆周率的近似值并对各种⽅法进⾏精确度的⽐较得出具体情况选择的⽅法，且通过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分别是⽤梯形共识和⾟普森公式计算出的π。

实验09数值微积分与方程数值解（第6章）《数学软件》课内实验王平（第6章MATLAB数值计算）一、实验目的1.掌握求数值导数和数值积分的方法。

2.掌握代数方程数值求解的方法。

3.掌握常微分方程数值求解的方法。

二、实验内容1.求函数在指定点的数值导数某f(某)1程序及运行结果：某2某36某2某3某2,某1,2,3 022.用数值方法求定积分(1)I120cot24in(2t)21dt的近似值。

程序及运行结果：(2)I220ln(1某)d某1某2程序及运行结果：3.分别用3种不同的数值方法解线性方程组6某5y2z5u49某y4zu133某4y2z2u13某9y2u11程序及运行结果：4.求非齐次线性方程组的通解2某17某23某3某463某15某22某32某449某4某某7某22341程序及运行结果（提示：要用教材中的函数程序line_olution）：5.求代数方程的数值解(1)3某+in某-e某=0在某0=1.5附近的根。

程序及运行结果（提示：要用教材中的函数程序line_olution）：(2)在给定的初值某0=1，y0=1，z0=1下，求方程组的数值解。

in某y2lnz70y33某2z10某yz50程序及运行结果：26.求函数在指定区间的极值某3co某某log某(1)f(某)在(0,1)内的最小值。

e某程序及运行结果：332(2)f(某1,某2)2某1在[0,0]附近的最小值点和最小值。

4某1某210某1某2某2程序及运行结果：7.求微分方程的数值解，并绘制解的曲线某d2ydy5y0d某2d某y(0)0y'(0)0程序及运行结果（注意：参数中不能取0，用足够小的正数代替）：令y2=y,y1=y'，将二阶方程转化为一阶方程组：1'5yy1某1某y2'y2y1y(0)0,y(0)0218.求微分方程组的数值解，并绘制解的曲线y'1y2y3y'yy213y'0.51yy123y1(0)0,y2(0)1,y3(0)1程序及运行结果：3三、实验提示四、教程：第6章MATLAB数值计算（2/2）6.2数值微积分p1556.2.1数值微分1.数值差分与差商对任意函数f(某)，假设h>0。

圆周率π圆周率是一个常数（约等于3.1415926），是代表圆周长和直径的比例。

它是一个无理数，即是一个无限不循环小数。

但在日常生活中，通常都用3. 14来代表圆周率去进行计算，即使是工程师或物理学家要进行较精密的计算，也只取值至小数点后约20位。

π(pai)是第十六个希腊字母，本来它是和圆周率没有关系的，但大数学家欧拉在一七三六年开始，在书信和论文中都用π来代表圆周率。

既然他是大数学家，所以人们也有样学样地用π来表示圆周率了。

但π除了表示圆周率外，也可以用来表示其他事物，在统计学中也能看到它的出现。

π=Pai（π=Pi）古希腊欧几里德《几何原本》（约公元前3世纪初）中提到圆周率是常数，中国古算书《周髀算经》（约公元前2世纪）中有“径一而周三”的记载，也认为圆周率是常数。

历史上曾采用过圆周率的多种近似值，早期大都是通过实验而得到的结果，如古埃及纸草书（约公元前1700）中取pi=（4/3）^4≒3.1604 。

第一个用科学方法寻求圆周率数值的人是阿基米德，他在《圆的度量》（公元前3世纪）中用圆内接和外切正多边形的周长确定圆周长的上下界，从正六边形开始，逐次加倍计算到正96边形，得到（3+（10/71））<π<（3+（1/7）），开创了圆周率计算的几何方法（亦称古典方法，或阿基米德方法），得出精确到小数点后两位的π值。

中国数学家刘徽在注释《九章算术》（263年）时只用圆内接正多边形就求得π的近似值，也得出精确到两位小数的π值，他的方法被后人称为割圆术。

他用割圆术一直算到圆内接正192边形，得出π≈根号10 （约为3.16）。

南北朝时代著名数学家祖冲之进一步得出精确到小数点后7位的π值（约5世纪下半叶），给出不足近似值3.1415926和过剩近似值3.141592 7，还得到两个近似分数值，密率355/113和约率22／7。

他的辉煌成就比欧洲至少早了1000年。

其中的密率在西方直到1573才由德国人奥托得到，1625年发表于荷兰工程师安托尼斯的著作中，欧洲不知道是祖冲之先知道密率的，将密率错误的称之为安托尼斯率。

蒲丰掷针法一． 实验目的通过蒲丰掷针法来计算pi 的近似值，其本质是蒙特卡罗法随机模拟，通过求概率来求得pi 的近似值。

二． 实验内容与要求（根据问题重新叙述）在白纸上画许多等距为d 的平行线，将一根长为的d/2的直针随机投掷向白纸，在进行n 次实验之后其中有m 次与平行线相交，当n 很大时，pi 的近似值可认为是n/m 。

三． 实验原理（问题假设，分析，模型建立）由于针投到纸上的时候，有各种不同的方向和位置，但是，每一次投针时，其位置和方向都可以由两个量唯一确定，即针的中点距平行线的距离（最近的平行线）d 和偏离水平的角度α。

因此只要随机生成n 对这样的 d 和 α，就可以模拟n 次的投针实验。

设平行线之间的距离为d ，以y 表示针的中点到最近的一条平行线的距离，α表示针与平行线间的夹角，其中0≤y ≤2d ,0≤α≤Pi 。

而直针与一平行线相交的充要条件：0≤y ≤4d sin α ，0≤α≤Pi 。

四． 实验过程（模型求解，模型结果）当n 的值为1000时，重复5次后，计算结果分别为：PI= 3.1546 PI= 3.2895 PI= 3.3223 PI= 3.2573 PI= 3.0303当n 的值为10000时，重复5次后，计算结果分别为：PI= 3.1546 PI= 3.1676 PI= 3.1066 PI= 3.2092 PI= 3.2010 当n 的值为50000时，重复5次后，计算结果分别为：PI= 3.1584 PI= 3.1283 PI= 3.1131 PI= 3.1273 PI= 3.1644五． 实验总结（结果分析）当n 的值越来越大时，可以发现通过随机模拟得到的PI 的值越来越接近真实的pi 的值（pi 的值的前8位为pi=3.1415926）虽然不是很接近pi 的值，但得到的PI 的值离真实的pi 的值越来越近。

且得到的PI 的值也越来越稳定。

六． 附录（源程序）。

估算和乘除法运算，这些数学知识点在解题过程中都要涉及到。

估算和乘除法运算是数学学科中非常重要的知识点，它们在现实生活中也有广泛的应用。

估算可以帮助我们快速地对数据进行大致的评估或预测，而乘除法运算则涉及到实际的计算过程。

在本文中，我们将具体介绍这两个知识点的定义、方法和举例。

一、估算1.定义估算是指在缺乏精确数据的情况下，通过一些近似方法或经验公式来估计某个量或数值。

估算过程通常包括确定一个适当的量级、对数据进行调整和推断，以得到一个大致的结果。

2.方法估算通常有以下几种方法：（1）保留有效数字：在进行估算时，通常只需要保存有效数字，即数字中最高位数和最后一位数。

例如，对于数字145.78，将其保留两位有效数字为146。

（2）舍去或进位：在进行估算时，可以根据需要去掉某些数字或将它们进位，来得到更加精确或近似的结果。

（3）利用比例和经验公式：估算还可以根据比例和经验公式来计算。

比例是指两个物体之间的数量比，例如人口密度、体积比等。

经验公式是科学家或数学家在长期研究中得出的数学公式，可以用来估算某些物理或化学量。

（4）利用近似值：当数据非常大或非常小时，可以利用近似值来进行估算。

这些近似值通常需要在多次实验和计算中确定出来，例如pi的值3.14、自然对数e的值2.71等。

3.举例利用估算的方法，我们可以快速地预测出一些结果。

例如，如果我们要计算100万块砖头铺满一块足球场的面积，我们可以使用估算方法来预估这个面积。

首先，我们可以估计一块砖头的面积大约为0.1平方米，然后将100万乘以0.1，得到一些大约为10万平方米的答案。

然而，这个估算结果可能不精确，因为砖头之间有缝隙，所以我们可以进一步调整，将结果乘以0.8（即80%），最终得到大约为8万平方米的答案。

二、乘除法运算1.乘法运算乘法是指将两个或多个数相乘的运算。

乘法的基本性质包括交换律、结合律和分配律。

例如，对于任意的a、b、c三个数，有以下基本性质：（1）交换律：a×b=b×a（2）结合律：a×(b×c)=(a×b)×c（3）分配律：a×(b+c)=a×b+a×c2.除法运算除法是指将一个数分成若干等份的运算。

实验四非线性方程近似解一、按揭还贷㈠问题描述（1）小张夫妇以按揭方式贷款买了一套价值20万元的房子，首付5万元，每月还款1000元，15年还清。

问贷款利率是多少？（2）某人想贷款50万元购房，他咨询了两家银行，第一家银行开出的条件是每月还4500元，15年还清；第二家银行开出的条件是每年还45000元，20年还清。

从利率方面看，哪家银行较优惠？（简单假设年利率=月利率*12）㈡简要分析初看本题，一个简单的思路是每次测试一个利率值，以这个值为基础计算15年后所剩还款数量，通过结果判断应将利率值增大或减小，从而实现迭代。

这其实是一个二重迭代的过程，之所以这样是因为不容易一眼看出本题的非线性方程。

事实上，转换思路后，可以利用一个简单的方程描述整个迭代过程。

这样就将二重迭代转化为了一层迭代。

使得处理更加简便。

㈢方法与公式1、解题方法(1)二次迭代给定总的本金，从每一次还款中扣去这段时间中增加的利息，再将其还到本金，使本金总量逐渐减少。

代码：for i = 1:time*12less = (repay-left*interest);left = left - less;(2)方程描述虽然并不是所有本金都在还款的整个期间中产生了相应的利息，但是可以设想成这样，与此同时，还款从在相应的还款时间开始产生利息，这样可以得出，两者最终的“本息和”相等，即nA(1+q)n=P(1+q)n−ii=1其中A为总还款金额，q为了利率，P为每次还款金额。

2、解方程方法(1)牛顿法x k+1=x k−f(x k) f′(x k)(2)直接使用公式fzero()㈣结果与分析1、第一问：(1)二次迭代[i,q]=iterate(150000,1000,15,2,0,1,100,10^-6); 公式表意为：总贷款量=200000-50000=150000；每月还款100元；还款期限15年；还款方式为按月还款；迭代区间设定为[0,1]；最大迭代次数为100次；精度要求为10^-6；最终结果为：迭代次数：45；使用时间0.003030989435705s；利率为0.002081163889457。