常数e和γ的Monte Carlo模拟估计.pdf

上海中学数学2023年第7 一8 期3常数e和的Monte Carlo模拟估计200234上海师范大学数理学院摘要：探讨了采用Monte Carlo(蒙特卡罗)方法估计常数e和的问题.通过概率论中的匹配问题对常数e进行估算，并应用概率论中的几何概型（面积法）使用软件编程估计常数e和的近似值，还就常数e和的估算提出了试验方案.关键词：自然对数底；欧拉常数；蒙特卡罗模拟完王蓉华0 引 言历史上著名的“蒲丰投针”问题具有深远的意义，它采用随机模拟的方法得到了圆周率元的近似值，最早体现统计学中应用MonteCarlo（蒙特卡罗)模拟的思想.“蒲丰投针”问题已经成为中学教学中一个激发学生好奇心的经典素材.例如，北师大版初中数学教材第七章“概率”的阅读材料就是“蒲丰投针”问题，教材将其作为频数与频率的补充阅读.值得指出的是，在人教社2 0 11年版义务教育教科书数学（七九年级)中，九年级上册第二十五章第14 9 页“试验与探究”栏目的内容就是“元的估计”，它也是利用Monte Carlo模拟的思想方法，设计一个类似于“蒲丰投针”的试验方案来估计圆周率元的值.其试验方案如下：在一个很大的房间里铺上边长为1的正方形地砖，每块地砖上画有内接圆，然后往房间里撒很多大米（个数为N），统计落在内接圆里的米粒数（记为no），以及落在圆上的米粒数（记为ni），则元的估计值允可以通过如下算式得到：no+n2元=4N常数元称为圆周率.其实，在中学数学教学中，学生还会遇到许许多多的常数，比较常用的有常数e 等。.常数e=2.71828=lim(1+1)n称为自然对数底，是无理数.关于e的较为详细的应用可见文献1、文献2 而=0.572.im(n m)称为欧拉常数,现+k=1k仍无法确定其是有理数还是无理数，欧拉常数有多种表达式，具体可见文献 3.下面的问题是如何利用MonteCarlo模拟（随机试验）的方法得到常数e,的近似值.1常数e的Monte Carlo模拟估计常数e的MonteCarlo模拟估计通常有两种方法，一种是匹配法，另一种是面积法.1.1匹配法假设某人撰写了内容不同的n封信，并在n个信封上写上不同的地址，接着随意地把这n封信放人n个信封中，现求至少有一封信刚好放进正确的信封的概率p（n).这就是古典概型中著名的“匹配问题”记事件A，表示 把第ii=1,2，n）封信放进正确的信封”,则易见p(n）=P(U A,).由于是随意把信装进信封的，那么把任意一封信放进正确的信封的概率P(A,)是,即有P(A,）=n=nn1.进而，因为可把第一封信放进n个信封中的任何一个,第二封信可放进剩下的n一1个信封中的任何一个，那么把第一封信和第二封信都放进正确的信封的概率P(A,A2）是-类似地，把任意第in(n-1):封信和第封信（ii）都放进正确的信封的概率P(A;A,)是1n(n-1)11类似地，把任何三封信i，i 和k（i n(n一1)2！n+80jk）都放进正确的信封的概率P(A,A;A,）是(n-1)(-2)因此ZP(AA.)=c:.11n(n-1)(n-2)n封信都正确放进各自信封的概率P（A A 2 A,）是工，从而可得至少把一封信放进正确的信封的概n！2！+3！4!1则ZP(A,A,)=Ciiiik1持续这个过程，计算把所有3！4n十的时候,p(n）的值趋向于下式的极限1+1-1+.=1-e-,即当lim p(n)=12！+3!4!n很大时，至少有一封信被放进正确的信封的概率p(n）的值接近1e-l=0.63212.注意到p(n)的值会随着 n的增加形成一个波动序列.当n以偶数2，4，6，增加时，p（n）的值会逐渐增大到极限值0.6 32 12；当n以奇数3，5，7，增加时，p（n)的值会逐渐减小到同样的极限值.p(n)的收敛速度很快，当n=7时，p(7)的值已经可以精确到小数点后第4 位了。利用上述的“匹配问题”的结果，可以通过实际试验得到e的近似值.记事件A表示“至少有一封信刚好放进正确的信封”,由于P(A)1一e-1,所以只要得到P(A)的近似值，就可以得到e的近似值（或称为估计值）为 1-P(A)崔宜兰在文献 4 中曾经操作过该试验.崔老师与某班级的4 0 名学生一起利用扑克进行匹配试验，并对e进行了估计.崔老师的方法是，取扑克中两种花色共2 6 张牌，每次随机取两张，若成对则认为是一个匹配.试验时，先将扑克牌充分洗匀，若出现对子就停止试验，洗匀后再进行下一轮试验，否则摸完26张牌.共进行了2 5 0 0 次试验，其中有对子出现的有15 7 8 次，则得=2.7 114 9 6 7 4 6，而e2.718281828.,其误差为le-el6.810-31.2面积法张乐成等在文献 5 中采用的是定积分=ln2,其算法是向矩形120y1中投点N次（N很大），统计落人阴影部分的点数n（如图1所示)，则e的估计值为=2岁.1X11图1dx=ln2的示意图N107e的近似值2.71827960864361上海中学数学2023年第7 一8 期1所示.毕艳辉和程希明在文献 6 中采用的是定积分1d=ln2，其算法是向矩形0 1,0yJoa+11中投点N次（N很大），统计落入阴影部分的点数n如图2 所示），则e的估计值为e=2%.0.51图2文献 5 与文献6 实质是一回事，文献 6 通过向图2 中的正方形投点计算得到e的近似值，如表2所示.值得指出的是，文献5 是通过VB6.0实现的，而文献 6 则是通过投点实现的.通过计算机软件来判断所投的点是否属于阴影部分是没有问题的，但如果要学生做投点试验，文献 6 在操作层面将会遇到很多问题.它提出在半米高的上方投点，如dc果不是半米，而是更高，试验中可能会出现许许多多Jia的点不落入正方形内的情形，所以投点的方向、所用的力量还有高度等都会影响试验结果，还将拉长试验时间.那么，在操作层面该如何处理这一问题呢？回想“蒲丰投针”问题，我们注意到“蒲丰为何要画一系列间隔相等的平行线而不是若干条平行线”，这正是“蒲丰投针”试验设计的巧妙之处.其试验方案设计如下：首先寻找N颗直径很小的“滚珠”，然后在一张很大的白纸上画满边长为1的正方形，每个正方形画有类似于图2 的图形，接着将这N颗“滚珠”撒在这张白纸上，统计落在阴影部分内部的“滚珠”数（记为no），落在阴影部分周边的四条线上的“滚珠”数（记为ni），则e的估计值可2X表1文献 5 的试验结果1082.71827943845056文献 5 采用VB6.0实现了算法，其结果如表dx=ln2的示意图+1以通过如下算式得到：=2%，1092.71829011738232x+IX110102.718282705345上海中学数学2023年第7 一8 期表2 文献 6 的试验结果N100n66e的近似值2.9048462欧拉常数的MonteCarlo模拟估计注意到有如下定积分1记函数=1n y1的图形如图3所示.1Inx1-x图3dx=的示意图如果通过软件很容易得到欧拉常数的估计值，即向图3上投点N次（N很大），统计落入阴影部分的点数n,则的估计值为=.下面采用R语言编写,其代码如下（以N=104为例）：n-le4#set the value of nm O#initialize the value of m#Perform the simulationfor(i in l:n)(X-runif(1)Y-runif(1)if(Y=1/(log(X,base=exp(1)+1/(1-X)(m-m+1N的近似值N的近似值550010003356822.8138342.7630771#Calculate Sdx=.S-m/n,0 x1,其在0 1,0#Print the resultscat(m=,m,“n)#m=5755cat(S=,S,“n)#S=0.5755分别取N=10,i=2,3，,9,得到的近似值如表3所示.如果不采用软件编程投点，为求欧拉常数的模拟估计值，试验方案可以设计如下：首先寻找N颗直径很小的“滚珠”，然后在一张很大的白纸上画满边长为1的正方形，每个正方形画有类似于图3的图形，接着将这N颗“滚珠”撒在这张白纸上，统计落在阴影部分内部的“滚珠”数（记为no），落在阴X影部分周边的四条线上的“滚珠”数（记为ni），则的估计值可以通过如下算式得到：=，N一2参考文献1张立新.从物理学中数学常数e思考数学之自然属性J.首都师范大学学报，2 0 15（3）：2 5-2 9.2徐怀，郭志军.常数e的几点说明.阴山学刊，2018(3):135-127.3谭卫平.欧拉常数的几种表示 门.数学理论与应用，2005(4):116-117.4崔宜兰.无理数e的统计估计法 J.工科数学，1998(3):116-119.5张乐成，邵梅，迟津愉，等.一种基于Monte Carlo方法计算数学常数e值的算法 J.计算机与现代化，2012(5):14-15;19.6毕艳辉，程希明.用随机试验法计算欧拉常数eJ.统计与管理,2 0 16(1):9-10.表3欧拉常数的模拟结果1021030.580.5661061070.5778650.5769956500035012.699449non11041050.57550.57906108100.57714710.5772272