Вычисление кратных интегралов методом Монте – Карло

Используя встроенный в систему программирования генератор случайных чисел, вычислить двойной интеграл вида

Вид функции f(x,y) и пределы интегрирования выбрать в соответствии с Вашим вариантом (номер варианта равен порядковому номеру в списке группы, после 12-го номера «номер варианта = номер в списке группы – 12»).

Сравнить результат вычисления с точным значением данного интеграла.

Варианты:

f(x,y)	x1	x2	y1	y2	Вариант
ex+y			-2	-1
xЧsin(y)				p
ex-y			-3	-2
cos(x)Чsin(y)		p/2	-p
cos(x)Чcos(y)	-p			p
sin(x)Чsin(y)	-p		p	2p
sin(x)Чcos(y)	-p/2	p/2		p
x2y	-1
xЧy2			-1
sin(x) Чy		p
cos(x) Чy2		p/2	-2
xЧy3

Отчет по работе

а) Постановка задачи.

б) Расчетные соотношения.

в) Алгоритм вычислений (в виде блок-схемы).

г) Результаты расчета (приближенное и точное значения интеграла).

Примечание. Поскольку результат вычисления имеет случайный характер, следует вычислить интеграл N раз и провести статистическую обработку результатов расчета (представить результат в виде I ± DI).

Тема 5. Генераторы случайных чисел. Различные алгоритмы генерации случайных чисел.

Тема 6. Простейшие системы массового обслуживания и их параметры.

Тема 7. Модель «очередь».

Тема 8. Применение имитационного моделирования для решения задач массового обслуживания.

Тема 9. Методы прогнозирования: методы экстраполяции, методы экспертных оценок, методы логического моделирования, эвристические методы. Временные ряды.

Тема 10. Точность прогнозов. Точность прогнозов. Проблемы определения точности методов прогнозирования. Анализ точности. Сравнительные характеристики.

Самостоятельная работа

Программой предусмотрены дополнительные часы для самостоятельной работы студентов по каждой из тем, изучаемых на теоретических и практических занятиях.

№ работы	Название работы	Расчетная трудоемкость, часов
	Моделирование движения тел в среде с учетом трения
	Моделирование движения небесных тел и заряженных частиц
	Моделирование колебательных процессов
	Моделирование физических процессов в приближении сплошной среды
	Моделирование динамики численности популяций
	Моделирование очередей в системах массового обслуживания
	Моделирование различных случайных процессов
Итого:

Возможные темы:

1. Моделирование физических процессов (броуновское движение и др.)

2. Моделирование биологических процессов (распространение эпидемий)

3. Моделирование биологических процессов (конкурирующие популяции)

4. Решение геометрических задач с использованием ПК

5. Использование математических пакетов в моделировании

6. Имитационное моделирование на основе программного продукта MathCad

7. Использования пакета Maple в моделировании

8. Нахождение площадей произвольных фигур методом Монте-Карло

9. Использование метода Монте-Карло для решения интегралов

10. Генераторы случайных чисел

11. Задача о взаимозачётах (два способа решения)

12. Задача о нефтепроводе максимальной пропускной способности

13. Методы прогнозирования и их применение

14. Прогнозирующее моделирование в экономике

15. Составление и решение двойственных задач

ПРИМЕРЫ РЕШЕНИЯ ЗАДАЧ

Задача 1

Сняты следующие значения входа X и выхода Y у модели типа черный ящик

X
Y

Предложите описание и найдите коэффициенты модели. Оцените погрешность.

Решение задачи 1

Нарисуйте в системе координат экспериментальные точки. Так как график функции практически является линией, то логично для начала принять гипотезу о линейной функции

Y =A₁ X + A₀

Теперь следует найти неизвестные коэффициенты A₁ , A₀, причем так, чтобы график функции подходил к экспериментальным точка как можно ближе. То есть, чтобы сумма F расстояний E от каждой точки до линии была как можно меньше. Стоит задача минимизации суммарной ошибки F(A₁, A₀). Чтобы найти минимум функции F(A₁ , A₀) по переменным A₁ , A₀ надо найти производные MF/MA₁ и MF/MA₀ и приравнять их нулю, разрешить систему двух уравнений с двумя неизвестными A₁ , A₀.

E_i = Y_i - A₀ - A₁ X_i, i = 1,n , где n - число снятых экспериментально точек.

Получим систему из двух уравнений

Надо найти коэффициенты A₀ и A₁, для этого решаем систему методом Крамера, построив предварительно определитель следующего вида:

Подставляя числа X и Y, число экспериментов n=6, имеем: A₁ = 0,71 , A₀ = 3,32.

Чтобы определить, принимается гипотеза или нет, нужно рассчитать ошибку между теоретической и экспериментальной зависимостями.

То есть F=31.84, F = 2,3 .

Строим параллельно графику функции Y = 0.71* X + 3.32две линии на расстоянии от него +2,3 и –2,3. В коридор [-2.3, +2.3] попадет 5 экспериментальных точек из шести, то есть 5/6*100%=83%, что удовлетворяет требованию нормального распределения ошибки (67%). Поэтому гипотеза о том, что данную функцию можно считать линейной выполняется.

Задача 2

Дана система из нескольких взаимодействующих тел (см. рис.). Рассматривается процесс теплопроводности. Предложите математическую модель системы. Определите, является ли система открытой? Рассчитайте движения системы методом Эйлера (5-8 шагов).

Решение задачи 2

Исследуется система двух материальных тел А и В с различными теплофизическими свойствами. Система контактирует с опорой Тн и помещена во внешнюю среду с температурой Тс. Интересует протекание процесса изменения температур тел.
Как видно, в процессе жизни в системе изменяются (могут измениться) четыре показателя - температуры тел A, B, Тс, Тн. Значит, мы имеем дела с четырьмя переменными от времени (поскольку переменные изменяют свои значения со временем). Введем эти переменные:X1(t), X2(t), X3(t), X4(t).

Для построения математической модели данной системы процесс теплопередачи изображен в виде графа зависимостей. Стрелка от А к В обозначает изменение температуры X2(t) объекта В под влиянием объекта А. Ряд стрелок отсутствует, то есть отсутствует влияние одних параметров на другие. Например, от B к Тс, так как тело В не в состоянии сколько-нибудь существенно нагреть открытую атмосферу. Строго говоря, такое влияние есть, но оно настолько ничтожно, что разумно им пренебречь.

Массивную опору тоже нагреть не удастся. Поскольку переменных четыре, то нам необходимо четыре закона, описывающих их изменение. В общем виде, учитывая, от каких переменных зависит каждый показатель, имеем:

- для тела A X1(t) зависит от температуры атмосферы Тс и температуры тела; В;
- для тела В X2(t) зависит от температуры атмосферы Тс, температуры тела A и температуры опоры Тн.

Стрелки, входящие в соответствующий элемент графа, указывают на количество влияющих параметров, а то, откуда они исходят, определяет конкретные названия переменных.

Для среды закон имеет вид X3(t)= const - температура атмосферы Тс не зависит от остальных составляющих данной системы и, соответственно, не изменяется. Для нагревателя закон имеет вид X4(t)= const - температура опоры Тн не зависит от остальных составляющих данной системы и, соответственно, не изменяется.

Основной динамический закон для описания изменения переменной имеет вид:
dX/dt = w( x(t), y(t), z(t), ...)

Рассмотрим первое уравнение dX1(t)/dt = f1( X1(t), X2(t), X3(t) ). Какие пары переменных взаимодействуют? Стрелки соединяют X1(t) с X2(t), X1(t) с X3(t). То есть имеет место два процесса, влияющих на темп. Мы рассматриваем процессы теплообмена тел. Известно, что тепло, переданное от одного тела, складывается с теплом, переданным от другого. Таким образом имеем:

dX1(t)/dt= g1(X1(t), X2(t)) + g2(X3(t), X1(t)).

Раскроем структуру оставшихся выражений g1 и g2. Очень удобно, что g1 ни как не зависит от g2 и может рассматриваться отдельно. Забудем на некоторое время о g2.

Какой знак нужно поставить между X1(t) и X2(t) в выражении g1?
Возможные варианты: X1(t) + X2(t) , X2(t) - X1(t) , X1(t) - X2(t), X1(t) * X2(t), X1(t) / X2(t) и другие.

Следует начать с наиболее простых - природа построена просто. И только, если простейшие не удовлетворяют, переходят к более сложным вариантам описания.

Попробуем вариант: dX1(t)/dt= X2(t)-X1(t)

Какие есть качественные варианты у этой физической системы ?

а).Х1>X2. Тело А теплее тела В. Теплопоток при контакте двух тел направлен от А к В. Тело А отдает тепло телу В. То есть в процессе контакта значение Х1 падает - уменьшается. Так ли это в уравнении? Посмотрим.

Х1>X2, значит Х2-X1<0, значит dX1(t)/dt<0, значит Х1 падает. Вывод не противоречит физической картине. Значит, пока данный вариант приемлем, и надо проверить его на остальных качественных ситуациях.

б).X1<X2. Тело А холоднее тела В. Теплопоток при контакте двух тел направлен от В к А. принимает тепло То есть в процессе контакта значение Х2 растет - увеличивается. Так ли это в уравнении? Посмотрим.

Х1<X2, значит Х2-X1>0, значит dX1(t)/dt>0, значит Х1 растет. Вывод не противоречит физической картине. Значит, пока данный вариант приемлем, и надо проверять его далее.

в). X1=X2. Температура тела А равна температуре тела В. Теплопоток при контакте двух тел равен нулю. То есть значение Х1 не изменяется - тело А не отдает и не принимает тепло. Так ли это в уравнении? Посмотрим.

Х1=X2, значит Х2-X12=0, значит dX1(t)/dt=0, значит Х1 не изменяется. Вывод не противоречит физической картине. Значит, данный вариант принимается, так как он правильно (пока только качественно!) отражает физическую картину во всех случаях.

Других вариантов существования системы нет, рассмотрение оканчивается.

Забыв на некоторое время о g1, также можно рассмотреть и g2. В итоге получаем:

dX1(t)/dt= (X2(t)-X1(t)) + (X3(t)-X1(t)).

Так как, во-первых, у разных материалов разность температур влияет на темп изменения температуры тела различным способом, во-вторых, скорости двух процессов (двух разных пар металлов) могут быть разными, то скорректируем модель коэффициентом теплопроводности, который играет роль усилителя (ослабителя) процессов. Это коэффициент влияния связи на объект. При К=0 влияние отсутствует, связь отключается. При К=0.0001 влияние слабое. При К=1000 влияние связи огромно. Понятно, что коэффициент стоит при выражении процесса К?(X2(t)-X1(t)), где знак ? означает знак некоторой операции. Какой? Умножение. Эта операция дает ЗАВИСИМОСТЬ одного члена от другого ( в нашем случае К от (X2(t)-X1(t))).

В итоге модель имеет вид: dX1(t)/dt= К21*(X2(t)-X1(t)) + К31*(X3(t)-X1(t)).

Теперь аналогично второе уравнение:

dX2(t)/dt= К12*(X1(t)-X2(t)) + К32*(X3(t)-X2(t)) + К42*(X4(t)-X2(t)).

Уравнение изменения температуры опоры: dX4(t)/dt=0.
Уравнение изменения температуры атмосферы: dX3(t)/dt=0.
Вся система уравнений в сборе имеет вид:
dX1(t)/dt=К21*(X2(t)-X1(t))+К31*(X3(t)-X1(t))
dX2(t)/dt=К12*(X1(t)-X2(t))+К32*(X3(t)-X2(t))+К42*(X4(t)-X2(t))
dX3(t)/dt=0
dX4(t)/dt=0

Ясно, что по физическим соображениям - сколько тепла вытекает из А в В, столько же и поступает в В из А, то есть К21=К12.

Далее заметим, что мы получили открытую систему, то есть такую, чье суммарное тепло не постоянно, а может изменяться. Это видно из асимметрии стрелок на графе. Проверим этот факт. Для этого сложим левые части всех уравнений и, отдельно, правые части.

вместе: dXсистемы/dt=К31*(X3(t)-X1(t))+К32*(X3(t)-X2(t))+К42*(X4(t)-X2(t)),
то есть dXсистемы/dt не равно нулю и есть утечка или приток тепла в систему извне.

Имеем слева: dX1(t)/dt+dX2(t)/dt+dX3(t)/dt+dX4(t)/dt или d(X1+X2+X3+X4)/dt или dXсистемы/dt

справа: К21*(X2(t)-X1(t)) + К31*(X3(t)-X1(t)) + + К12*(X1(t)-X2(t)) + К32*(X3(t)-X2(t)) + К42*(X4(t)-X2(t)) или К31*(X3(t)-X1(t)) + К32*(X3(t)-X2(t)) + К42*(X4(t)-X2(t))

Рассмотрим применение метода Эйлера для расчета траектории движения системы. Зададим значения коэффициентов модели:

Х3=22 град, Х4=15 град, k1=0.1 1/с, k2=0.1 1/с, k3=0.1 1/с, k4=0.05 1/с.

Подставим значения коэффициентов:
X1(t + ) = X1(t)+[0.1*(22-X1(t))+0.2*(X2(t)-X1(t))] *
X2(t + )=X2(t)+[0.1*(15-X2(t))+0.05*(22-X2(t))+0.2*(X1(t)-X2(t))]*
Зададим начальные условия системы
Х1(0)=30 град ;
Х2(0)=70 град.
Конечное значение времени моделирования - tk = 4 сек..
Выбираем шаг моделирования, например, =0,2с и приступаем к моделированию.
Процесс моделирования и численные значения отражены в таблице.

Задача 3

Сгенерируйте случайное трехразрядное число, распределенное по равномерному закону в интервале от 0 до 1, с помощью монеты.

Решение задачи 3

Подбросьте монету 10 раз, если монета упала решкой, то запишите 0, если орлом, то 1. Итак, допустим, что в результате эксперимента получили случайную последовательность 100110100.

Начертите интервал от 0 до 1. Считывая числа в последовательности слева направо, разбивайте интервал пополам и выбирайте каждый раз одну из частей (если 0, то левую, если 1, то правую). Таким образом, можно добраться до любой точки интервала, сколь угодно точно.

Итак, 1, интервал делится пополам [0 – 0,5] и [0,5 - 1]. Выбираем правую половину. Интервал сужается: [0.5 - 1].

Следующее число 0. Интервал [0.5 - 1] делится пополам: [0,5 – 0,75] и [0,75 - 1]. Выбираем левую половину [0,5 – 0,75]. Интервал сужается: [0,5 – 0,75].

Следующее число 0. Интервал [0.5 – 0,75] делится пополам: [0,5 – 0,625] и [0,625 - 0,75]. Выбираем левую половину [0,5 – 0,625]. Интервал сужается: [0,5 – 0,625].

Следующее число 1. Интервал [0.5 – 0,625] делится пополам: [0,5 – 0,5625] и [0,5625 – 0,625]. Выбираем правую половину [0,5625 – 0,6250]. Интервал сужается: [0,5625 – 0,6250].

По условию точности задачи решение найдено, им являются любое число из интервала [0,5625 – 0,6250], например, 0,625.

Способ 2.

Разобьем двоичное число на триады: 100, 110, 100. После перевода двоичных чисел в десятичные получаем 4, 6, 4. Проведем операцию конкатенации, подставив спереди «0,»: 0,464. Так как в принципе таким методом могут получаться только числа от 0,000 до 0,777, то следует получаемые числа масштабировать на интервал от 0 до 1. Коэффициент масштаба равен 1/0.777.

Итак, искомое число равно 0,464/0,777 = 0,597.

Задача 4

Предложите пример с возможным исходом четырех несовместных случайных событий. Задайте их вероятности. Промоделируйте выпадение последовательности событий.

Решение задачи 4

Выбрать из колоды карт наугад карту. Определить масть карты.

В колоде 36 карт четырех мастей по 9 карт каждой масти. Интервал от 0 до 1 разделим на равные четыре части: [0,0-0.25], [0,25-0,50], [0,5-0,75], [0,75-1,0]. Первая часть будет соответствовать картам масти Ч, вторая – картам масти П, третья – картам масти В, четвертая – Б. Взять случайное равномерно распределенное число в интервале от 0 до 1 (смотри задачу 3). Например, 0,597. Данное число попадает в третий интервал, соответствующий масти В.

Поскольку теперь в колоде 9 карт масти Ч, 9 карт масти П, 8 карт масти В, 9 карт масти Б, то интервал от 0 до 1 будет разбит на отрезки длиной: 9/35, 9/35, 8/35, 9/35.

То есть [0.00-0.257], [0,257-0,474], [0,474-0,743], [0,743-1,00]. Разыграем случайное равномерно распределенное число в интервале от 0 до 1. Например, 0,321. Данное число попадает во второй интервал, соответствующий масти П.

Продолжая процесс можно получить (в зависимости от конкретных случайных чисел), например, такую последовательность: В-П-В-Ч-Б-П-Ч-…

Задача 5

Сгенерируйте поток из 10 случайных событий с интенсивностью появления событий 5 шт/час.

Решение задачи 5

Возьмем случайные числа, равномерно распределенные в интервале от 0 до 1 (смотри задачу 3 или приложение) и вычислим их логарифмы:

0.0333 -3.4022

0.3557 -1.0337

0.2172 -1.5269

0.5370 -0.6218

Формула пуассоновского потока определяет расстояние между двумя случайными событиями как: t = -Ln(rnd)/8. Тогда, учитывая 8=5, имеем расстояния между двумя случайными соседними событиями: 0.68, 0.21, 0.31, 0.12.

То есть события наступают: первое в момент t=0, второе в t=0.68, третье в t=0.89, четвертое в t=1.20, пятое в t=1.32 и так далее.

Задача 6

Проимитируйте работу робота на примере вероятностной дискретной марковской цепи из 4 звеньев - ориентирование манипулятора, захват детали, потеря детали, укладка детали в гнездо. Определите среднее количество тактов, требуемое для проведения всей комплексной операции.

Решение задачи 6

Обозначим состояния марковской цепи: А-ориентирование манипулятора, В-захват детали, С-потеря детали, Д-укладка детали в гнездо. Введем вероятности переходов из i-го состояния в j-тое состояние таким образом, чтобы сумма вероятностей переходов из i-го состояния в остальные смежные с ним была равна 1 (переходы – полная группа несовместных события).

Таблица вероятностей переходов	в А	в В	в С	в Д	Сумма вероятностей переходов
из А	0,7	0,2		0,1	0,7+0,2+0,1=1
из В	0,7		0,3		0,7+0,3=1
из С
из Д

(Потренируйтесь рисовать на основании таблицы граф переходов).

Проимитируем, используя таблицу случайных чисел, процесс работы робота. Пусть начальное состояние будет А. Последовательность из таблицы случайных чисел:

0,51 (робот остается в состоянии А, так как 0.51<0,7),

0,89 (робот переходит из А в состояние В, так как 0,7<0.89<0.7+0,2),

0,63 (робот переходит из В в состояние А, так как 0.63<0,7),

0,31 (робот остается в состоянии А, так как 0.31<0,7),

0,92 (робот переходит из А в состояние Д, так как 0.7<0,7+0,2<0.92<0.7+0.2+0.1),

(робот переходит из Д в состояние А с вероятностью 1).

Полный цикл совершен за 6 тактов: А-А-В-А-А-Д-А.

Повторяя данную имитацию, можно получить, например, такие реализации (это зависит от того, какие конкретно случайные числа выпадут):

4 (АВАДА), 9 (АВСАВСААДА), 5 (АВСАДА), 6 (ААВСАДА), 12 (АВСААВСААВАДА), 6 (АВСААДА), 5 (АВААДА), 10 (АВСАВСАВАДА), 8 (ААВАВСАДА). Всего уложено 10 деталей. Среднее число циклов: (6+4+9+5+6+12+6+5+10+8)/10=7.1 или, округляя, 7.

Задача 7

Проимитируйте временной диаграммой работу системы массового обслуживания с интенсивностью прихода заявок 5 шт/час, двумя каналами обслуживания с интенсивностью обслуживания по 2 и 1 шт/час. и ограниченной двумя местами очередью к каждому из них. Определите пропускную способность системы массового обслуживания по временной диаграмме и вероятность обслуживания заявки.

Решение задачи 7

Аналогично тому, как это делалось в задаче 5, сгенерируем время наступления событий – прихода заявок в систему.

То есть события наступают: первое в момент t=0, второе в t=0.68, третье в t=0.89, четвертое в t=1.20, пятое в t=1.32 и так далее.

Первая заявка поступает сразу в канал обслуживания номер 1 и обслуживается там в течение времени t = -Ln(rnd)/3=-Ln(0.1958)/2=1.6307/2=0,815 час. В момент 0,815 заявка поступает на выход системы.

Вторая заявка поступает в систему в t=0.68 и, так как первый канал к этому моменту занят, поступает в канал обслуживания номер 2 и обслуживается там в течение времени t = -Ln(rnd)/1=-Ln(0.7003)/1=0.3562/1=0,3562 час. В момент 0,68+0,3562=1,0362 (час) заявка поступает на выход системы.

Третья заявка поступает в систему в t=0.89 и, так как первый канал к этому моменту освободился, то заявка поступает в канал обслуживания номер 1 и обслуживается там в течение времени t = -Ln(rnd)/2=-Ln(0.9499)/2=0.0514/2=0,0257 час. В момент 0,89+0,0257=0,9157 (час) заявка поступает на выход системы.

Наиболее удобно все построения и вычисления производить на временной диаграмме. Пример такой диаграммы с дальнейшими построениями приведен на рисунке.

Рис.1. Временная диаграмма моделирования системы массового обслуживания

Время моделирования – 3 часа. По диаграмме видно, что за это время в систему поступило 14 заявок, отказано – 1 заявке, обслужено – 10 заявок. Пропускная способность системы равна количеству обслуженных заявок за время моделирования, отнесенное к этому времени, то есть 10/3=3,3 [шт/час]. Вероятность обслуживания равна отношению количества обслуженных заявок к числу желавших обслужиться за тот же период, то есть 10/14=0,7.

Задача 8

Определите методом Монте-Карло площадь пятиугольника с координатами углов - (0,0), (0,10), (5,20), (10,10), (7,0).

Решение задачи 8

Нарисуйте в двухмерных координатах заданный пятиугольник. Наименьшее значение X=0, наибольшее значение X=10, наименьшее значение Y=0, наибольшее значение Y=20. Нарисуйте прямоугольник с координатами углов - (0,0), (0,20), (10,20), (10,0). Данный прямоугольник содержит пятиугольник. Площадь прямоугольника легко найти, и она равна: (10-0)*(20-0) = 200.

Используем таблицу случайных чисел для генерации пар чисел R,G, равномерно распределенных в интервале от 0 до 1. Поскольку одно число R будет имитировать координату X (0<=X<=10), то X=10*R. Поскольку второе число G будет имитировать координату Y (0<=Y<=20), то Y=20*G.

Сгенерируем 10 пар чисел R,G и отобразим 10 точек X,Y на графике.

Таблица 1