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

Министерство образования и науки Украины

Сумской Государственный Университет

Кафедра информатики

Курсовая работа

по численным методам

на тему

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

Сумы 2007

СОДЕРЖАНИЕ

1 История рождения метода Монте-Карло

1.1 Алгоритм Буффона для определения числа Пи

1.2 Связь случайных процессов и дифференциальных уравнений

1.3 Рождение метода Монте-Карло в Лос-Аламосе

1.4 Дальнейшее развитие и современность

2 Использование метода Монте-Карло в численном интегрировании

2.1 Численное интегрирование

2.2 Одномерный случай

2.3 Многомерный случай

3 Применение метода Монте-Карло для вычисления кратных интегралов (на примере двукратных интегралов)

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

3.2 Программная реализация метода

3.2.1 Математическое описание метода

3.2.2 Алгоритм метода в программе

3.2.3 Описание основных значений в программе

3.2.4 Инструкция по работе с программой

3.2.5 Результат программы

3.2.6 Реализация метода в пакете Mathcad

Вывод

Литература

Приложение А (код программы)

1 ИСТОРИЯ РОЖДЕНИЯ МЕТОДА МОНТЕ-КАРЛО

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

1.1 Алгоритм Буффона для определения числа Пи

Случайные величины использовались для решения различных прикладных задач достаточно давно. Примером может служить способ определения числа Пи, который был предложен Буффоном еще в 1777 году. Суть метода была в бросании иглы длиной N на плоскость, расчерченную параллельными прямыми, расположенными на расстоянии d друг от друга (см. Рис. 1).

Рисунок 1 - Метод Буффона

Вероятность того, что отрезок пересечет прямую связана с числом Пи:

,

где A - расстояние от начала иглы до ближайшей к ней прямой;

и - угол иглы относительно прямых.

Этот интеграл просто взять:

(при условии, что d > L), поэтому подсчитав долю отрезков, пересекающих прямые, можно приближенно определить это число. При увеличении количества попыток точность получаемого результата будет увеличиваться.

В 1864 году капитан Фокс, выздоравливая после ранения, чтобы как-то занять себя, реализовал эксперимент по бросанию иглы. Результаты представлены в следующей таблице:

	Число бросаний	Число пересечений	Длина иглы	Расстояние между прямыми	Вращение плоскости	Значение Пи
Первая попытка	500	236	3	4	отсутствует	3.1780
Вторая попытка	530	253	3	4	присутствует	3.1423
Третья попытка	590	939	5	2	присутствует	3.1416

Комментарии:

· Вращение плоскости применялось (и как показывают результаты - успешно) для того, чтобы уменьшить систематическую ошибку.

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

1.2 Связь случайных процессов и дифференциальных уравнений

Создание математического аппарата случайных методов началось в конце 19го века. В 1899 году лорд Релей показал, что одномерное случайное блуждание на бесконечной решетке может давать приближенное решение параболического дифференциального уравнения. Колмогоров в 1931 году дал большой толчок развитию случайных подходов к решению различных математических задач, поскольку он сумел доказать, что цепи Маркова связаны с некоторыми интегро-дифференциальными уравнениями. В 1933 году Петровский показал, что случайное блуждание, образующее Марковскую цепь асимптотически связано с решением эллиптического дифференциального уравнения в частных производных. После этих открытий стало понятно, что, случайные процессы можно описывать дифференциальными уравнениями и, соответственно, исследовать при помощи хорошо на тот момент разработанных математических методов решения этих уравнений.

1.3 Рождение метода Монте-Карло в Лос-Аламосе

Однако сначала Энрико Ферми в 1930х годах, а затем Джон фон Нейман и Станислав Улам в Лос-Аламосе предположили, что можно использовать связь между случайными процессами и дифференциальными уравнениями «в обратную сторону». Они предложили использовать случаный подход для аппроксимации многомерных интегралов в уравнениях переноса, возникших в связи с задачей о движении нейтрона в изотропной среде.

Идея была развита Станиславом Уламом, который, по иронии судьбы, также как и Фокс боролся с вынужденным бездельем во время выздоровления после болезни, и, раскладывая пасьянсы, задался вопросом, какова вероятность того, что пасьянс «сложится». Ему в голову пришла идея, что вместо того, чтобы использовать обычные для подобных задач соображения комбинаторики, можно просто поставить «эксперимент» большое число раз и, таким образом, подсчитав число удачных исходов, оценить их вероятность. Он так же в дальнейшем предложил использовать компьютеры для расчетов методом Монте-Карло.

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

Годом рождения метода Монте-Карло считается 1949 год, когда в свет выходит статья Метрополиса и Улама «Метод Монте-Карло». Название метода происходит от названия города в княжестве Монако, широко известного своими многочисленными казино, поскольку именно рулетка является одним из самых широко известных генераторов случайных чисел. Станислав Улам пишет в своей автобиографии «Приключения математика», что название было предложено Метрополисом в честь его дяди, который был азартным игроком.

1.4 Дальнейшее развитие и современность

В 1950х годах метод использовался для расчетов при разработке водородной бомбы. Основные заслуги в развитии метода в это время принадлежат сотрудникам лабораторий ВВС США и корпорации RAND.

В 1970х годах в новой области математики - теории вычислительной сложности было показано, что существует класс задач, сложность (количество вычислений, необходимых для получения точного ответа) которых растет с размерностью задачи экспоненциально. Иногда можно, пожертвовав точностью, найти алгоритм, сложность которого растет медленнее, но есть большое количество задач, для которого этого нельзя сделать (например, задача определения объема выпуклого тела в n-мерном евклидовом пространстве) и метод Монте-Карло является единственной возможностью для получения достаточно точного ответа за приемлемое время. В настоящее время основные усилия исследователей направлены на создание эффективных Монте-Карло алгоритмов различных физических (компьютерное моделирование играет в современной физике важную роль и метод Монте-Карло является одним из самых распространенных во многих областях от квантовой физики до физики твердого тела, физики плазмы и астрофизики.), химических и социальных процессов для параллельных вычислительных систем.

2 ИСПОЛЬЗОВАНИЕ МЕТОДА МОНТЕ-КАРЛО В ЧИСЛЕННОМ ИНТЕГРИРОВАНИИ

2.1 Численное интегрирование

Численное интегрирование - вычисление значения определённого интеграла (как правило, приближенное), основанное на том, что величина интеграла численно равна площади криволинейной трапеции, ограниченной осью абсцисс, графиком интегрируемой функции и отрезками прямых

и ,

где и - пределы интегрирования (см. рис.2).

Рисунок 2 - Определённый интеграл как площадь фигуры

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

2.2 Одномерный случай

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

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

Точки называются узлами метода, числа - весами узлов. При замене подынтегральной функции на полином нулевой, первой и второй степени получаются соответственно методы прямоугольников, трапеций и парабол (Симпсона). Часто формулы для оценки значения интеграла называют квадратурными формулами.

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

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

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

Рисунок 3 - Численное интегрирование функции детерминистическим методом

Предположим, что для функции, представленной на рис. 3, достаточно разбиения на 25 отрезков и, следовательно, вычисления 25 значений функции. Представим теперь, мы имеем дело с n-мерной функцией. Тогда нам необходимо 25ⁿ отрезков и столько же вычислений значения функции. При размерности функции больше 10 задача становится огромной. Поскольку пространства большой размерности встречаются, в частности, в задачах теории струн, а также многих других физических задачах, где имеются системы со многими степенями свободы, необходимо иметь метод решения, вычислительная сложность которого бы не столь сильно зависла от размерности. Именно таким свойством обладает метод Монте-Карло.

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

· ограничим функцию прямоугольником (n-мерным параллелепипедом в случае многих измерений), площадь которого S_par можно легко вычислить;

· «набросаем» в этот прямоугольник (параллелепипед) некоторое количество точек (N штук), координаты которых будем выбирать случайным образом;

· определим число точек (K штук), которые попадут под график функции;

· площадь функции S дается следующим выражением:

Для малого числа измерений интегрируемой функции производительность Монте-Карло интегрирования гораздо ниже, чем производительность детерминистических методов. Тем не менее, в некоторых случаях, когда функция задана не явно, а необходимо определить область, заданную в виде сложных неравенств, случайный метод более предпочтительный.

2.3 Многомерный случай

В небольших размерностях можно так же применять квадратурные формулы, основанные на многочленах Лагранжа. Однако в больших размерностях эти методы становятся неприемлемыми из-за быстрого возрастания числа точек сетки и/или сложной границы области. В этом случае применяется метод Монте-Карло. Генерируются случайные точки в нашей области и усредняются значения функции в них. Так же можно использовать смешанный подход - разбить область на несколько частей, в каждой из которых (или только в тех, где интеграл посчитать не удаётся из-за сложной границы) применить метод Монте-Карло.

3 ПРИМЕНЕНИЕ МЕТОДА МОНТЕ-КАРЛО ДЛЯ ВЫЧИСЛЕНИЯ КРАТНЫХ ИНТЕГРАЛОВ (НА ПРИМЕРЕ ДВУКРАТНЫХ ИНТЕГРАЛОВ)

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

Вычислить двукратный интеграл методом Монте-Карло

3.2 Программная реализация метода

3.2.1 Математическое описание метода

Метод Монте-Карло состоит в том, что рассматривается некоторая случайная величина о, математическое ожидание которой равно искомой величине z, т.е.

Осуществляется серия n независимых испытаний, в результате которых получается (генерируется) последовательность n случайных величин : , , …, и по совокупности этих значений приближенно определяется искомая величина, т.е.

.

Интеграл же может быть вычислен как математическое ожидание некоторой случайной величины о, которая определяется независимыми реализациями з_i случайной величины з с равномерным законом распределения. Двукратные интегралы вычисляются следующим образом:

, (1)

где , - независимые реализации равномерно распределенных на величин.

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

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

Поскольку эти числа генерируются по наперед заданному алгоритму, то они не совсем случайны (псевдослучайны), хотя и обладают свойственными случайным числам статистическими характеристиками. То есть случайные величины, равномерно распределенные на отрезке , в ЭВМ задаются с помощью специальных программ - генераторов псевдослучайных чисел.

3.2.2 Алгоритм метода в программе

1. Ввод количества испытаний n.

2. Запустить генератор случайных чисел.

3. Определение общих пределов интегрирования.

4. Случайным образом (с учетом общих пределов интегрирования) генерировать значения x и y.

5. Если значения x и y лежат в текущих пределах интегрирования, то к значению интеграла прибавить значение функции при этих значениях.

6. Пункты 4 и 5 повторить n раз.

7. Окончательно полученную сумму (из формулы 1) перемножить на оставшиеся части формулы (1) для получения значения интеграла.

3.2.3 Описание основных значений в программе

n	Количество испытаний
DInt	Основная подпрограмма
F	Функция для вычисления значений точек подынтегральной функции
x, y	Аргументы функции
S	Текущее значение суммы (интеграла)
Ax, Ay, Bx, By	Пределы интегрирования
i	Счетчик

В программе MonteKarlo используются стандартные модули WinCrt, Strings.

В программе MonteKarlo описан метод Монте-Карло, который осуществляется с помощью функций: function F(x,y: real): real; - расчет значения функции и function DInt(A: real; n: longint): real; - вычисление двойного интеграла.

3.2.4 Инструкция по работе с программой

Для начала работы программы запустить файл MonteKarlo.pas. Перед вами появится окно:



Запустите программу меню RunRun или Ctrl+F9, перед вами появится окно:

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

3.2.5 Результат программы

3.2.6 Реализация метода в пакете Mathcad

Mathcad позволяет вычислять кратные интегралы непосредственно, однако в большинстве случаев при кратности интегралов 3 и более применение самого метода Монте-Карло предпочтительнее. Дело в том, при одинаковой точности метод Монте-Карло дает существенный выигрыш во времени (в десятки и сотни раз), особенно при большой кратности интегралов. Идея метода состоит в том, что интеграл заменяется величиной Fср.·V, где V - объем области интегрирования, Fср. - среднее значение подынтегральной функции, вычисленное по нескольким случайно выбранным точкам.

Определим подынтегральную функцию.

И вычислим интеграл



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

,

где N - общее число точек, n - число точек, лежащих внутри области интегрирования, V - объем области, включающей область интегрирования.

ВЫВОД

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

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

Я реализовала метод Монте-Карло в приложении BPW, а также в пакете Mathcad. Программная реализация метода показывает, что при увеличении количества генерирования случайных чисел возрастает время подсчета интеграла. Mathcad позволяет вычислять кратные интегралы непосредственно, однако в большинстве случаев при кратности интегралов 3 и более применение самого метода Монте-Карло предпочтительнее. Дело в том, при одинаковой точности метод Монте-Карло дает существенный выигрыш во времени (в десятки и сотни раз), особенно при большой кратности интегралов.

В настоящее время основные усилия исследователей направлены на создание эффективных Монте-Карло алгоритмов различных физических химических и социальных процессов для параллельных вычислительных систем.

Метод Монте-Карло используется очень часто, порой некритично и неэффективным образом. Он имеет некоторые очевидные преимущества:

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

б) Он приводит к выполнимой процедуре даже в многомерном случае, когда численное интегрирование неприменимо, например, при числе измерений, больше 10.

в) Его легко применять при малых ограничениях или без предварительного анализа задачи.

Он обладает, однако, некоторыми недостатками, а именно:

а) Границы ошибки не определены точно, но включают некую случайность. Это, однако, более психологическая, чем реальная, трудность.

б) Статическая погрешность убывает медленно.

в) Необходимость иметь случайные числа.

ЛИТЕРАТУРА

1. Статья «Моделируя жизнь», автор Андрей Тепляков.

2. Книга «Fundamentals of the Monte Carlo method for neutral and charged particle transport», автор Alex F Bielajew (на английском).

3. Статья «Metopolis, Monte Carlo and the MANIAC».

4. Статья о Монте-Карло на www.riskglossary.com.

5. ISBN 5030033920 Д. Каханер, К. Моулер, С. Нэш. Численные методы и программное обеспечение (пер. с англ.). М.: Мир, 2001, 575 c.

6. Численные методы на Mathcad.