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

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

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

ТИМОФЕЕВ Константин Алексеевич

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

05.13.18 - математическое моделирование, численные методы и комплексы программ

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

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

003461055

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

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

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

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

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

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

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

ДЕМЬЯНОВИЧ Юрий Каземирович (Санкт-Петербургский государственный университет)

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

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

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

"Ж" cptjpa.t^_ 2009 г. в {б10

212.232.51 по защите докторских и канди

часов на за-

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

седании совета Д 212.232.51 по защите докторских и кандидатских диссертаций при Санкт-Петербургском государственном университете по адресу: 198504, Санкт-Петербург, Петродворец, Уциверситетский пр., 28, математико-механический факультет, ауд. .

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

Автореферат разослан "Ж" ХН.&0(р1 200! г.

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

диссертационного совета Даугавет И. К.

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

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

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

Актуальность тематики обусловлена, в частности, тем, что, как показано в работе [1], с ростом размерности задачи использование методов Монте-Карло для решения линейных уравнений в ряде случаев становится более предпочтительным, чем использование детерминистических методов.

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

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

В качестве другого подхода можно указать рандомизацию детерминистических методов. Примером применения такого подхода является рандомизированный метод Ньютона, предлагаемый в [4].

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

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

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

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

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

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

Теоретически и практически показывается возможность использования как мелкозернистого, так и крупнозернистого параллелизма при построении оценок.

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

Целью диссертационной работы является разработка новых мето-

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

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

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

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

1. Построенный метод Монте-Карло решения алгебраических уравнений с квадратичной нелинейностью является новым.

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

3. Разработанные алгоритмы применения предлагаемых методов на параллельных вычислительных системах являются новыми.

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

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

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

Апробация работы. Основные результаты диссертации докладывались и обсуждались на семинарах кафедры статистического моделирования математико-механического факультета СПбГУ, на всероссийской конференции по вычислительной математике КВМ-2007 (18-20 июня 2007 г. Академгородок, Новосибирск, Россия). Тезисы доклада опубликованы в программе конференции МС<ЗМС'08 (6-11 июля 2008, Монреаль, Канада).

Работа над диссертацией была поддержана грантом РФФИ (N01-05-00865а).

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

Структура и объем работы. Диссертация объемом 157 страниц состоит из четырех глав, разбитых на разделы и параграфы, списка иллюстраций, списка таблиц и списка литературы. Содержит 14 таблиц, 12 рисунков и список цитируемой литературы.

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

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

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

- метод сжимающих отображений;

- метод Ньютона и квази метод Ньютона;

- метод последовательных линеаризаций;

- рандомизированный метод Ньютона для решения нелинейных алгебраических уравнений;

- методы Эйлера и Рунге-Кутта.

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

Что касается метода Эйлера, показывается, что от его рандомизированного аналога следует ожидать сокращения времени расчетов при увеличении погрешности.

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

Рассматриваются следующие методы:

- рандоминизровапный метод Ньютона;

- рандомизированный метод последовательных линеаризаций;

- рандомизированный метод искусственного хаоса (см. далее в этом разделе).

Метод искусственного хаоса был предложен в работе [8] как аналог методов, основанных на распространении хаоса, используемых в газовой динамике.

Рассмотрим систему уравнений относительно х € Иг

г г

Ж» = /г+ ацхз + ^ г = 1,..., г. (1)

3 = 1 ],к=1

Пусть х = (жь ... ,хг) € К1" — решение этой системы уравнений. Домножая систему (1) слева и справа на XI, I = 1,..., г, получаем систему

г г

Уи = /¿£/ + ^ ачу51 + Х{ ЬцкУзк, Уи - щхи г, I = 1,... ,г.

}=1 j,k=l

В ряде случаев оказывается, что решению у = Цу^ ЦГ^! построенной системы соответствует решение х исходной системы. Можно строить следующий итерационный процесс нахождения решения:

1) выбирается начальное приближение х°, полагается п = 0;

2) решается линейное относительно уп+1 уравнение

У»+1 = Л*Г + ХХ1 +®Г £ (2)

з=1 з'Л=1

3) полагается ,тп+1 = ф{уп+г), где ф = ..., : К2г —> Мг — некоторое отображение;

4) значение п увеличивается на единицу н повторяются шаги 2-3 до срабатывания некоторого критерия остановки.

В диссертации предлагается несколько способов расчета хп+1 при известном значении уп+1. Например в случае, когда относительно решения х

г

известно, что ^ > 0, можно брать

¿=1

1 Т Г ~ 1/2 & (У) "" 2 + у^ ( ' ~ 1> ■ • • >п (3)

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

Теорема 1. (О сходимости метода искусственного хаоса) Пусть выполняются следующие условия:

1) для у € Е2г построим £ € Е2г следующим образом: — уц —

— 1, ...,г. Пусть существует такая константа со > 0, что для всех ||е|| < с0 отображение ф : К2г —» Кг (см. шаг 3 алгоритма) представимо в виде ф(у) = х + Ье + 0(||е|]2), где Ь : М2г —* Ег — некоторое линейное отображение;

2) существует оператор (I — Кх — К2) ~1, где операторы Кг : Е2г —» К2г, Кч : К2г —> К2г определяются выражениями

г г

(1<1У)ц = ~ £ ЬгкЩУки Ь3 = 1,

к-1 к,1=1

3) ги := \\(1-Кг -К2) 1||||Ь||||(/—Кз)ж|| < 1, где Ь — лилейный оператор из условия 1, оператор : Кг —> Мг определяется выражением

г

[К3х){ = (ЧкХк, г,3 = 1,... ,г. к=1

Тогда для любого 5, 0 < 6 < 1 — и), существует такое малое число р > 0, что для всех векторов х° е К.г таких, что ||х — ж°[| < р будет выполняться \\хп — ж|| —► О при п —> оо. При этом скорость сходимости будет линейной: \\хп — < —х\\.

При этом доказывается, что, в частности, преобразование (3) удовлетворяет первому условию теоремы.

Важно отметить, что размерность линейного уравнения (2), решаемого на втором шаге алгоритма, равна г2, что является препятствием на пути его решения детерминистическими методами. Тем не менее при решении этого линейного уравнения методами Монте-Карло и при использовании на третьем шаге преобразования (3) метод становится сравнимым по трудоемкости и требованиям к памяти с рандомизированным методом Ньютона, так как достаточно оценивать значения г+1 линейного функционала от решения.

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

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

В главе осуществляется экспериментальное сравнение указанных в на-

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

ТЕ 7

1 - рандомизировзнньи метод Ньютона;

2- прииейший мет од последовательных линеаризации;

3' рэндоиизироеэнньй метод ИХ;

4 ■ норма невязки начэльн ora приближения.

Количвство траэкторнА

Рис 1. Зависимость нормы невязки от числа траекторий на одном из

примеров

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

В третьей главе рассматриваются методы Монте-Карло решения эволюционных уравнений, основанные на методе Эйлера.

Одним из распространенных методов решения эволюционных уравнений является сведение их к системам обыкновенных дифференциальных уравнений вида

dv(t) dt

= f(t,v), v(0)=v°, t> 0,

(4)

где v{t) = {vi{t),... ,vm{t))T, f{t,v) = (fi(t,v),...,fm{t,v))T. Отображение / может быть конечно-разностной аппроксимацией дифференциального оператора в частных производных. Часто в практически интересных случаях число m превосходит 10s.

Если для решения системы (4) используются методы типа Эйлера или Рунге-Кутта, то каждый шаг по времени требует одного или нескольких

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

Рассмотрим систему дифференциальных уравнений (4). Для этого уравнения метод Эйлера заключается в построении оценок vn величин v(Atn), п £ N вида

vn+1 =vn + Atf(Atn,vn). (5)

Предлагается строить случайные оценки тР+1 в виде

=vn + Atf{Atn, v"), п > 0, (6)

где v:> = и0, f(t,v) — семейство случайных векторов, удовлетворяющих условию E/(i, v) да f(t, v) (смысл приближенного равенства уточняется в тексте диссертации).

Если удается выбрать семейство случайных векторов / так, что время моделирования правой части выражения (б) меньше времени расчета правой части (5), то достигается выигрыш по скорости расчетов. При этом, во-первых, случайные оценки ir оказываются смещенными относительно оценок vn и относительно решения v(Atn) исходного уравнения, и, во-вторых, возникает необходимость исследования поведения ковариационных матриц оценок.

Процесс построения оценок легко распараллеливается. Если для расчетов использовать L вычислителей с независимыми датчиками случайных чисел, то можно в L раз уменьшить дисперсии компонент оценки путем усреднения L независимых оценок г>", построенных для фиксированного п на этих вычислителях (крупнозернистый параллелизм). При этом можно строить доверительные интервалы для оценок. Можно проводить также усреднение оценок на каждом шаге п (мелкозернистый параллелизм). Очевидно, что при расчетах обоими способами можно добиться использования до 100% вычислительных ресурсов параллельной вычислительной системы.

Такой подход исследовался ранее в работах Некруткина В.В., Голян-диной Н.Э., Тура Н.И. и др (см., например, [5]). В этих работах приводятся, в том числе, достаточные условия сходимости оценок к точному решению и исследуется скорость сходимости для некоторых случаев. Исследования проводились для так называемых (п,к)-частичных процессов,

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

iW(ti„) = «Me<>i (7)

где вектор p(t,v) = (pi(t,v),... ,pm(t,v)) при фиксированных t € [О,Г] и v G Rm задает распределение на {1 ,...,m}, а — случайная величина с распределением p(t,v), ei — г-й единичный орт размерности то, функция дг(у) : [0, оо) —» [0,оо) при фиксированном z> 0 определяется выражением

( 1 , У<

gz(y)={ 2(y-z)3-3(y-z)2 + l ,z<y<z+1; [О ,y>2.

Доказана следующая теорема:

Теорема 2. Допустим, что на отрезке t € [О, Г] у уравнения (4J существует единственное решение г)(t). Обозначим с\ := ||f(i)||c[o,T]- Пусть при этом f(t, v) при t 6 [О, Т\ дважды непрерывно дифференцируема по второй переменной для всех v таких, что ||г>|| < с\.

Пусть существует такая константа сг, что для всех г — 1 ,...,тп выполняется Pi(t,v) > c2<?2ci (11^Ц)|/г(^jг')I или

т

Pi(t,v) > c2g2ci{\\v\\)\fiM\/ £ №,v)\.

Рассмотрим последовательность случайных векторов vn, п>0, строящихся по формулам

vn+l = i/1 + T£(2ci)(Ain,i/"),n > 0, г/° =гА

Тогда равномерно not G [О ,Т] при At стремящихся к нулю выполняется:

1) ЦЕг/LVAtJ = 0(At);

2) || cov = О (At) (здесь cov обозначает ковариационную матрицу случайного вектора);

3) ЦБ- v(t))(vWAti - t;(i))T|| = O(At).

То есть скорость сходимости предлагаемого варианта рандомизированного метода Эйлера в этом случае имеет тот же порядок, что и у детерминистического метода Эйлера.

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

Также приводится класс оценок, применимых для решения уравнения (4), функция /(£,«) в правой части которого квадратично нелинейна по V. Для этого класса оценок доказывается утверждение, аналогичное теореме 2.

Глава 4. Метод, результаты теоретического исследования которого приведены в главе 3, лег в основу программного комплекса, предназначенного для численного исследования его свойств на примере уравнения Навье-Стокса, описывающего движение вязкой несжимаемой жидкости в двумерном и трехмерном случаях. Рассматривался как общий случай уравнения, так и вариант с искусственной сжимаемостью. Решение уравнения осуществлялась в естественных координатах (см. [2] и [3]).

Комплекс состоит из нескольких компонент: средств расчета и средств визуализации. Их подробное описание приводится в тексте диссертации.

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

Если в случае дискретного по у уравнения (4) под сложностью расчетов понимать число компонент вектор-функции /, которое нужно рассчитать для нахождения оценки решения на очередном шаге но времени, то показано, что с уменьшением шага по времени можно, при некоторой потере точности, получить существенный выигрыш по трудоемкости (вплоть до т раз, где тп — размерность решаемого уравнения).

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

Наконец, приводятся примеры успешного расчета решения уравнения Навье-Стокса на достаточно мелкой сетке (размерность решаемой системы уравнений оказывается равной 3 • 106).

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

кавернах с одной подвижной стенкой.

Рис 3. Линии тока, проходящие сквозь плоскости 2 = 0.45 и г = 0.55

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

3. ЗАКЛЮЧЕНИЕ

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

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

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

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

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

СПИСОК ЛИТЕРАТУРЫ

[1] Ermakov S.M. Danilov D.L., Halton J. Asymptotic complexity of Monte-Carlo methods for solving linear systems // Journ. of statistical planning and inference. 85 (2000)N 1-2, 5-18. ed.V.Melas. Special issue. St.Petersburg State Univ. St. Petersburg, Russia.

[2] P. Темам. Уравнения Навье-Стокса. -M.: Мир. 1981.

[3] Андерсон Д., Таннехилл Дж., Плетчер Р. Вычисл. гидродинамика и теплообмен. В 2-х т. - М., Мир, 1990

[4] Halton Jh. Sequential Monte Carlo Techniques for Solving Non-Linear Systems // Monte Carlo Methods and Appl. 2006. Vol. 12. P. 113-141.

[5] В.В.Некруткин, Н.И Тур. К обоснованию схемы прямого статистического моделирования течений разреженных газов // Журнал вычислит. математики и матем. физики, 29, 9, стр. 1380-1391. 1989.

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

[6] Тимофеев К.А. Об одном классе методов Монте-Карло для решения уравнений с квадратичной нелинейностью.// Вестн. С.-Петербург, унта, СПб., Изд. СПбГУ, серия 10, вып. 3, с. 105 - 110, 2008.

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

[8] Калошин И.В., Ермаков С.М., Тимофеев К.А. Метод искусственного хаоса для решения методом Монте-Карло уравнений с квадратичной нелинейностью // Математические модели. Теория и приложение. Сб. науч. статей, под ред. Чиркова М.К., С-Пб., вып. 7, с. 3-20, 2006.

Отпечатано копировально-множительным участком отдела обслуживания учебного процесса физического факультета СПбГУ. Приказ № 571/1 от 14.05.03. Подписано в печать 25.12.08 с оригинал-макета заказчика. Ф-т 30x42/4, Усл. иеч. л. 1. Тираж 100 экз., Заказ № 911/с. 198504, СПб, Ст. Петергоф, ул. Ульяновская, д. 3, тел. 929-43-00.

Введение

1 Общие сведения

1.1 Введение.

1.2 Методы решения нелинейных алгебраических уравнений

1.2.1 Метод сжимающих отображений.

1.2.2 Метод Ньютона-Канторовича.

1.2.3 Метод последовательной линеаризации решения уравнений с квадратичной нелинейностью.

1.2.4 Рандомизированный метод Ньютона.

1.3 Класс методов Рунге-Кутта решения обыкновенных дифференциальных уравнений.

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

2 Решение алгебраических уравнений с квадратичной нелинейностью методами Монте-Карло

2.1 Введение.

2.2 Рандомизация итерационных методов для решения нелинейных уравнений.

2.3 Метод искусственного хаоса.

2.3.1 Общий случай

2.3.2 Случай алгебраических уравнений.

2.4 Рандомизированный метод искусственного хаоса.

2.4.1 Сходимость рандомизированного метода искусственного хаоса в смысле вторых моментов.

2.4.2 Замечания к практическому применению рандомизированного метода искусственного хаоса.

2.4.3 Алгоритм моделирования.

2.4.4 Требования к вычислительным ресурсам. Сравнение с методом Ньютона.

2.5 Вычислительные эксперименты.

2.5.1 Условия сходимости.

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

2.6 Выводы.

3 Решение обыкновенных дифференциальных уравнений методом Монте-Карло

3.1 Введение.

3.2 Рандомизированный метод Эйлера.

3.2.1 Примеры

3.3 Параллелизм.

3.4 Локально оптимальные параметры.

4 Вычислительные эксперименты

4.1 Введение.

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

4.2.1 Уравнение Навье-Стокса.

4.2.2 Уравнение Навье-Стокса с искусственной сжимаемостью

4.2.3 Построение разностной аппроксимации.

4.3 Алгоритмы решения.

4.4 Описание программного комплекса.

4.5 Результаты вычислительных расчетов.

4.5.1 Сходимость.

4.5.2 Мелкозернистый параллелизм

4.5.3 Крупнозернистый параллелизм.

4.5.4 Существенная и квазисущественная выборка.

4.5.5 Расчет оценки решения уравнения Навье-Стокса в переменных скорость-давление.

4.5.6 Сравнение трудоемкости с методом Эйлера.

4.5.7 Расчет задач большой размерности.

4.6 Использование квазислучайных чисел.

4.7 Выводы.

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

Сложность решаемых задач диктует требования к вычислительным ресурсам. В настоящее время все большую популярность приобретают многопроцессорные вычислительные системы. Для эффективного использования ресурсов таких систем алгоритмы расчета должны обладать свойством параллелизма. В то время как многие детерминистические методы не позволяют эффективно использовать вычислительные ресурсы многопроцессорных систем или требуют для этого применение сложных процедур (см, например, [11]), алгоритмы Монте-Карло, как правило, легко поддаются модификации, позволяющей распараллелить процесс расчетов (см. например [6], [12] и [13]). Исследования, проведенные в рамках диссертации, направлены па разработку методов Монте-Карло, использующих порядка 100% имеющихся вычислительных ресурсов. Такие методы можно применять при больших размерностях решаемых уравнений в случаях, когда высокая точность вычислений не требуется.

Актуальность тематики обусловлена, в частности, тем, что, как показано в работе [2], с ростом размерности задачи использование методов Монте-Карло для решения линейных уравнений в ряде случаев становится более предпочтительным, чем использование детерминистических методов.

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

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

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

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

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

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

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

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

Теоретически и практически показывается возможность использования как мелкозернистого, так и крупнозернистого параллелизма при построении оценок.

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

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

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

- метод сжимающих отображений;

- метод Ньютона и квази метод Ньютона;

- метод последовательных линеаризаций;

- рандомизированный метод Ньютона для решения нелинейных алгебраических уравнений;

- методы Эйлера и Рунге-Кутта.

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

Что касается метода Эйлера, показывается, что от его рандомизированного аналога следует ожидать сокращения времени расчетов при увеличении погрешности.

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

Рассматриваются следующие методы:

- рандоминизрованный метод Ньютона;

- рандомизированный метод последовательных линеаризаций;

- рандомизированный метод искусственного хаоса (см. далее).

Метод искусственного хаоса был предложен в работе [31] как аналог методов, основанных на распространении хаоса, используемых в газовой динамике.

Рассмотрим систему уравнений относительно ж Е1г г г

Х{ — -4" ^ ^ ^ ^ Ьц^х^хк) % ~ 1,., т. (1)

Пусть х = (^1,. ,хг) Е 1г - решение этой системы уравнений. Домножая систему (1) слева и справа на ж/, I = 1,., г, получаем систему г г

Уи = + Е сц^уц + XI Е ЪцкУук, Уи = Х1Х1, г, I = 1,., г.

1 3,к=1

В ряде случаев оказывается, что решению у = Ц^Ц^-х построенной системы соответствует решение х исходной системы. Можно строить следующий итерационный процесс нахождения решения:

1) выбирается начальное приближение ж0, полагается п = 0;

2) решается линейное относительно уп+1 уравнение г г уГ1 = №+Е +х1 Е (2)

3=1 ¿,*=1

3) полагается хп+1 = ф(уп+1), где ф = {ф\,., фг) : М2г Мг — некоторое отображение;

4) значение п увеличивается на единицу и повторяются шаги 2-3 до срабатывания некоторого критерия остановки.

В диссертации предлагается несколько способов расчета хп+1 при известном значении уп+1. Например в случае, когда относительно решения х известно, г что > 0, можно брать г=1

1 ^ ^ —1/2 ^ (у) = 2 + ^ ' * = 1, • ■ ■. г- (3) 1

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

Теорема 1. (О сходимости метода искусственного хаоса) Пусть выполняются следующие условия:

1) для у € М2г построим е 6 К г следующим образом: е^ — у^ — XiXj; 1 ,.,г. Пусть существует такая константа со > 0, что для всех Це|| < со отображение ф : М2г —> Жг (см. шаг 3 алгоритма) представимо в виде ф(у) = х + Ье + 0(||е:||2); где Ь : К2г —у Мг — некоторое линейное отображение;

2) существует оператор (/ — К\ — где операторы К\ : К2г —> Ж2г, К2 М2г —> М2г определяются выражениями г г

К1У)ч — <ЧкУк], у)ч — ]Г) Ьы№зУки Ь 3 — 1, • • •, Т\ к=1 к,1=1

3) V) := ||(1 - Кг- К2)~1\\\\Ь\\\\(1 - Кз)х\\ < 1, где Ь - линейный оператор из условия 1, оператор : Кг —> Кг определяется выражением г

К3х)г = £ а1кхк, г, 7 = 1,., г.

Тогда для любого 5, 0 < 5 < 1 — т, существует такое малое число р > О, что для всех векторов х° Е КГ таких, что ||ж — < р будет выполняться ||хп — ж|| —> 0 при п —¥ оо. При этом скорость сходимости будет линейной: ||жп-ж|| < +^Цж"-1 - ж||.

При этом доказывается, что, в частности, преобразование (3) удовлетворяет первому условию теоремы.

Важно отметить, что размерность линейного уравнения (2), решаемого на втором шаге алгоритма, равна г2, что является препятствием на пути его решения детерминистическими методами. Тем не менее при решении этого линейного уравнения методами Монте-Карло и при использовании на третьем шаге преобразования (3) метод становится сравнимым по трудоемкости и требованиям к памяти с рандомизированным методом Ньютона, так как достаточно оценивать значения г + 1 линейного функционала от решения.

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

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

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

В третьей главе рассматриваются методы Монте-Карло решения эволюционных уравнений, основанные на методе Эйлера.

Одним из распространенных методов решения эволюционных уравнений является сведение их к системам обыкновенных дифференциальных уравнений вида где = (г>1(£),., г>т(£))т, /(¿,17) = (/х(г, г>), ., /т(£, у))т. Отображение / может быть конечно-разностной аппроксимацией дифференциального оператора в частных производных. Часто в практически интересных случаях число т превосходит 106.

Если для решения системы (4) используются методы типа Эйлера или Рунге-Кутта, то каждый шаг по времени требует одного или нескольких вычислений значения вектора /(£, у). В случае функций / сложного вида и в случаях большой размерности т решаемой системы, это предъявляет большие требования к вычислительным системам и требует разработки специальных приемов распараллеливания. Прием, описанный ниже (рандомизация), позволяет в ряде случаев уменьшить вычислительную работу и указать простую процедуру распараллеливания алгоритма.

Рассмотрим систему дифференциальных уравнений (4). Для этого уравнения метод Эйлера заключается в построении оценок Vй величин г>(Д£п), ггЕМ вида „(0)=«°, (>0,

4) = „» +ДгДДйг,,,»).

5)

Предлагается строить случайные оценки vп+ в виде + Atf(Atn, Vй), п> О, (6) где = i;0, f(t,v) — семейство случайных векторов, удовлетворяющих уело -^ч вию E/(t, v) « /(t, v) (смысл приближенного равенства уточняется в тексте диссертации).

Если удается выбрать семейство случайных векторов / так, что время моделирования правой части выражения (6) меньше времени расчета правой части (5), то достигается выигрыш по скорости расчетов. При этом, во-первых, случайные оценки гт71 оказываются смещенными относительно оценок vn и относительно решения v(Atn) исходного уравнения, и, во-вторых, возникает необходимость исследования поведения ковариационных матриц оценок.

Процесс построения оценок легко распараллеливается. Если для расчетов использовать L вычислителей с независимыми датчиками случайных чисел, то можно в L раз уменьшить дисперсии компонент оценки путем усреднения L независимых оценок Vй, построенных для фиксированного п на этих вычислителях (крупнозернистый параллелизм). При этом можно строить доверительные интервалы для оценок. Можно проводить также усреднение оценок на каждом шаге п (мелкозернистый параллелизм). Очевидно, что при расчетах обоими способами можно добиться использования до 100% вычислительных ресурсов параллельной вычислительной системы.

Рассмотрим семейство случайных векторов с параметрами z: t и v

7) где вектор p(t,v) = (pi(t,v),. ,pm(t,v)) при фиксированных t G [0, Т] и v £ Mm задает распределение на {1, .,m}, а — случайная величина с распределением p(t,v), е* — г-й единичный орт размерности га, функция

9г(у) = <

9г{у) '■ [0, оо) —[0, оо) при фиксированном г > 0 определяется выражением

1 , У < я;

2 (у -г)*-Ъ{у-г)2 + 1 , * < у < * + 1; О

Доказана следующая теорема:

Теорема 2. Допустим, что на отрезке £ 6 [0,Т] у уравнения (4) существует единственное решение и{р). Обозначим С\ := И^СОЦсрДЧ- Пусть пРи этом г;) при t Е [О, Т] дважды непрерывно дифференцируема по второй переменной для всех V таких, что ||г>|| < с\.

Пусть далее существует такая константа С2, что для всех г — 1,., га выполняется р^, и) > С2д2С1{\\у\\)\^(£1У)\ или

ТП

3=1

Рассмотрим последовательность случайных векторов ип, п > О, строящихся по формулам ь>п+1 = ип + ге(2с1)(ЛЫ, 1Уп), п > О, I/0 = у

Тогда равномерно по £ Е [О, Г] при стремящихся к нулю выполняется:

1) — г>(£)|| — 0(Д£);

2) \\cov гЛ*/л*-1|| = 0(Д£) (здесь сои обозначает ковариационную матрицу случайного вектора);

3) ||Е(|Д'/Д^ - - ^))Г|| = 0(Д*).

То есть скорость сходимости предлагаемого варианта рандомизированного метода Эйлера в этом случае имеет тот же порядок, что и у детерминистического метода Эйлера.

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

Также приводится класс оценок, применимых для решения уравнения (4), функция f(t,v) в правой части которого квадратично нелинейна по v. Для этого класса оценок доказывается утверждение, аналогичное теореме 2.

Метод, результаты теоретического исследования которого приведены в главе 3, лег в основу программного комплекса, предназначенного для численного исследования его свойств на примере уравнения Навье-Стокса, описывающего движение вязкой несжимаемой жидкости в двумерном и трехмерпом случаях. Четвертая глава содержит некоторые полученные результаты расчетов. Рассматривался как общий случай уравнения, так и вариант с искусственной сжимаемостью. Решение уравнения осуществлялась в естественных координатах (см. [18] и [22]).

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

Если в случае дискретного по у уравнения (4) под сложностью расчетов понимать число компонент вектор-функции /, которое нужно рассчитать для нахождения оценки решения на очередном шаге по времени, то показано, что с уменьшением шага по времени можно при некоторой потере точности получить существенный выигрыш по трудоемкости (вплоть до т раз, где m — размерность решаемого уравнения).

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

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

Наконец, в конце главы приводится пример расчета уравнения на сетке из 3 • 106 узлов.

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

4.7 Выводы

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

Заключение

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

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

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

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

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

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

2. Halton Jh. Sequential Monte Carlo Techniques for Solving Non-Linear Systems // Monte Carlo Methods and Appl. 2006. Vol. 12. P. 113-141.

3. Ermakov S.M., Wagner W. Monte Carlo difference schemes for the wave equation // Monte Carlo Methods and Appl. 2002. Vol. 8, N 1. P. 1-29.

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

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

6. L. Devroye. Non-Uniform Random Variate Generation, Springer-Verlag, New York, 1986.

7. В.В.Некруткин, Н.И Тур. К обоснованию схемы прямого статистического моделирования течений разреженных газов // Журнал вычислит, математики и матем. физики, 29, 9, стр. 1380-1391. 1989.

8. V.Nekrutkin, N.Tur. Asymptotic expansions and estimators with small bias for Nanbu proceses // Monte Carlo Methods and Appl., V.3, No 1, 1-35. 1997

9. N.Golyandina. Central Limit Theorem for (n,k)-particle processes solving balance equations. Monte Carlo Methods and Appl., 9, 1, p.1-11. 2003.

10. В.П. Гергель, Р.Г. Стронгип. Основы параллельных вычислений для многопроцессорных вычислительных систем. Издательство Нижегородского госуниверситета. Нижний Новгород. 2003.

11. I. Dimov. A. Karaivanova. Monte-Carlo parallel algorithms.

12. J. Rosenthal. Parallel computing and Monte Carlo algorithms. // Far East Journal of Theoretical Statistics 2000. 4. p. 207-236.

13. Hairer E. Introduction a l'Analyse Numerique, Universite de Geneve, 2001.

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

15. А.А. Амосов, Ю.А. Дубинский, Н.В. Копченова. Вычислительные методы для инженеров. М.-Высшая школа. 1994.

16. Ландау Л.Д., Лифшиц Е.М. Гидродинамика. Москва. "Наука". 1986.

17. Темам Р. Уравнения Навье-Стокса. Теория и численный анализ. -М.: Мир, 1981.

18. Зоткевич А.А., Лаевский Ю.М. Численное моделирование распространения ламинарного пламени на основе двухуровневых явных разностных схем // Вычислительные технологии. 2006. Т. 11. N 6. С. 31-43

19. Роуч П. Вычислительная гидродинамика. М., Мир. 1980. 612 с.

20. Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен. Том 1. -М.: Мир, 1990.

21. Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен. Том 2. -М.: Мир, 1990.

22. J. Halton. Sequential Monte Carlo techniques for the solution of linear systems // Journal of Scientific Computing, vol. 9. issue 2. p. 213 257. 1994.

23. Prandtl L., Tietjens O. Hydro- und Aeromechanic. Berlin, 1929. P. 175-208.

24. Ламб Г. Гидромеханика: Пер. с англ./ Под ред. Н.А. Слезкина. М.-Л: Гостехтеориздат, 1947. С. 99, 304, 839.

25. Reynolds J. An Experimental Investigation of the Circulation which determine the Motion of Water shall be Direct or Sin. and the Law of Resistance in parallel Channels. Phil. Trans. CLXXIV, 935,1883. Papers 11,51.

26. Лаврентьев M.A., Шабат E.B. Проблемы гидродинамики и их математические модели. М.: Наука, 1973. С. 339, 340.

27. Ю.В. Лапин. Статистическая теория турбулентности // Научно технические ведомости 2, 2004. Проблемы турбулентности и вычислительная гидродинамика

28. Василевский Ю.В., Липников К. Error estimates for Hessian-based mesh adaptation algorithms with control of adaptivity //Proceedings of 13-th International Meshing Roundtable, September 19-22, Williamsburg, Virginia, 345-351.

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

30. Тимофеев К.А. Об одном классе методов Монте-Карло для решения уравнений с квадратичной нелинейностью.// Вестн. С.-Петербург, ун-та, СПб., Изд. СПбГУ, серия 10, вып. 3, с. 105 110, 2008.

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