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

УДК 519.622+519.868

АЛГОРИТМ ПРОВЕРКИ РАЗРЕШИМОСТИ ОДНОЙ СИНГУЛЯРНОЙ СИСТЕМЫ

КАБАЛЯНЦ П.С.____________________________

При математическом моделировании задач экономики, теории электрических цепей и теории управления возникают полулинейные дифференциальные уравнения с вырожденной матрицей при производной. В работе описывается алгоритм, проверяющий выполнение условий существования решения начальной задачи для случая, когда характеристический пучок ХА + B линейной части уравнения и сопряженный с ним пучок pA * +B * имеют нетривиальные аннуля-торы KerA П KerB и KerA* DKerB*.

1. Постановка задачи и обзор

Рассмотрим систему полулинейных дифференциально-алгебраических уравнений (DAEs), которая в векторной форме имеет вид

Au'(t) + Bu(t) = f(t,u), 0 < t <т , (1)

с начальным условием

u(0) = uo (2)

и вырожденной, вообще говоря, матрицей A размерности m х n . Здесь вектор-функция u(t) действует из [0, т] в Cn , а вектор-функция f (t, u) - из [0, т] х Cn в Cm .

Уравнение (1) возникает при математическом моделировании задач экономики, теории электрических цепей, теории управления: S.L.Campbell [4], W.Leontief [11], G.G.Yerrghese [9], T.Kailath [10], P.G.Luenberger [12], L.Pandolfi [13] и др. Известно, что если матрица A необратима, то задача (1),(2) может не иметь решения даже в случае линейной правой части уравнения (1) — функции f(t). Условия существования решения линейного уравнения (1) получены К.Вейерштрассом, Л.Кронеке-ром [2]. Стандартные математические пакеты (Maple 6, Ma^mati^a 4 и др.) применимы к задаче (1),(2) лишь в случае обратимой матрицы A . В 70-е годы XX ст. в работах Сибирской школы (Ю.Е.Боярин-цев, В.Ф.Чистяков и др.) разрабатывалось математическое обеспечение решения краевой задачи для системы линейных DAEs [1]. Позже K.D.Clark и L.R.Petzold [6] применили вычислительные методы стрельбы и конечных разностей для нахождения приближенных решений краевых задач этого класса. Начальная задача для некоторых классов уравнения (1) рассматривалась в работах S. L. Campbell [5], U.Ascher [3], A.Kvaern [8]. В данной статье описываются алгоритм и программное обеспечение проверки условий существования решения для одного класса сингулярных задач (1),(2).

Матричный пучок ХА + B называется регулярным, если матрицы A, B являются квадратными и най-

дется число X 0 є C, для которого существует обратная матрица (X0A + B)_1. В противном случае пучок ХА + B называется сингулярным.

Каноническая форма сингулярного матричного пучка, полученная Л.Кронекером [2], имеет блочно-диагональный вид. На диагонали могут располагаться блоки пяти типов: два типа регулярных блоков, три типа — сингулярных. Мы рассмотрим сингулярность одного из трех типов —нулевого.

Будем говорить, что матричный пучок ХА + B и отвечающее ему уравнение (1) имеют сингулярность только нулевого типа, если существует пара невырожденных матриц P, Q , преобразующих пучок к блочно-диагональной форме, в которой ненулевой блок является регулярным:

р0.А + B)Q = [°° 1Дг°+ Br} , de«XAr + Br) * 0(3)

Для геометрической трактовки перейдем к линейным операторам A, B : X ^ Y , индуцированным одноименными матрицами в координатных базисах пространств X = Cn, Y = Cm . Из матричного представления (3) следует, что пространства X,Y раскладываются в ортогональные суммы подпространств

X = Xs © Xr, Y = Ys © Yr , (4)

где Xs = KerA n KerB, Ys = KerA * ПKerB * — анну-ляторы пучков ХА + B , цА* + B* соответственно.

Пучок операторов kAs + Bs = ХА + b|xs , индуцированный в подпространствах Xs,Ys, является нулевым. Пара подпространств (Xr,Yr) инвариантна относительно пучка ХА+B :

A(Xr) с Yr, B(Xr) с Yr; пучок XAr + Br : Xr ^ Yr, индуцированный в подпространствах Xr , Yr, является регулярным. В дальнейшем предполагается дополнительно, что det Ar Ф 0 (для матрицы Ar) или KerAr = {0} для оператора Ar : Xr ^ Yr.

Пусть размерности сингулярных подпространств равны g = dimXs , h = dim Ys . Выберем ортонормированный базис x1,...,xn в пространстве х и ортонормированный базис y1,...,ym в пространстве Y так, чтобы первые части базисов относились к сингулярным подпространствам:

span{xi }g=1 = Xs , span{y j j = Ys (5)

(span — линейная оболочка).

2. Теоремы существования решения

Рассмотрим случай линейного уравнения (1). Предположим, характеристический пучок ХА + B системы (1) имеет сингулярность только нулевого типа, и пусть det Ar + 0 . Тогда для существования решения задачи (1),(2) необходимо и достаточно, чтобы линейная функция f (t) в правой части уравнения была ортогональна подпространству Ys при всех t:

РИ, 2001, № 4

81

(f(t),yi) = 0, i = Vt є [0,x] .

Для того чтобы решение было единственным, необходимо и достаточно, чтобы подпространство Xs было тривиально [2].

Перейдем к нелинейному уравнению (1).

Теоремаї. Пусть число уравнений в системе (1) не больше числа переменных (h < g). Тогда для существования решения u(t) на некотором нетривиальном интервале [ 0, т0 ] достаточно, чтобы выполнялись следующие условия:

I) функция f(t,u) непрерывна по совокупности переменных и непрерывно дифференцируема по ив некоторой области {t є [0, т], ||u - uo|| < r} ;

II) согласование начального вектора и0 и сингулярной компоненты функции f (t, и) :

1) (f(0,U0),yj) = 0, yj єbasYs (5),j = 1,...,h ;

* Sf(0,u0)

2) rang(L 0 M) = h, где

да

L = (yi,...,yh),M = (xi,...,xg),

Sf(0,u0) _(9fj(0,u0)'їn

9ui

yi=1j=i

Если переменных столько же, сколько и уравнений (g = h), то решение будет единственным. В противном случае решение начальной задачи (1),(2) будет неединственным.

т f Л(Г

Замечание. Матрицы L, —, M имеют размерности

от

h х m, m х n, n x g соответственно. Условие согласования 2) означает, что

(

Ші

rang

(f^y,,)

СЇЇ

от„

(f0yy,„

su

g

= h

y u=u0

где дифференцирование выполняется относительно проекций ui = (u,xi), i = 1,...,g вектора u на элементы базиса {xi }g=. (5).

Доказательство. Используем разложения по ортонормированным базисам {xi }n=. , {y j j : n m u

u(t) = E uixi , f(t,u) = E fj(t,u)yj , (6)

i=1 u j=1

где ~i = (u,xi), i = 1,...,n; uj(t,u) = (f (t, u),y j), j = 1,...m.

Тогда

n n n m

Bu(t) = EuiBxi = ЁuiBxi = E Eui(Bxi,yj)yj .

i=1

i=g+1

i=g+1 j=h+1

m

nm

Аналогично Au'(t) = E Eui(Axi,yj)yj .

i=g+1 j=h+1

Следовательно, система (1) распадается на две подсистемы:

1) fj(t,u) = 0, j = 1,...,h ; (7)

Е[(Axi,yj)u- + (Bxi,yj)~і ] = ~j (t,u),

2) i=g+1 (8) j = h + 1,...,m.

Из условий утверждения следует, что существует ненулевой минор порядка h . Перенумеруем векторы {xi }g=1 так, чтобы

9~1 (0, u) 9~1 (0, u)

6111 Cuh

det

Ф 0 .

dfh(0,u)

Cui

f (0,u) сїїь

2 а=аo

Подставим в подсистему (7) разложение u(t) по базису {xi}n=1 :

fj (t, u1,. ., un) = 0, j = 1, . ,h ,

fj(t,~1,...,un) = uj(t,E~ixi), j = 1, . ,h

i=1

Применяя теорему о неявной функции к подсистеме (7), находим функции

~j =Тj(t,uh+1,...,un) J = 1,...,h .

Они определены и имеют непрерывные частные

производные

ЁЁ1

і = h + 1,...,n,j = 1,...,h в неко-

торой области

Q = {t є [0, x1],||u - u0|| < r1}, r1 < r, x1 <x . Подсистема (8) в векторной форме имеет вид

Ar~r + Br~r = f(t,u1,...,Un) , (9)

u u u tr

где ur — (ug+1 ,..., un ) ,

f(t,u1,...,un) = (fh+1(t,u),...,fm(t,u))tr .

Умножим уравнение (9) на матрицу A^1 и подставим функции uj =у(t,uh+1,...,un),j = 1,...,h :

ur + A“1Brur = A“19(t,uh+1,...,Un). (10)

Здесь функция

9(t,Uh+1,...,un) = f(t, V1(t,Uh+1,...,un),...

..., W h(t,uh+1,...,un),uh+1,...,un)

имеет в области Q непрерывные производные

9ф

СЇЇі

f h f Эф j 9ui j=1j сїїі ’

i = g +1,

n

Применяем теорему Пикара [7] к системе (10) с начальными условиями щ(0) = (u0,xi), i = g + 1,...,n.

РИ, 2001, № 4

82

Получаем непрерывно дифференцируемые решения Uj(t,Uh+i,...,Ug), i = g + в некоторой обла-

сти Qi = {t Є [0,Tq],||u - Uq|| < Г2}, Г2 < ri, Т2 <Ti .

Они определяют непрерывно дифференцируемую функцию

U(t) = Vl(t,Uh+i,...,Un)Xi +... + yh(t,Uh+i,...,Un)Xh +

+ Uh+1xh+1 +... + Ugxg +

+ Ug+i(t,Uh+i,...,Ug)xg+i +... + Un(t,Uh+i,...,Ug)xn .

(11)

h+i,..., ug

Коэффициенты базисного разложения u^+i,..., Ug выбираются непрерывными на интервале [0,тq] с начальными условиями ~ц(0) = (uQ,xi) ,i = h + i,...,g. При таком выборе коэффициентов и функция (11) есть решение задачи (1),(2)

Теорема 2. Пусть число уравнений в системе (1) больше числа переменных (h > g). Тогда для существования единственного решения на некотором нетривиальном интервале [0, тq ] достаточно, чтобы выполнялись следующие условия [10] :

I) функция f (t, и) непрерывна по совокупности переменных и непрерывно дифференцируема по и в некоторой области {t є [0, т] , ||и - иО|| < r} ;

II) согласование начального вектора uq и сингулярной компоненты функции f(t,u): ортонормированный базис {yi }h=i подпространства Ys можно выбрать так, чтобы

1) (yi,f(t,u)) = 0, i = g + 1,g + 2 ,...,h (12)

при всех t из интервала [ 0, т ] и всех и из некоторой окрестности точки uQ;

2) (f(0,uo),yi) = 0, i = i,...,g; ;

* 6f(0,uo)

3) det(L ^M) * 0 , где

би

L = (yi,...,yg),M = (xi,...,xg)

X s n m

af(o,uQ) _f afj(o,uo)

Ji=ij=i

би

би;

Замечание. Матрицы L, —, M имеют размерности

би

h х m, m x n, n x g соответственно. Условие согласования 2) означает, что

Д|Ууі )

OU1

det

(^,Уі )'

OUg

(fQU>,y.)

би„

^ o,

/ U=

UQ

где дифференцирование выполняется относительно проекций и; = (u,x;), i = i,...,g вектора и на элементы базиса {x;}g=i (5).

Доказательство. Подставляем разложения (6) для функцийu(t) ,f(t,u) в уравнение (1). Учитывая (12), получаем две подсистемы :

1) ~j(t,u) = О, j = i,...,g ; (13)

2) Е [(Ax; ,yj )ui + (Bx;,yj )~і ] = fj (t,u),

2) i=g+i

(14)

j = h + i,...,m.

Подставим в подсистему (13) разложение u(t) по базису {x;}n=i:

fj(t,ui,...,un) = o, j = i,...,g ,

fj(t,ui,...,un) = UJ(t,EliiME j = i,...,g .

i=i

Применяя теорему о неявной функции к подсистеме (13), находим функции

Uj = Фj(t,Ug+i,...,Un),j = i,...,g .

Они определены и имеют непрерывные частные

5у j

производные --------, 1 = g + i,...,n,j = i,...,g в неко-

6Ui

торой области

Q = {t Є [О, Ті] , ||u - Uq || < ri}, ri < r, xi <T . Подсистема (14) в векторной форме имеет вид

Ar~r + Brur = f(t,ui,...,un), (15)

Ur = (Ug+1,...,Un)tr ,

ГдЄ f(t,ui,...,un) = (Ug+i(t,u),...,Um(t,u))tr .

Умножим уравнение (15) на матрицу A и подставим функции uj = y(t,Ug+i,...,iin), j = 1,...,g :

ur + A("1BrUr = A("19(t,ug+1,...,un) . (16)

Здесь функция

ф(й Ug+1,..., ~n ) — f (t, Фі (t, Ug+i,..., Un ),...

..., Ф g(t,Ug+1,...,Un),Ug+1,...,Un)

имеет в области Q непрерывные производные

5ф

6~i

6f h f 6y j

6Ui j=i бф j 6~i ’

i = g +1

n

Применяем теорему Пикара [7] к системе (16) с начальными условиями 1ц(О) = (uq,x;) , i = g + 1,...,n . Получаем непрерывно дифференцируемые решения ui (t), i = g +1,..., n в некоторой области

Пі = {t є [О, tq],||u - uq|| < r2}, r2 < ri, Т 2 <Ti .

Они определяют непрерывно дифференцируемую функцию

U(t) = yi(t,Ug+i,...,Un)xi +...

u u u u (17)

... +Vg(t,Ug+i,...,Un)xg + Ug+i(t)xg+i +... + Un(t)xn'v 7

РИ, 2001, № 4

83

Функция (17) есть решение задачи (1),(2).

Приведем алгоритм проверки выполнения условий существования решения задачи (1),(2) одновременно для случаев g > h и g < h .

3. Описание алгоритма проверки разрешимости

1) Решаем системы линейных уравнений

A

х - 0 и

A'

Б'

у = 0

(18)

Векторы фундаментальных систем решений для уравнений (18) образуют базисы в подпространствах аннуляторов пучков XA + Б и pA * + Б* -Xs = KerA П KerB и Ys = KerA* П KerB* , соответственно. Если обе системы имеют лишь тривиальное решения x(t) = 0, y(t) = 0 , то делаем вывод: уравнение не имеет сингулярности нулевого типа. Иначе переходим к следующему шагу алгоритма.

2) Используем процедуру ортогонализации для нахождения ортонормированных базисов xj,...,Xg в подпространстве KerA f| KerB и yb...,yh в подпространстве KerA * flKerB *.

3) Решаем системы линейных уравнений

< X1 ^

xg

V g

х = 0,

f Уі ^

I yh.

у = 0

(19)

Векторы фундаментальных систем решений для уравнений (19) образуют базисы в ортогональных дополнениях

Xr = (KerA П KerB)1, Yr = (KerA * fKerB*)1,

соответственно.

4) Используем процедуру ортогонализации для нахождения ортонормированных базисов Xg+i,..., хп в подпространстве (KerA П KerB)1 и yh+1,...,ym в подпространстве (KerA* fKerB*)1.

5) Опредеёяем Ar = A|(KerAnKerB)1 :

' (Axg+1,yh+1)

(AxnAh+1) 1

Ar =

v (Axg+1,ym)

(Axn,ym)

6) Вычисляем detAr . Если detAr = 0 или n - g Ф m - h, то делаем вывод, что уравнение (1) имеет сингулярность не только нулевого типа. Иначе переходим к следующему шагу алгоритма.

7) Проверяем условия гладкости I теоремы 1 для функций, заданных аналитически, с помощью стандартных функций пакета символьных вычислений Maple 6.

8) Если n > m, то проверяем условия согласования II 1) и 2) теоремы 1. Иначе проверяем условия согласования II 1), 2), 3) теоремы 2, соответственно.

9) Проверка условия II 3) теоремы 2.

Для нахождения ортонормированного базиса {y д!1^ подпространства Ys, удовлетворяющего условию I, решается функциональное уравнение

аДДри) +... + а hfh(t,u) = 0

относительно чисел {аД^ с C , где

fj(t,u) = (f (t, u),y j), j = 1,...h .

Если это уравнение имеет h - g линейно-независимых решений {ak}k=1,j = g + 1,...,h , то

yi =a1y1 +... + ahyh A = g + 1,...,h .

10) Если достаточные условия выполнены и n = m, то существует единственное решение. Если n ф m, то решение неединственно.

Программная реализация алгоритма проверки выполнена в виде расширения библиотеки пакета Maple 6.

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

В работах С.Кемпбелла [4] описываются численные методы решения начальной задачи для неявных уравнений вида

F(t,u,u') = 0 . (20)

С. Кемпбелл обобщает на случай неявного уравнения некоторые классические численные схемы (Рунге-Кута и др.). Наше уравнение (1) также является неявным, т.е. не разрешенным относительно вектора производных u'(t). Однако (1) имеет более специальный вид сравнительно с общим уравнением (20), производная u'(t) присутствует только в линейной части уравнения, сингулярно вырожденной. Относительно уравнения (1) можно надеяться получить более точные теоретические результаты и алгоритмы, учитывающие специфику задачи. В настоящей статье предложен один из таких алгоритмов, решающий проблему допустимости начальных данных задачи Коши для системы дифференциально -алгебраических уравнений, имеющих векторную форму (1). При этом учитывается лишь один из трех возможных типов сингулярностей Л.Кронекера, выявление которого не требует получения полной канонической формы характеристического пучка матриц XA + B . Алгоритмическое описание допустимых начальных данных и соответствующее вычислительное обеспечение для остальных двух типов сингулярностей связано с задачей приведения сингулярного пучка к канонической форме Л.Кронекера. Полное вычислительное обеспечение этой задачи, по-видимо -му, отсутствует, и преодоление связанных с этим трудностей будет изложено в следующей публикации. Заметим, что С.Кемпбелл предполагает, что начальные данные принадлежат классу допусти-

84

РИ, 2001, № 4

мых ‘consistent initial conditions’ для общего уравнения (20), не решая задачу описания этого класса.

Литература: 1. Бояринцев Ю.Е., Чистяков В.Ф. Численные методы решения сингулярных систем. АН СССР, Сиб. Отд., Иркут. В.Ц. — Новосиб. 1989. 2. Гантмахер Ф.Р. Теория матриц. М.: Наука, 1967. 576с. 3.Ascher U. On summetric schemes and differential-algebraic equations//SIAM J.Sci and Statist.Comput. 1989. 10,№5. Р.937-949. 4. Campbell S.L. Singular systems of differential equations/Research notes in mathematics; 40: 1980, Pitman Advanced Publishing program. San Francisco. London. Melbourne. 176p. 5. Campbell S.L.,R.Hollenbeck .Automatic differentiation and implicit differential equations // SIAM Computational Differentiation: Techniques, Applications, and Tools, Philadelphia, 1996. Р.215-2270. 6. Clark K.D, Petzold L.R. Numerical solution of boundary value problems in differential-algebraic systems//SIAM J.Sci and Statist.Comput. 1989. 10,№5. Р.915-936. 7. Dieudonn J. Foundations of Modern Analysis // Academic Press: New York-London, 1960. 8. Kvaern A. Runge-Kutta methods applied to fully imlisit differential-algebraic equations of index 1//Math.Comput. 1990. 54,№190.Р.583-625. 9. Verghese G. Infinity frequence behaviour in generalised dynamical systems, Ph.D.Thesis, Stanford Univ. U.S.A.(1978). 10. Verghese G. andKailath T., ‘Impulsive behavior in dynamical

systems, structure and significance’, in Proc. 4th Int. Symp. Math. Theory, Networks Syst., Delfi,The Netherlands, July 1979. 11. Leontief W. Essays in Economics, M.E. Sharpe, Inc., N.Y., 1977. 12. Luenberger D.G. Nonlinear descriptor systems, J. Economic Dynamic and Control (1979). Р. 219-242. 13. L.Pandolfi. On the regulator problem for linear degenerate control systems, J. Optimization Theory and Applications 33 (1981). Р. 241-254. 14. Каба-лянц П. С., Руткас А.Г. О квазилинейном дифференциальном уравнении с сингулярностью типа Кронекера // Вісник Харківського університету. Серія “Математика, прикладна математика і механіка”. 1999. №444. С.111-118.

Поступила в редколлегию 21.09.2001

Рецензент: д-р физ.-мат. наук, проф. Руткас А.Е.

Кабалянц Петр Степанович, старший преподаватель кафедры математического моделирования и обеспечения ЭВМ Харьковского национального университета им.В.Н.Каразина. Научные интересы: вырожденные дифференциальные уравнения, математическое моделирование и вычислительные методы. Увлечения: музыка, пение и сочинительство. Адрес: Украина, 61004, Харьков, ул. Окт. рев., 18, кв. 2-а, тел. 12-6388.

УДК 621.391.175

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ДИНАМИКИ

ВОССТАНАВЛИВАЕМЫХ СИСТЕМ

ПОДПЫКИН Н.С._______________________

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

1. Введение

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

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

РИ, 2001, № 4

Состояние системы свяжем с вероятностью отказа. Обычно состояние не наблюдаемо. Но всегда найдутся контролируемые параметры, которые содержат информацию о состоянии. Мы не будем рассматривать здесь задачу об идентификации состояния. Ее решение можно найти, например, в [5]. Заметим лишь, что в процессе идентификации состояния невозможно учесть все информативные параметры (признаки), определяющие его. Поэтому состояние системы в каждый момент времени следует считать случайной величиной, а процесс его изменения во времени — случайным процессом.

2. Модель отказа технической системы

Обозначим через x состояние системы и через X — множество всех возможных состояний. Исходя из принципа изоморфизма в моделировании [6], множество X можно считать отрезком [0,1].

Формальная модель динамики системы и ее отказа не должна противоречить естественным требованиям, которые подтверждены практикой:

а) вероятность отказа зависит от состояния, определяемого степенью износа;

б) с ухудшением состояния (увеличением степени износа) вероятность отказа увеличивается;

в) состояние системы монотонно не улучшается в процессе эксплуатации и хранения;

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

Выше мы предположили, что на множестве состояний X задано отношение порядка. Положим, что состоянию x=0 соответствует новая система. Вероятность отказа такой системы минимальна. ^стоянию x=1 соответствует полностью изношенная система, вероятность отказа ее на любом интервале времени равна 1.

85