Скалярные алгоритмы метода Монте-Карло для решения метагармонических уравнений

Официальной датой рождения метода Монте-Карло принято считать 1949 год, когда была опубликована статья [46) С. Улама и’Н. Мет-рополиса. Сам термин был предложен ещё во время Второй мировой войны выдающимися учеными XX века математиком Дж. фон Нейманом и физиком Энрико Ферми в Лос-Аламосе (США) в процессе работ по ядерной тематике. Хотя методы Монте-Карло были известны и до 40-х годов, интенсивное развитие статистическое моделирование получило несколько позже в связи с появлением компьютеров, что позволило проводить вычисления больших объемов. С другой стороны, всё большее распространение получает статистическое описание тех или иных сложных физических процессов в связи с чем методы Монте-Карло всё более активно используются во многих научных областях (теория переноса, теория массового обслуживания, теория надежности, статистическая физика и др.).

Основными преимуществами данных методов являются.

• физическая наглядность и простота реализации,.

• малая зависимость трудоемкости задачи от её размерности,.

• возможность решения задач со сложной геометрией,.

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

• вероятностные представления позволяют строить обобщённые решения уравнений,.

• одновременное оценивание вероятностной погрешности оценки искомого функционала,.

• простое распараллеливание методов.

Одна из схем решения краевых задач методом Монте-Карло заключается в сведении исходной дифференциальной задачи к некоторому интегральному уравнению, что даёт возможность использовать развитой аппарат методов Монте-Карло [10, 11] для решения интегральных уравнений второго рода. Среди подходов такого рода выделим следующие два. .

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

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

Другой способ состоит в использовании тесной связи эллиптических операторов с диффузионными процессами (см., например, [5, 40]) и основан на вероятностном представлении решения исходной дифференциальной задачи в виде функционала от решения соответствующей системы. стохастических дифференциальных уравнений (СДУ). В качестве вспомогательной задачи в диссертации рассматривается n-мерная задача Дирихле для эллиптического уравнения г 1 / л д2и f,. ди.

Lu + cu+g = - + bi® — + cu+g = Ot и ^ = <р (1) i, j=1 1 3 г=1 г в ограниченной области D с границей, которая предполагается односвяз-ной и регулярной. Дифференциальному оператору L из (1) можно поставить в соответствие систему СДУ в смысле Ито [38]: t t & = & + J b (Zi)dl + J v (b)dw (l), ' (2) о 0 где G Rra, fj® — нижняя треугольная матрица диффузии, определяемая разложением Холесского [13]:

А = а (г) • ат (г), a w (t) — стандартный винеровский процесс.

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

L + с) р+1и = —g, (L + с)*и|г = у*, А = 0,., р, (3) при L = Д. Здесь используется модификация винеровского процесса, называемая процессом «блуждания по сферам» или «сферическим процессом» .

Впервые эвристическое описание алгоритма «блуждания по сферам» было дано Брауном [4] для приближенного решения задачи Дирихле для уравнения Лапласа. Строгое обоснование алгоритма было предложено Мюллером в [45], где изучался также вопрос об эффективности такого алгоритма. Дальнейшее развитие алгоритмы «блуждания по сферам» получили в работах Михайлова Г. А. и Елепова Б. С. В работе [9] был предложен алгоритм для решения уравнения Гельмгольца, использующий вероятностное представление решения. Позднее в [8] был построен алгоритм для решения этой задачи, исходя из специального интегрального уравнения второго рода. Данный подход позволил разработать алгоритмы статистического моделирования для решения более общих уравнений и систем (см., например, [12, 32, 33, 41]). Алгоритм «блуждания по сферам» основан на моделировании точек последовательного выхода винеровского процесса из максимальных сфер, целиком лежащих в рассматриваемой области. В работах [7, 8, 12, 33] для широкого класса границ Г была получена логарифмическая оценка для среднего числа переходов q (ro, s) в цепи «блуждания по сферам» до момента попадания в^{-окрестность} границы:

Щ (го, е) < Сп 1пе|.

Кроме того, для больших значений п размерности пространства верна [45] линейная зависимость Сп ~ п.

В первой главе диссертации путём итерирования вероятностного представления решения задачи Дирихле для эллиптического уравнения (1) построено вероятностное представление для решения метаэллиптиче-ского уравнения (3). Это даёт возможность, в рамках третьего подхода, применить для нахождения функционалов от решения хорошо развитую теорию оценок метода Монте-Карло, использующих сферический процесс. Разработанные алгоритмы «блуждания по сферам» для решения как детерминированных, так и стохастических задач приводятся в третьей главе диссертации. Следует отметить, что в рамках первого и второго подходов векторные оценки для решения метагармонического уравнения были получены в [33]. Однако в этой работе было рассмотрено только однородное уравнение. Кроме того, полученные с помощью «двойной рандомизации» [41] оценки «блуждания по малым сферам» для ковариационной функции решения задачи А2и + ки — f со случайными параметрами имеют существенно большую дисперсию, сравнительно с новыми оценками из п. 3.4 и п. 3.5, полученными частичным осреднением скалярных оценок.

Для приближенного статистического моделирования решения системы СДУ (2) можно, в частности, использовать простой в реализации и физически наглядный дискретный метод Эйлера с постоянным шагом At [35]: rn+1 = Гп + b (rn)At + cr{rn)r)ny/Ai, (4) где гп — численная оценка решения (2) в узлах равномерной сетки по времени {пAt}, а {г)п} - последовательность независимых между собой случайных векторов с независимыми стандартными гауссовскими компонентами. Однако из-за невысокой скорости сходимости метода (4) приходится использовать методы более высоких порядков с более сложной реализацией (см., например, [43, 35]). В диссертации предложен дискретный метод Эйлера для решения бигармонического уравнения, который естественным образом распространяется для решения метаэллиптического уравнения.

Во второй главе диссертации использован ещё один подход к решению краевых задач, широко применяемый в практических расчётах в силу своей простоты и сравнительной универсальности. Речь идёт о методе «блуждания по решёткекоторое может интерпретироваться как дискретный вариант упомянутых выше подходов. Этот подход прост в реализации, не требует большой памяти ЭВМ, но является сравнительно трудоёмким из-за необходимости моделирования длинных случайных траекторий [32, 34]. Для глобальных оценок решения рассматриваемых задач приходится решать задачу выбора оптимального распределения начальной точки траектории, минимизирующую трудоёмкость метода.

Величиной, определяющей эффективность алгоритмов метода Монте-Карло при их практическом использовании, является трудоёмкость 5© случайной оценки? (см., например, [11]). В диссертационной работе под трудоёмкостью i5© будем понимать среднее количество вычислительной работы, необходимой для достижения заданной погрешности:

5© = DC • где t (() — среднее время, затрачиваемое на моделирование одного выборочного значения случайной величины С (например, для процесса «блуждания по сферам» величина ?(?) определяется средним числом переходов в цепи до момента попадания в? — окрестность границы Г). Поскольку общая погрешность оценки, полученной методом Монте-Карло, обычно состоит из двух частей: детерминированной и вероятностной, целесообразно выбирать параметры алгоритма таким образом, чтобы обе эти погрешности имели один и тот же порядок. Вероятностная погрешность возникает в связи с заменой математического ожидания Е£ ариф.

1 N метическим средним выборочных значений — Q и пропорциональна j=1 у/Щ величине —., где N — количество моделируемых выборочных значений. Детерминированная погрешность появляется, например, при введении е — окрестности границы (в алгоритмах «блуждания по сферам»), при замене диффузионного процесса кусочно-линейным процессом с шагом по времени At (в методе Эйлера) или при замене исходной дифференциальной задачи разностной (в алгоритмах «блуждания по решётке»).

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

Далее следует краткое содержание диссертационной работы по главам и параграфам.

В первой главе исходя из стандартного вероятностного представления решения задачи Дирихле для уравнения эллиптического построено вероятностное представление решения первой краевой задачи для ме-таэллиптического уравнения. Это представление показывает, что рассматриваемое решение может быть получено путём параметрического дифференцирования решения специальной краевой задачи для эллиптического уравнения. Таким образом, открывается возможность построения и обоснования требуемых алгоритмов метода Монте-Карло путём параметрического дифференцирования известных стандартных оценок, связанных с «блужданием по решётке» и «блужданием по сферам». Результаты главы 1 опубликованы в [17, 20, 24].

В параграфе 1.1 вводится вспомогательная задача для эллиптического уравнения в области D С Rn.

L + с 4-)и = — д, и г = ?>- формулируются условия, обеспечивающие существование и единственность решения данной задачивыписывается конкретный вид вероятностного представления для этого решения.

В параграфе 1.2 рассматривается краевая задача для метаэллипти-ческого уравнения.

L + cf+lu = д, (L + с) ки.

T = Wь, к = 0,. ., р.

Посредством итерирования вероятностного представления для решения эллиптического уравнения строится вероятностное представление решения и = Щр = Е т р I.

— 1Г1 [ ^-(^(eOA + ^t^c-^r'wfe) { Р' 1=0 для предварительно расщеплённой задачи настоящего параграфа. На основе данного представления показывается, что решение метаэллипти-ческого уравнения может быть получено путём параметрического дифференцирования стандартного решения специальной задачи для эллиптического уравнения. Доказывается конечность дисперсии случайной величины Ср

В параграфе 1.3 для полученного вероятностного представления реализована схема Эйлера с постоянным шагом по времени с целью нахождения решения бигармонического уравнения Л Ли = д. В конце параграфа 1.3 приведены результаты численных расчётов и сделано замечание об улучшении порядка погрешности в методе Эйлера.

Во второй главе построены и обоснованы новые весовые методы Монте-Карло для оценки решения задачи Дирихле для многомерного разностного бигармонического уравнения на основе моделирования «блужданйя по решетке». Векторные варианты построенных алгоритмов непосредственно распространяются на разностные метагармонические уравнения с сохранением вида условий несмещенности оценок и ограниченности их дисперсий. В связи с этим построен простой алгоритм для оценки первого собственного числа многомерного разностного оператора Лапласа. Кроме того, построены специальные алгоритмы «блуждания по решетке», позволяющие при определенных условиях оценивать решения задачи Дирихле для бигармонического уравнения со слабой нелинейностью и для задач со смешанными краевыми условиями, включающими условие Неймана. Результаты главы 2 опубликованы в [16, 18, 20, 22, 23].

В параграфе 2.1 формулируется разностная задача, аппроксимирующая вспомогательное уравнение Гельгольцавводятся соответствующие обозначенияопределяется случайный процесс «блуждания по решётке» .

В параграфе 2.2 с помощью приёма описанного в параграфе 1.2 строится оценка решения разностного метагармонического уравнения (Дл + cby+W = gh. Доказывается теорема о равномерной ограниченности полученной оценки.

В параграфе 2.3 рассматривается задача вида:

Д + с)(Д + 6) и = Аи + Ьи ср и Ч> и эквивалентная ей система уравнений.

А и + bu = v, и.

Av + cv — —д, v V.

В соответствующую оценку решения разностной задачи для первого уравнения данной системы входят неизвестные значения vно можно использовать их случайные оценки на этой же траектории. Таким образом, строится оценка для решения исходной задачи. На основе векторного представления доказываются конечность дисперсии данной оценки. Из рекуррентного представления оценки выводится равномерная по h ограниченность дисперсии.

В параграфе 2.4 строится глобальная оценка решения задачи из параграфа 2.3. Применительно к таким задачам известно несколько подходов к моделированию цепи [32, 30]. В этом параграфе было выбрано равномерное по сеточной границе области начальное распределение цепи, для которого среднее число шагов есть величина порядка 0{h~l). Для глобальной оценки, которая получается линейным восполнением оценок решения в узлах равномерной сетки, решается задача минимизации трудоёмкости.

В параграфе 2.5 рассматривается задача вида:

Л Аи + сАи + Ьи = — д, и.

Аи (р.

С помощью рандомизации соответствующей разностной задачи получается скалярная оценка решения. Доказывается конечность дисперсии представленной оценки.

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

ААи + сАи + Ъит — -д.

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

В параграфе 2.6 оценка «блуждания по решётке», построенная в п. 2.3, распространяется на следующую задачу Дирихле для метагар-монического уравнения:

Л + CW)(A +. (Д + с^)и = -д, т iJk г г с комплексными коэффициентами При этом обосновывается конечность и равномерная ограниченность по h дисперсии полученной оценки.

На основе соотношения (при с = const) p (ttA)(p-1) uh)(p) р->оо где uh — решение задачи.

A huh + cuh = 0, uh г = 1, изучается метод Монте-Карло для оценивания с* по формуле р[1 — ch2/(2п)] с* «Е с.

AT +р — l)[/i2/(2n)] j Здесь —с* - первое собственное число многомерного оператора Лапласа для области D.

В параграфе 2.7 приводятся результаты численных расчётов. В пп. 2.7.1 с помощью оценки, полученной в п. 2.3 численно решаются би-гармоническое и метагармоническое уравнения. В пп. 2.7.2 проводится сравнение двух алгоритмов, соответствующих «блужданию по решётке» и методу Эйлера.

В третьей главе при построении оценок параметрических производных в случае «блуждания по сферам» преодолены существенные трудности. Здесь дополнительно необходимо изучать смещение «еоценок» и обосновывать равномерную ограниченность их дисперсий. Это сделано в настоящей работе на основе использования «мартингального свойства» стандартного вероятностного представления, которое состоит в следующем: представление сохраняет свойство несмещённости, если траекторию процесса оборвать в некоторый момент времени, добавив к оценке значение стандартного решения с соответствующим весом. Таким образом построены и обоснованы оценки «блуждания по сферам» для решения первой краевой задачи для метагармонического уравнения.

В случае бигармонического уравненияА2и — —д полученные оценки записаны в явном виде. Они представляют собой линейные комбинации значений функции д и граничных функций. Это дало возможность построения просто реализуемых алгоритмов оценки ковариационных функций решений уравнений А2и = д, А2и + си = д со случайными функциональными параметрами. Показано, что неограниченность числа шагов в случайном блуждании не является существенной. Результаты главы 3 опубликованы в [17,19−21,24].

В параграфе 3.1 вводятся обозначения и определения, связанные с процессом «блуждания по сферам» .

В параграфе 3.2 на основе полученного в п. 1.2 вероятностного представления строятся оценки решения метагармонического уравнения, связанные с процессом «блуждания по сферам». Доказывается конечность дисперсий и обосновывается^{-смещённость} построенных оценок.

В параграфе 3.3 выписывается конкретный вид оценок решения бигармонического уравнения для случая двумерной и трёхмерной области D.

В параграфе 3.4 рассматривается задача Дирихле для бигармонического уравнения.

ДДи = —д, А и и г ¥->о, со случайными функциональными параметрами д, (fo, <р. Путём условного (для фиксированных траекторий) осреднения произведения соответствующих оценок решения в точках г и т¹) строится оценка ковариационной функции v (r, И1)) решения задачи, записанной выше. Для реализации полученного метода предлагается при превышении количества точек {г*} процесса «блуждания по окружностям» некоторого достаточно большого уровня М не сохранять информацию о следующих точках с номерами г = М+1,., N — 1, но, при попадании в ег-окрестность границы, использовать полученные вес Qm и координаты г мПри этом Qi = О для г = М+1,., N — 1, и можно заменить Qm на Qm{N-M— 1). Показывается, что порядок (по е) величины М можно эвристическиоценить по формуле.

М х 11пе|2.

В параграфе 3.5 рассматривается задача вида:

А 2и + си = д, и.

Ее = 0, Dc<�с*.

Решение и (г, с) аппроксимируется суммой двух первых членов ряда Лорана, по формуле и (г, с) «w (r, 0) + 0) с, с погрешностью приблизительно равной 0) с2. Выписывается конкретный вид оценок параметрических производных и^(г, 0), и^(г, 0). Путём условного осреднения (для фиксированных траекторий г*, г-1^ и значений случайных параметров д, сро, <рi) произведения полученных оценок решения в точках г* и г^ строится оценка ковариационной функции г-(г1,г^) решения исходной задачи этого параграфа. Оценивается погрешность предложенного метода Монте-Карло.

В параграфе 3.6 приводятся результаты числовых расчётов для оценок, полученных в третьей главе.

В пп. 3.6.1 решается неоднородное бигармоническое уравнение с помощью оценки, выписанной в п. 3.3.

В пп. 3.6.2 тестируются два алгоритма} предложенных в п. 3.4. В алгоритме А, после того как количество центров сфер из соответствующей оценки достигало уровня М, случайное «блуждание по кругам», исходящее из точки г, прекращалось. В алгоритме В в подобном случае процесс не обрывался, но соответствующие веса Qi и случайные точки не запоминались, а последний вес Qn0 домножался на величину (N — М + 1). Численно показывается, что эвристическая оценка для М (М х | In £г|2) является несколько завышенной.

В пп. 3.6.3 сравниваются алгоритм, полученный в п. 3.5, и рандомизированный алгоритм «блуждания по окружностям и кругам» [27]. Результаты расчётов подтверждают, что предложенный в настоящей работе алгоритм имеет значительно меньшую (приблизительно в 5 раз) дисперсию.

На защиту выносятся следующие основные результаты.

1. Построены вероятностное представление для решения метаэллип-тического уравнения (L + с) р = д и соответствующий метод Монте-Карло, связанный с дискретным методом Эйлера.

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

3. Разработан простой алгоритм для оценки первого собственного числа многомерного разностного оператора Лапласа.

4. Построены оценки для параметрических производных от решения задачи Дирихле для уравнения Гельмгольца.

5. Построена скалярная оценка метода Монте-Карло, связанная с-. процессом «блуждания по сферам», для решения многомерного полигармонического уравнения.

6. Предложены алгоритмы, позволяющие при определённых условиях относительно случайныхпараметров дне оценивать ковариационную функцию решения уравнений А2и = д, А2и + си = д.

7. Осуществлен ряд численных расчетов, подтверждающих теоретические результаты.

Результаты, изложенные в диссертации, докладывались и обсуждались на семинаре Отдела статистического моделирования в физике ИВМиМГ СО РАН, на конференциях молодых ученых ИВМиМГ СО РАН (1999 — 2005 гг.) — на международной конференции по вычислительной математике в г. Новосибирске (ICCM-2004) — на объединенных шестой международной конференции по методам Монте-Карло и квази Монте-Карло в научных вычислениях и второй международной конференции по методам Монте-Карло и вероятностным методам для решения дифференциальных уравнений в частных производных в г. Жуан-Ле-Пен, Франция (MC2QMC2004) — на четвертом (2001г.) и пятом (2005г.) международном семинаре по математическому моделирования в г. Санкт-Петербурге.

Пользуясь случаем, автор выражает искреннюю благодарность своему научному руководителю члену-корреспонденту РАН Геннадию Алексеевичу Михайлову за постоянное внимание и плодотворное руководство работой.

Заключение

.

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

1. Построены вероятностное представление для решения метаэллип-тического уравнения (L + с) р = д и соответствующий метод Монте-Карло, связанный с дискретным методом Эйлера.

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

3. Разработан простой алгоритм для оценки первого собственного числа многомерного разностного оператора Лапласа.

4. Построены оценки для параметрических производных от решения задачи Дирихле для уравнения Гельмгольца.

5. Построена скалярная оценка метода Монте-Карло, связанная с процессом «блуждания по сферам», для решения многомерного полигармонического уравнения.

6. Предложены алгоритмы, позволяющие при определённых условиях относительно случайных параметров д и с, оценивать ковариационную функцию решения уравнений А2и = д, А2и + си = д.

7. Осуществлен ряд численных расчетов, подтверждающих теоретические результаты.