Моделирование систем. Моделирование систем 2е издание, исправленное
 Скачать 5.72 Mb.
|
|
1. Что из себя представляет график функции плотности равномерного закона распределения случайных величин? 2. Как формируется поток Эрланга -го порядка? 3. Чему соответствует область определения функции распределения случайных величин экспоненциального закона? 4. Чему соответствует область определения функции распределения случайных величин равномерного закона? 5. Какая связь между стандартным отклонением и дисперсией случайных величин? 6. Что собой отображает гистограмма случайных величин? 7. Какая связь между функцией плотности и функцией распределения случайных непрерывных величин? 8. Что является параметром в экспоненциальном законе распределения случайных величин? Какой функциональный смысл имеет параметр экспоненциального распределения случайной величины? 9. Назовите допустимую область определения функции распределения потока Эрланга -го порядка. Выборочный метод Монте-Карло Цель работы: Изучить и практически освоить метод Монте-Карло на примерах расчета площадей плоских фигур, объемов пространственных тел, а также вычисления кратных интегралов. Среда программирования — MATLAB. Теоретическая часть Сущность метода Монте-Карло состоит в следующем: требуется найти значение некоторой изучаемой величины. Для этого выбирают такую случайную величину , математическое ожидание которой равно , т. е. [6] Практически же поступают следующим образом: производят испытаний, в результате которых получают возможных значений величины ; вычисляют их среднее арифметическое и принимают в качестве оценки (приближенного значения) искомого числа : Поскольку метод Монте-Карло требует произведения большого числа испытаний, его часто называют методом статистических испытаний. 1. Оценка погрешности метода Монте-Карло Как отмечалось, для получения оценки математического ожидания случайной величины необходимо произвести независимых испытаний и по ним найти выборочную среднюю, которая принимается в качестве искомой оценки. При каждой конечной серии испытаний будут получаться различные значения случайной величины и, следовательно, другая средняя, а значит, и другая оценка математического ожидания. Очевидно, что получить точную оценку математического ожидания невозможно. Поэтому возникает вопрос о допускаемой ошибке. Обычно ограничиваются отысканием лишь верхней границы допускаемой ошибки с заданной вероятностью , т. е. При этом возможны следующие случаи оценки числа испытаний: 1. Случайная величина распределена нормально и ее среднее квадратическое отклонение (стандартное отклонение) известно. В этом случае с заданной вероятностью верхняя граница ошибки определяется по формуле (4.1) где: — число испытаний (разыгранных значений случайной величины ); — значение аргумента функции Лапласа или интеграла вероятности, при котором она равна половине заданной вероятности; — известное среднее квадратическое отклонение. Из формулы (4.1) может быть найдено число испытаний. Один из вариантов интеграла вероятностей (функции Лапласа) имеет вид [6] (4.2) Значения табулированы и приведены в большинстве учебников по теории вероятностей и математической статистике. В зарубежной литературе большое распространение получила так называемая функция ошибок ( ) : (4.3) Связь между функцией ошибок (4.3) и интегралом вероятностей (4.2) выражается в виде (4.4) 2. Случайная величина распределена нормально, причем ее среднее квадратическое отклонение неизвестно. В этом случае с заданной вероятностью верхняя граница ошибки вычисляется по формуле (4.5) где: — число испытаний; — "исправленное" среднее квадратическое отклонение; находят по специальным таблицам, например, приведенной в [6] Из формулы (4.5) может быть найдено число испытаний для определения верхней границы ошибки. 3. Случайная величина распределена по закону, отличному от нормального. В этом случае при достаточно большом числе испытаний ( ), с вероятностью, приближенно равной (заданной вероятностью), верхняя граница ошибки может быть вычислена по формуле (4.1), если среднее квадратическое отклонение случайной величины известно; если же оно неизвестно, то можно подставить в формулу (4.1) его оценку — "исправленное" среднее квадратическое отклонение — либо воспользоваться формулой (4.5). При этом чем больше число испытаний, тем меньше различие между результатами, которые дают обе формулы. 2. Вычисление кратных интегралов методом Монте-Карло В случае когда, например, определенный интеграл не может быть вычислен в квадратурах, либо прибегают к численным методам интегрирования, либо расчет ведется с помощью метода Монте-Карло. Применение метода Монте-Карло становится оправданным при кратности интеграла больше трех. В данной лабораторной работе мы используем метод Монте-Карло для расчета интегралов с кратностью не более трех. Это позволит более ясно представить технику применения метода. Сначала рассмотрим вычисление простого определенного интеграла (4.6) где Введем под знак интеграла постоянный множитель (равный единице): (4.7) Вынесем из-под интеграла числитель дроби, получим (4.8) Как известно, если случайная величина распределена на заданном интервале (например, ) равномерно, то ее функция плотности обратно пропорциональна длине интервала, т. е. Кроме того, если известно распределение случайной величины , то функция от этой случайной величины будет иметь тот же самый закон распределения. В этом случае математическое ожидание непрерывной равномерно распределенной случайной величины рассчитывается по формуле Соответственно, математическое ожидание от функции случайной величины будет определяться следующим образом: (4.9) Сопоставляя (4.8) и (4.9), приходим к выводу, что определенный интеграл может быть рассчитан по формуле (4.10) Несмещенной оценкой математического ожидания случайной величины, как известно, является ее среднее арифметическое. Поэтому математическое ожидание можем приближенно найти по формуле (4.11) С учетом (4.11) получаем выражение для приближенного расчета определенного интеграла (4.12) Чем больше число испытаний , тем точнее будет расчет математического ожидания (4.11) и, следовательно, определенного интеграла, вычисляемого по формуле (4.12). Рассмотрим общий подход вычисления -кратного интеграла с помощью метода Монте-Карло [5] Пусть задан -кратный интеграл вида (4.13) где подынтегральная функция задана на замкнутой области Погрузим область интегрирования в -мерный промежуток (4.14) имеющий меру (4.15) Определим в промежутке (4.15) функцию (4.16) Тогда в соответствии с (4.13) и (4.16) получим (4.17) Введем в рассмотрение -мерную случайную величину , имеющую в замкнутой области равномерное распределение вероятностей с дифференциальной функцией плотности (4.18) Функция плотности равномерного распределения есть величина постоянная, поэтому введем ее под знак интеграла (4.17) следующим образом: (4.19) Вынесем числитель дроби за знак интеграла, т. е. (4.20) В (4.20) -кратный интеграл — это математическое ожидание от функции случайной величины в предположении, что случайная величина распределена равномерно с плотностью (4.18). Следовательно, можем записать (4.21) где В свою очередь математическое ожидание может быть оценено с помощью арифметического среднего. Тогда приближенное значение -кратного интеграла будет определяться приближенной формулой (4.22) где — значение случайной величины в -м испытании. Чтобы смоделировать выборку -мерной случайной величины , равномерно распределенной в -мерном промежутке , используются псевдослучайные числа. Для этого в каждом испытании с номером выбирают псевдослучайных чисел , и по ним определяют координаты случайной величины , псевдослучайной точки [5] Таким образом, техника применения метода Монте-Карло здесь будет заключаться в определении области , генерировании в ней псевдослучайных чисел, подсчета числа попаданий этих чисел в область и применении формулы (4.22). Расчет площадей и объемов можно рассматривать как частный случай вычисления кратных интегралов. Например, вычисление объема тел с помощью трехкратного интеграла сводится к взятию интеграла по области при подынтегральной функции, тождественно равной единице. Практическая часть Пример 1. Рассчитайте с помощью метода Монте-Карло площадь круга радиусом и с координатами центра . Число испытаний примите 500. Выполните графические построения, поясняющие расчет. Для определения площади плоской фигуры ее вписывают в соответствующую известную фигуру, площадь которой достаточно просто вычисляется. Вычисление площади круга может быть произведено через вычисление площади квадрата, в который вписывается круг. Выборочный метод Монте-Карло предполагает генерирование случайных (псевдослучайных) равномерно распределенных чисел. Этим числам сопоставляются координаты точек для рассматриваемой фигуры — квадрата, в которую вписывается круг. Если площадь квадрата подсчитана, то задается число испытаний , из которых исходов могут оказаться внутри круга или на его границе, т. е. на окружности. Тогда площадь круга будет определяться выражением где: — число сгенерированных случайных чисел, соответствующих количеству точек, находящихся в данном квадрате: — число случайных чисел (точек), которые попали в круг. Граница круга — это окружность, уравнение которой имеет вид где: , — координаты центра окружности; , — текущие значения переменных; — радиус окружности. Программный код решения примера: 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. Напишите программу расчета числа по методу Монте-Карло. Расчет числа произведите с точностью в четыре знака после десятичной точки. Пример 2. Рассчитайте с помощью метода Монте-Карло площадь фигуры, ограниченной эллипсом с большой осью и малой осью . Координаты пересечения осей эллипса , , оси коллинеарны декартовым осям. Теоретическая площадь эллипса вычисляется по формуле где: — большая ось эллипса; — малая ось эллипса. Каноническое уравнение эллипса (относительно центра координат): Программный код решения примера: 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. Координаты центра пересечения осей эллипса примите случайными при использовании функции 4. Оси эллипса примите случайными, равномерно распределенными из интервала , где — номер компьютера, за которым выполняется лабораторная работа ( .). 5. В программу включите расчет общего числа испытаний, при котором достигается заданная относительная погрешность в 2%. 6. В программу включите построение графика относительной погрешности от числа испытаний начиная с и более. Пример 3. Найдите оценку значения тройного интеграла методом Монте-Карло где — область, ограниченная круговым конусом и плоскостью , , Программный код решения примера: 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. Формирование случайных переменных , , вынесите за пределы цикла. 2. Объясните расстановку пределов интегрирования. 3. Вычислите по методу Монте-Карло объем тела для полученных пределов интегрирования. Сравните со значением, вычисленным аналитически. Задание 4 В нижеприводимых заданиях выполните графические построения заданных плоских фигур и напишите программы по расчету площадей этих фигур в системе MATLAB. 1. Рассчитайте с помощью метода Монте-Карло площадь фигуры, ограниченной гипоциклоидой с параметрами ; ; . Параметрическое уравнение гипоциклоиды: Подберите массив значений параметра . Относительная погрешность расчета не должна превышать 5%. 2. Рассчитайте с помощью метода Монте-Карло площадь фигуры, ограниченной кривой Штейнера с параметром Параметрическое уравнение кривой Штейнера: Подберите массив значений параметра . Относительная погрешность расчета не должна превышать 5%. 3. Рассчитайте с помощью метода Монте-Карло площадь фигуры, ограниченной астроидой с параметром . Параметрическое уравнение астроиды: Подберите массив значений параметра . Относительная погрешность расчета не должна превышать 5%. Задание 5 В нижеприводимых заданиях выполните графические построения заданных пространственных тел и напишите программы по вычислению объемов этих тел в системе MATLAB. 1. Вычислите с помощью метода Монте-Карло объем тела, ограниченного параболоидом и плоскостью Относительная погрешность расчета не должна превышать 5%. 2. Вычислите с помощью метода Монте-Карло объем тела, ограниченного поверхностями: Относительная погрешность расчета не должна превышать 5%. 3. С помощью метода Монте-Карло вычислите объем тела, ограниченного поверхностями и (рассмотреть ту из областей внутри конуса, для которой ). Относительная погрешность расчета не должна превышать 5%. Контрольные вопросы |