数学建模——π的计算

实验报告

实验目的

1. 掌握用MATLAB软件求π值的方法，并对结果作初步分析；

2.掌握用MATLAB软件实现蒙特卡洛法的过程。

实验内容

运用数值积分法、泰勒级数法、模拟蒲丰实验和随机整数互素法计算π值，取n =（k=3, 4, …,7），记录所得结果和程序

运行时间，并作简要分析。

实验过程

一、数值积分法

单位圆的面积等于以单位圆的圆心为原点建立直角坐标系，单位圆在第一象限内的部分G 是一个扇形，其面积S = π/4

只需计算出S 的近似值，它的4倍就是π的近似值

1、梯形公式法

将扇形G 分为n 个宽度相等的部分Gk（1 ≤k ≤n）

Gk：曲边梯形左右边界平行，上方边界为曲线n→∞Gk→梯

形

Gk 梯形面积：

等价于

π

取 n =10k （k=3, 4, …,7），按上述方法，通过计算扇形面积计算 π的近似值。

程序

Clear

tic

n=1e 3(3, 4, …,7); x=linspace(0,1,n+1); y=(1-x.^2).^0.5;

pai=(0.5*y(1)+0.5*y(n+1)+sum(y(2:n)))*4/n toc

G G

2、辛普森（Simpson）公式法

仍用分点xi = a + i(b-a)/n (1≤ i ≤ n-1) 将区间[a, b]分成n 等份，直线x = xi (1≤ i ≤ n-1) 将曲边梯形分成n 个小曲边梯形，再作每个小区间[xi-1, xi]的中点

将第i个小曲边梯形的上边界y=f(x) (xi-1≤x≤xi) 近似地看作经

过这三点的抛物线段，则可求得：

其中

于是得到S

;

即

程序：

clear

tic

n=1e3(3,4,5,6,7);

x=linspace(0,1,n+1);

x2=x(1:n)+1/(2*n);

y=4./(1+x.^2);

y2=4./(1+x2.^2);

pai=((y(1)+y(n+1))+2*sum(y(2:n))+4*sum(y2))/(6*n) toc

二、泰勒级数法

考虑反正切函数的泰勒级数

取x =1，n =1e3(3,4,5,6,7)

程序：

clear

tic

k=1e3(3,4,5,6,7);

x=1;

n=1:k;

pai=sum((-1).^(n-1).*x.^(2.*n-1)./(2.*n-1))

toc

结果显示：花费时间长！准确度差！！

原因：x =1时得到的arctan1的展开式收敛太慢！

Maqin公式法

提高收敛速度：

x << 1随着指数的增加，x的幂快速接近于0，泰勒级数就会快速收敛

取2k-1 = 63，其误差已经非常小，表明收敛速度非常快

有

再考虑收敛更快的：

/4=arctan 1/2+arctan 1/3程序：

clear

tic

k=1e3(3,4,5,6,7);

n=1:k;

x=1/2;

at_1=sum((-1).^(n-1).*x.^(2.*n-1)./(2.*n-1));

x=1/3;

at_2=sum((-1).^(n-1).*x.^(2.*n-1)./(2.*n-1));

pai_1=4*at_1+4*at_2

toc

/4=4arctan 1/5-arctan 1/239程序：

clear

tic

k=1e3(3,4,5,6,7);

n=1:k;

x=1/5;

at_5=sum((-1).^(n-1).*x.^(2.*n-1)./(2.*n-1));

x=1/239;

at_239=sum((-1).^(n-1).*x.^(2.*n-1)./(2.*n-1));

pai=16*at_5-4*at_239

toc

三、蒙特卡罗（Monte Carlo）法

在数值积分中，我们利用求1/4单位圆的面积来得到π/4，进而得到π 1/4单位圆是一个扇形G，它是边长为1的单位正方形G1

的一部分。G1的面积S1=1，只要能求出扇形

G 的面积S在正方形G1的面积S1中所占

的比例即可求得π值。4

问题转化为求扇形面积在正方形面积中所占比例k

在正方形中随机地投很多点，使所投的每

个点落在正方形中每一个位置的机会均

等，看其中有多少个点落在扇形内。记投

点总数为n，落在扇形内的点的个数为m，

则k

在计算机上实现随机投点，需要能产生[0, 1]上均匀分布的随机数。在Matlab中，使用rand函数实现这一功能。

rand 产生1个随机数

rand(n) 产生n×n的随机方阵

rand(m,n) 产生m×n的随机矩阵

需要产生两个[0, 1]上均匀分布的随机数x，y，表示一个点的坐标，这个点落在单位正方形内每个位置的机会均等，当x2 + y2 ≤ 1 时，该点落在扇形内。

程序：

clear

tic

n=1e7;

x=rand(1,n);

y=rand(1,n);

c=x.^2+y.^2;

m=sum(c<=1);

k=m/n;

pai=4*k

toc

四、蒲丰（Buffon）掷针试验

1777年法国数学家蒲丰提出了一种用于计算π的随机掷针试验，步骤如下：

（1）取一张白纸，在上面画许多间距为d 的等距平行线；（2）取一根长度为l ( l < d ) 的均匀直针，随机地向画有平行线的纸

上掷去，一共投掷n次（n是一个很大的整数），观察针和直线相交的次数m；

（3）针和直线相交的概率p = 2l / (πd )，取m/n为p的近似值，

则π

针和直线相交的概率为p=2l/(d) 的证明过程设

M 表示针的中点，x 表示针落下后M 到最近一条平

行线的距离，表示针与直线的夹角。显然有0≤x≤d/2

Ox平面上表示一个矩形Ω。针和某一条平

行线相交的充要条件是：Ox

平面上表示一个区域A，点M落在区域A内的概率为

程序：

clear

tic

n=1e7;

d=2;

l=1;

x=abs(d*rand(1,n)-1); y=0.5*pi*rand(1,n);

m=sum(x<=0.5*sin(y));

pai=n/m

toc

五、随机整数互素法

蒙特卡洛法随机整数互素的概率：

取一个大的整数N，在1到N 之间随机地取一对整数a，b，找出它们的最大公约数(a, b)，当(a, b)=1时称a，b互素。做n 次这样的实验，记录其中(a, b)=1的情况出现的次数m，算出p = m/n的值。注：随机整数互素的概率：

程序：

Clear

tic

n=1e3(3,4,5,6,7);

N=1e7;

a=round(1+N*rand(1,n));

b=round(1+N*rand(1,n));

m=sum(gcd(a,b)==1);

p=m/n

pai=(6*n/m)^0.5

toc

实验结论

实验结果统计表

π=3.141592653589793

3.141555466911023 3.141591477611321 3.14159261640199 3.141592653552502 3.1415926535525020.000345 s 0.001115 s 0.008649 s 0.113040 s 0.925687 s 3.141592653589790 3.141592653589802 3.141592653589803 3.141592653589796 3.1415926535898560.000345 s

0.001375 s 0.010334 s 0.128711 s 0.931969 s

3.140592653839794 3.141492653590035 3.141582653589780 3.141591653589692 3.1415925535897800.000409 s 0.001900 s 0.012589 s 0.113065 s 1.236871 s 3.141592653589792 3.141592653589792 3.141592653589792 3.141592653589792 3.1415926535897920.001705 s 0.008424 s 0.075838 s 0.669452 s 5.293582 s 3.141592653589794 3.141592653589794 3.141592653589794 3.141592653589794 3.1415926535897940.002017 s 0.009191 s 0.078560 s 0.694382 s 5.282906 s

3.152000000000000 3.142000000000000 3.149040000000000 3.141124000000000 3.1422036000000000.000293 s 0.001201 s 0.011782 s

0.110932 s

1.068707 s 3.174603174603174 3.089280197713933 3.123048094940662 3.135533432625226 3.1420130186167410.000337 s

0.001807 s

0.016089 s

0.150625 s

1.423036 s

3.216337604513385 3.144508978712735 3.144042691862183 3.140239********* 3.1411480710860620.039859 s

0.367654 s

3.461486 s

33.917678 s

239.278874 s

蒙特卡罗（M ont e C ar l o）法

模拟蒲松实验法随机整数互素法

方法 n 数值积分法

泰

勒级数法

4arct an 1/5-arct an 1/239

arct an 1/2+arct an 1/3

梯形公式

辛普森公式普通算法

Maq in 公式

结果分析：

由以上计算结果知：利用概率统计原理的三种方法(模拟蒲松实验法、蒙特卡罗（Monte Carlo）法、随机整数互素法)的计算结果耗时长，运算次数多，精确度差，但该原理适用性强；而数值积分法和普通的泰勒级数法的耗时短，运算次数较小，精确度相对较高；泰勒级数法中的Maqin公式法耗时很短，运算次数小，精确度高。

注：泰勒级数法中的

可将算到200位。