6Simulation systems Лекция 19 Монте-Карло (Материалы к лекциям), страница 5

Целесообразно обратить внимание на ряд возможных обобщений.

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

Поскольку прямоугольник АВED не является единичным, целесообразно изменить масштаб по осям координат.

X* = ; y* = Тогда квадрат АВED, будет единичным, а искомая площадь S=DB  AE  S^* , где S^* - площадь фигуры, выраженная в единицах измерения, соответствующих новым масштабам.

Величина S^* может быть определена методом статистических испытаний. Процедура решения задачи, в основном, совпадает с рассмотренной выше, за исключением того, что теперь для каждого x_j нужно вычислять два значения y₁(x_j) и y₂(x_j) ординат точек, лежащих на границе фигуры, и проверять справедливость неравенства: y₁(x_j)  y_j  y₂(x_j), причем y_j выбирается уже не из отрезка [0,1], а из преобразованного (сдвинутого) единичного отрезка.

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

S = f(x)dx

Поэтому процедура вычисления площади является процедурой вычисления интеграла. Напомним, что здесь 0  f(x)  1.

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

Например, пусть требуется вычислить интеграл S = f(x)dx,

где максимальное значение f(x) в отрезке [a,b] равно f_max.

Используя замену переменных x=a+(b-a)z и изменение масштаба по оси y, получим

S = = f_max(b-a)

где f^* (z)= f[a+(b-a)z]

Таким образом, рассматривается задача, сводимая к предыдущей.

9.3.5. Вторая модификация метода Монте-Карло для вычисления интегралов.

Кроме дискретных случайных величин рассматриваются случайные величины, возможные значения которых сплошь заполняют соответствующие интервалы на числовой оси. Такая случайная величина характеризуется функцией распределения, которая определяется так же, как и для дискретных случайных величин; при каждом x функция распределения

F(x)=P( < x)

Во многих случаях функция распределения F(x) оказывается дифференцируемой.

Если существует производная f(x)= , то случайная величина  называется непрерывной, а функция f(x) - функцией плотности вероятностей случайной величины .

Очевидно, что F(x)=

Рассмотрим случайную величину , возможные значения которой x принадлежат некоторому интервалу на оси 0Х. Закон распределения этой случайной величины задан функцией плотности f(x).

Вычислим вероятность попадания случайной величины внутрь полуинтервала [a,b], целиком лежащей в вышеуказанном интервале на оси ОХ.

Из теории вероятностей известно, что

Р(a   < b)=

Для решения этой задачи воспользуемся методом статистических испытаний.

Легко видеть, что в данном случае применима уже изложенная процедура вычисления интеграла вида

Из вышесказанного вытекает второй вариант метода статистических испытаний. Предположим, что в нашем распоряжении имеются случайные числа x_i, закон распределения которых описывается функцией плотности f(x).

Поскольку случайные числа x_i можно рассматривать как возможные значения случайной величины  , искомая вероятность Р(a   <b), приближенно равна (при достаточно большом числе испытаний N) частоте попадания случайных чисел x_i в полуинтервал a  x_i < b

Р 

где m - количество случайных чисел, удовлетворяющих неравенству

a  <b при N испытаниях. Другими словами, в условиях данной задачи

Процедура вычисления интеграла на ЭВМ следующая:

1. Из множества случайных чисел с законом распределения, заданным функцией плотности f(x), выбираем значения x_i.

2. Xi сравниваем с граничными значениями полуинтервала [a,b).

Выполнение неравенства ax_i<b отмечаем признаком =1, невыполнение =0;

т.е. = +1.

3. К содержимому счетчика количества испытаний добавляем единицу N=N+1.

4. Управление передается снова первой операции.

Проведя N таких испытаний, вычисляем значение частоты .

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

Так в первом варианте для реализации каждого испытания были необходимы два случайных числа Xi и Yi. Кроме того, требовалось вычисление подынтегральной функции для каждого значения Xi случайной величины , что, вообще говоря, может быть связано со значительным объемом вычислений. С другой стороны, первый способ является более общим случаем вычисления интеграла; этот способ пригоден для подынтегральных функций весьма широкого класса, тогда как второй способ предполагает, что подынтегральная функция является функцией плотности некоторой случайной величины, что накладывает на нее определенные ограничения. Первый способ оперирует с равномерно распределенной случайной величиной, которую легко получить на ЭВМ, тогда как рассмотренный здесь требует специального преобразования равномерного распределения в заданное. Перечисленные отличия разграничивают классы задач, которые удобно решать тем или другим способом. Так, первым способом широко пользуются для вычисления многократных интегралов, в то время как второй способ находит применение при моделировании сложных систем методом статистических испытаний.

9.3.6. Третья модификация метода Монте-Карло для вычисления интегралов.

Перейдем к рассмотрению третьей разновидности применения метода статистических испытаний для вычисления интегралов. Пусть  - случайная величина, возможные значения которой принадлежат отрезку [a,b],

a f_(X) - функция плотности этой случайной величины. Рассмотрим непрерывную функцию =g().Из теории вероятностей известно, что среднее значение или математическое ожидание функции g() равно:

М[ g()]=

Воспользуемся методом статистических испытаний для вычисления интегралов приведенного вида. Из совокупности возможных значений Х случайной величины  выберем последовательность Х₁, Х₂, ... , Х_N и вычислим значения g(Xi); i=1, 2, ... , N, причем Хi , не принадлежащие отрезку [a,b], в расчет не принимаются. При достаточно большом N среднее арифметическое весьма близко к математическому ожиданию случайной величины. Поэтому в качестве приближенного значения для интеграла приведенного вида может быть взято среднее арифметическое

М=

Рассмотрим практическую реализацию высказанного положения.

Пусть перед нами стоит задача вычислить

S =

Выберем некоторый закон распределения f(X), для которого имеется удобный способ получения случайных чисел, и областью определения которого является интервал (a,b). Преобразуем подынтегральную функцию следующим образом:

S = (1)

или S = , где g*(x)=

Последний интеграл вычисляется с помощью метода статистических испытаний. В частном случае, когда a и b конечны, или их можно считать конечными приближенно, в качестве f(X) целесообразно выбрать равномерный закон распределения.

Как известно, плотность вероятности равномерного закона распределения в интервале (a,b) равна: f(x)=

Подставим это значение в интеграл (1): S=(b-a)

и рассмотрим процедуру вычисления, описанную выше, для общего случая в связи с равенством: М=

Y

Из множества равномерно распределенных случайных чисел выбирается Xi. Для каждого значения Xi вычисляется U(Xi), затем вычисляется среднее значение: U(Xi)= функции U(X) в интервале (a,b). Таким образом, величина интеграла может быть представлена следующей формулой: S (b-a)U(Xi)

Из рисунка видно, что площадь под кривой заменяется площадью прямоугольника с основанием (b-a) и высотой, равной среднему значению функции U(Xi)на отрезке (a,b). Рассмотренный частный случай находит широкое применение для вычисления интегралов методом статистического моделирования в силу того, что границы области легко могут быть приведены к пределам интегрирования (a,b).

9.3.7. Кратные интегралы

Рассмотрим приемы вычисления интеграла методом статистических испытаний с точки зрения целесообразности их применения для вычисления многократных интегралов.

Рассмотрим интеграл S =

по ограниченной замкнутой области , заключенной внутри n-мерного параллелепипеда.

a_i  x_i  b_i; i=1, 2,...,n

Заменой переменных x_i=a_i+(b_i-a_i)y_i

приведем интеграл к виду:

S=(b₁-a₁)(b₂-a₂)...(b_n-a_n)S^*, где

S* =

Теперь область интегрирования целиком заключена внутри n-мерного куба. Для вычисления интеграла применим уже изложенный метод (первый). Изменением масштаба по оси z=f(x₁,x₂,...,x_n) преобразуем вычисляемый интеграл

S*=L+(M-L)S

где L - минимальное, а М - максимальное значение z по области 

S= и 0  z  1, где z=f(X₁,X₂,...,Xn)

Будем рассматривать вычисляемый интеграл как некоторый объем в

n- мерном пространстве.

S =

Координаты случайной точки в этом пространстве имеют вид совокупности n+1 случайных чисел R=(R₁,R₂,...,R_n) , равномерно распределенных в интервале (0,1).

Проведем N испытаний, состоящих в заполнении n-мерного единичного куба равномерно распределенными случайными точками. Будем проверять принадлежность этих точек объему V посредством неравенства:

V  f(R₁, R₂, ... , Rn)

По окончании эксперимента получим искомый объем

S 