=Paper=
{{Paper
|id=Vol-2064/paper04
|storemode=property
|title=
Анализ трудностей реализации стохастических численных методов Рунге—Кутты
(Implementation difficulties analysis of stochastic numerical Runge-Kutta methods)

|pdfUrl=https://ceur-ws.org/Vol-2064/paper04.pdf
|volume=Vol-2064
|authors=Dmitry Kulyabov,Migran Gevorkyan,Anastasiya Demidova,Anna Korolkova,Leonid Sevastianov,Mikhail Kotukov
}}
==
Анализ трудностей реализации стохастических численных методов Рунге—Кутты
(Implementation difficulties analysis of stochastic numerical Runge-Kutta methods)
==
<pdf width="1500px">https://ceur-ws.org/Vol-2064/paper04.pdf</pdf>
<pre>
УДК 519.216.2
Кулябов Д.С.1,2, Геворкян М.Н.1, Демидова А.В.1, Королькова А.В.1, Севастьянов Л.А.1,
Котюков М.М.1
1 Российский университет дружбы народов, г. Москва, Россия
2 Объединенный институт ядерных исследований, г. Дубна, Россия

АНАЛИЗ ТРУДНОСТЕЙ РЕАЛИЗАЦИИ СТОХАСТИЧЕСКИХ ЧИСЛЕННЫХ МЕТОДОВ
РУНГЕ—КУТТЫ*
Аннотация
В данной статье рассматриваются стохастические численные методы типа Рунге-Кутты
с слабой и сильной сходимостями для систем стохастических дифференциальных уравнений
в форме Ито. Дается краткий обзор основных работ по данной теме. Также кратко
излагается теоретический материал по стохастическим численным методам и сведения
из теории стохастических дифференциальных уравнений. В связи с тем, что
аналитическому решению поддается лишь малая доля стохастических дифференциальных
уравнений, стохастические численные методы являются единственным способом
получения решения. Стохастические методы Рунге-Кутта основываются на тех же идеях,
что и классические методы, однако их схемы существенно сложнее своих
детерминированных аналогов. Это обуславливает трудности при попытки программной
реализации. Авторами мотивируется подход к реализации данных методов с помощью
генерации исходного кода так как это позволяет добиться универсальности при сохранении
простоты кода. Обсуждаются детали реализации и используемые языки программирования
и библиотеки. Приводится ряд примеров стохастических дифференциальных уравнений,
используемых для тестирования программы и оценки погрешности вычислений. Дается
ссылка на репозиторий с исходным кодом программы на языке Python.
Ключевые слова
Стохастические дифференциальные уравнения; стохастические численные методы;
автоматическая генерация кода; Python; Julia; шаблонизатор.
Kulyabov D.S.1,2, Gevorkyan M.N.1, Demidova A.V.1, Korolkova A.V.1, Sevastianov L.A.1,
Kotukov M.M.1
1 Peoples’ Friendship University of Russia, Moscow, Russia
2 Joint Institute for Nuclear Research, Dubna, Moscow region, Russia

IMPLEMENTATION DIFFICULTIES ANALYSIS OF STOCHASTIC NUMERICAL RUNGE-
KUTTA METHODS
Abstract
This article discusses stochastic numerical methods of Runge-Kutta methods with weak and strong
convergences for systems of stochastic differential equations in the Ito form. At first section we give a
brief review of main works on the topic. In next section we make short introduction into theory of
stochastic numerical methods and facts from the theory of stochastic differential equations. Due to
the fact that the analytical solution exists only for a small subset of stochastic differential equations,
stochastic numerical methods are the only way to obtain the solution. Stochastic Runge-Kutta
methods are based on the same ideas as the classic methods, but their scheme is essentially more
complicated than their deterministic counterparts. This leads to difficulties when one attempts to
implement them. In this paper we motivate the approach to the implementation of these methods
using source code generation as it allows achieving universality while maintaining a simplicity of code.
We discuss some of the implementation details and the used programming languages and libraries.

* Труды II Международной научной конференции «Конвергентные когнитивно-

информационные технологии» (Convergent’2017), Москва, 24-26 ноября, 2017
Proceedings of the II International scientific conference "Convergent cognitive information
technologies" (Convergent’2017), Moscow, Russia, November 24-26, 2017

28
We give several examples of stochastic differential equations, used to test the program and evaluate
the error of the calculations. We provide the link to the repository with the source code of discussed
programs.
Keywords
Stochastic differential equations; stochastic numerical methods; automatic code generation; Python;
Julia; the template engine.

Введение
В статье авторов [1] описывалась реализация стохастических численных методов типа Рунге-Кутты на
языке Python [2] с использованием библиотек NumPy и SciPy [3]. Выбор языка был продиктован удобством
программирования и возможностью работать с многомерными массивами как с тензорами (функции
tensor_dot и einsum из библиотеки NumPy). Однако производительность созданных функций была на
низком уровне, и не столько из-за низкого быстродействия языка Python, сколько из-за большого числа
вложенных циклов (до семи циклов). В данной статье мы рассматриваем альтернативных подход к
реализации стохастических численных методов, основанный на генерации кода функций.
Данная статья состоит из трех разделов. В первом разделе дается обзор основных источников и
излагаются сведения из теории стохастических дифференциальных уравнений (СДУ) и численных
методов для их решения. Во втором разделе приведены стохастические численные схемы для скалярных
СДУ с сильной сходимостью и для систем СДУ с сильной и слабой сходимостью. Наряду с общими схемами
приводятся несколько таблиц коэффициентов, которые позволяют реализовать уже конкретный
численный метод. Наконец в третьем разделе мотивируется применение генерации кода для
стохастических численных методов и описываются некоторые детали реализованного нами генератора
(на языке Python с использованием шаблонизатора Jinja2).
Обзор основных источников
Кратко охарактеризуем основные доступные статьи, посвященные стохастическим методам. Особое
внимание уделим многостадийным методам. Первым стохастические броуновские процессы для
математического моделирования использовал французский математик Л. Ж.-Б. А. Башелье (1870-1946) в
1900 году в работе [4].
В качестве краткого введения в область стохастических численных методов можно использовать
обзорную статью А. В. Лукшина и С. Н. Смирнова [5], в которой тезисно изложены подходы к численному
решению стохастических уравнений по состоянию на 1990 год. В более новой статье [6] А. Ф. Ерешко и
Д.В. Филатовой также приведены некоторые численные методы (в том числе и один метод Рунге-Кутты
и методы Г. Н. Милштейна [7-9]), а также даны результаты численного эксперимента.
Используемые в данной статье численные схемы в основном были получены А. Росслером, который в
работе [10] последовательно изложил теорию стахостических методов Рунге-Кутты. Автор рассматривает
аппроксимацию систем СДУ Ито и Стратоновичав в слабом смысле для скалярного и многомерного
процесса Винера. После краткого обзора литературы автор переходит к изложению теории помеченных
деревьев, которая используется для вывода условий порядка в случае детерминированных методов
Рунге-Кутты [11] и которая была распространена на случай СДУ. Дальнейшие результаты Росслера
изложены в статьях [12,13]. В препринте [13] продолжена классификация стохастических методов Рунге-
Кутты со слабой сходимостью на основе условий порядка. Приведено несколько конкретных реализаций
и даны результаты численных экспериментов. В другом препринте [14] приведены таблицы методов
четвертой стадийности и сильного порядка сходимости 𝑝 = 2.0.
П.М. Бурраж и К. Бурраж в серии статей [15-19] изучили методы сильного порядка 𝑝 = 1.5, а также
распространили теорию помеченных деревьев на стохастический случай. В статье R. Soheili и M. Namjoo
[20] получены три метода с сильной сходимостью 𝑝 = 1.0 и проведено численное сравнение с методом из
книги [21]. Дальнейшее развитие методы со слабой сходимостью получили в статьях [22] и [23]. В статье
[24] получены два трехстадийных метода слабой сходимости 𝑝 = 2.0 , а также проведены численные
эксперименты.
В виду крайней сложности дальнейшего повышения порядка точности стохастических численных
схем, современные работы в основном концентрируются на получении численных схем для частных
случаев СДУ. Можно выделить работы, посвященные симплектическим стохастическим численным
методам Рунге–Кутта [27–29] и неявным схемам, в частности стохастическим аналогам метода
Розенброка [30].
В качестве монографии, посвященной целиком стохастическим численным методам, следует
упомянуть книгу [21], которая содержит подробную сводку разнообразных численных схем, а также
список СДУ, которые имеют аналитические решения.

29
Случайный процесс Винера и программное генерирование его траекторий
Случайный процесс 𝑊(𝑡), 𝑡 ≥ 0 называется скалярным процессом Винера, если выполняются
следующие условия [25, 32, 33]:
• 𝑃{𝑊(0) = 0} = 1 иначе говоря, 𝑊(0) = 0 почти наверное;
• 𝑊(𝑡) — процесс с независимыми приращениями, то есть приращения Δ𝑊𝑖 являются
независимыми случайными величинами
• ∆𝑊𝑖 =𝑊(𝑡𝑖+1 )-𝑊(𝑡𝑖 )∼𝒩(0,𝑡𝑖+1 -𝑡𝑖 ), где 0 ≤ 𝑡𝑖+1 < 𝑡𝑖 < 𝑇, 𝑖 = 0, 1, ⋯ , 𝑁 − 1.
Обозначение ∆𝑊𝑖 ~𝒩(0, Δ𝑡𝑖 ) говорит о том, что Δ𝑊𝑖 — нормально распределенная случайная величина
с математическим ожиданием 𝔼[∆𝑊𝑖 ] = 𝜇 = 0 и дисперсией 𝔻[Δ𝑊𝑖 ] = 𝜎 2 = Δ𝑡𝑖
Винеровский процесс является моделью броуновского движения (хаотического блуждания). Если
рассмотреть процесс 𝑊(𝑡) в те моменты времени 0 = 𝑡0 ≤ 𝑡1 < 𝑡2 <. . < 𝑡𝑁−1 < 𝑡𝑁 , когда он испытывает
случайные аддитивные изменения, то непосредственно из определения следует:
𝑊(𝑡1 ) = 𝑊(𝑡0 ) + Δ𝑊0 , 𝑊(𝑡2 ) = 𝑊(𝑡1 ) + Δ𝑊1 , … , 𝑊(𝑡𝑁 ) = 𝑊(𝑡𝑁−1 ) + Δ𝑊𝑁−1 ,
где по Δ𝑊𝑖 ~𝒩(0, Δ𝑡𝑖 ), ∀𝑖 = 0, … , 𝑁 − 1.
Запишем 𝑊(𝑡𝑁 ) в виде накопленной суммы:
𝑁−1

𝑊(𝑡𝑁 ) = 𝑊(𝑡0 ) + ∑ 𝑊𝑖
𝑖=0
и учтем, что 𝔼[Δ𝑊𝑖 ] = 0 и 𝔻[Δ𝑊𝑖 ] = Δ𝑡𝑖 можно показать, что сумма нормально распределенных случайных
чисел Δ𝑊𝑖 также является нормально распределенным случайным числом:
𝑁−1 𝑁−1 𝑁−1 𝑁−1

𝔼 ∑ Δ𝑊𝑖 = 0, 𝔻 ∑ Δ𝑊𝑖 = ∑ Δ𝑡𝑖 = 𝑇, ∑ Δ𝑊𝑖 ~𝒩(0, 𝑇).
𝑖=0 𝑖=0 𝑖=0 𝑖=0
Многомерный винеровский процесс 𝑊(𝑡): Ω × [𝑡0 , 𝑇] → ℝ𝑚 определяют как случайный процесс
составленный из совместно независимых одномерных винеровских процессов 𝑊1 (𝑡), … , 𝑊𝑚 (𝑡).
Приращения Δ𝑊𝑖𝛼 , ∀𝛼 = 1, … , 𝑚 являются совместно независимыми нормально распределенными
случайными величинами. С другой стороны? вектор Δ𝑊𝑖𝛼 можно представить как многомерную
нормально распределенную случайную величину с вектором математических ожиданий 𝜇 = 1 и
диагональной матрицей ковариаций. Рассмотрим теперь вопрос генерации траектории процесса Винера
на компьютере. Для симулирования одномерного винеровского процесса необходимо сгенерировать 𝑁
нормально распределенных случайных чисел 𝜇 = 𝜀1 , … , 𝜀𝑁 и построить их накопленные суммы 𝜀1 , 𝜀1 +
𝜀2 , 𝜀1 + 𝜀2 + 𝜀3 и так далее, то мы получим выборочную траекторию винеровского процесса 𝑊(𝑡)
(независимого случайного блуждания).
В случае многомерного случайного процесса следует сгенерировать уже 𝑚 последовательностей из 𝑁
нормально распределенных случайных величин.
Стохастические интегралы и СДУ для скалярного процесса Винера
Формулировка стохастического дифференциального уравнения начинается с введения операции
стохастического интегрирования.
Определение. Пусть 𝑔(𝑡, 𝑥(𝑡)) непрерывная скалярная функция, зависящая от случайного процесса
𝑥(𝑡) , 𝑊(𝑡) — процесс Винера, 𝑡 ∈ [𝑡0 , 𝑇], 0 ≤ 𝑡0 < 𝑡1 < 𝑡2 < ⋯ < 𝑡𝑛−1 < 𝑡𝑛 ≤ 𝑇 < ∞ — разбиение
отрезка [𝑡0 , 𝑇], тогда
𝑇 𝑛−1

𝐼(𝜃) =: ∫ 𝑔(𝑡, 𝑥(𝑡))𝑑𝑊(𝑡) ≝ lim ∑ 𝑔(𝜃𝑡𝑖+1 + (1 − 𝜃)𝑡𝑖 , 𝜃𝑥(𝑡𝑖+1 ) + (1 − 𝜃)𝑥(𝑡𝑖 ))(𝑊(𝑡𝑖+1 ) − 𝑊(𝑡𝑖 )),
𝑛→∞
𝑡0 𝑖=0
где 𝜃 ∈ [0,1] . 𝐼(𝜃) называется стохастическим интегралом. Строгое обоснование существования
данного интеграла для широкого класса функций см. [25, глава 3].
В физике и прикладной математике нашли применение два частных случая стохастического
интеграла:
• интеграл Ито (назван в честь японского математика К. Ито)
𝑛−1

𝐼 = 𝐼(0) = lim ∑ 𝑔(𝑡𝑖 , 𝑥(𝑡𝑖 ))(𝑊(𝑡𝑖+1 ) − 𝑊(𝑡𝑖 )),
𝑛→∞
𝑖=0
• интеграл Стратоновича (назван в честь советского физика Р. Л. Стратоновича)
𝑛−1
𝑡𝑖+1 + 𝑡𝑖 𝑥(𝑡𝑖+1 ) + 𝑥(𝑡𝑖 )
𝐼0.5 = 𝐼(0.5) = lim ∑ 𝑔 ( , ) (𝑊(𝑡𝑖+1 ) − 𝑊(𝑡𝑖 )).
𝑛→∞ 2 2
𝑖=0
После введения операции стохастического интегрирования, можно записать интегрального
уравнения для случайного процесса 𝑥(𝑡) [25].

30
𝑡 𝑡

𝑥(𝑡𝑘 ) = 𝑥(𝑡0 ) + ∫ 𝑓(𝜏, 𝑥(𝜏))d𝜏 + ∫ 𝑔(𝜏, 𝑥(𝜏))d𝑊.
0 ⏟
0
стохастический интеграл
Вышеприведенному интегральному уравнению соответствует стохастическое дифференциальное
уравнение:
d𝑥(𝑡) = 𝑓(𝑡, 𝑥(𝑡))d𝑡 + 𝑔(𝑡, 𝑥(𝑡))d𝑊.
При этом необходимо помнить, что дифференциал типа d𝑥 не является обычным «малым»
приращением функции x(t) это случайная величина. Фактически СДУ является символической записью
предела итерационной схемы.
Для численных схем необходимо уметь вычислять стохастические интегралы Ито следующего вида:
𝑡
𝑊 2 (𝑡) 𝑡
∫ 𝑊(𝜏)𝑑𝑊(𝜏) = − .
2 2
𝑡0
Мы сразу привели конечный результат так как полная процедура вычисления данного интеграла
громоздка.
СДУ в форме Ито для многомерного процесса Винера
Введем вероятностное пространство (Ω, 𝒜, ℙ), где Ω — пространство элементарных событий, 𝒜 —
сигма-алгебра подмножеств пространства Ω , ℙ — вероятностная мера. На отрезке [𝑡0 , 𝑇] ∈ ℝ1
определена переменная 𝑡 имеющая физический смысл времени.
𝑇
Рассмотрим случайный процесс 𝐱(𝑡) = (𝑥1 (𝑡), … , 𝑥 𝑑 (𝑡)) , где 𝐱(𝑡) принадлежит функциональному
пространству L2 (Ω) с нормой ‖ ∙ ‖. Будем считать, что случайный процесс 𝐱(𝑡) является решением для
стохастического дифференциального уравнения (СДУ) в форме Ито [1,29,30]:
𝐱(𝑡) = 𝐟(𝑡, 𝐱(𝑡))d𝑡 + 𝐆(𝑡, 𝐱(𝑡))d𝐖,
где 𝐖 = (𝑊 𝟏 , … , 𝑊 𝑚 )𝑇 — многомерный винеровский процесс, называемый ведущим (driver) процессом
СДУ. Функция 𝐟: [𝑡0 , 𝑇] × ℝ𝑑 → ℝ𝑑 называется вектором сноса, а матричнозначная функция 𝐆: [𝑡0 , 𝑇] ×
𝑇
ℝ𝑑 × ℝ𝑚 → ℝ𝑑 × ℝ𝑚 — матрицей диффузии. Кроме того, 𝐟(𝑡, 𝐱(𝐭)) = (𝑓1 (𝑡, 𝐱(𝐭)), … , 𝑓 𝑑 (𝑡, 𝐱(𝐭))) , а
матрица 𝐆 имеет вид:
𝑔11 (𝑡, 𝐱) ⋯ 𝑔𝑚 1
(𝑡, 𝐱)
𝐆(𝑡, 𝐱(𝑡)) = [ ⋮ ⋱ ⋮ ]
𝑔1𝑑 (𝑡, 𝐱) ⋯ 𝑔𝑚 𝑑
(𝑡, 𝐱)
То же уравнение можно переписать в индексном виде:
𝑚

𝑥 𝛼 (𝑡) = 𝑓 𝛼 (𝑡, 𝑥 𝛾 (𝑡))d𝑡 + ∑ 𝑔𝛽𝛼 (𝑡, 𝑥 𝛾 (𝑡))d𝑊𝛽 ,
𝛽=1
где 𝛼, 𝛾 = 1, … , 𝑑, 𝛽 = 1, … , 𝑚, а 𝑓 𝛼 (𝑡, 𝑥 𝛾 (𝑡)) = 𝑓 𝛼 (𝑡, 𝑥1 (𝑡), … , 𝑥 𝑑 (𝑡)).
На отрезке [𝑡0 , 𝑇] введем сетку 𝑡0 < 𝑡1 < ⋯ < 𝑡𝑁 = 𝑇 с шагом ℎ𝑛 = 𝑡𝑛+1 − 𝑡𝑛 , где 𝑛 = 0, … , 𝑁 − 1 и
максимальный шаг сетки ℎ = max{ℎ𝑛−1 }1𝑁 . Далее будем полагать, что сетка равномерная, то ℎ𝑛 = ℎ =
const. 𝐱𝑛 — сеточная функция, аппроксимирующая случайный процесс 𝐱(𝑡) так что 𝐱0 = 𝐱(𝑡0 ), 𝐱𝑛 ≈ 𝐱(𝑡𝑛 ),
∀ 𝑛 = 1, … , 𝑁.
Вычисление и аппроксимация кратных интегралов Ито
В общем случае для конструирования численных схем, порядок сильной сходимости которых был бы
больше чем 0,5, необходимо включение в формулы этих схем однократных и двукратных интегралов Ито:
𝑡𝑛+1 𝑡𝑛+1 𝜏1

𝐼𝑖 (𝑡𝑛 , 𝑡𝑛+1 ) = 𝐼𝑖 (ℎ𝑛 ) = ∫ d𝑊 , 𝑖
𝐼𝑖𝑗 (𝑡𝑛 , 𝑡𝑛+1 ) = 𝐼𝑖𝑗 (ℎ𝑛 ) = ∫ ∫ d𝑊 𝑖 (𝜏2 )d𝑊𝑗 (𝜏1 ),
𝑡𝑛 𝑡𝑛 𝑡𝑛
где 𝑖, 𝑗 = 1, … , 𝑚 и 𝑊𝑖 — компоненты многомерного винеровского процесса.
Возникает задача выражения этих интегралов через приращения Δ𝑊𝑛𝑖 = 𝑊 𝑖 (𝑡𝑛+1 ) − 𝑊 𝑖 (𝑡𝑛 ). В случае
однократного интеграла это можно сделать для любого индекса 𝑖: 𝐼𝑖 (ℎ𝑛 ) = 𝑊𝑛𝑖 В случае же двукратного
интеграла 𝐼𝑖𝑗 (ℎ𝑛 ) точная формула имеет место лишь при 𝑖 = 𝑗:
1
𝐼𝑖𝑖 (ℎ𝑛 ) = ((Δ𝑊𝑛𝑖 )2 − Δ𝑡𝑛 ),
2
в остальных же случаях при 𝑖 ≠ 𝑗 выразить 𝐼𝑖𝑗 (ℎ𝑛 ) через приращения Δ𝑊𝑛𝛼 и Δ𝑡𝑛 в конечном виде не
представляется возможным, поэтому остается лишь использовать численную аппроксимацию.

31
В книге [21] приведены следующие формулы для аппроксимации двукратного интеграла Ито 𝐼𝑖𝑗 :
𝑗
Δ𝑊𝑛𝑖 Δ𝑊𝑛 − ℎ𝑛 𝛿 𝑖𝑗
𝐼𝑖𝑗 (ℎ𝑛 ) = + 𝐴𝑖𝑗 (ℎ𝑛 ),
∞
2
ℎ 1 𝑗 𝑗 𝑗
𝐴𝑖𝑗 (ℎ𝑛 ) = ∑ (𝑉𝑘𝑖 (𝑈𝑘 + √2/ℎ𝑛 Δ𝑊𝑛 ) − 𝑉𝑘 (𝑈𝑘𝑖 + √2/ℎ𝑛 Δ𝑊𝑛𝑖 ) ) ,
2𝜋 𝑘
𝑘=1
Δ𝑊𝑛𝑖 ~𝒩(0, ℎ), 𝑉𝑘𝑖 ~𝒩(0, ℎ), 𝑈𝑘𝑖 ~𝒩(0, ℎ), 𝑖, 𝑗 = 1, … , 𝑚; 𝑘 = 1, … , ∞; 𝑗 = 1, … , 𝑁. Из формул видно, что в
случае 𝑖 = 𝑗 получаем конечное выражение для 𝐼𝑖𝑗 , но в случае 𝑖 ≠ 𝑗 приходится суммировать
бесконечный ряд 𝐴𝑖𝑗 . Данный алгоритм дает ошибку аппроксимации порядка 𝑂(ℎ2 /𝑛), где 𝑛 — число
оставленных слагаемых бесконечного ряда 𝐴𝑖𝑗 .
В статье [26] предложен матричный вид аппроксимирующих формул. Пусть 𝟏𝑚×𝑚 , 𝟎𝑚×𝑚 — единичная
и нулевая матрицы 𝑚 × 𝑚, тогда
Δ𝐖𝑛 Δ𝐖𝑛 − ℎ𝑛 𝟏𝑚×𝑚
𝐈(ℎ𝑛 ) = + 𝐀(ℎ𝑛 ),
∞
2
ℎ 1 𝑇
𝐀(ℎ𝑛 ) = ∑ (𝐕𝑘 (𝐔𝑘 + √2/ℎ𝑛 Δ𝐖𝑛 ) − (𝐔𝑘 + √2/ℎ𝑛 Δ𝐖𝑛 )𝐕𝑘𝑇 ),
2𝜋 𝑘
𝑘=1
где Δ𝐖𝑛 , 𝐕𝑘 , 𝐔𝑘 — независимые нормально распределенные многомерные случайные величины:
Δ𝐖𝑛 = (Δ𝑊𝑛1, Δ𝑊𝑛2 , … , Δ𝑊𝑛𝑚 )𝑇 ~𝒩(𝟎𝑚×𝑚 , ℎ, 𝟏𝑚×𝑚 ),
𝐕𝑘 = (𝑉𝑘 , 𝑉𝑘 , … , 𝑉𝑘𝑚 )𝑇 ~𝒩(𝟎𝑚×𝑚 , 𝟏𝑚×𝑚 ), 𝐔𝑘 = (𝑈𝑘1 , 𝑈𝑘2 , … , 𝑈𝑘𝑚 )𝑇 ~𝒩(𝟎𝑚×𝑚 , 𝟏𝑚×𝑚 ),
1 2

Слабая и сильная сходимости аппроксимирующей функции
Необходимо определить критерий точности аппроксимации процесса 𝐱(t) функцией 𝐱𝐧 . Таких
критериев принято выделять два: слабый и сильный [25, 32, 33].
Последовательность сеточных аппроксимирующих функций 𝐱1 , … 𝐱𝑁 сходится в сильном смысле с
порядком 𝑝 к решению 𝐱(𝑡) СДУ в момент времени 𝑇 если существует константа 𝐶 > 0 и 𝛿0 > 0 такая что
∀ℎ ∈ (0, 𝛿0 ] выполняется условие
𝔼(‖𝐱(𝑇) − 𝐱𝑁 ‖) ≤ 𝐶ℎ𝑝 .
Последовательность сеточных аппроксимирующих функций 𝐱1 , … 𝐱𝑁 сходится в слабом смысле с
порядком 𝑝 к решению 𝐱(t) СДУ в момент времени 𝑇 если существует константа 𝐶𝐹 > 0 и 𝛿0 > 0 такая что
∀ℎ ∈ (0, 𝛿0 ] выполняется условие
|𝔼[𝐹(𝐱(𝑇))] − 𝔼[𝐹(𝐱𝒏 )]| ≤ 𝐶ℎ𝑝 .
2(𝑝+1)
Здесь 𝐹 ∈ 𝐶𝑃 (ℝ, ℝ𝑑 ) — непрерывный дифференцируемый функционал с полиномиальным ростом.
Если матрица 𝐆 обращается в ноль, то условие сильной сходимости равносильно условию сходимости
для детерминированного случая. В отличие от детерминированного случая порядок сильной сходимости
не обязательно является натуральным числом и может принимать дробно-рациональные значения.
Важно отметить, что выбор типа сходимости зависит от решаемой задачи. Увеличение порядка
строгой сходимости ведет к увеличению точности аппроксимации траекторий решения 𝐱(𝑡) СДУ. Если же
требуется вычислить, например, момент случайного процесса 𝐱(𝑡) или обобщенный функционал вида
𝔼[𝐹(𝐱(𝑡))], то следует повышать порядок слабой сходимости.
Стохастические методы Рунге–Кутты
Перейдем к рассмотрению стохастичесских методов Рунге-Кутты. Мы приведем формулы для трех
численных схем, записав для каждой из них формулы, таблицу Бутчера в общем виде и несколько таблиц
с конкретным набором коэффициентов.
Простейшим численным методом для решения как скалярных уравнений, так и систем СДУ, является
метод Эйлера–Маруямы, названный так в честь Гиширо Маруямы (Gisiro Maruyama), который
распространил классический метод Эйлера для ОДУ на случай СДУ [27]. Метод легко обобщается на случай
многомерного винеровского процесса:
𝑥0𝛼 = 𝑥 𝛼 (𝑡0 ),
𝑑
𝛾 𝛽
𝛼 𝛼 𝛼 (𝑡 𝛼 )ℎ
𝑥𝑛+1 = 𝑥𝑛 + 𝑓 𝑛 , 𝑥𝑛 𝑛 + ∑ 𝐺𝛽𝛼 (𝑡𝑛 , 𝑥𝑛 )Δ𝑊𝑛 .
𝛾=1
Как видно из формул на каждом шаге для вычисления следующего значения аппроксимирующей
𝛽
функции требуется лишь соответствующее данному шагу приращение процесса 𝑊𝑛 Метод имеет
сильный порядок (𝑝𝑑 , 𝑝𝑠 ) = (1.0,0.5) . Величина 𝑝𝑑 обозначает порядок точности детерминированной
части численного метода, то есть той точности, которую будет давать численный метод при применении
его к уравнению функцией 𝐺(𝑡, 𝑥 𝛼 (𝑡)) ≡ 0 .Величина 𝑝𝑠 обозначает порядок приближения
стохастической части уравнения.

32
В случае скалярных уравнения и винеровского процесса справедлива следующая численная схема
Рунге–Кутты сильной сходимости 𝑝 = 1.5 [13]:
𝑗 𝑗 𝑗 𝑗 𝐼10 (ℎ𝑛 )
𝑋0𝑖 = 𝑥𝑛 + ∑𝑠𝑗=1 𝐴𝑖0 𝑓(𝑡𝑛 + 𝑐0 ℎ𝑛 , 𝑋0 )ℎ𝑛 + ∑𝑠𝑗=1 𝐵0𝑗
𝑖
𝑔(𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋1 ) ,
√ℎ𝑛

𝑗 𝑗 𝑗 𝑗
𝑋1𝑖 = 𝑥𝑛 + ∑𝑠𝑗=1 𝐴𝑖1𝑗 𝑓(𝑡𝑛 + 𝑐0 ℎ𝑛 , 𝑋0 )ℎ𝑛 + ∑𝑠𝑗=1 𝐵1𝑗
𝑖
𝑔(𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋1 ) √ℎ𝑛 ),
𝐼11 (ℎ𝑛 ) 𝐼10 (ℎ𝑛 ) 𝐼111 (ℎ𝑛 )
𝑥𝑛+1 = 𝑥𝑛 + ∑𝑠𝑖=1 𝑎𝑖 𝑓(𝑡𝑛 + 𝑐 𝑖 ℎ𝑛 , 𝑋0𝑖 )ℎ𝑛 + ∑𝑠𝑖=1 (𝑏𝑖1 𝐼1 (ℎ𝑛 ) + 𝑏𝑖2 + 𝑏𝑖3 ℎ𝑛
+ 𝑏𝑖4 ℎ𝑛
) 𝑔(𝑡𝑛 + 𝑐1𝑖 , 𝑋1𝑖 ),
√ℎ𝑛

обобщенная таблица Бутчера который имеет вид

В препринте Росслера [10] для данной схемы приведены две таблицы Бутчера для метода четвертой
стадийности 𝑠 = 4.

Численную схему, реализуемую первой таблицей, обозначим как SRK1W1, а вторую как SRK2W2. Метод
SRK1W1 имеет сильный порядок (𝑝𝑑 , 𝑝𝑠 ) = (2.0,1.5), а метод SRK2W1 сильный порядок (𝑝𝑑 , 𝑝𝑠 ) =
(3.0,1.5). Еще один метод с сильным порядком 𝑝𝑠 = 1.0 приведен в книге [21] и его таблица Бутчера имеет
вид:

Для системы СДУ Ито с многомерным винеровским процессом можно построить стохастическую
численную схему Рунге-Кутты сильного порядка 𝑝𝑠 = 1.0 с использованием однократных и двукратных
интегралов Ито [13]
𝑠 𝑚 𝑠
𝑖 𝑗 𝑖 𝑗
𝑋 0𝑖𝛼 𝛼 𝛼 𝑜𝑗𝛽
= 𝑥𝑛 + ∑ 𝐴0𝑗 𝑓 (𝑡𝑛 + 𝑐0 ℎ𝑛 , 𝑋 )ℎ𝑛 + ∑ ∑ 𝐵0𝑗 𝐺𝑙𝛼 (𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋𝑙𝑗𝛽 )𝐼𝑙 (ℎ𝑛 )
𝑗=1 𝑙=1 𝑗=1
𝑠 𝑚 𝑠
𝑖 𝑗 𝑖 𝑗 𝐼𝑙 (ℎ𝑛 )
𝑋 𝑘𝑖𝛼 𝛼 𝛼 𝑜𝑗𝛽
= 𝑥𝑛 + ∑ 𝐴1𝑗 𝑓 (𝑡𝑛 + 𝑐0 ℎ𝑛 , 𝑋 )ℎ𝑛 + ∑ ∑ 𝐵1𝑗 𝐺𝑙𝛼 (𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋𝑙𝑗𝛽 ) ,
𝑗=1 𝑙=1 𝑗=1
√ℎ𝑛
𝑠 𝑚 𝑠
𝑗
𝑥𝑛+1 = 𝑥𝑛 + ∑ 𝑎𝑖 𝑓 (𝑡𝑛 + 𝑐0 ℎ𝑛 , 𝑋 )ℎ𝑛 + ∑ ∑(𝑏𝑖1 𝐼𝑘 (ℎ𝑛 ) + 𝑏𝑖2 √ℎ𝑛 )𝐺𝑘𝛼 (𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋𝑘𝑖𝛽 ),
𝛼 𝛼 𝛼 𝑖 𝑜𝑗𝛽

𝑖=1 𝑘=1 𝑖=1
где 𝑛 = 0,1, ⋯ , 𝑁 − 1; 𝑖 = 1, ⋯ , 𝑠; 𝛽, 𝑘 = 1, ⋯ , 𝑚; 𝛼 = 1, ⋯ , 𝑑 и обобщенная таблица Бутчера имеет
следующий вид

33
В препринте Росслера [21] указаны две таблицы Бутчера для метода третей стадийности 𝑠 = 3:

Метод SRK1Wm имеет сильный порядок (𝑝𝑑 , 𝑝𝑠 ) = (1.0,1.0), а метод SRK2Wm сильный порядок
(𝑝𝑑 , 𝑝𝑠 ) = (1.0,2.0).
Стохастические методы Рунге–Кутты слабой сходимости 𝑝 = 2.0 Численные методы со слабой
сходимостью хорошо аппроксимируют характеристики распределения случайного процесса 𝑥 𝛼 (𝑡) .
Слабый численный метод не нуждается в информации о траектории винеровского процесса 𝑊𝑛𝛼 и
случайные величины для этих методов могут быть сгенерированы на другом вероятностном
пространстве. Поэтому данный метод можно использовать с легко моделируемом на компьютере
распределением. Формулы данной численной схемы самые громоздкие из трех рассматриваемых нами
схем [13]:
𝑠 𝑠 𝑚
𝑗 𝑗
𝑋 0𝑖𝛼 𝛼 𝑖 𝛼 𝑜𝑗𝛽
= 𝑥𝑛 + ∑ 𝐴0𝑗 𝑓 (𝑡𝑛 + 𝑐0 ℎ𝑛 , 𝑋 )ℎ𝑛 + ∑ ∑ 𝐵0𝑗 𝑖
𝐺𝑙𝛼 (𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋𝑙𝑗𝛽 )𝐼̂𝑙 ,
𝑗=1 𝑗=1 𝑙=1

𝑠 𝑠
𝑖 𝑗 𝑖 𝑗
𝑋 𝑘𝑖𝛼 𝛼 𝛼 𝑜𝑗𝛽
= 𝑥𝑛 + ∑ 𝐴1𝑗 𝑓 (𝑡𝑛 + 𝑐0 ℎ𝑛 , 𝑋 )ℎ𝑛 + ∑ 𝐵1𝑗 𝐺𝑙𝛼 (𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋𝑙𝑗𝛽 )√ℎ𝑛 ,
𝑗=1 𝑗=1

𝑠 𝑠 𝑚
𝑗 𝑗 𝐼̂𝑘𝑙
𝑋̂ 𝑘𝑖𝛼 = 𝑥𝑛𝛼 + ∑ 𝐴𝑖2𝑗 𝑓 𝛼 (𝑡𝑛 + 𝑐0 ℎ𝑛 , 𝑋 𝑜𝑗𝛽 )ℎ𝑛 + ∑ ∑ 𝐵2𝑗
𝑖
𝐺𝑙𝛼 (𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋𝑙𝑗𝛽 ) ,
𝑗=1 𝑗=1 𝑙=1,𝑙≠𝑘
√ℎ𝑛
𝑠
𝛼
𝑥𝑛+1 = 𝑥𝑛 + ∑ 𝑎𝑖 𝑓 𝛼 (𝑡𝑛 + 𝑐1𝑖 , 𝑋𝑘𝑖𝛽 )ℎ𝑛
𝛼

𝑖=1
𝑠 𝑚 𝑠 𝑚
̂𝐼𝑘𝑘
+ ∑ ∑ (𝑏𝑖1 𝐼̂𝑘 + 𝑏𝑖2 ) 𝐺𝑘𝛼 (𝑡𝑛 + 𝑐1𝑖 ℎ𝑛 , 𝑋𝑘𝑖𝛽 ) + ∑ ∑(𝑏𝑖3 𝐼̂𝑘 + 𝑏𝑖4 √ℎ𝑛 )𝐺𝑘𝛼 (𝑡𝑛 + 𝑐2𝑖 ℎ𝑛 , 𝑋̂𝑘𝑖𝛽 )
𝑖=1 𝑘=1
√ℎ𝑛 𝑖=1 𝑘=1
где 𝑛 = 0,1, ⋯ , 𝑁 − 1; 𝑖 = 1, ⋯ , 𝑠; 𝛽, 𝑘 = 1, ⋯ , 𝑚; 𝛼 = 1, ⋯ , 𝑑.
Обобщенная таблица Бутчера имеет вид

В численной схеме слабого стахостического метода Рунге-Кутты используются следующие случайные
величины [13]:

34
1 𝑘 𝑙
(𝐼̂ 𝐼̂ − √ℎ𝑛 𝐼̃𝑘 ), 𝑘 < 𝑙,
2
1 𝑘 𝑙
𝐼̂𝑘𝑙 = (𝐼̂ 𝐼̂ − √ℎ𝑛 𝐼̃𝑙 ), 𝑙 < 𝑘,
2
1 2
̂𝑘
{ 2 ((𝐼 ) − ℎ𝑛 ) , 𝑘 = 𝑙,
Случайная величина 𝐼̂𝑘𝑙 имеет трехточечное распределение. Это означает, что 𝐼̂𝑘𝑙 может принимать
три значения: {−√3ℎ𝑛 ,0, √3ℎ𝑛 } с вероятностями 1/6, 2/3 и 1/6 соответственно. В свою очередь
случайная величина 𝐼̂𝑘 имеет двухточечное распределение и принимает два значения {−√ℎ𝑛 , √ℎ𝑛 } с
вероятностями 1/2 и 1/2.
Анализ трудностей реализации стохастических численных методов Рунге-Кутты
Как видно уже из формул, стохастические методы Рунге-Кутты значительно сложнее своих
классических аналогов. Кроме громоздкости формул, можно выделить еще следующие факторы,
усложняющие реализацию стохастических методов в программном виде, а также их применение для
численного решения СДУ.
• При выборе конкретного метода надо учитывать, какой тип сходимости необходимо обеспечить
для данной конкретной задачи, а также какое из стохастических уравнений необходимо решать
— в форме Ито или в форме Стратоновича. Это увеличивает количество алгоритмов, которые
нужно запрограммировать.
• Для методов с сильной сходимостью большей единицы на каждом шаге необходимо решать
ресурсоемкую задачу по аппроксимации двукратных стохастических интегралов.
• В численной схеме присутствуют не только матрицы и векторы, но и тензоры (четырехмерные
массивы), с которыми необходимо совершать операцию свертки по нескольким индексам.
Реализация свертки через суммирование с помощью обычных циклов приводит к существенному
падению производительности.
• Для использования слабых методов необходимо применять метод Монте Карло, проводя
несколько серий по множеству испытаний в каждой. Так как метод Монте-Карло сходится
приблизительно как 1/√𝑁, где 𝑁 — число испытаний, то для достижения точности хотя бы 10-3,
необходимо провести минимум 106 испытаний.
Наиболее существенное падение производительности происходит при реализации универсального
алгоритма, то есть такой программы, которая может произвести расчет, используя произвольную
таблицу коэффициентов. В этом случае приходится использовать большое количество вложенных циклов,
для того чтобы организовать суммирование. Наличие в схемах для систем СДУ двойных сумм и сложной
комбинации индексов в множителях, находящихся под знаком этих сумм, еще более усложняет задачу и
число вложенных циклов вырастает до шести.
Кроме этих специфических особенностей стохастических численных методов, стоит упомянуть также
несколько причин падения производительности, присущих также и классическим схемам, которые в
стохастическом случае также играют роль. Очевидным способом хранения коэффициентов методов
является использование массивов. Однако у явных методов, которые мы рассматриваем, матрица
является нижне-диагональной и хранение ее в виде двумерного массива приводит к тому, что более чем
половина выделенной для массива памяти тратится на хранение нулей.
Если изучить исходные коды популярных подпрограмм, реализующих классические явные вложенные
методы Рунге—Кутты [28], то можно обнаружит, что в этих программах для хранения коэффициентов
метода используется набор именованных констант, а не массивов. Это вызвано также и тем, что операции
со скалярными величинами в большинстве языков программирования проводятся быстрее, чем операции
с массивами.
При сохранении требования универсальности создаваемого кода и вместе с тем желание увеличить
скорости вычислений и уменьшить расход памяти, привели нас к решению использовать автоматическую
генерацию кода по одному шаблону для каждого отдельного метода.
Кроме выигрыша в производительности, автоматическая генерации кода, позволяет добавлять или
изменять все функции за раз путем редактирования одного лишь шаблона, а не каждой отдельной
функции. Это позволяет как уменьшить количество ошибок, так и генерировать различные варианты
функций для разных целей.
Краткое описание автоматической генерации кода
В качестве языка для написания генератора кода используется язык Python. Исходный код и
документация созданных в нашем коллективе программ, доступны по ссылке

35
https://bitbucket.org/mngev/sde_num_generation. Данный репозиторий содержит модуль stochastic, в
котором реализован винеровский случайный процесс и численные методы типа Рунге-Кутты для
скалярных и многомерных стохастических дифференциальных уравнений в форме Ито. Большая часть
кода для этого модуля генерируется набором скриптов из директории generator. Вручную написан класс,
реализующий многомерный винеровский процес и методы Эйлера-Маруямы.
Для работы генератора кода мы использовали библиотеку для обработки шаблонов Jinja2 [29]. Данный
шаблонизатор (template engine) разрабатывался изначально для генерации HTML страниц, однако он
обладает очень гибким синтаксисом и может использоваться как универсальное средство для генерации
текстовых файлов любого вида, в том числе и исходных кодов на любых языках программирования. Кроме
Jinja2 мы также использовали библиотеку NumPy для работы с массивами и ускорения вычислений.
Кроме перечисленных выше двух внешних библиотек, был использован стандартный модуль fraction,
который позволяет задать коэффициенты метода в виде рациональных дробей, а потом уже
преобразовывать их в вещественный вид с нужным порядком точности, а также модуль typing для
аннотации типов аргументов функций.
Шаблоны для генерации представляют собой файлы с исходным кодом на языке Python со вставками
специальных команд шаблонизатора Jinja2. Информация о коэффициентах методов хранится отдельно, в
структурированном виде в формате JSON. Это позволяет легко добавлять новые методы и изменять
старые путем редактирования и JSON файлов. В настоящее время генерируются методы с
коэффициентами, представленными в работах [12,13] и [21].
В качестве языка для уже сгенерированных функций используется сам python с активным
применением библиотеки NumPy, которая позволяет получить приемлемую производительность. Однако
генерируемый код легко можно переформатировать так, чтобы он соответствовал синтаксису любого
другого языка программирования. Мы планируем доработать программу, для генерации кода на языке
Julia (julialang.org). Данный язык был представлен в 2012 году и изначально ориентирован на научные
вычисления. В настоящее время он интенсивно развивается и набирает популярность. На сегодняшний
день актуальной является версия 0.6.1. Julia позволяет получить производительность, сравнимую с C++ и
Fortran, но при этом является динамическим языком с возможностью интерактивной командной строки
(REPL) наподобие IPython и может интегрироваться в интерактивную среду Jupyter. Текущая версия
библиотеки превосходит описанную авторами в статье [1]. Применение автогенерации позволю
избавится от большого количества вложенных циклов, выделения лишней памяти для хранения нулевых
коэффициентов, а также значительно упростило код.
Тестирование созданных программ
Для тестирования созданных программ мы использовали стохастические численные методы с
известными аналитическими решениями. Большой список скалярных СДУ с аналитическими решениями
представлен в книге [21], а для многомерного случая мы использовали уравнение Блека-Шоулза в форме
Ито, представленное в книге [36]. Для проверки корректности работы генератора, был создан шаблон для
генерации формул в формате LaTeX, которые затем проверялись на корректность визуально.
В качестве одномерного уравнения для тестирования методов с сильной сходимостью выбрано
уравнение логарифмического блуждания (применяется в финансовой математике для моделирования
рынков). Оно имеет следующий вид:
𝑑𝑥(𝑡) = 𝜇𝑥(𝑡)𝑑𝑡 + 𝜎𝑥(𝑡)𝑑𝑊(𝑡), 𝑥(0) = 𝑥0,
где W(t) — скалярный процесс Винера, а x0 —начальное значение. Его аналитическое решение
выражается через процесс Винера в конечном виде:
𝑥(𝑡) = 𝑥0 exp((𝜇 − 0.5𝜎 2 ) + 𝜎𝑊(𝑡)).

Рис. 1 Абсолютная локальная погрешность скалярных стохастических методов с сильной сходимостью

36
При вычислениях будем использовать коэффициенты µ = 2, σ = 1 и x0 = 1. Результаты тестирования
погрешности представлены на рисунке 1. Аббревиатурой EM обозначен метод Эйлера-Маруямы,
коэффициенты трех других методов представлены выше в пункте «Стохастические методы Рунге–Кутты».
Шаг дискретизации h равен 0.001. Результаты согласуются с теоретической оценкой сильной
погрешности.
Уравнение Блека-Шоулза в форме Ито имеет следующие вектор сдвига и матрицы диффузии:
𝑎1 𝑥1 𝑏1 𝑥1 0
𝑓(𝑥1 , 𝑥2 ) = (𝑎 𝑥 ) , 𝐺(𝑥1 , 𝑥2 ) = [ ],
2 2 𝑏2𝜌𝑥2 𝑏2 √1 − 𝜌 2
Значения параметров выберем следующие a 1 = a2 = 0.1, b1 = b2 = 0.2, ρ = 0.8. Начальные значения: x1(0)
= x1(0) = 1. Аналитическое решение имеет следующий вид:
𝑥1 (𝑡) = 𝑥01 𝑒𝑥𝑝((𝑎1 − 0.5𝑏12 )𝑡 + 𝑏1 𝑊1 (𝑡))
{ 2 .
𝑥 (𝑡) = 𝑥0 𝑒𝑥𝑝((𝑎2 − 0.5𝑏22 )𝑡 + 𝑏2 (𝜌𝑊 1(𝑡) + √1 − 𝜌𝑊 2 (𝑡)))
2

Погрешность для метода, обозначенного нами выше как SRK1Wm изображена на рисунке 2. Для
второго метода (SRK2Wm) погрешность вычислений аналогична.

Рис. 2 Абсолютная локальная погрешность векторного стохастического метода SRK1Wm
Проведенные тесты позволяют убедится в корректности генерируемого кода. С помощью
«магической» команды %timeit оболочки iPython мы также оценили производительность векторных
численных стохастических методов. Для случая уравнения Блека-Шоулза время работы функции
strong_srk1w2 (реализует метод SRK2Wm для размерности 2) на сравнительно медленном процессоре
Intel i3-4100 составило около секунды, что вполне приемлемо в случае разовых вычислений.
Заключение
В данной работе рассмотрены стохастические численные схемы с порядком сходимости выше 0.5.
Показано, что такие методы существенно сложнее эквивалентных численных методов для систем
обыкновенных дифференциальных уравнений. Выделены их особенности, делающие эффективную
программную реализацию таких методов не тривиальной задачей. В заключительной части статьи
обсуждается подход, основанный на автоматической генерации программного кода, который позволяет
получить эффективную реализацию рассматриваемых методов, давая возможность использовать любую
таблицу коэффициентов. Дано краткое описание программы, созданной авторами и ссылка на
репозиторий с исходными кодами.
Благодарности
Публикация подготовлена при поддержке программы РУДН «5-100» и при финансовой поддержке
РФФИ в рамках научных проектов РФФИ № 15-07- 08795, 16-07-00556.
Литература
1. Gevorkyan M. N., Velieva T. R., Korolkova A. V., Kulyabov D.S., Sevastyanov L.A. Stochastic Runge–Kutta Software Package for
Stochastic Differential Equations // Dependability Engineering and Complex Systems. –– Springer International Publishing, 2016. –
– Vol. 470.––P. 169–179.–– 1606.06604
2. Python Reference Manual: Rep. ; Executor: Guido Rossum. — Amsterdam, The Netherlands, The Netherlands: 1995.
3. Jones Eric, Oliphant Travis, Peterson Pearu et al. SciPy: Open source scientific tools for Python. — 2001–. — [Online; accessed
19.01.2017]. Access mode: http://www.scipy.org/
4. Bachelier L. Théorie de la spéculation // Annales Scientifiques de l’École Normale Supérieure.—1900. —Vol. 3, no. 17.—P. 21–86.
5. Лукшин А. В., Смирнов С. Н. Численные методы решения стохастических дифференциальных уравнений //
Вычислительные алгоритмы и методы. –– 1990. –– Т. 2, № 11. ––С. 108–121.

37
6. Ерешко А. Ф., Филатова Д. В. Анализ явных численных методов решения стохастических дифференциальных уравнений //
Труды ИСА РАН. Динамика неоднородных систем.––2008.––Т. 32, № 2.––С. 164–173.
7. Milstein G.N.Approximate Integration of Stochastic Differential Equations //TheoryProbab.Appl. —1974.—no. 19.—P. 557–562.
8. Milstein G. N. A Method of Second-Order Accuracy Integration of Stochastic Differential Equations // Theory Probab. Appl.— 1979.—
no. 23. — P. 396–401.
9. Milstein G. N. Weak Approximation of Solutions of Systems of Stochastic Differential Equations // Theory Probab. Appl.—1986.—
no. 30.— P. 750–766.
10. Rößler Andreas. Runge-Kutta Methods for the Numerical Solution of Stochastic Differential Equations : Ph.D. thesis / Andreas Rößler ;
Technischen Universität Darmstadt. — Darmstadt, 2003.— februar.
11. Hairer E., Nørsett S.P., Wanner G. Solving Ordinary Differential Equations I. — 2edition. — Berlin : Springer, 2008.—ISBN: 978-3-
540-56670-0 .
12. Debrabant K., Rößler A. Continuous weak approximation for stochastic differential equations // Journal of Computational and
Applied Mathematics. — 2008. — no. 214. — P. 259–273.
13. Debrabant K., Rößler A. Classification of Stochastic Runge–Kutta Methods for the Weak Approximation of Stochastic Differential
Equations / Technische Universität Darmstadt, Fachbereich Mathematik.—2013.— Mar.— arXiv:1303.4510v1.
14. Rößler A. Strong and Weak Approximation Methods for Stochastic Differential Equations – Some Recent Developments /
Department Mathematik. Schwerpunkt Mathematische Statistik und Stochastische Prozesse.—2010
15. Burrage K., Burrage P. M. High strong order explicit Runge-Kutta methods for stochastic ordinary differential equations // Appl.
Numer. Math. — 1996. — no. 22. — P. 81–101.
16. Burrage K., Burrage P. M., Belward J. A. A bound on the maximum strong order of stochastic Runge-Kutta methods for stochastic
ordinary differential equations. // BIT. — 1997. — no. 37. — P. 771–780.
17. Burrage K., Burrage P. M. General order conditions for stochastic Runge-Kutta methods for both commuting and non-commuting
stochastic ordinary differential equation systems // Appl. Numer. Math.—1998.—no. 28.— P. 161–177.
18. Burrage P. M. Runge-Kutta Methods for Stochastic Differential Equations: Ph.D. thesis /P. M. Burrage ; University of Qeensland.—
Australia, 1999.
19. Burrage K., Burrage P. M. Order conditions of stochastic Runge-Kutta methods by B-series // SIAM J. Numer. Anal.—2000.— no.
38.— P. 1626–1646.
20. Soheili A. R., Namjoo M. Strong approximation of stochastic differential equations with Runge–Kutta methods // World Journal of
Modelling and Simulation. — 2008. — Vol. 4, no. 2. — P. 83–93.
21. Kloeden Peter E., Platen Eckhard. Numerical Solution of Stochastic Differential Equations.— 2 edition.—Berlin Heidelberg New York :
Springer, 1995. — ISBN: 3-540-54062-8
22. Komori Y.,Mitsuri T.Stable ROW-Type Weak Scheme for Stochastic Differential Equations// RIMS Kokyuroku.—1995.—no. 932.—
P. 29–45.
23. Mackevičius V. Second-order weak approximations for stratonovich stochastic differential equations // Lithuanian Mathematical
Journal. — 1994. — Vol. 34, no. 2. — P. 183–200. — Access mode: http://dx.doi.org/10.1007/BF02333416
24. Tocino A., Ardanuy R. Runge–Kutta methods for numerical solution of stochastic differential equations // Journal of Computational
and Applied Mathematics. — 2002. — No. 138. — P. 219–241.
25. Øksendal Bernt. Stochastic differential equations. An introduction with applications. — 6 edition.—Berlin Heidelberg New York:
Springer, 2003. — ISBN: 3-540-04758-1
26. Wiktorsson M. Joint characteristic function and simultaneous simulation of iterated Itô integrals for multiple independent Brownian
motions // The Annals of Applied Probability. —2001. — Vol. 11, no. 2.—P. 470–487.
27. Maruyama G. Continuous Markov processes and stochastic equations // Rendiconti del Circolo Matematico.—1955.—no. 4. —P. 48–
90.
28. Burrage Kevin, Burrage Pamela M. Low rank Runge–Kutta methods, symplecticity and stochastic Hamiltonian problems with
additive noise // Journal of Computational and Applied Mathematics. — 2012. — Vol. 236, no. 16. — P. 3920–3930. — Access mode:
http://www.sciencedirect.com/science/article/pii/S0377042712001240.
29. Ma Qiang, Ding Xiaohua. Stochastic symplectic partitioned Runge–Kutta methods for stochastic Hamiltonian systems with
multiplicative noise // Applied Mathematics and Computation. — 2015. — Vol. 252, no. Supplement C. — P. 520–534. — Access
mode: http://www.sciencedirect.com/science/article/pii/S0096300314016993.
30. Stochastic symplectic Runge–Kutta methods for the strong approximation of Hamiltonian systems with additive noise / Weien Zhou,
Jingjing Zhang, Jialin Hong, SongheSong// Journal of Computational and Applied Mathematics. — 2017. — Vol. 325, no. Supplement
C. — P. 134– 148. — Access mode: http://www.sciencedirect.com/science/article/pii/S0377042717302285.
31. Amiri Sadegh, Hosseini S. Mohammad. Stochastic Runge–Kutta Rosenbrock type methods for SDE systems // Applied Numerical
Mathematics. — 2017. — Vol. 115, no. Supplement C. — P. 1–15. — Access mode:
http://www.sciencedirect.com/science/article/pii/S0168927416302355.
32. Булинский А.В., Ширяев А. Н. Теория случайных процессов. ––5 изд.––Москва: ФИЗМАТЛИТ, 2003.––ISBN: 5-9221-0335-0
33. Кузнецов Д. Ф. Стохастические Дифференциальные Уравнения: Теория и Практика Численного Решения. –– 4 изд. –– Санкт-
Петербург : Издательство Политехнического Университета, 2010.
34. Fortran and Matlab Codes. — Access mode: https://www.unige.ch/~hairer/software.html.
35. Jinja2 official site.—Access mode: http://http://jinja.pocoo.org.
36. Platen Eckhard, Bruti-Liberati Nicola. Numerical Solution of Stochastic Differential Equations with Jumps in Finance . — 2 edition.
— Heidelberg Dordrecht London New York : Springer, 2010.—ISBN: 978-3-642-12057-2

References
1. Gevorkyan, M. N.; Velieva, T. R.; Korolkova, A. V.; Kulyabov, D. S. & Sevastyanov, L. A. (2016), Stocha stic Runge–Kutta Software
Package for Stochastic Differential Equations'Dependability Engineering and Complex Systems', Springer International Publishing, ,
pp. 169--179.
2. Rossum, G. (1995), 'Python Reference Manual', CWI (Centre for Mathematics and Computer Science), Amsterdam, The Netherlands,
The Netherlands.
3. Jones, E.; Oliphant, T.; Peterson, P. & others (2001), 'SciPy: Open source scientific tools for Python', [Online; accessed 19.01.2017].
4. Bachelier, L. (1900), 'Théorie de la spéculation', Annales Scientifiques de l’École Normale Supérieure 3(17), 1986.
5. Lukshin A., V. & Smirnov S., N. (1990), 'Chislennye metody reshenija stohasticheskih differencial'nyh uravnenij', Vychislitel'nye
algoritmy i metody 2(11), 108–121.
6. Ereshko A., F. & Filatova D., V. (2008), 'Analiz javnyh chislennyh metodov reshenija stohasticheskih differencial'nyh uravnenij',

38
Trudy ISA RAN. Dinamika neodnorodnyh sistem 32(2), 164–173.
7. Milstein, G. N. (1986), 'Weak Approximation of Solutions of Systems of Stochastic Differential Equations', Theory Probab. Appl.(30),
750--766.
8. Milstein, G. N. (1979), 'A Method of Second-Order Accuracy Integration of Stochastic Differential Equations', Theory Probab.
Appl.(23), 396--401.
9. Milstein, G. N. (1974), 'Approximate Integration of Stochastic Differential Equations', Theory Probab. Appl.(19), 557–562.
10. Rößsler, A. (2003), 'Runge-Kutta Methods for the Numerical Solution of Stochastic Differential Equations', PhD thesis, Technischen
Universität Darmstadt.
11. Hairer, E.; Nørsett, S. P. & G.Wanner (2008), Solving Ordinary Differential Equations I, Springer, Berlin.
12. Debrabant, K. & Rößsler, A. (2013), 'Classification of Stochastic Runge–Kutta Methods for the Weak Approximation of Stochastic
Differential Equations', Technische Universität Darmstadt, Fachbereich Mathematik, arXiv:1303.4510v1.
13. Debrabant, K. & Rößsler, A. (2008), 'Continuous weak approximation for stochastic differential equations', Journal of Computational
and Applied Mathematics(214), 259--273.
14. Burrage, K. & Burrage, P. M. (2000), 'Order conditions of stochastic Runge-Kutta methods by B-series', SIAM J. Numer. Anal.(38),
1626--1646.
15. Rößsler, A. (2010), 'Strong and Weak Approximation Methods for Stochastic Differential Equations – Some Recent Developments',
Department Mathematik. Schwerpunkt Mathematische Statistik und Stochastische Prozesse.
16. Burrage, K. & Burrage, P. M. (1998), 'General order conditions for stochastic Runge-Kutta methods for both commuting and non-
commuting stochastic ordinary differential equation systems', Appl. Numer. Math.(28), 161--177.
17. Burrage, K. & Burrage, P. M. (1996), 'High strong order explicit Runge-Kutta methods for stochastic ordinary differential equations',
Appl. Numer. Math.(22), 81--101.
18. Burrage, K.; Burrage, P. M. & Belward, J. A. (1997), 'A bound on the maximum strong order of stochastic Runge-Kutta methods for
stochastic ordinary differential equations.', BIT(37), 771-780.
19. Burrage, P. M. (1999), 'Runge-Kutta Methods for Stochastic Differential Equations', PhD thesis, University of Qeensland.
20. Soheili, A. R. & Namjoo, M. (2008), 'Strong approximation of stochastic differential equations with Runge–Kutta methods', World
Journal of Modelling and Simulation 4(2), 83--93.
21. Kloeden, P. E. & Platen, E. (1995), Numerical Solution of Stochastic Differential Equations, Springer, Berlin Heidelberg New York.
22. Komori, Y. & Mitsuri, T. (1995), 'Stable ROW-Type Weak Scheme for Stochastic Differential Equations', RIMS Kokyuroku(932), 29--
45.
23. Mackevičius, V. (1994), 'Second-order weak approximations for stratonovich stochastic differential equations', Lithuanian
Mathematical Journal 34(2), 183-200.
24. Tocino, A. & Ardanuy, R. (2002), 'Runge–Kutta methods for numerical solution of stochastic differential equations', Journal of
Computational and Applied Mathematics(138), 219--241.
25. Øksendal, B. (2003), Stochastic differential equations. An introduction with applications, Springer, Berlin Heidelberg New York.
26. Wiktorsson, M. (2001), 'Joint characteristic function and simultaneous simulation of iterated Itô integrals for multiple independent
Brownian motions', The Annals of Applied Probability 11(2), 470-487.
27. Maruyama G. Continuous Markov processes and stochastic equations // Rendiconti del Circolo Matematico.—1955.—no. 4. —P. 48–
90.
28. Burrage Kevin, Burrage Pamela M. Low rank Runge–Kutta methods, symplecticity and stochastic Hamiltonian problems with
additive noise // Journal of Computational and Applied Mathematics. — 2012. — Vol. 236, no. 16. — P. 3920–3930. — Access mode:
http://www.sciencedirect.com/science/article/pii/S0377042712001240.
29. Ma Qiang, Ding Xiaohua. Stochastic symplectic partitioned Runge–Kutta methods for stochastic Hamiltonian systems with
multiplicative noise // Applied Mathematics and Computation. — 2015. — Vol. 252, no. Supplement C. — P. 520–534. — Access
mode: http://www.sciencedirect.com/science/article/pii/S0096300314016993.
30. Stochastic symplectic Runge–Kutta methods for the strong approximation of Hamiltonian systems with additive noise / Weien Zhou,
Jingjing Zhang, Jialin Hong, SongheSong// Journal of Computational and Applied Mathematics. — 2017. — Vol. 325, no. Supplement
C. — P. 134– 148. — Access mode: http://www.sciencedirect.com/science/article/pii/S0377042717302285.
31. Amiri Sadegh, Hosseini S. Mohammad. Stochastic Runge–Kutta Rosenbrock type methods for SDE systems // Applied Numerical
Mathematics. — 2017. — Vol. 115, no. Supplement C. — P. 1–15. — Access mode:
http://www.sciencedirect.com/science/article/pii/S0168927416302355.
32. Bulinskij A.V., Shirjaev A. N. Teorija sluchajnyh processov. ––5 izd. ––Moskva: FIZMATLIT, 2003.––ISBN: 5-9221-0335-0
33. Kuznecov D. F. Stohasticheskie Differencial'nye Uravnenija: Teorija i Praktika Chislennogo Reshenija. –– 4 izd. –– Sankt-Peterburg:
Izdatel'stvo Politehnicheskogo Universiteta, 2010.
34. Maruyama, G. (1955), 'Continuous Markov processes and stochastic equations', Rendiconti del Circolo Matematico (4), 48--90.
35. Fortran and Matlab Codes http://jinja.pocoo.org
36. Platen Eckhard, Bruti-Liberati Nicola. Numerical Solution of Stochastic Differential Equations with Jumps in Finance. — 2 edition. —
Heidelberg Dordrecht London New York: Springer, 2010.—ISBN: 978-3-642-12057-2

Об авторах:
Кулябов Дмитрий Сергеевич, кандидат физико-математических наук, доцент, доцент кафедры
прикладной информатики и теории вероятностей, Российский университет дружбы народов;
сотрудник Лаборатории информационных технологий, Объединенный институт ядерных
исследований, kulyabov_ds@rudn.university
Геворкян Мигран Нельсонович, кандидат физико-математических наук, доцент кафедры прикладной
информатики и теории вероятностей, Российский университет дружбы народов,
gevorkyan_mn@rudn.university
Демидова Анастасия Вячеславовна, кандидат физико-математических наук, старший преподаватель
кафедры прикладной информатики и теории вероятностей, Российский университет дружбы
народов, demidova_av@rudn.university
Королькова Анна Владиславовна, кандидат физико-математических наук, доцент, доцент кафедры
прикладной информатики и теории вероятностей, Российский университет дружбы народов,

39
korolkova_av@rudn.university
Севастьянов Леонид Антонович, доктор физико-математических наук, профессор, профессор кафедры
прикладной информатики и теории вероятностей, Российский университет дружбы народов,
sevastianov_la@rudn.university
Котюков Михаил Михайлович, аспирант кафедры прикладной информатики и теории вероятностей,
Российский университет дружбы народов, kotukov_mm@rudn.university

Note on the authors:
Kulyabov Dmitry S., candidate of physics and mathematics, associate professor of Department of Applied
Probability and Informatics, Peoples’ Friendship University of Russia; Employee of the Laboratory of
Information Technologies, Joint Institute for Nuclear Research, kulyabov_ds@rudn.university
Gevorkyan Migran N., candidate of physics and mathematics, associate professor of Department of Applied
Probability and Informatics, Peoples’ Friendship University of Russia, gevorkyan_mn@rudn.university
Demidova Anastasiya V., candidate of physics and mathematics, assistant professor of Department of Applied
Probability and Informatics, Peoples’ Friendship University of Russia, demidova_av@rudn.university
Korolkova Anna V., candidate of physics and mathematics, associate professor of Department of Applied
Probability and Informatics, Peoples’ Friendship University of Russia, korolkova_av@rudn.university
Sevastianov Leonid A., doctor of physics and mathematics, full professor of Department of Applied Probability
and Informatics, Peoples’ Friendship University of Russia, sevastianov_la@rudn.university
Kotukov Mikhail M., postgraduate student, professor of Department of Applied Probability and Informatics,
Peoples’ Friendship University of Russia, kotukov_mm@rudn.university

</pre>