Содержание
Рассматриваемые ранее методы называются Детерминированными, то есть лишенными элемента случайности.
Методы Монте–Карло (ММК) – это численные методы решения математических задач с помощью моделирования случайных величин. ММК позволяют успешно решать математические задачи, обусловленные вероятностными процессами. Более того, при решении задач, не связанных с какими-либо вероятностями, можно искусственно придумать вероятностную модель (и даже не одну), позволяющую решать эти задачи. Рассмотрим вычисление определенного интеграла
(2.20)
При вычислении этого интеграла по формуле прямоугольников интервал [A, B] разбиваем на N одинаковых интервалов, в серединах которых вычислялись значения подынтегральной функции. Вычисляя значения функции в случайных узлах, можно получить более точный результат:
(2.21)
(2.22)
Здесь γi — случайное число, равномерно распределенное на интервале
[0, 1]. Погрешность вычисления интеграла ММК
, что значительно больше, чем у ранее изученных детерминированных методов.
На рис. 2.5 представлена графическая реализация метода Монте-Карло вычисления однократного интеграла со случайными узлами (2.21) и (2.22).
 |
Рис. 2.5. Интегрирование методом Монте-Карло (1-й случай)
Однако при вычислении кратных интегралов детерминированными методами оценка погрешности перерастает в задачу порой более сложную, чем вычисление интеграла. В то же время погрешность вычисления кратных интегралов ММК слабо зависит от кратности и легко вычисляется в каждом конкретном случае практически без дополнительных затрат.
Рассмотрим еще один метод Монте-Карло на примере вычисления однократного интеграла:
 |
(2.23)
Рис. 2.6. Интегрирование методом Монте-Карло (2-й случай)
Как видно на рис. 2.6, интегральная кривая лежит в единичном квадрате, и если мы сумеем получать пары случайных чисел, равномерно распределенных на интервале [0, 1], то полученные значения (γ1, γ2) можно интерпретировать как координаты точки в единичном квадрате. Тогда, если этих пар чисел получено достаточно много, можно приблизительно считать, что 
. Здесь S – число пар точек, попавших под кривую, а N – общее число пар чисел.
Пример 2.1. Вычислить следующий интеграл:
Поставленная задача была решена различными методами. Полученные результаты сведены в табл. 2.1.
Решение определенного интеграла методом Монте-Карло
Метод Монте-Карло занимает особое положение среди методов вычисления определенных интегралов по двум причинам. Во-первых, это единственный метод, позволяющий вычислять интегралы высокой кратности. И, во-вторых, это метод, который дает лишь вероятностные гарантии степени точности вычисления интегралов.
Метод Монте-Карло — это статистический метод, его используют при вычислении сложных интегралов, решении систем алгебраических уравнений высокого порядка, моделировании поведения элементарных частиц, в теориях передачи информации и массового обслуживания, при исследовании сложных систем (экономических, биологических и т. д.).
Сущность метода состоит в том, что в решаемую задачу вводят случайную величину о, изменяющуюся по какому-то закону p(о). Как правило, случайную величину выбирают таким образом, чтобы искомая в задаче величина стала математическим ожиданием от случайной величины.
Таким образом, искомая величина определяется лишь теоретически. А вот чтобы найти ее численно, пользуются статистическими методами: берут выборку случайной величины о объемом элементов. В результате получают вариант случайной величины оi, для которых вычисляют их среднее арифметическое (выборочное среднее)
которое и принимают в качестве приближенного значения (оценки) искомой величины :
Для получения результата приемлемой точности по методу Монте-Карло требуется большое число статистических испытаний. Именно поэтому этот метод иногда так и называют: метод статистических испытаний.
Теория метода Монте-Карло изучает способы выбора случайных величин о для решения различных задач, а также способы уменьшения дисперсии используемых случайных величин. Уменьшение дисперсии играет большую роль, поскольку при равных объемах выборок, выборка с меньшей дисперсией имеет меньшую погрешность.
Итак, для вычисления однократного интеграла методом Монте-Карло может быть применена формула
где xi равномерно распределенное на интервале [a, b] случайное число. Справедлива следующая оценка точности вычисления интеграла по формуле (5.1.3) с вероятностью p=1-з выполняется неравенство
Например, если положить p=99%, тогда з = 0.01 и можно утверждать, что с вероятностью 99% справедливо неравенство
Все, что нужно для вычисления интегральной суммы по формуле (5.1.1) — это научиться получать случайные числа, равномерно распределенные на интервале [a, b]. Для этой цели можно использовать генератор случайных чисел, входящий в состав стандартных библиотек, поставляемых с компилятором. С помощью функции random легко получить случайное вещественное число, равномерно распределенное на интервале [0, 1] — например, результатом выполнения оператора x=random является случайное вещественное число из интервала [0, 1]. Имея случайное вещественное число из интервала [0, 1] легко получить случайное число из любого интервала. Например, если z — случайное число из интервала [0, 1],тогда x=a+(b-a)*z — случайное число из интервала [a, b].
Как видно из приведенных выше оценок погрешности формулы (5.1.3) точность вычисления интеграла и в методе Монте-Карло определяется числом слагаемых N в интегральной сумме — чем больше слагаемых, тем точнее результат. Ниже приведен пример и программа, вычисляющая определенный интеграл методом Монте-Карло.
Метод Монте-Карло легко обобщается на интегралы произвольной кратности. Например, двукратный интеграл может быть вычислен по формуле
где xi, yi — случайные числа, равномерно распределенные на интервалах [a, b] и [c, d] соответственно. Оценка точности вычисления интеграла по формуле (5.1.6) совершенно аналогично приведенной выше для случая однократного интеграла и поэтому здесь не приводится.
Основные понятия. Численные методы линейной алгебры и анализа
Навигация
Оглавление
Быстрый поиск
Метод прямоугольников для двойного интеграла¶
Рассмотрим вопросы численного нахождения двойного интеграла по прямоугольной области ([a,b] imes [c, d]) :
Рассмотрим два способа вывода аналога формулы прямоугольников для численного решения данного интеграла.
Вывод с помощью одномерных интегралов¶
Так как мы уже знаем формулы прямоугольников для интегралов от функций одной переменной, то удобно представить двойной интеграл через два интеграла, каждый из которых будет вычисляться от функции одной переменной и может быть численно найден с помощью уже известной нам одномерной формулы прямоугольников. С этой целью введем вспомогательную функцию (g(x)) :
Каждый из интегралов
можно вычислить с помощью формул численного интегрирования для одномерных интегралов. Воспользуемся формулой прямоугольников ( 3.1.5 ) и начнем с (g(x) = int_c^d f(x, y) dy) . На отрезке ([c, d]) введем равномерную сетку с шагом (h_y) :
Тогда для рассматриваемого интеграла составная формула прямоугольников будет иметь вид:
Для численного нахождения интеграла по переменной (x) на отрезке ([a, b]) введем равномерную сетку с шагом (h_y) :
на которой для интеграла (int_a^b g(x) dx) получим формулу прямоугольников:
Объединяя формулы по каждому направлению, получим составную формулу прямоугольников для двойного интеграла:
Прямой вывод¶
Формула (1) может быть получена непосредственно с использованием основной идеи построения одномерной формулы прямоугольников. Разобьем прямоугольник ([a, b] imes [c, d]) на (N_x imes N_y) ячеек (прямоугольников). Идея метода прямоугольников заключается приближении функции (f) кусочнопостоянной функцией: с постоянными значениями на ячейках равными значению (f) в центре ячейки. Ячейка ((i, j)) занимает область
центром которой является точка ((x_i, y_j) in omega_h = omega_
Программная реализация двойной суммы¶
Формула (1) содержит двойную сумму, которую обычно реализуют как двойной цикл for . Функция, реализующая формулу (1) на языке Python, может выглядеть следующим образом:
Сохранив эту функцию в файле-модуле rectangular_double.py , мы можем вычислить интеграл, например, (int_2^3int_0^2 (2x + y) dxdy) в интерактивной оболочке и показать, что функция возвращает правильное значение:
Повторное использование кода для одномерного интеграла¶
Вполне естественно реализовать метод прямоугольников для численного нахождения двойного интеграла так, как мы сделали в функции rectangular_double1 при условии, что мы знаем формулу (1). Однако, возникает вопрос: можем ли мы (как это делали при выводе формулы через формулу для одномерного интеграла) повторно использовать реализацию для одномерного интеграла при вычислении двойного? То есть, можем ли мы использовать функцию rectangular из раздела Реализация “дважды”? Ответ: “Да, можем”. Соответствующая функция будет очень короткой:
Важное преимущество такого подхода состоит в том, что мы используем уже протестированную функцию для одномерного метода прямоугольников и используем одномерный интеграл для вычисления двумерного в точности как и делали это математически.
Проверка с помощью тестовых функций¶
Как можно протестировать, что функция для двойного интеграла работает правильно? Лучший юнит тест — найти задачу, где бы погрешность была равна нулю, так как в этом случае мы будем точно знать численный ответ. Одномерная формула прямоугольников точна на линейных функциях не зависимо от количества частичных отрезков. Также и двумерная линейная функция (f(x,y) = px + qy + r) будет интегрироваться точно с помощью двумерной формулы прямоугольников. Можно выбрать (f(x,y) = 2 x + y) и создать соответствующую тестовую функцию, которая будет автоматически проверять наши две реализации двумерной формулы прямоугольников. Для вычисления интеграла от (f(x,y)) мы воспользуемся SymPy, чтобы избежать ошибок при ручном вычислении. Тестовая функция будет такой:
Должна ли тестовая функция что-то выводить? Если мы просто вызовем приведенную выше функцию и ничего не произойдет, то наша реализация корректна. Тем не менее, несколько раздражительно иметь функцию, которая полностью бесшумна, когда работает — мы можем быть уверены, что все правильно вычислено? Поэтому в процессе разработки настоятельно рекомендуется добавлять оператор печати так, чтобы мы могли контролировать расчеты и быть уверены в том, что тестовая функция делает то, что мы хотим. Так как в тестовой функции не должно быть никаких операторов печати, мы просто закомментировали его, как показано во фрагменте кода выше.
В качестве альтернативы методу прямоугольников можно выбрать формулу трапеций. Вывод формулы для двойного интеграла и реализация использует точно те же идеи, которые были приведены для формулы прямоугольников. Поэтому Упражнение 3.13: Вывести и реализовать формулу трапеций для двойного интеграла должно быть выполнено обязательно.
Формула прямоугольников для тройного интеграла¶
Теория¶
Так как нам удалось обобщить метод работающий для одномерного интеграла на двойной, то, естественно, довольно просто распространить его на три измерения. Рассмотрим тройной интеграл
Получим формулу прямоугольников для этого интеграла. Введем вспомогательные функции:
Для каждого из этих одномерных интегралов применим формулу прямоугольников:
Таким образом, получим прямоугольников для тройного интеграла:
Реализация¶
Аналогично случаю двойного интеграла, создадим файл rectangular_triple.py с реализацией соответствующих функций:
Метод Монте—Карло для численного интегрирования по областям сложной формы¶
Повторное использование формул для одномерного интегрирования при численном нахождении двойных и тройных интегралов возможно только, если область интегрирования прямоугольная (прямоугольник или параллелепипед). Для других форм областей нужно использовать совершенно другой подход. Наиболее общий способ для двумерных и трехмерных областей заключается в разбиении их на большое число маленьких треугольников или тетраэдров и использовании численного интегрирования для каждого треугольника или тетраэдра. Полный алгоритм и его реализация слишком сложны для данного спецкурса. Вместо этого мы будем использовать альтернативный, достаточно простой и общий метод, называемый методом интегрирования Монте—Карло (ММК). Код его реализации можно вместить на полстраницы, но эта реализация требует на порядок больше выполнений функций при вычислении двойных интегралов в сравнении с методом прямоугольников.
Однако ММК является более эффективным чем формула прямоугольников при вычислении многомерных интегралов (более трех переменных) по гиперкубическим областям. Представленные идеи для двойных и тройных интегралов могут быть легко обобщены на численное интегрирование по (m) переменным. В этом случае формула прямоугольников требует вычисления (m) сумм. При (N) частичных отрезках по каждой переменной формула прямоугольников потребует (N^m) вычислений функций. Это значит, что вычислительная работа возрастает как экспонента от числа размерностей. ММК не зависит таким образом от числа размерностей и является более предпочтительным при численном нахождении многомерных интегралов. Поэтому имеет смысл в главе, посвященной численному интегрированию, рассмотреть применение ММК как для сложных областей, так и для больших размерностей.
Алгоритм интегрирования по ММК¶
Идея ММК для численного интегрирования (int_a^b f(x) dx) заключается в использовании теоремы о среднем из математического анализа, в которой утверждается, что интеграл (int_a^b f(x) dx) равен произведению длины отрезка (здесь (b — a) ) и среднего значения (ar
Введем некоторые величины, которые позволят нам формализовать алгоритм численного интегрирования. Пусть дан двумерный интеграл:
где (Omega) — двумерная область заданная посредством вспомогательной функции (g(x,y)) :
Таким образом, граница области (partialOmega) задана неявной функцией (кривой) (g(x,y) = 0) . Такое описание областей распространено в последние десятилетия, при этом (g) называется функцией уровня, а граница (g = 0) — нулевым контуром функции уровня. Для простых областей можно легко построить функцию (g) вручную, но в более сложных промышленных приложениях следует обратиться к математическим моделям построения (g) .
Пусть (A(Omega)) — площадь области (Omega) . Мы можем численно найти интеграл по следующему ММК:
- помещаем область (Omega) внутрь прямоугольника (R) ;
- генерируем большое число случайных точек на (R) ;
- вычисляем долю (q) точек, которые попали в область (Omega) ;
- приближаем (A(Omega)/A(R)) числом q , т.е., полагаем (A(Omega) = q A(R)) ;
- вычисляем среднее значение (ar
) функции (f) на области (Omega) ; - вычисляем приближенное значение интеграла как (A(Omega)ar
) .
Отметим, что площадь (A(R)) прямоугольника (R) легко вычислить, при том что площадь (A(Omega)) нам не известна. Однако, если предположить, что доля площади (A(R)) занимаемой областью (Omega) такая же как доля случайных точек, попавших внутрь (Omega) , можно получить простое приближение для (A(Omega)) .
Для того, чтобы получить представление о методе, рассмотрим область, ограниченную окружностью и помещенную в прямоугольник, как показано на рисунке ниже. Набор случайных точек изображен синим цветом.
Реализация¶
Функция на Python, реализующая численное интегрирование (int_Omega f(x, y) dxdy) , может выглядеть следующим образом (файл mc_double.py ):
Проверка¶
Возьмем простейший тестовый пример: найти площадь прямоугольника ([0, 2] imes [3, 4.5]) , помещенного в прямоугольник ([0, 3] imes [2, 5]) . Точное значение площади — 3, но ММК, к сожалению, никогда не является точным, поэтому не возможно предсказать результат алгоритма. Все, что мы знаем — это то, что приближенное значение интеграла должно стремиться к 3 при числе случайных точек стремящемся к бесконечности. Также при фиксированном числе точек мы можем запускать алгоритм несколько раз и получать разные значения, которые колеблются около точного значения, так как разные выборки точек используются при разных вызовах алгоритма.
Площадь прямоугольника можно вычислить с помощью интеграла (int_0^2int_3^ <4.5>dxdy) , поэтому в этом случае положим (f(x,y) = 1) , а функцию (g) можно определить, например, следующим образом: 1 если точка ((x, y)) лежит внутри прямоугольника ([0, 2] imes [3, 4.5]) иначе -1. Ниже приведен пример использования функции MonteCarlo_double в интерактивной сессии IPython при вычислении площади с разными мощностями выборки:
Мы видим, что значения колеблются около 3, что подтверждает правильную реализацию, но в принципе, за неточным ответом могут скрываться ошибки.
Известно, что стандартное отклонение приближенного значения интеграла сходится как (n^<-1/2>) , где (n) — мощность выборки (число случайных точек). Такой вид оценки скорости сходимости может использоваться при тестировании реализации метода, но мы этот способ рассматривать не будем.
Тестовая функция для функций со случайными числами¶
Чтобы написать тестовую функцию, нам нужен юнит тест, который имеет идентичное поведение при разных запусках. Это представляется сложным при использовании случайных чисел, так как эти числа меняются каждый раз при запуске алгоритма, и каждый запуск дает разные результаты. Стандартный способ тестирования алгоритмов со случайными числами заключается в фиксировании начального числа генератора случайных чисел. Тогда последовательность чисел будет одной и той же при каждом запуске алгоритма. Предположим, что функция MonteCarlo_double работает. Мы фиксируем начальную точку, наблюдаем некторый результат, и берем этот результат в качестве правильного. При условии, что тестовая функция всегда использует эту начальную точку, мы должны получать в точности тот же результат каждый раз, когда вызывается функция MonteCarlo_double . Наша тестовая функция может выглядеть так (файл mc_double.py ):
Интеграл по кругу¶
Тест представленный выше включал простейшую функцию (f(x,y) = 1) . Нам следует выполнить также тесты для не постоянной функции (f) и для более сложной области. Пусть (Omega) — круг с центром в начале координат и радиусом 2, и пусть (f(x,y) = sqrt) . Такой выбор позволяет вычислить точное значение интеграла: в полярных координатах (int_Omega f(x,y) dxdy) преобразуется в (2piint_0^2 r^2dr = 16pi/3) . Мы должны быть готовы к достаточно грубым приближениям, которые колеблются около этого результата. Мы знаем по опыту предыдущего тестового примера, что более точные результаты можно получить на большем числе случайных точек. Ниже дается пример тестовой функции: