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

Министерство образования и науки Российской Федерации
Ивановский государственный университет
Математический факультет
Кафедра вычислительной и прикладной математики
Бакалаврская работа по основной образовательной программе бакалавриата направления
«МАТЕМАТИКА.КОМПЬЮТЕРНЫЕ НАУКИ»
на тему
«Численное интегрирование с использованием метода Монте-Карло»
Выполнила:
Терзи Марина
Петровна студентка
4 курса дневного отделения
Научный руководитель:
Щаницина Светлана
Валерьевна
Иваново 2010 г.

Оглавление
Введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4 1 Теоретические основы метода Монте-Карло
6 1.1
Сведения из теории вероятностей и математической статистики
6 1.2
Описание метода Монте-Карло . . . . . . . . . . . . . . . . .
11 1.2.1
История метода Монте-Карло . . . . . . . . . . . . . .
11 1.2.2
Сущность метода Монте-Карло . . . . . . . . . . . . .
13 1.2.3
Погрешность метода Монте-Карло . . . . . . . . . . .
15 1.3
Генераторы случайных чисел . . . . . . . . . . . . . . . . . .
17 1.4
Интегрирование с помощью метода Монте-Карло
22 1.4.1
Способ усреднения подынтегральной функции . . . .
22 1.4.2
Способ существенной выборки, использующей вспомо- гательную плотность распределения . . . . . . . . . .
26 1.4.3
Способ, основанный на истолковании интеграла как площади . . . . . . . . . . . . . . . . . . . . . . . . . .
28 1.4.4
Способ выделения главной части . . . . . . . . . . . .
30 1.4.5
Вычисление кратных интегралов . . . . . . . . . . . .
32 1.5
Погрешность интегрирования методом Монте-Карло . . . . .
36 1.5.1
Формула оценки погрешности интегрирования мето- дом Монте-Карло . . . . . . . . . . . . . . . . . . . . .
36 1.5.2
Некоторые способы уменьшения дисперсии . . . . . .
39 2 Численный эксперимент
40 2.1
Алгоритм вычисления интегралов методом Монте-Карло . .
40 2.2
Программа для вычисления интегралов методом Монте-Карло 42 2.3
Пример 2.1. Интеграл от функции, не имеющей первообраз- ной в классе элементарных функций . . . . . . . . . . . . . .
44 2.4
Пример 2.2. Двойной интеграл . . . . . . . . . . . . . . . . .
46 2.5
Пример 2.3. Тройной интеграл. Объем восьмой части шара
47 2

2.6
Пример 2.4. Тройной интеграл. Объем тела, ограниченного поверхностями . . . . . . . . . . . . . . . . . . . . . . . . . .
48 2.7
Пример 2.5. Шестикратные интегралы в задаче о взаимном притяжение двух материальных тел . . . . . . . . . . . . . .
51
Заключение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
54
Литература . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
Приложение 1. Таблица равномерно распределенных слу- чайных цифр . . . . . . . . . . . . . . . . . . . . . . . . . .
56
Приложение 2. Таблица значений функции Φ(t) . . . . . .
57
Приложение 3. Код основной программы
58
Приложение 4. Код модифицированной программы
61 3

Введение
Методами Монте-Карло называют численные методы решения матема- тических задач при помощи моделирования случайных величин. Однако,
решать методами Монте-Карло можно любые математические задачи, а не только задачи вероятностного происхождения, связанные со случайными величинами.
Важнейшим приемом построения методов Монте-Карло явлеяется све- дение задачи к расчету математических ожиданий. Так как математиче- ские ожидания чаще всего представляют собой обычные интегралы, то центральное положение в теории метода Монте-Карло занимают методы вычисления интегралов.
Преимущества недетерминированных методов особенно ярко проявля- ются при решении задач большой размерности, когда применение тради- ционных детерминированных методов затруднено или совсем невозможно.
Границы между простым и сложным, возможным и невозможным суще- ствуют всегда, но с развитием вычислительной техники сдвигаются вдаль.
До появления электронных вычислительных машин (ЭВМ) методы Монте-
Карло не могли стать универсальными численными методами, ибо модели- рование случайных величин вручную — весьма трудоемкий процесс. Раз- витию методов Монте-Карло способствовало бурное развитие ЭВМ. Ал- горитмы Монте-Карло сравнительно легко программируются и позволя- ют производить расчеты во многих задачах, недоступных для классиче- ских численных методов. Так как совершенствование ЭВМ продолжается,
есть все основания ожидать дальнейшего развития методов Монте-Карло и дальнейшего расширения области их применения.
Цель работы – изучение особенностей применения метода Монте-Карло в решении задач численного интегрирования.
Для достижения цели были поставлены следующие задачи:
1) изучение теоретических основ метода Монте-Карло (идея метода,
алгоритмы расчета, способы оценки погрешности);
2) создание электронного информационного ресурса по теме Численное интегрирование с использованием метода Монте-Карло;
3) создание программы для приближенного вычисления интегралов с использованием метода Монте-Карло;
4

4) подбор примеров для проведения численного эксперимента.
При изучении теоретических основ метода Монте-Карло были проана- лизированы одиннадцать источников.
Основные сведения из теории вероятностей и математической статисти- ки, необходимые для понимания метода Монте-Карло, содержатся в [3], [7],
[9].
Для моделирования случайных процессов, связанных с применением ме- тода Монте-Карло, необходимы случайные числа; информацию о генера- торах случайных чисел можно найти в [4], [5], [6], [7].
Использование метода Монте-Карло для вычисления одномерных инте- гралов подробно обсуждается в [8].
Наиболее типичный из возможных способов вычисления многократных интегралов методом Монте-Карло описан в [6].
Некоторые способы уменьшения дисперсии, а, следовательно, ошибки интегрирования достаточно развернуто изложены в [4], [6].
Численные методы и программирование в системе MATLAB подроно описаны в [10].
Для проведения численных экспериментов методом Монте-Карло ис- пользуются некоторые примеры из [1], [2], [11].
В используемой литературе не приводится оценка погрешности вычис- ления кратных интегралов, в связи с чем в процессе выполнения работы возникла необходимость вывода формулы для оценки этой погрешности.
Изученный материал составляет основу теоретической главы электрон- ного информационного ресурса.
Результаты применения изложенной теории на практике получены с по- мощью разработанной программы и обсуждаются в главе 2.
Нумерация формул в работе организована следующим образом. Номер состоит из двух чисел, разделенных точкой, первое из которых соответству- ет номеру главы, а второе – номеру формулы внутри главы. Обозначения в тексте вводятся по мере их возникновения.
5

Глава 1
Теоретические основы метода
Монте-Карло
1.1
Сведения из теории вероятностей и математической стати- стики
Вероятностным пространством называется тройка (Ω, S, P (S)), где
1) Ω – некоторое заданное множество, называемое пространством элемен- тарных исходов; его элементы называются также точками;
2) S – непустое множество подмножеств Ω, являющееся σ-алгеброй, т.е.
замкнутым относительно операций суммирования (объединения) счетного числа своих элементов, счетного пересечения и дополнения множеством;
предполагается, что Ω (и, следовательно, пустое множество) является эле- ментом S;
3) P (S) – вероятностная мера, т.е. неотрицательная функция множеств из S такая, что P (Ω) = 1 и P

∞
P
i=1
A
i

=
∞
P
i=1
P (A
i
), где A
i
– последователь- ность попарно непересекающихся множеств из S.
Множества A из S называют событиями, а величины P (A) – их веро- ятностями.
Вероятность удовлетворяет следующим условиям:
1) P (A) ≥ 0;
2) P (Ω) = 1;
3) A, B ∈ S, A ∩ B = ∅ ⇒ P (A ∩ B) = P (A) + P (B);
4) A
1
⊇ A
2
⊇ . . . , ∩
∞
n=1
A
n
= ∅ ⇒ lim n→∞
P (A
n
) = 0.
Функция X(ω), ω ∈ Ω называют S-измеримой, если множество {ω :
X(ω) < x} ∈ S для любого вещественного числа x.
6

Случайной величиной, определенной на (Ω, S, P (S)), называется произ- вольная S-измеримая функция X(ω), принимающая конечные значения.
Для каждой случайной величины X(ω) однозначно определяется функ- ция распределения
F
X
(x) = P {ω : X(ω) < x} (F
X
(x) = P (X < x)),
(1.1)
представляющая собой монотонно неубывающую, неотрицательную, непре- рывную слева функцию такую, что lim x→−∞
F
X
(x) = 0, lim x→∞
F
X
(x) = 1.
Если случайная величина может принимать конечное или счетное число значений, то она называется дискретной.
Случайная величина называется абсолютно непрерывной, если ее функ- ция распределения представляется в виде
F (x) =
x
Z
−∞
f (t)dt,
при этом f (x) = F
0
(x) называется плотностью вероятности.
Плотность должна удовлетворять свойствам:
1) f (x) > 0;
2)
∞
R
−∞
f (x)dx = 1.
Математическим ожиданием дискретной случайной величины X на- зывается число
M [X] =
n
X
i=1
x i
p i
,
(1.2)
где x i
– возможные значения случайной величины X, а p i
– соответствую- щие им вероятности.
Вероятности p i
(i = 1, ..., n) должны удовлетворять условиям:
1) ∀i p i
> 0;
2)
n
P
i=1
p i
= 1.
7

Чтобы выяснить смысл величины M [X], запишем ее в следующем виде:
M [X] =
n
P
i=1
x i
p i
n
P
i=1
p i
Отсюда видно, что M [X] – это среднее значение величины X, причем более вероятные значения x i
входят в сумму с большими весами.
Математическим ожиданием непрерывной случайной величины назы- вается число
M [X] =
∞
Z
−∞
xf (x)dx.
(1.3)
Укажем формулу для математического ожидания случайной функции.
Пусть по-прежнему случайная величина X имеет плотность вероятностей f (x). Выберем произвольную непрерывную функцию g(x) и рассмотрим случайную величину Y = g(X), которую называют случайной функцией.
Тогда
M [g(X)] =
∞
Z
−∞
g(x)f (x)dx.
(1.4)
Отметим основные свойства математического ожидания:
если c – какая-нибудь не случайная величина, X и Y — две любые случайные величины, то
1) M [X + c] = M [X] + c;
2) M [c · X] = c · M [X];
для независимых случайных величин X и Y справедливо:
3) M [X · Y ] = M [X] · M [Y ].
Дисперсия D[X] – это математическое ожидание квадрата отклонения случайной величины X от ее среднего значения M [X]:
D[X] = M [(X − M [X])
2
].
(1.5)
Очевидно, D[X] > 0.
В случае, если случайная величина дискретная, то
D[X] =
n
X
k=1
(x k
− M [X])
2
P (X = x k
).
(1.6)
8

Если случайная величина непрерывна, то дисперсия
D[X] =
∞
Z
−∞
(x − M [X])
2
f (x)dx.
(1.7)
Формулу (1.5) для дисперсии можно преобразовать так:
D[X] = M [X
2
−2·X ·M [X]+(M [X])
2
] = M [X
2
]−2·M [X]·M [X]+(M [X])
2
,
откуда
D[X] = M [X
2
] − (M [X])
2
(1.8)
Вычислять дисперсию по этой формуле обычно проще.
Отметим основные свойства дисперсии:
если c — какая-нибудь не случайная величина, то
1) D[X + c] = D[X];
2) D[c · X] = c
2
M [X];
для независимых случайных величин X и Y справедливо:
3) D[X + Y ] = D[X] + D[Y ].
Случайная величина R, определенная в интервале (a, b) и имеющая плотность f (x) =
1
b−a
, называется равномерно распределенной в (a, b).
Известно, что для этой случайной величины математическое ожидание
M [R] =
b+a
2
, а дисперсия D[R] =
(b−a)
2 12
Нормальной (гауссовской) случайной величиной называется случайная величина ζ, определенная на всей оси (−∞, ∞) и имеющая плотность f (x) =
1
√
2πσ
e
−
(x−m)2 2σ2
, где m и σ > 0 - числовые параметры, причем математическое ожидание M [ζ] = m, дисперсия D[ζ] = σ
2
Правило трех сигм. Каковы бы ни были m и σ в f (x)
P {m − 3σ < X < m + 3σ} = 0.997,
или m+3σ
Z
m−3σ
f (x)dx = 0.997.
Вероятность 0,997 настолько близка к 1, что иногда последнюю формулу интерпретируют так: при одном испытании практически невозможно получить значение X, отличающееся от M [X] больше чем на 3σ.
9

Центральная предельная теорема. Если случайные величины
X
1
, ..., X
n
– независимые и одинаково распределенные, математическое ожидание M [X
i
] = m и дисперсия D[X
i
] = σ
2
, то
∀t
1
< t
2
lim n→∞
P







t
1
<
n
P
i=1
(X
i
− m)
√
nσ
< t
2







=
1 2π
t
2
Z
t
1
e
−
z2 2
dz
(1.9)
Выберем t
2
= −t
1
= t :
lim
N →∞
P
(
1
N
N
X
i=1
(X
i
− m)
< t r
D[X]
N
)
= 2Φ(t),
(1.10)
где
Φ(t) =
1 2π
t
Z
0
e
−
z2 2
dz.
(1.11)
Таким образом, при достаточно больших значениях N
P
(
¯
X − m
< t r
D[X]
N
)
≈ 2Φ(t).
(1.12)
Эта формула содержит целое семейство оценок, зависящих от параметра t. Если задать любой коэффициент доверия β, то можно найти (по табли- це значений функции Лапласа из приложения) корень t = t
β
уравнения
2Φ(t) = β. Тогда вероятность неравенства
¯
X − m
< t
β
r
D[X]
N
(1.13)
приблизительно равна β.
Чаще других используют коэффициенты доверия β = 0.95, которому отвечает x
β
= 1.96; или β = 0.997, которому отвечает x
β
= 3(«правило трех сигм»).
10

1.2
Описание метода Монте-Карло
1.2.1
История метода Монте-Карло
Метод Монте-Карло можно определить как метод моделирования слу- чайных величин с целью вычисления характеристик их распределений.
Определение метода можно сделать более общим, говоря не о случайных величинах, а о случайных функциях.
Идея моделирования случайных явлений известна давно и, по мнению некоторых авторов, восходит ко временам Древнего Вавилона и Ветхого
Завета. Что же касается использования такого рода явлений для целей приближенных вычислений, то первой работой в этой области принято считать работу Холла 1873 г. о вычислении числа π с помощью случай- ных бросаний иглы на разграфленную параллельными линиями бумагу.
Существуют также ряд более поздних работ, в которых до появления элек- тронных вычислительных машин (ЭВМ) использовались по существу идеи метода Монте-Карло. Однако, до появления ЭВМ этот метод не мог найти сколько-нибудь широкого применения, ибо моделировать случайные вели- чины вручную – очень трудоемкий процесс. Поэтому, идеи эти не получили заметного развития вплоть до 1944 г., когда в связи с работами по созданию атомной бомбы Дж. фон Нейман предложил широко использовать аппа- рат теории вероятностей для решения прикладных задач с помощью ЭВМ.
Первая работа, где этот вопрос систематически излагался, принадлежит
Метрополису и Уламу. В Советском Союзе первые статьи о методе Монте-
Карло были опубликованы в 1955—1956 годах (Чавчанидзе, Шрейдер).
В настоящее время имеется более двух тысяч работ, где исследуются теоретические основы метода или рассматриваются его применения к кон- кретным задачам. Метод используется для решения прикладных задач из большого числа областей науки и техники.
Первоначально метод Монте-Карло использовался главным образом для решения задач нейтронной физики, где традиционные численные методы оказались малопригодными. Далее его влияние распространилось на ши- рокий круг задач статистической физики, очень разных по своему содер- жанию. К разделам науки, где все в большей мере используется метод
11

Монте-Карло, следует отнести задачи теории массового обслуживания, за- дачи теории игр и математической экономики, задачи теории передачи со- общений при наличии помех и ряд других.
Метод Монте-Карло оказал и продолжает оказывать существенное влия- ние на развитие методов вычислительной математики (например, развитие методов численного интегрирования) и при решении многих задач успеш- но сочетается с другими вычислительными методами и дополняет их. Его применение оправдано в первую очередь в тех задачах, которые допускают теоретико-вероятностное описание. Это объясняется как естественностью получения ответа с некоторой заданной вероятностью в задачах с вероят- ностным содержанием, так и существенным упрощением процедуры реше- ния.
Наиболее сложными этапами решения той или иной задачи на ЭВМ в настоящее время следует считать математическое описание исследуемого явления, необходимые упрощения задачи, выбор подходящего численного метода, исследование его погрешности и запись алгоритма. В тех случа- ях, когда имеется теоретико-вероятностное описание задачи, использование метода Монте-Карло может существенно упростить упомянутые промежу- точные этапы. Во многих случаях полезно и для задач детерминированных строить вероятностную модель (рандомизовать исходную задачу) с тем,
чтобы далее использовать метод Монте-Карло.
Само название «Монте-Карло» происходит от города Монте-Карло в княжестве Монако, знаменитого своим игорным домом. Дело в том, что одним из простейших механических приборов для получения случайных ве- личин является рулетка. Термин «метод Монте-Карло» равнозначен терми- ну «метод статистических испытаний», также принятому в отечественной литературе. В зарубежной литературе обычно говорят о методах Монте-
Карло, имея в виду то обстоятельство, что подлежащие вычислению вели- чины могут быть оценены, исходя из различных вероятностных моделей
(например, случайные величины с различными законами распределения могут иметь одинаковое среднее значение). Часто, когда речь идет об изу- чении некоторых реальных явлений, то моделирование связанных с ними случайных величин процессов называют имитацией (simulation).