蒙特卡洛方法求圆周率

基于Julia语言的蒙特卡洛方法求圆周率 $\pi$

原理：$(x,y)$是平面上的一点，其中 $x,y$ 服从均匀分布 $U(-1,1)$，假设得到2000个样本点，如下图所示

圆的面积 $S_1=\pi r^2=\pi$，正方形的面积 $S_2=(2r{)}^2=4r^2=4$。若落在圆内的点为$N_1$个，落在正方形内的样本点数是$N_2$个，$N_2$ 也是所有的样本点数。因此可以得到下面的近似关系
$$\frac{S_2}{S_1}=\frac{4}{\pi }\approx \frac{N_2}{N_1}$$
即 $S_1=\pi \approx \frac{4N_1}{N_2}$。

用Julia语言实现的代码如下所示：

#导入包
using Plots
using Distributions
#用蒙特卡洛方法求pi值
sample_num = 2000;
dist = Uniform(-1,1);

sample_x = rand(dist,sample_num);
sample_y = rand(dist,sample_num);

count = 0;
plot()

for k in range(1,sample_num,step=1)
    if (sample_x[k]^2+sample_y[k]^2)<=1
    #落入圆内的点
        count = count + 1
        plot!([sample_x[k]],[sample_y[k]],shape=:circle,markersize=2,
            markercolor=:black,mscolour=:black)
    else
    #落入圆外的点
        plot!([sample_x[k]],[sample_y[k]],shape=:circle,markersize=2,
            markercolor=:green,mscolour=:green)
    end
end

print("sample_num = ",sample_num,"\n","count = ",count,"\n","pi=",count/sample_num * 4.0)


#画出圆心为（x0,y0）,半径为r的圆
x0 = 0
y0 = 0
r = 1

tmin = 0
tmax = 2pi

tvec = range(tmin,tmax,length=100)

x = r.*sin.(tvec).+ x0
y = r.*cos.(tvec).+ y0

plot!(x,y,aspect_ratio = 1,linecolor=:red,legend=false)

# 画圆心
plot!([x0],[y0],shape=:circle,markercolor=:black,mscolour=:black)

#画出矩形 (x_left_up,y_left_up)、(x_right_up,y_right_up)、
#       (x_right_down,y_right_down)、(x_left_down,y_left_down)

x_left_up = -1; x_right_up = 1; x_right_down = 1; x_left_down = -1;
y_left_up = 1; y_right_up = 1; y_right_down = -1; y_left_down = -1;

rectangle_x = [x_left_up, x_right_up, x_right_down, x_left_down, x_left_up]
rectangle_y = [y_left_up, y_right_up, y_right_down, y_left_down, y_left_up]
plot!(rectangle_x,rectangle_y,linecolor=:blue,legend=false,linewidth=2)

#保存图片
savefig("montecarluo")

仿真结果：