Алгоритмы статистического моделирования решений уравнений эллиптического и параболического типа

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

СИМОНОВ Николай Александрович АЛГОРИТМЫ СТАТИСТИЧЕСКОГО МОДЕЛИРОВАНИЯ РЕШЕНИЙ УРАВНЕНИЙ ЭЛЛИПТИЧЕСКОГО И ПАРАБОЛИЧЕСКОГО ТИПА 01.01.07 – Вычислительная математика

АВТОРЕФЕРАТ

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

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

Работа выполнена в учреждении Российской академии наук, Институте вычислительной математики и математической геофизики СО РАН

Официальные оппоненты: доктор физико-математических наук, профессор Андросенко Пётр Александрович (Обнинский институт атомной энергетики) доктор физико-математических наук, профессор Григорьев Юрий Николаевич (Институт вычислительных технологий СО РАН) доктор физико-математических наук, профессор Рябов Виктор Михайлович (Санкт-Петербургский государственный университет)

Ведущая организация: Новосибирский государственный университет

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

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

Автореферат разослан 2010 года.

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

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

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

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

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

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

Во многих случаях вероятностные модели для описания какого-либо феномена используются в совокупности с различными другими моделями, отличающимися друг от друга масштабами, степенью детализации и, как следствие, применяемым математическим аппаратом. Согласование этих моделей, а также алгоритмов, создаваемых для решения поставленных в рамках этих моделей задач, является естественным требованием, обеспе чивающим их обоснованность и состоятельность. Для природных явлений, которые на макроуровне описываются дифференциальными уравнениями в частных производных эллиптического и параболического типа, условия согласования состоят в том, что макропараметр, удовлетворяющий урав нению, представляется в виде функционала от случайного диффузионного процесса. Построение оценок метода Монте-Карло для такого функциона ла обычно основано на моделировании траекторий некоторой марковской цепи. Начало разработки и применения алгоритмов статистического моде лирования к решению уравнений эллиптического и параболического типа восходит к пятидесятым годам двадцатого века. К этому времени мето ды Монте-Карло уже являлись основным вычислительным инструментом решения задач, связанных с переносом излучения. Основоположниками этих методов были J.v.Neumann и S.Ulam, а в СССР развитие и практи ческое применение алгоритмов статистического моделирования для реше ния уравнения переноса связано с именами Г.И.Марчука, В.Г.Золотухина, С.М.Ермакова, Г.А.Михайлова, И.М.Соболя, А.И.Хисамутдинова, Л.В.Ма йорова, В.В.Учайкина и многих других. Развитие методов Монте-Карло в применении к решению уравнений в частных производных эллиптического и параболического типа восходит к работам W.Wasow, J.Curtiss, M.Muller, G.Brown, A.Haji-Sheikh. Интенсивные исследования в этом направлении были инициированы С.М.Ермаковым, Г.А.Михайловым и продолжают ве стись ими и их учениками: К.К.Сабельфельдом, А.С.Расуловым, А.С.Си пиным, О.Курбанмурадовым, А.А.Кронбергом, Б.С.Елеповым, В.Вагнером и другими, а также Г.Н.Мильштейном, D.Talay, S.Maire, A.Lejay, M.Masca gni, J.Given, I.Dimov и многими другими. Данная диссертация продолжает традиции новосибирской школы методов Монте-Карло. Работа над ней ве лась в рамках исследований, проводившихся группой Стохастических за дач математической физики Института вычислительной математики и ма тематической геофизики Сибирского отделения Российской академии на ук. Начиная с 1985 года, группу возглавляет профессор К.К.Сабельфельд, чьи взгляды и идеи оказали существенное влияние на формирование на учных интересов автора и направление его собственных исследований.

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

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

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

Основные цели и задачи работы.

• Построение и обоснование новых алгоритмов статистического моде лирования решений уравнений эллиптического и параболического ти па.

• Эффективная компьютерная реализация построенных вычислитель ных методов.

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

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

Научная новизна.

Все основные результаты работы являются новыми и заключаются в сле дующем.

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

Апробация работы. Результаты, изложенные в диссертации, регу лярно, начиная с 1979 года, докладывались и обсуждались на семинарах отдела Статистического моделирования в физике Вычислительного центра (Института вычислительной математики и математической геофизики) СО РАН под руководством чл.-корр. РАН Г.А.Михайлова.

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

VII всесоюзная конференция ‘Методы Монте-Карло в вычислитель ной математике и математической физике’, Новосибирск, 1985;

Всесоюзная конференция ‘Актуальные проблемы вычислительной и прикладной мате матики’, Новосибирск, 1987;

Третья республиканская конференция ‘Инте гральные уравнения в прикладном моделировании’, Одесса, 1989;

Актуаль ные проблемы статистического моделирования и его приложения, Ташкент, 1989;

VIII всесоюзная конференция ‘Методы Монте-Карло в вычислитель ной математике и математической физике’, Новосибирск, 1991;

Междуна родная конференция AMCA-95, Новосибирск, 1995;

Математические мо дели и численные методы механики сплошной среды, Новосибирск, 1996;

The 2nd St. Petersburg Workshop on simulation, St. Petersburg, 1996;

GAMM Annual Meeting, Regensburg, Germany, 1997;

First IMACS Seminar on Monte Carlo Methods, Brussels, Belgium, 1997;

15th IMACS World Congress, Berlin, Germany, 1997;

Munchener Stochastik-Tage, Munich, Germany, 1998;

SibIN PRIM-98, Новосибирск, 1998;

The 3rd St. Petersburg Workshop on Simulation, St. Petersburg, 1998;

SibINPRIM-2000, Новосибирск, 2000;

Algorithms and Complexity for Continuous Problems, Schloss Dagstuhl, Wadern, Germany, 2000;

The 4th St. Petersburg Workshop on Simulation, St. Petersburg, 2001;

Международная конференция по вычислительной математике ICCM-2002, Новосибирск, 2002;

The International Conference on Computational Science ICCS-2003, St. Petersburg, Russia, 2003;

4th International Conference on ‘Large Scale Scientic Computations’, Sozopol, Bulgaria, 2003;

IVth IMACS Seminar on Monte Carlo Methods, Berlin, Germany, 2003;

AMS 2004 Spring South eastern Section Meeting, Tallahassee, USA, 2004;

NATO Advanced Research Workshop: Advances in Air Pollution Modelling for Environmental Security, Borovetz, Bulgaria, 2004;

Международная конференция по вычислительной математике ICCM-2004, Новосибирск, 2004;

SIAM Conference on Computa tional Science and Engineering, Orlando, USA, 2005;

The 5th IMACS Seminar on Monte Carlo Methods, Tallahassee, USA, 2005;

The International Conference on Computational Science ICCS-2005, Atlanta, USA, 2005;

5th International Conference on ‘Large Scale Scientic Computations’, Sozopol, Bulgaria, 2005;

17th IMACS World Congress, Paris, France, 2005;

7th International Conference on Monte Carlo and Quasi-Monte Carlo Methods in Scientic Computing, Ulm, Germany, 2006;

6th Conference on Numerical Mathematics and Applications, Borovets, Bulgaria, 2006;

Workshop on Quantitative Computational Biophysics, Tallahassee, USA, 2007;

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

The 6th IMACS Seminar on Monte Carlo Methods, Reading, UK, 2007;

Международная конференция ‘Диффе ренциальные уравнения. Функциональные пространства. Теория прибли жений’ посвященная 100-летию со дня рождения С.Л. Соболева, Новоси бирск, 2008.

Публикации. Результаты диссертации изложены в 41 опубликован ной работе, в том числе в двух монографиях и 12 статьях, напечатанных в журналах, рекомендованных ВАК для опубликования основных научных результатов диссертаций на соискание учёной степени доктора наук. Пе речень публикаций приведён в конце автореферата. Издания [3, 4, 6, 12] входят в список ВАК, а издания [5, 7 – 10, 11, 13, 14] входят в систему цитирования Web of Science.

В совместных работах [1, 2] Н.А.Симонову принадлежит детальная раз работка алгоритмов случайного блуждания для решения первой, второй и третьей краевых задач для уравнения Лапласа и уравнений Ламе. В работе [4] автору диссертации принадлежит разработка алгоритма для вычисле ния итераций сингулярного интегрального оператора. В работах [10, 11, 13, 30 – 38] Н.А.Симонову принадлежит разработка и обоснование алгоритма, проведение численных экспериментов и анализ результатов. В работе [14] расчёты проводились с использованием программы, созданной автором, он принимал участие в анализе результатов и формулировке выводов. В ра боте [19] автору принадлежит алгоритм блуждания в подобластях и его реализация в виде подпрограмм. В работе [27] Н.А.Симонову принадле жит разработка и реализация алгоритма метода Монте-Карло и проведе ние численных экспериментов.

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

Структура и объём работы. Диссертация состоит из введения, трёх глав, разбитых на параграфы, дополнения и заключения. Результаты ис следований изложены на 286 страницах с использованием 26 рисунков и таблиц. Библиографический список состоит из 282 наименований.

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

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

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

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

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

Пусть область, в которой ищется решение, представляется в виде объ единения пересекающихся подобластей Gj, в каждой из которых оценка решения задачи Дирихле в точке x имеет вид [uj ](x) = gj (t )Q(x;

t ;

t1,..., tN ), (1) где ti Gj, i 1 – марковская цепь случайных блужданий, t – точка первого выхода на границу этой подобласти, а gj – значение функции uj на этой границе. Тогда оценка для решения во всей области представляется в виде (K ) K k1 k k k [u](x) = g(t ) Q(t ;

t ;

t1,..., tNk ) (2) k= и обладает всеми свойствами оценок (1). Здесь tk – последовательность точек выхода из подобластей, t1 x, g – функция, заданная на границе.

Теорема 1. 1. Если оценка [uj ](x), определяемая равенством (1), явля ется несмещённой оценкой решения задачи Дирихле в каждой из подоб ластей Gj, рассматриваемых отдельно, то оценка [u](x), построенная в соответствии с формулой (2), будет несмещённой оценкой решения за M дачи Дирихле во всей области G = Gj.

j= 2. Если оценка (1), обладает некоторым смещением:

E [uj ](x) = uj (x) + j (x), (3) где |j (x)| Bj, а Bj = sup |uj (x)|, то оценка (x) решения задачи Ди xGj рихле в области G также смещённая и | E (x) u(x)| B, 1 ( + q) где q = max qi,k, B = sup |u(x)|. Здесь qj,k 1 – константа из леммы i,k xG Шварца в применении к паре подобластей Gj и Gk.

3. Дисперсия оценки (x) конечна.

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

Эффективные методы статистического моделирования строятся на основе рандомизации этой формулы.

Пусть функция u является решением линейного уравнения Пуассона Больцмана (Гельмгольца):

u 2 u = 0, = const 0 (4) в области G R3. Обозначим 1 sinh((a |y x|)),a (y) = 4 |y x| sinh(a) функцию Грина задачи Дирихле для уравнения (4) в шаре B(x, a). То гда значение решения в точке x G удовлетворяет при любом a d(x) формуле среднего (d(x) – расстояние до границы области ) 1 a u(x) = u ds. (5) S(x,a) 4a sinh(a) Рассмотрим смешанную краевую задачу для уравнения (4):

u (y) + (y)u(y) = g(y), y, (y) (6) n где = 1, = 0 на 0 и наоборот = 0, = 1 на 1 = \ 0. Для построе ния алгоритма случайного блуждания выведено и используется следующее интегральное соотношение,a u(x) = 2 u ds n S\{x} u + [2,a ] ds, (7) n B(x,a)\{x} где S – граница области G B(x, a). Здесь x 0. Рандомизация соотно шений (5), (7) позволяет построить марковскую цепь случайного блуж дания {xk, k = 0, 1,..., N } и определить оценку решения краевой за дачи u(x0 ) в виде функционала от траектории этой цепи. Внутри обла сти используется блуждание по сферам, основанное на рандомизации (5) при a = d(x), вплоть до первого попадания точки xk в -окрестность границы. При выходе цепи в окрестность границы с условиями Неймана u (y) = g(y) к оценке решения прибавляется величина, равная второму n интегралу в (7) при x = x, x 0 – ближайшая к xk. При этом модели k k рование марковской цепи продолжается, а её следующая точка выбирается на S в соответствии с равномерным распределением по углу видимости из x. При x 1 используется граничное условие Дирихле и моделиро N k вание цепи обрывается.

Теорема 2. Функция [u](x0 ), построенная на траектории случайного блуждания по сферам, которое используется как внутри области, так и при моделировании отражения от границы, является оценкой для реше ния смешанной краевой задачи (4), (6). Дисперсия этой оценки конечна, а смещение есть величина того же порядка, что и ширина полосы вбли зи границы, то есть совпадает по порядку со смещением оценки решения задачи Дирихле. Трудоёмкость построенной оценки (для выпуклой грани цы) есть O(log() 2 ). Здесь, требуемая точность решения, совпадает с шириной приграничной полосы.

Если 0, то утверждение остаётся верным и для задачи Неймана.

Для невыпуклой границы используется аппроксимация интегрального представления (7), которая позволяет преодолеть алгоритмические трудно сти, вызванные знакопеременностью ядра интегрального оператора. Вво дится касательная плоскость в точке x и следующая точка марковской k цепи выбирается равномерно на полусфере S(x, a), отсекаемой этой плос k костью. Выбирая = O(a/2R)3 при малых a/R, где R – минимальный радиус кривизны поверхности в данной точке, получаем оценку с трудо ёмкостью O(log() 5/2 ).

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

Пусть функция ue удовлетворяет линейному уравнению Пуассона-Больц мана (4) во внешности ограниченной области G, а функция ui удовлетво ряет внутри G уравнению Пуассона i ui = 4. (8) На границе заданы условия непрерывности решения и потока:

ui ue u e = ui, i (y) = e (y). (9) n n Для построения алгоритма блуждания по сферам в G и её внешности G1 используется рандомизация соотношения о среднем (5). Во внешней области при 0 весовой множитель в этой интегральной формуле ин терпретируется как вероятность выживания. В случае = 0 используется модификация алгоритма блуждания по сферам с прямым моделированием ухода диффузионной траектории на бесконечность.

Правая часть уравнения Пуассона учитывается с помощью представ 1 ления ui в виде суммы объёмного потенциала g(x) = (y)dy, i |x y| G и регулярного слагаемого. Учёт условий на границе раздела приводит к интегральной формуле (x ):

e 1 a u(x) = ue ds e + i S(x,a) G1 2a2 sinh(a) i 1 a + u ds 2 sinh(a) i e + i S(x,a) G 2a (e i ) 0,a 2 Q,a u ds (10) e + i B(x,a)\{x} n i [22,a ]ui dy + e + i B(x,a) G + [2,a ] dy, e + i B(x,a) G sinh((a r)) + r cosh((a r) 1, r = |x y|.

где Q,a (r) = sinh(a) Для построения марковской цепи используется приближённая рандоми зация этой формулы, основанная на разделении шара B(x, a) касательной e плоскостью в точке x = x. С вероятностью pe = следующая точка k i + e марковской цепи xk+1 выбирается равномерно на той полусфере S + (x, a), k в которую направлен вектор нормали, а с дополнительной вероятностью pi – в замкнутом полушаре, B (x, a). Значение объёмного потенциала с k плотностью оценивается по одной случайной точке, которая выбирается с плотностью, согласованной с.

При 0 коэффициент, стоящий под интегралом, рассматривается при xk+1 S + (x, a) как вероятность выживания. Если же направле k ние вектора выбрано так, что (xk+1 x, n(x )) 0, то с вероятностью k k a следующая точка цепи выбирается на S (x, a), а с до q(, a) = k sinh(a) полнительной вероятностью – внутри полушара B (xk, a).

Лемма 1. Среднее число выходов построенной марковской цепи на гра 2v ницу есть EN = (1 + dv). Здесь v – ограниченное решение задачи (8), a vi ve (y) + 1, а dv = O(a) при a 0.

(4) с условием на границе i (y) = e n n Пусть {x, j = 1,..., Ni } – последовательность точек выхода по k,j строенной марковской цепи из области G на границу, а {xk+1,j G, j = 1,..., Ni 1} – последовательность точек возврата цепи в эту область.

Тогда Теорема 3.

Ni [ ] g(x ) g(xk+1,j ) g(x [u](x0 ) = g(x0 ) + k,j+1 ) (11) k, j= является оценкой для решения краевой задачи (8), (4), (9). Для = (a/2R) смещение оценки равно O(a2 ) при a 0. Дисперсия оценки конечна, а трудоёмкость равна O(log() 5/2 ) при заданной точности.

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

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

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

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

Для задачи с условиями непрерывности и = 0 во всём пространстве используется представление решения в виде суммы регулярной части и объёмного потенциала g. Регулярная часть ищется в виде поверхностного 1 Сабельфельд, К. К. Векторные алгоритмы метода Монте-Карло для решения эл липтических уравнений второго порядка и уравнений Ламе / К. К. Сабельфельд // Доклады АН СССР.— 1982.— Т. 262, № 5.— С. 1076–1080.

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

µ = 0 Kµ + f. (12) e i 0. Ядро интегрального оператора K имеет Здесь e i, 0 = e + i 1 cos yy вид, где yy есть угол между внешней нормалью в точке y и 2 |y y | g вектором y y. Свободный член уравнения равен f = 0 и вычисля n(y) ется аналитически. При e i ряд Неймана для (12) сходится медленно, а в невыпуклом случае сходимость при замене ядра оператора на его мо дуль не сохраняется. Применяется аналитическое преобразование ряда с помощью замены спектрального параметра и используется конечное чис ло членов. В результате, O(q n+1 )-смещённая оценка для регулярной части решения есть (n) n li (0 )i Qi hx (yi ).

= (13) i= 1 1 1 1 (n) Здесь hx (y) =,q=, коэффициенты li определя 2 |x y| 1 + 2|0 | ются способом аналитического продолжения, Q0 = p0(y00)). Для построения f (y оценки используется марковская цепь изотропного случайного блуждания по границе {y0, y1,..., yn }. В невыпуклом случае используется либо равно вероятный выбор одного из пересечений, либо моделирование ветвящейся марковской цепи. В последнем случае Qi+1 = Qi sign(cos yi+1 yi ).

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

u = ur + us. Регулярная часть удовлетворяет уравнению Лапласа внутри G, а в G1 – уравнению e ur e 2 (ur + us ) = 0. Функция us является решением задачи с условиями непрерывности при = 0. Для его оценки используется либо алгоритм блуждания по границе, либо блуждание по сферам. Оценка для ur строится с использованием блуждания по сферам.

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

В параграфе 2.1 рассматривается уравнение общего вида wt = w (u · )w + cw + f, (x, t) H Rm (0, T ], (14) где u(x, t) = (u1, u2,..., um ) – ограниченное поле: |u| Cu, 0 – неко торый постоянный коэффициент, а c(x, t), f (x, t) – ограниченные функции.

Для решения уравнения (14) выписывается интегральное представление в виде суммы тепловых потенциалов, ядрами которых является фундамен тальное решение уравнения теплопроводности Z = Z(x x, t t ). По казано, что если поставлена задача Коши w(x, 0) = w0 (x), x Rm, то функция w и её производные по пространственным переменным удовле творяют системе интегральных уравнений Вольтерра:

W = KW + F, (15) где W = (W0, W1,..., Wm )T, K – матрично-интегральный оператор с мат рицей ядер K = {kij }m. Здесь i,j= k00 (x, t;

x t ) = Z c(x, t ), k0j (x, t;

x t ) = Z uj (x, t ), Z Z ki0 (x, t;

x t ) = c(x, t ), kij (x, t;

x t ) = uj (x, t ), xi xi i, j = 1,..., m, (16) а F0 и Fi – тепловые потенциалы, определяемые свободным членом урав нения и начальными данными.

Теорема 4. Пусть функции uj (x, t), c(x, t), f (x, t) a) непрерывны в H;

b) ограничены;

c) непрерывны по Гёльдеру относительно x равномерно по t, а функция w0 (x) непрерывна и ограничена в Rm.

Тогда 1) ряд Неймана для уравнения (15) сходится к вектор-функции W с ком понентами в L (H);

2) w(x, t) = W0 (x, t) является решением задачи Коши для уравнения (14);

3) Wi (x, t) = w/xi, т.е. i-я компонента суммы ряда Неймана совпадает с соответствующей производной нулевой компоненты.

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

Скалярная сопряжённая оценка (0, 1,..., m )T для вектора W (x0, t0 ) строится на основе рекуррентных соотношений m a0 cj j (x, t ) + F0 (xn, tn ), 0 (xn, tn ) = (17) n+1 n+ j= m ai cj j (x, t ) + Fi (xn, tn ), i (xn, tn ) = (18) n+1 n+ j= Здесь x 1/, tn+1 = tn, – выборочное n+1 = xn + 2((tn tn+1 ) m ) значение равномерно распределённой на отрезке [0, 1] случайной величи ны, m – -распределённая с параметром m/2 случайная величина, а tn – единичный вектор с изотропным случайным направлением. a0 = с q вероятностью q0 и a0 = 0 с вероятностью 1 q0.

x = xn + 2((tn t ) )1/2, t = (1 ( )2 )tn, при этом, n+1 n+1 n+ m+, и, m, попарно независимы, а весовой множитель ai = m+ ( ) 2 m+1 tn (2) с вероятностью q1 и ai = 0 с вероятностью 1 q1. Та q1 m 1 i ким образом, оценка, несмотря на линейность интегрального уравнения, строится на траекториях ветвящейся марковской цепи, каждая точка ко торой (xn, tn ) с вероятностями q0 и q1 порождает в следующем поколении независимые точки (x, t ) и (x, t ).

n+1 n+1 n+1 n+ Теорема 5. Рекуррентные соотношения (17), (18) определяют несме щённую оценку с конечной дисперсией для решения задачи Коши для урав нения (14) и для его первых производных по пространственным перемен ным в точке (x0, t0 ).

При c = 0 строится векторная оценка для градиента решения W (x, t).

При этом её можно вычислять, не оценивая собственно решение задачи.

Несмещённая сопряжённая оценка определяется соотношением N (x0, t0 ) = Q F (xn, tn ), (19) n n= где Q n+1 = An Qn, – матричные веса, Q0 = E – единичная матрица, An = {ai cj }i,j=1.

m N a0 (x0, t0 ) · Q F (xn, tn ) + F0 (x, t) 0 (x, t) = u n n= N F (xn, tn ) · Q + F0 (x, t), = (20) n n= где Q = a0 (x0, t0 ), Q = AT Q – векторные веса. K = {kij }m, n u 0 n+1 n i,j= F = (F1,..., Fm )T.

Несмещённая локально-векторная оценка для W (x, t) определяется со отношениями N K (x, t;

xn, tn )Qn + F (x, t).

(x, t) = (x, t;

x0, t0 )Q0 + F (x, t) = n= Здесь весовые матрицы Qn и векторные веса Qn определяются соотноше ниями Q0 = E, Qn+1 = AT Qn, n F (x, t ) Qn+1 = An Qn, Q0 =, p0,2 (x, t ) для согласованной плотности p0,2, ( ) 2 m+1 (t t ) j ui (x, t ) (m 2) n (An )ij = nn q2 2 с вероятностью q2 и (An )ij = 0 с вероятностью 1 q2.

Матрица вторых моментов локально-векторной оценки полностью опре деляется массивом матриц m ip = {Eij pq }i,q=1, i, p = 1,..., m.

Каждая из этих матриц есть сумма сходящегося ряда Неймана для соот ветствующего интегрального уравнения:

(K )T ip K T + Li kp + i Lp + i p, i, p = 1,..., m.

kT ip kT k = p X Здесь Li = ((W K )i1,..., (W K )im ).

Далее рассматривается первая начально-краевая задача для уравнения (14): w(y, t) = g(y, t), y, t (0, T ]. Решение этой задачи представ ляется в виде суммы тепловых потенциалов, в которую добавляется по t Z ds(y )2 µ(y, t ) с неизвестной плот верхностный потенциал dt n(y ) 0 ностью. Получена система интегральных уравнений Вольтерра со слабо полярными ядрами, которой удовлетворяет вектор-функция с компонен тами w и µ. Показана сходимость ряда Неймана для этой системы. Её ран домизация определяет ветвящуюся марковскую цепь X = {(xn, tn ), n = 0, 1,...} с фазовым пространством, включающим в себя как внутреннюю часть области решения, так и её боковую цилиндрическую поверхность.

Оценка решения начально-краевой задачи [w](x, t) строится на траекто риях цепи X в соответствии с рекуррентными соотношениями [w](xn, tn ) = a0 G (x )(div u(x, t ) + c(x, t )) [w](x, t ) n+1 n+1 n+1 n+1 n+1 n+1 n+ G (x )(u(x, t ) · a) [w](x, t ) n+1 n+1 n+1 n+1 n+ 2k j bj [µ](yn+1, tj ) + F (xn, tn ).

j + (21) n+ j= [µ](yn, tn ) = a0 G (x )(div u(x, t ) + c(x, t )) [w](x, t ) n+1 n+1 n+1 n+1 n+1 n+1 n+ G (x )(u(x, t ) · a) [w](x, t ) n+1 n+1 n+1 n+0 n+ 2k j bj [µ](yn+1, tj ) + F (yn, tn ) g(yn, tn ), j + (22) n+ j= Здесь G и – индикаторные функции области G и поверхности соот ветственно.

Теорема 6. Рекуррентные соотношения (21), (22) определяют несме щённую оценку с конечной дисперcией для решения w(x, t) первой начально краевой задачи для уравнения (14).

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

Для уравнения конвекции-диффузии во внешности ограниченной одно связной области G0 с односвязной границей wt = w (u · )w, (23) (x, t) X = G [0, T ), G = R \ G0, m m рассматривается первая краевая задача. Её решение представляется в виде суммы объёмных и поверхностного тепловых потенциалов, что позволяет получить систему интегральных уравнений Вольтерра, которой удовлетво ряют решение и плотность поверхностного потенциала:

w = Kw + Kµ + F0 (24) µ = BKw + BKµ + BF0 g. (25) Здесь ядро интегрального оператора K определено соотношением 2 k(x, t;

x, t ) · = (4(t t ))m/2+ m/ ( )[ ] |x x | · u(x, t ) · (x x ), exp (26) 4(t t ) dx Z(x x, t)w0 (x ). Оператор K с ядром F0 = G ( ) |x y | · cos y x |x y | k(x, t;

y, t ) = exp (27) 4(t t ) (4)m/2 ((t t ))m/2+ переводит функции из L (Y) (Y = [0, T )) в функции, бесконечно дифференцируемые в любой точке, отделённой от границы. Здесь y x = (n(y ), xy ), а вектор нормали n(y ) считаем для определённости направ ленным внутрь G0 (то есть n(y ) – внешняя по отношению к G нормаль).

Предполагаем также, что поверхность – ляпуновская из класса L1 (B, ).

BK – интегральный оператор с ядром k(y, t;

x, t ), определённым в (26), действующий из L (X ) в C(Y), BK – интегральный оператор с ядром k(y, t;

y, t ), действующий из L (Y) в C(Y), BF0 – граничное значение функции F0, а g – краевое условие.

Определяется ветвящаяся марковская цепь, на траекториях которой строится оценка решения:

[w](X0 ) ( )1/ ( m+1 ) [ ] t = 1 · 2 u(X1 ) · e,1 · [w](X1 ) ( m ) q 2k 1 j · sign(cos yj x0 ) [µ](Y1j ) + F0 (X0 ).

+ (28) q j= [µ](Y1 ) ( )1/ ( m+1 ) [ ] t = 2 · 2 u(X2 ) · e,2 · [w](X2 ) m ( 2 ) q 2k 1 j · sign(cos yj y1 ) · [µ](Y2j ) + F0 (Y1 ) g(Y1 ). (29) + q j= Здесь, j – индикаторные функции области G и поверхности.

Лемма 2. Обозначим n(G) (x, t) и n() (x, t) – среднее число точек вет вящейся марковской цепи, начинающейся в точке (x, t) X, в G и на соответственно. В условиях поставленной начально-краевой задачи n(G) (x, t) и n() (x, t) ограничены.

Теорема 7. Случайная величина [w](X), построенная на основе рекур рентных соотношений (28), (29), является несмещённой оценкой с ко нечной дисперсией для решения w(x, t) первой начально-краевой задачи.

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

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

Рассмотрим первую краевую задачу для уравнения теплопроводности.

Для оценки решения применяется алгоритм случайного блуждания по гра нице. В невыпуклом случае оценка строится на траекториях ветвящейся марковской цепи, переходная плотность которой есть [ ] 2| cos yi yi+1 | p(yi, ti yi+1, ti+1 ) = (30) m |yi yi+1 |m [ ( )] 4|yi yi+1 |m |yi yi+1 | · exp.

4(ti ti+1 ) ( 2 )(4(ti ti+1 ))m/2+ m Доказана следующая лемма, позволяющая обосновать реализуемость оце нок метода Монте-Карло при решении уравнений параболического типа в невыпуклых областях.

Лемма 3. Среднее число точек нестационарного ветвящегося случайного блуждания по границе равно t n= dt0 d(y0 )m(y0, t0 )p0 (y0, t0 ), 0 где m – решение интегрального уравнения t dt d(y )p(y, t y, t )m(y, t ) + 1, m(y, t) = (31) 0 p0 – плотность распределения начальной точки цепи.

В параграфе 2.5 рассматривается двумерное уравнение Прандтля для завихренности, построено интегральное уравнение, которому удовлетворя ет эта функция. Для решения этого уравнения определяется итерационный алгоритм и доказывается его сходимость. На основе рандомизации итера ционного процесса построена ветвящаяся марковская цепь и конструктив но определена монте-карловская оценка решения.

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

Построена оценка решения и доказаны её вероятностные свойства.

Рассмотрим уравнение wt = w + cw + f, (32) где коэффициент c(x, t;

), f (x, t;

) и начальные данные w0 (x;

) являются случайными полями над некоторым подходящим вероятностным простран ством (, A, P ), (x, t) DT = Rm (0, T ). Решение представляется в виде суммы тепловых потенциалов w = KZ0 c w + F0, где Z0 – фундаментальное dx Z0 (xx, t)w0 (x ;

). Это решение уравнения теплопроводности, F0 = Rm представление рассматривается как интегральное уравнение, а его решение ищется в виде суммы ряда n w= KZ0 c F0, (33) n= где предел рассматривается либо почти наверное (п.н.), либо в среднеквад ратичном (с.к.) смысле.

Марковская цепь X = {(xn, tn ), n = 0, 1,...} с фазовым пространством DT определяется начальной точкой и переходной плотностью p(xn, tn q xn+1, tn+1 ) = Z0 (xn xn+1, tn tn+1 ).

tn Cn = E(cn |X) (C0 1) – моментные функции случайного поля c(x, t;

), n C n = E(|cn | |X) и Cn = sup |Cn |, где cn (X;

) есть произведение c(xj, tj ;

) X j= значений случайного коэффициента рассматриваемого параболического урав нения в точках последовательности X.

Оценка решения определяется рекуррентными соотношениями tn [w](xn, tn ;

) = c(xn+1, tn+1 ;

) [w](xn+1, tn+1 ;

) sn q + [F0 ](xn, tn ;

), (34) где k(xn, tn ;

xn+1, tn+1 ) Q [c] = Q [c] p(xn, tn xn+1, tn+1 ) n+1 n tn = Q [c] c(xn+1, tn+1 ;

), Q [c] = 1, n q q – вероятность выживания на переходе (xn, tn ) (xn+1, tn+1 ), а {sn, n = 0, 1,...} – последовательность случайных величин, равных единице с веро ятностью q и нулю с дополнительной вероятностью 1 q, N – случайное число членов в выборочной траектории марковской цепи X.

1) если для почти всех фиксированных функции c, Теорема 8.

f и w0 ограничены почти всюду в DT, то п.н. существует условное математическое ожидание E( [w]|), равное сумме ряда Неймана из (33);

2) условное математическое ожидание E( [w]|) является п.н. реше нием случайной задачи Коши и условная дисперсия этой оценки ко нечна;

3) для произвольного конечного t 0, [w](x, t;

) является несмещён ной оценкой для математического ожидания случайного решения w(x, t) Ew(x, t;

·) ;

4) если кроме того f и w0 имеют ограниченные вторые моменты, а мо ментные функции случайного коэффициента c удовлетворяют усло виям C n n!!q n/2 tn a(n) для некоторого сходящегося ряда a(n), то дисперсия оценки [w](x, t;

) конечна.

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

Глава 3 посвящена применению построенных методов статистического моделирования к решению модельных и прикладных задач.

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

Константа реакции макромолекулы G с броуновскими частицами опреде ляется как интегральный поток этих частиц на поверхности молекулы.

Задача вычисления диффузионно обусловленной константы реакции K сводится в стационарном случае к решению внешней краевой задачи для уравнения Лапласа:

u(x) = 0, x G1 = R3 \ G, u s (y)u(y) D (y) = s (y), y G, (35) n lim u(x) = 0.

|x| Здесь u = 1. (x) – концентрация частиц, D – коэффициент диф фузии. Вычисления основаны на формуле K = 4RD u(R), где u(r) = d 4 u(r, ), а R таково, что G B(x, R) для некоторого x. Рандомиза ция этой формулы приводит к следующей монте-карловской оценке:

K = 4RD E (x0 ). (36) Здесь случайная точка x0 равномерно распределена на S(x, R), а – неза висимая от неё оценка решения задачи (35). Эта оценка строится на тра екториях случайного блуждания по сферам с прямым моделированием ухода на бесконечность (метод А.С.Сипина). При первом попадании в окрестность границы используется рандомизация конечно-разностного при ближения к граничному условию третьего рода. Для частного случая сме шанной краевой задачи (модель Solc-Stockmayer) или задачи Неймана ис пользуется алгоритм, основанный на рандомизация формулы среднего, по строенный в первой главе. Результаты тестовых расчётов показывают эф фективность предложенных алгоритмов в сравнении с использовавшими ся ранее оценками, основанными на прямом моделировании броуновского движения.

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

Рис.1).

Рис. 1: Зависимость от (ширины приграничной полосы, в лог-лог мас штабе) среднего числа переходов в блуждании по сферам, в котором от ражение от границы с условиями Неймана производилось в соответствии с соотношением о сферическом среднем (7) (ромбы). Результаты расчётов аппроксимируются зависимостью EN = 4.419 ln 9.505, то есть EN име ет тот же порядок, что и при решении задачи Дирихле.

Для сравнения приведена зависимость среднего числа переходов при ис пользовании конечноразностного приближения к нормальной производ ной (кресты). Результаты расчётов аппроксимируются степенной зависи мостью EN = 0.691 h1.059. Здесь h = 1/2 – шаг аппроксимации.

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

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

Общепринятым подходом, позволяющим вычислять электростатические свойства макромолекул в растворе, является использование физической модели, в которой только структура самой молекулы (область Gi R3 ) описывается точно, а окружающая её вода с растворёнными в ней солями рассматривается как сплошная среда. При такой постановке задачи элек тростатический потенциал i внутри молекулы удовлетворяет уравнению Пуассона M Qm (x x(m) ) = 0, x Gi, i i (x) + 4 (37) m= где Qm – величины зарядов, расположенных в точках x(m), i – диэлек трическая постоянная в Gi. В растворе за пределами молекулы (область Ge = R3 \ Gi ) электростатический потенциал e удовлетворяет уравнению Пуассона-Больцмана. При малых значениях потенциала это уравнение ли неаризуется:

e (x) 2 e (x) = 0. (38) Значения функций i и e согласованы условиями непрерывности ви да (9). Для вычисления значений потенциала и электростатической энер гии используются оценки вида (11), описанные в первой главе. При этом g = c (x) – сингулярная составляющая потенциала, представимая в ана литическом виде. Здесь i (x) = rf (x) + c (x). Для описания геометрии молекулы рассматривается модель, определяемая поверхностью Ван-дер Ваальса.

Постановка задачи определения электростатических свойств макромо лекул такова, что применение методов Монте-Карло для её решения явля ется естественным. Использование построенных в данной работе алгорит мов позволяет, не находя решения во всём пространстве, вычислять зна чения потенциала в отдельных точках. Не возникает проблем, вызванных неограниченностью расчётной области, в том числе, проблем с компьютер ной памятью. Статистический характер погрешности позволяет адекватно оценивать точность полученного решения. Разработанные алгоритмы ста тистического моделирования позволяют вычислять значения потенциала диэлектрической реакции и энергии одновременно для конечного набора значений параметра, квадрат которого линейно зависит от концентрации соли типа N aCl:

Nins [rf ](x(m) ) = c (x ) + Fj () (c (xins ) c (x )). (39) 1 j j,ins j= Здесь весовая функция Fj () 1 при = 1, а при других значениях параметра домножается на каждом шаге во внешней области на отноше ние q(i, dk )/q(1, dk ). Оценка для свободной электростатической энергии диэлектрической реакции есть линейная комбинация оценок точечных зна чений rf :

M Qm [rf ](x(m) ).

[Wrf ] = (40) 2 m= При этом значения потенциала оцениваются непосредственно в местах рас положения зарядов, что невозможно при использовании конечноразност ных и конечноэлементных методов. Тестовые расчёты для модельной зада чи, имеющей аналитическое решение, показали, что данное обстоятельство приводит к относительно большим ошибкам в результатах. Например, для концентрации соли 0.1 моль расчёт по широко используемой программе APBS с шагом сетки 0.1 даёт результат с ошибкой 0.856%, в то время как оценка (39) приводит к результату со смещением 0.028%. При i = 1 энер гия диэлектрической реакции Wrf совпадает с электростатической сво бодной энергией растворения. На Рис. 2 показана зависимость разности между энергией при заданной концентрации соли и её значением при кон центрации 0.5 моль для пептидов которые идентифицируются в PDB (бел ковой базе данных) как (A) 1A4T и (B) 1QFQ. Для сравнения приведены результаты решения, полученные на основе метода граничных элементов с ускорением по мультипольному методу. Существенно, что каждое зна чение, полученное с использованием детерминированного метода, требует отдельного расчёта, в то время как алгоритм статистического моделиро вания позволяет получить всю кривую на основе одного вычислительного эксперимента. Результаты вычислений по методу Монте-Карло сильно кор релированы, и дисперсия разности результатов, полученных для 1 и стремится к нулю при стремлении к нулю разности |1 2 |.

В параграфе 3.4 описан комплекс программ, построенный с использо ванием разработанных алгоритмов и применяющийся для решения задач молекулярной биофизики. Программа написана на языке Фортран. При её создании использован опыт работы над предыдущими пакетами программ, реализующими алгоритмы блуждания по сферам в сочетании с блужда нием в подобластях, описанным в первой главе. Применялись все идеи, Рис. 2: Разность между электростатической энергией растворения и кон трольным значением энергии при концентрации NaCl в 0.5 M для двух пептидов: 1A4T и 1QFQ. Диэлектрические постоянные равны 1 внутри и 78.5 вовне молекулы. В качестве границы используется поверхность Ван дер-Ваальса. Все детерминированные расчёты проводились с использова нием метода граничных элементов с ускорением по мультипольному ме тоду. Результаты обозначены треугольниками (1A4T) и ромбами (1QFQ).

Показаны доверительные интервалы (±2) для результатов, полученных методом Монте-Карло.

реализованные ранее при разработке программного обеспечения для бо лее ранних версий операционных систем. Были учтены конкретные осо бенности геометрических моделей, применяемых в задачах определения электростатических свойств макромолекул в растворе. Внутри молекулы используется блуждание в подобластях, а для нахождения расстояния до границы во внешней для молекулы области – алгоритм, основанный на построении многомерных двоичных деревьев. Созданный на основе раз работанных алгоритмов программный комплекс используется в настоящее время для исследовательских расчётов в Институте молекулярной биофи зики Флоридского государственного университета.

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

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

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

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

СПИСОК РАБОТ ПО ТЕМЕ ДИССЕРТАЦИИ Монографии 1. Курбанмурадов О. А., Сабельфельд К. К., Симонов Н. А. Алгорит мы случайного блуждания по границе. — Новосибирск: ВЦ СО АН СССР, 1989.

2. Sabelfeld K. K., Simonov N. A. Random Walks on Boundary for solving PDEs. — Utrecht, The Netherlands: VSP, 1994.

Публикации в журналах, рекомендованных ВАК 3. Симонов Н. А. Модификация алгоритма блуждания по сферам, удоб ная в практическом применении // Журн. вычисл. матем. и матем.

физ. — 1984. — Т. 24, № 6. — С. 936–938.

4. Сабельфельд К. К., Симонов Н. А. Решение пространственных задач теории упругости в детерминированной и стохастической постанов ках методом Монте-Карло // Доклады АН СССР. — 1984. — Т. 275, № 4. — С. 802–805.

5. Simonov N. A. Monte Carlo solution of the nonlinear integral equation system in the boundary layer theory // Russian Journal of Numerical Analysis and Mathematical Modelling. — 1993. — Vol. 8, no. 3. — Pp. 265– 274.

6. Симонов Н. А. Стохастические итерационные методы решения урав нений параболического типа // Сиб. матем. журн. — 1997. — Т. 38, № 5. — С. 1146–1162.

7. Simonov N. A. Monte Carlo methods for convective diusion equation // Russian Journal of Numerical Analysis and Mathematical Modelling. — 1997. — Vol. 12, no. 1. — Pp. 67–81.

8. Simonov N. A. Monte Carlo solution to three-dimensional vorticity equa tion // Mathematics and Computers in Simulation. — 1998. — Vol. 47, no.

2-5. — Pp. 455–459.

9. Simonov N. A. Stochastic iterative method for a system of parabolic equations // Zeit. Angew. Math. Mech. — 1998. — Vol. 78, Suppl.1. — Pp. S185–S188.

10. Mascagni M., Simonov N. A. Monte Carlo methods for calculating some physical properties of large molecules // SIAM Journal on Scientic Computing. — 2004. — Vol. 26, no. 1. — Pp. 339–357.

11. Mascagni M., Simonov N. A. The random walk on the boundary method for calculating capacitance // The Journal of Computational Physics. — 2004. — Vol. 195, no. 2. — Pp. 465–473.

12. Симонов Н. А. Методы Монте-Карло для решения эллиптических уравнений с граничными условиями, включающими в себя нормаль ную производную // Доклады Академии наук. — 2006. — Т. 410, № 2. — С. 164–167.

13. Simonov N. A., Mascagni M., Fenley M. O. Monte Carlo-based linear Poisson-Boltzmann approach makes accurate salt-dependent solvation free energy predictions possible // Journal of Chemical Physics. — 2007. — Vol. 127. — P. 185105 (6 pages).

14. Using Correlated Monte Carlo Sampling for Eciently Solving the Linear ized Poisson-Boltzmann Equation Over a Broad Range of Salt Concentr ation / M. O. Fenley, M. Mascagni, J. McClain et al. // J. Chem. Theory Comput. — 2010. — Vol. 6, no. 1.— Pp. 300–314. Publication Date (Web):

December 23, 2009. http://dx.DOI.org/10.1021/ct9003806.

Прочие публикации 15. Симонов Н. А. Алгоритмы случайного блуждания для решения крае вых задач с разбиением на подобласти // Методы и алгоритмы стати стического моделирования. — Новосибирск: ВЦ СО АН СССР, 1983. — С. 48–58.

16. Симонов Н. А. О рандомизации итерационных процессов решения уравнений первого рода // Теория и приложения статистического мо делирования. — Новосибирск: Вычислительный центр СО АН СССР, 1985. — С. 31–37.

17. Симонов Н. А. Решение граничного интегрального уравнения для третьей краевой задачи методом Монте-Карло // Методы статисти ческого моделирования. — Новосибирск: Вычислительный центр СО АН СССР, 1986. — С. 53–59.

18. Симонов Н. А. Случайные оценки производных от решения задачи Неймана вблизи границы области // Теория и приложения статисти ческого моделирования. — Новосибирск: Вычислительный центр СО АН СССР, 1988. — С. 76–87.

19. Назначение и архитектура пакета прикладных программ НЕКТОН 1 / М. М. Бежанова, И. И. Кабанихина, С. Е. Макаров и др. // VIII Всесоюзное совещание "Методы Монте-Карло в вычислительной ма тематике и математической физике 19-21 февраля 1991г. Часть 1. — Новосибирск: ВЦ СО АН СССР, 1991. — С. 142–145.

20. Симонов Н. А. Монтекарловская оценка решения системы интеграль ных уравнений, порожденной уравнениями Прандтля // Труды ВЦ СО РАН. — Новосибирск: Вычислительный центр СО РАН, 1993. — Вычислительная математика. Выпуск 1. — С. 107–112.

21. Simonov N. A. Solution of two-dimensional Prandtl equations by the Monte Carlo method // Bulletin of the Novosibirsk Computing Center. — Novosibirsk: Computing Center SB RAS, 1993. — Numerical Analysis.

Issue 4. — Pp. 83–102.

22. Simonov N. A. Stochastic algorithm for solving two-dimensional boundary layer equations // Siberian Journal of Computer Mathematics. — 1995. — Vol. 2, no. 1. — Pp. 41–56.

23. Simonov N. A. Boundary value problem and stochastic algorithm for two-dimensional Navier-Stokes equations // Monte Carlo Methods and Applications. — 1995. — Vol. 1, no. 1. — Pp. 59–70.

24. Симонов Н. А. Решение методом Монте-Карло первой краевой задачи для параболических уравнений с граничной функцией, зависящей от решения // Труды ВЦ СО РАН. — Новосибирск: Вычислительный центр СО РАН, 1996. — Вычислительная математика. Выпуск 4. — С. 129–145.

25. Simonov N. A. Monte Carlo solution of the rst boundary value problem for parabolic equation // Mathematical methods in stochastic simulation and experimental design. Proc. of the 2nd St. Petersburg Workshop on simulation, St. Petersburg, June 18-21, 1996 / Ed. by S. M. Ermakov, V. B. Melas. — St. Petersburg: Pub. House of St. Petersburg University, 1996. — Pp. 109–111.

26. Simonov N. A. Stochastic computational methods for parabolic equations with random data // IMACS World Congress. Berlin, Germany, August 24-29, 1997. Proceedings. Volume 1. Computational Mathematics. — 1997. — Pp. 449–452.

27. Iterative procedure for multidimensional Euler equations / W. Dreyer, M. Kunik, K. Sabelfeld et al. // Monte Carlo Methods and Applications. — 1998. — Vol. 4, no. 3. — Pp. 253–271.

28. Simonov N. A. Nonlinear Monte Carlo estimators for parabolic equation // Simulation 2001, Proc. 4th St. Petersburg Workshop on Simulation / Ed.

by S. M. Ermakov, Y. N. Kashtanov, V. B. Melas. — St. Petersburg: St.

Petersburg State University, 2001. — Pp. 453–456.

29. Simonov N. A. Monte Carlo solution of a parabolic equation with a random coecient // Сибирский журнал вычислительной матема тики. — 2001. — Т. 4, № 4. — С. 389–402.

30. Mascagni M., Simonov N. A. Random walk on the boundary methods for computing reaction rate and capacitance // The International Conference on Computational Mathematics. Novosibirsk, Russia, June 24-28, 2002, Proceedings / Ed. by G. A. Mikhailov, V. P. Il’in, Y. M. Laevsky. — Novosibirsk: ICM&MG; Publisher, 2002. — Pp. 238–242.

31. Mascagni M., Simonov N. A. Monte Carlo method for calculating the electrostatic energy of a molecule // Lecture Notes in Computer Science. — 2003. — Vol. 2657. — Pp. 63–74.

32. Karaivanova A., Mascagni M., Simonov N. A. Solving BVPs using quasi random walks on the boundary // Lecture Notes in Computer Science. — 2004. — Vol. 2907. — Pp. 162–169.

33. Simonov N. A., Mascagni M. Random walk algorithms for estimating eective properties of digitized porous media // Monte Carlo Methods and Applications. — 2004. — Vol. 10, no. 3-4. — Pp. 599–608.

34. Karaivanova A., Mascagni M., Simonov N. A. Parallel quasirandom walks on the boundary // Monte Carlo Methods and Applications. — 2004. — Vol. 10, no. 3-4. — Pp. 311–320.

35. Simonov N. A., Mascagni M. Random walk algorithms for estimating electrostatic properties of large molecules // The International Conference on Computational Mathematics. Novosibirsk, Russia, June 21-25, 2004, Proceedings / Ed. by G. Mikhailov, V.P.Il’in, Y.M.Laevsky. — Novosibirsk:

ICM&MG; Publisher, 2004. — Pp. 352–358.

36. Fleming C., Mascagni M., Simonov N. An ecient Monte Carlo approach for solving linear problems in biomolecular electrostatics // Lecture Notes in Computer Science. — 2005. — Vol. 3516. — Pp. 760–765.

37. Simonov N., Mascagni M. The method of random walk on spheres for solving boundary-value problems from molecular electrostatics // 17th IMACS World Congress, Paris, France, July 11-15, 2005. Proceedings. — 2005. — Pp. T1–I–62–0945. — ISBN 2-915913-02-1.

38. Karaivanova A., Simonov N. A. Quasi-Monte Carlo methods for investi gating electrostatic properties of organic pollutant molecules in solvent // Lecture Notes in Computer Science. — 2006. — Vol. 3743. — Pp. 172–180.

39. Симонов Н. А. Алгоритмы случайного блуждания по сферам для решения смешанной краевой задачи и задачи Неймана // Сибир ский журнал вычислительной математики. — 2007. — Т. 10, № 2. — С. 209–220.

40. Simonov N. A. Random walks for solving boundary-value problems with ux conditions // Lecture Notes in Computer Science. — 2007. — Vol.

4310. — Pp. 181–188.

41. Simonov N. A. Walk-on-spheres algorithm for solving boundary-value problems with continuity ux conditions // Monte Carlo and Quasi Monte Carlo Methods 2006 / Ed. by A. Keller, S. Heinrich, H. Niederreiter. — Heidelberg: Springer-Verlag, 2008. — Pp. 633–644.

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

Отпечатано в типографии OOO Омега Принт 630090 г.Новосибирск, пр. Академика Лаврентьева, 6, оф. 3- тел./факс. (383) 335 65 23, e-mail: [email protected]