Matlab (Монте-Карло)

clear all, clc, close all
%%%%%%% Расчет площади круга методом Монте-Карло
r = 5;
x0 = 1;
y0 = 2;
%%% Построение окружности
t = 0 : r/500 : 2*pi;
x = r*cos(t) + x0;
y = r*sin(t) + y0;
line(x, y, 'linew', 2, 'color','r')
%%% Центр круга
line(x0,y0,'marker','o','markerfacecolor','r','color', 'r' )
%%% Построение квадрата, в который вписан круг
line([x0-r, x0+r],[y0+r, y0+r], 'lines','-.')
line([x0-r, x0+r],[y0-r, y0-r], 'lines','-.')
line([x0+r, x0+r],[y0-r, y0+r], 'lines','-.')
line([x0-r, x0-r],[y0-r, y0+r], 'lines','-.')
%%% Генерирование случайных чисел для области D*
N = 500; %%% число испытаний
%%% Генерация чисел по горизонтальной стороне квадрата
rx = min(x) + (max(x) - min(x))*rand(N,1);
%%% Генерация чисел по вертикальной стороне квадрата
ry = min(y) + (max(y) - min(y))*rand(N,1);
%%%% Заполнение случайными числами квадрата с кругом
for J = 1 : N
line(rx(J),ry(J),'marker','o','markersize',2,'color', 'k','markerfacecolor', 'k' )
end
%%% Подсчет количества случайных чисел, попавших в круг
m = 0;
for J = 1 : N
if (rx(J) - x0)^2 + (ry(J) - y0)^2 <= r^2
m = m + 1;
end
end
%%%%% Расчет площади квадрата
S = ((x0+r) - (x0-r))^2;
%%%%% Расчет числа ПИ
pi1 = m*4/N
%%% Расчет площади круга методом Монте-Карло
S0 = m/N*S
%%% Проверка расчета
Scontrol = pi*r^2
%%% Заголовок для диаграммы
str = sprintf('%s Радиус круга r = %g, координаты центра x_0 = %g, y_0 = %g.%s%g%s%g. %s%g%s%g', '\bf',r,x0,y0, ...
'\newline Теоретическая площадь круга S = ', Scontrol, '\newline Значение числа pi = ', pi1, 'Площадь круга по методу Монте-Карло S0 = ', S0, ...
'\newlineЧисло испытаний N = ', N);
title(str)
grid on
xlabel('\bf - - - - - - - X - - - - - - - ')
ylabel('\bf - - - - - - - Y - - - - - - - ')
%%% Ограничения по осям
xlim([x0 - r - r/10, x0 + r + r/10])
ylim([y0 - r - r/10, y0 + r + r/10])
axis equal