开放性数学实验报告

(2016 / 2017学年 第 2学期)

题 目：基于MATLAB的圆周率近似计算

专 业 通信工程

学生姓名 杨 坤 冯著豪 周李鑫

班级学号 B16011115 B16011110 B16011124

指导教师 赵礼峰

指导单位 南京邮电大学理学院

日 期 2017/5/20

MATLAB圆周率的近似计算

B16011115 杨 坤 B16011110 冯著豪 B16011124周李鑫

摘要：圆周率(Pi)是圆的周长与直径的比值，一般用希腊字母π表示，是一个在数学及物理学中普遍存在的数学常数。π也等于圆形之面积与半径平方之比。是精确计算圆周长、圆面积、球体积等几何形状的关键值。 在分析学里，π可以严格地定义为满足sin x = 0的最小正实数x。

计算圆周率一直是很多人的追求。在电子计算机还没有发明的时候就有很多先贤用各种方法计算了圆周率的近似值最著名的应该是祖冲之，他计算出了圆周率的位数达到了小数点后七位。该记录在世界范围内保持了八百年。

之后圆周率的计算进入了分析法时期，这一时期人们开始利用无穷级数或无穷连乘积求π，摆脱可割圆术的繁复计算。无穷乘积式、无穷连分数、无穷级数等各种π值表达式纷纷出现，使得π值计算精度迅速增加。

在分析法的基础上，电子计算机的出现使得圆周率的计算精度大幅提高。计算圆周率已经成为评判超级计算机的性能指标的项目之一。

如今个人计算机的性能也达到了一个极高的程度。学习使用计算机计算圆周率可以帮助我们更好地学习matlab同时对数学也会有更深的理解。

关键词：圆周率计算；投点法；定积分计分法；幂级数；韦达公式

一、 问题分析

计算圆周率有很多方法，不同方法之间自然也有好坏之分。在强大的计算机性能的支持下，我们能使用不同的方法计算圆周率并且感受不同方法孰优孰劣。首先我们需要了解不同的计算方法是怎么计算圆周率的，然后使用matlab编写代码帮助我们实现算法，计算出圆周率。

二、 实验方法

1. 投点法：

投点法，顾名思义就是通过投点计算圆周率。在一个边长为1的正方形里以1为半径画一个四分之一圆，再向正方形里投点，在概率的学习中我们知道，大量地向这个正方形中投点时，在投的点足够多的前提下，落在四分之一扇形里的点与投的所有点的个数之比应该为扇形与整个正方形的面积之比。扇形的面积为四分之一圆，即1/4*pi,正方形的面积为1.设投n个点，落在扇形里的点的个数为count即可推出pi=4*(count/n)。

代码如下：

count=0;

ezplot('x^2+y^2=1',[0,1,0,1]),hold on ,grid on

n=10000;

for i=1:1:n

x=rand(1,1);

y=rand(1,1);

plot(x,y,'*'),hold on

pause(0.001)

if x^2+y^2<=1

count=count+1;

end

end

p=4*(count/n)

投200个点时多次运行分别获得以下结果：

3.3600 3.0600 3.1800 3.1400 3.1800

可见这时所得的结果并不稳定，且结果并不精确于是继续实验投一千个点。

多次实验后得出以下结果

3.1040 3.1120 3.2200 3.1520 3.0560 3.1280

此时的精确度仍然不尽人意，于是我们直接投了一万个点

此时得到了以下结果

3.1376 3.1355 3.1413 3.1415 3.1490 3.1457

此时基本达到了两位小数的精确度

于是我们更改投点个数得到以下结果

投100000个点时：3.1415 3.1414 3.1414 3.1413

当我们准备投更多点时发现投点法消耗太多计算机资源，已经无法投太多的点了。投点法计算圆周率利用了计算机的高性能，但是性能的使用效率并不太高。

2. 定积分积分法

0111+x2dx=π4 π=40111+x2dx

将区间[0,1]分成n等份，在每个小区间上，选中点为ξi，使用积分的方法计算π的近似值。

代码如下：

n=50;%等分积分区间数。

i=0:1/n:1;

s=0;

for k=1:length(i)-1

s=s+(1/(1+((i(k)+i(k+1))/2)^2))*1/n;

end

vpa(4*s，20)

设置不同的区间数，获得不同的近似值，对应关系如下：

区间数

50

100

200

300

400

500

600

近似值

3.1416

3.1416

3.1415947

3.14159357

3.14159317

3.14159298

3.14159288

相比于投点法，该方法以较高的效率使用了计算机资源达到了更高的精度，当n=100，000，000时，精度达到了小数点后十二位。

3. 幂级数

幂级数的方法有很多，我们学习了使用arctan(x)的Maclaurin展开式计算π的近似值，π=4arctan1=4(1-13+15-…+-1n+12n-1+…

代码如下

n=10; %展开次数

s=0;

digits(50) %计算精度

for k=1:n

s=s+4*(-1)^(k+1)/(2*k-1);

end

vpa(s,20)

改变n的值，计算不同情况下得出的圆周率近似值

n

10

100

1000

10000

100000

近似值

3.0418

3.131592

3.14059

3.1414926535

3.1415826535897

该计算方法能比较准确地得出较高精度的圆周率的近似值

4. 韦达公式

根据韦达在1593年给出的公式 2π=22*2+22*…

代码如下：

a=sqrt(2);

s=1;

n=100

for i=1:n

s=(s*a)/2;

a1=sqrt(2+a);

a=a1;

end

vpa(2/s,20)

使用该算法计算圆周率的近似值时，n与所得的近似值如下：

n

10

100

1000

Pi

3.1415914

3.14159265358979385

3.1415926535897932385

事实上，当n为100时所算出的圆周率精度已经超过了小数点后100位，可见其计算精度之高，效率之快。一下列出n=100时所得结果精确到100位的输出值：

3.1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164 0628620899 8628034825 342117068

三、实验感想

圆周率的近似计算一直是一种挑战，经过这次研究，我们更加深刻地体会到了数学的魅力，同时也对matlab更加熟悉。在不断的学习中，我们成长了很多。

参考文献：

【1】 杨振华，郦志新《数学实验》北京：科学技术出版社，2010.2

【2】 赵洪牛《高等数学》北京：高等教育出版社

【3】 百度百科：圆周率

6 / 6

展开阅读全文