Метод Монте-Карло (ММ РУиС)

Навеяно введением в курс Statistical Mechanics Algorithms and Computations with Werner Krauth

% S of square: d^2
% S of cicle: pi*r^2 = pi/4 * d^2
% Ss / Sc = 4/pi; pi = Sc / Ss * 4 = Nc / Ns * 4;

N = 10000000;
M = 1000000;
K = N/M;

Nc = 0;
rng('shuffle');

doplot = 0;

for k = 1:K
    x = 2*rand(M,1) - 1;
    y = 2*rand(M,1) - 1;
    inc = find(x.*x + y.*y <= 1);
    Nc = Nc + length(inc);
    newpi = Nc / (M*k) * 4;
   
    if doplot
        plot(x, y, '*b');
        hold on
        plot(x(inc), y(inc), '*r');
        drawnow
    end

    fprintf('pi = %.6f\n', newpi);
end

См. также

• http://habrahabr.ru/post/128454/