автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Разработка адаптивно-статистических методов вычисления определенных интегралов

На правах рукописи УДК 519.245 /

РГ Б ОД

- з май иод

Кореневский Максим Львович

РАЗРАБОТКА АДАПТИВНО-СТАТИСТИЧЕСКИХ МЕТОДОВ ВЫЧИСЛЕНИЯ ОПРЕДЕЛЕННЫХ ИНТЕГРАЛОВ.

Специальность 05.13.18. Теоретические основы математического моделирования, численные методы и комплексы программ

Автореферат диссертации на соискание ученой степени кандидата физико-математических наук

САНКТ-ПЕТЕРБУРГ 2000

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

университете.

Научный руководитель: доктор фнз.-мат. наук,

профессор В.М.Иванов, Санкт-Петербургский государственный технический университет;

Официальные оппоненты: доктор физ.-мат. наук,

профессор В.Б.Мелас, '' Санкт-Петербургский государственный университет;

кандидат техн. наук, доцент А.А.Суханов, Санкт-Петербургский государственный технический университет;

Ведущая организации — ГНЦ РФ - ЦНИИ " Электроприбор",

Санкт-Петербург

Защита состоится " И " 2000 г. в /^ пасов на заседай

диссертационного совета К 063.38.30 СПбГТУ.

По адресу: 195220, Санкт-Петербург, Гражданский пр., д.28, ауд. 4о5>.

С диссертацией молено ознакомиться в фундаментальной библиотеке унив' ситета.

Автореферат разослан " 3 " ¡МхС^!^ 2000 г.

Ученый секретарь диссертационного совета, докт. техн. наук

А . О Я

(^Жг-

Д.Г.Арсен.

Актуальность темы. Одна из самых распространенных и важных за-.ч численного анализа -— приближенное вычисление определенных интегралов, .дачи такого рода возникают в многочисленных областях науки, таких, напри-!р, как теория управления, физика высоких энергий, финансовая математика, оретическая астрономия, и т.д. Вопросам разработки, исследования, тести-вания и применения методов численного интегрирования посвящено огромное личество литературы, включая обширные монографии, справочники и мно-?ство журнальных статей. Тем не менее, достаточно полно разработанной 1Жно считать лишь теорию методов одномерного интегрирования. Задачи приближенного вычисления многомерных интегралов являются зна-тельно более сложными по целому ряду причин, основные из которых — рез-е возрастание трудоемкости традиционных численных методов с ростом раз-рности задачи (так называемое "проклятие размерности'"), а также необхо-мость учета геометрических особенностей возможных областей интегрирова-я. Если вторая из указанных трудностей обычно обходится путем разбиения ласти интегрирования на подобласти, близкие к стандартным (шар. куб, пшене), то первая может оказаться серьезным препятствием при необходимости лисления интеграла достаточно большой размерности с высокой точностью щачи такого рода очень часто возникают как промежуточные этапы при ре-изации сложных вычислительных алгоритмов). Наиболее часто в подобных учаях используются статистические методы интегрирования, основанные на тоде Монте-Карло, скорость сходимости которых не зависит от размерности, [нако, эта скорость является весьма низкой и значительное внимание уделя-поискам путей ускорения сходимости статистических методов. В связи с ;ш совершенствование и разработка статистических методов многомерного тегрирования, обладающих способностью к адаптации, представляется весь-актуальным направлением исследований.

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

тод выделения главной части и обобщенные адаптивные методы.

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

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

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

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

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

Публикации и аппробация работы. Содержание работы отражено в печатных трудах. Результаты диссертации докладывались и обсуждались к семинаре кафедры информатики Санкт-Петербургского государственного те: нического университета (СПбГТУ) , на студенческой конференции в рамках 28-Недели Науки СПбГТУ (1999), на Всероссийской научно-технической конфере!

ни "Фундаментальные исследования в технических университетах" (Санкт-етербург, 1998), на международных научных конференциях "3-rd St.Petersburg 'orkshop on Simulation" (Санкт-Петербург, 1998) и "Second IMACS Seminar on [onte Carlo Methods" (Варна, Болгария, 1999).

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

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

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

Кубатурные формулы представляют собой наиболее широко исследованную >уппу методов многомерного интегрирования. Большой вклад в современную орию кубатурных формул внесли В.И.Крылов, И.П.Мысовских, В.И.Половин-:н, С.Л.Соболев, P.J.Davis, J.Radon. A.Sard, A.H.Stroud, V.Tchakaloff и дру-e ученые. Узлы и веса кубатурных формул выбирают таким образом, что-[ обеспечить точное интегрирование всех функций определенного класса или илучшую на классе погрешность интегрирования. Погрешность при этом -адиционно оценивают на классах Ст функций, имеющих в области интегрп-вания всевозможные ограниченные производные до порядка т включитель-, и близких к ним. Как показал Н.С.Бахвалов, на этих классах погрешность )бой кубатурной формулы с п узлами (и, вообще, любого детерминированно-метода интегрирования, использующего значения подынтегральной функции более, чем в п точках) не может убывать быстрее, чем 0(n~m/'s), где s — змерность вычисляемого интеграла. Таким образом, с ростом размерности теграла скорость сходимости быстро снижается и для достижения хорошей чности необходимо выбирать очень большие значения п, что в свою очередь зывает резкий рост временных затрат на реализацию алгоритма интегрова-я.

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

которым удалось для важных классов функций построить методы интегр! рования, скорости сходимости которых практически не зависят от размернс сти вычисляемого интеграла. Эти методы являются представителями груг пы теоретико-числовых методов (методов квази-Монте-Карло) вычисления m тегралов, в основе которых лежит использование сеток и последоватсльнс стей узлов, равномерно распределенных (в смысле Вейля) в области инт« грирования. Квадратурные формулы с равными весами и узлами, образ} юшими равномерно распределенную последовательность, имеют для многи классов функций (например, функции с ограниченной полной вариацией) скс рость убывания погрешности, близкую к 0(п"1). Основной в теории методо квази-Монте-Карло является проблема построения сеток и последовательнс стей точек, обладающих хорошими свойствами равномерности. Среди мат« матиков, которые внесли серьезный вклад в ее исследование, следует, помим Н.М.Коробова и И.М.Соболя назвать таких ученых как H.Fanre, J.H.Haltor E.Hlawka, J.F.Koksma, H.Niederreiter, K.F.Ptoth, Y.Wang, S.K.Zaremba.

Другой путь построения методов интегрирования, скорость сходимость ко торых не ухудшается с ростом размерности, связан с отказом от получения га рантированной оценки погрешности и переходе к вероятностным оценкам, t.t таким, которые имеют место лишь с некоторым наперед заданным уровне; достоверности. Численные методы, основанные на статистическом моделирс ванпп и использующие такую трактовку сходимости образуют большую груп пу статистических методов (методов Монте-Карло). Большой вклад в развити статистических методов многомерного интегрирования внесли Б.Л.Грановский

C.М.Ермаков, В.М.Иванов. О.Ю.Кульчицкий, Г.А.Михайлов, В.В.Некруткш: И.М.Соболь. А.С.Фролов, Н.Н.Ченцов, S.Haber, J.H.Haltoii, J.M.Hammerslej

D.C.IIandscomb, H.Kahn, J.Spanier и другие ученые.

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

эдынтегральной функции, выборка по группам, случайные квадратурные фор-улы и т.д. Эффективность указанных приемов существенным образом зависит г количества используемой априорной информации, а их скорость сходимости. 1К правило, является той же, что и у простейшего метода Монте-Карло, т.е. (я-1/2) (известны отдельные варианты этих методов с повышенной скоростью содимости, однако они не являются последовательными, что весьма неудобно). В работах О.Ю.Кульчицкого, В.М.Иванова и Д.Г.Арсеньева было установле-что если априорная информация о подынтегральной функции отсутствует, тагожна ее "реконструкция" в в процессе моделирования и учет этой рекон-грукшш в алгоритме. Более того, на примере простейшей задачи одномерного зтегрирования было показано, что такая адаптация алгоритма к свойствам щынтегральной функции позволяет получить повышенную скорость сходимо-:и алгоритмов интегрирования при сохранении основных положительных черт эадиционных методов Монте-Карло.

Настоящее исследование посвяшено дальнейшему развитию и обобщению ■их идей. Его целью является разаботка методики построения адаптивных [горитмов интегрирования, обладающих основными достоинствами традици-[ных статистических методов, а также построение конкретных адаптивных юцедур многомерного интегрирования, имеющих повышенную скорость схо-1М0СТИ по сравнению с известными аналогами.

В первой главе разрабатываются общие принципы построения адаптив-»IX алгоритмов интегрирования на основе "последовательного метода Монте-грло". Эта вычислительная схема, изначально предложенная Дж.Холтоном .Н.НаНоп), предоставляет весьма удобный инструмент для описания различ-■ix адаптивных процедур, что и было использовано в данной главе. При этом зультаты Холтона существенно обобщены и расширены. Итак, пусть требуется оценить неизвестную величину ■/ по последователь-сти некоторым образом смоделированных реализаций случайных величин ,.С2, • ■ и имеется последовательность несмещенных оценок Бк = Зк{%1,..., личины ,/, называемых первичными оценками. Построим вторичные оценки (также несмещенные) по правилу

п

■Гп = + (1 - или Л = 53«(1)

с 0 ^ с*к,Рк 1 — некоторые коэффициенты алгоритма, причем Д — 1, 2=1 = 1 (связь между а^"® и Зь легко устанавливается). В эту схему эследовательный метод Монте-Карло) укладываются многочисленные вари-ты метода Монте-Карло для вычисления интеграла. При этом, как правило.

(п) _1 _

= а\. — п , точки Хк суть независимые реализации одной и той лее случаи-й величины и вк = Бк{хк), т.е. первичные оценки независимы. Этот случай

исследовался самим Холтоном.

В настоящей работе рассмотрен более общий случай, когда первичные оцен ки зависят от всех реализаций Х1,...,Хк и являются некоррелированнымъ Тогда дисперсии вторичных оценок имеют вид

О Ш = £ 4")2 А, Вк = О {5,} = М {\эк - А2} . (2

к=1

где математическое ожидание берется по совокупности величин х\,... ,х/,. Пр] этих условиях удалось доказать следующие утверждения о сходимости после довательного метода Монте-Карло.

Теорема 1. Пусть первичные оценки некоррелированны и их дисперси Бк имеют порядок 0{к~11п° к) для некот.орых у, 6 ^ 0. Пусть коэффициенты последовательного метода Монте-Карло выбраны в виде

Ри = ——7? соответственно а\ ' = +1)——-—, , (3

п + 7 Г{к)Г(п -г + 1)

при произвольном 7' > (7 — 1)/2. Тогда дисперсии вторичных оценок О {1п имеют порядок 0(п-1-71п<5 п), и этот порядок невозможно улучшить другш выбором коэффициентов. Если же дисперсии /Л- имеют, порядок о(к~1\пд к) то дисперсии вторичных оценок И {Л} — о(п~"1_т п).

Теорема определяет скорость сходимости метода в среднеквадратичном смысле Из нес в силу неравенства Чебышева вытекает следующая оценка скорости схо димости по вероятности:

Следствие. Пусть выполнены условия теоремы 1. Тогда для сколь угоды малого 6 > 0 найдется постоянная С > 0, С = С{6). такая что

Р

||Л - 7| ^ Сп"<7+1)/21пг гг} > 1 - 6. (4

Следующие два утверждения относятся к сходимости почти наверное: Теорема 2. Пусть выполнены следующие условия:

1) величины кп = тахог[п' стремятся к нулю при п —> оо;

2) дисперсии первичных оценок Дь есть величины порядка 0(к~~1);

со

3) найдется такое £ > 0, + 1 < 7, для которого сходится ряд Вкк~1 ■

к={

Тогда вторичные оценки сходятся к значению 3 с вероятностью 1.

Следствие. Пусть выполнены условия теоремы 1. Тогда вторичные оцьнк1 сходятся к значению 7 с вероятностью 1 при любом 7 > 1. Кроме этого, для случая некоррелированных первичных оценок в работе по ■пучина чиппричогиш пгггшуя гтттгпрргтттт пппугуятотттяя простой рекуррентны!

ресчет при увеличении п и позволяющая контролировать погрешность в ходе :численнй. Она имеет вид

=-(5)

1 + £ а[п) к=1 к= 1

шляется асимптотически несмещенной при п —> ос оценкой дисперсии D {./„}■ Пусть теперь искомая величина J представляет собой интеграл

J = I" f(x)dx (6)

h

ограниченной замкнутой области О евклидова пространства М5 от некоторой нкции функции f(x) € LПредположим, что параллельно с моделирова-ем случайных точек х\,х-2,... € П производится построение последовательно-i аппроксимаций f\(x), /2(2;),. • ■ подынтегральной функции /(.г), таких что г) - Л-(ж; жь ..., Xk-i) и интегралы

h = J fk{x)dx (7)

п

известны. В этом случае первичные оценки

Sk = Ik + /u(n)[f(xk)- Ыхк)}. (8)

хк — независимые реализации равномерно распределенной в Q случайной шчины, определяют адаптивный метод выделения главной части, а первич-е оценки

с / \ / I \ Л('г') <rv

Sk = —7—г, Рк(х) = рЦхг......rjt_i) = —1.9J

1 Хк распределены в Г2 с условными плотностями Рк{х), определяют адаптив-й мет.од существенной выборки (здесь необходимо дополнительно предпола-■ь неотрицательность f(x) и всех аппроксимаций Л-(.с)). Нетрудно показать. ) такие первичные оценки некоррелированы, т.е. к ним полностью примени-результаты теорем 1 и 2. В случае, когда fk(x) = g(x) для всех к, указан-е адаптивные методы переходят в традиционные методы выделения главной :ти (с главной частью д{х)) и существенной выборки (с плотностью, пропор-шальной д(х)). Если же с ростом к имеет место сходимость fk{x) к /(:г), можно ожидать. что первичные оценки Sk с ростом к будут стремиться к J в соответствии с теоремой 1, обеспечит ускоренную сходимость вторичных нок J, 1 к интегралу J.

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

Я* = МХ1„.....11/(1) - 1к{х)?с1х - - /*|21 . (ю;

ii

= .....,, I/-/«Л. (и:

соответственно, т.е. они тем меньше, чем ближе /¡¿(х) к /(.г) в норме Таким образом, скорость сходимости предложенных адаптивных методов определяется скоростью сходимости последовательности аппроксимаций /*(х) к по дынтегральной функции /(я). В второй и третьей главах рассматривается дв;: различных подхода к построению таких аппроксимаций. В первой главе помимо двух основных адаптивных методов, описанных выше, предложен ряд и> модификаций, которые в некоторых случаях могут оказаться более удобными.

Во второй главе изучен процесс построения аппроксимаций, которые являются кусочными на последовательности разбиении области интегрирования, I: доказаны утверждения о сходимости адаптивных методов интегрирования с такими аппроксимациями.

Пусть по точкам Х\,Х2, - ■ ■ построена последовательность разбиений П1,02,.. области интегрирования О. Каждое разбиение С1к состоит из Л^. подобласте! П* п служит основой для построения кусочной аппроксимации Л(з:). Будел говорить, чго последовательность функций приближает, функцию /{х]

на последовательности разбиений О4 с порядком аппроксимации, I > 0, еслг суу^ествует, постоянная С > 0, такая что

\Нх)~Мх)\^с[ц(п))]1, при хел), (12;

для любых к > 0 и ] — 1,2Аппроксимации, удовлетворяющие этом} условию, могут при заданной последовательности разбиений быть построень для многих классов подынтегральных функций. Одним из наиболее общих таких классов является, по-видимому, класс #(т, а,А), состоящий из функций имеющих всевозможные частные производные до порядка (т — 1), ограничен ные постоянной а, и т-с производные, удовлетворяющие условию Гельдера < постоянной а и показателем А. Для функций этого класса аппроксимация ш данной подобласти может выбираться, например, в виде отрезка ряда Тейлора до порядка т вблизи одной из точек подобласти, что обеспечивает порядог

(Щцриксимацнц I — {т I \)/з.____

Наложение на аппроксимации условия (12) позволяет свести задачу оцени-нпя дисперсий Ик первичных оценок предложенных адаптивных методов к сниванию величины

торую можно назвать моментом порядка 7 разбиения П* относительно ра-хранных точек я^а'г,... .¿ч— ]_ (в адаптивном методе существенной выборки •ебуется дополнительно предположить, что все аппроксимации равно-

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

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

При рассмотрении многомерного случая внимание было ограничено интегри-ванием по гиперпараллелепипеду или, без ограничения общности, по единич-му гиперкубу пространствам5. Однако, даже в этом частном случае попытки тосредствснного переноса описанного выше способа построения разбиения на огомерный случай сталкиваются с серьезными трудностями. В связи с этим л разработан альтернативный способ построения последовательности раз-зний — .метод последовательной бисекции, одинаково применимый как для юмерных, так и для многомерных задач и обеспечивающий оптимальные эядки убывания моментов разбиений. Его сущность молено сформулировать ;дующим рекуррентным правилом: каждое последующее разбиение области тегрирования получается из предыдущего бисе.кцией (делением на 2 равные сти) той из подобластей кусочной аппроксимации, в которую попадает ередная разыгранная случайная точка (при этом подразумевается, что все мбласти — гиперпараллелепипеды).

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

(13)

подобластей Nk разбиения Qk при использовании метода бисекции в точност! равно к, а для моментов разбиений в обоих предложенных адаптивных мето дах получены оценки вида = 0{к1_7), откуда вытекают следующие оцены дисперсий Dk\

Dk = 0(к-71).. (14,

где I - порядок аппроксимации. В результате имеет место следующая теорема Теорема 3. Если выполнено условие (12), то при использовании метода по следователъной бисекции погрешность адаптивных методов выделения глав ной части и существенной выборки (при заданном уровне достоверности имеет порядок 0{п~1~1^2).

Следствие: Если f(x) £ Н(т,а,\). то можно построить адаптивные алгоритмы со скоростью сходимости О (n~(m+A)/s~1/2).

Как показал Н.С.Бахвалов, данная скорость сходимости является оптималь ной по порядку для всех недетерминированных методов интегрирования на классе Н(т, а, А). Она превосходит как скорость сходимости традиционных ыетодог Монте-Карло — О (п-1/2), так и наилучшую скорость сходимости детерминированных методов —

Помимо описанных адаптивных методов, во второй главе предложены также последовательные варианты метода выборки по группам (расслоенной выборки), основанные на методе последовательной бисекции и имеющие оптимальные порядки скорости сходимости на классах Н(т, а, 1) при т = 0,1.

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

Итак, пусть 9i(.t), <ръ{х),... — полная ортонормнрованная в ¿¡¡(П) система функций (вообще говоря, комплексная). Будем строить аппроксимации в виде

Гк

h(x) = Y^CkWix), (15)

¡=i

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

>эффидиенты линейной комбинации совпадают с коэффициентами Фурье при-шжаемой функции. Поэтому наиболее разумным было бы положить

сы = С{ = У /(х)<р{(х)с1х.

16)

[е черта сверху означает комплексное сопряжение, однако величины с, суть нечестные интегралы, вычисление которых равносильно по сложности исходной даче. В связи с этим приходится использовать их приближенные значения. Основная идея, предложенная в главе, состоит в том, чтобы вычислять оцен-[ интегралов С{ параллельно с вычислением исходного интеграла тем лее са->ш адаптивным методом выделения главной части. В самом деле, к моменту строения аппроксимации Д(х-) уже разыграны случайные точки х),..., построены аппроксимации /дх),..., Д._1(ж), и, значит, функции /,(х)<л(£) [я ] < к могут быть использованы в качестве аппроксимаций при построения рвичных оценок адаптивного метода вычисления интеграла с,-. 'Таким обра-ч. предлагается вычислять см по формулам адаптивного метода выделения авной части:

г _ _

с*,- - £ = / + ц{ЩПх,) - (17)

а

торые легко преобразовать в рекуррентную форму:

си =

Е а^цтПхЛ - г > г,_,

7=1

г и /4-1 — некоторым образом выбранные коэффициенты последова-

пьного метода Монте-Карло. Таким образом, набор коэффициентов разложе-я может пересчитываться рекуррентно, с каждой новой разыгранной точкой, ш хранится последовательность невязок [/(^) — /Д^,)]. Погрешность предложенного способа аппроксимации складывается из двух ;тавляющих — остатка ряда Фурье и погрешности вычисления коэффициен-з о,-. Если величины г к с ростом к возрастают медленно, то преобладающим дет, очевидно, первый фактор, а если быстро — второй. Поэтому надо найти зумные условия, при которых темпы убывания остатка ряда Фурье сравнимы скоростью уточнения используемых коэффициентов разложения. Эти усло-I даются следующей теоремой:

орема 4. Пусть выплнены следующие предположения: Существуют, постоянные Сх, 62,71,72 > 0, 6 ^ 0 такие, что для любого

г > 0 имеют место оценки:

Т оо

V < С\гъ для любого г е П: М'2 ^ С2г~ъ 1пд' г, (1Ё

1=1 1=Г+|

2. Коэффициенты последовательного метода Монте-Карло для вычисле.ни величин с*,- выбираются в виде (3) при 2-/ + 1 > 7 = 72/713. Величины гь выбираются в соответствии с формулой

Гк =

(Х(к - 1))^ , х,

д = С'1//(П)Л

где [х] означает целую часть от х, а А > 0 таково, что

(У + 1)' ^

(21' -7 + 1)

Тогда существует постоянная А > 0, такая чт.о для любого к > 0:

Т(к)Ь6 к

А- < -4

Г(А- + 7)

Г2С

(21

(22

т.е. Вк =0 (к 71п к). Если вместо второго условия (19) выполняется

ос

]Г |с;|2 = о(г-'Чп^г), (23

¡=>■+1

то ик = о (к-1 1п* к).

Замечание 1: На самом деле в теореме доказана более сильная оценка — о иен ка величины ||/ - }к\\ь2-, из которой в силу (10) уже вытекает оценка дисперси:

А-

Замечание 2: В случае, когда разложение /(х) по функциям ^¡(х) содержи лишь конечное число слагаемых, скорость сходимости зависит лишь от выбор, коэффициентов последовательного метода Монте-Карло и может быть скол, угодно высокой, т.к. условия (19) будут выполнены при любых положительны: 71.72-

На основании результатов теоремы 4, в главе проведено исследование пред ложенного алгоритма построения глобальных аппроксимаций (и основанноп на нем метода интегрирования) для двух важных классов функций. Пер вый из них - класс функций 5р с быстро сходящимися рядами Фурье-Хаара (Детерминированные методы интегрирования на нем подробно исследовалиа И.М.Соболем.) Для этого класса предложен способ нумерации ортонормиро ванных функций, при котором погрешность интегрирования с заданной досто верностью имеет порядок

IЛ

о I п 11о£

(т » ' '

о2

р < 2.

(24

шная скорость сходимости превосходит скорость сходимости детерминиро-нных методов интегрирования, предложенных И.М.Соболем для функций tacca Sp. Налучшая для данной конкретной функции / 6 Sp оценка погрешно-и, полученная И.М.Соболем, имеет порядок о

Второй рассмотренный класс функций — класс Е°, состоящий из функций Зыстро сходящимися тригонометрическими рядами Фурье. (Класс Е!? и ме-ды интегрирования на нем подробно изучались Н.М.Коробовым.) Для этого асса предложен способ нумерации ортонормированных функций, при котором грешность интегрирования с заданной достоверностью имеет порядок

е а > 1 определяет темпы убывания модулей коэффициентов ряда Фурье, яная скорость сходимости совпадает (с точностью до некоторой степени 1п п) скоростью сходимости детерминированного метода оптимальных коэффици-хов Н.М.Коробова, которая является наилучшей по порядку для детермини-ванных методов интегрирования на классе Е". Уступая методу оптимальных эффициентов в простоте реализации, предложенный метод обладает, однако, умя важными преимуществами, состоящими в последовательном характере толненпя и возможности контроля точности по эмпирической дисперсии в :е вычислений.

Кроме этого, на классе Е° получена оценка погрешности аппроксимации в же пространства непрерывных функций:

■>ядок убывания которой (с точностью до логарифмического множителя) со-щает с порядком, оптимальным для детерминированных методов аппроксн-ции. Следует, однако, отметить, что наилучший среди известных детермини-шшых методов аппроксимации на классе Е* (предложенный С.А.Смоляком и Г.Рябеньким), обеспечивает порядок скорости сходимости вдвое хуже, чем в сношении (26). Таким образом, предложенный метод построения аппрокси-дий является существенным улучшением по сравнению с известными ранее. В четвертой, заключительной главе описаны результаты численных экспе-.гентов по применению предложенных адаптивных методов интегрирования >яду тестовых задач. Проведен сравнительный анализ результатов рабо-адаптивных методов, кубатурных формул, квазимонтекарловскнх методов гегрирования и различных вариантов метода Монте-Карло. На его осно-гформулированы практические рекомендации по применению предложенных хшов для вычисления многомерных интегралов.

Эсновная тестовая задача, рассмотренная в главе, состоит в вычислении ин-

(25)

(26)

теграла по единичному s-мсрному гиперкубу от функции

f(x) = il2-^{7rXk) + a\ *ек\ *=Л l + ak

где Uk > 0 параметры, определяющие степень зависимости подынтегральнс функции от переменной хк (при больших ак зависимость слабая, при ак бли' ких к нулю — существенная). Подынтегральная функция, аналогичная f(x рассматривалась в работах И.М.Соболя при исследовании эффективности kbó зимонтекарловских методов интегрирования. Функция f(x) бесконечно дифф< ренцируема и принадлежит классам Sp с любым р > 1 и классу ЕПоэтому ней могут быть применены практически все адаптивные методы, предложенны в настоящей работе.

Интеграл от f(x) равен единице при любых значениях a¡¡. В то же врем) поведение f(x) на области интегрирования существенно зависит от выбора эти величин. При проведении численных экспериментов было рассмотрено 4 ра: личных набора параметров ак:

1. ai = a-2 = ... = a-s = 0.01;

2. ai = ü2 — . ■. = ah = 1;

3. aje = fc, 1 ^ k ^ s;

4. ak - k2, 1 k < s.

В первых двух случаях наибольшее значение f(x). принимаемое в вершина гиперкуба, экспоненциально зависит от размерности s, т.е. с ростом s вблиз вершин появляются области с большими градиентами функции /(х). В Tpt тьем случае наибольшее значение линейно по отношению к размерности, а четвертом — ограничено. Поэтому можно ожидать постепенного улучшени сходимости методов интегрирования от первого случая к четвертому.

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

В ходе экспериментов были проведены расчеты по всем четырем случаям вы бора параметров ак для различных значений s. Результаты расчетов позволяю-сделать следующие выводы:

• В задачах малой размерности (s = 1-4) адаптивные методы с кусочным) аппроксимациями mol у i у сненша^гошасрировать как с детеоминированны ми методами, так и с различными вариантами метода Монте-Карло:--

• В задачах с размерностью s = 5-15 адаптивные методы с кусочными аппроксимациями являются эффективными при достаточно высоких требованиях к точности вычисления интеграла, особенно в случае отсутствия априорной информации о подынтегральной функции. Также эффективными представляется использование адаптивных методов в задачах аппроксимации.

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

• В задачах большой размерности (s > 15) адаптивные методы с кусочными аппроксимациями будут предпочтительней простейшего метода Монте-Карло и его модификаций в случае высокой гладкости подынтегральной функции, высокого порядка аппроксимации и высоких требований к точности.

Вторая тестовая задача взята из работы С.М.Ермакова и Г.А.Илюшиной и стоит в интегрировании по пятимерному единичному гиперкубу функции

/(.г) = ехр {.сгАг} ,

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

Результаты диссертации отражены в следующих публикациях.

1. Иванов В.М., Кореневский M.JL, Кульчицкий О.Ю. // Адаптивно-зтистические методы численного интегрирования с повышенной скоростью ушмости. Фундаментальные исследования в технических университетах анкт-Петербург, 25-26 июня 1998 г.): тез. докл., СПб, 1998, с. 87-88.

2. Ivanov V.M., Korenevsky M.L., Kul'chitsky O.Yn. // Adaptive Monte Carlo orithms with advanced accuracy. 3-rd St.Petersburg Workshop on Simulation анкт-Петербург, 28 июня - 3 июля 1998 г.): тез. докл., СПб, СПбГУ, 1998, 36-39.

3. Иванов В.М., Кореневский М.Л., Кульчицкий О.Ю. Адаптивные схемы года Монте-Карло повышенного порядка точности. // Доклады Академии ук России, 1999, т. 367, №5, с. 590- 593.

4. Ivanov V.M., Korenevsky M.L. Superccmvergent adaptive stochastic metho for numerical integration. // Second IMACS Seminar on Monte Carlo Method (Варна, Болгария, 7-11 июня 1999 г.): тез. докл., Sofia, CLPP-BAS, 199Í с. 16-17.

5. Арсеньев Д.Г., Кореневский M.JL, Кульчицкий О.Ю. О повышении то*; ности статистических методов расчета интегральных характеристик сложны систем. // Известия высших учебных заведений России. Радиоэлектронике Принята к печати.

6. Кульчицкий О.Ю., Арсеньев Д.Г., Кореневский M.JI. Методы адаптивне го управления сеткой в задаче расчета интегральных характеристик сложны: систем. // Автоматика и телемеханика. Принята к печати.

7. Кореневский М.Л. О последовательном методе расслоенной выборки Научно-техническая конференция студентов (в рамках 28-й Недели наую СПбГТУ): тез. докл., СПб, СПбГТУ, 1999.

8. Ivanov V.M., Korencvski M.L. Sequential Monte Carlo and adaptive multi dimensional integration. // International Conference on Monte Carlo Simulatioj (Монте-Карло, Монако, 18-21 нюня 2000 г.). Включена в программу конферен ции.

Введение . . о

Глава 1. Последовательный метод Монте-Карло и адаптивное интегрирование.

1.1 Постановка задачи.

1.2 Метод Монте-Карло вычисления определенных интегралов

1.2.1 Простейший метод Монте-Карло.

1.2.2 Метод существенной выборки.

1.2.3 Метод выделения главной части.

1.2.4 Достоинства и соотношение методов существенной выборки и выделения главной части

1.3 Последовательные схемы метода Монте-Карло

1.3.1 Основные соотношения.

1.3.2 Сходимость в среднеквадратичном смысле

1.3.3 Сходимость почти наверное.

1.3.4 Оценивание погрешности в процессе вычислений

1.4 Адаптивные методы интегрирования.

1.4.1 Простейший адаптивный метод одномерного интегрирования.

1.4.2 Адаптивный метод существенной выборки

1.4.3 Адаптивный метод выделения главной части

1.4.4 Обобщенные адаптивные методы существенной выборки и выделения главной части .

1.4.5 О затратах времени и оперативной памяти

Глава 2. Методы адаптивного интегрирования на основе кусочных аппроксимаций.

2.1 Кусочные аппроксимации по подобластям.

2.1.1 Кусочные приближения и порядок аппроксимации

2.1.2 Построение аппроксимаций для конкретных классов функций

2.1.3 Моменты разбиения и оценки дисперсий • •

2.1.4 Обобщенные адаптивные методы.

2.2 Простейший одномерный метод.

2.2.1 Метод выделения главной части.

2.2.2 Метод существенной выборки.

2.2.3 Выводы и замечания.

2.3 Последовательная бисекция.

2.3.1 Описание способа бисекции.

2.3.2 Метод выделения главной части.

2.3.3 Метод существенной выборки.

2.3.4 Временные затраты при использовании метода бисекции.

2.4 Интегрирование функций с локальным нарушением гладкости.

2.5 Последовательный метод расслоенной выборки

2.5.1 Основные соотношения метода расслоенной выборки.

2.5.2 Метод расслоенной выборки с последовательной бисекцией.

Глава 3. Методы адаптивного интегрирования на основе глобальных аппроксимаций.

3.1 Глобальные аппроксимации.

3.1.1 Аппроксимации по ортонормированиым функциям. Основные соотношения.

3.1.2 Условия сходимости алгоритма.

3.2 Адаптивное интегрирование на классе Sp.

3.2.1 Система функций Хаара и классы функций Sp одной переменной.

3.2.2 Адаптивное интегрирование на классе Sp. Одномерный случай.

3.2.3 Разложение на разноразмерные слагаемые. Многомерные классы Sp.Ill

3.2.4 Адаптивное интегрирование на классе Sp. Многомерный случай.

3.3 Адаптивное интегрирование на классе

3.3.1 Классы функций Е®

3.3.2 Адаптивное интегрирование с тригонометрическими аппроксимациями.

Глава 4. Численные эксперименты.

4.1 Постановка тестовых задач.

4.1.1 Первая задача.

4.1.2 Вторая задача.

4.2 Описание программы.

4.3 Результаты экспериментов

4.3.1 Первая тестовая задача

4.3.2 Вторая тестовая задача

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

Вопросам разработки, исследования, тестирования и применения численных методов приближенного вычисления определенных интегралов посвящено огромное количество литературы, включая обширные монографии [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 variâtes 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. Эффективность предложенных адаптивных алгоритмов при вычислении многомерных интегралов подтверждена проведенными численными экспериментами.

1. Арсеньев Д.Г., Иванов В.М., Кульчицкий О.Ю. Адаптивные методы вычислительной математики и механики. Стохастический вариант. СПб.: Наука, 1996.

2. Бабаев А. О способах ускорения сходимости при вычислению интегралов по методу Монте- Карло. // Вопросы вычислительной и прикладной математики. Ташкент, 1978, №53, с. 8-17.

3. Бахвалов Н.С. О приближенном вычислении кратных интегралов. // Вестник МГУ, серия мех., матем., астрон., физ., 1959, №4, с. 3-18.

4. Бахвалов Н.С. Численные методы. М.: Наука, 1973.

5. Березин И.С., Жидков Н.П. Методы вычислений. В 2-х томах. М.: Наука, 1966.

6. Вальд А. Последовательный анализ. М.: Физматгиз, 1960.

7. Войтишек A.B., Колюхин Д.Р. Использование аппроксимации Стренга-Фикса при вычислении многократных интегралов методом Монте-Карло. // Препринт, Новосибирск, ВЦ СО РАН, 1996.

8. Воронцов Ю.В. О вычислении многомерных интегралов методом Монте-Карло с использованием ортонормированных функций. // Тр. Радиотехн. инст-та АН СССР, 1976, №24, с. 77-84.

9. Грановский Б.Л. О случайных квадратурах гауссовского типа. // Ж. вычисл. мат. и матем. физ., 1968, т. 8, №4, с. 879-884.

10. Грановский Б.Л. Условия допустимости одного класса квадратурных формул со случайными узлами. // Математ. заметки, 1974, т. 16, №2, с. 293-304.

11. Грановский Б.Л., Ермаков С.М. Случайные квадратуры с частично фиксированными узлами. // В сб.: Методы вычислений, вып.6, 1970, с. 79-88.

12. Дейвид Г. Порядковые статистики. М.: Наука, 1979.

13. Демидович Б.П., Марон И.А. Основы вычислительной математики. М.: Наука, 1970.

14. Епихова Н.В. Адаптивный метод Монте-Карло п его применение в задачах нелинейной идентификации. // Автоматика и телемеханика, 1981, N210, с. 98-104.

15. Ермаков С.М. Случайные квадратуры повышенной точности. // Ж. вычисл. мат. и матем. физ., 1964, т. 4, №3, с. 550-554.

16. Ермаков С.М. О понятии допустимости в методе Монте-Карло. //В сб.: Методы вычислений, вып.6, 1970, с. 72-79.

17. Ермаков С.М. Метод Монте-Карло и смежные вопросы, изд. 2-е. М.: Наука, 1975.

18. Ермаков С.М., Золотухин В.Г. Полиномиальные приближения и метод Монте-Карло. // Теор. вер. и ее прим., 1960, т. 5, №4, с. 473 476.

19. Ермаков С.M., Илюшина Г.А. О вычислении одного класса интегралов с помощью случайных интерполяционно- квадратурных формул. //В сб.: Методы Монте-Карло в вычислительной математике и математической физике. Новосибирск, 1974, с. 67-71.

20. Ермаков С.М. Михайлов Г.А. Курс статистического моделирования. М.: Наука. 1976

21. Иванов В.М., Кореневский М.Л., Кульчицкий О.Ю. Адаптивные схемы метода Монте-Карло повышенного порядка точности. // Докл. АН России, 1999, т. 367, №5, с. 590-593.

22. Коробов Н.М. Теоретико-числовые методы в приближенном анализе. М.: Физматгиз, 1963.

23. Крамер Г. Математические методы статистики. М.: ГИИЛ, 1948.

24. Крылов В.И. Приближенное вычисление интегралов. М.: Наука, 1967.

25. Крылов В.Pl., Шульгина Л.Т. Справочная книга по численному интегрированию. М.: Наука, 1966.

26. Кульчицкий О.Ю., С'кроботов C.B. Адаптивный алгоритм метода Монте-Карло для расчета интегральных характеристик сложных систем. // Автоматика и телемеханика, 1986, №6, с. 88-95.

27. Ланцош К. Практические методы прикладного анализа. М.: Физматгиз, 1961.3132 333438