Содержание к диссертации
Введение
1 Статистические свойства квазислучайных чисел, оценка погрешности метода КМК 9
1.1 Основные определения и факты 9
1.2 Экспериментальное исследование стохастических свойств квазислучайных чисел 17
1.3 Оценка погрешности метода КМК 26
1.3.1 Пример 1 31
1.3.2 Пример 2 34
1.3.3 Пример 3 36
2 Решение СЛАУ модифицированным методом МК, оптимизация оценок модифицированного метода 39
2.1 Решение СЛАУ методом МК 39
2.2 Применение марковских цепей для решения СЛАУ 41
2.3 Параллелизм метода МК и его сравнительная сложность . 49
2.4 Предлагаемая модификация метода МК 53
2.5 Выбор оптимальных параметров оценок при использовании модифицированного МК 56
2.6 Устойчивость модифицированного МК 60
3 Решение СЛАУ модифицированным методом КМК 65
3.1 Применение модифицированного КМК и оценка его погрешности 65
3.2 Эксперименты по сравнению методов МК, модифицированного МК, КМК и модифицированного КМК 70
4 Численные эксперименты 90
4.1 Метод последовательной верхней релаксации 90
4.2 Применение к уравнению Лапласа 93
4.3 Устойчивость модифицированного МК, примененного совместно с методом верхней релаксации 95
4.4 Применение модифицированного КМК при ослаблении условия мажорантной сходимости 96
4.5 Решение уравнения Навье-Стокса 104
Заключение 110
Список литературы 112
- Экспериментальное исследование стохастических свойств квазислучайных чисел
- Параллелизм метода МК и его сравнительная сложность
- Эксперименты по сравнению методов МК, модифицированного МК, КМК и модифицированного КМК
- Применение модифицированного КМК при ослаблении условия мажорантной сходимости
Введение к работе
Методы Монте-Карло(МК) и Квази Монте-Карло(КМК) являются одним из важнейших инструментов для решения задач в многочисленных областях физики, математики, экономики, оптимизации, теории управления и др. В том числе, см. [1]-[4],[7], [11]-[15], [19]-[22], [33]-[34], они применяются для решения основной задачи линейной алгебры - нахождения решения системы линейных алгебраических уравнений (СЛАУ).
Интерес к этим методам обусловлен, в частности, тем, что как показано в работах Ермакова СМ. и Данилова Д.Л.([2],[3]), при большой размерности СЛАУ метод МК обладает меньшей трудоемкостью (в статистическом смысле), чем детерминированные итерационные методы. Также, большое преимущество перед другими методами дает методу МК его естественное свойство параллелизма. В некоторых случаях методы МК и КМК допускают неограниченный параллелизм и, следовательно, полную загрузку оборудования, что особенно важно для решения задач большой размерности, см. [1], [5].
В диссертации исследовались способы применения методов МК и КМК для решения систем линейных алгебраических уравнений (СЛАУ). Процесс решения СЛАУ методами МК и КМК может быть связан с моделированием марковских цепей с дискретным временем. В последнее время такого рода подход применялся во многих работах, например [2], [3], [7], [11], [12], [15]. В этих работах результаты были получены при наличии ряда ограничений, накладываемых на класс уравнений, которые могут быть решены с помощью
методов МК и КМК.
Одно из ограничений связано с необходимостью выполнения условия мажорантной сходимости. Оно было частично преодолено в работах Ермакова СМ. и Вагнера В. ([4],[5]) и в работе Сабельфельда К.К. и Шалимовой И.А.([13]).
Другое ограничение связано с невысоким порядком убывания погрешности метода МК при линейном увеличении времени расчета. Для достижения точности порядка є требуется проделать порядка І/є2 операций. Для преодоления этой проблемы применяют методы КМК, поскольку при оценке s-кратных интегралов они в 0(^) точнее методов МК (7V - число траекторий, промоделированных для получения одной оценки).
Применение методов КМК для решения СЛАУ приводит к возникновению новой проблемы, а именно, к необходимости использовать квазислучайные векторы большой размерности s для оценки высокой степени матрицы оператора. Это приводит к ухудшению сходимости из-за наличия множителя lns N в числителе оценки скорости убывания остатка КМК. Проблема отмечается многими авторами ([1], [27], [28], [34]).
В диссертации предлагается новый вариант модификации метода КМК (в дальнейшем МКМК), основанный на идее снижения конструктивной размерности оцениваемых интегралов. Предложенная модификация позволяет, во-первых, повысить скорость сходимости оценок КМК, и, во-вторых, позволяет ослабить условие мажорантной сходимости, так как суммирование отрезка ряда Неймана происходит в определенном порядке.
Приводится способ оценки ошибки метода КМК. Это способ имеет важное практическое применение, так как теоретическая оценка погрешности КМК существует для класса функций ограниченной вариации в виде неравенства Коксмы-Хлавки, но на практике получение численного значения верхней оценки является трудной задачей, как это отмечено в работах [1], [12], [13], [14].
Основными целями диссертации являются:
разработка модификации методов МК и КМК для решения СЛАУ. Разработанная модификация должна преодолеть некоторые трудности, возникающие при применении КМК, в частности, увеличить скорость убывания погрешности метода КМК и ослабить условие мажорантной сходимости, накладываемое на решаемое уравнение;
оптимизация параметров предложенных методов;
численное исследование предлагаемых методов, в том числе в случае, когда условие мажорантной сходимости для МК и КМК не выполняется;
исследование вопроса применения модифицированных методов МК и КМК для решения СЛАУ совместно с одним из методов увеличения скорости сходимости итерационного процесса;
разработка стохастического метода для экспериментальной оценки погрешности метода КМК.
Первая глава диссертации посвящена исследованию статистических свойств квазислучайных чисел и оценке погрешности метода КМК. Приводится краткое сравнение свойств псевдослучайных и квазислучайных чисел. Для квазислучайных последовательностей формально применяется терминология случайных чисел. Исследуются выборочные функции распределения, коэффициенты корреляции, кросскорреляционые функции для последовательностей Холтона и Соболя.
Описывается одна из сложностей в применении метода КМК, а именно тот факт, что трудно получить теоретическое значение оценки погрешности
метода i2jv[/], где
Г N
RnU}= J f(X)dX-Y^AjfiX-).
[одр j=1
Предлагается и детально описывается стохастический метод для экспериментальной оценки погрешности метода КМК. Теоретические предположения подтверждаются серией проведенных экспериментов.
Вторая глава посвящена модифицированному методу МК для решения СЛАУ. Как отмечалось ранее, существуют определенные ограничения, налагаемые на класс задач, которые могут быть решены путем моделирования траекторий случайных процессов. Во-первых, это естественные ограничения, связанные с возможностью представления решения в виде интеграла по траекториям. Во-вторых, это ограничения, связанные с известным невысоким порядком убывания погрешности стохастических методов - для достижения точности порядка є требуется проследить порядка І/є2 траекторий.
Чтобы ослабить первое ограничение, в диссертации предлагается модификация метода МК. Идея этой модификации впервые появилась в работе Ермакова СМ. и Вагнера В. Предложенная модификация позволяет ослабить условие мажорантной сходимости, так как суммирование отрезка ряда Неймана происходит в определенном порядке.
Автору диссертации принадлежит применение этой модификации для решения СЛАУ, выбор оптимальных параметров цепи Маркова и исследование стохастической устойчивости модифицированного метода при решении СЛАУ
В третьей главе предлагается применить описанную во второй главе модификацию метода МК для метода КМК. Исследуется оценка погрешности модифицированного метода КМК. Приводятся результаты экспериментов по сравнению методов МК, КМК, модифицированного МК и модифицированного КМК для решения СЛАУ Делается вывод, что по сравнению с методами
МК и КМК, предложенная модификация позволяет повысить скорость сходимости оценок КМК, то есть частично преодолеть второй недостаток МК.
В четвертой главе приводятся численные эксперименты.
Модифицированные методы МК и КМК успешно применяются совместно с одним из методов улучшения сходимости итерационного процесса - аналогом метода верхней релаксации.
Описанный в диссертации метод МКМК наряду с методом КМК применяются для решения уравнения Навье-Стокса, описывающего в декартовых координатах движение вязкой несжимаемой жидкости в трехмерном случае. Делается вывод об эффективности предложенного в диссертации МКМК.
Экспериментальное исследование стохастических свойств квазислучайных чисел
Тогда функция 1/N ]Г) &(ХІ—t) по определению является эмпирической функ 1=1 цией распределения случайной величины, имеющей своими реализациями название "discrepancy" и принимается в теории чисел в качестве меры неравномерности распределения множества точек х\%..., х . Предполагается 0 Xi 1. Рассмотрим квадратурную формулу с равными коэффициентами. Легко проверить, что, если функция / имеет суммируемую производную / , то для нее справедлива точная оценка В многомерном случае имеет место похожий результат, называемый неравенством Коксмы-Хлавки, который показывает, что величина Dj /N, где является нормой остатка в классе функций ограниченной вариации в смысле Харди-Краузе. Неравенство Коксмы-Хлавки означает, что при оценке интеграла с помощью набора квазислучайных последовательностей, ошибка не будет превосходить усредненной суммы "discrepancy" этих наборов, умноженной на вариацию подынтегральной функции в смысле Харди-Краузе. Итак, подобрав последовательность векторов X\t..., Xjy с известным значением "dicsrepancy", можно судить о норме остатка кубатурной формулы с равными коэффициентами. Понятие "discrepancy" часто встречается в литературе. Например, в математической статистике величина "discrepancy"(в одномерном случае), деленная на N, фигурирует под названием критерия согласия Колмогорова (для равномерного распределения). В теории чисел Соболя, см. [25], величина "discrepancy" называется отклонением и вводится через теорию сеток. Определение 1.1.2 Равномерной сеткой So называется сетка, состоящая из N = Мп точек с координатами (Гі , -, г"м ) где гЪ— гп пробегают независимо друг от друга все целые значения от 0 до М — 1, величины Pv, i/ = l,..., п фиксированы, 0 Д, 1. Определение 1.1.3 Отклонением сетки {XQ,XI,...,XN-I} называется величина Ду := D(XQ, ..., Величина DM В каком-то смысле характеризует отклонение распределения точек сетки от равномерного. Рассмотрим так называемую последовательность векторов Холтона — Xi, ... ,.XJV G M.s, для которой RN[I] 0(lnsN/N). Впервые эта последовательность была проанализирована Дж. Холтоном в 1960, см. [28]. Приведем известный алгоритм получения чисел Холтона [28]. Пусть п, г -натуральные числа. Тогда число п может быть единственным образом пред-3 ставлено в виде где для V і Мг — 1 выполняется 0 щ г — 1, а также 1 пмг г — 1, Mr = \logrn\ = [log п/ log г J С помощью, так называемого, "отражения" представления п, получаем число рг{п) — (О щпіп2 ... nMr)r = Щг 1 + щг 2 + п2г 3 + ... + пМгг Мт г. Одномерная последовательность Холтона длины N имеет вид pr(n\),...,рг(пн). При фиксированном г и щ числа Рг(ті), п = 1,2,... сначала возрастают и достигают максимума при п = г — 1. Далее при i — m имеет место резкий спад и монотонный подъем до очередного максимума с длиной серии равной г. Для заполнения интервала [0,1] число моделируемых точек N должно быть не менее г. Свойства последовательности проявляются лучше, если N т\ Многомерной последовательностью Холтона называется последовательность векторов в единичном гиперкубе Ds: Х\ = Х(п),Х2 = Х(п+1),.. ,,XN — Х(п + N — 1),..., такая, что при каждом фиксированном п = 1,2,... где гі, r2,..., rs - попарно взаимно простые числа. В работе [28] для "discrepancy" последовательностей Холтона доказана следующая теорема. Рассмотрим сетку в гиперкубе Ds, состоящую из начальных участков последовательностей Холтона (при построении ХІ,І 1 полагаем п — 0) Е = {Xi, ...,XN}. Приведем кратко теорию чисел Соболя [25]. Изначально вводится понятие Ро-сетки и LPQ последовательности. Определение 1.1.4 Одномерная сетка, состоящая из N = 2V точек, называется Ро-сеткой, если каждому двоичному отрезку с длиной 1/N ([ 7, г], к = О,..., 2" — 11, принадлежит одна точка сетки. Определение 1.1.5 Последовательность XQ, Ж І, жг,..., хп,... называется LPQ последовательностью, если любой ее двоичный участок есть Ро-сетка. Где двоичный участок - множество членов с номерами {г : &2J г (к -f 1)2» , AT = 0,1,2,...; j = 1,2,...}
Параллелизм метода МК и его сравнительная сложность
В результате проведенных экспериментов, был замечен интересный факт сохранения моментов четных порядков при моделировании с использованием чисел Холтона нормального распределения. Проиллюстрируем этот эффект на следующем примере.
Ниже приводится алгоритм моделирования s-мерного нормального распределения при заданных ковариационной матрице С и векторе средних at. Многомерная плотность нормального распределения, которую необходимо промоделировать, выглядит следующим образом: Представим матрицу ковариаций в виде произведения двух треугольных матриц С = GG матрицу G можно найти методом квадратных корней из уравнения GGT = С. Сделаем замену X — GY+ ct, якобиан замены равен det G\. Получаем новое выражение для моделируемой плотности ф(у) = ЩФ ехР(--5уТу) = щщ exp(--5yi) - ехр(--5 2) так как YTY = у\ + ... + у%- Таким образом, плотность s-мерного нормального распределеіїия разбивается на произведение плотностей независимых одномерных нормальных величин. Для моделирования одномерного нормального распределения применялся полярный метод. То есть моделируются одновременно две независимые нормальные величины 1, 2 по следующему алгоритму: Алгоритм 1.1. Формальное моделирование нормального распределения при помощи чисел Холтона. (a) Получаем два независимых числа Холтона pr(ni) = ai,pr(ri2) аг, где г = r(s),n = n(s) - параметры, зависящие от размерности, и строим вспомогательную величину d = (2ai — І)2 + (2ci2 — І)2 (b) Если выполняется d 1, то fi = (2ai - l)y/-2logd/d,2 = (2ai - 1)y/-2logd/d. (с) Иначе, (если d 1), то изменяем параметры чисел Холтона щ щ + 1,П2 = П2 + 1 и возвращаемся к пункту 1. Промоделируем двумерное нормальное распределение с нулевым средним и заданной ковариационной матрицей ( 77.5 3.853103 \ у 3.853103 0.411343 у подставляя в качестве псевдослучайных чисел числа Холтона. Каждый нормальный вектор (1, 2) моделируется с помощью фиксированных для этого вектора чисел щ, пі и набора простых чисел, начиная сЗ. Считаем, что щедт = 200, nend = 30199, при этом щЄдіп nend + 1 = 30000. Получаем выборку из 30000 векторов. Оказалось, что для построенной выборки нет оснований отвергать гипотезу о нормальности на основании х2 теста с доверительным уровнем р = 0.99967 для первой компоненты вектора кр — 0.95940 для второй. После подсчета выборочной ковариационной матрицы, получаем (ИЗ) - /77.274138 3.852429 С = У 3.853103 0.411868 Сравним (1.12) и (1.13). Видно, что моменты совпадают с большой точностью (до второго знака после запятой). Аналогичный результат можно получить для четвертого выборочного момента. Посмотрим сравнительную таблицу значений четвертого выборочного момента для двух выборок нормально распределенных векторов: первая выборка получена моделированием с помощью чисел Холтона, вторая с использованием превдослучайных чисел. Символ Е(...) означает, что рассматривается не истинное значение момента, его выборочная оценка. На самом деле, сохранение четных моментов произойдет и в том случае, если при моделировании выборки случайных векторов вместо чисел Холтона подставить детерминированные значения из промежутка [0,1]. Например, покроем единичный квадрат сеткой, которая разбивает каждую из граней на п частей. Построим выборку из п2 величин -, используя детерминированные значения из промежутка [0,1] в алгоритме 1.1. Другую выборку из п2 величин ? построим с помощью псевдослучайных чисел. Выберем п = 1000 и приведем сравнительные таблицы значений четырех первых центральных моментов для обеих выборок (таб. 1.3-1.6).
Эксперименты по сравнению методов МК, модифицированного МК, КМК и модифицированного КМК
Потребность в решении систем линейных алгебраических уравнений (СЛАУ) часто возникает как в теоретических так и в инженерных прикладных задачах, причем решение может быть ценно само по себе, либо в качестве промежуточного шага для решения более глобальной задачи. В практических задачах размер СЛАУ может быть очень велик, поэтому огромное значение играют как точность метода решения СЛАУ, так и его трудоемкость. Рассмотрим СЛАУ вида ВХ = С, где В = \\bij\\fj i Є Rnxri заданная матрица, а С = (/і,..., fn)T Є Rn заданный вектор. Известно (например см. [1], [36],[37]), что если решение X существует, то оно может быть найдено с помощью: прямых методов, таких как метод исключения Гаусса, метод разложения Холесского, сложность которых есть 0(п3); детерминистических итерационных методов, такими как метод Якоби, Гаусса и Зейделя и различными релаксационными методами, которые преобразуют систему к виду и применяют итерации вида до достижения необходимой точности; при этом сложность есть 0(кп2), где к число шагов; многосеточных методов, которые используют последовательность вложенных сеток, первоначально ищут решение на грубой сетки, затем применяют операторы проектирования и интерполирования для уточнения решения СЛАУ на более мелкой сетке; многосеточные методы обладают вычислительной сложностью 0(п); методов Монте-Карло(МК), которые также требуют преобразования системы к виду и используют случайные траектории для оценки суммы вида Методы МК эффективны для оценки отдельных функционалов от решения СЛАУ, в частности скалярного произведения (Н,Х): где Н -произвольный заданный вектор, в частности они применяются для нахождения отдельной компоненты решения СЛАУ жг- = (ЄІ,Х), где ег- і-й единичный орт. Методы МК требуют 0(п) времени для нахождения п компонент решения на произвольном шаге. Учитывая, что скорость сходимости методов МК есть 0(JV-1/2), где JV число случайных траекторий, для достижения высокой точности необходимо очень большое число шагов МК. Однако, в практических задачах зачастую высокая точность либо не требуется, либо неточное приближение, полученное МК, может быть использовано как начальное приближение для последующих вычислений. При этом метод МК обладает рядом преимуществ. Во-первых, устойчивостью к росту размерности, во-вторых, практически неограниченным параллелизмом ([1], [5], [33]). В работах Ермакова СМ. и Данилова Д.Л. ([2], [3]) показано, что при большой размерности СЛАУ метод МК может быть предпочтительнее детерминистических итерационных методов. Методы МК для решения СЛАУ подразумевают использование цепей Маркова с произвольным числом состояний. Рассмотрим СЛАУ вида где Dn есть n-мерный гиперкуб, принадлежащий n-мерному евклидову пространству, A = {aij}fj-i матрица размерности п х п, F — {/І}"=І- Пусть сходится метод последовательных приближений для решения СЛАУ (2.1), то есть решение X может быть представлено в виде сходящегося ряда Неймана под сходимостью понимается абсолютная сходимость, то есть jAfcFj оо, где \А\ = {ЫЦ=Ъ \F\ = {Ш=1. Рассматривается задача о вычислении методом МК скалярного произведения векторов J = (Х,Н), X искомое решение СЛАУ (2.1), Н Є Dn произвольный заданный вектор. Функционал J можно также выразить через решение X сопряженной к (2.1) системы линейных уравнений Легко показать, что {Х,Н) = (X ,F), то есть исходная задача является вполне симметричной относительно исходной СЛАУ (2.1) и сопряженной к ней (2.3). С системами уравнений (2.1) и (2.3) можно связать однородную цепь Маркова, заданную начальным дискретным распределением р и матрицей переходных вероятностей V = {pij}fj i- Относительно однородной цепи {р, V} делаются естественные предположения, называемые, также, условиями согласования цепи Маркова 1. для таких г, что fa ф О, начальное распределение р\ О 2. для таких (ij), что ( ф 0, переходная вероятность тоже pij ф 0. Обычно предполагают, что величину gi трактуют как вероятность поглощения при переходе из состояния і в состояние j. Так как ді не зависит от і, можно считать, что если известно состояние і, то с вероятностью gi траектория цепи Маркова обрывается, а с вероятностью 1 — gi траектория продолжается и определяется следующее состояние. Выборочная траектория цепи Маркова строится следующим образом. Сначала моделируется случайная величина с дискретной плотностью распределения р и находится состояние го- Затем с вероятностью #г-0 происходит поглощение. Если поглощение не произошло, то моделируется случайная величина с дискретной плотностью
Применение модифицированного КМК при ослаблении условия мажорантной сходимости
Общая идея описываемой модификации метода МК обсуждалась в работах С.М.Ермакова и В.Вагнера ([4]). Автору принадлежит конкретная реализация модификации МК для решения СЛАУ, оптимизация параметров цепи Маркова и утвержденяи об устойчивости описанной модификации.
Как указано выше метод МК для решения СЛАУ X = АХ + F может применяться только при выполнении условия мажорантной сходимости Ai(A) 1, где Ai максимальное собственное число. Это условие может быть очень ограничительным. В частности, разностный аналог краевой задачи для уравнения Пуассона AU = F в двумерном случае (г, &)-номер узла сетки, h - шаг сетки, удовлетворяет условию мажорантной сходимости в случае метода простых итераций, но не удовлетворяет ему в случае введения параметра релаксации, обеспечивающего оптимальную сходимость метода верхней релаксации.
Пусть п достаточно велико, Аі(Л) 1, но Аі(А) 1. Тогда для оценки решения СЛАУ X = АХ + F по-прежнему можно применить метод МК, но с учетом предлагаемой модификации. Применение метода МК предполагает оценку суммы ряда Неймана X = оо ]Г) AlF как интеграла по считающей мере, см., например, [[9]]. г=0 Идея модификации заключается в том, что вместо того, чтобы оценивать сумму ряда Неймана целиком, можно выделить и оценить конечное число слагаемых этой суммы. Зафиксируем достаточно большое число к Є N вычисляемых слагаемых ряда Неймана, тогда где Єк - систематическая погрешность, такая, что при к ь- со, если в качестве нормы матрицы взять спектральный радиус. Пусть имеется начальное приближение решения СЛАУ Хо, тогда исходная система X — АХ + F эквивалентна X — Хо — А{Х — Хо) + (AXQ — Хо + F). Используя метод МК, можно последовательно вычислить оценки вида AY: Z ч- АХ0 + F; на первом шаге вычисляется оценка 1 = AY, где Y = АХо — Хо + F, Z+-Z + &; на втором шаге оценивается 2 = АУ, где У = f і - оценка, полученная на первом шаге, Z 4— Z + &; на шаге к оценивается = АУ, где У = &_і - оценка, полученная на шаге k-1, Z - Z -\-(к В качестве оценки X берется Z. Компоненты всех оценок вида ( = AY) вычисляются не поочередно, а так, что номер оцениваемой компоненты выбирается в соответствии с некоторым распределением Р = рр j = 1,..., п. Подробнее о виде используемых оценок и о выборе оптимального распределения говорится ниже. Вопрос о сходимости модифицированного метода МК рассматривался в работах Ермакова СМ. и Вагнера В. Более подробно устойчивость ММК, примененного совместно с методом верхней релаксации рассматривается в главе 4. В целом для описанной модификации МК можно сделать следующие выводы: 1. Условие Лі(Л) 1 может быть заменено условием Лі(А) 1 за счет асинхронизации алгоритма на определенном этапе, это условие является необходимым условием сходимости метода последовательных приближений для произвольного свободного члена. 2. Свойства параллелизма алгоритма сохраняются за счет введения "круп-нозернистости", так как для оценки на каждом шаге справедливо ESTO = Хт. Это означает, что преимущества метода МК сохраняются при решении гораздо более широкого круга задач, чем предполагалось раньше, в частности в работе ([8]), где обсуждается вопрос о решении сеточного аналога гиперболических уравнений. Таким образом, модифицированный метод МК может быть значимым для решения широкого класса задач математического моделирования. 3. Модифицированный метод МК имеет интересную особенность. При его использовании кратности вычисляемых "интегралов" равны единице, что позволяет использовать алгоритм вместе с квазислучайными числами, избавляясь от неприятного множителя lns N в числителе оценки погрешности.
