Моделирование систем
[+]
Мордовский государственный университет имени Н.П. Огарева
Опубликован: 30.11.2010 | Доступ: свободный | Студентов: 3376 / 2068 | Оценка: 4.12 / 4.13 | Длительность: 14:37:00
ISBN: 978-5-9963-0352-6
Темы: Программирование, Алгоритмы и дискретные структуры
Специальности: Программист, Системный архитектор
Лекция 4:

Выборочный метод Монте-Карло
A
|
версия для печати
< Лекция 3 || Лекция 4: 12 || Лекция 5 >

Практическая часть

Пример 1. Рассчитайте с помощью метода Монте-Карло площадь круга радиусом r = 5 и с координатами центра (1; 2). Число испытаний примите 500. Выполните графические построения, поясняющие расчет.

Для определения площади плоской фигуры ее вписывают в соответствующую известную фигуру, площадь которой достаточно просто вычисляется. Вычисление площади круга может быть произведено через вычисление площади квадрата, в который вписывается круг. Выборочный метод Монте-Карло предполагает генерирование случайных (псевдослучайных) равномерно распределенных чисел. Этим числам сопоставляются координаты точек для рассматриваемой фигуры — квадрата, в которую вписывается круг. Если площадь квадрата S подсчитана, то задается число испытаний N, из которых m исходов могут оказаться внутри круга или на его границе, т. е. на окружности. Тогда площадь круга S0 будет определяться выражением

S0=\frac{m}{N}S,

где:

N — число сгенерированных случайных чисел, соответствующих количеству точек, находящихся в данном квадрате:

m — число случайных чисел (точек), которые попали в круг. Граница круга — это окружность, уравнение которой имеет вид

(x-x0)^2+(y-y0)^2=r^2,

где:

х0, y0 — координаты центра окружности;

х, y — текущие значения переменных;

r — радиус окружности.

Программный код решения примера:

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;
 
%%% Расчет площади круга методом Монте-Карло
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', '\bf',r,x0,y0, ... 
'\newline Теоретическая площадь круга S = ', Scontrol, ...
'\newline Площадь круга по методу Монте-Карло S0 = ', S0, ...
 'Число испытаний 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

Задание 1

  1. Видоизмените программу так, чтобы выполнялось условие непревышения относительной погрешности заданной величины (например, 1%) и чтобы производился необходимый подсчет числа испытаний.
  2. Предусмотрите, чтобы точки, попавшие в круг, были красного цвета, а не попавшие в него – синего цвета.
  3. Напишите программу расчета числа \pi по методу Монте-Карло. Расчет числа \pi произведите с точностью в четыре знака после десятичной точки.

Пример 2. Рассчитайте с помощью метода Монте-Карло площадь фигуры, ограниченной эллипсом с большой осью a = 6 и малой осью b = 4. Координаты пересечения осей эллипса x0 = –2, y0 = 3, оси коллинеарны декартовым осям.

Теоретическая площадь S эллипса вычисляется по формуле

S=\pi ab,

где:

a — большая ось эллипса;

b — малая ось эллипса.

Каноническое уравнение эллипса (относительно центра координат):

\frac{x^2}{a^2}+\frac{y^2}{b^2}=1.

Программный код решения примера:

clear all, clc, close all
%%%% Расчет площади эллипса методом Монте-Карло
%%% Параметры эллипса
a = 6;
b = 4;
x0 = -2;
y0 = 3;

%%% Построение эллипса
t = 0 : min([a, b])/500 : 2*pi;
x = a*cos(t) + x0;
y = b*sin(t) + y0;
line(x,y, 'linew', 2, 'color', [1,0,0])
%%% Центр пересечения осей эллипса
line(x0,y0,'marker','o','markerfacecolor','r','color','r' )
 
%%% Построение прямоугольника, в который вписан эллипс
line([x0-a, x0+a],[y0+b,y0+b], 'lines','-.')
line([x0-a, x0+a],[y0-b,y0-b], 'lines','-.')
line([x0-a, x0-a],[y0-b,y0+b], 'lines','-.')
line([x0+a, x0+a],[y0-b,y0+b], 'lines','-.')
 
%%% Гененрирование случайных чисел для области D* - области прямоугольника
N = 500; %%% число испытаний
%%% Генерация чисел по горизонтальной стороне прямоугольника
rx = (x0-a) + (x0+a - (x0-a))*rand(N,1);
%%% Генерация чисел по вертикальной стороне прямоугольника
ry = (y0-b) + (y0+b - (y0-b))*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)/a)^2 + ((ry(J) - y0)/b)^2 <= 1
        m = m + 1;
    end
end
 
%%%%% Расчет площади прямоугольника
S = ((x0+a) - (x0-a))*((y0+b)-(y0-b));
 
%%% Расчет площади эллипса методом Монте-Карло
S0 = m/N*S
 
%%% Проверка расчета 
Scontrol = pi*a*b
 
%%% Заголовок для диаграммы
str = sprintf('%s Оси эллипса a = %g, b = %g, координаты центра x_0 = %g,  y_0 = %g.%s%g%s%g. %s%g', '\bf',a, b,x0,y0, ... 
'\newline Теоретическая площадь эллипса S = ', Scontrol, ...
'\newline Площадь эллипса по методу Монте-Карло S0 = ', S0, ...
 'Число испытаний N = ', N);
title(str)
grid off %% выключение сетки на диаграмме
xlabel('\bf - - - - - - - X  - - - - - - - ')
ylabel('\bf - - - - - - - Y  - - - - - - - ')
%%% Установление размеров и свойств графического окна
gfig = get(0, 'screensize');
set(gcf, 'color', 'w', 'position', [gfig(1) + 100, gfig(2) + 100, gfig(3)*0.8, gfig(4)*0.7])
%%% Ограничения по осям
xlim([x0-a-a/10, x0+a+a/10])
ylim([y0-b-b/10, y0+b+b/10])
axis equal

Задание 2

  1. Для построения эллипса используйте его каноническое уравнение.
  2. Предусмотрите, чтобы точки, попавшие в эллипс, были красного цвета, а не попавшие в него — синего цвета.
  3. Координаты центра пересечения осей эллипса примите случайными при использовании функции randn.
  4. Оси эллипса примите случайными, равномерно распределенными из интервала [Х; 5Х], где Х — номер компьютера, за которым выполняется лабораторная работа ( 1, 2, 3, ...).
  5. В программу включите расчет общего числа испытаний, при котором достигается заданная относительная погрешность в 2%.
  6. В программу включите построение графика относительной погрешности от числа испытаний начиная с N = 100 и более.

Пример 3. Найдите оценку значения тройного интеграла методом Монте-Карло

I=\mathop{{\int\int ... \int}}\limits_{V}zdxdydz,

где V — область, ограниченная круговым конусом z^2=\frac{h^2}{R^2}(x^2+y^2) и плоскостью z = h, h = 5, R = 2.

Программный код решения примера:

clear all, clc, close all
%%% Вычисление тройного интеграла методом Монте-Карло, 
R = 2;
h = 5;

%%% Для построения конуса
[x,y] = meshgrid(-1.2*R:0.02:1.2*R, -1.2*R:0.02:1.2*R);
z = h/R*sqrt(x.^2 + y.^2);
xmax = 1.2*R;
%% Радиус окружности в сечении конуса на высоте max(z)
rmax = sqrt(xmax^2 + xmax^2);
%%% Максимальное значение аппликаты
zmax = h/R*sqrt(xmax^2 + xmax^2);
%%% Тангенс угла фронтального сечения конуса
koef = (rmax/zmax);
%%% Радиус окружности в сечении конуса на высоте h
rh = h*(koef); 
 
%%% Круговой конус
mesh(x,y,z) %% функция построения пространственных фигур
hold on
%%% Плоскость на высоте h
line(x, y, h*ones(size(x)))
 
%%% Параллелепипед 
line([-rh, rh],[-rh,-rh],[0,0])
line([-rh, -rh],[-rh,-rh],[0,h])
line([rh, rh],[-rh,-rh],[0,h])
 
line([-rh, rh],[rh,rh],[0,0])
line([-rh, -rh],[rh,rh],[0,h])
line([rh, rh],[rh,rh],[0,h])
 
line([-rh, -rh],[-rh,rh],[0,0])
line([rh, rh],[-rh,rh],[0,0])
 
%%% Вычисление тройного интеграла по встроенным функциям int
syms x y z
In0 = int(int(int(z, z, h/R*sqrt(x^2+y^2), h), ...
    y, -sqrt(rh^2 - x^2), sqrt(rh^2 - x^2)), x, -rh, rh);
In = double(In0)
%%% Расчет тройного интеграла методом Монте-Карло
N = 1000; %%% Число испытаний
S = 0;
for J = 1 : N
    xr = -rh + (rh - (-rh))*rand;
    yr = -rh + (rh - (-rh))*rand ;
    zr = h*rand;
if  zr^2<=h^2/R^2*(xr^2 + yr^2) & h/R*sqrt(xr^2 + yr^2)>=0 & ... 
         h/R*sqrt(xr^2 + yr^2) <= h 
 
        S = S + zr; 
      end
end
 
mD = [2*rh]*[2*rh]*[h]; %%% Объем параллелепипеда, мера
In_Car = (mD/N)*S
 
grid on
str = sprintf('%s%g;%s Точное значение тройного интеграла: %g',...
    '\bf Значение интеграла по методу Монте-Карло: ', In_Car, '\newline', In);
title(str)
 
view(-30,20)
axis equal

Задание 3

  1. Формирование случайных переменных xr, yr, zr вынесите за пределы цикла.
  2. Объясните расстановку пределов интегрирования.
  3. Вычислите по методу Монте-Карло объем тела для полученных пределов интегрирования. Сравните со значением, вычисленным аналитически.

Задание 4

В нижеприводимых заданиях выполните графические построения заданных плоских фигур и напишите программы по расчету площадей этих фигур в системе MATLAB.

  1. Рассчитайте с помощью метода Монте-Карло площадь фигуры, ограниченной гипоциклоидой с параметрами R = 3 ; r = 2.7 ; m = r/R. Параметрическое уравнение гипоциклоиды:

    x =(R-m*R)*cos(m*t)+m*R*cos(t-m*t),\\ y = (R-m*R)*sin(m*t)-m*R*sin(t-m*t).

    Подберите массив значений параметра t.

    Относительная погрешность расчета не должна превышать 5%.

  2. Рассчитайте с помощью метода Монте-Карло площадь фигуры, ограниченной кривой Штейнера с параметром r = 5. Параметрическое уравнение кривой Штейнера:

    x = 2*r*cos(t/3) + r*cos(2*t/3);\\ y = 2*r*sin(t/3) - r*sin(2*t/3);

    Подберите массив значений параметра t.

    Относительная погрешность расчета не должна превышать 5%.

  3. Рассчитайте с помощью метода Монте-Карло площадь фигуры, ограниченной астроидой с параметром R = 4. Параметрическое уравнение астроиды:

    x = R*(cos(t/4)).^{\wedge} 3;\\ y = R*(sin(t/4)).^{\wedge} 3;

    Подберите массив значений параметра t.

    Относительная погрешность расчета не должна превышать 5%.

Задание 5

В нижеприводимых заданиях выполните графические построения заданных пространственных тел и напишите программы по вычислению объемов этих тел в системе MATLAB.

  1. Вычислите с помощью метода Монте-Карло объем тела, ограниченного параболоидом \frac{y^2}{4}+\frac{z^2}{9}=2\frac{x}{3} и плоскостью х = 2.

    Относительная погрешность расчета не должна превышать 5%.

  2. Вычислите с помощью метода Монте-Карло объем тела, ограниченного поверхностями:

    y^2=16-6x,\qquad y^2=2x,\qquad z=\pm 5

    Относительная погрешность расчета не должна превышать 5%.

  3. С помощью метода Монте-Карло вычислите объем тела, ограниченного поверхностями x^{2}+y^{2}=a^{2}, x^{2}+z^{2}=b^{2} и b=2, a=3 (рассмотреть ту из областей внутри конуса, для которой х \ge 0 ).

    Относительная погрешность расчета не должна превышать 5%.

Контрольные вопросы

  1. Какой закон распределения случайной величины используется при вычислении кратных интегралов методом Монте-Карло?
  2. Чему равен интегрант трехкратного интеграла при вычислении объема тел методом Монте-Карло?
  3. От чего зависит точность метода Монте-Карло?
  4. Как определяется оценка "исправленного" среднеквадратического отклонения?
Дальше >>
< Лекция 3 || Лекция 4: 12 || Лекция 5 >

Вопросы и ответы

вопросов: 2
Мария Ястребинская
Мария Ястребинская

Добрый день. Я приступила сегодня к самостоятельному изучению курса "Моделирование систем". Хочу понять - необходимо ли отсылать мои решения практических заданий на сайт, (и если да - то где найти волшебную кнопку "Загрузить...") или практические задания остаются полностью на моей совести? (никто не проверяет, и отчётности по ним я предоставлять не обязана?)

P.S.: тьютора я не брала

ответить
алена зянтерекова
алена зянтерекова
ответить