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

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

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

РГБ ОД

1 5 ДЬК 1996 БЕЛЯЕВА

Аэлита Аркадьевна

МЕТОДЫ МОНТЕ-КАРЛО ДЛЯ РЕАЛИЗАЦИИ ИТЕРАЦИОННЫХ ПРОЦЕССОВ

Специальность 01.01.07. - вычислительная математика

АВТОРЕФЕРАТ

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

САНКТ-ПЕТЕРБУРГ 1996

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

Научный руководитель- доктор физико-математических наук,

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

Официальные оппоненты- доктор физико-математических наук,

вед. научный сотрудник КАРКЖИН Владимир Владимирович доктор физихо-математических наук, профессор СЕДУНОВ Евгений Витальевич

Ведущая организация- Санкт-Петербургский Технический

Университет

Зашита состоится 1996 г. в _ часов на заседании

диссертационного совета Д 063.57.30 по защите диссертаций на соискание учёной степени доктора физико-математических наук в Санкт-Петербургском Государственном Университете по адресу: 198904 Санкт-Петербург, Старый Петергоф, Библиотечная площадь, 2.

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

ОН

Автореферат разослан МХШ 1996 г.

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

диссертационного совета Д 063.57.30 кандидат физ.-мат. наук

Ю. А. Сушков

Общая характеристика работы.

Актуальность темы. Как известно, метод Монте-Карло применяется для решения широкого круга задач. Он успешно конкурирует с другими вычислительными методами, а в ряде случаев этот метод является единственно возможным для получения ответа за разумное машинное время. Наиболее изученной и часто применяемой схемой решения уравнений методом Монте-Карло является классическая схема, предложенная Дж. Фон Нейманом и Уламом. Её применение, однако, возможно лишь при выполнении довольно ограничительных условий. Очевидно, что построение новых алгоритмов метода Монте-Карло, при применении которых возможно ослабить эти условия является актуальной задачей. Имеет смысл также модифицировать классические схемы метода Монте-Карло, реализующие различные итерационные процессы, с целью расширения класса решаемых с его помощью задач. Наиболее перспективные схемы, отличные от схемы Неймана-Улама носят -1 последовательный характер. Примеры использования таких схем имеются у Холтона и учёных новосибирской школы метода Монте-Карло. В последовательных схемах важное значение имеет стохастическая устойчивость алгоритма. Впервые понятие стохастической устойчивости введено С.М. Ермаковым. Ивановой A.B. получены условия применимости алгоритмов "послойного" счёта для эволюционных уравнений математической физики. Для стационарных задач эта проблема не изучалась.

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

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

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

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

Апробация работы. Результаты работы обсуждались на семинарах кафедры статистического моделирования и лаборатории моделирования систем и статистических методов Санкт-Петербургского Го^дарственного Университета, на семинаре "Дифференциальных и интегральные уравнения" Новгородского Государственного Университета. Основные результаты докладывались на конференции по теории вероятностей и математической статистике ( сентябрь 1995 г., Фергана), на Второй Международной конференции по моделированию (июнь 1996 г., Санкт-Петербург), Публикации. Основные результаты работы отражены в работах [1-4]. Структура и объём работы. Диссертация состоит из введения и трёх глав. Объём основной части диссертации составляет ^^^страниц. Библиография содержит 43 наименования.

Содержание работы.

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

Суть задачи в самом общей виде заключается в следующем. Пусть задано некоторое, вообще говоря, нелинейное уравнение вида <р = Р(ф), где ф-функция переменной хеХсЗГ из класса функцийФ, заданный на Ф оператор. Предположим. что /"(р) = у/) и сходится метод последовательных приближений

?>„! 4 = 0,1,к, р0 - задано (1)

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

= />Д„«)> £0 = <р0еФ, (2)

где ^„„¡(я),-оценки соответственно (п + 1)-го и (п)-го приближения задачи (1) а-случайная точка некоторого вероятностного пространства (0,3, Р) (например, траектория цепи Маркова), /„{<&,-)- измеримые случайные отображения этого пространства в пространство Х„, с X с сг- алгеброй .

Считаем, что при каждон п мы имеем выборку независимых случайных величин /,(е>,,5,), /я(о2,4).К. ,/„(<";/■£,) такую, что распределение выборочного среднего этой выборки, определенного как

. (4)

п ы

можно считать достаточно близким к нормальному.

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

Определи яле. Алгоритм (2) назовём статистически устойчивым, если существует ограниченная функция гдел еХсШ'.такая, что

£%„(х) < Ч/(х) для всех п=1,2,... (5)

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

Х=АХ+Р , (б)

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

Предположим, что сходится метод последовательных приближений

Х„ = АХ^+Р. (7)

Во втором к четвёртом параграфах первой главы описываются алгоритмы построения случайных оценок для решения задачи (6). На каждом шаге итерационного процесса (7) моделируется случайный вектор Е"(Лгп) = (^"(Л7л),4г(Лг«)>"-^"(Лг»)) являющийся оценкой Ха итерационного решения уравнения (6). Последовательность случайных векторов Е" построена так, чтобы условное математическое ожидание Н" при фиксированных .а, было равно

Е{Н„ /Ел.1,Н^>...,51) = АЕ^ +/' . (8)

Тогда для полного математического ожидания имеем

= = = ^ . (9)

Расчетные формулы первого алгоритма, задающие последовательность рекуррентных соотношений типа (2) следующие:

I

. СО)

где, =

^иаИГЧЛ,-,) ^па^'^)

4 =

А А А

Из формулы (10) видно, что вычисление этой оценки требует моделирования дискретной случайной величины и приписывания соответствующего знака величине ¡^(Л1^,), гДе З-случайный номер полученный в результате моделирования.

Расчётные формулы второго алгоритма, также удовлетворяющего требованиям (8), (9) следующие:

. пи

'»Л 1=1

„ . , „ гГ'№-,> <4 »-I где ^>=7/ +/, , ц = —• . Рд -начальная вероятность,/^,-вероятность Л,

перехода из ^ в /.

Построение вектора Х„(№„) требует моделирования на каждом шаге "двухзвенной цепи Маркова" с начальным состоянием, распределение которого ме!яется от шага к шагу, так как зависит от оценки, полученной на предыдущем шаге и матрицей переходных вероятностей Р, согласованной с матрицей А.

Поведение дисперсии оценки 1-го алгоритма описывается следующим рекуррентным соотношением, связывающим дисперсию оценки п-ой итерации с дисперсией оценки предыдущей итерации :

ввчл.) = в£ (Ю+Е^-'Щ-,)

. , (12) ■

+ 2(1соу^К-^-Ч^.,))

где -обозначает дисперсию оценки для А', итерации, получаемой при

т

условии, если бы Хп_х было известно без случайной погрешности, А, = '^Г1\ац\ •

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

Теорема 1. Если спектральная норма матрицы А удовлетворяет условию ¡¡А||<ц<1, то оценка итерационного решения уравнения (6), получаемая по схеме (10) или (11) будет стохастически устойчива при выполнении не менее ^

испытаний на каждом шаге игераций, где У. > ^ - ^^ а оператор А,

НИИ

зависит от вида используемой оценки((10) или (11)).

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

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

В том случае, когда метод последовательных приближений (7) расходится или сходится медленно, рассматривают методы отличные от (7), сходящиеся достаточно быстро к решению. Оказывается можно построить и вероятностные методы решения (6), основанные на методах последовательных приближений, отличных от (7). В главе 2 подробно рассмотрена и исследована реализация метода Монте-Карло для решения (6), основанная на методе Зейделя. Метод Зейделя предполагает представление исходной матрицы А задачи (б) в виде суммы двух матриц А = ^ + А2< где Д-нижняя треугольная матрица с нулями на диагонали, Д-верхняя треугольная матрица, и рассмотрение итерационного процесса вида:

(13)

Метод последовательных приближений (15) эквивалентен методу

Непосредственное применение схемы Неймана-Улама к (16) требует, чтобы спектральный радиус модуля матрицы Г=(1 -Д)""'Д был меньше единицы: р([Л) < ] . Это довольно жесткое требование, а также необходимость его проверки, можно снять, если рассмотреть некоторые модификации схемы Неймака-Улама, основанные на идее пошагового запоминания промежуточных результатов вычислений . А именно, на каждом шаге итераций строится оценка решения уравнения = ДА^, + , гдеф„й =АгХп .

Это означает, что методом Монте-Карло обращается оператор(/-Д)~' , а в свободном члене используется предыдущее итерационное приближение. Такая процедура применения метода Монте-Карло предполагает выполненным условие

Рв40<1- <15)

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

Теорема 2. Пусть система (6) решается методом Монте-Карло, с запоминанием промежуточных результатов вычислений, реализующим метод Зейделя (13). Если выполнено условие (15) и матрица Р2 - является стохастической и согласованной с Аг ,то оценки

,,, <?,г(х), г./ с6)

Р р-\ Ра, —Р',-,!, Р Рй, -Л

гд^ = (!„ 1) >а = | р£> — , а*о при фГ'("/)*о

V р ' ( «я, ^ = 2

являются несмещенными оценками решения уравнения (6).

Данный алгоритм построения оценок решения с учётом структуры матриц А1;Аг означает

во-первых, моделирование типа частицы (случайная величина %); во-вторых, случайное блуждание этой частицы влево до положения гк, в котором происходит либо поглощение ( если % = 2), либо ещё один переход , но вправо, и поглощение (когда % = \).

Теорема 3. Пусть система (б) решается методом Монте-Карло, с запоминанием промежуточных результатов вычислений, реализующим метод Зейделя. Тогда, если спектральная норма оператора Т удовлетворяет условию ||(/—Д)~'у42!|<ц<1, и на каждом шаге производится не менее ^ испытаний, где

,, ЛЛ1НИНИ'1Н1Л1* - -

- ^ -—, то данный алгоритм стохастически устойчив. Здесь

оператор/? зависит от вида используемой оценки. Далее в главе 2 рассмотрен и исследован алгоритм метода Монте-Карло для решения линейных интегральных уравнений 2-го рода

й

ф) = /К(х, х') Р(х') йх + /(х) , (17)

а

реализующий абстрактный метод Зейделя

л Ь '

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

= ^/в^^ЛсЧв-'. (19)

№

, ест 2 = 1

где я-вероятность поглощения частицы. Метод Зейделя для уравнения (19) имеет вид

Показано, что норма оператора шага при уе[е,1], гдеО<Б<1 строго меньше единицы. Следовательно, имеет место стохастическая устойчивость алгоритма. Количество испытаний на каждом шаге оценено снизу величиной

N > ^ + ^ ,Из этой оценки следует, что чем меньше вероятность поглощения, 1-9

тем больше испытаний необходимо провести при каждом фиксированном п.

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

У = Л(У,У) + ВУ+Р , (20)

где У-Б-мсрный вектор евклидова пространства З?1, ^-линейный оператор, действующий из Я' еЗТ'.Хх X—ЭТ5, Хс К5 -билинейный непрерывный оператор. Для решения уравнения (20) рассмотрен метод последовательных приближений

+ У0=0, (21)

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

С(У„) = (1-В-А(У„,-)Г'Г. (22)

Теорема 4. Пусть метод последовательных приближений (21) сходится и на каждом итерационном шаге задачу (21) можно решать методом Монте-Карло, используя предыдущее итерационное приближение. Тогда данный алгоритм стохастически устойчив, если норма производной Фреше оператора шага (22) удовлетворяет условию ЦО'(Х)||<|д.<1, и при каждом фиксированном и выполняется не менее испытаний, где

||Л"(Х)/2|Н|С'(Х)Ц|1С'(Х)-Ц

Н|б'(Х)|Ц|0'(Х)'|| ' 1 '

При дискретизации нелинейных задач математической физики часто приходится решать конечно-разностные уравнения вида :

Xi+i = + BYr+l +CY„ + F, (24)

где операторы А,В определены как и выше , С-линейный оператор, a -

могут иметь смысл значений искомой величины на n-ом и n+l-ом временных слоях соответственно. Для каждого п уравнение (24) при выполнении ряда условий, эквивалентно бесконечной системе линейных уравнений, оператор которой и свободный член будут зависеть от Y„. К данной системе применяется метод Монте-Карло и для него также получено аналогичное теореме 4 условие стохастической устойчивости.

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

Публикации по теме диссертации :

[1] Ермаков С.М., Беляева А.А. О стохастической устойчивости метода Монте-Карло с запоминанием промежуточных результатов. // Тезисы докладов конференции по теории вероятностей и мат. статистике Фергана, 27-29 сентября 1995 г.,с.34-35.

[2] Belyaeva А.А., Ermakov S.M. On the stochastical stability of the Monte-Cario algorithms for solving non-linear equations. // 2-nd St.-Petersburg Workshop on Simulation, June 18-21, 1996, p.36-37

[3] Беляева А.А., Ермаков C.M. О методе Монте-Карло с запоминанием промежуточных результатов И Вестник Санкт-Петербургского Университета, серия: мат., мех., астр., 1996 г., вып.З., № 15,с.8-14.

[4] Беляева А.А. Метод Монте-Карло и метод Зейделя.// Деп. в ВИНИТИ, ífo 2673-В96, от 13.08.96.