ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2007, том 47, < 7, с. 1221-1228
УДК 519.676
ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЯ КОЛМОГОРОВА-ФЕЛЛЕРА1-*
© 2007 г. Н. А. Баранов, Л. И. Турчак
(119991 Москва, ул. Вавилова, 40, ВЦ РАН) e-mail: banial@yandex.ru; turchak@ccas.ru Поступила в редакцию 27.11.2006 г.
Предлагается конечно-разностный метод для решения интегродифференциального уравнения Колмогорова - Феллера. Построенная расчетная схема является безусловно устойчивой маршевой схемой, причем граничные условия определяются на основе явного решения исходного уравнения в граничных точках. Библ. 4. Фиг. 6.
Ключевые слова: интегродифференциальное уравнение Колмогорова - Феллера, метод конечных разностей, устойчивость конечно-разностной схемы.
введение
Традиционной задачей теории безопасности и надежности является определение вероятности нахождения системы в том или ином состоянии для заданного момента времени функционирования системы. При этом часто приходится рассматривать системы с непрерывным ограниченным множеством состояний 5, например 5 е [0, 1], которые изменяются под действием некоторого марковского разрывного процесса, характеризующего внешние воздействия на систему. Последние приводят к изменению ее состояния, в частности к повреждениям системы.
При предположении о том, что состояния системы упорядочены, а восстановление системы отсутствует, переходы системы из состояния 5Х возможны только в состояния 52 такие, что 52 > 5Х. При этом состояние 5 = 1 является поглощающим (см. [1]) и соответствует ситуации прекращения функционирования системы (в результате исчерпания ресурса, разрушения или других факторов).
Таким образом, в этом случае задача состоит в определении для произвольного момента времени г плотности распределения состояний системы р(5, г), которая удовлетворяет интегродиф-ференциальному уравнению Колмогорова - Феллера (см. [1], [2]). В настоящей работе рассматривается численный метод решения уравнений такого типа, основанный на конечно-разностной аппроксимации исходного уравнения.
постановка задачи
Рассмотрим уравнение вида
5 1
дрд—^ = | р (а, г) К (а, 5, г) йа - р (5, г )| К (5, а, г) йа. (1)
0 5
Функцияр(5, г) определена на области [0, 1] х [0, <»), а К(а, 5, г) - на области [0, 1] х [0, 1] х [0, <»). Ядро интегральных членов в правой части уравнения (1) удовлетворяет следующим условиям:
1) К(а, 5, г) > 0;
2) К(а, 5, г) = 0 при а > 5;
3) К(а, 5, г) - непрерывная по каждому аргументу функция;
4) К(1, 5, г) = 0.
5) I1 К (5, а, г)йа < С(г) для любого 5 < 1, где С(г) > 0 - некоторая функция времени.
Работа выполнена при финансовой поддержке РФФИ (коды проектов 06-07-89275, 07-07-00118) и программы фундаментальных исследований ОМН РАН № 3.
Требуется найти решение уравнения (1), удовлетворяющее начальному условию
p(0) = po(^), (2)
где p0(s) > 0 - непрерывная функция.
Прежде чем приступить к численному решению уравнения (1), выясним некоторые свойства решений этого уравнения.
некоторые свойства решения
Рассмотрим некоторые свойства решения уравнения (1). Для этого проинтегрируем это уравнение по s:
1 1 s 11
1 ^P(dt ds = ЦPt)K(°'s, t)d°ds - Jp(s, t)Jk(s,o, t)dods.
0 0 0 0 s
Меняя в левой части порядок интегрирования и дифференцирования, а в правой части - порядок интегрирования в первом члене, получаем
1 11 11 ddtJp(s, t)ds = Иp(o, t)K(o, s, t)dsdo - Jp(s, t)JK(s,o, t)dods.
0 0 s 0 s
Нетрудно видеть, что правая часть полученного уравнения тождественно равна нулю. Таким образом, мы показали, что
1
ifp (s,t) ds s 0,
0
т.е. решение уравнения удовлетворяет условию нормировки вида
1
Jp(s, t)ds = P = const, (3)
0
где константа P, в силу начального условия (2), определяется по формуле
1
P = Jp0 (s) ds.
0
Далее заметим, что при s = 0 уравнение (1) принимает вид
1
^pM = -p(0, t)JK(0, o, t)do.
0
Относительно функции p(0, t) это уравнение является обыкновенным дифференциальным уравнением, решение которого с учетом начального условия (2) можно записать в явном виде:
, 11
p(0, t) = p0(0) exp
-JJK(0, o, Z)dodZ
00
(4)
конечно-разностная аппроксимация уравнения (1) Для построения конечно-разностной аппроксимации уравнения (1) введем сетку узлов
= кИ, ^ = ] т,
где И, т - шаги сетки по пространству и времени соответственно, а к, - номера узлов.
ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЯ КОЛМОГОРОВА-ФЕЛЛЕРА 1223
Введем в рассмотрение сеточные функции
Рк = Р(), Ккт = К(аЬ ) .
Интегралы в правой части уравнения (1) заменим конечно-разностным приближением, используя метод трапеций численного интегрирования определенных интегралов из [3]:
| Р (а, ) К (а, 8т, t] ) й<5 = 2 рОКм + 2 £ РКт + РтК
-]
тт
к = 1
N - 1
\ К(8 т, С, t] )йа = И Ктт + 2 £ Ктк + К
мN
? V к= т +1
причем при т = 1 аппроксимационное выражение для первого интеграла имеет вид
|Р(а, tj)К(а, tj)йа = 2(Рт-1 км-1, т + РтКтт),
о
а при т = N - 1 аппроксимация для второго интеграла записывается в виде
| Р (а, ^) К (а, 81, t]) йа = И (К^ + К]ш).
о
Таким образом, конечно-разностная аппроксимация правой части уравнения при 8 = 8м, 1 < т < < N - 1, t = tj имеет вид
I Р(а, ^) К (а, 8 т, ^) йа - р (8 т, ^) | К( 8 т, а, ^) йа
8м
\ _ ( N-1
<> 2 £ К„к + К
к
к = 1
2 ^
1
mN
V к= м+1
V к = м + 1
рОк0м + 2 £ ркК
при 8 = 8т, т = 1, I = ^ - вид
|Р(а 0)К(а, 8m, ^)йа - Р(8m, ^^К(8m, а ^)йа = И>рОК0м - И>2 £ К^к + К
О
при 8 = 8м, т = N - 1, t = ^ - вид
IР (а, tj) К (а, 8 т, tj) йа - р (8 т, tj )| К( 8 т,а, tj) йа = ^ рОкОм + 2 £ ркК
О 8м V к =1
При т = N с учетом того, что уравнение (1) при 8 = 1 имеет вид
1
| р (а, t) К (а, 1, t) йа,
mN
т - 1
И
^кт
_ИЛ К
2* т^mN•
др ( 1 , t) = д t
конечно-разностная аппроксимация его правой части запишется в виде
1 / N-1
\р (а, ^) К (а, 1, ^) йа = И рО^ + 2 £ рА + РмК
NN
к = 1
8
т
О
8
т
О
О
Таким образом, конечно-разностная аппроксимация уравнения при 1 < т < N - 1 запишется в виде
/ 1 + 1 I 1 + 1 1 I 1
1 (Рт + Рт - 1 Рт + Рт-1
при т = 1 - в виде
Л . ( N-1
1 +1
р0+1 Кот1+2 £ рк+1 Кь+1 -;трт+12 £ 11+к
+-
4
р0К0т + 2 £ р]кК)
к = 1
Л . í N-1 1
И +1
^mN
к = 1
4 рт
V к = т + 1 \
2 £ 1 + К
1
mN
\ к = т +1
1 (р1и + р1щ - 1 + р1щ - 1
2 2 )
+ Ир0К1т -
- в виде
р 1 + 1 + р 1 + 1 рт р т- 1 р /я + р т - _
/ N - 1
^ _ И + 1+ 1 И р + 1 2 V К + 1 + К
I 4 р0 К 0т 4 рт 2 £ Ктк + Кт
V •
N - 1
2 £ Ктк + KmN
4рт
V к _ т+1 \
V к _ т+1
р0+1 К0+1 +2 £ рк +1
к _ 1
пт + 1 Кт +1
4Р т J^^mN
+ -
4
р0К0т + 2 £ рК
V к _ 1
После преобразований получаем при 1 < т < N - 1
- И р К1
4Р mГi^mN■
{ { N -1
2
2 £ Кт+к: +Кт+/
V к _ т + 1
р1 +1 _ р1 р т р т
{ { N -1
1 -И.Т
- 2
2 £ 1 + К
1
mN
т - 2
+ 1Т £ Кк+ 1 рк + 1 - ( 1 - 1Т К^, т ) р1+-\ + у р0К0т + 1Т £ К1трк + ( 1 + ИТ К1 - 1, т ) р1 - 1,
V к _ т + 1
т-2
р' +1К+1
2 р0 0т
к _ 1
при т = 1
/ / N -1
1 + 1ит 2
2 £ Кт+к
V к _ т + 1
к _ 1
( ( N-1
и т
И т ,
I 1 щ Кт -1, т 1- 1 + I 1 + щ Кт -1, т |р1 -1,
1-2
И т ,
2 £ Ктк+к
1
mN
V к _ т + 1
при т = N - 1
1 + иТКтл/ ^р;т 1 _ (1 - иТ^т^р'т + 1 К0т1+ Ит £ Ккт^ рк ^ ( I-ИтК'т^1, т ХрЯ-1 +
к _ 1
-2
+ у р0К0т + ИТ £ Кктрк + ( 1 + ИТК]т - 1, т)р1 - 1 .
к _ 1
ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЯ КОЛМОГОРОВА-ФЕЛЛЕРА Для узла т = N имеем
1 (Рл^+РуЛ Р1N + РJN -Л = XV 2 2
+
РО +1 кО^ + 2 £ Рк +1 к^1+ pN+1 К к = 1
N - 1
Р]КШ + 2 £ Р]КШ + РК
т+1
NN
к = 1
откуда вытекает равенство
и^К^)Р1;1 = Р1N[ 1+ + Иг ро +1 КО^+ Ит£ р[+1 к^1- р^+-11( 1-Их к^д N) +
к = 1
+ у p0к0N + и х£ р[к]ш+Р1N -1 (1 + Их к^ -1, N).
к = 1
Построенная конечно-разностная схема для исходного уравнения (1) является маршевой схемой, однако она требует задания значения сеточной функции в узле т = 0 на каждом временном слое, т.е. для нее требуется задание граничного условия на левом конце интервала. В качестве граничного условия можно использовать полученное ранее явное решение (4) для функции р(0, 0.
Заметим, что для вычисления значений сеточной функции в узле т = N на каждом временном слое вместо конечно-разностного уравнения может быть использовано условие нормировки (3).
Запишем численное значение для интеграла от решения:
|Р( 8, ^ +1 )й8
т + 1
РО
■2 £ Рк +1+ РN1
к = 1
откуда, принимая во внимание условие нормировки, находим
т+1 = 2 р -
РN -РО
N - 1 Л
т + 1
РО +1 + 2 £ Рк
к = 1
анализ устойчивости конечно-разностной схемы
Для оценки устойчивости построенной расчетной схемы воспользуемся спектральным признаком устойчивости из [4]. Представим сеточную функцию в виде
Имеем
рМ = V? ехр (¡ту).
1 + у С]т+1)^] + 1ехр (¡ту) = ^ 1-И2 С]т)^] ехр (¡ту) +
+ V
т + 1
+ ^
у ком1 + Их£ К]к+т1 ехр (¡ку) - (1 - И х Км+-11, т) ехр [ i (т - 1) у ]
к = 1 т - 2
И2хКОм + Их£Ккмтехр(¡ку) + (1 + ИхК]т-1,т)ехр[i(т - 1 )у]
к = 1
N-2
О
1226 где
Ст 2 £ Ктк + KmN.
к _ т + 1
Поделив левую и правую части полученного соотношения на |т ехр(/т ф), получим
т-2
ИТ
|[1 + ехр (-/ф)] + | 2
Ст+1 - К0т2 ехр(-/тф) - 2 £ К^1 ехр [-/(т - к)ф] + 2КЯ+Д техр(-/ф)
_ 1 + ехр (-/ф)
к _ 1 т - 2
С1 - К0техр(-/тф) - 2 £1 ехр [-/(т - к)ф] + 2К'т-1,техр(-/ф)
к _ 1
откуда найдем выражение для спектрального множителя
I _
1 + ехр(-/ф) - ^Рт 1 + 008ф - уРт - / 8Шф
1 + ехр (-/ф) + И 1 1 + 008 ф + 1 - / 8Ш ф
Ип^т + 1 2
где
т-2
Р]т _ С]т - К0т ехр (-/тф)-2 £1 ехр [-/(т - к)ф] + 2Кт-1, т ехр (-/ф).
к _ 1
Для оценки устойчивости построенной конечно-разностной схемы рассмотрим частный случай, когда ядра интегральных членов являются стационарными функциями, т.е. не зависят от времени и, следовательно,
р _ р + 1
тт
Тогда
( 2- ИТ Рт )( 1 + 008 ф) + ( Рт )2
|2 _-т1— < 1
(2 + И т )(1 + 008 ф) + (И4- ()2
что означает безусловную устойчивость построенной расчетной схемы при 1 < т < N - 1. Аналогично устойчивость расчетной схемы может быть показана при т = N - 1.
Таким образом, при 1 < т < N расчетная схема является безусловно устойчивой.
результаты численных расчетов
Рассмотрим некоторые результаты численных экспериментов. Для численного анализа возьмем ядро уравнения и начальное условие в виде
К (а, 5, г) _
X
а< 5, а< 1,
1- а'
0, а > 5 или а _ 1,
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.