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

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

Вопросам разработки, исследования, тестирования и применения численных методов приближенного вычисления определенных интегралов посвящено огромное количество литературы, включая обширные монографии [26], [32], [41], [43], [65], справочники [27], [66] и множество журнальных статей. Эти вопросы в той или иной мере отражены и в большинстве общих курсов численного анализа [5], [6], [14], [29].

Тем не менее, достаточно полно разработанной можно считать лишь теорию методов одномерного интегрирования (т.е. интегрирования по подмножествам вещественной оси М1). Для таких задач построено и исследовано множество квадратурных формул и других способов интегрирования (методы с автоматическим выбором шага, статистические методы). Эти способы различаются по своей сложности, эффективности и области приложения, однако для подавляющего большинства встречающихся на практике одномерных задач может быть выбран подходящий метод, позволяющий произвести интегрирование с заданной точностью при разумных затратах машинного времени.

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

Большинство известных в настоящее время методов вычисления многомерных интегралов можно отнести к одной из трех больших групп — кубатурные формулы, теоретико-числовые методы (методы квази-Монте-Карло) и методы Монте-Карло.

Кубатурные формулы представляют приближения к интегралу в виде взвешенной суммы значений подынтегральной функции в некоторых точках области интегрирования. Узлы и веса квадратурных формул подбираются так, чтобы обеспечить точное интегрирование функций определенного класса (как правило, алгебраических или тригонометрических полиномов степени не выше заданной). Основы теории кубатурных формул изложены в монографиях С. Л. Соболева [34], И. П. Мысовских [31], А. Страуда (А.Н^гоисГ) [65].

Наиболее простой и до определенного момента наиболее распространенный способ построения кубатурных формул состоит в сведении кратного интеграла к повторному и последовательному применению одномерных квадратурных формул по каждой из переменных. Известно, что на классе Ст функций з переменных, имеющих в области интегрирования все ограниченные частные производные до порядка т включительно, погрешность такого способа интегрирования имеет порядок О (п~тгде п — число точек области, в которых вычисляются значения подынтегральной функции или ее производных. Видно, что с ростом размерности интеграла «$ указанный метод становится практически неприменимым из-за слишком медленной сходимости. Оказывается, однако, что на данном классе функций порядок убывания погрешности невозможно улучшить — для любого детерминированного метода интегрирования погрешность на классе Ст не может убывать быстрее, чем 0(п~Ш//5). Впервые это было установлено Н. С. Бахваловым [3], [4].

Анализ причин подобной ситуации, привел к появлению следующего соображения: хотя указанный класс подынтегральных функций и является весьма важным в теоретическом плане, он не вполне адекватно описывает реальные классы функций, встречающиеся в многочисленных приложениях (например, в математической физике), а для таких классов положение может оказаться несколько иным. Так, Н. М. Коробовым [24] были введены в рассмотрение классы функций Щ и Е", в определенном смысле близкие к рассмотреному выше классу Ст, оценка снизу погрешности интегрирования на которых имеет порядок 0{п" а) при, а > 1, т. е. не зависит от размерности. Используя методы теории чисел, Н. М. Коробов построил методы интегрирования на единичном кубе (а точнее, сетки специального вида, так называемые параллелепипедальные сетки), на которых погрешность интегрирования отличается от оценки снизу лишь множителем вида In7 п, где 7 не зависит от п. Построение параллелепипедальных сеток проводится с помощью набора специальным образом выбираемых целых чисел, называнных оптимальными коэффициентами. Результаты подобного же характера были получены И. М. Соболем [35] для классов Sp функций с быстро сходящимися рядами Хаара. И. М. Соболю удалось построить сетки и последовательности точек единичного куба, обеспечивающие почти наилучший порядок убывания погрешности на этих классах, равный где р > 1 не зависит от размерности.

Все веса квадратурных формул Коробова и Соболя совпадают и равны гГ1, а сетки узлов обладают медленно растущими с ростом п значениями отклонения (discrepancy) — функции, характеризующей степень равномерности расположения узлов в единичном кубе. Таким образом, эти методы оказываются тесным образом связанными с теорией равномерного распределения (в смысле Г. Вейля). Равномерно распределенные последовательности являются детерминированным аналогом последовательностей статистически независимых реализаций равномерно распределенной случайной величины и могут (при определенных условиях) подменять последние в задачах статистического моделирования — поэтому их часто называют квазислучайными. Среднее арифметическое значений функции в первых п узлах такой последовательности при п —> оо стремится к интегралу от нее (для любой функции, интегрируемой по Риману), что позволяет использовать их в задачах численного интегрирования. Методы, основанные на использовании равномерно распределенных последовательностей, а также сеток с медленно растущими отклонениями, образуют вторую группу способов многомерного интегрирования, которые называют методами квази-Монте-Карло (или теоретико-числовыми методами). Основными в теории методов квази-Монте-Карло являются вопросы построения последовательностей, отклонения которых возрастают не быстрее, чем некоторая степень.

In77 (low-discrepancy sequences). Их использование позволяет для многих классов подынтегральных функций (например, функции с ограниченной полной вариацией в смысле Хардп) получить скорость убывания погрешности интегрирования порядка 0(n~1+d), где 6 > 0 сколь угодно мало, вне зависимости от размерности.

Первый класс последовательностей с медленно растущими отклонениями был построен в I960 году Холтоном (J.H.Halton) [46]. В 1966 году И. М. Соболем [3−5] был построен новый класс равномерно распределенных последовательностей (ЛПг-последовательности), обладающих очень хорошими свойствами равномерности и быстрым алгоритмом расчета. Обобщение ЛПг-последовательностей производилось Форе (H.Faure) [44] и Нидер-рейтером (H.Niederreiter) [-59]. Также имеется целый ряд работ (например, [60], [64]), посвященных построению теоретико-числовых сеток с малыми значения отклонения — так называемые решеточные формулы интегрирования (lattice rules).

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

Другой путь, по которому удается преодолеть резкое ухудшение сходимости методов интегрирования с ростом размерности задачи, состоит в отказе от получения гарантированной оценки погрешности метода и переходе к вероятностным оценкам, т. е. таким, которые имеют место не обязательно, а лишь с некоторым наперед заданным уровнем достоверности. Н. С. Бахвалов, рассматривая в работах [3], [4] различные недетерминированные способы интегрирования, показал, что для очень многих классов подынтегральных функций оценка снизу вероятностной погрешности таких способов оказывается по порядку лучше, нежели оценка снизу гарантированной погрешности всевозможных детерминированных методов, причем это улучшение не зависит от размерности. Методы, использующие такую трактовку сходимости и имеющие вероятностную основу, образуют третью большую группу способов многомерного интегрирования, объединяемых под общим названием «метод Монте-Карло» (метод статистических испытаний). Эти методы также начали активно развиваться с внедрением вычислительных машин (их возникновение относят к 1949 году) и область их применения отнюдь не ограничивается задачами вычисления интегралов. Среди книг, посвященных теории и различным приложениям метода Монте-Карло необходимо назвать монографии С. М. Ермакова и Г. А. Михайлова [18], [21], [30], И. М. Соболя [36], Дж. Хэммерсли (J.M.Hammersley) и Д. Хэндскомба (D.C.Handscomb) [51], а также Дж. Спанье (J.Spanier) и Э. Гелбарда (Е.М.Gelbard) [38].

Применительно к задаче вычисления интегралов метод Монте-Карло обнаруживает ряд несомненных достоинств. Это, в первую очередь, чрезвычайная простота его реализации, последовательный характер исполнения, а также возможность контролировать погрешность вычислений апостери-орно, во время счета. Немаловажным является и то, что даже простейший метод Монте-Карло сходится (с любым наперед заданным уровнем достоверности) для всех функций из L-2 со скоростью 0{п~1/'2)1 не зависящей от размерности вычисляемого интеграла. Применение же специальных приемов уменьшения дисперсии (variance reduction methods) позволяет в некоторых случаях значительно улучшить сходимость метода Монте-Карло.

Подавляющее большинство работ, касающихся интегрирования методом Монте-Карло, посвящено именно вопросам улучшения его сходимости, поскольку асимптотическая скорость убывания погрешности 0(п~1//2) является довольно низкой. В монографии С. М. Ермакова [18] дано подробное изложение основных приемов понижения дисперсии оценки интеграла. Одним из наиболее общих приемов такого рода является использование случайных квадратурных формул. Случайные интерполяционно-квадратурные формулы были введены в рассмотрение С. М. Ермаковым и В. Г. Золотухиным [19] и представляют собой обобщение детерминированных интерполяционно-квадратурных формул, в которых узлы представляют собой случайные величины с заданным законом распределения, а веса являются некоторыми функциями узлов. Дальнейшая разработка случайных квадратурных формул проводилась в работах С. М. Ермакова и Б. Л. Грановского [10]-[12], [16], [17]. В частности ими было введено понятие допустимости (неулучшаемости) случайной квадратуры на классе подынтегральных функций, получен ряд критериев допустимости и построены конкретные допустимые случайные квадратуры. Примеры использования случайных квадратурных формул, приведенные в [18], показывают их значительное превосходство перед простейшим методом Монте-Карло.

Одним из наиболее распространенных приемов уменьшения дисперсии при интегрировании методом Монте-Карло является метод существенной выборки или выборки по важности (importance sampling). Этот метод был предложен Каном (H.Kalm) [56], [57] и суть его заключается в том, чтобы выбирать большее количество случайных точек в существенных областях, т. е. там, где подынтегральная функция принимает большие по модулю значения. Математическим основанием такого подхода служит тот факт, что плотность распределения случайных точек, доставляющая минимум функционалу дисперсии метода Монте-Карло, пропорциональна модулю подынтегральной функции. Более того, в случае знакопостоянной подынтегральной функции, эта плотность обеспечивает нулевое значение дисперсии оценки. Такая плотность не может непосредственно использоваться в расчетах, поскольку поиск коэффициента пропорциональности не менее сложен, чем исходная задача, однако она может быть достаточно хорошо приближена, если имеется аппроксимация модуля подынтегральной функции, интеграл от которой легко вычисляется.

Весьма близким к методу существенной выборки является метод выделения главной части (в англоязычной литературе его называют correlated или control variate sampling). Он также основан на построении легко интегрируемой аналитически аппроксимации подынтегральной функции, однако, в отличие от метода существенной выборки, эта аппроксимация используется непосредственно — интеграл от нее выбирается в качестве первого приближения к искомому интегралу. После этого метод Монте-Карло применяется к вычислению интеграла лишь от погрешности этой аппроксимации, что позволяет значительно сократить дисперсию. Достоинством такого метода по сравнению с методом существенной выборки является отсутствие необходимости моделирования случайных величин, имеющих заданное «хорошее» распределение. Однако и для его эффективного использования необходимо иметь достаточно хорошую, легко интегрируемую аппроксимацию подынтегральной функции. Некоторые авторы называют подобные аппроксимации «легкими» (easy) функциями.

Различные варианты методов существенной выборки и выделения главной части определяются способами построения «легких» аппроксимаций. Так, например. Розенберг (L.Rosenberg) [62] предложил использовать в качестве аппроксимации (на единичном кубе) полином Бернштейна подынтегральной функции. Использование полинома Сю (L.C.Hsu), близкого к полиному Бернштейна, исследовалось в работе А. Бабаева [2]. В статье А. В. Войтишека и Д. Р. Колюхина [8] предложено приближать подынтегральную функцию с помощью аппроксимации Стренга-Фикса, а в работе.

Блага (P.Blaga) [40] — с помощью B-сплайна. Во всех указанных случаях построение аппроксимаций проводится по значениям подынтегральной функции в узлах некоторой регулярной сетки, что несколько снижает практическую ценность предложенных подходов, поскольку число узлов регулярных сеток экспоненциально растет с ростом размерности. С вычислительной точки зрения использование таких аппроксимаций в методах существенной выборки и выделения главной части можно трактовать как статистическое уточнение соответствующих этим аппроксимациям интерполяционных кубатурных формул. В работе Сасаки (T.Sasaki) [G3] предлагается приближенно факторизовать подынтегральную функцию s переменных в произведение (s — 1) функций, каждая из которых зависит только от двух переменных. В работе Ю. В. Воронцова [9] аппроксимация строится в виде отрезка разложения подынтегральной функции в ряд Фурье по какой-либо ортонормированной системе. В качестве такой системы предложено использовать, например, систему функций Хаара. Применение подобного разложения по полиномам Эрмита рассмотрено в работе Хитцл (D.L.Hitzl) и Зеле (F.Zele) [54].

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

Метод симметризации подынтегральной функции (в англоязычной литературе antithetic variates method) берет свое начало в работах Хэммерсли (J.M.Hammersley) и Мортона (K.W.Morton) [52], [58]. Суть метода состоит в замене подынтегральной функции / на некоторую функцию /о, интеграл от которой равен исходному и которая обладает лучшими чем / свойствами с точки зрения применения метода Монте-Карло к ее интегрированию, т. е. обеспечивает меньшее значение дисперсии. Большое число способов симметризации, применимых к достаточно гладким одномерным функциям предложено Холтоном и Хэндскомбом (J.H.Halton, D. C'.Handscomb) [50]. В определенных случаях использование метода симметризации приводит к ускорению сходимости метода Монте-Карло (т.е. увеличению порядка убывания дисперсии). Необходимо заметить, однако, что при переходе в многомерный случай методы симметризации становятся весьма громоздкими, что усложняет их практическое использование.

Метод выборки по группам или расслоенной выборки (stratified sampling) состоит в разбиении области интегрирования на ряд подобластей, интеграл по каждой из которых вычисляется простейшим методом Монте-Карло с числом разыгрываемых точек, зависящим от подобласти. На основе метода выборки по группам также можно построить методы интегрирования с повышенной скоростью сходимости. Метод такого типа был впервые предложен и обоснован Дупачем (V.Dupac) [42], дальнейшее исследование и обобщение проводилось Хабером (S.Haber) [45]. В частности, путем объединения метода выборки по группам с методом симметризации, Хаберу удалось построить метод интегрирования, имеющий на классе дваждыдиф-ференцируемых функций С2 вероятную ошибку с порядком 0(n~1//2~2//s), оптимальным на данном классе. Это ускорение сходимости, однако, оборачивается потерей одного из важнейших достоинств метода Монте-Карло — последовательного характера исполнения, т.к. при изменении п (в случае, когда требуемая точность еще не достигнута) требуется заново пересчитывать как приближенное значение интеграла, так и оценку погрешности. Кроме того, в силу использования регулярного разбиения п не может изменяться произвольно — оно должно иметь вид п = ms, где т — некоторое натуральное число.

Те же недостатки присущи другому методу с повышенной скоростью сходимости, предложенному Хайнриком (S.Heinrich) [53]. Этот метод представляет собой специальный случай метода выделения главной части, в котором «легкая» аппроксимация строится в конечно-элементной форме (аналогично [8]) на регулярной сетке с числом узлов равным числу используемых случайных точек. Такой подход позволяет получить на классах функций Ст вероятную погрешность, убывающую со скоростью оптимальной на данном классе среди всех недетерминированных методов интегрирования. Однако, в силу указанных недостатков данный метод представляет интерес в основном с точки зрения теории сложности [67].

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

Основы последовательного подхода в рамках общей теории статистического оценивания были заложены Вальдом (A.Wald) [7]. Последовательная вычислительная схема, удобная для исследования методов статистического моделирования, была предложена в 1962 году Холтоном [47] и названа последовательным методом Монте-Карло. В ней рассматриваются несмещенные оценки интеграла двух типов — первичные и вторичные. При этом вторичные оценки представляют собой выпуклые комбинации ранее полученных первичных оценок, коэффициенты которых подбираются так, чтобы обеспечить наилучшую точность. Первичные же оценки могут строиться, вообще говоря, с учетом всей информации о подынтегральной функции, полученной ранее в процессе вычислений. Простейший метод Монте-Карло, методы существенной выборки и выделения главной части представляют собой простейшие варианты последовательной схемы, в которых первичные оценки зависят от текущей случайной точки. Из априорных соображений ясно, однако, что использование более сложных первичных оценок может привести к улучшению сходимости рассматриваемой процедуры. Холтон предложил несколько подобных оценок общего вида, однако теоретическим их исследованием не занимался. Сущность его предложений состоит в том, что по значениям подынтегральной функции в ранее разыгранных точках строится «легкая» функция, моделирующая ее поведение, и по ней уже вычисляется очередная первичная оценка, а также корректируется дальнейший ход алгоритма (например, плотность распределения очередной точки в методе существенной выборки). Важным достоинством подобного подхода является то. что для его реализации не требуется наличия априорной информации о поведении подынтегральной функции на области интегрирования. При заданном способе построения «легких» аппроксимаций алгоритм сам учитывает поведение подынтегральной функции на основе накапливаемой информации о ней, т. е. адаптируется, подстраивается под особенности конкретной задачи.

В течение длительного времени указанные идеи Холтона оставались невостребованными — по всей видимости, из-за предполагаемых сложностей теории. Первые конкретные адаптивные процедуры интегрирования, близкие к предложенным Холтоном, появились в литературе лишь спустя более пятнадцати лет. Среди них в первую очередь следует отметить статью Н. В. Епиховой [15], в которой впервые конкретизирован способ адаптации плотности распределения случайных точек в методе существенной выборки, хотя и отсутствуют оценки погрешности предложенного алгоритма. Первое исследование сходимости адаптивного метода интегрирования было проведено в работе О. Ю. Кульчицкого и С. В. Скроботова [28]. В ней, применительно к исследованию интегральных характеристик сложных систем управления, рассматривается и изучается одномерный адаптивный метод интегрирования, основанный на методе существенной выборки, в котором плотность распределения адаптируется к виду подынтегральной функции, представляя собой кусочно-постоянную ее аппроксимацию по ранее разыгранным случайным точкам. В работе показано, что на классе функций с кусочно-непрерывной ограниченной первой производной, для предложенного метода можно получить оценку дисперсии, убывающую как 0(п-2), т. е. метод имеет ускоренную сходимость. В то же время, он сохраняет простоту классических методов Монте-Карло, имеет последовательный характер и обладает удобным механизмом оценивания погрешности по эмпирической дисперсии в процессе вычислений. Дальнейшее исследование и развитие данного метода, проводимые в работах О. Ю. Кульчицкого, В. М. Иванова, Д. Г. Арсеньева и М. Л. Кореневского [1], [22], [23], [-55] показали, что он может быть обобщен на случай использования более точных аппроксимаций, что приводит к дальнейшему повышению скорости сходимости.

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

Структура и содержание диссертации. Диссертация состоит из введения, четырех глав и заключения. Она изложена на 161 странице, включает 3 таблицы и 12 рисунков.

Список литературы

содержит 67 наименований.

Заключение

.

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

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

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

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

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

5. Предложены адаптивные варианты метода расслоенной выборки и исследована их сходимость.

6. Предложен и исследован способ адаптивного оценивания коэффициентов Фурье разложения функции по заданной ортонормированной системе и построения глобальных аппроксимаций в виде отрезков этого разложения. На его основе построены конкретные адаптивные методы интегрирования и аппроксимации для некоторых классов функций с быстро сходящимися рядами Фурье.

7. Эффективность предложенных адаптивных алгоритмов при вычислении многомерных интегралов подтверждена проведенными численными экспериментами.