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

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

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

(9.1)

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

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

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

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

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

(9.2)

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

.

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

. (9.3)

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

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

. (9.4)

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

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

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

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

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

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

.

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

??

Рис. 9.1

Рассмотрим две случайные величины ξ — равномерно распределенную на отрезке [a, b] и η — равномерно распределенную на отрезке [0, f_max], где .

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

(9.5)

Проведем серию из N испытаний и получим N случайных точек , принадлежащих прямоугольной области . Подсчитаем число N_f точек, удовлетворяющих условию . Тогда вероятность попадания случайной точки (ξ, η) в криволинейную трапецию приближенно равна относительной частоте попадания в трапецию, т.е. и интеграл приближенно вычисляется по формуле:

(9.6)

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

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

(9.7)

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

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

(9.8)

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

Пример 9.1. Вычислить интеграл

.

Решение в программе Excel. Точное значение интеграла составляет .

В данном примере возьмем N = 100 и

Введем формулы, как показано в табл. 9.1.

Таблица 9.1.

Выделим ячейки A2:C2 и маркером заполнения протянем их вниз до строки 101 включительно. В ячейки C102, C103 соответственно введем формулы

=СУММ(C2:C101), =(C102/100)*(3-1)*9.

В ячейке C103 получим приближенное значение интеграла. Это значение будет изменяться при каждом пересчете листа, так как используется функция СЛЧИС(), которая автоматически дает новое случайное число из интервала при каждом пересчете. Приведем некоторые из полученных значений: 7,32; 9,72; 9,18.

Создадим макрос — функцию для вычисления интеграла по формуле (9.6). С помощью команды меню «Сервис — Макросы — Редактор Visual Basic» откроем окно редактора, выполним команду «Insert — Module» и введем тексты описаний функций:

Function f(ByVal x)

f = x ^ 2

End Function

Function int_st(ByVal a, ByVal b, ByVal Fmax, ByVal N)

Randomize

Nf = 0

For i = 1 To N

x = Rnd * (b - a) + a

y = Rnd * Fmax

If y < f(x) Then Nf = Nf + 1

Next i

int_st = Nf * (b - a) * Fmax / N

End Function

Function int_MK1(ByVal a, ByVal b, ByVal N)

Randomize

S = 0

For i = 1 To N

x = Rnd * (b - a) + a

S = S + f(x)

Next i

int_st1 = (b - a) * S / N

End Function

Здесь функция int_st() реализует формулу (9.6), а функция int_MK1() — формулу (9.7). Обратим внимание на то, что датчик случайных чисел Rnd языка Visual Basic дает случайные числа, равномерно распределенные в интервале , которые с помощью преобразования отображаются на отрезок [a, b].

Введем в ячейку D103 формулу =int_st(1;3;9;100).

При каждом пересчете формулы в ячейке D103 получим значения, имеющие погрешность того же порядка, что и значения в ячейке C103. Изменим в формуле D103 число испытаний N = 100 на N = 1000000: =int_st(1;3;9;1000000). Теперь при каждом пересчете будем получать существенно более точные значения интеграла.

Ниже в таблице 9.2 приведены результаты вычисления интеграла по формуле (9.6) при различных значениях числа испытаний N. Использование формулы (9.7) дает таблицу (9.3).

Таблица 9.2

Таблица 9.3

Таблицы 9.2, 9.3 наглядно иллюстрируют оценку порядка погрешности: чтобы уточнить одну десятичную цифру (т.е. уменьшить погрешность в 10 раз), надо увеличить число испытаний в 100 раз!