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

Введение.

1. Оценка по времени для вычисления линейных функционалов от концентрации траекторий многомерных диффузионных процессов

1.1. Оценка по времени.

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

1.3. Реализация оценки по времени в случае однородной диффузии

1.4. Способы повышения порядка детерминированной погрешности

2. Оценка градиента решения стационарного диффузионного уравнения методом Монте-Карло на основе его вероятностного представления

2.1. Построение вспомогательного веса

2.2. Алгоритмы оценки производной ^г—.

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

3.1. Оценки на основе "блуждания по сферам и в шарах"

3.2. Асимптотическая несмещенность и равномерная ограниченность дисперсии оценок.

3.3. Комбинированная оценка для случая невырождающегося конвективного слагаемого.

3.4. Трехмерный случай.

4. Численные результаты

4.1. Оценка по времени.

4.2. Численное исследование порядка детерминированной погрешности

4.3. Вычисление производных.

4.3.1. Блуждания по сферам и в шарах.

4.3.2. Метод Эйлера.

4.3.3. Один контрпример.

Методы Монте-Карло (методы статистического моделирования) -"численные методы решения математических задач при помощи моделирования случайных величин и статистической оценки их характеристик" [27]. Их применение для решения разнообразных задач в различных областях физики, математики и техники имеет богатую историю [17, 27, 18]. Официальной датой рождения метода Монте-Карло принято считать 1949 год, когда была опубликована статья [39] С. Улама и Н. Метрополиса. Впрочем, сам термин был предложен еще во время Второй мировой войны выдающимися учеными XX века математиком Дж.фон Нейманом и физиком Энрико Ферми в Лос-Аламосе (США) в процессе работ по ядерной тематике. В 50-х годах состоялись несколько симпозиумов по методам Монте-Карло, которые продемонстрировали широкие перспективы их применения. Хотя методы Монте-Карло были известны и до 40-х годов, они не получили широкого распространения из-за больших объемов вычислений. Появление компьютеров сделало их практически применимыми и за последние 50 лет вместе с развитием компьютерных технологий методы Монте-Карло всё более активно используются во многих научных областях (теория переноса, теория массового обслуживания, теория надежности, статистическая физика и др.).

Основными преимуществами данных методов являются

• физическая наглядность и простота реализации,

• малая зависимость трудоемкости задачи от её размерности,

• возможность решения задач со сложной геометрией,

• одновременное оценивание вероятностной погрешности оценки искомого функционала,

• простое распараллеливание методов.

В диссертационной работе рассматривается ¿-мерная задача Дирихле для эллиптического уравнения

1 к д^и к ди !

2 и=1 ду1дуз дщ |г в ограниченной области П с границей, которая предполагается одно-связной и регулярной. Граница будет регулярной, если, например, для каждой точки а Е Г существует шар Вр(го) = {г Е В,к : | г* — у о | < р} такой, что а — Г2 П Вр(го). Отметим, что уравнение (0.1) является стационарным сопряжённым уравнением.

Теорема 0.1 [21]. Пусть Ь - эллиптический оператор из (0.1) такой, что к к - £ Ьцгщ + г ]Г VjZj + с 1,3=1 з=1

1/(1 + |г|2), = г2 = — 1. (0.2)

7=1

Пусть, кроме того, функции с(-) Е С6(Ик), то есть удовлетворяют условию Гёлъдера с показателем 5 в Як, а д(-) Е С*(П); Г - регулярная граница ограниченной области О, С Як, функция ф(-) непрерывна на Г. Тогда существует единственное решение и Е С?+г{П) П С(П) задачи (0.1). □

Здесь - пространство всех функций $ таких, что произведение д ■ (р принадлежит пространству Гёльдера С2^6(11к) для любой бесконечно дифференцируемой функции (р равной нулю вне Пив окрестности дО, = Г.

Один из способов решения первой краевой задачи (0.1) методами статистического моделирования состоит в сведении ее к интегральному уравнению, удовлетворяющему некоторым условиям, и в построении несмещенных оценок решения на траекториях сходящейся марковской цепи, связанной с полученным интегральным уравнением естественным способом (см., например, [18]). Например, алгоритм "блуждания по сферам" основан на моделировании точек последовательных выходов винеровского процесса на максимальные сферы, целиком лежащие в рассматриваемой области; он был впервые предложен в [14] для решения задачи Дирихле для уравнения Лапласа. Этот алгоритм является частным случаем использования формул Грина для стандартных областей, целиком лежащих в данной области (шар, сфера, эллипсоид и т. д.). В диссертационной работе (с помощью центральной и нецентральной функций Грина для оператора Лапласа в шаре) данная схема была реализована для записи локальных интегральных уравнений не только на само решение исходной дифференциальной задачи, но также и на его градиент, что позволило построить несколько алгоритмов "блуждания по сферам и в шарах" для одновременной оценки как решения, так и его производных по пространственным координатам.

Другой способ состоит в использовании тесной связи эллиптических операторов с диффузионными процессами (см., например, [37]) и основан на вероятностном представлении решения исходной дифференциальной задачи в виде функционала от решения соответствующей системы стохастических дифференциальных уравнений (СДУ). Дифференциальному оператору Ь из (0.1) можно поставить в соответствие систему СДУ в смысле Ито [33]: г = 6 + / + / (7й)<мо, (о.з) о о где £ Як,ст(г) - нижняя треугольная матрица диффузии, определяемая разложением Холесского [19]:

В = а(г)-ат(г), а ги(£) - стандартный винеровский процесс.

Теорема 0.2 [33]. Пусть О С Як - ограниченная область с регулярной границей, - момент первого выхода решения & системы СДУ (0.3) на границу 80, = Г; ф Е С(дО), функции удовлетворяют условию Гёлъдера в Нк, ад удовлетворяет условию Гёльдера в Q, кроме того с(-) < 0. Тогда функция г / * \ /г у г0) = Б / ехр / сЙ)Л + ^(60 ехР / , (0.4)

0 \о / \о /. где ^ - решение системы СДУ (0.3) с начальным условием £о = го является единственным дважды непрерывно дифференцируемым решением (0.1) в точке гд. □

Известно, что величине и(го) в случае однородного краевого условия Дирихле можно придать смысл осреднённой с весом д(-) полной концентрации диффузионных траекторий для источника, сосредоточенного в точке го, то есть и(го) = (#, Ц), где «о(-) - решение однородной краевой задачи для "прямого" диффузионного уравнения со свободным элементом /о (г) = 5 (г — г о). Фактически при этом используется двойственное представление функционала (д,иц). В диссертационной работе не только прослежена связь диффузионных процессов с вероятностным представлением решений краевых задач, но и построен весовой алгоритм, позволяющий, кроме оценки решения исходной дифференциальной задачи, одновременно оценивать его градиент.

Для приближенного статистического моделирования решения системы СДУ (0.3) можно, в частности, использовать простой в реализации и физически наглядный дискретный метод Эйлера с постоянным шагом г [28]: гп+1 = гп + ъ(гп)т + а(гп)г)пУ/т. (0.5) где гп - численная оценка решения (0.3) в узлах равномерной сетки по времени {тг ■ г}, а {г)п} - последовательность независимых между собой случайных векторов с независимыми стандартными гауссов-скими компонентами. Однако из-за невысокой скорости сходимости метода (0.5) приходится использовать методы более высоких порядков с более сложной реализацией (см., например, [38, 28]). В диссертации предложен новый метод повышения порядка сходимости дискретного метода Эйлера. Он основан на одновременной оценке искомого функционала методом (0.5) с двумя разными шагами по времени и последующей экстраполяции и является простым и наглядным.

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

S (С) = v< • £(С), где t(Q - среднее время, затрачиваемое на моделирование одного выборочного значения случайной величины £ (например, для процесса "блуждания по сферам" величина t(Q определяется средним числом переходов в цепи до момента попадания в е - окресность границы Г, причём для широкого класса границ известна логарифмическая оценка [17, 18]: t(Q < Сг| 1п(е)|). Поскольку общая погрешность оценки, полученной методом Монте-Карло, обычно состоит из двух частей: детерминированной и вероятностной, целесообразно выбирать параметры алгоритма таким образом, чтобы обе эти погрешности имели один и тот же порядок. Вероятностная погрешность возникает в связи с заменой математического ожидания арифметическим

1 N средним выборочных значений — ^ Cj и пропорциональна величине

N j=i л/vr где N - количество моделируемых выборочных значений. Де-v N терминированная погрешность появляется, например, при введении е - окресности границы (в алгоритмах "блуждания по сферам") или при замене диффузионного процесса кусочно-линейным процессом с шагом по времени т (в методе Эйлера).

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

В диссертационной работе предполагаются выполненными все условия регулярности данных задачи (коэффициенты, граничная поверхность области, краевые условия и т.п.) для существования и единственности некоторого обобщенного решения задачи (0.1), а также его вероятностного и интегрального представление с помощью шаровой функции Грина. При этом условия теорем 0.2 и 0.1 могут видоизменяться, в частности, вместо условия с(г) < 0 полагаем с(г) < с*, где —с* — первое собственное число оператора Ь в области а вместо условия (0.2) считаем матрицу В(г) = {Ьц(г))^=1 к равномерно положительно определенной в О (см., например, [20, 34]).

Далее следует краткое содержание диссертационной работы по главам и параграфам.

Заключение

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

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

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

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

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

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

6. Осуществлен ряд численных расчетов, подтверждающих теоретические результаты.

1. Михайлов Г.А., Бурмистров А.В. Оценка "по времени" в методе Монте - Карло / / Доклады Академии Наук. 1999, том 367, № 1, с.7 - 10.

2. А. V. Burmistrov and G.A.Mikhailov. "Path" and "Time" estimates in the Monte Carlo method // Russian Journal of Numerical Analysis and Mathematical Modelling.1999, vol.14, № 3, pp.221 236.

3. Burmistrov A.V. "Time" estimate in the Monte Carlo method.Bull. NCC. 1999. -Spec. Iss. - pp.25 - 34.

4. Михайлов Г.А., Бурмистров А.В. Решение стационарного диффузионного уравнения методом Монте Карло с вычислением производных / / Доклады Академии Наук. 2000, том 372, № 6, с.736 - 739.

5. Михайлов Г.А., Бурмистров А.В. Новый метод Монте Каро для решения стационарного диффузионного уравнения // Сибирский математический журнал.2000, том 41. № 5, с.1098 1105.

6. Burmistrov А. V. The Monte Carlo method for conjugate stationary diffusion equation with special item // Bull. NCC. Ser. Num. Anal. 2001. - Iss. 10. - pp.11 -18.

7. Бурмистров А.В. Численное решение однородного стационарного диффузионного уравнения новым методом Монте-Карло / / Труды конференции молодых ученых, Новосибирск, 2001, с.29 36.

8. А. V.Burmistrov. Numerical solution of homogeneous stationary diffusion equation by new Monte Carlo method // Proceedings of the 4th St.Petersburg Workshop on Simulations. N11 Chemistry St.Petersburg University Publishers, 2001, pp.172 177.

9. G.A.Mikhailov, A.V.Burmistrov. The Monte Carlo method for conjugate stationary diffusion equation // Proceedings of the 4th St.Petersburg Workshop on Simulations. Nil Chemistry St.Petersburg University Publishers, 2001, pp.343 348.

10. Бурмистров А.В. Оценка производных от решения стационарного диффузионного уравнения методом Монте-Карло // Труды конференции молодых ученых, Новосибирск, 2002, с.ЗЗ 40.

11. A.V.Burmistrov and G.A.Mikhailov. Estimation of the gradient of the solution of an adjoint diffusion equation by the Monte Carlo method // Russian Journal of Numerical Analysis and Mathematical Modelling. 2002, vol.17, № 4, pp.367 380.

12. Михайлов Г.А., Бурмистров А.В. Весовой метод Монте-Карло для оценки производных от решения сопряженного диффузионного уравнения // Доклады Академии Наук. 2002, том 386, № 1, с. 18 21.

13. Браун Дж. Методы Монте-Карло // Современная математика для инженеров. М.: Изд-во иностр. лит., 1959. - С.275-301.

14. Вентцелъ А.Д. Курс теории случайных процессов. М.:Наука, 1975.

15. Данфорд Н., Шварц Дж.Т. Линейные операторы (общая теория.) М.: йзд-во иностр. лит., 1962.

16. Ермаков С.М., Михайлов Г.А. Статистическое моделирование. М, Наука, 1982.

17. Ермаков С.М., Некруткин В.В., Сипин А.С. Случайные процессы для решения класических уравнений математической физики. М, Наука, 1984.

18. Икрамов Х.Д. Численное решение матричных уравнений. М.: Наука, 1984.

19. Кошеле в А.И. Регулярность решений эллиптических уравнений и систем. М.: Наука, 1986.

20. Крылов И.В. Лекции по эллиптическим и параболическим уравнениям в пространствах Гёльдера. Новосибирск, Научная книга, 1998.

21. Макаров Р.Н. Комбинированные оценки метода Монте-Карло для решения третьей краевой задачи для уравнения параболического типа / / Russian Journal of Numerical Analysis and Mathematical Modelling. 2002, vol.17, № 6, pp.501 515.

22. Математический энциклопедический словарь, Москва, "Советская энциклопедия", 1988

23. Михайлов Г.А. Новые весовые методы Монте-Карло для решения уравнения Гельмгольца // Доклады Академии Наук. 1992, том 326, № 6, с.943 947.

24. Михайлов Г.А. Весовые методы Монте-Карло. Новосибирск. Изд-во СО РАН, 2000.

25. Михайлов Г.А., Макаров Р.Н. Решение краевых задач второго и третьего рода // Сиб. матем. журн. 1997. Т.38, № 3, С.603-614.

26. Соболь И.М. Численные методы Монте-Карло. М, Наука, 1973.

27. Artemiev S.S., Averina Т.A. Numerical Analysis of Systems of Ordinary and Stochastic Differential Equations. Utrecht, Tokyo:VSP, 1997.

28. Baldi P. Exact asymptotics for the probability of exit from a domain and applications to simulation // Ann. Prob. (1995) 23(4), pp.1644 1670.

29. Borodin A.N., Salminen P. Handbook of Brownian Motion. Facts and Formulae. Basel-Boston-Berlin: Birkhauser Verlag, 1996.(Русское издание:Бородин A.H., Сал-минен П. Справочник по броуновскому движению. Факты и формулы. Санкт-Петербург, 2000.)

30. Cheshkova A.F. "Walk on spheres" algorithms for solving Helmholtz equation in the n-dimensional Euclidean space // Bull. NCC. Ser. Num. Anal. 1993. - Iss. 4. - pp.7- 18.

31. Cvesielski Z., Tayler S.J. First passage times and sojorn times for Brownian motion in space and the exact Hausdorff measure of the sample path // Trans. Amer. Math. Soc., 1962, v.103, № 3, pp.434 450.

32. Dynkin E.B. Markov Processes. Academic Press, New York, 1965.

33. Freidlin M. Functional Integration and Partial Differential Equations. Princeton Univ. Press., Princeton, 1985 (Ann. of Math. Stud.,109)

34. Gobet E. Weak approximation of killed diffusion using Euler schemes // Stochastic Processes and their Applications 87 (2000), pp.167 197.

35. Hammersley I.M., Morton K.W. A new Monte Carlo technique antithetic variates // Proc. Cambrige Philos. Soc. 1956. v.52. pp.449 - 475.

36. Kushner H.J. Probability Methods for Approximation in Stochastic Control and Elliptic Equations. Academic Press, New York, 1991.

37. Milsht ein G.N. The solving of boundary value problems by numerical integration of stochastic equations // Mathematic and Computers in Simulation (1995) 38, pp.77 85.

38. S. Ulam and N. Metropolis. The Monte Carlo method, Journal of American Statistical Association, 44, 335, 1949.