如何计算π的值(MATLAB)

求解圆周率π的几种方法,MATLAB代码

如何计算 的值

1、蒙特卡罗（Monte Carlo）法

思想：

取一正方形A，以A的一个顶点为圆心，A的边长为半径画圆，取四分之一圆（正方形内的四分之一圆）为扇形B。已知A的面积，只要求出B的面积与A 的面积之比k

SB

，就能得出SB，再由B的面SA

积为圆面积的四分之一，利用公式S圆= R2即可求出 的值。因此，我们的目的就是要找出k的值。

可以把A和B看成是由无限多个点组成，而B内的所有点都在A内。随机产生n个点，若落在B内的有m个点（假定A的边长为1，以扇形圆心为坐标系原点。则只要使随机产生横纵坐标x、y满足，则可近似得出k的值，即k x2 y2 1的点，就是落在B内的点）由此就可以求出 的值。 程序（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)

m

，n

求解圆周率π的几种方法,MATLAB代码

程序运行结果： p =

3.14000000000000000000000000000

2、泰勒级数法

思想：

反正切函数arctanx的泰勒级数展开式为：

2k 1

x3x5k 1xarctanx x ( 1) 352k 1

将x 1代入上式有

111

arctan1 1 ( 1)n 1. 4352n 1

利用这个式子就可以求出 的值了。 程序（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时，就会花费很长时间，而且精度也不是很高。原因是x 1时，arctan1的展开式收敛太慢。因此就需要找出一个x使得arctanx收敛更快。

求解圆周率π的几种方法,MATLAB代码

若取x ，则我们只有找出 与的关系，才能求出 的值。 令 arctan，

12

12 4

4

，

根据公式tan( )

tan tan 1 11

有tan ，则有 arctan arctan。

34231 tan tan

1

2

13

所以可以用 arctan arctan来计算 的值。

4

程序（2）：

i=1;n=1000;s=0;s1=0;s2=0; for i=1:n

s1=s1+(-1)^(i-1)*(1/2)^(2*i-1)/(2*i-1); s2=s2+(-1)^(i-1)*(1/3)^(2*i-1)/(2*i-1); end s=s1+s2; p=vpa(4*s,30) 程序运行结果： p =

3.14159265358979323846264338328

显然，级数收敛越快，取同样的n值可以得到更高的精度。以同样的方法，能得出 4arctan arctan

4

'

1

51

，程序和上面的一样。这样 239

的近似值可以精确到几百位。

3、数值积分法

思想：

半径为1的圆的面积是 。以圆心为原点建立直角坐标系，则

求解圆周率π的几种方法,MATLAB代码

圆在第一象限的扇形是由y 与x轴，y轴所围成的图形，扇形的面积

s

4

。只要求出扇形的面积，就可得出

的值。而扇形面积可近似等于定积

分 的值。

对于定积分 af(x)dx的值，可以看做成曲线与x轴，x a，x b所围的曲边梯形的面积S。把 a,b 分成n等分，既得n 1个点x a，x1，

xn 1，x b组成n个小区间，每一个小区间与x轴，x a，x b所围

b

成的图形是一个小曲边梯形。而梯形的面积计算公式是

(上底+下底) 高 2，对于第i个小曲边梯形有上底为f(xi 1)，下底为f(xi)。所有小梯形的高都为h (b a)/n。所以第i个小曲边梯形的面

积为(f(xi 1)+f(xi)) h 2。曲边梯形的总面积即定积分 af(x)dx的值就是所有小梯形的面积总和。

为了避免根号，我们也可以利用积分

1

01 x2

4

1

b

得出 的值。

我们可以利用对求曲边梯形的面积来得出定积分 0

从而得出 的值。 程序（3）：

a=0;b=1;s=0;n=1000;i=0; h=(b-a)/n; for i=0:(n-1) xi=a+i*h; yi=1/(1+(xi)^2);

1

1

的值，2

1 x

求解圆周率π的几种方法,MATLAB代码

xj=a+(i+1)*h; yj=1/(1+(xj)^2); s=s+(yi+yj)*h/2; end

p=vpa(4*s,30) 程序运行结果： p =

3.14159248692312775830259852228

对于数值积分法求 值，以上程序简洁明了。我们也可以以x做循环，用一条语句求出 值。 程序（3’）： s=0;n=1000; for x=0:(1/n):(1-(1/n))

s=s+(1/(1+x^2)+1/(1+(x+(1/n))^2))*(1/n)/2; end

p=vpa(4*s,30) 程序运行结果： p =

3.14159248692312775830259852228

用以上三种方法求 ，n都取1000时，泰勒级数法求 ，得到的近似值精度最高。

求解圆周率π的几种方法,MATLAB代码