Методы Монте-Карло и Квази Монте-Карло для решения систем линейных алгебраических уравнений тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

На правах рукописи

Рукавишникова Анна Игоревна

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

01.01.07 - вычислительная математика

АВТОРЕФЕРАТ диссертации на соискание уданпй степени кандидата физико-матема

и0346Ю56

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

003461056

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

Научный руководитель:

Официальные оппоненты:

Ведущая организация:

доктор физико-математических наук, профессор

ЕРМАКОВ Сергей Михайлович

доктор физико-математических наук, профессор

КОРНЕЕВ Вадим Глебович (Санкт-Петербургский государственный политехнический университет)

кандидат физико-математических наук, инженер

КУЧКОВА Ирина Николаевна (Sun Microsystems)

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

Защита состоится

аил

2009 г. в ^ часов на заседании совета Д 212.232. 4 9 по защите докторских и кандидатских диссертаций при Санкт-Петербургском государственном университете по адресу: 199048, Санкт-Петербург, Васильевский остров, Ц-я линия, д. 29, математико-механический факультет, ауд. 13.

С диссертацией можно ознакомиться в библиотеке им. Горького Санкт-Петербургского государственного университета по адресу: 199034, Санкт-

Петербург, Университетская наб. д.7/9.

»)

Автореферат разослан

. 200.2 г.

Ученый секретарь

диссертационного совета

Архипова А.А.

1. ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

Интерес к этим методам обусловлен, в том числе, тем, что как показано в работах Ермакова С.М. и Данилова Д.Л., при большой размерности СЛАУ метод МК обладает меньшей трудоемкостью (в статистическом смысле), чем детерминированные итерационные методы. Также, большое преимущество перед другими методами дает методу МК его естественное свойство параллелизма. В некоторых случаях методы МК и КМК допускают неограниченный параллелизм и, следовательно, полную загрузку оборудования, что особенно важно для решения задач большой размерности.

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

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

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

Применение методов КМК для решения СЛАУ приводит к возникновению новой проблемы, а именно, к необходимости использовать квазислучайные точки большой размерности в для оценки высокой степени матрицы оператора. Это приводит к ухудшению сходимости из-за наличия множителя 1п3 N в числителе оценки скорости убывания остатка КМК. Про-

блема отмечается многими авторами (Ермаков С.М., Лекот К. и Лекюйе П.). Она решается различными предложенными ими вариантами модификации КМК, но лишь частично.

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

Далее в диссертации приводится способ оценки ошибки метода КМК. Это способ имеет важное практическое значение, так как теоретическая оценка погрешности КМК существует для класса функций ограниченной вариации в виде неравенства Коксмы-Хлавки, но на практике получение численного значения верхней оценки является трудной задачей, как это отмечено в работах Ермакова С.М., Йоханга Ли, Маскани М. и Караваи-новой А.

Целями диссертационной работы являются:

1) разработка модификации методов МК и КМК для решения СЛАУ. Разработанная модификация должна преодолеть некоторые трудности, возникающие при применении методов МК и КМК, в частности, увеличить скорость убывания погрешности и ослабить условие мажорантной сходимости, накладываемое на решаемое уравнение;

2) оптимизация параметров предложенных методов;

3) численное исследование предлагаемых методов, в том числе в случае, когда условие мажорантной сходимости для МК и КМК не выполняется;

4) исследование вопроса применения модифицированных методов МК и КМК для решения СЛАУ совместно с одним из методов увеличения скорости сходимости итерационного процесса (методом верхней релаксации);

5) разработка стохастического метода для экспериментальной оценки погрешности метода КМК.

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

В частности, применялись итерационные методы решения уравнений, аппроксимация уравнений методом конечных разностей, метод верхней релаксации, методы статистического анализа: оценка параметров выборки, проверка статистических гипотез. Использовался пакет Statistica7. Программирование осуществлялось в среде MS Visual С++.

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

Результаты, выносимые на защиту:

1) предложен и обоснован стохастический метод для экспериментальной оценки погрешности методов КМК;

2) предложена и изучена модификация методов МК и КМК для решения СЛАУ, применимая в том числе при нарушении условия мажорантной сходимости для методов МК; указаны методы оптимизации предложенных оценок;

3) доказано, что предложенная модификация уменьшает конструктивную размерность, что упрощает применение квазислучайных:чисел и увеличивает скорость сходимости по сравнению с КМК при решении СЛАУ;

4) показана эффективность применения методов МК и КМК совместно с методом верхней релаксации в ряде случаев;

5) реализованы программы и просчитано большое число примеров, подтверждающих теоретические положения диссертации.

Научная новизна. Все результаты диссертации являются новыми. В частности:

1) предложенный стохастический метод для экспериментальной оценки погрешности методов КМК является новым;

2) предложенная модификация метода КМК для решения СЛАУ является новым методом, применимым в том числе при нарушении условия мажорантной сходимости для методов МК;

3) успешное применение предложенной модификации МК и КМК совместно с методом верхней релаксации при нарушении условия мажорантной сходимости для МК также является новым результатом;

4) впервые новыми методами решены системы линейных уравнений большой размерности.

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

Как показано в работах Ермакова С.М. и Данилова Д.Л., применение методов МК эффективно для решения СЛАУ больших размерностей. Следовательно, именно в этом случае эффективно применение описанных модификаций методов МК и КМК, как альтернатива методу МК.

При этом выгода применения модификаций заключается в следующем:

- ослабляется условие мажорантной сходимости, то есть при решении уравнения X = АХ 4- F условие |Ài(|A|)| < 1 заменяется менее строгим условием |Ai(A)| < 1, где Ai - максимальное собственное число матрицы системы, |— матрица, составленная из модулей компонент матрицы А\

- остаток модифицированного метода КМК убывает быстрее, чем остаток КМК;

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

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

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

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

- 5th St.Petersburg Workshop on Simulation, June 26-July 2, 2005;

- MCQMC 2006 Seventh International Conference on Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing Ulm, Germany Aug 14-18 August, 2006;

- всероссийская конференция по вычислительной математике КВМ-2007, Новосибирск, с 18 по 20 июня 2007;

и опубликованы в статьях [1]-[4] (см. список публикаций автора в конце автореферата) и в материалах конференций.

Структура и объем работы. Диссертация объемом 115 страниц состоит из введения, заключения и четырех глав. Содержит 28 рисунков, 42 таблицы и список цитируемой литературы из 37 наименований.

2. СОДЕРЖАНИЕ РАБОТЫ

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

Описывается научная новизна и практическая значимость работы, сформулированы цели и задачи. Указаны методы исследования, приведена краткая аннотация всех глав диссертации, перечислены основные тезисы, выносимые на защиту.

Основным объектом исследования диссертации является применение МК и КМК для решения систем линейных алгебраических уравнений (СЛАУ).

Описывается актуальность предлагаемой в диссертации модификации МК и КМК для решения СЛАУ.

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

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

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

Проблема оценки погрешности КМК отмечается многими авторами, в частности, она обсуждается в работах Ермакова С.М., Йоханга Ли, Мас-кани М. и Караваиновой А.

Эффективность применения метода КМК зависит от "равномерности" распределения квазислучайных чисел, поэтому квазислучайные числа выбираются из множества точек с малым значением "discrepancy".

Определение 1. Для множества s-мерных точек Y = {Yq,..., Yjv-i}; Yi £ [О, l]s, i = 0,..., N — 1, и для измеримого по Лебегу подмножества Q С [О, l]s локальная "discrepancy" определяется как

где xq есть индикатор подмножества Q.

Определение 2. "Star discrepancy" множества s-мерных точек Y определяется как D*(Y) := sup\D(Q*,Y)\, где Q* пробегает все интер-

Q*

валы вида [0, щ) х ... х [0, us), 0 < щ < 1, i = l,...,s. Имеет место теорема Коксмы-Хлавки:

Теорема 1. Для любой последовательности s-мерных точек типа Y = {Уо,...,Улг—i} С [О, l]s и для любой функции f, имеющей ограниченную вариацию в смысле Харди-Краузе, ошибка метода КМК оценивается следующим образом

/ f(Y)dY\ < V(f)D*(Y),

n=0 [ОД]3

где V(f ) -вариация функции f в смысле Харди-Краузе.

На практике вычислить правую часть неравенства Коксмы-Хлавки в большинстве случаев очень сложно, особенно множитель V(f).

Предлагается следующим стохастический метод оценки погрешности. Пусть в качестве Jo берется оценка КМК

j=1

где Xj - квазислучайные вектора.

Для оценки погрешности строится набор стохастических оценок искомого интеграла, использующих узлы метода КМК (вектора Xj). Для получения одной оценки Ji все узлы Xj, j = 1,...,N сдвигаются на некоторый случайный вектор Ег = (а|, ...,а*), компоненты которого независимы и равномерно распределены на множестве [0,1], и, чтобы остаться в пределах куба [0, l]s, берется дробная часть от полученного результата. Имеет место лемма

N

Лемма 1. Оценка вида J = ^ ^ f{{Xj + S}), где S = (cti,...,as) €

з=1

[0, l]s - случайный вектор, компоненты которого независимы и равномерно распределены на [0,1], Xj = (x{,...,xi) G [0, l]s произвольный вектор, является несмещенной оценка интеграла J f(X)dX, то есть EJ = J.

[одр

В данном случае следует ожидать, что дисперсия а2 оценок J из леммы 1 невелика, так как исходные квазислучайные узлы Xj, j = 1 ,...,N, обладают высокой равномерной распределенностью по своей природе, а операция сдвига на случайный вектор и взятия дробной части {xj + а^} не меняет расстояние между узлами ( в метрике тора).

Таким образом следует ожидать, что построив выборку Ji,..., Jjv из N

независимых реализаций случайной величины J, можно с большой точно__л п

стью оценить величину J = f f(X)dX при помощи оценки J = ^ ^ ,7г,

[0,1]* г= 1

а следовательно, оценить = / f(X)dX — Jo- Также можно no-

lo, i]8

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

Была проведена серия экспериментов для выяснения точности предлагаемой оценки для остатка метода КМК. С помощью метода КМК вычислялись интегралы от функции f(x) = П exp(sin(xfe)cos(:rfc)) для раз-

k=1

мерностей s — 2 и 4 при N = 1296 и 9261000 соответственно. В обоих случаях полученная стохастическая оценка остатка КМК уже при п = 10 имела оценку дисперсии меньшую, чем 6% от оцениваемой погрешности, что говорит о применимости предлагаемой оценки.

Во второй главе описывается модифицированный метод МК для ре-

шения СЛАУ.

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

Чтобы ослабить первое ограничение в диссертации предлагается модификация метода МК. Идея этой модификации впервые появилась в работе Ермакова С.М. и Вагнера В.

Модификация метода МК для решения СЛАУ. Рассматривается задача решения СЛАУ X = АХ + Р, где А = ||ау - матрица СЛАУ, Р = {/¿}"=1 ~ вектор правой части СЛАУ. Допустим, что р(А) < 1 (р -спектральный радиус)

Применение метода МК предполагает оценку суммы ряда Неймана

оо

X — А1^ как интеграла по считающей мере, см. работу Форсит Ж.Е. ¿=о

и Лейблер Р.З.

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

Зафиксируем достаточно большое число КМ вычисляемых слагаемых ряда Неймана, тогда

X = Р + АР + ... + АкР + ек, (1)

где £к - систематическая погрешность, такая, что ЦеьЦ < )] А||ь+1 0

при к I—► оо.

Пусть имеется начальное приближение решения СЛАУ Хо, исходная системах = АХ+Р эквивалентна (Х-Х0) = А(Х-Х0) + (АХ0-Аг0 + Г). Используя метод МК, последовательно вычисляем оценки векторов вида АУ (будем обозначать их АУ):

- г <- АХ0 + Г;

- на первом шаге вычисляется оценка = АУ, где У = АХо — Хо+Р,

- на втором шаге оценивается £2 = АУ, где У = £ 1 - оценка, полученная на первом шаге, Z <— 2 +

- на шаге к оценивается £ = АУ, где У = - оценка, полученная на шаге к-1, 2 2 +

В качестве оценки X берется 2. Один из вариантов оценки для вычисления АУ приводится ниже.

Выбор оптимальных параметров цепи Маркова. Существуют различные способы решения СЛАУ X = АХ+Р методом МК, которые основаны на моделировании траекторий цепей Маркова. При этом возникает вопрос об оптимизации параметров этих цепей Маркова в зависимости от параметров исходной системы. В зависимости от задачи, оптимальность может пониматься как минимизация суммы дисперсий оценок компонент

п ______л

^ ВХг или как минимизация дисперсии скалярного произведения (X, Н),

г=1

где Н заданный вектор. Оптимальные параметры для прямых и сопряженных оценок по поглощению и столкновениям могут быть найдены в работах С.М.Ермакова (например, "Методы Монте-Карло и смежные вопросы"). Предпочтительный вид используемой оценки зависит от конкретной задачи.

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

На каждом шаге метода оценивается величина вида АУ, У = {у±,..., уп). Фиксируем некоторое распределение Р° = ||р°||"=1 и стохастическую матрицу Р = \\pij (ее строки задают распределения). Пусть а — случайная величина с распределением Р°, ¡3 — случайная величина с распределением Р{(3 = ]\а = г) — р^, ¿,] = 1,...,п. Справедливо утверждение: Утверждение 1. Рассмотрим оценку задаваемую выражением

с _ У а ара

? = "ТГ-е/Э> I2)

Ра

где е^ — 1-й единичный орт размерности п. Пусть выполняются условия

п

согласования: если а^у^ ф 0, то р^ > 0 и если у3 ^ аи Ф О, то р® > 0. Тогда Е£ = АУ.

Возникает вопрос выбора оптимальных параметров Р° и Р. Если ПО-

71

нимать оптимальность в смысле минимизации выражения т0 тео~

¿=1

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

Теорема 2. При выполнении условий согласования выражение ^

¿=1

достигает своего минимума при выборе параметров

п п п

p°j = 1С \aikVk\, Pji = Wkjl

¿=1 i,k=1 fe=l

Вопрос о сходимости модифицированного метода МК рассматривался в работах Ермакова С.М. и Вагнера В. Более подробно для случая решения СЛАУ вида X = АХ + F этот вопрос рассматривается в диссертации. В частности, сформулирована и доказана теорема о достаточных условиях ограниченности вторых моментов строящихся оценок модифицированного метода Монте-Карло. Она показывает, что при выборе параметров из приведенной теоремы 2 и при усреднении на каждом шаге ММК достаточно большого числа оценок вида (2) имеет место ограниченность вторых моментов оценки решения (а следовательно, процесс расчета стохастически устойчив).

Примеры применения модифицированного метода МК с оптимизированными параметрами для решения СЛАУ, в том числе при нарушении мажорантной сходимости для метода МК, приводятся в главе 4.

В третьей главе предлагается применить описанную выше модификацию метода МК к методу КМК.

В некоторых случаях использование методов КМК, для которых порядок убывания остатка есть logs N/N, с ростом N оказывается более эффективным по сравнению с методами МК, где убывание остатка имеет порядок 1/v^V, где N - число моделируемых траекторий для получения одной оценки, s - конструктивная размерность оцениваемого интеграла. Этот факт описывается в работах Дж. Холтона, С.М.Ермакова, Маскани М. (при этом рассматривались примеры, когда значение s мало).

Несмотря на эффективность КМК из-за относительно высокой скорости убывания его погрешности, применение КМК для решение СЛАУ сталкивается с некоторыми проблемами. Во-первых, существует необходимость выполнения условия мажорантной сходимости, унаследованная от метода МК. Во-вторых, оценка большого количества слагаемых ряда Неймана или суммы ряда приводит к необходимости использовать квазислучайные точки большой размерности s. С ростом s сильно увеличивается значение дГ N), что приводит к медленной сходимости метода. В главе 3 предлагается модифицированный метод КМК (в дальнейшем МКМК). Идея модификации описана в главе 2 для метода МК. В случае КМК отличие состоит в использовании квазислучайных чисел вместо псевдослучайных.

Применительно к методу КМК предложенная модификация позволяет,

во-первых, повысить скорость сходимости оценок КМК, и, во-вторых, как и для метода МК позволяет ослабить условие мажорантной сходимости, так как суммирование отрезка ряда Неймана происходит в определенном порядке. Конструктивная размерность оцениваемых интегралов при этом постоянна и равна 1 (в = 1).

Оценка погрешности метода МКМК. При применении МКМК для решения СЛАУ X = АХ + Р с начальным приближением Хо используем числа Холтона вида рг(п), г - простое число, п ~ натуральное.

Используя метод КМК, последовательно вычисляем оценки вида АУ:

- 2<-АХйЛ-Р\ _

на первом шаге вычисляется оценка £1 = АУ, где У = АХо — Хо + Р, 2 2 Л-1, параметры используемых чисел Холтона гбйь п £ Д; на втором шаге оценивается £2 = АУ, где У = £1, 2 2 + £2> параметры используемых чисел Холтона г £ Я2, п & /2; на шаге к оценивается £ = АУ, где У = ^-ъ 2 2 + параметры используемых чисел Холтона г € п € В качестве оценки X берется 2. Д1, ...,Як обозначают непересекающиеся между собой группы простых чисел, Д,...,/^ - непересекающиеся между собой группы натуральных чисел.

Рассмотрим погрешность модифицированного метода КМК при Х0 = 0. Метод предполагает вычисление последовательности оценок £п+х = АЕп, где Еп есть оценка Оп = Хп - Хп_ь Хп+1 - АХп +Р,Ег = Г.

Положим Е„ — £)„ + еп, тогда £1 = 0.

Утверждение 2. Для ошибки одного шага модифицированного метода Квази Монте-Карло справедливо рекуррентное соотношение = {А + 5п) еп + 6пОп, где 5п - детерминированная матрица, норма которой, в силу характера убывания погрешности в методе КМК, имеет порядок ||<У = причем выполняется оценка

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

Утверждение 3. Справедливо рекуррентное соотношение:

||Д*-Хк||

Сравнивая модифицированный КМК с другими методами, можно утверждать, что он точнее, чем модифицированный метод МК в О (ц^) и точнее, чем

обычный КМК в 0( 1п3-1 Ы) раз, где й средняя длина траектории в методе КМК.

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

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

В четвертой главе модифицированные методы МК и КМК успешно применяются совместно с одним из методов улучшения сходимости итерационного процесса - методом верхней релаксации.

В качестве одного из примеров решается СЛАУ, полученная после аппроксимации на сетке двумерного уравнения Пуассона AU = х{1 —х)у{ 1 — у); х,у £ [0,1] с нулевыми граничными условиями. Квадрат [0,1] х [0,1] покрыт прямоугольной сеткой размерности Nnet х Nnet, h = l/(Nnet — 1) - шаг сетки. Применяется метод верхней релаксации:

Un+1 = Un-T{Un-AUn-F), (3)

где Un,F - векторы размерности N%et, А - матрица размерности N%et х N%et и выбирается оптимальное значение параметра г* = На этом примере демонстрируется, что применение метода МКМК совместно с аналогом метода верхней релаксации эффективно. К тому же, построенная стохастическая схема позволяет эффективно распараллеливать алгоритм.

В качестве другого примера в четвертой главе рассматривается разностный аналог уравнения Навье-Стокса, описывающего в декартовых координатах движение вязкой несжимаемой жидкости в трехмерном случае:

dv 1 d2ViVk ...

— + (vV)v = — gradp + vAv, А (4)

i,k

Здесь v — векторное поле скорости, р — давление (скаляр). Константы р и v определяют плотность и кинематическую вязкость соответственно.

Выбранный для расчетов пример относится к задаче движения жидкости в кубической каверне с подвижной стенкой с прилипанием жидкости

на границах. Скорость стенки равна 1 — 1/(1 + ¿)5, £ > 0. В начальный момент жидкость находится в состоянии покоя, и = 0.0025.

Рассматривалась прямоугольная ортогональная сетка [15x15x15] и неявная разностная схема второго порядка точности во внутренних точках. На каждом шаге правая часть сеточного уравнения линеаризовыва-лась и решалось СЛАУ вида X — АХ + Р. Размерность системы п = 153 • 3 = 10125 (множитель 3 - из-за трехмерной размерности).

На интервале времени £ = [0,0.8) оценка решения сеточного уравнения находилась детерминистическим образом, а далее, для расчета очередного шага неявного метода Эйлера, решение системы линейных алгебраических уравнений осуществлялось различными методами КМК.

Сначала для решения этой системы применялся метод КМК с использованием прямой оценки по поглощению. Вероятности перехода в используемой цепи Маркова ру, г, j = 1 ,...,п выбирались пропорциональными модулям коэффициентов ац системы. Вероятность поглощения в каждой точке равна 0.1. Средняя длина траектории равнялась 10.

Путем усреднения некоторого числа оценок метода КМК строилась оценка X решения СЛАУ. Для исследования ее характеристик строилась выборка из 10 независимых реализаций Х{, г = 1,..., 10. Затем рассчитывалась выборка — Х\\2, г = 1,..., 10, где || • ||г обозначает евклидову норму, X — точное решение СЛАУ. Далее для решения той же СЛАУ применялся модифицированный метод КМК и строились аналогичные оценки. Результат приводится в таблице 1, которая показывает зависимость нормы ошибки решения СЛАУ (строки КМК и МКМК содержат оценки математического ожидания ||Х — Х||г для соответствующего метода) от числа N использованных чисел Холтона. Строку N можно считать характеристикой трудоемкости метода.

Таблица 1 - погрешности оценок решения СЛАУ при решении уравнения Навье-Стокса_

repeats 104 2 • 104 4 • 104 8 • 104 16 • 104 32 • 104 64 • 104

КМК 8.29 5.93 4.05 2.97 2.11 1.46 1.04

МКМК 0.84 0.83 0.65 0.52 0.41 0.3 0.24

Из таблицы 1 видно, что в данном примере при фиксированной трудоемкости метода, понимаемой как число использованных чисел Холтона, метод МКМК является более точным, чем метод КМК.

Третий пример посвящен решению различных СЛАУ со случайно заполненными коэффициентами. В диссертации приводятся результаты расчета большого числа примеров. Например, рисунок 1 иллюстрирует зависимость оценки евклидовой нормы погрешности от числа обращений к датчику псевдо или квазислучайных чисел для различных методов при

решении СЛАУ размерности 100. Масштаб осей логарифмический.

I.S0Î

¡Ж

1

3

1-МК 2 - MMR" 4j3- КМК

ж 1 seas гит <тж 4 - МКМК

АЯИ Säüds 1024505

Рис. 1 — сравнение погрешностей при решении СЛАУ размерности 100

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

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

[1] Ермаков С.М.,Рукавишникова А.И. Квази Монте-Карло алгоритмы решения систем оинейных алгебраических уравнений// Мат. модели. Теория и приложения под ред. Чиркова 2005. Выпуск 6, стр.3-27.

[2] Ermakov S.M., Rukavishnikova A.I. Quasi Monte-Carlo Algorithms for Solving Linear Algebraic Equation// MC Methods and Applications. 2006. Vol.12, №5, pp.363-384.

[3] Рукавишникова А.И. Оптимизация алгоритма Квази Монте-Карло решения систем линейных алгебраических уравнений// Вестник С-Петербургского университета, сер.1, вып.1, стр.73-78, 2008.

[4] Ермаков С.М., Тимофеев К.А., Рукавишникова А.И. О некоторых стохастических и квазистохастических методах решения уравнений// Вестник С-Петербургского университета, сер.1, вып.4, стр.75-85, 2008.

СПИСОК ПУБЛИКАЦИИ АВТОРА

Подписано в печать 25.12.2008г. Формат 60x84 1/16. Бумага офсетная. Печать офсетная. Усл. печ. л. 1,0. Тираж 100 экз. Заказ № 1025.

Отпечатано в ООО «Издательство "JIEMA"»

199004, Россия, Санкт-Петербург, В.О., Средний пр., д.24, тел./факс: 323-67-74 e-mail: izd_lema@mail.ru http://www.lemaprint.ru

Введение

1 Статистические свойства квазислучайных чисел, оценка погрешности метода КМК

1.1 Основные определения и факты.

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

1.3 Оценка погрешности метода КМК.

1.3.1 Пример

1.3.2 Пример 2.

1.3.3 Пример 3.

2 Решение СЛАУ модифицированным методом МК, оптимизация оценок модифицированного метода

2.1 Решение СЛАУ методом МК.

2.2 Применение марковских цепей для решения СЛАУ.

2.3 Параллелизм метода МК и его сравнительная сложность

2.4 Предлагаемая модификация метода МК.

2.5 Выбор оптимальных параметров оценок при использовании модифицированного МК.

2.6 Устойчивость модифицированного МК

3 Решение СЛАУ модифицированным методом КМК

3.1 Применение модифицированного КМК и оценка его погрешности.

3.2 Эксперименты по сравнению методов МК, модифицированного МК, КМК и модифицированного КМК.

4 Численные эксперименты

4.1 Метод последовательной верхней релаксации.

4.2 Применение к уравнению Лапласа.

4.3 Устойчивость модифицированного МК, примененного совместно с методом верхней релаксации.

4.4 Применение модифицированного КМК при ослаблении условия мажорантной сходимости.

4.5 Решение уравнения Навье-Стокса.

Методы Монте-Карло(МК) и Квази Монте-Карло(КМК) являются одним из важнейших инструментов для решения задач в многочисленных областях физики, математики, экономики, оптимизации, теории управления и др. В том числе, см. [1]-[4],[7], [11]-[15], [19]-[22], [33]-[34], они применяются для решения основной задачи линейной алгебры - нахождения решения системы линейных алгебраических уравнений (СЛАУ).

Интерес к этим методам обусловлен, в частности, тем, что как показано в работах Ермакова С.М. и Данилова Д.Л.([2],[3]), при большой размерности СЛАУ метод МК обладает меньшей трудоемкостью (в статистическом смысле), чем детерминированные итерационные методы. Также, большое преимущество перед другими методами дает методу МК его естественное свойство параллелизма. В некоторых случаях методы МК и КМК допускают неограниченный параллелизм и, следовательно, полную загрузку оборудования, что особенно важно для решения задач большой размерности, см. [1], [5].

В диссертации исследовались способы применения методов МК и КМК для решения систем линейных алгебраических уравнений (СЛАУ). Процесс решения СЛАУ методами МК и КМК может быть связан с моделированием марковских цепей с дискретным временем. В последнее время такого рода подход применялся во многих работах, например [2], [3], [7], [11], [12], [15]. В этих работах результаты были получены при наличии ряда ограничений, накладываемых на класс уравнений, которые могут быть решены с помощью методов МК и КМК.

Одно из ограничений связано с необходимостью выполнения условия мажорантной сходимости. Оно было частично преодолено в работах Ермакова С.М. и Вагнера В. ([4],[5]) и в работе Сабельфельда К.К. и Шалимовой И.А.([13]).

Другое ограничение связано с невысоким порядком убывания погрешности метода МК при линейном увеличении времени расчета. Для достижения точности порядка е требуется проделать порядка 1/е2 операций. Для преодоления этой проблемы применяют методы КМК, поскольку при оценке э-кратных интегралов они в 0( точнее методов МК (ТУ - число траекторий, промоделированных для получения одной оценки).

Применение методов КМК для решения СЛАУ приводит к возникновению новой проблемы, а именно, к необходимости использовать квазислучайные векторы большой размерности 5 для оценки высокой степени матрицы оператора. Это приводит к ухудшению сходимости из-за наличия множителя 1п5 N в числителе оценки скорости убывания остатка КМК. Проблема отмечается многими авторами ([1], [27], [28], [34]).

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

Приводится способ оценки ошибки метода КМК. Это способ имеет важное практическое применение, так как теоретическая оценка погрешности КМК существует для класса функций ограниченной вариации в виде неравенства Коксмы-Хлавки, но на практике получение численного значения верхней оценки является трудной задачей, как это отмечено в работах [1], [12], [13], [14].

1) разработка модификации методов МК и КМК для решения СЛАУ. Разработанная модификация должна преодолеть некоторые трудности, возникающие при применении КМК, в частности, увеличить скорость убывания погрешности метода КМК и ослабить условие мажорантной сходимости, накладываемое на решаемое уравнение;

2) оптимизация параметров предложенных методов;

3) численное исследование предлагаемых методов, в том числе в случае, когда условие мажорантной сходимости для МК и КМК не выполняется;

4) исследование вопроса применения модифицированных методов МК и КМК для решения СЛАУ совместно с одним из методов увеличения скорости сходимости итерационного процесса;

5) разработка стохастического метода для экспериментальной оценки погрешности метода КМК.

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

Описывается одна из сложностей в применении метода КМК, а именно тот факт, что трудно получить теоретическое значение оценки погрешности метода Длг[/], где

Предлагается и детально описывается стохастический метод для экспериментальной оценки погрешности метода КМК. Теоретические предположения подтверждаются серией проведенных экспериментов.

Вторая глава посвящена модифицированному методу МК для решения СЛАУ. Как отмечалось ранее, существуют определенные ограничения, налагаемые на класс задач, которые могут быть решены путем моделирования траекторий случайных процессов. Во-первых, это естественные ограничения, связанные с возможностью представления решения в виде интеграла по траекториям. Во-вторых, это ограничения, связанные с известным невысоким порядком убывания погрешности стохастических методов - для достижения точности порядка е требуется проследить порядка 1/е2 траекторий.

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

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

В третьей главе предлагается применить описанную во второй главе модификацию метода МК для метода КМК. Исследуется оценка погрешности модифицированного метода КМК. Приводятся результаты экспериментов по сравнению методов МК, КМК, модифицированного МК и модифицированного КМК для решения СЛАУ. Делается вывод, что по сравнению с методами

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

В четвертой главе приводятся численные эксперименты.

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

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

Заключение

В результате проведенных исследований в диссертации были получены следующие важные результаты:

1) Предложен и обоснован стохастический метод для экспериментальной оценки погрешности методов КМК

2) Предложена и изучена модификация методов МК и КМК для решения СЛАУ, применимая в том числе при нарушении условия мажорантной сходимости для методов МК. Указаны методы оптимизации предложенных оценок.

3) Доказано, что предложенная модификация уменьшает конструктивную размерность, что упрощает применение квазислучайных чисел и увеличивает скорость сходимости по сравнению с КМК при решении СЛАУ.

4) Показана эффективность применения методов МК и КМК совместно с аналогом метода верхней релаксации в ряде случаев.

5) Реализованы программы и просчитано большое число примеров, подтверждающих теоретические положения диссертации.

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

1) предложенный стохастический метод для экспериментальной оценки погрешности методов КМК является новым;

2) предложенная модификация метода КМК для решения СЛАУ является новым методом, применимым в том числе при нарушении условия мажорантной сходимости для методов МК;

3) успешное применение предложенной модификации МК и КМК совместно с аналогом метода верхней релаксации при нарушении условия мажорантной сходимости для МК также является новым результатом;

4) впервые новыми методами решены системы линейных уравнений большой размерности.

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

2. Ermakov S.M., Danilov D.L. On the comparative complexity of the Monte-Carlo method for solving systems of linear equations// Comput. Math. Mat. Phys., №5, pp.35-40, 1995.

3. Ермаков С.M., Данилов Д.JI. Асимптотическая сложность оценки по столкновениям для решения линейных систем// Журнал Вычислительной Математики и Мат. Физики, №5, т.37, стр.515-523, 1997.

4. Вагнер В., Ермаков С.М. Monte-carlo difference schemes for the wave equations // Monte-Carlo Methods and Appl., Vol.18, pp. 1-19, 2002.

5. Вагнер В., Ермаков С.М. Стохастическая устойчивость и параллелизм метода Монте-Карло // ДАН. 2001, 179 №4., стр.438-441, 2001.

6. Д. Кнут, Э.Яо Сложность моделирования неравномерных распределений // Кибернетический сборник, Новая серия, вып. 19, М., Мир, 1983, стр.97-158

7. Danilov D.L., Ermakov S.M., Haïton J.H. Asymptotic complexity of Monte Carlo methods for solving linear systems. //J.Statistical Planning and Inference, 2000, p.p. 5-18.

8. Ермаков С.М. Дополнение к одной работе по методу Монте-Карло, Ж.Выч. Математики и Мат. Физики, 2001, т.41, н.б.

9. Адамов А.В., Ермаков С.М. О стохастической устойчивости метода Монте-Карло(случай операторов)// Вестник С-Петербургского университета, сер.1, вып.2, стр.3-7, 2004.

10. Ermakov S.M. Neumann-Ulam scheme and particle methods // IV-th IMACS Seminar on Monte-Carlo Methods: Abstracts(SEptember 15-19, 2003, Berlin), Berlin, WIAS, p.55,2003.

11. Forsith G.E. Leibler R.Z. Matrix inversion by a Monte-Carlo methods // Math tables and other aids to computation,4, 1950, pp.127-120.

12. Maskagni M, Karavainova A. A Parallel QMC Method for Solving of Linear Equations// Amsterdam, Lecture Notes in Computer Science, vol.2330, pp.598-608, 2008.

13. Maskagni M, Karavainova A. Matrix Computation Using Quasirandom Sequences// Amsterdam, Lecture Notes in Computer Science, Springer, vol.1988, pp.552-559, 2001.

14. Maskagni M, Yaohang Li. Grid Based Quasi-Monte-Carlo Applications// Monte-Carlo Methods and Appl., Vol.11, pp.39-55, 2005.

15. Dimov I., Aleksandrov V., Karaivanova Resolvent Monte-Carlo Methods for Linear Algebra Problems// Mathematics and Computers in Simulation, Vol.55, pp.25-36, 2001.

16. Sabelfeld, K, Shalimova I. Random walk on spheres methods for iterative solution of elastisity problems// Monte-Carlo Methods and Appl., Vol.8, pp.171-202, 2002.

17. Сабельфелъд К.К. Методы Монте-Карло в краевых задачах// Новосибирск, "Наука", 280с., 1989 The Jons Hopkins Univ. Press,Baltimore, 1996

18. Sabelfeld K, Shalimova I., Levikin A. Random walk on fixed spheres methods for Laplace and Lame equations// Monte-Carlo Methods and Appl., Vol.12, pp.55-93, 2006.

19. Ермаков С.M.,Рукавишникова А.И. Квази Монте-Карло алгоритмы решения систем линейных алгебраических уравнений// Мат. модели. Теория и приложения под ред. Чиркова 2005. Выпуск 6, стр.3-27.

20. Ermakov S.M., Rukavishnikova A.I. Quasi Monte-Carlo Algorithms for Solving Linear Algebraic Equation// MC Methods and Applications. 2006. Vol.12, №5, pp.363-384.

21. Рукавишникова А.И. Оптимизация алгоритма Квази Монте-Карло решения систем линейных алгебраических уравнений// Вестник С-Петербургского университета, сер.1, вып.1, стр.73-78, 2008.

22. Ермаков С.М., Тимофеев К.А., Рукавишникова А.И. О некоторых стохастических и квазистохастических методах решения уравнений// Вестник С-Петербургского университета, сер.1, вып.4, стр.75-85, 2008.

23. Товстик Т.М. Сравнение некоторых статистических свойств квазислучайных и псевдослучайных последовательностей// Вестник С-Петербургского университета, сер.1, вып.2, стр.62-70, 2006.

24. G.H. Golub, C.F. Van Loon Matrix Computations// The Jons Hopkins Univ. Press, Baltimore, 1996.

25. Соболь И.М. Многомерные квадратурные формулы и функции Хаара// М., Наука, 288с., 1969

26. J.H. Halton Sequential Monte-Carlo Techniques for the Solution of Linear Systems // SIAM Journal of Scientific Computing, Vol.9,pp213-257,1994

27. Halton J.H. On the efficiency of certain quasi-random sequences of points in evaluating multidimensional integrals// Numer. Math., vol.2, pp.84-90, 1960

28. Halton J.H. Quasi-probability // Monte-Carlo Methods and Appl., vol.11, №3, pp.203-350.

29. Михайлов Г. А. Войтишек А.В. Численное статистическое моделирование. Методы Монте-Карло // М., Академия, 367с, 2006.

30. Ripley B.D. Stochastic Simulation // Wiley, New-York, 1987

31. Ермаков C.M. Михайлов Г.А. Статистическое моделирование // M., Наука, 2-е изд., доп., 1982.

32. Hammersley Т.М. Handscomb D.C. Monte-Carlo methods // John Wiley and Sons, N.Y., London, Sydney, Methuen, 1964.

33. Niederreiter H. Random Number Generation and Quasi-Monte Carlo Methods // SIAM, Phil., Pensilvania, 1992.

34. Темам P. Уравнения Навье-Стокса // M., Мир, 1981.

35. Марчук Г.И. Методы вычислительной математики // М., Наука, 456с., 1977.

36. G.H.Golub, C.F. Van Loon Matrix Computations// The Jons Hopkins Univ. Press,Baltimore, 1996