现在的位置: 首页 > 综合 > 正文

采用蒙特卡洛方法计算PI的值

2014年02月20日 ⁄ 综合 ⁄ 共 1151字 ⁄ 字号 评论关闭

蒙地卡罗方法(Mente Carlo method)是一种通过产生大量随机数并结合相关统计和数值计算来解决问题的一类方法。本文介绍如何利用蒙地卡罗方法来计算pi的值。


基本思想: 利用圆与其外接正方形面积之比为pi/4的关系,通过产生大量均匀分布的二维点,计算落在单位圆和单位正方形的数量之比再乘以4便得到pi的近似值。样本点越多,计算出的数据将会越接近真识的pi(前提时样本是“真正的”随机分布)。


在Matlab中,只需以下几行代码:

N = 10000;%样本点
rng('shuffle');%初始化随机发生器
x = rand(N,1);%二维样本的x坐标,x: [-1,1]
rng('shuffle');%再次初始化随机发生器,与上次不同,将产生独立的随机数
y = rand(N,1);%二维样本的y坐标, y:[-1:1]
s = sum(x.^2+y.^2 <= 1); %计算落在单位圆内的点数
my_pi = s/N*4; %计算pi值

在我的电脑上其中一次计算出的pi值为3.1672, 与真正的pi还有一段矩离。原因是 1) 样本点不够大, 和 2) matlab产生数是伪随机的,而不是真的随机。原因2暂时不在此深究,有兴趣的同学可以点此参考. 关于原因1, 可以通过下面的程序看得更清楚:

%%
clc;
clear;
close all;
%%
sample_size = 10.^(2:1:7);
my_pi = zeros(length(sample_size),1);
for i = 1:length(sample_size)
    N = sample_size(i);
    rng('shuffle');
    x = rand(N,1);
   rng('shuffle');
    y = rand(N,1);
    s = sum(x.^2+y.^2 <= 1);
    my_pi(i) = s/N*4;
end
%%
plot(abs(my_pi-pi),'g-s','LineWidth',1.5);
set(gca,'XScale','log');
xlabel('number of points');
ylabel('abs(my\_pi-pi)');
grid on;


以上的程序通过增加样本数(从100到10,000,000)来计算pi的值并与真实值进行比较,最后得到下图:



由图可见,样本数越大,产生的pi值越接近真实值。


小结:

本文简单演示了如何通过蒙地卡罗的方法计算pi的值,旨在介绍蒙地卡罗的思想。关于pi的值如何计算,可以参考http://en.wikipedia.org/wiki/Pi。 事实上,没有人会用这样的方法来计算pi,因为机数生成算法(random
number generator, RNG)实在是“太不随机了”,在我的电脑上进行了上亿次采样也只能到3.141592左右变动,再也无法更精确了。



抱歉!评论已关闭.