Адаптивные дискретно-стохастические алгоритмы численного интегрирования

0.1. Детерминированные и стохастические кубатурные формулы. С развитием вычислительной техники возрастает интерес к численным методам решения прикладных задач. Математическое описание исследуемого процесса сводится, как правило, к рассмотрению неизвестной функции многих переменных, для которой записывается система дифференциальных уравнений (см., например, [1]). Одним из способов приближенного решения этой системы является построение разностной схемы и сведение задачи к решению системы линейных уравнений для значений неизвестной функции в узлах сетки (см., например, [2−4]). Часто также оказывается целесообразным сведение (или постановка) задачи к интегральной форме, когда исследуемая функция представляет собой многократный интеграл, зависящий от параметра, или является решением интегрального уравнения (см., например, [1]). В последнем случае при возникновении интегрального уравнения Фредгольма второго рода соответствующий интегральный оператор предполагается сжимающим, и решение записывается в виде ряда Неймана, представляющего собой сумму параметрических интегралов бесконечно возрастающей кратности (см., например, [5]). Далее используются численные методы приближенного вычисления получаемых многократных интегралов на ЭВМ. При разработке этих методов решается

ЗАДАЧА 0.1. Построить алгоритм вычисления интеграла здесь X — замыкание тех х € Д', для которых <�у (х) ф 0.

Для интегралов I малых размерностей I с гладкими (в обычном или обобщенном смыслах) подынтегральными функциями д и относительно простыми областями интегрирования X развита теория квадратурных (для случая I = 1) и кубатурных (для I > 1) формул (см., например, [4, 6]). Кубатурная формула в общем случае имеет вид где {х1,., хп} - заданные детерминированные (и, как правило, регулярные) узлы сетки в В1, а {с1,., с71} - веса. Вычисление интеграла (0.1) по формуле (0.2) будем называть АЛГОРИТМОМ 0.1. В качестве X в данной работе будем использовать, как правило, простые компактные подмножества в В} (чаще всего — /-мерный единичный куб = {х = (жь. .. , ?-): 0 < Хг < 1- 2 = 1,., /}).

Оптимальный выбор узлов и весов связан с минимизацией погрешности 5п = 1 — 5П| и основан (явно или неявно) на использовании аппроксимаций подынтегральной функции д. Главным достоинством кубатурных формул является возможность получения гарантированной и сравнительно быстрой сходимости 5п к нулю при п —У оо для классов гладких подынтегральных функций д.

К недостаткам «классических» (детерминированных) кубатурных формул на классах подынтегральных функций следует отнести:

— слабый учет специфики той или иной подынтегральной функции;

— необходимость разработки специальных численных алгоритмов поиска оптимальных весов и (или) узлов;

— чувствительность к росту размерности и к гладкости начальных данных (подынтегральной функции д и области X);

0.1) п

0.2)

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

Для существенно многомерных задач 0.1 (т.е. для случая ¿->1и даже I —V сю) достаточно эффективным оказывается стандартный метод Монте-Карло (см., например, [7−9]). Этот алгоритм основан на представлении

1 = 1 д (х)Дх) ¿-х = ЕС, С = <7®, д (х) = </(х)//(х). (0.3)

Для выполнения (0.3) достаточно потребовать /(х) ф 0 при х 6 X. Весовая функция / является вероятностной плотностью в Я1, т. е. /(х) > 0 и //(х)сЬс = 1, а /-мерный случайный вектор? распределен согласно плотности /.

АЛГОРИТМ 0.2 [7−9]. 1. Реализуем п значений случайного вектора. , £п.

2. Вычисляем приближенное значение интеграла (0.1):

3=1 3=1

Если в (0.4) трактовать как набор независимых одинаково распределенных случайных векторов с плотностью распределения /, то случайные величины (1 = 9(^1), -■-■-)(п ~ я{€п) будут также независимыми одинаково распределенными с математическим ожиданием I (см. соотношение (0.3)) и дисперсией

БС = а2 = Бд (0 — I ?2(х) Дх) ¿-х — /2 = I ¿-х ~ 1 (0.5)

Если величина (0.5) конечна, то в силу закона больших чисел (см., например, [10]) формула (0.4) верна для достаточно больших п.

Сразу заметим, что для /(х) = 1 (т. е. для равномерного распределения = формула (0.4) является частным случаем формулы (0.2) для с =. — сп = 1/п и ^ = Таким образом, разница между алгоритмами 0.1 и 0.2 с постоянными весами связана с выбором узлов — детерминированным или стохастическим. Отметим также, что и для «детерминированного» алгоритма 0.1 введение весовой функции типа / и связанного с ней набора узлов позволяет улучшать качество кубатурных формул (0.2).

Алгоритм 0.2 имеет следующие положительные свойства:

— возможность уменьшения трудоемкости алгоритма за счет, уда.'чного (согласованного с видом подынтегральной функции д) выбора весовой функции / (алгоритм выборки по важности) или преобразования исходного интеграла (выделение главной части) и весовой функции (методы математического ожидания и расщепления, выборка по группам и т. п.) — см. далее подразд. 0.5;

— веса имеют простой вид, а узлы реализуются согласно выбираемому вероятностному распределению с плотностью /;

— относительно слабая чувствительность к росту размерности I и к гладкости начальных данных (подынтегральной функции д и области X);

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

Последнее свойство обусловлено тем, что затраты 5 алгоритма 0.2 можно подсчитать по простой формуле 5 = пЬ, где t — среднее время получения одного выборочного значения случайной величины С (это время, в свою очередь, связано со «сложностью» вычисления функции д и со средними затратами на реализацию одного выборочного значения вектора ?). Далее, в силу центральной предельной теоремы (см., например, [10]) для достаточно больших п имеет место соотношение

При фиксированном уровне погрешности число случайных узлов п прямо пропорционально дисперсии а2 случайной величины и вместо я в качестве величины, отражающей затраты стандартного метода Монте-Карло, можно ввести число которое называется трудоемкостью алгоритма 0.2 [7−9]. Рассмотрение величины (0.7) вместо 5 является более удобным при оптимальном выборе весовой функции /, т. к. величина 51 в явном виде содержит вероятностную характеристику (дисперсию) случайной величины С, которая, в свою очередь, определяется через /. Неизвестное значение (0.5) можно приближенно вычислить по набору выборочных значений. [7−9]:

1 V 1 (у^У (ой а8)

Формула (0.6) отражает также и главный недостаток алгоритма 0.2 — относительно низкую (порядка 1/у/п) скорость сходимости погрешности к нулю при возрастании числа узлов п. К примеру, для / = 1 и д € С2(Х) простейшая формула прямоугольников имеет порядок погрешности 1/п2 (см., например, [4]). Это обуславливает использование алгоритма 0.2 только для достаточно больших размерностей I.

0.2. Дискретно-стохастические методы численного интегрирования. Формулы Бахвалова и теория сложности. Эффективными могут оказаться и смешанные, комбинированные процедуры численного интегрирования, сочетающие в себе элементы алгоритмов 0.1 и 0.2. Как правило, эти схемы содержат «детерминированную» составляющую, связанную с регулярной дискретизацией области X, а также «стохастическую» составляющую, связанную с применением метода Монте-Карло. Поэтому, вслед за А. В. Войтишеком [11], мы будем называть такие алгоритмы дискретно-стохастическими.

По-видимому, одна из первых численных схем подобного рода была представлена в работе Н. С. Бахвалова [12] (см. также [4] и подраздел 1.2.4 данной работы). Исходной областью X являлся ¿—мерный единичный куб (?1, который разбивался на п равных кубов с вершинами в точках равномерной сетки. В каждом ^'-ом элементарном кубе выбирался узел кубатурной формулы (0.2) случайным образом (согласно равномерному распределению). Представленный алгоритм можно считать как кубатурной формулой (0.2) со случайными узлами (такой термин для подобных конструкций имеется, например, в [8]), так и предельным случаем выборки по группам в методе Монте-Карло [7−9]. Н. С. Бахвалов показал в [12], что его алгоритм является оптимальным (по скорости сходимости к нулю погрешности 5п при п —> со) в пространстве непрерывно дифференцируемых подынтегральных функций С1(Х). Для этого ему потребовалось получать оценки сверху и снизу для 5п.

Методология работы [12] явилась основой для интенсивного развития так называемой теории сложности (см. [13] и сопутствующие этой монографии работы). Эту теорию можно назвать «философией вычислительных алгоритмов». В ней изучается вопрос о том, каков максимальный порядок стремления к нулю погрешности 5й для

0.6)

0.7) данного класса вычислительных алгоритмов с количеством вычислительных операций п.

Проблема численного интегрирования (задача 0.1) является наиболее удобной и часто используемой иллюстрацией в теории сложности. В работах по теории сложности рассмотрены вопросы оптимальности кубатурных формул (0.2) для различных пространств В (Х) подынтегральных функций. При построении соответствующих нижних границ для 5п требуется строить конкретную численную схему типа (0.2). В некоторых (достаточно редких) случаях эти алгоритмы имеют вид, вполне пригодный для практических вычислений. Например, для В (Х) = С2(Х) в работе [12] Н. С. Бахвалов построил оптимальный алгоритм, в котором, в дополнение к алгоритму для д € С1{Х), кроме узла х^ выбирается точка х7, симметричная х^ относительно центра ?-го элементарного куба (см. подраздел 1.2.4 настоящей работы).

С точки зрения теории сложности алгоритм 0.2 решения задачи 0.1 является не слишком интересным, так как, в силу соотношения (0.6), имеет неулучшаемый порядок сходимости 1 /у/п. Поэтому в связи с задачей 0.1 в теории сложности распространено обращение к случаю использования так называемых квази-случайных узлов кубатур-ной формулы (0.2) [13]. Здесь теория кубатурных формул смыкается со специальным нетривиальным разделом теории чисел. Наиболее эффективный алгоритм построения квази-случайных узлов представлен в [7] (это ЛПГ-последовательность Соболя). Вопросы использования квази-случайных чисел в численном интегрировании в данной работе не рассматриваются.

В общей теории сложности (в том числе и для задачи 0.1) принимается ряд допущений, которые не всегда выполняются на практике. Например, приравниваются затраты на вычисление значения подынтегральной функции д в точке и одно арифметическое действие. В связи с этим открытым остается вопрос, будут ли теоретически оптимальные схемы иметь на практике минимальную трудоемкость. Забегая вперед, отметим, что в представляемой рабоге указаны примеры ситуаций, когда оптимальные (с точки зрения теории сложности) кубатурные формулы не являются наилучшими на практике (см. раздел 1.2).

0.3. Цель и структура диссертации. Как уже было сказано выше, возможность эффективной практической реализации оптимальных кубатурных формул является большой редкостью (особенно для многомерных случаев). Поэтому вполне разумными видятся подходы, связанные с построением и исследованием эффективно реализуемых дискретно-стохастических численных алгоритмов, дающих относительно небольшие значения трудоемкости (0.7).

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

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

Заключение

В работе получены следующие результаты.

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

• Разработан двусторонний геометрический метод Монте-Карло. Исследованы вопросы оптимизации этого метода для кусочно-постоянных мажорант и минорант подынтегральной функции.

• Исследована эффективность дискретно-стохастических состоятельных оценок метода Монте-Карло: взвешенной равномерной выборки и оценки с поправочным множителем. Получены верхние оценки дисперсий и оценены затраты соответствующих дискретно-стохастических численных схем.

• Проведено сравнение численных алгоритмов аппроксимации решения интегрального уравнения второго рода на основе рандомизации конечных и бесконечных отрезков ряда Неймана. Исследованы возможности использования отрезков однородных цепей Маркова конечной и случайной длицы, а также эффективных дискретно-стохастических метод<�эв численного интегрирования. Представлецы примеры, подтверждающие целесообразность использования смещенных детермини-рованно-стохастических оценок решения.

• Проведен сравнительный анализ известных аппроксимационных базисов с точки зрения использования их в алгоритмах численного статистического моделирования. Показана целесообразность использования приближений Стренга-Фикса и Бернштейна в дискретно-стохастических численных схемах.

• Рассмотрены возможности применения траекторий спектральных (гауссовских и негауссовских) моделей случайных функций при тестировании численных алгоритмов интегрирования. Такой подход позволил добиться независимости тестирования, получить требуемые свойства подынтегральных функций (гладкость, «сложность» вычисления и др.), вывести аналитические выражения для средних погрешностей. Исследован вопрос о слабой сходимости используемых негауссовских численных моделей.

Проведенные численные эксперименты (в том числе, с использованием стохастической тестовой системы функций) показали, что разрабатываемые здесь дискретно-стохастические алгоритмы эффективны, как правило, для задач «умеренно большой» кратности (конкретнее, для размерностей от трех до десяти).