数学实验报告
实
验
二
学院：数学与统计学院
班级：信息与计算科学(1)班
姓名：郝玉霞
学号：201171020107
实验二
一、实验名称：π的计算
二、实验目的：首先在Mathematica环境中用多种方法计算圆周率π的值，通过
实验来体会各种方法的区别，比较各种方法的优劣，接着尝试自己提出新的
方法来计算圆周率π的值。

三、实验环境：学校机房，Mathematica软件。

四、实验的基本理论和方法
1、用Mathematica绘图函数Plot绘制圆周率π；
2、计算圆周率π的数值积分法、泰勒级数法、蒙特卡罗法，并且利用特定
的公式来计算圆周率π。

五、实验的内容和步骤及实验的结果和结果分析
步骤一、数值积分法计算π
因为单位圆的半径为1，它的面积等于π，所以只要计算出单位圆的面积，就算出了π。

在坐标轴上画出以圆点为圆心，以1为半径的单位圆，则这个单位圆在第一象限的部分是一个扇形，而且面积是单位圆的1/4，于是，我们只要算出此扇形的面积，便可以计算出π。

当n=5000时；
语句：
n=5000;y[x_]:=4/(1+x*x);
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,30],N[Pi,30]}];
实验结果：
3.1415926469231265718,3.14159265358979323846264334
3.14159265358979323846264338328
当n=10000时；
语句：
n=10000;y[x_]:=4/(1+x*x);
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,30],N[Pi,30]}];
Plot[{4(1-x*x)},{x,0,1}]
实验结果：
3.1415926519231265718,3.14159265358979323846264338
3.14159265358979323846264338328
图1 1/4个单位圆
结果分析：当数值积分法得到π的近似值为3.14159265358979323846264338328， 可以看出，用这种方法计算所得到的π值是相当精确的，n 越大，计算出来的扇形面积的近似值就越接近π的准确值。

步骤二、泰勒级数法计算π 利用反正切函数的泰勒级数
+--+-+-=--1
2)1(53a r c t a n 1
2153k x x x x x
k k 来计算π。

语句：T[x_,n_]:=Sum[(-1)^k*x^(2k+1)/(2k+1),{k,0,n}]; N[4*T[1,20000],20]//Timing
T[x_,n_]:=Sum[(-1)^k*x^(2k+1)/(2k+1),{k,0,n}]; Print[N[4*(T[1/2,260]+T[1/3,170]),150]]; Print[N[16*(T[1/5,110]-4*T[1/239,30]),150]]; Print[N[Pi,150]]
实验结果：
9.14Second,3.1416426510898869
3.14159265358979323846264338327950288419716939937510582494459230781640628620899862803482534211706798214808651230664709384460955058223172535940813
2.89054809346530980659035048572237571973428548091718877376781907690970580083540220107847652474250068362104652048128394634092219187032819003167814
3.14159265358979323846264338327950288419716939937510582494459230781640628620899862803482534211706798214808651230664709384460955058223172535940813
结果分析：从实验过程可以看出，这种方法花费的时间很长。

原因是当x=1时得到的arctan1的展开式收敛太慢。

要使泰勒级数收敛得快，容易想到，应当使x
的绝对值小于1，最好是远比1小。

例如，因为11
arctan1arctan arctan 23
=+，所
以我们可以计算出11
arctan ,arctan 23
的值，从而得到arctan1的值。

这样，就使得
收敛速度加快。

改进后可以看出，泰勒级数法得到的结果比数值分析法精确到小数点后更多位。

步骤三、蒙特卡罗法计算π
在数值分析法中，我们利用求单位圆的1/4面积来得到/4π，从而得到π。

单位圆的1/4是一个扇形，它是边长为1的单位正方形的一部分，单位正方形的面积11S =。

只要能够求出扇形的面积S 在正方形的面积中所占的比例1/k S S =，就能立即得到S ，从而得到π的值。

下面的问题归结为如何求k 的值，这就用到了一种利用随机数来解决此种问题的蒙特卡罗法，其原理就是
在正方形中随机的投入很多点，是所投的每个点落在正方形中每一个位置的机会均等，看其中有多少个点落在扇形内。

降落在扇形内的点的个数m 与所投店的总数n 的比可以近似的作为k 的近似值。

语句：
n=10000;p={}; Do[m=0;
Do[x=Random[];y=Random[]; If[x^2+y^2<=1,m++],{k,1,n}];
AppendTo[p,N[4m/n]],{t,1,10}];
Print[p];
Sum[p[[t]],{t,1,10}]/10 实验结果：
3.1528,3.1472,3.1276,3.134,3.1384,3.1516,3.1424,3.1664,3.1436,3.1
3.14668
结果分析：
从运行结果来看，蒙特卡罗法的计算结果为3.14668，虽然精确度不太高，但运行时间短，在很多场合下，特别是在对精确度要求不高的情况下很有用的。

步骤四、针对步骤三提出疑问：步骤三中我们发现当n=10000时，蒙特卡罗法的计算结果为3.14668，精确度不太高，那么对n 取不同的值，所得结果的精确度
会不会有变化？假如有变化，会有什么变化呢？
猜想：对n 取不同的值，所得结果的精确度应该会有变化，且当n 值越大，所得结果越精确。

现令n=1000；
语句：
n=1000;p={}; Do[m=0;
Do[x=Random[];y=Random[]; If[x^2+y^2<=1,m++],{k,1,n}];
AppendTo[p,N[4m/n]],{t,1,10}];
Print[p];
Sum[p[[t]],{t,1,10}]/10 实验结果：
3.16,3.132,3.08,3.156,3.144,3.184,3.156,3.116,3.092,3.
3.
令n=100000； 语句：
n=100000;p={}; Do[m=0;
Do[x=Random[];y=Random[]; If[x^2+y^2<=1,m++],{k,1,n}];
AppendTo[p,N[4m/n]],{t,1,10}];
Print[p];
Sum[p[[t]],{t,1,10}]/10 实验结果：
3.14,3.13172,3.13692,3.13752,3.140923.13852,3.13976,3.14572,3.14028,3.14
3.1
结果分析：
从运行结果来看，虽然蒙特卡罗法的计算结果的精确度不太高，但对n 取不同的值，所得结果的精确度有变化，且当n 值越大，所得结果越精确，这与我们的猜想完全一致。

步骤五、利用麦琴给出
239
1
arctan 51arctan
44
-=π
，推出π
=4(239
1
arctan
51arctan
4 )。

对比以上方法，这种简单的直接用公式求的π的方法要简单得多，所以用处更广。