Научная статья на тему 'О МЕТОДЕ МОНТЕ-КАРЛО ДЛЯ РЕШЕНИЯ БОЛЬШИХ СИСТЕМ ЛИНЕЙНЫХ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ'

О МЕТОДЕ МОНТЕ-КАРЛО ДЛЯ РЕШЕНИЯ БОЛЬШИХ СИСТЕМ ЛИНЕЙНЫХ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ Текст научной статьи по специальности «Математика»

CC BY
163
28
Читать
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
МЕТОД МОНТЕ-КАРЛО / СИСТЕМЫ ОДУ / ИНТЕГРАЛЬНОЕ УРАВНЕНИЕ / ЗАДАЧИ МАССОВОГО ОБСЛУЖИВАНИЯ / ОПТИМАЛЬНАЯ ПЛОТНОСТЬ / НЕСМЕЩЕННАЯ ОЦЕНКА / СТАТИСТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

Аннотация научной статьи по математике, автор научной работы — Ермаков Сергей Михайлович, Смиловицкий Максим Григорьевич

В статье рассматривается применение метода Монте-Карло к решению задачи Коши для больших систем линейных дифференциальных уравнений. В первой части статьи дается краткий обзор уже известных результатов применения метода для решения интегральных уравнений Фредгольма. В основной части статьи разбирается применение подхода к системе линейных ОДУ, которая приводится к эквивиалентной системе интегральных уравнений Вольтерра. Это позволяет снять ограничения, связанные со сходимостью мажорантного процесса. Формулируются следующие ключевые теоремы. Теорема 1 указывает требуемые условия согласования, которым должны отвечать переходная и начальная плотности распределения, инициирующие соответствующую цепь Маркова, для которой выполняется равенство между математическим ожиданием оценки и интересующим нас функционалом. Теорема 2 формулирует выражение для дисперсии оценки, в то время как теорема 3 указывает параметры цепи Маркова, минимизирующие значение дисперсии для оценки функционала. В работе приводятся доказательства всех трех теорем. В практической части предложенный метод применяется к системе линейных ОДУ, описывающих замкнутую систему массового обслуживания из десяти условных машин и семи условных рабочих. Решение приводится как для системы с постоянной матрицей коэффициентов, так и для системы с переменной матрицей, где в зависимости от времени меняется интенсивноcть выхода машин из строя. Также произведено сравнение решения методом Монте-Карло с решением методом Рунге - Кутта. Все результаты отражены в таблицах.

i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
Предварительный просмотр
i Надоели баннеры? Вы всегда можете отключить рекламу.

MONTE-CARLO FOR SOLVING LARGE LINEAR SYSTEMS OF ORDINARY DIffERENTIAL EQUATIONS

Monte-Carlo approach towards solving Cauchy problem for large systems of linear differential equations is being proposed in this paper. Firstly, a quick overlook of previously obtained results from applying the approach towards Fredholm-type integral equations is being made. In the main part of the paper, a similar method is being applied towards a linear system of ODE. It is transformed into an equivalent system of Volterra-type integral equations, which relaxes certain limitations being present due to necessary conditions for convergence of majorant series. The following theorems are being stated. Theorem 1 provides necessary compliance conditions that need to be imposed upon initial and transition distributions of a required Markov chain, for which an equality between estimate’s expectation and a desirable vector product would hold. Theorem 2 formulates an equation that governs estimate’s variance, while theorem 3 states a form for Markov chain parameters that minimise the variance. Proofs are given, following the statements. A system of linear ODEs that describe a closed queue made up of ten virtual machines and seven virtual service hubs is then solved using the proposed approach. Solutions are being obtained both for a system with constant coefficients and time-variable coefficients, where breakdown intensity is dependent on t. Comparison is being made between Monte-Carlo and Rungge Kutta obtained solutions. The results can be found in corresponding tables.

Текст научной работы на тему «О МЕТОДЕ МОНТЕ-КАРЛО ДЛЯ РЕШЕНИЯ БОЛЬШИХ СИСТЕМ ЛИНЕЙНЫХ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ»

УДК 519.245 Вестник СПбГУ. Математика. Механика. Астрономия. 2021. Т. 8 (66). Вып. 1 МБС 65С05

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

С. М. Ермаков, М. Г. Смиловицкий

Санкт-Петербургский государственный университет,

Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9

Для цитирования: Ермаков С. М., Смиловицкий М. Г. О методе Монте-Карло для решения больших систем линейных обыкновенных дифференциальных уравнений // Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия. 2021. Т. 8(66). Вып. 1. С. 37-48. https://doi.org/10.21638/spbu01.2021.104

В статье рассматривается применение метода Монте-Карло к решению задачи Коши для больших систем линейных дифференциальных уравнений. В первой части статьи дается краткий обзор уже известных результатов применения метода для решения интегральных уравнений Фредгольма. В основной части статьи разбирается применение подхода к системе линейных ОДУ, которая приводится к эквивиалентной системе интегральных уравнений Вольтерра. Это позволяет снять ограничения, связанные со сходимостью мажорантного процесса. Формулируются следующие ключевые теоремы. Теорема 1 указывает требуемые условия согласования, которым должны отвечать переходная и начальная плотности распределения, инициирующие соответствующую цепь Маркова, для которой выполняется равенство между математическим ожиданием оценки и интересующим нас функционалом. Теорема 2 формулирует выражение для дисперсии оценки, в то время как теорема 3 указывает параметры цепи Маркова, минимизирующие значение дисперсии для оценки функционала. В работе приводятся доказательства всех трех теорем. В практической части предложенный метод применяется к системе линейных ОДУ, описывающих замкнутую систему массового обслуживания из десяти условных машин и семи условных рабочих. Решение приводится как для системы с постоянной матрицей коэффициентов, так и для системы с переменной матрицей, где в зависимости от времени меняется интенсивность выхода машин из строя. Также произведено сравнение решения методом Монте-Карло с решением методом Рунге — Кутта. Все результаты отражены в таблицах. Ключевые слова: метод Монте-Карло, системы ОДУ, интегральное уравнение, задачи массового обслуживания, оптимальная плотность, несмещенная оценка, статистическое моделирование.

1. Введение. Метод Монте-Карло для интегральных уравнений и для больших систем линейных алгебраических уравнений достаточно подробно описан в литературе, например [1, 2], а также [3] и [4]. Иным образом дело обстоит с задачами Коши для больших систем обыкновенных дифференциальных уравнений. Вместе с тем даже некоторые классы линейных систем представляют собой значительный интерес для приложения. В качестве примера можно привести задачи массового обслуживания, которые могут описываться системами с очень большим, или даже бесконечным, числом уравнений. В простейшем случае пуассоновских систем легко

(¡5 Санкт-Петербургский государственный университет, 2021

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

Этих соображений достаточно чтобы рассмотреть отличные от имитационных методы Монте-Карло для решения системы

^У(*)=Л(*)У(*) + ^),У( 0) = Уо,

т = к (*щ=1, у (£) = ...,уп(тт, ^ (г) = ит, ..,ытт,

для значений г из промежутка [0, Т].

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

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

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

2. Решение интегральных уравнений Фредгольма методом Монте-

Карло. Пусть задано уравнение относительно ((х)

Ах) = I к(х,у)(р(у)/№у) + /(x), (2)

выполняющееся на носителе X вероятностной! меры и для множества Н функций К(х) таких, что существует интеграл ((, К) = § ((х)К(х)/(¿х). Пусть сходится метод последовательных приближений (имеется в виду сходимость по норме в некотором нормированном пространстве, см. [2, с. 272-275]):

<рт+1(х) = J 11к(х,у)1^рт(у)/^у) + /(х)1 (о = /(х)1; т = 0, l,.... (3) Тогда справедливо представление

((,К) = (/,К) / /о... щКо, &од...&г-1,г/г, (4)

1=1

где обозначаем = /(¿хг), Кг = К(хг), /г = /(хг), кг-1,г = к(хг-1, хг),г = 0,1,..., I.

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

J p(x,y)/(dy) = 1 - д(х), 0 < д(х) < 1, (5)

•Ут(Жо, ..., хт) — ————— —- —— (I)

и условию согласования

ро(х) > 0, если к(х) = 0,

д(х) > 0, если /(х) = 0, (6)

р(х, у) > 0, если к(х, у) = 0.

Здесь и далее подразумеваем х, у € X.

Моделируем цепь Маркова, и на ее траекториях хо ^ х\ ^ ... ^ хт вычисляем функцию

к(хо)к(хо, х\)...к(хт-ххт)/(хт) р°(х0)р(х0, х1)...р(хт-1хт)д(хт)'

Предполагается, что почти все траектории конечные. Это будет так, во всяком случае, если д(х) > 0.

Легко показать (см. [1, с. 94-102]), что при выполнении условий сходимости (3) и условий согласования (6) имеет место равенство

Е ,1т (хо,...,хт ) = (Н,ф). (8)

Для этого достаточно учесть, что плотность вероятности случайной траектории есть р0(хо)р(хо,х\), ...,р(хт_1,хт)д(хт), и воспользоваться определением математического ожидания.

Существует ряд других оценок на траекториях процесса (5), но оценка (7) — одна из простейших в вычислительном отношении оценок, и для нее очень просто получить аналитическое выражение для дисперсии и задать цепь Маркова, удовлетворяющую (3) и (6), для которой значение дисперсии имеет минимальное значение. Эта оптимальная цепь Маркова является однородной и задается

~0, N \КХ)\'ТР{Х) начальной плотностью р 1х) = — — и

(НИ) _

переходной плотностью р(х, у) =-_ -,

<р(х)

где Тр есть предел (3) при то —>• оо в (3).

Интересно отметить (это важно для дальнейшего), что в случае неотрицательных к, / и к оптимальное значение дисперсии оценки (8) оказывается нулем. И хотя выражения (9), определяющие оптимальную цепь Маркова, зависят от точного значения Тр, они позволяют сделать некоторые выводы относительно выбора цепи Маркова при решении реальных задач. Так, мы видим что в состав плотности р(х, у) следует включать сомножителем |к(х, у)|.

3. Решение задачи Коши для системы обыкновенных дифференциальных уравнений методом Монте-Карло. Вернемся к уравнению (1). Элементы матрицы А(Ь) и вектора Р(Ь) будем предполагать ограниченными функциями Ь на вещественной оси. Отсюда следует выполнение условий Липшица для правой части (1) и, следовательно, существование и единственность решения на любом конечном промежутке [0, Т]. Тогда можем переписать (1) в форме

У(Ь) = У(0) + / ^(т)йт + [ А(т)У(т)йт (10)

■1 о ■! о

Y(t) = i А(в)Е(t - 0)Y(в) de + F(t), (11)

Jo

где F(t) = Yo + ft F(r)dr, а E(x) = Л' еСЛИ X > 0 t G [0,T]. 0 [О, если x < 0,

По аналогии, мы можем выбрать цепь Маркова, удовлетворяющую условиям

согласования

p0(i,t) > О, если hi(t) = 0,

g(i,t) > 0, если fi(t)=0, (12)

p(i,t; j,e) > 0, если aij (в) =0,

и построить аналог оценки (7) функционала H(r)Y(r)dr для некоторого фиксированного вектора H(t) = (hi(t),hn(t)). Но далее мы замечаем, что уравнение (10) является уравнением вольтерровского типа, и мажорантный процесс

Ym+1(t)= [ \A(e)\E(t-e)Ym(e)d,e+\F(t)\ (13)

o

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

Определим скалярное произведение двух n-мерных векторов как (A, B) = aibi + ... + anbn. Для построения несмещенной оценки (H,Y) единственным требованием к цепи Маркова остается выполнение условий согласования (12). Если эти условия выполнены, то алгоритм решения задачи состоит в моделировании траектории цепи

io,to ^ ii,ti ^ ... ^ iT,tT

и вычислении оценки

J = _hio (¿oKo,»i ,iT (U)fiT (tT)__^^

Pi0 (to)p(io, to', 4,ti)...p(iT-i,tT-i', iT, tT)g(iT,tT)'

to = T, для которой выполняется равенство

EJt = (H,Y).

Оценка JT равна нулю вне симплекса to > ti > ... > tT. Траектория разыгрывается внутри этого симплекса.

Эти предварительные рассмотрения приводят нас к следующим результатам:

Теорема 1. Пусть задан случайный процесс Маркова с начальной и переходной плотностями

p0(t0)=p°(t0):{p°(e0)=p0(t0,e0)}- i0 = T~n-, t0 = т,

__15

P{t) = \\PiMlj=i--{PiM=Pi4,ek]ik+i,ek+i)}] к = о,т; te [0,Т]

соответственно, для которых выполнено условие согласования (12) и 'равенство

n г t

n rt

V p(i,t;j,e)de = 1 - g(i,t), 0 < g(i,t) < 1. (16)

i=iJo

Тогда для оценки (14) выполнено 'равенство

ЕЛ = (Н,У), (17)

где Н(Ь) — заданный вектор, а У(Ь) — решение уравнения (1).

Доказательство. В соответствии с теоремой Пикара [7, с. 53] имеем

^ г оТ

то г-г г-о т-1

(.Н,У)(Ь) = (Н,Р)(Ь) +V/ а91... 3,втН(в1)А(в2)...А(вт)Р(вт) = т =1-1о -)о

то п п . г . от-1

= ...^ а9о... а9т К (Оо)<Н0м (в1)...04т-1 лт (9т )/т (9т).

Имея заданную цепь Маркова, получаем N траекторий 'о,Ьо ^ 4^1 ^ ... ^ 'т,Ьт с плотностью

ро('о,Ьо)р('о,ьо;'1,Ь1)...р('т-1,Ьт-1; 'т,Ьт)д(''тМ). (18)

Совокупность этих плотностей определяет вероятностную меру Р(¿9) на множестве траекторий. По условию (16) почти все траектории являются конечными, поэтому выполняется равенство

то п п г рвТ-1

/ а9о... а9тро('о,9о)р('о^о;'1,9{)...р(1т-1,Ьт-1;'т ,9т )д('т,Ьт) = 1.

т=о ¿о = о ¿т = о о •/о

Пользуясь определением математического ожидания, при выполнении условий согласования (12) для оценки (14) получаем

то п п .г „вТ-

-1

У) У)...У) а9о... а9т Зтр0('о,9о)р('оМ;'1,9{)...р(1т-1Ьт-1; ¡т ,9т )д('т ,Ьт ) =

т=о ¿0 =о ¿Т =о}о ^

то п п г ,-вТ-1

= ¿9о... ¿9тk¿о (9о)^1 (91)-^т-1 ¿Т (9т)к (9т)■

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

т=о ¿0=о ¿Т=о-]о -]о

Очевидным отличием от рассматриваемого в начале статьи классического случая является:

1) выполнение условия (16);

2) автоматическое выполнение мажорантного условия (13). При выполнении условий теоремы 1 справедлива также

Теорема 2. Дисперсия оценки (14) конечна и равна

'Н? Р*

(19)

т А%(9)Е2(Ь - 9)_ р2(Ь) где 2(4) есть решение уравнения 2(4) = —--——-2(6>)сЙ? -|--а пота~

ция А2%; Р\; ро; p¿j; g¿ подразумевает покомпонентные операции над матрицами и векторами (операции Адамара).

Доказательство. Воспользуемся стандартным выражением для дисперсии:

ВЗТ = ЕО2 - [ЕОТ]2.

Из результатов теоремы 1 мы знаем, что [Е,1Т]2 = (Н, У)2. Обратимся к выражению Ео/Т2. Возводя в квадрат, получаем выражение

ТО П П г. £

т=0 го = 0 2Т

^ л-1 ^ ^(^оКа^О-Ч-л^^) (20)

.]О тр°{ъ0, в0)р{%0, ¿о; «1, 6»1)...р(гт_1, гт, вт)д{гт, 6Т)

(Н2 Л ,т А2а(е)Е2(г - в)„ Р2(г)

для скалярного произведения I —2 I, где 2(4) = _/0 —---Е(в)с1в-\—г—гт>

V р ) Рг3 (") д2(г)

или равносильную систему дифференциальных уравнений а А2 р2

т = + 5(0)= Но.

аг Ргз дг

В соответствии с теоремой Пикара, решение этой системы уравнений существует и единственно:

Н2 = А ГШо Г^-1 т Н2(в0)А2](в1) А23(вТ)Р^вТ) __1 -'о -'о

Н22

¿е.

р°; Т=^о Jо Р0(ео) Р23№)'рц(ет) дг(ет)'

Таким образом, = □

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

Теорема 3. Дисперсия оценки ОТ принимает минимальное значение при выборе следующих параметров цепи Маркова:

.о „

Чопт(г,г; е)

(|Я|,У)

\^{9)\Е{Ь-9)щ{9) [ }

Ш

где У — решение уравнения У(4) = /0Т \А{в)\Е{Ь — в)Ут(в)<1в + |-Р(4)|; а Уг — г-й компонент У.

Доказательство. Обратимся к оценке (14). Насколько мы видели, с помощью этой оценки для каждой цепи Маркова, удовлетворяющей условиям согласования, имеем

(Н,У)(г) = ЕОт(г) = ( ОтКв)р(ае),

о

где Р — мера на траектории марковской цепи, такая что плотность распределения траектории по отношению к мере ^Т+1 есть выражение (18). Предположим, что мы

¿р

хотим ввести другую меру (^(¿6), такую что производная Радона — Никодима —

аЦ

существует. Для этого мера Р должна быть абсолютно непрерывна по отношению к мере Ц. Другими словами, она должна иметь сходную с Р структуру: траектория шт,в имеет плотность чт(г0,Ьо;...; гт, Ьт) по отношению к ^т +1, и чт должна быть больше нуля для тех траекторий шт,в, для которых выражение плотности (18) принимает строго положительные значения.

В соответствии с теоремой о существенной выборке [1, с. 106], для оптимальной плотности имеем

Чт, опт(^т,в) = С\Х\р0(г,Ь)р(го,Ьо; г1,Ь1)...р(гт-1и-1; г^т), что аналогично

Чт, опт (^т,в) = С\Нг0 (Ьо )йг0 (^...й^ — (Ьт )/гт (и )\,

где С — константа нормировки,

то п п . дт-1

С-1 = ]Т Е ... Е / ¿°о... авк\Ьг° (°о)агоМ (О1)...сцт- (От)к (От)\,

т=0 го = 0 гт =0 0 Л

с-1 = (\Н\,Ту,

а У удовлетворяет уравнению (13).

Следует заметить, что необязательно должна существовать цепь Маркова, индуцируемая Ц. Однако в рассматриваемом случае представляется возможным указать цепь Макрова, индуцирующую такую оптимальную меру. Из уравнения (13) следует

]Г Г Мв)\Е(г - о)у0{в)<1в = ш - \Ш\, ВД < 1,

¿Л Ш ~ УМ •

тт ел

Принимая допт(г,1) = _ , получаем выражение

УМ

п р г

У2 Чопт(г,Ь;з, в)¿в = 1 - допт(г,Ь).

0

3 = 1

Другими словами, условие (16) выполнено для вновь индуцируемой меры Ц. Выражение плотности приобретает вид

*7т, ОПТ )

(|Я|,У)

Из уравнения (13) очевидно, что Ч<0пт(г,Ь) — неотрицательная функция, и она нормирована (то есть является плотностью), и функция чопт(г,Ь;з,О) неотрицательна. Таким образом, для меры Ц выполняются условия согласования (12).

Используя неравенство Коши — Буняковского, находим оценку снизу для EJT2:

Sl (^)2ди) - =№y)2(t)-

Подставляя qT, опт(^т,в) в выражение для EJT2, получаем гt Л dP t't „ 1

Е J2T(t) = J2T(1FM)2Q(de)= J2T-—p--C\JT\PT(oJT,e) = Jo dQ J0 C2\J.T \

= fo Щ-рт = (\HiY)(t). j*\jT\p(de) = mxm

и убеждаемся в справедливости теоремы 3. □

Полученные выше результаты позволяют нам вывести

Следствие. Если aitj(t) > 0 и fi(t) > 0 для всех i,j и t, то оптимальное значение дисперсии оценки JT есть нуль.

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

(Я,У)2 = (|Я|,У)2, и тогда в выражении (19) получаем

DJT = (|Я|,У)2 - (Я, У)2 = 0.

Заметим, что для справедливости результатов аналогичных теореме 3, зависимость вероятности поглощения от t является существенной. Поясним это на следующем простом примере. Пусть F(t) = 0, и A(t) = \\aijj \\ij=i не зависит от t. Тогда решение может быть записано в замкнутой форме:

то Ak tk

y(i)=exP(Ai) = ]r— (22)

k=0 '

Скалярное произведение (Y, H) можно вычислить с помощью метода Монте-Карло, моделируя однородную цепь Маркова р°, Р, для которой выполнены условия согласования (12), и вычисляя оценку

т _ К aip,il ■■■ aiT-l,iT УъТ

JT — —g {¿6)

Pi0 Pia ,ii ■■■Pir-1 ,ir 3ir • T ■

Легко проверить (см., например, [6]), что при aitj > 0, у0т > 0 и hi0 > 0 нельзя выбрать цепь fP; Р, обращающую дисперсию в нуль.

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

4. Постоянная матрица коэффициентов. Приведем иллюстративный пример. Возьмем уравнение системы массового обслуживания, описывающее замкнутую систему из M =10 условных станков, которые выходят из строя с интенсивностью Л, и S = 7 рабочих, которые чинят станки с интенсивностью ц. Пользуясь [8, с. 102], получаем следующую систему уравнений:

3/0 (t) = -10Ay0(t)+ m/7 3i(t),

31 (t) = - [9A + m/7] yi (t) + 10A yo (t) + 2m/7 /2 (t),

32 (t) = - [8A + 2^/7] /2 (t) + 9A yi (t) + 3m/7 33 (t),

3/3 (t) = - [7A + 3m/7] уз (t) + 8A /2 (t) + 4m/7 34 (t),

34 (t) = - [6 A + 4m/7] 34 (t) + 7A уз (t) + 5m/7 35 (t),

//5 (t) = - [5A + 5m/7] 35 (t) + 6A 34 (t) + 6m/7 36 (t), (24)

//6 (t) = - [4A + 6m/7] 36 (t) + 5A 35 (t) + m 37 (t),

37 (t) = - [3A + m] 37 (t) + 4A 36 (t) + m 38 (t),

38 (t) = - [2A + m] 38 (t) + 3A 37 (t) + m 39 (t),

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

39 (t) = - [A + m] 39 (t) + 2A 38 (t) + m 310 (t),

310(t) = -M 3io(t) + A 39(t).

В рамках данной статьи был рассмотрен конкретный пример со следующими параметрами: ^ = 0.0202 и Л = 0.02, а также ^ = 0.0202 и Л = 0.01. Такая комбинация параметров дает нам значения коэффициента зарузки (см. [8, с. 104]) Ф = 0.99 и Ф = 0.49 соответственно. При коэффициенте загрузки Ф = 0.99 имитационные методы работают плохо. Второй вариант значения коэффициента загрузки демонстрирует работоспособность предлагаемого метода в сравнении с методом Рунге — Кутта в более конвенциональных условиях. Выберем начальную позицию y0 = 1; y0 = 0, Vi = 2; T = 1.5; количество траекторий установим равным N = 1 800 000.

Для сравнения результатов вычислений используем библиотеку SciPy программного языка Python [9]. Воспользуемся рутиной scipy.integrate.ode("dop853"), которая реализует явный алгоритм Дорманда — Принса порядка 8(5, 3) для метода Рунге — Кутта. Данный алгоритм является адаптивным алгоритмом, автоматически изменяющим шаг интегрирования в зависимости от того, попадает ли локальная ошибка в требуемое значение погрешности. Погрешность рассчитывается, исходя из формулы AbsTol + RelTol * |y step |, где AbsTol — абсолютная погрешность, RelTol — относительная погрешность, а y step — рассчитанное решение на каждом шаге (см. [10, с. 188-204]). Нами были установлены значения AbsTol = 1e-06; RelTol = 1e-04.

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

s„ = _ и значения абсолютной фактической ошибки х = \y(t) ~ yRK(t) IN

Результаты вычислений приведены в табл. 1 и 2, а также в табл. 3 и 4. При построении таблиц было принято решение не округлять результаты расчетов, полученных методом Рунге — Кутта, но исключить из расчета х те yRK(t), чья значащая цифра по порядку меньше либо равна значению AbsTol. На этих местах в таблице х стоит прочерк. Нетрудно также увидеть, что значения, полученные методом Монте-Карло, в некоторых состояниях меньше или равны sn, что объясняется их близостью к нулю.

Таблица 1. Сравнение решения методом Монте-Карло и методом Рунге — Кутта при Т = 1.5, Ф = 0.99, _N = I 800 ООО_

Монте-Карло Рунге — Кутта

У0 6.584041е-06 1.425017е-05

У1 6.602467е-03 6.672860е-03

У2 7.843779е-01 7.819781е-01

УЗ 1.900002е-01 1.898743е-01

У4 2.019594е-02 2.018617е-02

УЪ 1.276377е-03 1.226557е-03

Ув 4.864802е-05 4.658389е-05

У7 3.336787е-07 1.132339е-06

У8 1.867346е-08 1.721394е-08

У9 9.384152е-11 1.49590бе-10

У10 < «п 5.691357е-13

Таблица 2. Сравнение доверительного интервала метода Монте-Карло и абсолютной погрешности между решениями МК и РК при Т = 1.5, _Ф = 0.99, N = 1 800 000_

X

У0 8.884365е-06 7.666127е-06

У1 3.728903е-04 7.039296е-05

У2 4.111465е-03 2.399780е-03

УЗ 1.706013е-03 1.258374е-04

У4 3.551892е-04 9.779131е-06

УЪ 5.111231е-05 4.981978е-05

Ув 4.754939е-06 2.064131е-06

У7 1.866752е-07 -

У8 1.055185е-08 -

У9 8.037395е-11 -

У10 1.036435е-16 -

Таблица 3. Сравнение решения методом Монте-Карло и методом Рунге — Кутта при Т = 1.5, Ф = 0.49, _N = 1 800 000_

Монте-Карло Рунге — Кутта

У0 1.400258е-05 1.630445е-05

У1 7.620840е-03 7.575809е-03

У2 8.787064е-01 8.804753е-01

УЗ 1.058588е-01 1.061585е-01

У4 5.537256е-03 5.601951е-03

УЪ 1.644043е-04 1.689383е-04

У6 2.892959е-06 3.184300е-06

У7 4.260589е-08 3.841370е-08

У8 < «п 2.898 ЮОе-10

У9 < «п 1.249898е-12

У10 < «п 2.359078е-15

Таблица 4- Сравнение доверительного интервала метода Монте-Карло и абсолютной погрешности между решениями МК и РК при Т = 1.5, _Ф = 0.49, N = 1 800 000_

X

Уо 4.204214е-06 2.301872е-06

У1 2.519817е-04 4.503130е-05

У2 3.694861е-03 1.768928е-03

Уз 8.335724е-04 2.996921е-04

У4 9.236051е-05 6.469456е-05

УЪ 6.750112е-06 4.533990е-06

У6 3.615292е-07 -

У7 1.237473е-08 -

У8 2.560555е-11 -

У9 5.164668е-14 -

У10 4.225728е-24 -

Таблица 5. Сравнение решения методом Монте-Карло и методом Рунге — Кутта при Т = 1.5,

Щ) = Ар/(1 +£2), Ар = 0.02, N = 2 200 000

Монте-Карло Рунге — Кутта

У0 1.571291е-05 1.569973е-05

У1 7.285691е-03 7.297974е-03

У2 8.508264е-01 8.486886е-01

УЗ 1.342014е-01 1.343154е-01

У4 9.442907е-03 9.304705е-03

УЪ 3.745529е-04 3.683801е-04

У6 9.152182е-06 9.115712е-06

У7 1.262589е-07 1.443672е-07

У8 1.206253е-09 1.430467е-09

У9 1.160412е-11 8.104518е-12

У10 1.902090е-16 2.009527е-14

Таблица 6. Сравнение доверительного интервала метода Монте-Карло и абсолютной погрешности между решениями МК и РК при Т = 1.5,

\{Ь) = Ар/(1 +£2), Ар = 0.02, N = 2 200 000

X

У0 1.327831е-07 1.318437е-08

У1 3.981266е-05 1.228340е-05

У2 3.025420е-03 2.137891е-03

Уз 1.408177е-03 1.139761е-04

У4 2.221353е-04 1.382022е-04

УЪ 1.686217е-05 6.172786е-06

Ув 8.709463е-07 -

У7 2.493599е-08 -

У8 5.125696е-10 -

У9 7.929071е-12 -

У10 8.923702е-17 -

5. Переменная матрица коэффициентов. Продемонстрируем возможность предложенного метода для решения системы уравнений с матрицей коэффициентов, зависимой от времени. Возьмем систему (24) и сделаем параметр Л зависимым от Ь: Л(Ь) = Ло/(1 + Ь2). Сохраним при этом начальный параметр Ло = 0.02. Алгоритм вычисления функционала общего вида (И,У) повторяет алгоритм, описанный выше, для случая с постоянной матрицей коэффициентов, однако с необходимостью учитывать изменения матрицы коэффициентов А(Ь) с течением времени.

Оценка примет вид (14) при условии, что мы оставим прежние параметры

для переходной и начальной плотностей распределения, индуцирующих цепь Мар-

Ръ

кова: р° = р^ для начального распределения; р(в) : {р(1,вк'^,0к+1) = —для

Ьг

плотности перехода. Расчет проводим с параметрами ^ = 0.0202, Л0 = 0.02, Т = 1.5, у20 = 1; у0 = 0, Уг = 2, и N = 2200000. Результаты вычислений продемонстрированы в табл. 5 и 6.

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

Литература

1. Ермаков C.M. Метод Монте-Карло в вычислительной математике. Санкт-Петербург, Невский Диалект (2009).

2. Ермаков C.M. Метод Монте-Карло и смежные вопросы. Москва, Наука (1975).

3. Ермаков C.M., Сипин А. С. Метод Монте-Карло и параметрическая разделимость алгоритмов. Санкт-Петербург, Изд-во С.-Петерб. ун-та (2014).

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

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

5. Ermakov S., Pogosian A. On solving stochastic differential equations. Monte Carlo Methods and Applications 25, iss. 2, 155-161 (2019). https://doi.org/10.1515/mcma-2019-2038

6. Ермаков С.М., Товстик Т. М. Метод Монте-Карло для решения систем ОДУ. Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия 6 (64), вып. 3, 411421 (2019). https://doi.org/10.21638/11701/spbu01.2019.306

7. Бибиков Ю.Н. Курс обыкновенных дифференциальных уравнений. Санкт-Петербург, Лань (2011).

8. Кофман А., Крюон Р. Массовое обслуживание: теория и приложения, пер. с франц. Москва, Мир (1965).

9. Virtanen P., Gommers R., Oliphant T. E. et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat. Methods 17, 261-272 (2020). https://doi.org/10.1038/s41592-019-0686-2

10. Hairer E. Solving ordinary differential equations I: Nonstiff problems. Springer (1993).

Статья поступила в редакцию 3 июня 2020 г.;

после доработки 27 июля 2020 г.; рекомендована в печать 17 сентября 2020 г.

Контактная информация:

Ермаков Сергей Михайлович — д-р физ.-мат. наук, проф.; sergej.ermakov@gmail.com Смиловицкий Максим Григорьевич — аспирант; smmxgr@gmail.com

Monte-Carlo for solving large linear systems of ordinary differential equations

S. M. Ermakov, M. G. Smilovitskiy

St. Petersburg State University, 7—9, Universitetskaya nab., St. Petersburg, 199034, Russian Federation

For citation: Ermakov S. M., Smilovitskiy M. G. Monte-Carlo for solving large linear systems of ordinary differential equations. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy, 2021, vol. 8(66), issue 1, pp. 37-48. https://doi.org/10.21638/spbu01.2021.104 (In Russian)

Monte-Carlo approach towards solving Cauchy problem for large systems of linear differential equations is being proposed in this paper. Firstly, a quick overlook of previously obtained results from applying the approach towards Fredholm-type integral equations is being made. In the main part of the paper, a similar method is being applied towards a linear system of ODE. It is transformed into an equivalent system of Volterra-type integral equations, which relaxes certain limitations being present due to necessary conditions for convergence of majorant series. The following theorems are being stated. Theorem 1 provides necessary compliance conditions that need to be imposed upon initial and transition distributions of a required Markov chain, for which an equality between estimate's expectation and a desirable vector product would hold. Theorem 2 formulates an equation that governs estimate's variance, while theorem 3 states a form for Markov chain parameters that minimise the variance. Proofs are given, following the statements. A system of linear ODEs that describe a closed queue made up of ten virtual machines and seven virtual service hubs is then solved using the proposed approach. Solutions are being obtained both for a system with constant coefficients and time-variable coefficients, where breakdown intensity is dependent on t. Comparison is being made between Monte-Carlo and Rungge — Kutta — obtained solutions. The results can be found in corresponding tables.

Keywords: Monte-Carlo, ODE system, integral equation, queuing theory, optimal density, unbiased estimate, statistical modelling.

References

1. Ermakov S.M. Monte-Carlo for quantitative mathematics. St Petersburg, Nevsky Dialect Publ. (2009). (In Russian)

2. Ermakov S.M. Monte-Carlo and related topics. Moscow, Nauka Publ. (1975). (In Russian)

3. Ermakov S. M., Sipin A. S. Monte-Carlo and parametrical divisibility of algorithms. St. Petersburg, St. Petersburg Univ. Press (2014). (In Russian)

4. Ermakov S. M., Mikhailov G. A. Statistical modelling. Moscow, Nauka Publ. (1982). (In Russian)

5. Ermakov S., Pogosian A. On solving stochastic differential equations. Monte Carlo Methods and Applications 25, iss. 2, 155-161 (2019). https://doi.org/10.1515/mcma-2019-2038

6. Ermakov S. M., Tovstik T. M. Monte-Carlo method for solution of systems ODE. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy 6 (64), iss. 3, 411-421 (2019). https://doi.org/10.21638/11701/spbu01.2019.306 (In Russian) [Engl. transl.: Vestnik St. Petersb. Univ. Math. 52, iss. 3, 272-280 (2019). https://doi.org/10.1134/S1063454119030087].

7. Bibikov J.N. A course in ordinary differential equations. St. Petersburg, Lan' Publ. (2011). (In Russian)

8. Kaufmann A., Cruon R. Los fenomenos de espera; teoria y aplicaciones. Continental (1966). [Russ. ed.: Massovoe obsluzhivanie: teorija i prilozhenija. Moscow, Mir Publ. (1965)].

9. Virtanen P., Gommers R., Oliphant T. E. et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat. Methods 17, 261-272 (2020). https://doi.org/10.1038/s41592-019-0686-2

10. Hairer E. Solving ordinary differential equations I: Nonstiff problems. Springer (1993).

Received: June 3, 2020 Revised: July 27, 2020 Accepted: September 17, 2020

A u t h o r s' i n fo r m a t i o n:

Sergey M. Ermakov —sergej.ermakov@gmail.com Maxim G. Smilovitskiy — smmxgr@gmail.com

i Надоели баннеры? Вы всегда можете отключить рекламу.