Вероятностные методы. Метод Монте-Карло

В отличие от традиционных численных методов реше­ния задач, заключающихся в разработке алгоритма пост­роения решения, для исследования некоторых классов задач оказывается более целесообразным моделирование их сущности с использованием законов больших чисел теории вероятностей. Оценки f1, f2, ... , fn искомой ве­личины f в этом случае строят на основании статистичес­кой обработки материала, полученного в результате многократных случайных испытаний. Основным требо­ванием при этом является сходимость по вероятности рассматриваемой случайной величины к искомому реше­нию задачи:

(9.1)

Здесь Р – вероятность, ε – сколь угодно малое поло­жительное число.

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

 

9.1. Вычисление определенного интеграла

Пусть задана непрерывная случайная величина ξ (кси), с плотностью вероятности р(х), значения которой распре­делены на интервале (а, b)

Плотность вероятности р(х) обладает следующими свойствами:

 

1) p(x) ³ 0;

2) . (9.2)

 

Тогда математическое ожидание случайной величины ξ равно интегралу

M(ξ) = . (χ – хи)

Для функции f(ξ), аргументом которой является слу­чайная величина ξ, т. е. для случайной функции, матема­тическое ожидание равно

M(f((ξ)) = . (9.3)

Отсюда следует, что любой интеграл вида

где функция р(х) обладает свойствами (9.2), можно при­нять за математическое ожидание некоторой случайной функции f(ξ).

Но математическое ожидание случайной величины f(ξ) можно приближенно определить с помощью статистичес­ких оценок, т. е. на основе выборки {ξ1, ξ2,..., ξN} объема N как среднее арифметическое значений {f(ξ1), f(ξ2),..., f(ξN)} поэтому интеграл (9.3) можно вычислить прибли­женно по формуле

» (9.4)

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

Среднее арифметическое N испытаний случайной ве­личины ξ

ξN =

 

является случайной величиной с тем же математическим ожиданием

M(ξN) = М(ξ)

и дисперсией, равной

D(ξN) =

и при N ® ¥ закон распределения слу­чайной величины ξN стремится к нормальному закону.

Очевидно, что чем больше N, тем меньше дисперсия D(ξN) = . Величину погрешности формулы (9.4)можно оценить по вероятности. Например, при боль­ших N можно утверждать, что с вероятностью 0,997 ошибка не превосходит величину 3σ = 3× (прави­ло «трех сигм» для нормально распределенной случайной величины).

Можно считать, что погрешность формулы (9.4) есть величина порядка O , но для повышения точности в данном случае нельзя применять правило Рунге-Ромберга.

Приведем другой способ статистической оценки для одномерного интеграла

.

Для этого вспомним его геометрический смысл. Предпо­ложим, что функция f(x) положительна на отрезке [a, b]. Тогда интеграл равен площади криволинейной трапеции, ограниченной графиком функции f(x), осью абсцисс и прямыми х = а и х = b (рис. 9.1).

Рассмотрим две случайные величины ξ – равномерно распределенную на отрезке [а, b] и η (эта) – равномерно рас­пределенную на отрезке [0, fmax] где fmax = .

Рис. 9.1

Вероятность попадания случайной точки (ξ, η) в криво­линейную трапецию равна отношению площади трапеции к площади прямоугольника

{(х, у), a £ х £ b, 0 £ у £ fmax }:

(9.5)

 

Проведем серию из N испытаний и получим N случай­ных точек (х, у), принадлежащих прямоугольной области {a £ х £ b, 0 £ у £ fmax }. Подсчитаем число Nf точек, удов­летворяющих условию у £ f(x). Тогда вероятность попада­ния случайной точки (ξ, η) в криволинейную трапецию приближенно равна относительной частоте попадания в криволинейную трапецию, т. е. р » и интеграл приближенно вычисляется по формуле

» (9.6)

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

Проведем серию из N испытаний случайной величины, равномерно распределенной на отрезке [а, b], и полу­чим N случайных чисел xi, принадлежащих отрезку [а, b]. Вычислим интеграл по формуле

» (9.7)

Отметим, что в (9.7) подынтегральная функция может принимать положительные и отрицательные значения, тогда как формула (9.6) применима только для неотрица­тельной функции f(x).

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

(9.8)

и применить формулу (9.7).

9.2. Вычисление кратных интегралов

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

При использовании квадратурных или кубатурных формул, число операций быстро возрастает с ростом раз­мерности интеграла. Например, если для вычисления од­номерного интеграла методом трапеций с заданной точностью необходимо вычислить сумму порядка N слагаемых, то для вычисления двойного интеграла тем же методом необходимо сложить порядка N2 слагаемых, а для тройно­го интеграла число слагаемых составляет порядка N3.

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

Применяется следующий критерий выбора между кубатурной форму­лой р-го порядка точности и методом Монте-Карло для вычисления с точностью ε кратного интеграла функции m переменных:

1) если число измерений m < 2р, лучше использовать кубатурные или квадратурные формулы;

2) если m > 2рме­тод Монте-Карло.

Например, если р = 1, тройные интегралы выгоднее вычислять методом Монте-Карло, а одномерные – квад­ратурными формулами.

Если р = 2, лучше вычислять методом Монте-Карло пя­тимерные интегралы, а одномерные, двойные и трой­ные – квадратурными или кубатурными формулами.

Рассмотрим конкретные формулы метода Монте-Кар­ло для вычисления кратных интегралов, получающиеся способом, который применялся для вывода формулы (9.7).

Пусть требуется вычислить двойной интеграл

. (9.9)

Проведем серию из N испытаний случайной точки (xi, yi), где xi равномерно распределены на отрезке [a, b], a yi равномерно распределены на отрезке [с, d]. Вычислим интеграл (9.9) по формуле

= (9.10)

 

Для тройного интеграла аналогично получим формулу

 

= (9.11)

 

где xi равномерно распределены на отрезке [a, b], yi – на отрезке [с, d], a zi – на отрезке [р, q]; N – число испы­таний.

Для m-кратного интеграла формула метода Монте-Карло имеет вид

 

. (9.12)

9.3. Решение систем линейных уравнений

Рассмотрим применение метода Монте-Карло для ре­шения системы линейных уравнений

 

, i = 1, ..., n. (9.13)

Пусть система (9.13) приведена к виду

, i = 1, ..., n. (9.14)

где норма матрицы α = || aij || удовлетворяет условию ||α|| < 1. Тогда система (9.14) имеет единственное решение.

Приведем без доказательства схематическое описание метода Монте-Карло для решения системы линейных уравнений (9.14).

Подберем такие множители vij, чтобы решения pij уравнений αij = pij×vij удовлетворяли условиям:

1) pij ³ 0, причем pij > 0, если αij ¹ 0;

2) i =1, ..., n.

Положим pi, n+1 = l – , i = 1, ..., n;

p n+1, j = 0, если j < n + l; pn+1, n+1 = l.

Будем трактовать число pij как вероятность перехода некоторой частицы из состояния Si в состояние Sj. Всего состояний п + 1:

S1, S2, ... , Sn+1,

причем граничное состояние Sn+1 = Г соответствует ос­тановке частицы, так как рп+1, j = 0, т. е. частица не мо­жет выйти из состояния Sn+1.

Матрица с элементами pij называется матрицей пере­хода состояний {Si}:

П . (9.15)

Назовем траекторией Ti совокупность состояний, че­рез которые проходит блуждающая частица, начиная с состояния Si и заканчивая граничным состоянием Г. Про­межуточные состояния частицы обозначим ;

Ti ={ = Г }. (9.16)

Поставим в соответствие траектории Ti случайное чис­ло Xi:

Xi = βi + × + × + ××× + × , (9.17)

где – свободные члены системы (9.14) с индексами, совпадающими с индексами состояний, обра­зующих траекторию (9.16).

Теорема 9.1. Математические ожидания случайных величин (9.17) равны корням системы (9.14): M(Xi) = xi, i = 1, ... , п.

 

Сформулируем алгоритм решения системы линейных уравнений (9.13) методом Монте-Карло.

1. Привести систему (9.13) к виду (9.14):

, i = 1, ..., n, ||α|| < 1. (9.18)

2. Подобрать такие числа vij, чтобы решения pij урав­нений αij = pij×vij удовлетворяли условиям:

1) pij ³ 0, причем pij > 0, если αij ¹ 0;

2) i =1, ..., n. (9.19)

Вычислить (п + 1)-й столбец и (п + 1)-ю строку матрицы перехода:

pi, n+1 = l – , i = 1, ..., n;

p n+1, j = 0, если j < n + l; pn+1, n+1 = l. (9.20)

 

3. Выполнить для каждого i = 1, ... , п:

3.0. присвоить значение xi = 0;

3.1. для каждого k = 1, ... , N (где N число испытаний случайной величины Xi) выполнить статистические ис­пытания по алгоритму:

3.1.0. присвоить m = 0; im = i; v = 1; xi = xi + βi;

3.1.1. присвоим переменной очередное значение: m = m + 1;

3.1.2. возьмем случайное число ξ из интервала (0; 1), и выберем значение im по следующему правилу:

im = , 1 £ j £ n

3.1.3. если j < n + 1, перейти к 3.1.1;

3.1.4. для каждого l = 0, 1, ... , m – 2 вычислить

v = v× , xi = xi + v× ;

3.2. вычислить xi = xi / N.

 

<== предыдущая лекция | следующая лекция ==>
Уравнения в частных производных | Признаки и переменные

Дата добавления: 2020-08-31; просмотров: 600;


Поиск по сайту:

Воспользовавшись поиском можно найти нужную информацию на сайте.

Поделитесь с друзьями:

Считаете данную информацию полезной, тогда расскажите друзьям в соц. сетях.
Poznayka.org - Познайка.Орг - 2016-2024 год. Материал предоставляется для ознакомительных и учебных целей.
Генерация страницы за: 0.023 сек.