CC BY
99
УДК 519.245 Вестник СПбГУ. Математика. Механика. Астрономия. 2019. Т. 6 (64). Вып. 3 МБС 65С05
Метод Монте-Карло для решения систем ОДУ*
С. М. Ермаков, Т. М. Товстик
Санкт-Петербургский государственный университет,
Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9
Для цитирования: Ермаков С. М., Товстик Т. М. Метод Монте-Карло для решения систем ОДУ // Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия. 2019. Т. 6(64). Вып. 3. С. 411-421. https://doi.org/10.21638/11701/spbu01.2019.306
Рассматривается применение метода Монте-Карло для решения задач Коши для систем линейных и нелинейных обыкновенных дифференциальных уравнений. Метод Монте-Карло актуален для решения больших систем уравнений и в случае малой гладкости исходных функций. Система приводится к эквивалентной системе интегральных уравнений Вольтерра. Для линейных систем это преобразование позволяет снять ограничения, связанные со сходимостью мажорантного процесса. Приводятся примеры оценок функционалов от решения и обсуждается поведение их дисперсий. В общем случае промежуток интегрирования разбивается на конечные интервалы, на которых нелинейная функция аппроксимируется полиномом. Полученное интегральное уравнение решается с помощью ветвящихся цепей Маркова с поглощением. Обсуждаются возникающие при этом задачи распараллеливания алгоритмов. В качестве примера рассматривается одномерное уравнение с кубической нелинейностью. Обсуждается выбор переходных плотностей ветвящегося процесса. Подробно описывается метод поколений. Дается сравнение численных результатов с решением, полученным методом Рунге — Кутта.
Ключевые слова: метод Монте-Карло, системы ОДУ, интегральные уравнения, статистическое моделирование.
1. Введение. Теория методов Монте-Карло для решения линейных краевых задач математической физики в настоящее время достаточно хорошо развита [1, 2]. Обычно задачу сводят к эквивалентному интегральному уравнению. Как известно, интегральные уравнения 2-го рода тесно связаны с цепями Маркова и наиболее распространенные алгоритмы основаны на моделировании N независимых траекторий цепи Маркова, связанной с соответствующим интегральным уравнением условием согласования. При этом предполагается существование итерационного решения мажорантного интегрального уравнения. Последнее условие может быть довольно ограничительным, и имеется ряд работ [3, 4], где предложены специальные методы для снятия этих ограничений. Для нелинейных уравнений используются различные методы последовательной аппроксимации, среди которых особое место занимает приближенная замена нелинейности общего вида полиномиальной нелинейностью, что позволяет строить алгоритмы с естественным параллелизмом.
Что касается задач с начальными условиями, то алгоритмы методов Монте-Карло для них мало исследованы даже в случае ОДУ. Можно сослаться лишь на работу [5], где имеются примеры решения методом Монте-Карло систем ОДУ. В то же время имеются важные прикладные задачи, требующие решения больших систем
* Работа выполнена при финансовой поддержке РФФИ (грант № 17-01-00267-а). (¡5 Санкт-Петербургский государственный университет, 2019
такого рода или систем определенного вида (например, стохастических дифференциальных уравнений). В указанных случаях метод Монте-Карло может быть более эффективным, чем классические вычислительные методы.
Начальная краевая задача для систем ОДУ приводится к системам интегральных уравнений типа Вольтерра. Это позволяет для линейных систем снять ограничения, связанные со сходимостью мажорантного процесса теоретически для всех времен t G [0, то). Относительно некоторых классов нелинейных систем также возможен подход, основанный на приближении полиномами, но теорема Пикара диктует ограничения на длину промежутка, где существует итерационное решение, и метод Монте-Карло работает без синхронизации только на этом промежутке. При этом возникает стохастическая погрешность, которая может накапливаться в процессе вычислений, поэтому необходим анализ ее поведения.
2. Метод, основанный на цепях Маркова. Наиболее распространенным алгоритмом метода Монте-Карло является алгоритм, основанный на моделировании цепей Маркова.
Пусть / — вероятностная мера, и требуется решить интегральное уравнение
tp(x) = J k(x,y)tp(y) /(dy) + f (x), (mod/л), (1)
где k и f — заданные функции, (mod/) обозначает, что равенство выполняется на носителе меры
Предположим, что сходится последовательность (фп, h) при n ^ то для всех h из некоторого множества H функций, где Фрп, n = 0,1,..., — последовательность функций, определяемых равенством
<рп(х) = J Ik(x,y)l фп-1 (y) /(dy) + lf (x)\), (mod/), ф0 = |f 1 (2)
и ((рп, h) = J (ftn (x) h(x) /(dx). Алгоритм (2) приводит к решению уравнения
<p(x) = J \k(x,y)l<f(y) /(dy) + \f(x)1, (mod/), (3)
которое является мажорантным по отношению к (1). Алгоритм (2) сходится, если сходится ряд Неймана, и в этом случае выполняется равенство
то „
(p(x) = If(x)| / |k(x,xi)||k(xi,x2)| • •• |k(xm-i,xm)||f(xm)| <g>™i /(dxi). (4)
m=1
Выберем цепь Маркова, определяемую начальной плотностью распределения p0(x), переходной плотностью p(x,y) (по отношению к мере /) и вероятностью поглощения g(x) так, что выполняются равенства
jp(x, y) /(dy) = 1 — g(x), (mod/), 0 < g(x) < 1, (5)
и условия согласования
p0(x) > 0, если h(x) = 0,
p(x,y) > 0, если k(x,y)=0, (6)
g(x) > 0, если f(x) = 0.
Моделируя с помощью известных алгоритмов траектории этой цепи
х0 ^ х1 ^ ■ ■ ■ ^ хт_1 ^ хт (7)
в предположении, что случайный момент поглощения т конечен, с вероятностью 1 можем оценить (вычислить) величину (к, у), так как справедливо следующее утверждение:
17 7 /и \ Т к0к0,1 •••кт _1,т /т
ЕЛ = (/1,у>), Л =-п—:-:-• (8)
р0ро,1 • • .Рт _1,т дт
Здесь ко = к(хо), к^+1 = к(х^, х^+1), г = 0,1, • • •, т — 1, и аналогичные обозначения имеют место для /т, р0, Рг,г+1, дт. Случайную величину Ja называют оценкой «по поглощениям».
Основанные на методе Монте-Карло оценки функционала (к, у) используют представление решения уравнения (1) в виде ряда Неймана и предполагают его абсолютную сходимость, т.е. сходимость ряда (4).
Моделируем N независимых траекторий, для каждой из которых вычисляются значения Ja(1) (1 — номер траектории) и (к, у) оценивается с помощью среднего 1/NJ2Ja (1). Соответствующие доказательства и модификации алгоритма содержатся в [4]. Мы привели самую простую из известных оценок Монте-Карло для (к, у). Для нас важно отметить параллелизм алгоритма.
Известен ряд алгоритмов Монте-Карло для решения некоторых типов нелинейных интегральных уравнений [6]. Как правило, они основаны на последовательной приближенной линеаризации задачи. Это обычно ухудшает свойства параллелизма, так как на каждом этапе линеаризации требуется синхронизация. Исключением являются уравнения с полиномиальной нелинейностью (уравнения Ляпунова — Шмидта) вида
то „
у(х)=Е / К;(х, х1, • • • ,х1)П1т=1У(хт) ®т=1 ) + / (х) • (9)
г=1
При соответствующих условиях (см. [7,8]) функционал (к, у) от решения этого уравнения можно оценить, моделируя ветвящийся марковский процесс. На его траекториях вычисляется аналог оценки Ja и процедура не требует дополнительной синхронизации. Практический интерес представляют сравнительно небольшие степени нелинейности.
Перечисленные выше вычислительные методы и их модификации представляют интерес при решении многомерных задач, если функция у(х) зависит от большого числа переменных и (или) требуется решить большие системы уравнений.
Отметим, что важные в прикладном отношении краевые задачи для уравнений в частных производных сводятся к задачам, перечисленным выше. В этой области имеется обширная литература.
Вместе с тем стохастические методы для решения задачи Коши для больших систем обыкновенных дифференциальных уравнений исследованы мало. Первые результаты в этом направлении получены в [5], где намечен общий подход, но не обсуждалась сравнительная эффективность метода и его параллелизм. Задача Коши для ОДУ, однако, обладает своеобразием, которое можно уже проследить на примере задачи Коши для системы линейных уравнений.
3. Системы линейных дифференциальных уравнений. Пусть г € [0, то), у (г) = (у1(г),уп(г))т, ^ (г) = (¡1(г),и(г))т и
¿У
— =Аа)¥ + РЦ1 АЦ) = \Мтъ=1, У(0)=Уо. (10)
Эквивалентной системой интегральных уравнений будет
г г
У (г)=I А(т)У (т) ¿т ^(т) ат+Уо. (11)
о о
У (г) может быть найдено с помощью метода последовательных приближений, что эквивалентно вычислению суммы ряда
ТО г Т1 Тк — 1
У (г)=Р(г)+^ I А(т1 )вт-1! А(т2) ¿т2... I А(тк )Р(тк) ¿тк,
г
где Р(г)= ^(т) ат+Уо. (12)
о
В рассматриваемом случае решение задачи (10) может быть записано в виде интегральной экспоненты. Ее вычисления также приведут нас к ряду, аналогичному (12). Ниже об этом сказано подробнее.
Ряд (12) абсолютно сходится при любом г и для его вычисления методом Монте-Карло можно использовать подход, аналогичный описанному в п. 2, т.е. построить несмещенную оценку решения на траекториях цепи Маркова. Случайные дискретные моменты времени г к будем предполагать распределенными с плотностями г к (т), т€[0,гк-1], к=1, 2,... Для наших целей достаточно ограничиться равномерным распределением гк (т )=1/гк-1.
Введем в рассмотрение цепь Маркова (рР°(г), V(г)) с начальной плотностью распределения р°(г), ^П=1 р°(г)=1 и переходной матрицей V(г)=\\рг^(г)\\П]=1, Е;=1 Рг(г)=1-дг(г), г=1,...,п, связанную с матрицей А(г) условиями
согласования вида (6). Определим для функционала (Н, У)= Н(т)У (т) ¿т оценку
у = Ы0 (¿о)а»0,н (¿1) • • •а»т_1,»т(^т)-Р1(^т) -гт^ , ,
~ 7=0
Здесь гк при к>0 — случайные величины, равномерно распределенные на отрезках (0,гк-1), а г_1=г. Случайные величины гк имеют дискретные распределения [чгк-1,1(гк),..., Чгк-1 ,п(гк)], Чг,] (г)=Ръ,з (г)/(1-дг(г), к=1, 2,... Последнее произведение — результат нормировки равномерных распределений.
Имеет место следующее утверждение.
Утверждение. При сделанных выше предположениях случайная величина 7 является несмещенной оценкой функционала (Н,У).
Утверждение легко доказывается непосредственными вычислениями.
Нетрудно также получить формальное выражение для дисперсии J и даже записать систему дифференциальных уравнений, через решение которой дисперсия выражается. Однако задача минимизации дисперсии оказывается здесь существенно более сложной, чем в случае уравнений Фредгольма, и мы рассмотрим самый простой пример. Обсудим вопрос о выборе цепи Маркова, обращающей в нуль дисперсию оценки.
Для решения уравнения х = ах+6 при |а| < 1 имеем «цепь Маркова» р+д=1, р>0, д>0. Случайное блуждание: с вероятностью р продолжить процесс, с вероятностью </ — оборвать. Оценка (8) имеет вид £1 = где г — номер последнего перехода перед поглощением. Легко получить Е^1=х при а>0, р=а, д=1—а, дисперсия оценки обращается в нуль.
Для уравнения у'=ау, у(0)=6 имеем, применяя метод последовательных приближений,
* то * ткг1
Уй+1^)= / (т) ¿т + 6, у(£)=6+^ / а^Т1 ■■■ а^тдД
0 Й=10 0
Используя траектории той же «цепи Маркова», что и в предыдущем случае,
т
можем построить несмещенную оценку для у(£) в виде £2 = ат П а^6/(ртд), где а^,
к =1, 2, • • •, равномерно распределены на промежутке от 0 до а&_1, а а0 =
Теорема о существенной выборке [4, теорема 3.1] указывает меру, минимизирующую дисперсию, но она не обязана быть мерой, индуцируемой однородным марковским процессом. Задача ее моделирования в общем случае — это отдельная сложная задача.
В более простом случае, если речь идет о вычислении суммы ряда еа= Х^о Ж' дисперсию оценки £з=—г—г при а>0 с помощью выбора р обратить в нуль нельзя, но, выбирая вместо геометрического распределения траекторий {ртд} распределение Пуассона рт=^е~а, мы очевидным образом обратим дисперсию в нуль. Это замечание указывает на существенные отличия, возникающие при анализе дисперсий оценок решений уравнений фредгольмовского и вольтерровского типов.
4. Общий случай и дополнительные замечания. В предыдущих пунктах мы провели параллель между методами Монте-Карло для решения уравнений типа Фредгольма и методами решения задачи Коши для систем линейных ОДУ. Можно утверждать, что похожие результаты получаются при решении следующих задач:
1) построение двойственной сопряженной к (13) оценки;
2) построение аналога оценки по поглощению для уравнений с полиномиальной нелинейностью.
Последняя проблема подробно обсуждается ниже при решении численного примера.
Таким образом, для решения задачи Коши для систем обыкновенных дифференциальных уравнений (ОДУ) можно предложить следующую общую схему приближенного алгоритма.
Задача
Х=(х1,...,ха), /=(/!,..., /Д Х(Ьо)=Хо=(х1,---,х°), (14)
аппроксимируется задачей
где
Р (¿,Х )=(Р (¿,Х Д^Р^х)) и Р^,Х )= (¿)
¿1 +...+¿5 = 1
^ 1 ¿8 с7п )х1 • • • х3 •
Функции с^1 подбираются из условия минимизации ||/—Р|| — некоторой выбранной исследователем нормы разности /—Р на промежутке (¿о —Н,¿о+Н). Если в линейном случае при выполнении условий теоремы Пикара [9] о существовании решения системы (14) Н может быть практически неограниченным, то в случае нелинейном — и условия теоремы Пикара, и требование малости нормы /—Р могут налагать на Н довольно жесткие условия. Теоретические соображения и, возможно, пробные вычисления должны дать возможность выбрать подходящие т и Д^<Н для вычисления решения системы в точке ¿о+Д^.
Если ставится задача оценки X(Т) при Т>(£о+ДД то промежуток (¿о, ¿о+Т) делится на части, имеющие длину порядка Д^: и решаются последовательно задачи с начальными данными Хо,Х(¿о+Д:),Х(¿о+2Д:), • • • Волна обозначает, что соответствующее значение X имеет систематическую ошибку, вызванную полиномиальной аппроксимацией и случайную ошибку метода Монте-Карло. Анализ поведения погрешности в самом общем случае вряд ли может быть содержательным. Его следует проводить для конкретных классов уравнений. В детерминированном случае имеются фундаментальные результаты относительно существования решения и устойчивости по отношению к погрешностям. Если устойчивость имеет место, то при достаточно большом числе испытаний можно гарантировать близость с заданной вероятностью оценки метода Монте-Карло к точному решению.
Следующее замечание относительно параллелизма алгоритма должно быть сделано.
В пределах одного промежутка длины Д: алгоритм метода Монте-Карло состоит в моделировании независимых (может быть ветвящихся) траекторий и обладает неограниченным параллелизмом. Для перехода на следующий интервал требуется синхронизация. Это соображение также может быть важным при выборе Д:, т и точности полиномиальной аппроксимации.
Подробному анализу поведения погрешности метода и повышению его эффективности для определенных классов задач авторы предполагают посвятить дальнейшие публикации. Наиболее перспективным представляется применение его к решению стохастических дифференциальных уравнений, где он может служить источником новых алгоритмов.
Ниже приводится иллюстративный пример.
5. Пример. Пусть
х = вхп — ах, х(:о) = хо, а > 0, ¿о < t < Т (15)
Это уравнение при Т = £ эквивалентно интегральному
(¿)= в е_а(4_т Цт )п ¿т + e_а(í_ío)xo• (16)
Ло
Нелинейное интегральное уравнение (16) запишем в виде
х(*)=/' к(^т)х(т)п¿г + /(¿), к(^т)= ве-а('-т/(¿) = е-а(4-4о)жо. (17) Ло
Для вычисления решения уравнения (17) методом Монте-Карло, согласно изложенным выше алгоритмам, нужно использовать ветвящийся процесс, у которого в каждой точке ветвления рождается ровно п новых частиц. Весь интервал (0, Т) интегрирования нужно разбить на малые подинтервалы и искать решение на них. Будем вычислять
), = кд^ к = 1,..., £ = [т/дг]. (18)
Для нормы ядра к^^), ¿о < ¿2 < < г, должно выполняться неравенство
п||А;(*1,*2)||= вир = (19)
*о<*1 <еЛо |а|
из которого следует найти ограничения на параметры в, а, ¿,¿0.
Воспользуемся методом «поколений». В рассматриваемой задаче в начальный момент времени в точке £ имеется одна частица, которая может погибнуть, и тогда алгоритм заканчивается. Если частица не гибнет, то она превращается в п новых частиц, которые образуют первое поколение.
Далее каждая частица независимо от других ведет себя так же, как исходная частица. Выжившие частицы образуют следующее поколение. Алгоритм заканчивается, если на очередное поколение не переходит ни одной частицы.
Каждая частица, появившись в новом поколении в момент в, может либо погибнуть с вероятностью д(в), либо превратиться в п новых частиц в случайные моменты времени
т8(к), к =1,...,п, (20)
которые одинаково распределены и независимы. Плотность распределения г(в, ¿2), ¿о < ¿2 < в, этих случайных величин, которая называется условной плотностью
(к)
перехода частицы из точки в в точку т8 , удовлетворяет равенству
[ г(8,*2)^2 = 1.
Ло
По
Функции $(в) и г (в, ¿2) нужно задать и найти переходную плотность, которая равна
р(*1,*2) = (1 - з(*1)) г(^ 1, ¿2). (21)
Моделирующие формулы случайных величин тя(к), к = 1,...,п, при любом к получаются из уравнений
/ г(в, ¿2)^2 = а(к), (22)
в которых а(к) — независимые случайные величины, равномерно распределенные на (0,1).
На случайных моментах (20), в которых происходит изменение состава частиц, строится весовая функция
Посмотрим, какой вклад в оценку вносит исходная частица. Она с вероятностью д(г) гибнет, тогда алгоритм заканчивается. Если частица не гибнет, то появляются п новых частиц. Вклад начальной частицы таков:
, —с вероятностью шг),
= ^ д(г)
Пп=1 С(к), С(к) = №(г,т(к)) с вероятностью 1 - д(г), (к)
здесь тг — независимые случайные величины, которые моделируются с помощью формулы (22).
Пусть {г[к\...
, гПк} — координаты множества частиц Т(к), образующих к-е поколение. В Т(к) = С(к) и Н(к) входят частицы С(к), из которых образуется к + 1 поколение и частицы
Н(к), которые гибнут. Частицы к-го поколения вносят с соответствующими вероятностями следующий вклад в оценку:
п (к)
Л= п п Ш
г<м)сСцк) г=1 3 г(к)си(к) д(г1 )
а оценка, учитывающая состояние системы с нулевого по к-е поколения, равна Пк=о ^г. Алгоритм заканчивается на К поколении, если множество 0(К+1 оказывается пустым, а 0(К) — нет. В результате, одна реализация, оценивающая решение уравнения (17), принимает вид
К
(1)
V(1) =Ц Л, (23)
к=о
здесь К(1) — число поколений в первой реализации. Этот алгоритм нужно проделать N раз и при больших значених N получаем
N
х(г) « x(г)N = £ V(3)/к
3 = 1
Условие (19) приводит к гибели ветвящегося процесса с вероятностью единица и, следовательно, к конечности числа поколений
К(3)
при любом ].
Рассмотрим конкретную задачу, которая имеет вид (15),
Х = 0.4х3 - 0.6х, х(0) = 0.7, в = 0.4, а = 0.6, 0 < г < 4. (24)
Учитывая вид ядра к(г1,г2) = 0.4е о б(г1 г2), выбирае
м
Ье_Ь(г1_г2)
Г(*ЬЪ)= 1_е-ь(г1-г0)> Ь-а' ^о < ¿2 < ¿1-
В частности, можно выбрать Ь = а. Так как переходная плотность удовлетворяет равенству (21), то в качестве функции д(в) — вероятности поглощения в точке в, берем д(в) = е_с(в_го). Заметим, что с ростом номера поколения вероятности поглощения в нем возрастают.
Таблица 1. Таблица 2.
А4 = 0.25, Ъ = 0.6
№ с 5 ^тах К X
1 4.0 0.368 0.035 и 3 1
2 6.6 0.192 0.048 13 4 0
3 9.0 0.105 0.097 17 5 6
А4 = 0.15, 6 = 0.5
№ с 5 ^тах К X
1 2.5 0.678 0.020 6 1 3
2 4.0 0.549 0.017 8 2 0
3 6.5 0.377 0.036 13 2 0
Для сходимости алгоритма достаточно выполнения условия (^ 1, ¿2) | < 1, которое приводит к следующему неравенству:
|Д|р-(а-Ь)(41 -«2)(1 р-Ь(«1-М)
«I- " т-Л-.))-14 >• (25>
Неравенство (25) накладывает условия: а > 6, с > 6.
Нужно задать параметры 6 и с. Основываясь на идее существенной выборки при задании величины параметра 6, можно взять 6 = а или число, близкое к а. Параметр с определяет поглощение в начальный момент времени и, так как с этой вероятностью алгоритм не дойдет до первого поколения, следовательно, д(2) не должно быть слишком большим.
С числом испытаний N = 1000 были промоделированы два варианта при шаге Д4 = 0.15 и Д4 = 0.25. Для проверки точности были вычислены выборочные среднеквадратические отклонения оценок ж^)^ и найдены величины
ь
X = У^ I(И^) - | > 2<Г^), 5 = тах |ж(^) - |. (26)
к=1 _ _
В каждом из вариантов найдены величины максимального Ктах и среднего К (округленная величина) чисел поколений, которые проявились при вычислении оценок (18).
Исходные данные для трех вариантов параметров и результаты моделирования приведены, соответственно, при Д^; = 0.25 в табл. 1 (рис. 1) и при Д^; = 0.15 в табл. 2 (рис. 2). На каждом из рисунков построены графики ж^в виде синей, зеленой и фиолетовой линий, которые соответствуют данным с номерами 1-3 в таблицах. На каждом графике красной линией отмечено решение, полученное методом Рунге — Кутта.
Однако не при всех значениях параметра с алгоритм приводит к хорошему приближению. Чтобы устранить это препятствие, можно, используя описанный метод, найти оценки xn при разных значениях c и затем, усреднив их, получить xn .
В нашем примере, чтобы получить Xn, найдем и усредним приближения при со = 2, Ci = со + i1.5,1 < i < 4. При At = 0.15 величина S, полученная по формуле (26) с заменой XN на XN, дает S(b = 0.6) = 0.056 и S(b = 0.5) = 0.024, а при At = 0/25 получено S(b = 0.6) = 0.056 и S(b = 0.5) = 0.018. Отсюда и из анализа таблиц делаем вывод, что лучшая оценка получается, если взять b = 0.5.
Литература
1. Ермаков С. М., Сипин А. С. Метод Монте-Карло и параметрическая разделимость алгоритмов. СПб.: Изд-во С.-Петерб. ун-та, 2014.
2. Михайлов Г. А., Войтишек А. В. Численное статистическое моделирование. Методы Монте-Карло. М.: Академия, 2006.
3. Ermakov S.M., Wagner W. Monte Carlo difference schemes for the wave equation // Monte Carlo Methods and Appl. 2002. Vol.8, no. 1. P. 1-30.
4. Ермаков С. М. Метод Монте-Карло //В кн. Метод Монте-Карло в вычислительной математике (вводный курс). СПб.: Невский Диалект, Бином. Лаборатория знаний, 2009.
5. Akhtar M.N., Durad M.H., Ahmed A. Solving initial value ordinary differential equations by Monte Carlo method // Proc. of IAM. 2015. Vol.4, no. 2. P. 149-174.
6. Halton J. H. Sequential Monte Carlo techniques for solving non-linear systems // Monte Carlo Methods and Applications. 2006. Vol. 12, no. 2. P. 113-141.
7. Ермаков С. М. Об аналоге схемы Неймана — Улама в нелинейном случае // ЖВМ и МФ. 1973. Vol. 13, no. 3. P. 564-573.
8. Ермаков С. М. Метод Монте-Карло и смежные вопросы. М.: Наука, 1975.
9. Понтрягин Л. С. Обыкновенные дифференциальные уравнения. М.: Наука, 1974.
Статья поступила в редакцию 29 декабря 2018 г.;
после доработки 27 февраля 2019 г.; рекомендована в печать 21 марта 2019 г.
Контактная информация:
Ермаков Сергей Михайлович — д-р физ.-мат. наук, проф.; [email protected] Товстик Татьяна Михайловна — канд. физ.-мат. наук, доц.; [email protected]
Monte Carlo method for solution of systems ODE*
S. M. Ermakov, T. M. Tovstik
St. Petersburg State University, Universitetskaya nab., 7-9, St. Petersburg, 199034, Russian Federation
For citation: Ermakov S. M., Tovstik T. M. Monte Carlo method for solution of systems ODE. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy, 2019, vol. 6(64), issue 3, pp. 411-421. https://doi.org/10.21638/11701/spbu01.2019.306 (In Russian)
An application of the Monte Carlo method to solution of problems Cauchy for system of linear and nonlinear ordinary differential equations is considered. the Monte Carlo method is topical for solution of large systems of equations and in the cases with given functions are not smooth enough. A system is reduced to an equivalent system of integral equations of the Volterra type. For linear systems this transformation allows to avoid restrictions, connected with a convergence of a majorizing process. Examples of estimates of functionals of solutions
*The work is supported by Russian Foundation for Basic Research (grand N 17-01-00267-а).
are given, and a behavior of their dispersions are discussed. In the general case an interval of solution is divided into finite subintervals in that nonlinear functions are approximated by polynomials. An obtained integral equations are solved by using branching Markov chains with an absorption. Appearing problems of parallel algorithms are discussed. As an example an one dimensional cubic equation is considered. A method of levels of branching process is discussed in details. a comparison of numerical results with a solution, obtained by the Runge — Kutta method, is presented.
Keywords: Monte Carlo method, systems ODE, integral equations, statistic simulation. References
1. Ermakov S.M., Sipin A.S., Monte Carlo method and parametric separation of algorithms (St. Petersburg Univ. Press, St. Petersburg, 2014). (In Russian)
2. Mikhailov G.A., Voytishek A. V., Numerical statistical modeling. Monte Carlo methods (Akademiya Publ., Moscow, 2006). (In Russian)
3. Ermakov S. M., Wagner W., "Monte Carlo difference schemes for the wave equation", Monte-Carlo Methods and Appl. 8(1), 1-30 (2002).
4. Ermakov S.M., Monte Carlo method in quantitative mathematics (Introduction) (Nevskii Dialekt, Binom. Laboratoria znaniy Publ., St. Petersburg, 2009). (In Russian)
5. Akhtar M.N., Durad M.H., Ahmed A., "Solving initial value ordinary differential equations by Monte Carlo method", Proceeding of IAM 4(2), 149-174 (2015).
6. Halton J.H., "Sequential Monte Carlo techniques for solving non-linear systems", Monte Carlo Methods and Applications 12(2), 113-141 (2006).
7. Ermakov S.M., "On analogiue of Neuman-Ulam scheme in the nonlinear case", Computational Mathematics and Mathematical Physics 13(3), 564-573 (1973). (In Russian)
8. Ermakov S. M., Monte Carlo method and related topics (Nauka Publ., Moscow, 1975). (In Russian)
9. Pontryagin L. S., Ordinary Differential Equations (Addison-Wesley, 1962).
Received: December 29, 2018 Revised: February 27, 2019 Accepted: March 21, 2019
Author's information:
Sergey M. Ermakov —[email protected] Tatiana M. Tovstik — [email protected]
