【MATLAB】蒙特卡洛估计pi值

x =-1+ 2*rand(1e4,1);
y =-1+ 2*rand(1e4,1);

%搜寻落到半径为1的圆内的个数
z = x((x.^2 + y.^2 )<=1);%其实是先判断，后索引符合条件的x值

pitest = 4*length(z)/length(x); %原理: 所有随机数落到半径为1的圆内的概率等
err = abs(pitest - pi)/pi;
fprintf( '2.梦塔卡罗模拟π的值为%3.4f,相对误差为%3.4f\n\n ' ,pitest,err)
figure(1)
hold on
axis equal
set(gca, 'XLim',[-1 1]);
set(gca, 'YLim',[-1 1]);
for k = 1:1000 %length(x)
    if x(k)^2 + y(k)^2 <= 1
        %the point is inside the sphere, show it in red
        plot(x(k),y(k), '*b')
    else
        %the point is outside the sphere, show it in blue .
        plot(x(k), y(k),'.g')
    end
    drawnow update
end