一、实验目的探索精确计算π值的方法，并且比较不同方法之间的不同之处和优缺点。

掌握数值积分的辛普森公式。

二、问题描述1. 任务11) 用反正切函数的幂级数展开式结合有关公式求π，若要精确到40位、50位数字，试比较简单公式和Machin 公式所用的项数。

2) 验证公式111=arctan arctan arctan 4258π++ 试试此公式右端做幂级数展开完成任务1所需要的步数。

2. 任务2用数值积分计算π，分别用梯形法和Simpson 法精确到10位数字，用Simpson 法精确到15位数字。

3. 任务3用Monte Carlo 法计算π，除了加大随机数，在随机数一定时可重复算若干次后求平均值，看能否求得5位精确数字？设计方案用计算机模拟Buffon 实验4. 任务4利用积分20(1)!!sin !!2n n xdx n ππ-=⎰ ，n 为奇数 推导公式224422213352121n n n n π=-+ ……… 用此公式计算π的近似值，效果如何？5. 任务5利用学过的知识（或查阅资料），提出其他计算π的方法（先用你学过的知识证明），然后实践这种方法。

对你在实验中应用的计算π的方法进行比较分析。

6. 任务6e 是一个重要的超越数1e lim 1)n n n→∞=+（ 1111...2!!(1)!e e n n θ=++++++ 试用上述公式或其他方法近似计算e 。

三、问题解法1. 任务11) 根据幂级数展开的相关知识，易知：24122211(1)1n n x x x x--=-+-+-++……… 因为21(arctan )'1x x =+，故可以求得arctan x 的幂级数展开式为： 35211arctan (1)3521n n x x x x x n --=-+-+-+-……… 当x=1时，-11111--(-1)4352-1n n π=+⋯⋯++⋯ 当叠加了十万次以后得到结果π=3.141582654…只有五位有效数字，可见其精度与效率极低。

如果想要精确计算π的数值的话，非常有必要寻找改进以后的方法，这就引出了两个能够提高计算效率的公式——简单公式：11=arctan arctan 423π+ Machin 公式： 11=4arctan arctan 45239π- 对以上两式进行arctan 的幂级数展开可以非常快速的求得π比较精确的数值。

下面比较π精确到40位和50位数字时两个公式各需要计算多少项。

用简单公式得到40位有效数字需要叠加62项：3.141592653589793238462643383279502884197用简单公式得到50位有效数字需要叠加79项：3.1415926535897932384626433832795028841971693993751用Machin 公式得到40位有效数字需要叠加27项：3.141592653589793238462643383279502884197用Machin 公式得到50位有效数字需要叠加35项：3.1415926535897932384626433832795028841971693993751从上面简单的对比可以看出Machin 公式要优于简单公式，简单公式要优于不用公式的arctan 幂级数展开。

在得到相同精度的条件下，Machin 公式所需要的叠加步数要显著少于简单公式，并且在计算精度越高的情况下，优势越明显。

道理很简单，因为Machin 公式计算的收敛速度要显著快于普通公式。

简单公式决定收敛速度的是12n⎛⎫ ⎪⎝⎭，而Machin 公式决定收敛速度的是15n ⎛⎫ ⎪⎝⎭，因为15n ⎛⎫ ⎪⎝⎭的收敛速度快于12n ⎛⎫ ⎪⎝⎭，故Machin 公式计算pi 的时候收敛速度要快于普通公式。

所以Machin 公式比普通公式更加精确，并且在计算高精度的时候有更大优势。

2) 根据三角函数公式有：111tan(arctan arctan arctan )258111tan(arctan arctan )2581111tan(arctan arctan )2581123111123++++=-++==- 故111=arctan arctan arctan 4258π++ 用这个公式计算π值得到的结果与使用简单公式得到的结果完全相同：在62项得到40位有效数字，在79项得到50位有效数字。

2. 任务2 因为已知1201=41+x A dx π=⎰，故可以应用数值积分来计算π的近似值。

梯形法：12102[2(...))]n n y y y y y nπ-=+++++ Simpson 法：022********[()2(...)4(...)]3m n m y y y y y y y y m π--=+++++++++ 其中：211k ky x =+ 用matlab 计算后发现梯形法要将积分区间16811等分以后可以得到10位有效数字，而simpson 法只要22等分即可，精确到15位数字则需要做152等分。

可见Simpson 法要比梯形法更优越。

下面来分析一下Simpson 法优越的原因。

Simpson 公式法在达到相同精度的时候所用的区间分划远小于梯形法，这主要是Simpson 法所生成的曲线与原来函数的曲线吻合程度比较高。

3. 任务3Monte Carlo 法：在正方形 0<x <1，0<y<1上随机投大量的点，那么落在四分之一园内的点数数m 与在正方形内的点数n 之比m/n 应为π/4,故：π＝4 m/n 。

计算机模拟：产生区间[0,1]上数目为n 的一组随机数(x,y)，计算满足22y x +<1的数m 。

用matlab 模拟了910随机投点以后得到π的近似值为3.1417，只有四位有效数字，但是却需要很可观的计算时间，因此，想通过增加投点次数来实现π值的精确计算是不可行的。

下面尝试一下重复若干次以后取平均值能否提高精度。

我尝试了一下每一百次取一下平均，计算一百万次，结果依然不理想.有的时候能够达到5位精度，但是绝大多尝试都无法达到希望的5位精度。

要用计算机模拟Buffon 实验，可利用matlab 生成两个随机数列，分别模拟针中心到较近平行线距离y 及针与平行线夹角θ，将θsin ≤y 的情况的次数除以所有情况的总数就是2/π。

Buffon 投针实验精度依然不高，模拟一千万次以后得到π的近似值为3.1419。

4. 任务4首先要更正一下，原题中有一个小错误：公式成立的条件“n 为奇数”应该改为“n 为偶数”。

下面证明一下这个公式并导出一个新的近似计算π的公式。

根据分部积分公式，有：222200sin (1)cos sin nn xdx n x xdx ππ-=-⎰⎰ 将cosx 替换为sinx ，然后移项化简得：222001sin sin n n n xdx xdx n ππ--=⎰⎰ 这样就得到了一个递推公式。

易知：20sin 1xdx π=⎰2012dx ππ=⎰ 因此2020n-1)!!sin x n !!n-1)!!sin x n !!2n n dx n dx n πππ⎧=⎪⎪⎨⎪=⎪⎩⎰⎰ （为奇数（为偶数 下面推导π的近似值计算式：由n 为偶数时有：220(2)!!=sin 2(21)!!k k xdx k ππ-⎰ 又当k 的值很大时：lim 2lim 21k k k k →∞→∞=+，故可以用与之相邻的奇数代替这个偶数。

有：2120(2)!!sin (21)!!k k xdx k π+=+⎰ 这样就得到了π值的近似计算式：(2)!!(2)!!2(21)!!(21)!!k k k k π=-+ 也就是：224422213352121n n n n π=-+ ………用matlab 计算了一下n=10000的情况，结果是3.14151。

又计算了一下n=100000的情况，这时已经能感觉到程序明显的停顿，但是结果是3.14158，有效数字并没有增加。

从以上的结果可以看出，这种算法并不具有优越性，不仅不如辛普森积分公式，甚至不如Monte Carlo 法。

之所以得到这样的结果，主要是因为当n 取得很大的值的时候，2n/(2n+1)项已经十分接近1，这样对π的计算值逼近真实值贡献很小。

从这可以看出这种算法的特点是在开始的时候可以比较快速的逼近π的真实值，越到后来逼近真实值越困难。

也就是说在后面计算过程中，效率越来越低。

总之，这不是一个好的算法。

5. 任务5我在网上找到了一个计算π非常强的公式——拉马努金公式这个公式每迭代一次可以增加π计算值的14位精度，这是非常强大的。

从前几个任务的分析我们很容易得出，评价一个公式计算π近似值的优劣主要是看得到相同精度的情况下，谁用的循环次数更少。

如果想要使循环次数少的话就要让计算式的收敛速度变快，所以寻找一个好的计算π值的算法关键是要找到一个收敛速度快的公式。

Machin 公式的收敛是15指数收敛，收敛速度已经非常快了。

但是再看一下拉马努金公式，他的主要收敛项是41396，这个收敛速度是非常惊人的，因此用这个公式计算π的值效率是惊人的，可以很快逼近π的真实值。

在需要精确计算π值的场合应该优先考虑使用拉马努金公式。

在以上计算π值的方法中精确度最差的可以算是Monte Carlo 法了，如果想精确到5位数字还要在运气好的时候才能实现，但是这并不意味着在实际应用中Monte Carlo 法没有优点。

Monte Carlo 在金融等实际问题中有很广泛的应用，金融中的问题可能是一个很复杂的问题。

如果一个问题有500个甚至是更多的变量，求出其解析解是不可能的，这时候Monte Carlo 法就发挥了它巨大的优越性。

在用Monte Carlo 法计算π值的过程中，我们发现随机过程次数的增加似乎并不能非常显著的增加π值结果的精确度。

在matlab 的环境下，即使有无穷次的随机过程依然无法得到π的精确值。

因为matlab 生成的随机数的位数有限，因此生成的随机数中无理数是无法包含的，这样也就是说生成的随机数列其实是离散的。

也就是在一个平面网格上投点，而不是在平面圆域上投点，这样得到的π值精度是有限的。

6. 任务6超越数e 的算法比较简单，就是指数函数的泰勒级数展开。

因为展开以后阶乘的倒数，因此收敛速度也非常快，可以很快得到e 的近似值。

四、附录——实验中用到的matlab 源码1. Arctan 幂级数展开法function y=tpi(k)for n=1:ka(n)=(-1)^(n-1)/(2*n-1);endy=vpa(4*sum(a));2. 简单公式function f=simple(k,j)for n=1:ka(n)=(-1)^(n-1)*(1/2)^(2*n-1)/(2*n-1)+(-1)^(n-1)*(1/3)^(2*n-1)/(2*n-1);end f=vpa(4*sum(a),j);3.Machin公式function y=machin(k)for n=1:ka(n)=((-1)^(n-1)*(1/5)^(2*n-1)/(2*n-1))*4-(-1)^(n-1)*(1/239)^(2*n-1)/ (2*n-1);endy=vpa(4*sum(a));4.梯形法function f=trapezia(n,j)for k=0:nx(k+1)=k/n;y(k+1)=1/(1+x(k+1)^2);endvpa(f=2/n*(2*sum(y)-y(1)-y(n)),j)5.Simpson法function f=simpson(n,j)for k=0:2*n;x(k+1)=k/(2*n);y(k+1)=1/(1+x(k+1)^2);endfor i=1:nif i<=n-1a(i)=y(2*i+1);endb(i)=y(2*i);endf=vpa(2/(3*n)*(y(1)+y(2*n+1)+2*sum(a)+4*sum(b)),j)6.Monte Carlo法function y=monte(k)m=0;for n=1:kif rand(1)^2+rand(1)^2<=1m=m+1;end;end;4*m/k7.平均Monte Carlo法function y=monte2(k,t)for r=1:tm=0;for n=1:kif rand(1)^2+rand(1)^2<=1m=m+1;end;end;g(r)=4*m/k;endsum(g)/t8.模拟Buffon投针实验function f=buffon(n)m=0;for k=1:nif rand(1)<=sin(rand(1)*pi/2) m=m+1;endendf=n/m*29.Willas公式法function f=fpi(n)t=1;for i=1:na(i)=2*i/(2*i-1);b(i)=2*i/(2*i+1);t=t*a(i)*b(i);endf=2*t;10.超越数e的泰勒展开function f=normale(n)e=0;a(1)=1;for i=1:na(i+1)=1/(prod(1:i));endf=sum(a);。