Ерешко Арт. Ф.,
Вычислительный центр им. РАН,
Свентокшиская Академия в Кельцах, Польша
АНАЛИЗ ЯВНЫХ ЧИСЛЕННЫХ МЕТОДОВ
РЕШЕНИЯ СТОХАСТИЧЕСКИХ
ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
Рассматриваются основные принципы построения численных методов решения стохастических дифференциальных уравнений (СДУ). Анализируется проблема жесткости систем СДУ. Для одномерного СДУ Ито сравнивается точность аппроксимации существующих явных численных методов.
1. Введение
Анализ и синтез стохастических динамических систем часто связан с использованием численного решения СДУ. Для ряда задач таких, как фильтрация, идентификация, прогнозирование и оптимальное управление, интегрирование численного решения СДУ должно выполняться в реальном времени и, кроме того, с определенной точностью и устойчивостью. В связи с этим возникает ряд проблем. С одной стороны, очень не многие СДУ имеют аналитические решения (в основном это – линейные СДУ с аддитивным или мультипликативным шумом или нелинейные СДУ, сводимые к линейным [1]), а с другой – физические особенности реальных динамических систем [2] приводят к проявлению жесткости, что неудовлетворительно влияет на получаемое численное решение. Поэтому особо важным этапов при проектировании стохастической динамической системы является выбор схемы численного решения СДУ.
2. Принципы построения численных методов решения
стохастических дифференциальных уравнений
В настоящее время существует несколько подходов создания численных схем решения СДУ. Одной из возможностей является адаптация существующих для обыкновенных дифференциальных схем (ОДУ) схем с учетом свойств стохастических интегралов, другой – разработка специальных методов решения СДУ [3]. Большинство исследователей использует первый подход, поскольку теория численного решения ОДУ хорошо разработана и достаточно легко можно провести аналогии между ОДУ и СДУ.
Самым простым методом аппроксимации численного решения СДУ (с вычислительной точки зрения) является метод Эйлера, разработанный Маруямой в 1955г [3]. Эта схема удовлетворяет многим необходимым свойствам, предъявляемым к численным методам (она имеет порядок сходимости 
), но в тоже время обладает рядом ограничений (не всегда устойчива, ошибка аппроксимации достаточно высока и т. п.). Для устранения этих недостатков, а также повышения порядка сходимости численных схем решения СДУ были проведены и до сих пор ведутся исследования, направления которых можно представить в виде схемы (см. рис. 1).
По аналогии с разработкой схем численного решения ОДУ для повышения порядка сходимости, точности аппроксимации и устойчивости можно использовать разложение в ряд в точке аппроксимации, т. е. использования производных различных порядков, как переменной, так и коэффициентов дрейфа и диффузии. В литературе этот подход получил название метода Тейлора [3]. Однако недостатком схем Тейлора является то, что на каждом шаге аппроксимации требуется вычислять кратные стохастические интегралы, связанные с вышеуказанными производными. Для того, чтобы избежать вычислительные трудности можно использовать многократное деление шага аппроксимации (методы Рунге-Кутта [4]) или результаты аппроксимации предыдущих шагов (многошаговые методы [5 – 7]).
Как обыкновенные, так и стохастические системы дифференциальных уравнений, описывающие многие физические, биологические или экономические явления, при компьютерном моделировании с использованием обычных численных схем демонстрируют «нежелательное» поведение и могут быть отнесены к классу некорректных задач. В большинстве случаев под «нежелательным» поведением понимается очень высокая нестабильность численного решения, связанная с так называемым явлением жесткости. Существует несколько возможных объяснений этого явления.
Первая причина ассоциируется с техническими возможностями компьютера. Так для достижения желаемой точности можно применить многократное деление шага интегрирования. С одной стороны, это приводит к накоплению ошибки округления, и как следствие, возникает переполнение регистров компьютера. С другой стороны использование очень малых значений шага интегрирования требует огромных ресурсов времени и также приводит к накоплению ошибки округления. Вторая причина связана с физической стороной рассматриваемой системы. Это означает, что система описывает процессы различных скоростей или градиентов (прежде всего это характерно для некорректных задач). Такое явление обычно выступает в задачах пограничного слоя (гидродинамика), скин-эффекта (электромагнетизм), реакции химической кинетики и т. п. Наконец, жесткость может быть вызвана обеими причинами. Поэтому при разработке стабильных численных методов требуется учитывать вышеуказанные ситуации.
Анализ современной литературы показал, что создание численных методов решения жестких систем в большинстве случаев основано на идеях, представленных Хайрером и Ваннером [8]. В своей работе они постулировали, что жесткие системы не могут быть решены явными методами, и представили подходы, основанные только на использовании неявных методов. Однако следует отметить, что непосредственное применение этих методов всегда связано с крайне сложной процедурой определения параметров схемы, основанной на заранее выделенной области устойчивости только для рассматриваемой системы. Это обстоятельство делает предложенные подходы не приемлемыми для большинства вышеуказанных приложений, но позволяет выделить два важных математических свойства жесткости. Во-первых, все жесткие системы обладают очень широким спектром (или присутствием очень разных экспонент Ляпунова). Во-вторых, согласно теореме единственности и существования решения, для жестких систем характерны большие значения константы Липшица.
Итак, анализ принципов создания численных схем решения СДУ показал необходимость тщательного исследования существующих и, возможно, поиска новых методов, при решении конкретных задач.
3. Явные сильные численные схемы
Запишем СДУ в представлении Ито в общем виде
, 
, (3.1)
где 
- 
-мерный вектор состояния системы с начальным условиями 
; 
-
мерный процесс Винера; 
и 
- непрерывно дважды дифференцируемые функции дрейфа и диффузии; 
- 
мерный вектор параметров.
Получение сильного решения СДУ (3.1) является важным моментом во многих практических задачах, целью работы является сравнительный анализ существующих сильных явных численных методов решения СДУ.
Рассмотрим наиболее распространенный в финансовой литературе случай [9, 10] – случай одномерного уравнения (3.1), используя схемы: Эйлера, Мильштейна, Тейлора, Рунге-Кутта и двухшаговую [3]. В одномерном случае схема Эйлера имеет вид:
, (3.2)
где 
и 
(
).
Схема Мильштейна, обладающая порядком сильной сходимости 
, представляется как
, (3.3)
схема Тейлора порядка 
:
(3.4)
а двухшаговая схема порядка 
:
, (3.5)
где
Схема Рунге-Кутта, где порядок сходимости 
, может быть задана как
(3.6)
где
, 
, 
,
,
.
4. Сравнение схем по критерию абсолютной ошибки
Оценка качества аппроксимации той или иной численной схемы обычно связана с главной целью проводимого исследования или интерпретацией получаемого решения. Поскольку решение СДУ может быть решением в сильном или слабом смысле, то введение критерия качества аппроксимации должно учитывать этот факт. В работе рассматривается случай сильного решения СДУ, поэтому в такого критерия можно использовать критерий «абсолютной ошибки» или математическое ожидание абсолютного значения меры близости между результатом аппроксимации 
и аналитическое решением СДУ (3.1) 
на конце интервала интегрирования 
:
, (4.1)
где 
- оператор математического ожидания.
Заменим теоретическое значение критерия «абсолютной ошибки» (4.1) его статистическим аналогом, основываясь на моделировании Монте-Карло. Для этого предположим, что имеются по 
траекторий аналитического и численного решения процесса, описываемого СДУ, на конце интервала интегрирования . Обозначим их соответственно как 
и 
, тогда статистический аналог критерия (4.1) есть
(4.2)
Сравним вышеописанные схемы по критерию абсолютной ошибки. В качестве первого тестового примера исследуем линейное СДУ с постоянными однородными коэффициентами
, (4.3)
аналитическое решение которого имеет вид
.
Вторым тестовым примером является нелинейное СДУ Ито вида
с дифференцируемой функцией 
и общим решением
, где 
.
В частности, для уравнения
(4.4)
аналитическое решение есть
.
Выполним вычислительные эксперименты с использованием компьютера типа IBM PC с процессором AMD Athlon (TM) XP 3000+ для тестовых уравнений (4.3) и (4.4), используя численные схемы (3.2) – (3.6) и исследуя зависимость между длиной шага интегрирования 
, количеством траекторий и точностью аппроксимации (4.2). Результаты вычислений приведены в таблицах 1 – 3, проанализируем их, используя усредняющий критерий (4.2).
Для первого и второго тестовых уравнений (см. табл.1 и табл.2) при уменьшении длины шага интегрирования и увеличении порядка сходимости численной схемы возрастает точность аппроксимации для всех исследуемых численных схем.
Однако этого нельзя утверждать в третьем случае, который представлял жесткое СДУ [11] (см. табл.3). Удалось рассчитать значение абсолютной ошибки для всех комбинаций длины шага интегрирования и количества траекторий только для схемы Эйлера и двухшаговой схемы.
Таблица 1. Точность аппроксимации численного решения уравнения (4.3) (
, 
, 
, 
)
Схема | Длина шага интегрирования, | ||||||
|
|
|
|
|
|
| |
| |||||||
Эйлера | 0.2131 | 0.4025 | 0.2903 | 0.0239 | 0.2540 | 0.0228 | 0.0051 |
Мильштейна | 0.0368 | 0.0976 | 0.7012 | 0.0658 | 0.0201 | 0.0201 | 0.0055 |
Тейлора | 0.1513 | 0.3566 | 1.1820 | 0.0282 | 0.0749 | 0.0125 | 0.0007 |
двухшаговая | 0.2131 | 0.0486 | 0.0846 | 0.0117 | 0.0056 | 0.0001 | 0.0015 |
Рунге-Кутта | 0.0906 | 0.0296 | 0.0039 | 0.0015 | 0.0043 | 0.0010 | 0.0001 |
| |||||||
Эйлера | 0.1917 | 0.4300 | 0.2606 | 0.4059 | 0.3476 | 0.1705 | 0.0899 |
Мильштейна | 0.2039 | 0.3218 | 0.1000 | 0.1050 | 0.4573 | 0.0273 | 0.0123 |
Тейлора | 0.1446 | 0.4783 | 0.1026 | 0.1302 | 1.1038 | 0.0351 | 0.0123 |
двухшаговая | 0.1917 | 0.2896 | 0.0527 | 0.0499 | 0.0987 | 0.0079 | 0.0036 |
Рунге-Кутта | 0.0768 | 0.0709 | 0.0326 | 0.0181 | 0.0114 | 0.0026 | 0.0006 |
| |||||||
Эйлера | 0.2999 | 0.3999 | 0.6878 | 0.3729 | 0.2091 | 0.2180 | 0.1382 |
Мильштейна | 0.2728 | 0.2800 | 0.3892 | 0.1170 | 0.0608 | 0.0368 | 0.0220 |
Тейлора | 0.1994 | 0.3079 | 0.5748 | 0.3086 | 0.0737 | 0.0445 | 0.0389 |
двухшаговая | 0.2999 | 0.1941 | 0.1792 | 0.0496 | 0.0197 | 0.0132 | 0.0052 |
Рунге-Кутта | 0.1119 | 0.0757 | 0.0862 | 0.0234 | 0.0064 | 0.0036 | 0.0010 |
Для схем Мильштейна, Тейлора и Рунге-Кутта при 
, 
, 
значения абсолютной ошибки было гораздо выше, чем для схемы Эйлера или двухшаговой схемы, а при уменьшении длины шага интегрирования (
, 
, 
, 
) происходило переполнение регистров, что приводило к невозможности проведения дальнейших вычислений.
Таким образом, можно отметить, что в отличии от ОДУ, при численном интегрировании решения жестких СДУ следует использовать «простые» явные методы решения, т. е. избегать методов, использующих многократного деления шага аппроксимации или производных функций дрейфа и диффузии. В случае потребности численного решения СДУ в таких задачах, как фильтрация или идентификация параметров СДУ с использованием процедуры Монте-Карло [12 – 13], предпочтительной длиной шага является 
.
Таблица 2. Точность аппроксимации численного решения уравнения (4.4) (, 
)
Схема | Длина шага интегрирования, | ||||||
|
|
|
|
|
|
| |
| |||||||
Эйлера | 0.1618 | 0.2674 | 0.0448 | 0.1437 | 0.0964 | 0.0017 | 0.0830 |
Мильштейна | 0.0001 | 0.0483 | 0.0194 | 0.0073 | 0.0029 | 0.0018 | 0.0001 |
Тейлора | 0.0102 | 0.0074 | 0.1385 | 0.1513 | 0.0977 | 0.2649 | 0.1172 |
двухшаговая | 0.1618 | 0.1181 | 0.2250 | 0.1733 | 0.0537 | 0.2545 | 0.1243 |
Рунге-Кутта | 0.0571 | 0.0502 | 0.0001 | 0.0390 | 0.0638 | 0.0040 | 0.0073 |
| |||||||
Эйлера | 0.1847 | 0.1754 | 0.1164 | 0.1317 | 0.0512 | 0.0316 | 0.0253 |
Мильштейна | 0.0298 | 0.0268 | 0.0088 | 0.0069 | 0.0045 | 0.0017 | 0.0008 |
Тейлора | 0.0013 | 0.0496 | 0.0564 | 0.0847 | 0.0871 | 0.0647 | 0.1415 |
двухшаговая | 0.1847 | 0.1162 | 0.0655 | 0.0912 | 0.0863 | 0.0633 | 0.1432 |
Рунге-Кутта | 0.0710 | 0.0721 | 0.0401 | 0.0669 | 0.0194 | 0.0100 | 0.0168 |
| |||||||
Эйлера | 0.1883 | 0.1992 | 0.1380 | 0.0818 | 0.0727 | 0.0505 | 0.0374 |
Мильштейна | 0.0228 | 0.0248 | 0.0116 | 0.0060 | 0.0032 | 0.0017 | 0.0010 |
Тейлора | 0.0001 | 0.0726 | 0.1040 | 0.1022 | 0.0859 | 0.0967 | 0.1198 |
двухшаговая | 0.1883 | 0.1003 | 0.1159 | 0.1092 | 0.0898 | 0.0981 | 0.1207 |
Рунге-Кутта | 0.1106 | 0.1443 | 0.1051 | 0.0477 | 0.0329 | 0.0280 | 0.0217 |
Таблица 3. Точность аппроксимации численного решения уравнения
(4.3)(, 
, , )
схема | длина шага интегрирования, | ||||||
|
|
|
|
|
|
| |
| |||||||
Эйлера | 4.1449 | 3.4965 | 8142.3 | 0.4759 | 0.0181 | 0.0001 | 0.0001 |
Мильштейна | 8.3574 | 184.31 | 37.620 | 0.0020 | 0.0001 | 0.0001 | 0.0001 |
Тейлора | 48.501 | 24278 | ∞ | ∞ | ∞ | ∞ | ∞ |
двухшаговая | 4.1449 | 51.147 | 33723 | 36.716 | 0.2408 | 0.0718 | 0.0829 |
Рунге-Кутта | 46.144 | 1156.1 | ∞ | 21.454 | 0.0001 | 0.0001 | 0.0001 |
| |||||||
Эйлера | 6.3458 | 43.712 | 1065.8 | 65.253 | 133.38 | 0.0022 | 0.0017 |
Мильштейна | 23.326 | 178.74 | 1843.7 | 2695.4 | ∞ | ∞ | ∞ |
Тейлора | 73.265 | ∞ | ∞ | ∞ | ∞ | ∞ | ∞ |
двухшаговая | 6.3458 | 380.59 | 17858 | 25.88 | 0.7215 | 0.3548 | 0.1109 |
Рунге-Кутта | 64.016 | 1835.8 | 28289 | 19.238 | ∞ | ∞ | ∞ |
| |||||||
Эйлера | 4.0971 | 39.017 | 373.01 | 10934 | 3379.2 | 1.9114 | 0.0015 |
Мильштейна | 14.133 | 389.1 | 16903 | 2829.6 | 0.0546 | ∞ | ∞ |
Тейлора | 40.655 | ∞ | ∞ | ∞ | ∞ | ∞ | 0.0001 |
двухшаговая | 4.0971 | 636.98 | 26002 | 2848.3 | 7.1647 | 2.1052 | 0.1743 |
Рунге-Кутта | 38.225 | 2636 | 44269 | 831.95 | ∞ | ∞ | ∞ |
Рис. 1. Методы и направления преобразования схемы Эйлера
Литература
1. Oksendal B. Stochastic differential equations. Berlin: Springer, 2000.
2. , Методы решения некорректных задач. М.: Наука, 1986.
3. Kloden P. E., Platen E. Numerical solution of Stochastic Differential Equations. Berlin: Springer, 1999.
4. Burrage K., Tian T. Stiffly accurate Runge-Kutta methods for stiff stochastic differential equations // Computer Physics Communications. 2001. V. 142. P. 186 – 190.
5. Burrage K., Burrage P., Mitsui T. Numerical solutions of stochastic differential equations – implementation and stability issues // Journal of computational and applied mathematics. 2000. V. 125. P. 171 – 182.
6. Kuznetsov D. F. The three-step strong numerical methods of the orders of accuracy 1.0 and 1.5 for Ito Stochastic differential equations // Journal of Automation and Information Sciences. 2002. V. 34. № 12. P. 22 – 35.
7. Gaines J. G., Lyons T. J. Variable step size control in the numerical solution of stochastic differential equations // SIAM Journal of Applied Mathematics. 1997. V. 57. № 5. P. 1455 – 1484.
8. Hairer E., Wanner G. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. Berlin: Springer-Verlag, 1996.
9. Lamberton D., Lapeyre B. Introduction to stochastic calculus applied to finance. London: Chapman and Hall. 2000.
10. Н. Основы стохастической финансовой математики. М.: ФАЗИС, 1998.
11. Milstein G. N., Platen E., Schurz H. Balanced implicit methods for stiff stochastic systems // SIAM Journal of Numerical Analysis. 1998. V. 35. P. 1010 – 1019.
12. Filatova D., Grzywaczewski M., McDonald D. Estimating parameters of stochastic differential equations using a criterion function based on the Kolmogorov-Smirnov statistics // ACTA ET COMMENTATIONES UNIVERSITATIS TARTUENSIS DE MATHEMATICA. 2004. V. 8. – P. 93 – 99.
13. Nielsen J. N., Madsen H. Applying the EKF to stochastic differential equations with level effects // Automatica. 2001. V. 37. P. 107 – 112.
14. Nielsen J. N., Madsen H., Young P. C. Parametric estimation in stochastic differential equations: an overview // Annual Reviews in Control. 2000. V. 24. P. 83 – 94.
