Просмотр кода страницы Участник:Grachovamc/Вычисление кратного интеграла методом Монте - Карло.

== Свойства и структура алгоритма ==

=== Общее описание алгоритма ===
Метод Монте-Карло — общее название группы численных методов, основанных на получении большого числа реализаций стохастического процесса, который формируется таким образом, чтобы его вероятностные характеристики совпадали с аналогичными величинами решаемой задачи. Метод Монте-Карло имеет множество различных приложений. Он применяется для решения задач в различных областях математики, физики, химии, биологии, экономики, оптимизации, теории управления и др.

Годом рождения метода Монте-Карло считается 1949 год, когда в свет выходит статья Метрополиса и Улама «Метод Монте-Карло»<ref>''N. Metropolis, S. Ulam The Monte Carlo Method, — J. Amer. statistical assoc. 1949 44 № 247 335—341.''</ref>. Название метода происходит от названия города в княжестве Монако, широко известного своими многочисленными казино, поскольку именно рулетка является одним из самых широко известных генераторов случайных чисел.

Особенностью данного метода является простая структура вычислительного алгоритма. Как правило, составляется программа для осуществления одного случайного испытания, затем это испытание повторяется N раз, причем каждый опыт не зависит от остальных. Для получения итогового значения усредняются результаты всех опытов. Поэтому метод Монте - Карло также иногда называют методом статистических испытаний<ref>''Соболь И. М. Метод Монте-Карло. — М.: Наука, 1968. — 64 с. — (Популярные лекции по математике). — 79 000 экз.''</ref>.

Большим преимуществом метода Монте-Карло является то, что он позволяет учесть в модели элемент случайности и сложность реального мира. Кроме того, метод является робастным по отношению к изменению различных параметров, таких как распределение случайной величины.

Рассмотрим использование метода Монте — Карло для задачи вычисления кратных интегралов.

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

Для определения этой площади можно разбить отрезок, по которому производится интегрирование, на подотрезки, подсчитать площадь под графиком функции на каждом из них и сложить. И если предположить, что в случае функции одной переменной достаточно разбиения на 50 подотрезков и, соответственно, 50  вычислений значений функции, то в случае функции <math>n</math> переменных количество подотрезков и вычислений значений функции возрастает до <math>50^n</math>. При размерности функции больше 10 задача становится очень громоздкой. Поскольку пространства большой размерности встречаются во многих задачах, необходимо иметь метод решения, вычислительная сложность которого бы не столь сильно зависела от размерности. Именно таким свойством обладает метод Монте-Карло.
=== Математическое описание алгоритма ===
==== Интегрирование функции одной переменной ====
Рассмотрим два наиболее распространенных варианта использования метода Монте - Карло для вычисления интегралов.
===== Метод усреднения подынтегральной функции =====
Рассмотрим интеграл вида

<math>I = \int_a^b \phi(x) dx</math> 

Пусть <math>X</math>  - случайная величина, равномерно распределенная на отрезке интегрирования <math>[a, b]</math> с плотностью <math>f(x) = \frac{1}{b - a}</math>. Тогда <math>\phi(X)</math> также является случайной величиной, и ее математическое ожидание выражается как

<math>\mathbb{E}\phi(X)=\int\limits_a^b \phi(x)f(x) dx = \frac{1}{b - a}\int\limits_a^b \phi(x) dx</math>

Отсюда

<math>\int_a^b \phi(x) dx = (b - a)\mathbb{E}\phi(X)</math>

Заменив математическое ожидание <math>\mathbb{E}\phi(X)</math> его оценкой – выборочным средним, получим оценку интеграла

<math>\int\limits_{a}^{b}\phi(x)\,dx\approx\frac{b-a}{N}\sum^{N}_{i=1}\phi(x_i)</math>,

где <math>x_i</math> – возможные значения случайной величины <math>X</math>, <math>N</math> – число испытаний.

===== Геометрический метод =====
Рассмотрим интеграл вида

<math>I = \int_a^b \phi(x) dx</math>, где подынтегральная функция неотрицательна и ограничена <math>0 \le \phi(x) \le c</math> для любого <math>x \in [a, b]</math>.

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

* функция ограничивается прямоугольником площадью <math>S</math>;
* в данном прямоугольнике случайным образом выбираются координаты <math>N</math> точек;
* определяется число <math>n</math> точек, координаты которых попали под график функции; тогда выполняется соотношение <math>\frac{\int\limits_a^b \phi(x) dx}{S} \approx \frac{n}{N}</math>
* искомая оценка интеграла равна <math>S\frac{n}{N}</math>.

==== Вычисление кратного интеграла ====
Случай вычисления <math>n</math> - мерного интеграла аналогичен одномерному случаю с учетом некоторых корректировок.
===== Метод усреднения подынтегральной функции =====

Пусть функция <math>y = f(x_1, x_2,...,x_n)</math> непрерывна в ограниченной замкнутой области D и требуется вычислить n-кратный интеграл I по области D:

<math>I = {\int\limits_{D}{\ldots \int{f\left( {{x}_{1}},\ldots ,{{x}_{n}} \right)d{{x}_{1}}\ldots d{{x}_{n}}}}}</math>

Геометрически число <math>I</math> представляет собой <math>(m + 1)</math> - мерный объем вертикального цилиндрического тела в пространстве <math>Ox_1x_2...x_ny</math>, построенного на основании D и ограниченного сверху соответствующим куском поверхности <math>y = f(\vec x)</math>,где <math>\vec x = (x_1, ..., x_n)</math>.

Пусть область интегрирования D расположена в <math>n</math>-мерном параллелепипеде <math>a_k \le x_k \le b_k</math>, <math>k = 1,...,n</math>.
Преобразуем интеграл таким образом, чтобы новая область интегрирования целиком находилась внутри единичного <math>n</math> - мерного куба. Для этого сделаем замену переменных <math>x_k = a_k + (b_k - a_k) y_k</math>, <math>k = 1,...,n</math>. Тогда <math>0 \le y_k \le 1</math>, <math>k = 1,...,n</math>. Следовательно новая область интегрирования G целиком расположена в единичном <math>n</math>-мерном кубе.

Таким образом, получаем выражение <math>I = \frac{D(x_1,...x_n)}{D(y_1,...y_n)}I'</math>, где якобиан преобразования имеет вид 

<math>\frac{D(x_1,...x_n)}{D(y_1,...y_n)} = (b_1 - a_1)...(b_n - a_n)</math>,

<math>I' = {\int\limits_{G}{\ldots \int{g\left( {{y}_{1}},\ldots ,{{y}_{n}} \right)d{{y}_{1}}\ldots d{{y}_{n}}}}} = {\int\limits_{G}{\ldots \int{f\left( {a_1 + (b_1 - a_1)y_1},\ldots ,{a_n + (b_n - a_n)y_n} \right)d{{y}_{1}}\ldots d{{y}_{n}}}}}</math>,

а область интегрирования G расположена внутри <math>n</math> - мерного единичного куба <math>0 \le y_k \le 1</math>, <math>k = 1,...,n</math>.

Выберем N равномерно распределенных на отрезке <math>[0, 1]</math> последовательностей случайных величин <math>\vec y^{(i)} = ({y^{(i)}}_{1}, ..., {y^{(i)}}_{n})</math>, <math>i = 1,...,N</math>. Эти последовательности являются координатами случайных точек <math>K_i({y^{(i)}}_{1}, ..., {y^{(i)}}_{n})</math>. Далее определим, какие из этих точек принадлежат области интегрирования G. Пусть таких точек всего m. Тогда при достаточно большом m можно считать, что <math>\mathbb{E} g(K) = \frac{1}{m}\sum^{m}_{i=1}g(K_i) </math>. Следовательно оценка интеграла <math>I'</math> имеет вид 
<math>I' = \mathbb{E} g(K) * V_G</math>, где <math>V_G</math> - <math>n</math> - мерный объем области интегрирования G.
Таким образом, оценка искомого интеграла имеет вид <math>I = V_{par} * \mathbb{E} g(K) * V_G = \frac{V_{par}}{N}\sum^{m}_{i=1}g(K_i)</math>, где <math>V_{par} = (b_1 - a_1)...(b_n - a_n)</math> - объем параллепипеда, ограничивающего исходную область интегрирования D.

===== Геометрический метод =====
В случае использования геометрического метода осуществляется та же последовательность действий, что и при использовании данного метода в случае функции одной переменной с учетом того, что функция ограничивается <math>n</math> - мерным параллелепипедом, и, соответственно, в результирующей формуле вместо площади прямоугольника берется объем ограничивающего функцию параллелепипеда.
=== Вычислительное ядро алгоритма ===
Таким образом, вычислительным ядром последовательной реализации любого из представленных вариантов метода Монте - Карло интегрирования является многократное вычисление значения функции в случайной точке.
=== Макроструктура алгоритма ===
Согласно описанию вычислительного ядра алгоритма в качестве макровершины можно взять вычисление подынтегральной функции в случайной точке. Данное действие выполняется в алгоритме N раз.
=== Схема реализации последовательного алгоритма ===
==== Интегрирование функции одной переменной ====
===== Метод усреднения подынтегральной функции =====
Последовательность исполнения метода следующая:

1. Генерация последовательности <math>{{\{x_i\}}_{i = 1}}^{N}</math> значений случайной величины, равномерно распределенной на отрезке интегрирования <math>[a, b]</math>. Для экономии памяти лучше генерировать не всю случайную последовательность целиком, а по одному значению в цикле. 

2. Вычисление значений функции <math>f(x_i)</math> в каждой из точек последовательности <math>{{\{x_i\}}_{i = 1}}^{N}</math>.

3. Вычисление выборочного среднего <math>\frac{1}{N}\sum^{N}_{i=1}f(x_i)</math>.

4. Получение искомого значения интеграла согласно формуле <math>\int\limits_{a}^{b}\phi(x)\,dx\approx\frac{b-a}{N}\sum^{N}_{i=1}\phi(x_i)</math>.

===== Геометрический метод =====
Пусть интегрируемая функция имеет следующий вид

[[Файл:0_step.png]]

Последовательность исполнения метода следующая. В цикле осуществляется выполнение последовательности действий:

1. Функция ограничивается прямоугольником площадью <math>S</math>

[[Файл:1_step.png]]

2. В прямоугольнике случайным образом выбираются координаты <math>N</math> точек

[[Файл:3_step.png]]

3. Определяется число <math>n</math> точек, координаты которых попали под график функции

[[Файл:4_step.png]]

4. Искомая площадь вычисляется согласно формуле <math>S\frac{n}{N}</math>

==== Вычисление кратного интеграла ====
===== Метод усреднения подынтегральной функции =====
Последовательность выполнения метода следующая:

1. Преобразование интеграла таким образом, чтобы новая область интегрирования G целиком находилась внутри единичного <math>n</math> - мерного куба с использованием замены переменных <math>x_k = a_k + (b_k - a_k) y_k</math>, <math>k = 1,...,n</math>.

2. Вычисление якобиана преобразования и получение следующего выражения получаем выражение <math>I = \frac{D(x_1,...x_n)}{D(y_1,...y_n)}I'</math>, где якобиан преобразования имеет вид 

<math>\frac{D(x_1,...x_n)}{D(y_1,...y_n)} = (b_1 - a_1)...(b_n - a_n)</math>, <math>I' = {\int\limits_{G}{\ldots \int{g\left( {{y}_{1}},\ldots ,{{y}_{n}} \right)d{{y}_{1}}\ldots d{{y}_{n}}}}} = {\int\limits_{G}{\ldots \int{f\left( {a_1 + (b_1 - a_1)y_1},\ldots ,{a_n + (b_n - a_n)y_n} \right)d{{y}_{1}}\ldots d{{y}_{n}}}}}</math>.

3.Выбор N равномерно распределенных на отрезке <math>[0, 1]</math> последовательностей случайных величин <math>\vec y^{(i)} = ({y^{(i)}}_{1}, ..., {y^{(i)}}_{n})</math>, <math>i = 1,...,N</math>, которые являются координатами случайных точек <math>K_i({y^{(i)}}_{1}, ..., {y^{(i)}}_{n})</math>. 

4. Определение того, какие из этих точек принадлежат области интегрирования G. Пусть таких точек всего m.

5. Вычисление <math>\mathbb{E} g(K) = \frac{1}{m}\sum^{m}_{i=1}g(K_i) </math>.

6. Получение оценки интеграла <math>I'</math>: <math>I' = \mathbb{E} g(K) * V_G</math>, где <math>V_G</math> - <math>n</math> - мерный объем области интегрирования G.

7. Получение оценки искомого интеграла <math>I = V_{par} * \mathbb{E} g(K) * V_G = \frac{V_{par}}{N}\sum^{m}_{i=1}g(K_i)</math>, где <math>V_{par} = (b_1 - a_1)...(b_n - a_n)</math> - объем параллепипеда, ограничивающего исходную область интегрирования D.
==== Геометрический метод ====
Схема реализации последовательного алгоритма аналогична одномерному случаю.

=== Последовательная сложность алгоритма ===
Для вычисления определенного интеграла методом Монте - Карло требуется осуществить:

1) при использовании метода усреднения подынтегральной функции

* N генераций случайных величин;
* N вычислений значений функции;
* N - 1 сложение.

2) при использовании геометрического метода

* N генераций случайных величин;
* N вычислений значений функции;
* N сравнений.

Таким образом, алгоритм является линейным относительно количества случайных точек.
=== Информационный граф ===
[[Файл:My perfect graph.png]]
=== Ресурс параллелизма алгоритма ===
=== Входные и выходные данные алгоритма ===
=== Свойства алгоритма ===
== Программная реализация алгоритма ==
=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
=== Масштабируемость алгоритма и его реализации ===
=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===
== Литература ==