计算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和阶数n

r= taylor(f,n,v,a),指定自变量v、阶数n，计算f在a的级数

(3)级数法

由于利用arctanx的幂级数展开法的收敛较慢，可采用公式

的计算来求pi值。

(4)特殊公式（BBP）

四、设计算法、编程求解

4.1数值积分法

梯形公式Matlab代码：

format long

x=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)

MAtlab结果：

数值积分法（梯形公式）Pi值

n=10 3.104518326248318

n=100 3.140417031779045

n=500 3.141487477002141

n=1000 3.141555466911028

n=10000 3.141591477611324

4.2幂级数法Matlab代码：

(1)

format long

syms x

f=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)

结果：

级数法（arctanx）Pi值

n=10 3.339682539682540

n=100 3.121594652591011

n=500 3.137592669589472

n=1000 3.139592655589785

n=10000 3.141392653591792

(2)

format long

syms x

f=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.141592682404399

级数法N=10的pi值

3.339682539682540

3.141592682404399

4.3迭代法Matlab代码：

format long

a=1;b=1/sqrt(2);s=1/2;

for n=1:1:10

n,

y=a;

a=(a+b)/2;

b=sqrt(b*y);

c=a^2-b^2;

s=s-2^n*c;

pi=2*a^2/s

end

结果：

迭代法Pi值

n=1 3.187672642712109

n=2 3.141680293297657

n=3 3.141592653895460

n=4 3.141592653589831

n=5 3.141592653589880

n=6 3.141592653589978

n=7 3.141592653590173

n=8 3.141592653590564

n=9 3.141592653591346

n=10 3.141592653592909

4.4蒙特卡罗方法Matlab代码：

format long

s=0;

n=10; % n=100; n=1000; n=10000; n=100000; n=1000000

for i=1:n

a=rand(1,2);

if a(1)^2+a(2)^2<=1

s=s+1;

end

end

pi=4*s/n

结果：

随机模拟法Pi值

n=10 2

n=100 3.240000000000000

n=1000 3.132000000000000

n=10000 3.135600000000000

n=100000 3.138480000000000

n=1000000 3.140024000000000 4.5 BBP公式

format long

syms x

y=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:10

s1=eval(y);

x,s=s+s1

end

结果：

BBP公式Pi值

n=0 3.133333333333333

n=1 3.141422466422466

n=2 3.141587390346582

n=3 3.141592457567436

n=4 3.141592645460337

n=5 3.141592653228088

n=6 3.141592653572881

n=7 3.141592653588973

n=8 3.141592653589752

n=9 3.141592653589791

n=10 3.141592653589793

五、分析求解结果

由上表可知，蒙特卡罗方法计算出的pi值与真实值的误差相差较大并且收敛速度很慢；对于级数法，但由于所选择的的级数方法、公式不同，得到的结果也就不同，收敛速度较慢，而的收敛速度就较快；数值积分法和迭代法准确度较高，但数值积分法的收敛速度没有迭代法快、精度高，所以一般情况下采用迭代法求近似值较准确。对于特殊BBP公式，其打破了传统的计算方法，可以直接计算pi的任何第n位数，而不是先计算前面的n-1位数。随着学习研究，利用特殊的公式计算明显地提高了pi的计算值的精确度。

六、参考文献

[1]刘慧颖.MTALAB R2007基础教程[M].北京：清华大学出版社，2008.7,200

[2]李庆扬，王能超，易大义.数值分析[M]. 北京：清华大学出版社，2008.7，187-190

[3]数学实验课件

PI值计算方法1

计算圆周率的一些公式 -|waruqi 发表于 2005-12-8 9:24:00 Machin公式 这个公式由英国天文学教授John Machin于1706年发现。他利用这个公式计算到了100位的圆周率。Machin公式每计算一项可以得到1.4位的十进制精度。因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。 还有很多类似于Machin公式的反正切公式: pi/4=arctg(1/2)+arctg(1/5)+ arctg(1/8) 1844.达塞利 = arctg(1/2)+ arctg(1/3) =2 arctg(1/3)+ arctg(1/7) =12 arctg(1/18)+8 arctg(1/57)-5 arctg(1/239) 在所有这些公式中，Machin公式似乎是最快的了。虽然如此，如果要计算更多的位数，比如几千万位，Machin公式就力不从心了。下面介绍的算法，在PC机上计算大约一天时间，就可以得到圆周率的过亿位的精度。这些算法用程序实现起来比较复杂。因为计算过程中涉及两个大数的乘除运算，要用FFT(Fast Four ier Transform)算法。FFT可以将两个大数的乘除运算时间由O(n2)缩短为O(nl og(n))。 （FFT算法不在此文讲诉）

Ramanujan公式 1914年，印度数学家Srinivasa Ramanujan在他的论文里发表了一系列共1 4条圆周率的计算公式，这是其中之一。这个公式每计算一项可以得到8位的十进制精度。1985年Gosper用这个公式计算到了圆周率的17,500,000位。 1989年，David & Gregory Chudnovsky兄弟将Ramanujan公式改良成为： 这个公式被称为Chudnovsky公式，每计算一项可以得到15位的十进制精度。1994年Chudnovsky兄弟利用这个公式计算到了4,044,000,000位。Chudnovsky 公式的另一个更方便于计算机编程的形式是： AGM(Arithmetic-Geometric Mean)算法 Gauss-Legendre公式： 初值： 重复计算：

计算pi

一、实验目的 探索精确计算π值的方法，并且比较不同方法之间的不同之处和优缺点。掌握数值积分的辛普森公式。 二、问题描述 1. 任务1 1) 用反正切函数的幂级数展开式结合有关公式求π，若要精确到40位、50位数 字，试比较简单公式和Machin 公式所用的项数。 2) 验证公式 111=arctan arctan arctan 4 258π ++ 试试此公式右端做幂级数展开完成任务1所需要的步数。 2. 任务2 用数值积分计算π，分别用梯形法和Simpson 法精确到10位数字，用Simpson 法精确到15位数字。 3. 任务3 用Monte Carlo 法计算π，除了加大随机数，在随机数一定时可重复算若干次后求平均值，看能否求得5位精确数字？ 设计方案用计算机模拟Buffon 实验 4. 任务4 利用积分 2 0(1)!!sin !!2 n n xdx n π π-=? ，n 为奇数 推导公式 224422213352121 n n n n π=-+ ……… 用此公式计算π的近似值，效果如何？ 5. 任务5 利用学过的知识（或查阅资料），提出其他计算π的方法（先用你学过的知识证明），然后实践这种方法。 对你在实验中应用的计算π的方法进行比较分析。 6. 任务6 e 是一个重要的超越数 1e lim 1)n n n →∞=+（ 1111...2!!(1)! e e n n θ =++++++ 试用上述公式或其他方法近似计算e 。

三、问题解法 1. 任务1 1) 根据幂级数展开的相关知识，易知： 24122211(1)1n n x x x x --=-+-+-++……… 因为2 1(arctan )'1x x =+，故可以求得arctan x 的幂级数展开式为： 35 211arctan (1)3521 n n x x x x x n --=-+-+-+-……… 当x=1时， -11111--(-1)4352-1 n n π=+??++? 当叠加了十万次以后得到结果π=3.141582654…只有五位有效数字，可见其精度与效率极低。如果想要精确计算π的数值的话，非常有必要寻找改进以后的方法，这就引出了两个能够提高计算效率的公式—— 简单公式： 11=arctan arctan 4 23 π + 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) 根据三角函数公式有：

PI的计算算法

上机7：文件操作 一、问题： 古人计算pi （3.1415926535 89793 23846……），一般是用割圆法。即用圆的内接或外切正多边形来逼近圆的周长。大约公元前1200年，中国人计算到小数点后1位；Archimedes 在公元前250年用正96边形得到pi 小数点后3位的精度；公元263年，刘徽用正3072边形计算到小数点后5位 ； 公元480年，祖冲之计算到小数点后7位；Ludolph Van Ceulen 在1609年用正262边形得到了35位 精度；1706年，John Machin 计算到小数点后100位。1874年，William Shanks 计算到小数点后707位（527位正确）。这种基于几何的算法计算量大，速度慢，吃力不讨好。 1973年Guilloud & Bouyer 用 CDC 7600计算机计算到1,001,250位；1989年Kanada & Tamurar 用 HITACHI S-820/80计算机计算到1,073,741,799位；1999年Takahashi & Kanada 用 HITACHI SR8000计算机计算到206,158,430,000位；pi 的最新计算纪录由两位日本人Daisuke Takahashi 和Yasumasa Kanada 所创造。他们在日本东京大学的IT 中心，以Gauss-Legendre 算法编写程序，利用一台每秒可执行一万亿次浮点运算的超级计算机，从日本时间1999年9月18日19:00:52起，计算了37小时21分04秒，得到了pi 的 206,158,430,208(3*236)位十进制精度，之后和他们于1999年6月27日以Borwein 四次迭代式计算 了46小时得到的结果相比，发现最后45位小数有差异，因此他们取小数点后206,158,430,000位的p 值为本次计算结果。这一结果打破了他们于1999年4月创造的68,719,470,000位的世界纪录。随着数学的发展，数学家们在进行数学研究时有意无意地发现了许多计算pi 的公式。 Machin 公式 1 15239164arg arctg tg π=- 3 5 7 211()(1)35721 n n x x x x a rctg x x n --=-+-++-- 这个公式由英国天文学教授John Machin 于1706年发现。他利用这个公式计算到了100位的pi 值。Machin 公式每计算一项可以得到1.4位的十进制精度。因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。 Bailey-Borwein-Plouffe 算法 0142 1 1 ()1681848586n n n n n n π∞ ==---++++∑ 这个公式简称BBP 公式，由David Bailey, Peter Borwein 和Simon Plouffe 于1995年共同发表。它打破了传统的p 的算法，可以计算p 的任意第n 位，而不用计算前面的n-1位。这为p 的分布式计算提供了可行性。 因此使用这个公式计算，与计算常数e （2.718282……）的公式非常类似，编程很容易实现。 二、要求： 编写程序计算常数π到小数点后任意位，并将数据写入文件中。 参考程序： #include #include #include #include

数学实验：怎样计算圆周率

怎样计算 姓名： 学号 班级：数学与应用数学4班

实验报告 实验目的：自己尝试利用Mathematica软件计算的近似值，并学会计算的近似值的方法。 实验环境：Mathematica软件 实验基本理论和方法： 方法一：数值积分法（单位圆的面积是，只要计算出单位圆的面积也就计算出了的值） 其具体内容是：以单位圆的圆心为原点建立直角坐标系，则单位圆在第一象限内的部分G是一个扇形， 由曲线（）及坐标轴围成，它的面积是，算出了S的近似值，它的4倍就是的近似值。而怎样计算扇形G的面积S的近似值呢?如图

图一 扇形G中，作平行于y轴的直线将x轴上的区间[0，1]（也就是扇形在x轴上的半径）分成n等份（n=20），相应的将扇形G分成n个同样宽度1/n的部分（）。每部分是一个曲边梯形：它的左方、右方的边界是相互平行的直线段，类似于梯形的两底；上方边界是一段曲线，因此称为曲边梯形。如果n很大，每个曲边梯形的上边界可以近似的看成直线段，从而将近似的看成一个梯形来计算它的面积；梯形的高（也就是它的宽度）h=1/n，两条底边的长分别是和，于是这个梯形面积可以作为曲边梯形面积的近似值。所有这些梯形面积的和T就可以作为扇形面积S的近似值： n越大，计算出来的梯形面积之和T就越接近扇形面积S，而4T就越接近的准确值。 方法二：泰勒级数法 其具体内容是：利用反正切函数的泰勒级数 计算。 方法三：蒙特卡罗法

其具体内容是：单位正方形的面积=1，只要能够求出扇形G 的面积S在正方形的面积中所占的比例，就能立即得到S，从而得到的值。而求扇形面积在正方形面积中所占的比例k的值，方法是在正方形中随机地投入很多点，使所投的每个点落在正方形中每一个位置的机会均等，看其中有多少个点落在扇形内。将落在扇形内的点的个数m与所投的点的总数n的比可以作为k的近似值。能够产生在区间[0，1]内均匀分布的随机数，在Mathematica中语句是 Random[ ] 产生两个这样的随机数x，y，则以（x，y）为坐标的点就是单位正方形内的一点P，它落在正方形内每一个位置的机会均等。P落在扇形内的充分必要条件是。这样利用随机数来解决数学问题的方法叫蒙特卡罗法。 实验内容、步骤及其结果分析： 问题1：在方法一中，取n=1000，通过计算图一中扇形面积计算的的近似值。 分析：图一中的扇形面积S实际上就是定积分。 与有关的定积分很多，比如的定积分

计算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和阶数n r= taylor(f,n,v,a),指定自变量v、阶数n，计算f在a的级数 (3)级数法

由于利用arctanx的幂级数展开法的收敛较慢，可采用公式 的计算来求pi值。 (4)特殊公式（BBP） 四、设计算法、编程求解 4.1数值积分法 梯形公式Matlab代码： format long x=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) MAtlab结果： 数值积分法（梯形公式）Pi值 n=10 3.104518326248318 n=100 3.140417031779045 n=500 3.141487477002141 n=1000 3.141555466911028 n=10000 3.141591477611324 4.2幂级数法Matlab代码： (1) format long syms x f=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) 结果： 级数法（arctanx）Pi值 n=10 3.339682539682540 n=100 3.121594652591011 n=500 3.137592669589472 n=1000 3.139592655589785

如何计算π的值

1、蒙特卡罗（Monte Carlo ）法 思想： 取一正方形A ，以A 的一个顶点为圆心，A 的边长为半径画圆，取四分之一圆（正方形内的四分之一圆）为扇形B 。已知A 的面积，只要求出B 的面积与A 的面积之比B A S k S =，就能得出B S ，再由B 的面积为圆面积的四分之一，利用公式2=S R π圆即可求出π的值。因此，我 们的目的就是要找出k 的值。 可以把A 和B 看成是由无限多个点组成，而B 内的所有点都在A 内。随机产生n 个点，若落在B 内的有m 个点（假定A 的边长为1，以扇形圆心为坐标系原点。则只要使随机产生横纵坐标x 、y 满足221x y +≤的点，就是落在B 内的点），则可近似得出k 的值，即m k n =，由此就可以求出π的值。 程序（1）： i=1;m=0;n=1000; for i=1:n a=rand(1,2); if a(1)^2+a(2)^2<=1 m=m+1; end end p=vpa(4*m/n,30) 程序运行结果：

p = 2、泰勒级数法 思想： 反正切函数arctan x 的泰勒级数展开式为： 35 211arctan (1)3521 k k x x x x x k --=-+-???+-+???- 将1x =代入上式有 1 111arctan11(1)43521 n n π-==-+-???+--. 利用这个式子就可以求出π的值了。 程序（2）： i=1;n=1000;s=0; for i=1:n s=s+(-1)^(i-1)/(2*i-1); end p=vpa(4*s,30) 程序运行结果： p = 当取n 的值为10000时，就会花费很长时间，而且精度也不是很高。原因是1x =时，arctan1的展开式收敛太慢。因此就需要找出一个x 使得arctan x 收敛更快。

pi的计算 数值分析论文

古人计算圆周率，一般是用割圆法。即用圆的内接或外切正多边形来逼近圆的周长。Archimedes用正96边形得到圆周率小数点后3位的精度；刘徽用正3072边形得到5位精度；Ludolph Van Ceulen用正262边形得到了35位精度。这种基于几何的算法计算量大，速度慢，吃力不讨好。随着数学的发展，数学家们在进行数学研究时有意无意地发现了许多计算圆周率的公式。下面挑选一些经典的常用公式加以介绍。除了这些经典公式外，还有很多其它公式和由这些经典公式衍生出来的公式，就不一一列举了。 1、 Machin公式 [这个公式由英国天文学教授John Machin于1706年发现。他利用这个公式计算到了100位的圆周率。Machin公式每计算一项可以得到1.4位的十进制精度。因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。 Machin.c 源程序 还有很多类似于Machin公式的反正切公式。在所有这些公式中，Machin公式似乎是最快的了。虽然如此，如果要计算更多的位数，比如几千万位，Machin公式就力不从心了。下面介绍的算法，在PC机上计算大约一天时间，就可以得到圆周率的过亿位的精度。这些算法用程序实现起来比较复杂。因为计算过程中涉及两个大数的乘除运算，要用FFT(Fast Fourier Transform)算法。FFT可以将两个大数的乘除运算时间由O(n2)缩短为O(nlog(n))。 2、 Ramanujan公式 1914年，印度数学家Srinivasa Ramanujan在他的论文里发表了一系列共14条圆周率的计算公式，这是其中之一。这个公式每计算一项可以得到8位的十进制精度。1985年Gosper用这个公式计算到了圆周率的17,500,000位。 1989年，David & Gregory Chudnovsky兄弟将Ramanujan公式改良成为： 这个公式被称为Chudnovsky公式，每计算一项可以得到15位的十进制精度。1994年Chudnovsky兄弟利用这个公式计算到了4,044,000,000位。Chudnovsky公式的另一个更方便于计算机编程的形式是： 3、AGM(Arithmetic-Geometric Mean)算法 Gauss-Legendre公式： 这个公式每迭代一次将得到双倍的十进制精度，比如要计算100万位，迭代20次就够了。1999年9月Takahashi和Kanada用这个算法计算到了圆周率的206,158,430,000位，创出新的世界纪录。 4、Borwein四次迭代式： 这个公式由Jonathan Borwein和Peter Borwein于1985年发表，它四次收敛于圆周率。

pi的无穷级数算法

3π的计算模型—π的无穷级数算法

思想：利用一些特殊函数的幂级数展开式计算圆周率。 根据幂级数展开的相关知识，易知： 241222 1 1(1)1n n x x x x --=-+-+-++……… 因为2 1 (arctan )'1x x =+，故可以求得arctan x 的幂级数展开式为： 3521 1arctan (1)3521n n x x x x x n --=-+-+-+-……… 当x=1时， -11111--(-1)4352-1 n n π =+??++? 用Matlab 计算： 创建m 文件calpi1.m ，内容如下： function y=calpi1(k) for n=1:k a(n)=(-1).^(n-1)./(2*n-1); end ; 4*sum(a) 在命令窗口中输入如下命令：

它是一个与π有关的无穷级数，实际计算时，我们只能使用有限项。如果取级数前n 项之和作为π的近似值，其误差为 ｜｜≤ 为了保证误差不超过 ，就要取级数的前20000项进行计算，计算 量之大可想而知。现在看来，计算π的级数有明显的缺点：级数收敛太慢，计算量过大。其原因是｜x ｜偏大。如果想要精确计算π的数值的话，非常有必要寻找改进以后的方法，这就引出了两个能够提高计算效率的公式—— （1）欧拉公式：11=arctan arctan 4 23 π + 现取，令α=，显见0＜α＜， 记β=，而 = =， 所以β= ，就是11=arctan arctan 4 2 3 π + （2）马庭公式：11=4arctan arctan 4 5239 π - 可取，令，则 ， ，故4α≈，再令 ， 即，而，就有11=4arctan arctan 4 5 239 π - 将欧拉公式和马庭公式与arctanx 的泰勒级数相结合，会加快该级数的收敛速度，具有很强的实用性。

圆周率的计算数学实验报告

数学实验报告（二） 一、实验题目：圆周率的计算 二、实验目的： 1.用多种方法计算圆周率的值； 2.通过实验来说明各种方法的优劣； 3.尝试提出新的计算方法。 三、实验内容和方法： 1.古典方法： 用圆内接正多边形和圆外切正多边形来逼近 以阿基米德的圆内接96边形和圆外切96边形逼近为例 已知：sin <

y=4*x 得出当k=10时，π≈3.232315809405593 编写程序 syms k x=symsum((-1)^k/(2*k+1),k,0,20) y=4*x 得出当k=20时，π≈3.189184782277595 依次，加大k 的值 K=50，π≈3.161198612987050 K=100，π≈3.151493401070990 K=200，π≈3.146567747182986e+159 … (2).沃里斯(Wallis)方法 ∏∞=??? ??+?-?=??? ??????? ??? ???? ????=1122122276565 43432122k k k k k π 编写程序： format long x=1; for k=1:10 x=x*(2*k/(2*k-1)*2*k/(2*k+1)); end y=2*x 得k=10时，π≈3.067703806643498 增加k 的值 K=20，π≈3.103516961539230 K=50，π≈3.126078********* K=100，π≈3.133787********* K=10000，π≈3.141514118681864 K=1000000，π≈3.141591868191880 … (3).利用公式3 1arctan 21arctan 1arctan 4+==π 推出π=4(31 arctan 21 arctan +) 编写程序： syms n; f1=(-1)^(n-1)*(1/2)^(2*n-1)/(2*n-1); f2=(-1)^(n-1)*(1/3)^(2*n-1)/(2*n-1); ans1=symsum(f1,n,1,79);

积分法计算pi值

Pi值计算(积分法) 程序清单： #include #include #include static long num_step; const int Threadnum=10; double step,pi=0,sum=0.0; CRITICAL_SECTION g; DWORD WINAPI threadFunc(LPVOID pParam) { int num=*((int *)pParam); double x,sum1=0; step=1.0/(double)num_step; for(int i=num;i>num_step; LARGE_INTEGER temp; double dFreq; //系统时钟 QueryPerformanceFrequency(&temp); dFreq=(double)temp.QuadPart;// 获得计数器的时钟频率LONGLONG Start,End; QueryPerformanceCounter(&temp); Start=temp.QuadPart;// 获得初始计数器数值 HANDLE hthread[Threadnum]; InitializeCriticalSection(&g); for(int i=0;i

数学实验怎样计算圆周率

怎样计算 姓名: 学号 班级:数学与应用数学4班 实验报告 实验目的:自己尝试利用Mathematica软件计算的近似值,并学会计算的近似值的方法。 实验环境:Mathematica软件 实验基本理论与方法: 方法一:数值积分法(单位圆的面积就是,只要计算出单位圆的面积也就计算出了的值) 其具体内容就是:以单位圆的圆心为原点建立直角坐标系,则单位圆在第一象限内的部分G就是一个扇 形, 由曲线()及坐标轴围成,它的面积就是,算出了S的近似值,它的4倍就就是的近似值。而怎样计算扇形G的面积S的近似值呢?如图

图一 扇形G中,作平行于y轴的直线将x轴上的区间[0,1](也就就是扇形在x轴上的半径)分成n等份(n=20),相应的将扇形G分成n个同样宽度1/n的部分()。每部分就是一个曲边梯形:它的左方、右方的边界就是相互平行的直线段,类似于梯形的两底;上方边界就是一段曲线,因此称为曲边梯形。如果n很大,每个曲边梯形的上边界可以近似的瞧成直线段,从而将近似的瞧成一个梯形来计算它的面积;梯形的高(也就就是它的宽度)h=1/n,两条底边的长分别就 是与,于就是这个梯形面积 可以作为曲边梯形面积的近似值。所有这些梯形面积的与T就可以作为扇形面积S的近似值: n越大,计算出来的梯形面积之与T就越接近扇形面积S,而4T就越接近的准确值。 方法二:泰勒级数法

其具体内容就是:利用反正切函数的泰勒级数 计算。 方法三:蒙特卡罗法 其具体内容就是:单位正方形的面积=1,只要能够求出扇形G 的面积S在正方形的面积中所占的比例,就能立即得到S,从而得到的值。而求扇形面积在正方形面积中所占的比例k的值,方法就是在正方形中随机地投入很多点,使所投的每个点落在正方形中每一个位置的机会均等,瞧其中有多少个点落在扇形内。将落在扇形内的点的个数m与所投的点的总数n的比可以作为k 的近似值。能够产生在区间[0,1]内均匀分布的随机数,在Mathematica 中语句就是 Random[ ] 产生两个这样的随机数x,y,则以(x,y)为坐标的点就就是单位正方形内的一点P,它落在正方形内每一个位置的机会均等。P落在扇形内的充分必要条件就是。这样利用随机数来解决数学问题的方法叫蒙特卡罗法。 实验内容、步骤及其结果分析: 问题1:在方法一中,取n=1000,通过计算图一中扇形面积计算的 的近似值。

Mathematica计算Pi的值

Mathematica计算π的值 姓名: 学号: 班级: 实验目的 学习使用Mathematica软件的一些基本功能计算π的值，以下通过三种不同的方法求解π： 1.数值积分法 2.泰勒级数法 3.蒙特卡洛（Monte Carlo）方法 实验的基本原理和方法 1.Mathematica中常用绘图函数Plot在绘制高次函数时的方法； 2.计算圆周率π的数值积分法、泰勒级数法、蒙特卡洛（Monte Carlo）方法，并且利用特定的公式来计算圆周率π。 实验的内容和步骤 （1）数值积分法计算π 半径为1的圆称为单位圆，它的面积等于π。只要计算出单位圆的面积，就算出了π。在坐标轴上画出以圆点为圆心，以1为半径的单位圆（如下图），则这个单位圆在第一象限的部分是一个扇形，而且面积是单位圆的1/4，于是，我们只要算出此扇形的面积，便可以计算出π。 1.1．Mathematica输入如下： Plot[{4(1-x*x)},{x,0,1}] 图1 在计算扇形面积时，很容易想到使用数学分析中积分的方法，第一象限中的

扇形由曲线])1,0[(12∈-=x x y 及两条坐标轴围成，实际操作中，我们不能准确地计算它的面积，于是就通过分割的方法，将其划分为许多小的梯形，通过 利用梯形的面积近似于扇形面积来计算4 1102π =-=?dx x S 。 利用Mathematica 编程计算上式: 运行结果如下： 图2 从而得到 的近似值为3.14159265358979323846264338328，可以看出，用这种方法计算所得到的值π是相当精确的。n 越大，计算出来的扇形面积的近似值就越接近π的准确值。 2.泰勒级数法计算π 反正切函数的泰勒级数 +--+-+-=--1 2)1(53a r c t a n 1 2153k x x x x x k k 计算π，实验运行如下：

计算PI值

计算PI值 （作者：草原飞狼 2015年10月28日）声明：仅供学习与交流使用，高手请飘过，谢谢！ 布局 运行界面（1）

运行界面（2） 源代码如下： 主程序代码： Private Sub Command1_Click() Rem 计算PI值，核心算法是累加求和，注意循环结构的执行条件Dim sum As Double Dim n As Double Dim temp As Double Dim my_mubiao As Double sum = 0 temp = 1 n = 1 k = 0 '记录循环次数 my_mubiao = 0.000000000001 Do While temp >= my_mubiao '循环执行条件

sum = sum + temp n = n + 1 temp = 1 / (n * n) '得到下一个temp值 k = k + 1 Loop Text1.Text = Str(k) Text2.Text = Str(Sqr(6 * sum)) End Sub Private Sub Command2_Click() Rem 退出 Dim int_msg As Integer int_msg = MsgBox("单击“是”退出程序，单击“否”返回程序！", vbYesNo + vbQuestion + vbDefaultButton1, "退出提示") If int_msg = vbYes Then Unload Me Else MsgBox "返回程序，继续进行！", vbOKOnly + vbInformation, "程序返回" End If End Sub Private Sub Command3_Click() Rem 清空 Text1.Text = "" Text2.Text = "" End Sub

如何计算π的值(MATLAB)

如何计算π的值 1、蒙特卡罗（Monte Carlo ）法 思想： 取一正方形A ，以A 的一个顶点为圆心，A 的边长为半径画圆，取四分之一圆（正方形内的四分之一圆）为扇形B 。已知A 的面积，只要求出B 的面积与A 的面积之比B A S k S = ，就能得出B S ，再由B 的面积为圆面积的四分之一，利用公式2=S R π圆即可求出π的值。因此，我们的目的就是要找出k 的值。 可以把A 和B 看成是由无限多个点组成，而B 内的所有点都在A 内。随机产生n 个点，若落在B 内的有m 个点（假定A 的边长为1，以扇形圆心为坐标系原点。则只要使随机产生横纵坐标x 、y 满足221x y +≤的点，就是落在B 内的点），则可近似得出k 的值，即m k n =，由此就可以求出π的值。 程序（1）： i=1;m=0;n=1000; for i=1:n a=rand(1,2); if a(1)^2+a(2)^2<=1 m=m+1; end end p=vpa(4*m/n,30)

程序运行结果： p = 3.14000000000000000000000000000 2、泰勒级数法 思想： 反正切函数arctan x 的泰勒级数展开式为： 3521 1arctan (1)3521 k k x x x x x k --=-+-???+-+???- 将1x =代入上式有 1111 arctan11(1)43521 n n π -==-+-???+--. 利用这个式子就可以求出π的值了。 程序（2）： i=1;n=1000;s=0; for i=1:n s=s+(-1)^(i-1)/(2*i-1); end p=vpa(4*s,30) 程序运行结果： p = 3.14059265383979413499559996126 当取n 的值为10000时，就会花费很长时间，而且精度也不是很高。原因是1x =时，arctan1的展开式收敛太慢。因此就需要找出一个x 使得arctan x 收敛更快。

实验二 怎样计算Pi

数学实验实验

报 告 学院：数学与统计学院 班级：数学与应用数学3班 学号：201370010314 姓名：康萍 时间：2016.04.05 实验二怎样计算π 一、实验目的 分别用下列三种方法计算π的近似值，并比较三种方法的精确度： 数值积分法：通过使用Mathematica7.0编写梯形公式和辛普森公式的程序语言计算π。 泰勒级数法：利用反正切函数泰勒级数计算π。 蒙特卡罗（Monte Carlo）法：通过使用Mathematica7.0编写蒙特卡罗公式的程序语言来计算π。 二、实验环境 基于Windows环境下的Mathematica7.0软件。 三、实验的基本理论和方法 1、数值积分法

以单位圆的圆心为原点建立直角坐标系，则单位圆在第一象限内的部分G 是一个扇形，由曲线])1,0[(12∈-=x x y 及两条坐标轴围成，它的面积4 π = S 。算 出了S 的近似值，它的4倍就是π的近似值。而扇形面积S 实际上就是定积分 4 11 2π = -? dx x 。 与π有关的定积分有很多，比如211x +的定积分4 11102π =+?dx x 就比21x -的定积分更容易计算，更适合于用来计算π。 一般地，要计算定积分()dx x f b a ?，也就是计算曲线()x f y =与直线 b x a x y ===,,0所围成的曲边梯形G 的面积S 。为此，用一组平行于y 轴的直线 ()b x x x x x a n i x x n n i =<<<<<=-≤≤=-1210,11 将曲边梯形T 分成n 个小曲边梯形，总面积S 分成这些小曲边梯形的面积之和。如果取n 很大，使每个小曲边梯形的宽度都很小，可以将它上方的边界()()i i x x x x f ≤≤-1近似的看作直线段，将每个小曲边梯形近似的看作梯形来求面积，就得到梯形公式。如果更准确些，将每个小曲边梯形的上边界近似的看作抛物线段，就得到辛普森公式。具体公式如下： 梯形公式 设分点11,,-n x x 将积分区间],[b a 分成n 等份，即 ()n i n a b i a x i ≤≤-+=0,/。 所有的曲边梯形的宽度都是()n a b h /-=。记()i i x f y = 则第i 个曲边梯形的面积i S 近似的等于梯形面积()h y y i i +-12 1 。将所有这些梯形面积加起来就得到 ?? ???? +++++-= -20121n n y y y y y n b a S 这就是梯形公式。 辛普森公式 仍用分点()11-≤≤-+=n i n a b i a x i 将区间[] b a ,分成n 等份，直线()11-≤≤=n i x x i 将曲边梯形分成n 个小曲边梯形。再做每个小区间[]i i x x ,1-的中点n a b i a x i -? ?? ??-+=- 212 1。将第i 个小曲边梯形的上边界()()i i x x x x f y ≤≤=-1近

π的计算方法有哪些

圆周率π的计算方法 圆周率π的计算方法，是一个饶有趣味，值得探讨的问题。最直观的计算方法自然是从几何上着手，历史上也正是如此，这便是割圆法。设一半径为1的圆，作这个圆的内接正n边形，用此正n边形的周长去近似圆的周长。显然当n→∞时，正n边形的周长就无限趋近于圆周长，求得正n边形周长后除以直径便求出了圆周率。 I. 从几何上观察，可知：正n边形周长随n递增而递增，但始终是个有限值。割法如图1： 向左转|向右转 图1 割圆法 设圆半径为1，令半弦长AB=2a，AC=2c，OG和OD分别是等腰△OAB和△OAC 的中线。则我们要做的只是求出c关于a的表达式c=c（a）.令GC=b，根据勾股定理有: 向左转|向右转 (1) 进而有 向左转|向右转 (2)

得到此式后，编写计算机程序就很容易了，C语言程序如下： #include #include main() {double a,b,c,d,pi;double sqrt(double);int i,j,n;a=0.5;b=0;c=0;d=0.5;scanf("%d",&n); for(i=1;i<=n;i++) {b=sqrt(1-a*a);c=(1-b)*0.5;d=sqrt(c);a=d;} j=pow(2,n)*3;pi=2*d*j;printf("%d\n",j);printf("%f\n",pi)} 这里有一个问题就是a的初值如何选择？显然越简单直观越好，而已知对于圆内接正六边形的每一条边长等于圆的半径。所以取a=0.5，程序中参数n是对正六边形分割的次数，d的作用是当输入n=0（正六边形）的时候，得到π=3，此所谓的“径圆一三”。将这个文件保存为文本，在linux下用“gcc -lm”命令编译后，打开编译后得到的文件就能执行。 在古代可没有电子计算机，而祖冲之利用割圆法算得圆周率介于3.1415926和 3.1415927之间，可见古人之伟大！ II. 上面的方法简单直观，但是缺点也很明显。计算机在底层只能做“加减乘除四则整数运算”，显然开根号运算还是要通过转化为整数运算（级数展开等）才最后到硬件级计算。那么我们能否直接用整数的四则运算得到π的值？有！而且方法是多样的，其中一种叫作“Wallis公式”，有几种表达方式。如下： 向左转|向右转 (3) 或 向左转|向右转 (4) 或 (5)

实验2π的计算

. 22||||2 12188 . 224 4242))21 (11()21(|||||| .22 . . 11121 111221212 2 21122N N ∈?=?=????===?=∈--=--=--+=+====++-+n a S a OE AB S S n a a a a a EF AF AE a S a n S a n n n OAE n n n n n n ， ，故又注意到，相应的边数，故， 注意到，边长和面积，则步单位圆内接多边形的第分别是 ，设值也就越接近圆周率近单位圆的面积，其数边形的面积就越来越接那么内接正多 正多边形的边数开始，逐步成倍地增加从单位圆的内接正方形Δ实验2：π的计算 学院：机械与动力工程学院 姓名：唐子彦 学号：515020910137 一、实验目的 通过求π的近似值，了解历史上计算π值的一些方法，包括刘徽割圆术、级数展开、数值积分和Monte Carlo 法等. 复习微积分中相关知识，比较它们的差异，了解计算方法对提高计算效率的意义. 二、实验内容 （一）利用鲁道夫割圆术计算π值 【问题】德国人鲁道夫一生计算圆周率，他同样是用圆的内接多边形逼近圆周. 不过，他是从正方形开始成倍增加边数. 试推导出他计算所采用的递推公式，然后求π的近似值到10位和20位有效数字. 【解】 图1.1

. 1 ,21 ,21 421 , 2 1112???>∈?==?????>∈--==---n n a n S n n a n a n n n n n 且，且， 的递推公式： 由此可得鲁道夫割圆术N N 为通过“鲁道夫法”计算π的近似值，现编写M 文件如下： 为求出指定有效数字位数的π值，还需编写M 文件如下： 运行结果如下： 图1.2 【答】由图1.2不难发现，“鲁道夫法”计算π值效率较高，n=16时即可得到π的10位有效数字，n=32时即可得到π的20位有效数字. 然而，由于中途过程涉及开方等运算，因此计算较为复杂繁琐 .

数学实验：怎样计算圆周率

怎样计算 ： 学号 班级：数学与应用数学4班

实验报告 实验目的：自己尝试利用Mathematica软件计算的近似值，并学会计算的近似值的方法。 实验环境：Mathematica软件 实验基本理论和方法： 方法一：数值积分法（单位圆的面积是，只要计算出单位圆的面积也就计算出了的值） 其具体容是：以单位圆的圆心为原点建立直角坐标系，则单位圆在第一象限的部分G是一个扇形， 由曲线（）及坐标轴围成，它的面积是，算出了S的近似值，它的4倍就是的近似值。而怎样计算扇形G的面积S的近似值呢?如图 图一

扇形G中，作平行于y轴的直线将x轴上的区间[0，1]（也就是扇形在x轴上的半径）分成n等份（n=20），相应的将扇形G分成n个同样宽度1/n的部分（）。每部分是一个曲边梯形：它的左方、右方的边界是相互平行的直线段，类似于梯形的两底；上方边界是一段曲线，因此称为曲边梯形。如果n很大，每个曲边梯形的上边界可以近似的看成直线段，从而将近似的看成一个梯形来计算它的面积；梯形的高（也就是它的宽度）h=1/n，两条底边的长分别是和，于是这个梯形面积可以作为曲边梯形面积的近似值。所有这些梯形面积的和T就可以作为扇形面积S的近似值： n越大，计算出来的梯形面积之和T就越接近扇形面积S，而4T就越接近的准确值。 方法二：泰勒级数法 其具体容是：利用反正切函数的泰勒级数 计算。 方法三：蒙特卡罗法 其具体容是：单位正方形的面积=1，只要能够求出扇形G的面积S 在正方形的面积中所占的比例，就能立即得到S，从而得到的值。而求扇形面积在正方形面积中所占的比例k的值，方法是在正方形中随机地投入很多点，使所投的每个点落在正方形中每一个位置的机会均等，看其中有多少个点落在扇形。将落在扇形的点的个数m与所投