Алгоритмы численного решения стохастических дифференциальных систем с переключаемой диффузией Текст научной статьи по специальности «Математика»

УДК 519.6+004.94 ББК 22.193

АЛГОРИТМЫ ЧИСЛЕННОГО РЕШЕНИЯ СТОХАСТИЧЕСКИХ ДИФФЕРЕНЦИАЛЬНЫХ СИСТЕМ С ПЕРЕКЛЮЧАЕМОЙ ДИФФУЗИЕЙ

1 2 Черных Н. В. , Пакшин П. В.

(Арзамасский политехнический институт (филиал) Нижегородского государственного технического университета им Р.Е. Алексеева)

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

Ключевые слова: стохастические системы, диффузионные процессы, марковские переключения, схемы Тейлора, сходимость, устойчивость.

1. Введение

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

1 Надежда Валентиновна Черных, аспирантка ([email protected]).

2 Павел Владимирович Пакшин, доктор физико-математических наук, профессор (pakshinpv@gmail. com).

няется также их тесной связью с уравнениями математической физики.

Один из известных и перспективных подходов к численному интегрированию стохастических дифференциальных уравнений Ито основан на стохастических аналогах формулы Тейлора. Теоретически это разложение позволяет строить методы как угодно высокого порядка точности (при соответствующих предположениях о коэффициентах системы уравнений).

Важнейшей отличительной особенностью стохастических аналогов формулы Тейлора для решения стохастических дифференциальных уравнений Ито является присутствие в них, так называемых, повторных стохастических интегралов в форме Ито или Стратоновича, которые являются функционалами сложной структуры относительно компонент векторного винеровского процесса [3]. Аналогом обычных степеней винеровского процесса служат в теории Ито многочлены Эрмита, а обычную экспоненту от стохастического интеграла дублирует экспоненциальный супермартингал. Проблема совместного численного моделирования (исходя из среднеквадратического критерия сходимости) совокупностей повторных стохастических интегралов Ито является сложной как с теоретической, так и с вычислительной точки зрения.

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

Исследование скачкообразных систем началось с работ I. YA. Kats, NN. Krasovskii, Е.А. Lidskii (1960-1961), а затем Т. Kazangey и D.D. Sworder (1971) представили скачкообразную систему, где макроэкономическая модель народного хозяйства использовалась, чтобы изучить эффект федеральной политики по

стабилизации жилищного сектора. Известна работа Г.Н. Мильштейна (1972) «Mean square stability of linear systems driven by Markov chain». A.S. Willsky и B.C. Levy рассматривали скачкообразные системы для моделирования систем электроэнергии и управления солнечным тепловым центральным приемником. M. Mariton в своей работе «Jump Linear Systems in Automatic Control» поясняет, что скачкообразные системы появились как удобная математическая структура для моделирования разнообразных проблем в различных областях, таких как задачи слежения, контроля допустимых ошибок и производственных процессов. [11]

Как правило, рассматриваются переключения диффузии по закону марковской цепи с конечным числом состояний - марковские переключения.

При компьютерном моделировании таких систем возникает необходимость численного решения СДУ с марковскими переключениями.

Известны работы по изучению управляемости, устойчивости, стабилизации таких систем. Вопросы численного решения СДУ с марковскими переключениями изучены мало.

В работах [9, 10] рассмотрена известная численная схема Эйлера (для стохастических систем) в применении к решению стохастического дифференциального уравнения с марковскими переключениями, приведены доказательства некоторых результатов частного характера.

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

2. Предварительные сведения

Пусть (Q, F, P) - вероятностное пространство, Ft, t0 < t <t0 + T

- неубывающее семейство о - подалгебр F, (rnr(t), Ft), r = 1,...,d -независимые винеровские процессы.

Рассмотрим стохастическое дифференциальное уравнение

Ито

й

(1) йХ = a(t, X )Л + ^&г (t, X )йаг ^),

г=1

где X, а, аг - векторы размерности п. [4]

Предполагается, что функции а(^ х) и аг(1, х) определены и непрерывны при t e[t0, t0 + Т], х еЧЯ” и удовлетворяют условию Липшица при всех t е [t0, t0 + Т], х е^ ” , у е^ ” :

(2) |а^,х) - a(t,у)| + |ст(t,х) -ст(t,у) < К|х - у| . [4]

Будем использовать далее следующие обозначения:

|х| означает евклидову норму вектора х, ху - скалярное произведение векторов х и у; К - положительная константа.

Х4х(0 или просто Х(0 - решение уравнения (1), удовлетворяющее начальным данным Х4х(0) = х. Разобьем промежуток ^0, ^ + Т] точками деления 4 на N равных частей, так что tk+\ - tk = h, к = 0,

1, ..., N - 1, ^ + Т = ^, h = T/N. Приближение к Х(Д) будем обозначать X (^), где Х0 = X (^). Далее пусть X - -измеримая

случайная величина и Е\Х1 < да; Х( х ^) означает решение уравнения (1) для 4 < t < t0 + Т, удовлетворяющее начальным данным при t = tk.

Пусть (X(t), F(t)), t0 < t < t0 + Т - некоторое решение уравнения (1), с начальным условием, удовлетворяющим неравенству

Е|Х(^)|2 < да.

Определим одношаговую аппроксимацию Х(х + h),

t0 < t < t + h < t0 + Т, которая формируется в зависимости от х, t, h, и {ю(9) - ),...,(9) - <юй^): t < 9 < t + к\:

(3) Х1х ^ + h) = х + / (t, х, h;юi (9) - а>1 (t), i = 1,...й, t < 9< t + h),

где функция / определяется конкретным методом нахождения решения.

На основании одношаговой аппроксимации рекуррентно построим приближение (Хк,Ft ), к = 0, ..., N, ^+\ - tk = ^+\,

tN = to + Т:

(4) X0 = X0 = X (^ ), Хк+! = Х^,Хк (tk+l) =

= Хк + / (tk, Хк, К+1; ю (9) - ю (tk), i = и й, tk < 9 < tk+1).

Для простоты считаем, что tk+\ - tk = h = Т / N. [4]

Теорема 1. Пусть одношаговая аппроксимация Х(х ^ + И)

имеет порядок точности рх для математического ожидания отклонения и порядок точности р2 для среднеквадратичного отклонения, т. е., при любых t0 < t < t0 + Т - h, х е ^” выполняются неравенства

(5) ЕХ,,х (t + h) - X,,х (t + И)) < К(1 + |х|2)1/2 И*

Г — ^ “|1/2 2

(6) е|х,х(, + И)-Х,х(, + h)| < К(1 + |х|)1/2ИР2

и пусть

(7) р2 > 1/2, р1 > р2 + 1.

Тогда при любых N и k = 0, 1, ..., N выполняется неравенство

|2"

(8)

EXto, Xo(tk ) Xt„, Xo(tk )|

й К{1 + E|X0|2 )12 h

т.е. порядок точности метода, построенного с использованием одношаговой аппроксимации Xtx (t + h), равен р = р2 - -2. Постоянные ^ не зависят от X0 и N. [4]

3. Постановка задачи

Пусть (Q, F, P) - вероятностное пространство, Ft,

t0 < t < t0 + T - неубывающее семейство о - подалгебр F, rnr(), r = 1, ..., d - независимые винеровские процессы. Пусть M = {1, ..., m} - конечное множество.

Рассмотрим стохастическое дифференциальное уравнение с марковскими переключениями в форме

(9) dX(t) = a(0(t),X(t))dt + ]T ar (0(t),X(t))dar (t),

r=1

где P(t) - однородный марковский процесс со счетным множест-

вом состояний M, 0(0) = u0, X (0) = x0.

(10) р(р(, + И) = 1|^(0 = и,х(^),Р(,$),s < qul(,)И + о(и), и ФI,

где х(0 е^ п, а(-,-): ^ п х М ^ ^ п и ст(-,-): ^ п х М ^ ^ пхп.

Переходная функция такого процесса определяется набором функций Р(,,и,I) = ри1 (t), образующих стохастическую матрицу Р(0 переходных вероятностей (ри1 ^) > 0, Е ри1 ^) = 1), где числа

I

ри1 имеют смысл вероятностей перехода из и в I за интервал И при условии, что процесс в(0 на этом интервале покинул состояние и.

Эволюция процесса Р(,) целиком описывается с помощью чисел qul, ql, образующих матрицу интенсивности переходов

Q(t) = (ди1 ^)) е^тхт, Q(t) = Нт РР(-)—Р(0), где для каждого t

^0 t

т

qul(t) > 0 при и ФI, quu =-qu, Е qul (0 = 0для кажцого и е М ; [1]

I=1

Предполагается, что функции а(в(0, х(,)) и ^(А^), х(0) определены и непрерывны при t е ^0, t0 + Т], х еЧЯп и удовлетворяют условию Липшица при всех t е [^, t0 + Т], х е ^ п , у е ^ п , и е М:

(11) |а(и, х) - а(и,у) + |ст(и,х) -и(и,у) < К|х - у|, а также условию:

(12) |а(и, х)| + |ст(и,х)| < К(1 + |х|),

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

При выполнении всех описанных выше условий уравнение (9) имеет единственное, непрерывное решение Хи, х на интервале t > 0 для каждого начального значения. [8]

Поставим задачу численного решения уравнения (9), т. е. нахождения аппроксимации его решения на рассматриваемом временном интервале.

При такой постановке проблемы, когда матрица переходов Q не зависит от х, а зависит только от t, ее можно свести к известной задаче решения стохастического дифференциального уравнения вида (1) с помощью схем, описанных в [2, 4, 7] на совокуп-

ности случайных подынтервалов [0,^Х^,^ + t2),...,где ^ - случайные моменты переключения марковской цепи, которые можно определить заранее. Так как следующий момент преключения 4 и состояние в(4) есть случайные величины, распределение которых зависит лишь от (в^к), 4), можно предварительно сгенерировать совокупность (в^к), 4) правых частей уравнения (1) и решать уравнение (1) известными методами. При этом точность аппроксимации, которая зависит от выбора шага на отдельных подынтервалах, будет достижима и на всем рассматриваемом временном интервале.

Но можно предложить для решения уравнения (9) другой подход.

Для стохастического дифференциального уравнения вида (1) Е. Платен предложил очень простой вывод (использующий лишь формулу Ито) разложения решения X, х(, + И) по степеням И и интегралам, зависящим от приращений сог (3) - сог ^), где t < 3 < t + И, г = 1, ..., d. Это разложение в детерминированном случае представляет собой формулу Тейлора для Х{, х(, + И) по степеням И в окрестности точки ^, х). [2, 4, 7]

Будем использовать разложение Платена, подробно рассмотренное в [4], для решения уравнения (9), заменив функцию // х) функцией /в(0, х(0) с переключаемой компонентой.

Пусть Хи, х^) = Х^) - решение уравнения (9) для всех и е М ; /(в, х), где в = в(0, X = Х(0 - достаточно гладкая функция (скалярная или векторная). Согласно формуле Ито имеем для ^ < t < 3 < ^ + Т

(13) /(Р(3),Х(3)) = /08, х) +

С 3 3

+Цл / (8(3,), х & ))сфг 3)+} ь/ (8(3!), х & )>3,

г=1 t t

где операторы Лг, г = 1, ..., С и Ь определены как:

а

(14) Лг = |^ — |,

их

( а\ 1ч п п а2

(15) Ь = | а,— \ +1 ----,

к ’ I ах J 2 г=1 к и г г ах1 бх]

Применим формулу (13) к функциям Лг/ и Ь/, а затем полученные выражения для Лг/(8(3),X(3)) и Ь/(8(3),Х(3)) подставим в (13). Имеем

(16) /8), х^) ) = / + £ Л г/\Саг (3) + Ь/\С3 +

г=1 t t

С 5 ( d 3

+11 11Л5Л/(8(3:),X(3 ))С®, (3 )С©г (3) +

г=1 t V 5=1 t

+1 | (| ЬЛ г/(8(31), X (3 ))С3 ^ (3)+

г=1 t V t

+1 | (| Л Ь/ (8(31), X (31 )У©г (31 ))С3 +

г=1 t V t

+} Г| ь2/(8(31), X (31))С31 )С3,

t V t

где, например, Лг/ вычисляется в точке (в, х).

Поступая так дальше можно получить разложения для т + И), Х(, + И)), где роль степеней выполняют случайные величины вида (которые не зависят от Ft)

,+и 3 3 3]-2

(17) 1^ (И) = |(3)|С^_13)|... |С®11(3;-1),

t t II

где /1, ..., 1] принимают значения из множества чисел 0, 1, ..., С, t < 3< 31 < 32 < ... < 3]--1 < t + И и под Сю0(3г) понимается С3г.

Интеграл (17) соответствует /-степени параметра И в детерминированном разложении.

Очевидно, что Е1: : = 0, если хотя бы одно из 1к ф 0,

'1 ^^

к = 1, ..., ], и Е1^ ь = 0(И]), если все 1к = 0, к = 1, ..., /

Рассмотрим примеры вычисления нескольких интегралов Ито. Для этого будем использовать формулы (18) из теории интегралов Ито ([4, 5, 6]):

г „2 (^) _ + г г

(18) |ю(З^ю(З) =-; JЗdю(З) = ію(і) _|ю(З^З;

0 2 0 0

Ісо2(З^ю(З) = _ |„(З^З ;

0 3 0

г + h г + h

(19) І0(И) = |dЗ = Л ; ІДИ) = |da(З) = „(г + И) _„(г) = А„(И);

гг

г + И З г + И З г+и и ^

100(И) = | da0(З)|da0(З1) = |dЗJdЗ1 = JЗdЗ = —;

г г г г г 2

г+и З г+и Дгл2(И) _ И

Іи(И) = | da(З)|dю1(З1) = | со(З^а(З) =-------------—-;

г г г 2

г + И З г + И З г+И

І10(И) = | da0(З)|da(З1) = | dЗJda(З1) = |ю(З^З;

г г г г г

г + И З г+И З г+И

101(И) = | d„(З)|da0(З1) = | d„(З)|dЗ = JЗdю(З) =

г г г г г

г+И

= ИАю _ |ю(З^З = ИА„(И) _І10;

г

г+И З З

І1,1,1 (И) = | dю(З){dЮl(Зl) |dю2(З2) =

г г г

г+И З г+И „ 2(0) _ з

= | dю(З)jю1(З1)dю1(З1) = | ——---------da(З) =

г г г 2

Аю3(И) 1

=----------- ИА„(И);

6 2

В формулах (19) все рассмотренные повторные интегралы

г + И

Ито выражаются через случайные величины А„(И), |ю(З^З .

г

Вернемся к разложению (16) и выпишем формулу, которая получается путем ряда непосредственных подстановок:

t+И

t+И

(20) /8 + И),X^ + И)) = / + £Лг/1 с®г(3) + ь/1С3 +

г=1 t

с с

t+И

3

+ ЦЛ,Л г/ | Саг (3){ (31) +

г=11=1 t t

с с с ,+и 3 3

+ 1ЦЛ 5 Л/Л г/ | стг (3)| с®, (31) | с®, (32) +

г=1 / =1 5=1

с t+И 3

t

t

г

с t + И 3

+£ль/1с3|с®г(31)+^ьл/ | с®г(3)|с31 +

г=1 t t

1+И 3

г=1

t

г

+ Ь2/ |с3{ с31 + р, где / - /8), х(1)),

t t с с с с (+ И (3 3 32

(21) р = £ИЦ 1 (1 (/Л ]Л, Л,Л,/ х

г=1 1 =1 5=1 ] =1 t V t t t

х (8 (33), х (33 ))с®, (33 )У®; (33 ))с®5 (32) )с®, (31 )^® г (3) +

с с t+И(»(3 ,

+И1 1 1 ЬЛ, л г/(8(32), X (32))с32 ^3) К (3)

г=1 , =1 t

+

Vt V

d й (+И (3(3

+И1 1 1л,ЬЛг/(8(32), Х(32))с® (32)>3 )с®г (3) +

г=1 , =1 t

d с t+ И(3(3

+И1 1 1л,ЛгЬ/(8(32), Х(32))с® 3)^ (31))с3 +

г=1 , =1 t

с с с t+И (3(3 32

+ Ш1 1 1 (1 ЬЛ Л Л г/(8(33), X (33)^33^ (32)):

г=1 , =1 5=1 t V t V t t

х с®, (31))4®г (3) +

с t+и(3(3

+Ц 1 I Ь2Лг/(8(32), Х(32))с32 >3 )с®г (3) +

г=1 t

с t+ И(3(31

+11 1 1 ьл гЬ/ (8(32), Х(32))с32 )с®г (31))с3 +

г=1 t

Й г + ИҐЗҐЗ

+и І \лд2/(р(з2).Х(0))й„г(0)И )о +

г=1 г

!+И ( 3(31

+ 1 1 1 ь3/(8(32), х (32) ^ >3 )с3.

t V t V t

Рассмотрим леммы, доказанные в [4], применительно к разложению (20).

Лемма 1. Справедливо соотношение

»2) (е > )Ц

= О

( 1 2-і, \ И4= 2

V

Другими словами, при подсчете порядка малости интеграла (17) нужно руководствоваться правилом: йЗ вносит в порядок малости единицу, а (З), г = 1, ..., й- одну вторую.

Лемма 2. Пусть

< К (і + |х|2 У2.

(23)

Л г1 ...Л Ч/ (Р(г ). Х(г ))

Тогда величина

(24) I (/. И) =

г+И з З З_2 , .

: І (З)| йф З) І... |Л, ...Л 11/(р(З]_1). X (З^фЗ)

удовлетворяет неравенству

(25) Е

,(/. И)

< К (1 + Е|Х (г )|2 И=

2(2_гі )

т. е., в частности, ее порядок малости тот же, что и у величины 1,г , (И). Далее, если хотя бы один из индексов 1к,

к = 1, ...,] отличен от нуля, то

(26) Е1 (/,И) = 0, £/2 * 0.

к=1

Из этой леммы вытекает, что каждое слагаемое, составляющее р, имеет не более чем второй порядок малости. Более того, математическое ожидание всех членов второго порядка малости

2

и порядка малости 2,5 из р равно нулю согласно (26). Поэтому \Ер\ = 0(И ) Разумеется, это верно, если для всех подынтегральных функций из р выполняется, например, (23).

Леммы 1, 2 справедливы для разложения (20) в силу того, что в формуле (20) интегралы Ито остались в том же виде, в котором они используются в разложении Платена (см. [4, 7]).

На основании этих лемм непосредственно подсчитывается порядок каждого слагаемого в разложениях вида (20).

Легко увидеть [4], что можно получить как разложения, содержащие все члены включительно до какого-то полуцелого порядка малости, так и до какого-то целого порядка.

Возьмем теперь в (20), (21) в качестве /в, х) вектор х. В этом случае Лг/ = аг, Ь/= а. Поэтому d И

(27) Хи, х ^ + И) = х + Хстг 1 (3) + аИ +

Г=1

г

d d 1+И

+ ЕЕЛг^г 1(юг (3) - а. (t))dЮг (3) +

Г=1 I =1 t

d 1+И d

+ Х Ьег 1 (3-1 )<я?юг (3) + £л га 1(юг (3) - аг ^ ))d3 +

I г=1 t

t+и ( 3 И ^

+ЁЁЁЛ,Л г°г 1 !(®, (31) - (t ))аЧ- (3))яЧ(3)+Ьау+р.

г=1 d d d

В формуле (27) все коэффициенты аг, а, Л. аг, Ьог, Лг а, ЛS Л. аг, Ь а вычисляются в точке (в, х), а остаток р равен

d d d d t+И (3 31 32

(28) р = ЦЦ 1 1 (1 ({Л;ЛЛ.СТг (83), Х(3з)) х

г=1 1=1 £=1 ] =1 t V t t t

х da : (33 )^а3 (32 )^юг (31))<3юг (3) +

| +

1

d d t + И (3 3

+ И 1 1 (1ЬЛ г (83), X(32))d32)d®г (3))^ (3)-

г=1 г =1 t V t t

d d t + И (3 31

+ И 1 1 ( {Лг-Ь^г (8(32), Х^)^. (3)^3 ^ (3) +

г=1 I =1 ! V t t

| +

й d t + Н С$ $

+ И І I(|ЛгЛга(8(32),Х(32))й®г(32)й®г(3,))й3-

Г=1І=1 t V t t

й й й t + Н С $ $ $2

+ Ш І І(І(ILЛsЛгar(8($),Х(33))й33)йа6.($2)) х

Г=1 І =1 Я=1 t V t t t

х й©г ($1))й©г ($) +

й t+ Н С $ $1

X 1 1 (1ЬЧ (8(32), Х(32)М32М3 )d®г (3) +

г=1 t V t t

d t+ И ( 3 31

+ 1 1 1 (1ЬЛ га(8(32), Хф)^^ (31)) +

г=1 t V t t

d t+ И ( 3 31

+ 1 1 1 (\ЛгЬа(Р(32), Хф)^ ф)^ )3 +

г=1 ! V t t

1+И (3 3

+ 1 1 (11}а(Р(32), X (32)^32^3 )d3.

t V ^ t

В связи с формулами (27), (28) рассмотрим следующую одношаговую аппроксимацию:

(29) Хих((+ Н) = х + Естг(©г(ґ + Н)-©г(0) + аН +

Г=1

й й t+Н й t+

+ ££Лі*г I (©г ($) - ©г (t))й©г ($) + ЁLCTr I($-Ой©г ($)

Г=1 і =1 t

й t+Н

+ ЁЛга |(©г ($) - ©г (0)й$

Г=1 t

й й й

+

t + Н

г=1 t

+

+

t + Н[ $

Н2

Ё££Л£Л 1СТг 1 I }(®£ (31) - ©£ ^(31) ^г (3) +^"2 ,

г=1 I =1 £=1 t V t / 2

где а = а(в^), х(1)), стг = стг(в(0, х(1)).

Пользуясь леммами 1 и 2, можно установить (при выполнении условия (23) для соответствующих функций), что для одношаговой аппроксимации (29)

|Ер| = 0(И3), Е| р2 = 0(И4 ),

т. е. р1 = 3, р2 = 2, и этот метод имеет порядок сходимости равный 3 / 2 (на основании Теоремы 1).

Общий результат построения одношаговых аппроксимаций обоснован и сформулирован в [4].

Все многократные интегралы Ито вида (17) имеют множителем Аш(Н) (см. формулы (19)), т. е. величины, способные принимать сколь угодно большие значения. Тем не менее, схемы, основанные на разложении Платена, обеспечивают хорошую аппроксимацию решений стохастических дифференциальных уравнений, что подтверждается компьютерными экспериментами, например в [7].

Обобщим теорему Платена (приведенную в [4]) на рассматриваемый случай решения стохастического дифференциального уравнения с марковскими переключениями (9). Все выводы теоремы Платена остаются в силе для разложения (27), так как выполняются леммы 1, 2.

Теорема 2 (обобщение теоремы Платена). Пусть

Хих (ґ + Н) содержит все слагаемые вида Л г ...Л ^ • I^ ^,

где / = х, до целого порядка т включительно. Пусть все функции

] 2 - і к

Лг ...Лі f (Р(ґ),х(ґ)), где / = х, ^--- < т +1, удовлетворяют

7 1 к=1 2

неравенству (23). Тогда среднеквадратичный порядок точности метода, основанный на такой аппроксимации, равен т.

Пусть Хих (ґ + Н) содержит все слагаемые вида

Лг ...Л/ ь, где / = х, до полуцелого порядка т + -2 вклю-

Ґ+Н $ $т-1 Нт+1

чительно и слагаемое і!"а І й$... | й$т = і!" а----------------.

І І І т (т +1)!

Пусть все функции Лг ...Лі^(8(0,х(ґ)), где / = х,

7 2 - і и

^------- < т + 2, удовлетворяют неравенству (23). Тогда сред-

к=1 2

неквадратичный порядок точности метода, основанный на такой аппроксимации, равен т + 2.

Для нахождения решения уравнения (9) предложим следующую численную схему, построенную на одношаговой аппроксимации (29), т. е. среднеквадратичного порядка точности 3 / 2.

где ( = Р(ґ), х = х(г), Ь2к = | |dasds2 , Аак - есть N(0, И) гаус-

совски распределенных приращений винеровского процесса ш на подынтервалах гк < г < гк+1, Х( = х0.

Рассмотрим теперь модель переключения, зависящего от фазового состояния.

Вернемся к уравнению (9), где 8(0 - однородный марковский процесс со счетным множеством состояний М, в(0) = и0, Х(0) = х0,

(31) Рр(г + h) = /|Р(0 = и,x(s),Р(s),s < г) = qul(x)h + о(к), и Фї, где х(0 п, а(у): ^ п х М ^ ^ п и ст(-,-): ^ п х М ^ ^ пхп.

Такая состояние - зависимая модель находит значительно большее применение в управлении и оптимизации. Если матрица Q зависит от состояния х(0, то случай становится значительно

гк+1 52

гк гк

(32) Пусть б(-): ^п ^ ^тхт - ограниченная и непрерывная функ-

ция. Q(x )=(Яиі(х)) тхт для каждого х, qul (х) — 0 при и Фї,

т

X Чиї (х) = 0 для каждого и е М .

1=1

более сложным. Одна из главных трудностей - та, что из-за непрерывной состояние - обусловленной связи, ув(0 и х(0 становятся зависимы. в(0 - является цепью Маркова только для фиксированного х, то есть фактически не является марковской. [9, 10]

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

Суть задачи в следующем - рассмотреть пару процессов в(0 и х(0 совместно; т. е. процесс с двумя компонентами рассматривается как марковский. [9, 10]

Пусть переключаемая диффузия (в(0, х(0) задана уравнением (9). Пусть первое начальное значение (в(0), х(0)) = (и0, х0), также задано другое начальное значение (в(0), у(0)) = (и0, у0) для у ф х. Для состояние - зависимого случая ри°,х°^) Ф ри°’Уо^) бесконечно часто, даже при том, что первоначально ри0,х°(0) = ри0,У0(0) = Р , где верхний индекс показывает начальное значение зависимости.

Будем считать, что в(?) - стохастически непрерывный процесс, для которого + И) ^Р(0 при И ^ 0. В соответствии с

требованием сепарабельности будем считать, что в(?) не может менять свое состояние за «нулевое время» более одного раза. В этом случае траектории процесса будут кусочно - постоянными, т. е. временной интервал разбивается на подынтервалы [0, t\), [^, t1 + t2),..., на которых Р(?) постоянно, tk - случайные моменты переключения марковской цепи. [1]

Далее будем рассматривать последовательность Р(^) как дискретно - временной стохастический процесс, аппроксимирующий в(?) в соответствующем значении.

Используя функцию Q(x) можно построить матрицу переходных вероятностей Р = е^х)И . Используя разложение функции в ряд Тейлора I + ЛQ( х) + 0(Л2)и отбросив 0(Д2) ввиду ограниченности и непрерывности Q(x), будем моделировать матрицу

переходных вероятностей как Р = I + ДQ(х), где I - единичная матрица.

Будем рассматривать пару процессов в(0 и х(0 совместно как марковскую цепь следующим образом.

Значения Х^ генерируются рекурсивно согласно (30), используя предыдущие значение Хг , и одновременно генерируются значения Р( ^, также используя значение Хг (х = Хг в матрице переходных вероятностей Р = I + ДQ( х)).

При всех вышеописанных условиях Леммы 1, 2 и обобщенная теорема Платена выполняются для случая марковского процесса, зависящего от фазового состояния, т. е. схема (30) будет обеспечивать нахождение приближенного решения уравнения (9) с матрицей переходов Q(x).

Подобный подход используется в [9, 10], где представлены леммы и теоремы частного характера (рассматриваемая в этих работах схема Эйлера является частным случаем схемы (30)).

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

4. Пример

Рассмотрим процесс Ито X = {X, t > 0}, удовлетворяющий линейному стохастическому уравнению

(33) dXt = / (Pt, Xt) Xtdt + g (Pt, Xt) Xtdnt

на временном интервале [0, T] при T = 1, X0 = 1.

Пусть pt - однородная марковская цепь, M = {1, 2} - число состояний марковской цепи, Q - матрица интенсивности переходов, P - матрица переходных вероятностей.

/(fit, Xt) и g(fit, Xt) - непрерывные, ограниченные, дважды дифференцируемые функции, удовлетворяющие условиям (11),

(12).

(34) Xt=x о exp f/(Pt, xt) - 2 g 2(pt, xt ^+g(Pt, xt )nt

есть решение уравнения (33) для t є [0,T] и данного винеровского процесса n = {nt, t > 0}.

Перепишем схему (30) в виде

(35) Yn+l = Yn + a(Pn, Yn )Л n + a(Pn, Yn )I^ +

+ a(Pn, Yn )a'(Pn, Yn )10,„ + a'(Pn ,Yn )a(Pn, Yn )ї0,о) +

+ f a(Pn ,Yn )a'(Pn ,Yn) +1 ^2(Pn, Yn )a"(Pn ,Yn) \ +

2

1 — <

2

+ (a(Pn ,Yn )a\pn ,Yn) +1 a2(Pn ,Yn )a"(Pn, Yn ))I (0,1) +

+ а(Рп, Уп )(<т(Р„ ,Уп )а"(Рп ,Уп) + (а'(Рп, Уп ))2)/ (111),

где а(в, х^ = fl.Pt, Xt) Xt, о($, Xt) = #($, Xt) Xt для рассматриваемого уравнения (33).

Для численного моделирования интегралов Ито введем в рассмотрение независимые N(0, А) - распределенные случайные величины

_ 1 ,-- _ з (т»+1 1 ^

(36) & =Лп2Л®п и ^2 =^Г2Дп2 | (а(3) _ а(т))^_ -ДпЛ®п

у *п

С помощью этих величин получим:

1 V _1 (1 1 Л

(37) Дап = Д2п^!, | (ю(3) _@(т'У)$3 = Дп2 “£ + ^= ^2 .

тп у2 -V12 у

Тогда формулы (19) можно переписать в виде:

1 _____ ( 1 Л

1 (1,0) — "2 л/Ап Ап ■ + ^3 ^2 ; 1 (0,1) — А®п Ап - 1 (1,0);

1 (1,1,1)— -2 А®п2 - А п ^А®п; где А®п— 4^;

1 (1

2 ^ 3

Зададим начальные значения 70 = Х0, р0 = и0 и будем рекурсивно генерировать 100 значений Уп с равным значением шага А согласно (35), где Дп - есть длина временного интервала дискретизации ґ0 = т0 < т1 < ... < тп < ... < = Т на временном интервале

Т];

Для сравнения будем использовать (34), чтобы определить соответствующие значения точного решения, используя ту же примерную траекторию винеровского процесса юг на подынтервалах тп < ґ < тп + 1, а именно

(( 1 Л п Л

(38) Хп+1 — X 0 ехр^ /(Рп, Хп) - - g \рп, Хп )^п + g (Рп, Хп )£А®г

Рассмотрим результаты численного решения уравнения (33), выбирая различные варианты задания матрицы переходов Р, шага дискретизации А, значений функций/Д, х), g(вt, х)

Случай 1.

Пусть /(рь х) принимает два значения - (а1, а2}, соответствующие первому и второму состоянию марковской цепи. Положим а1 = 1,5, значение а2 будет изменяться в ходе эксперимента. g(вt, х) принимает два значения - (Аь А2}. Положим А1 = 0,1, значение А2 будет изменяться в ходе эксперимента.

Будем полагать, что начальному моменту времени соответствуют значения а1 и А1.

1. Р =

( 0,2 0,8Л

ч 0,4 0,6 ,

а2 = 10а1; А2 = 50Я1;

File Tools Edit

|0|6ЕР[ tj\

1.5 strong Taylor scheme

Рис. 1. A = 0,01 (t e [0;0,1])

Scilab Graphic {2) [ <—> | I°1

File lools Edit

I | O | GEP| k\

--------------- xt

------- yt

Рис. 2. A = 0,0001 (t e [0;0,01])

В данном случае значения коэффициентов меняются не очень значительно, поэтому устойчивая аппроксимация решения уравнения наблюдается с размером шага дискретизации Д = 0,01.

(0,3 0,7'

2. Р =

0,6 0,4

; а2 = 100а1; А2 = 200А1;

В данном примере оба коэффициента увеличиваются значительно (в 100 и 200 раз), поэтому выбор меньшего шага дискретизации заметно повышает точность вычислений.

Случай 2.

fifit, x) принимает два значения - {ai, «2}, соответствующие первому и второму состоянию марковской цепи. g(fit, Xt) принимает два значения - {Аь Я2}. Будем полагать, что начальному моменту времени соответствуют значения ai и Х\.

( 0,2 0,8 ^

3. P = ; ai = cos x + sin x;

0,4 0,6

a2 = 1 + cos x; X2 = 0,2x;

& Э - .Е. Ь

Рис. 6. А = 0,01 (t є [0;0,1])

э|эНЧИ

■ ..............--------11"-----М4ММЧ.,и,.------

Рис. 7. Марковская цепь

4. Р =

а1 = cos х + sin х; а2 = 1 + cos х;

А1 = 0,5х; Х2 = 2х;

О 10 20 30 40 50 60 70 80 90 100

----------

------ У*

__________________________________________ /а

Рис. 8. А = 0,01 (t е [0;0,1])

Я.1 ЗсПаЬ СгарЬк (0) 1 ■=■ I Ш 1иД»Г

£|1е ТооЬ Е<)й

;э|э|о|сев| 1»|

--------------У1

Рис. 9. А = 0,001 (I е [0;0,01])

Рис. 10. Марковская цепь

Варианты 3, 4 численного моделирования отличаются матрицей переходов Р и значением Л,2.

Случай 3.

5. б =

^ - 50cos3 х 50cos3 х ^ 10sin х - 10sin х

; Р = I + б А;

а = 2 + sin х; Л1 = 0,5;

а2 = 1 + sin х cos х;

Ь = о, 1;

О 10 20 30 40 50 60 70 80 90 100

— хГ

- уг

Рис. 11. А = 0,01 (t е [0;0,1])

0 10 20 30 40 50 60 70 80 90 100

-------- ЪеПа

А

^ БсЛаЬ ОгарЫс (2) ^ ЕНе 1оо1б Ес1И:

-------------- х!

------------- У*

Рис. 13. А = 0,001 (t е [0;0,01])

У ЕсПаЬ СгарЫс {2) [ = 1 е

БЕо| 1-г |

2.4-

2.2-

2.0

1.8-

1.6

1.4-

1.2-

1.0

0.8-

0.6

10 20 30 40 50 60 Ьега 70 80 90 100

6. б =

а! = 20 + 0082 кх = 0,2;

- 5 008 2 х

10 008 2 х - 10 008 2 х J

; Р = I + б А;

х; а2 = 8т х + 5008 2х;

^2 = 3;

Ейе 1оо|5 Ейй

^>|^3>|о|бЕр| Ь I

* I

Рис. 15. А = 0,001 (t е [0;0,01])

Я' 5а1аЬ йгарЫс (3) ~ 1-=. | И |м£Ы

П1е 1оо15 №

£>1э1оЫ ь!

^ Scilab Graphic (1)

File lools Edit

^1'

1. 5strong Taylof scheme

Рис. 17. А = 0,0001 (t є [0;0,001])

^ Scilab Graphic [1)

File Tools Edit

Me\o\ 6ED |

2.4-

2.2-

2.0-

1.8-

1.6-

1.4-

1.2-

1.0-

0.8-

0.6

10 20 30 40 50 60 70 betta 0 90 100

/a

Эксперимент показывает, что в случае состояние - зависимой модели переключения предложенная численная схема также позволяет найти приближенное решение уравнения (33).

а1 = 20 + cos2 х; а2 = sin х + 5 cos 2х;

Л-i = 0,2; I2 = 0,3;

Этот вариант данных отличается от предыдущего значением Я,2.

Scilab Graphic (7)

ImU

File lools Edit

|э|о|«ер| k\

Рис. 19. A = 0,01 (t e [0;0,1])

JU

10 20 30 40 50 60 70

^|э|о|«| Е»|

*£ ЗсПоЬ СгзрЬнс (10) 11=1 | В |»£5ЙГ

ЕИе 1оо1'5 Е(1^

Рис. 23. А = 0,0001 (t е [0;0,001])

БсИаЬ вга рЫс{10) Ь-Иэиа*

РНе Тоок Е(й

$>\Р г. 2.4- 2.2- 2.0 в|

1.6- 1.4- 1.2- 1.0 0.8 Об

10 20 30 ЬеЯа 3 50 60 70 80 90 100

7. Р =

Г - 5^2 х

10cos х - 10cos х

аі = 3 + cos х; Л = 0,5;

а2 = 20 - sin х; ^2 = 0,1;

Рис. 25. А = 0,01 (t є[0;0,1])

Scilab Gra phic (0) 1 = 1 в мы

File Tools Edit

bed| k\

2.4- 2.2- 2.0-

1.6-

1.4-

1.2-

0.8- 0.6-

10 20 30 40 50 60 70 80 90 100

File lools Edit

д>|э|с|ч

Рис. 29. А = 0,0001 (t є [0;0,001])

*£ Scilab Graphic {2)

File Tools Edit

eb|

2.4-

2.2-

2.0-

1.6- 1.4-

1.2-

1.0

0.8-

0.6-

10 20 30 40 50 60 70 0 90 100

betta

Л

5. Заключение

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

Возникающие вопросы устойчивости схем можно решить, учитывая, что во все слагаемые в той или иной степени входит множителем h. При достаточно малом h устойчивость схем будет обеспечена. Эксперименты показали, что для многих случаев подходящим является шаг дискретизации h = 0,01, но для случаев сильного колебания значений функций а(в, X), о(в, X) в уравнении (9) устойчивость будет обеспечивать выбор более малого шага h = 0,001, h = 0,0001.

В случае состояние - зависимой модели переключений остается открытым вопрос - что считать точным решением и как определять меру отклонения от получаемых приближенных значений.

Также необходимо отметить, что для практических применений важнее оказывается не поиск решения, как реализации единственной траектории, хотя бы и с заданной точностью, а скорее нахождение значений функционалов от случайных траекторий, таких как математическое ожидание интегральных функционалов и/или вероятностей пребывания в заданном множестве. Нахождение этих величин с заданной точностью методом Монте-Карло требует предварительной оценки числа итераций, которое также должно быть связано с точностью аппроксимации.

Литература

1. БОРОВКОВ А.А. Теория вероятностей: Учеб. пособие для вузов. - 2-е изд, перераб. и доп. - М.: Наука. Гл. ред. физ.-мат. лит., 1986. - 432 с.

2. КУЗНЕЦОВ Д.Ф. Стохастические дифференциальные уравнения: теория и практика численного решения. - Спб.: Издательство Политехнического университета, 2007. - 800 с.

3. КУЗНЕЦОВ Д.Ф. Стохастические дифференциальные уравнения: теория и практика численного решения // Дифференциальные уравнения и процессы управления. -2008. - № 1.

4. МИЛЬШТЕИН Г.Н. Численное интегрирование стохастических дифференциальных уравнений. - Свердловск: Изд-во Уральского университета, 1988. - 224 с.

5. МИЛЛЕР Б.М., ПАНКОВ А.Р. Теория случайных процессов в примерах и задачах. - М.: ФИЗМАТЛИТ, 2002. - 320 с.

6. ОКСЕНДАЛЬ Б. Стохастические дифференциальные уравнения. Введение в теорию и приложения: Пер. с англ. - М.: Мир, ООО «Издательство АСТ», 2003. - 408 с.

7. KLOEDEN P. E., PLATEN E., SCHURZ H. Numerical Solution of SDE Through Computer Experiments. - Berlin: Springer -Verlag, 1994.

8. MAO X. Stability of stochastic differential equations with Markovian switching // Stoch. Proce. Appl. - 1999. - Vol. 79. - 45-67.

9. YIN G., MAO X., YUAN C., AND CAO D. Approximation methods for hybrid diffusion systems with state-dependent switching processes: numerical algorithms and existence and uniqueness of solutions// SIAM Journal on Mathematical Analysis - 2010. - Vol. 41, №6. - p. 2335-2352.

10. YIN G., ZHU C. Stochastic modeling and applied probability. Hybrid switching diffusions. Properties and applications. // Springer Science + Business Media, LLC 2010.

11. YUAN C., LYGEROS J. Stochastic markovian switching hybrid processes // Project IST-2001-38314, COLUMBUS. Design of Embedded Controllers for Safety Critical Systems, June 17, 2004.

ALGORITHMS FOR NUMERICAL SOLUTION OF STOCHASTIC DIFFERENTAL SYSTEMS WITH SWITCHING DIFFUSION

Nadezda Chernykh, post-graduate student, Arzamas Polytechnic Institute of R.E. Alekseev, Nizhny Novgorod State Technical University, 19, Kalinina Street, Arzamas, 607227, Russia, ([email protected]).

Pavel Pakshin, Dr.Sci., professor, Arzamas Polytechnic Institute of R.E. Alekseev, Nizhny Novgorod State Technical University, 19, Kalinina Street, Arzamas, 607227, Russia, ([email protected])

Abstract: Mathematical models are considered of hybrid systems in the form of stochastic differential equations with Markovian switching of the diffusion component. An extension of Taylor schemes for numerical approximation of their solutions is proposed. Modeling results in SCILAB are presented to demonstrate efficiency of the obtained algorithms.

Keywords: stochastic systems, diffusion process, Markovian switching, Taylor schemes, convergence, stability.

Статья представлена к публикации членом редакционной коллегии А. В. Добровидовым