利用Matlab近似计算圆周率的若干方法

叶提芳

摘 要 本文分别采取幂级数展开式的方法、随机数的方法、数值积分的方法和公式法结合Matlab程序实现对圆周率的近似计算，分析实验结果，比较每种方法的近似程度的高低，实现了Matlab实验和数学理论的很好结合。

关键词 Matlab实验 圆周率 幂级数 随机数 数值积分

中图分类号：TP312 文献标识码：A

Some Methods of Approximatively Calculating Using Matlab Experiments

YE Tifang

(Industrial and Commercial College, Wuhan Polytechnic University, Wuhan, Hubei 430065)

Abstract In this paper, we used expanding power series, random number, numerical integration and formula methods combining matlab experiments to achieve the approximate value of. Then, we analyzed the experiments results, and compared the degree approximation of every method. It achieved the satisfying results of combination withMatlab experiment and mathematics theory.

Key words Matlab experiments; ; power series; random number; numerical integration

古今中外，历史上有许多人积极致力于圆周率的研究与计算。我国的刘徽用正3072边形得到小数点后的5位精度值,Ludolph Van Ceulen用262正边形得到了小数点后的35位精度值。这种方法虽然经典，但相当耗时。20世纪，很多数学家采取级数来近似计算圆周率的方法，已经能把圆周率近似计算精确到了上亿位，可以说，我们对圆周率的近似计算研究已经相当成熟。本文试在利用Matlab实验和高等数学中的知识有机结合起来，分别采取幂级数展开式的方法、随机数的方法、数值积分的方法和公式法结合Matlab程序实现对圆周率的近似计算。

1 利用幂级数展开式的计算方法

设是以2 为周期的周期函数，在[- , )上的表达式如下：

显然，为奇函数，利用我们在高等数学幂级数一章知识①，可以将其展开为正弦级数为：= 4

当我们取 = 时，得到一数项级数4，因为，则有4 =，求其部分和为 = 4，显然有 =，下面我们分别给出 = 1000,10000,50000时，利用Matlab程序求得的 的近似值。程序如下：

>> s=0;

>> n=50000;

>> digits(22)

>> for k=1:n

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

end

>> s=vpa(4*s,20)

s =

3.1415726535897814387

= 1000,10000,50000时， 的近似值如下：

2 采取随机数的近似计算方法

设一制作均匀的冰激凌可以看做是由圆锥面 = 和球面 ++ （）2= 1围成②。我们利用积分知识求得它的体积为：

其中。

我们还可以采取随机数的方法，由于所求锥形体可表示为：= {（）| + ≤ （）,+ ≤}，

它位于长方体 = {（）|-1≤≤1}内部，该长方体体积为8，往长方体内部随机投点个，然后统计锥形体内的随机点数，则≈ ≈8祝因为计算结果带有随机性，我们用十次重复的计算作比较，下面为N取100000时，用Matlab近似 的程序：

>> for k=1:10

r=rand(100000,3);

x=2*r(:,1)-1;

y=2*r(:,2)-1;

z=2*r(:,3);

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

p(k)=8*sum(fl-z.*z<=0&fl-z.*(2-z)<=0)/100000;

end

>> p

运行程序，得到 的近似值如下：

3.11363.12003.10003.16963.15603.12483.14643.09203.20003.1752

当取5000000时，运行程序，得到的 的近似值如下，可见近似程度是较好的：

3.14293.14293.14373.14053.14183.13973.13983.14093.14323.1430

从计算结果看：这种方法虽然简单可行，但收敛的速度慢，距离真实值误差较大。

3 利用定积分 = 近似计算

根据定积分的定义，（积分的结果和区间的分法及的取法无关，现在采取特殊分割和特殊取法不影响结果）将区间[0,1]分成等份，在每个小区间上，选取中点为，有

程度如下：

>> n=1000;

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

运行程序，得到的 的近似值为3.1415927369231306798，可见，近似程度已经很可观了。

4 其它方法

我们知道 = ,= ,…= （重根号），即 = ，由此我们得到韦达公式 = …，据此= 2/（…），

我们编写Matlab程序如下（下面为计算式分母中取前10项乘积的程序）:

a=sqrt(2);

>> s=1;

>> for i=1:10

s=(s*a)/2;

a1=sqrt(2+a);

a=a1;

end

>> vpa(2/s,20)

运行程序，得到结果：3.1415914215111997443。当计算式分母中取前100项乘积时，运行程序，得到结果：3.1415926535897932385，可见，这种方法收敛速度快，近似精度高。

注释

① 高纯一,周勇.高等数学[M].上海:复旦大学出版社,2006:249-259.

② 李继成.数学实验[M].北京:高等教育出版社,2006:92.