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 收敛更快。
若取12x =,则我们只有找出α与4π的关系,才能求出π的值。
令1
arctan 2α=,4πβα=-, 根据公式tan tan tan()1tan tan αβαβαβ++=
-有1tan 3β=,则有11arctan arctan 423π=+。
所以可以用1
1arctan arctan 4
23π=+来计算π的值。
程序('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.
显然,级数收敛越快,取同样的n 值可以得到更高的精度。
以同样的方法,能得出114arctan arctan 45239
π=+,程序和上面的一样。
这样π的近似值可以精确到几百位。
3、数值积分法
思想:
半径为1的圆的面积是π。
以圆心为原点建立直角坐标系,则
圆在第一象限的扇形是由y =与x 轴,y 轴所围成的图形,扇形的面积
4
s π
=。
只要求出扇形的面积,就可得出π的值。
而扇形面积可近似等于定积
分0⎰的值。
对于定积分()b
a f x dx ⎰的值,可以看做成曲线与x 轴,x a =,x
b =所
围的曲边梯形的面积S 。
把[],a b 分成n 等分,既得1n +个点x a =,1x ,⋅⋅⋅1n x -,x b =组成n 个小区间,每一个小区间与x 轴,x a =,x b =所围成的图形是一个小曲边梯形。
而梯形的面积计算公式是(+)2⨯÷上底下底高,对于第i 个小曲边梯形有上底为1()i f x -,下底为()i f x 。
所有小梯形的高都为()/h b a n =-。
所以第i 个小曲边梯形的面积为1(()+())2i i f x f x h -⨯÷。
曲边梯形的总面积即定积分()b
a f x dx ⎰的值就是
所有小梯形的面积总和。
为了避免根号,我们也可以利用积分
120114dx x π=+⎰ 得出π的值。
我们可以利用对求曲边梯形的面积来得出定积分120
11dx x
+⎰的值,从而得出π的值。
程序(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);
xj=a+(i+1)*h;
yj=1/(1+(xj)^2);
s=s+(yi+yj)*h/2;
end
p=vpa(4*s,30)
程序运行结果:
p =
3.
对于数值积分法求π值,以上程序简洁明了。
我们也可以以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.
用以上三种方法求π,n都取1000时,泰勒级数法求π,得到的近似值精度最高。