montecarlo_pi.m

来自「monte carlo方法在各个方面中的小应用」· M 代码 · 共 23 行

%求pi的程序
n=input('n=?');%总次数
m=input('m=?');%步长
k=1;
j=1;
h=0;
for i=1:n
    x=rand();y=rand();
    rxy=x.^2+y.^2;
    if rxy<1
        h=h+1;
    end
    if(i==k*m)%记录单位步长的pi值
        aps(j)=m*k;
        k=k+1;
        pie(j)=(4./aps(j))*h;
        j=j+1;    
    end
end
plot(aps,pie,'r')
pie(j-1)*(1-pie(j-1)/4)%最后一个pi值的方差
4*sqrt(pie(j-1)/4*(1-pie(j-1)/4)/n)%最后一个pi值的标准差