CC BY
22
Раздел IV. Информатика С.Г. Буланов
СХЕМА КОМПЬЮТЕРНОГО АНАЛИЗА УСТОЙЧИВОСТИ СИСТЕМ ЛИНЕЙНЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С НЕЛИНЕЙНОЙ ДОБАВКОЙ
Исследования устойчивости по Ляпунову составляют актуальное направление в качественной теории дифференциальных уравнений. Для технических приложений представляется важным обеспечить возможность компьютерного анализа устойчивости в режиме оперативного контроля за протеканием моделируемого процесса. Существенно обеспечить мониторинг состояния устойчивости системы в режиме реального времени.
Рассматривается задача Коши для системы дифференциальных уравнений
dY
— = A(t)Y + F(t, Y), dt (1)
[Y(t0)=Y0,
где A(t) - матрица nxn, F(J,Y) вектор-функция. Все коэффициенты A(I) функции (7) i,j = \,...,n и элементы F(J,Y) определены, непрерывны и непрерывно дифференцируемы на отрезке , Т при любом выборе Т = const, Те |(1. х
n
Всюду ниже используется каноническая норма матрицы II.1II = max / a J и согласован" " 1 <i<n^ 1 1 к=1
ная с ней норма вектора || Y || = max \
\<к<п
Предполагается, что 3 <5 > 0, такое, что для каждого невозмущенного и возмущенного решения (1) с начальными данными, удовлетворяющими неравенству
|Y0- Y0|<£, (2)
выполнены все условия существования и единственности решения на J0, со .
Требуется исследовать решение задачи (1) на устойчивость в смысле Ляпунова в области
R :{ to< t< оо; Y(t), Y (t): \\% - Yo||<£}. (3)
Традиционное определение устойчивости заимствовано из [3] с упрощениями [1, 2] на случай рассматриваемой области R .
В данных условиях для VT е [(|. х на отрезке / _ выполнено [2]:
1) (/) || < с , где с = const зависит от Т ;
2) пусть gk (t,, Y, ) = aki(t,) уио-к..+ akn (t,) yn(i), тогда \g'k (t, Y( t ))|< Ci, \g'k (t, Y(t ))|< Ci, где
с, = const зависит от T и Y0, к = 1.....п .
Исследование системы (1) проводится в предположении, что устойчива соответствующая ей линейная система
= A{t)Y,
dt (4)
Y (to) = Yo. Имеет место
Теорема 1 [1, 2]. Чтобы в рассматриваемых условиях решение задачи (4) было устойчиво, необходимо и достаточно выполнение неравенства
limQ^ + A^,))
' 1=0
< сЕ = const
(5)
для V/ е [0, оо . Решение асимптотически устойчиво тогда и только тогда, когда выполнено (5) и. кроме того, имеет место соотношение
I
Ига \\{Е + кА{1г_,))
' е=о
->0
(6)
при t —» СЮ .
Конструируемый метод анализа устойчивости опирается на разностную схему решения системы (1). В качестве разностной схемы выбран метод Эйлера, который для системы (1) имеет вид
YI+1 =Y1 +hA(tj)Yj + hF(t„ Yj) , i = 0,1,...,
или
Yi+\ = (E+ hA(ti))Yi +hF{ tt, Yi), (7)
где h - шаг разностной схемы, узлы которой здесь и всюду ниже предполагаются равноотстоящими; при любом выборе t = const. / е [0. х h . i и / всегда предполагаются связанными соотношениями:
t = tM, h = tl+l , / = 0,1,....
/+1
Для возмущённого решения системы (1) разностная схема (7) примет вид
(8) (9)
где У(/0) = У0, Y0 из (2).
Соответствующие приближённым решениям (7), (9) точные решения могут быть представлены в форме метода Эйлера с остаточным членом на каждом шаге, остаточный член понимается как векторная разность между точным решением и его приближением по методу Эйлера.
Иными словами,
Ъ+1=(Е + ИА(^1+ИР(Г1,^) + ®Е1,%+1=(Е + кА(0)%+11Е((1,%) + ©Е1, (10) где & ¡г,. (-) ¡,:, - векторы остаточных членов
®Е1={вш,...,вЕт), ©Я1.=(3
E\i>
дЕш ) >
компоненты которых описываются ниже.
Для определения и оценки компонент остаточного члена ®Е{ запишем систему (1) подробно:
% = (U, Y), dt
(11)
Ф* (tг, Y ) = ak 1 (ti ) Уи+ - + akn (ti ) Упг + fk (ti, Y X * = 1,..., Замечание 1. Условия накладываемые на систему (1) в совокупности с ограничениями, накладываемые на элементы матрицы А (I), обеспечивают то, что правая часть (11) определена, непрерывна, непрерывно дифференцируема на отрезке [0, t _ при любом выборе t - const,
te Io,aC
Точное решение (11) определяется равенством
Ук(,+1) =Укг +h<5>k(ti,Yi) + 6Eki, к = \,...,п, (12)
где 0Eki - остаточный член формулы Тейлора. Для уравнения с номером к остаточный член можно представить в виде
0EkI=-^k(tI+SklKY{tl+Sklh))h\ 0<Ski<l,k = \,...,n. (13)
Аналогично
6Ekl=-®k(.tI+SklKY(tl+Sklh))h2, О < 6ki < 1, к = \ ,...,п.
В рассматриваемых ограничениях с учетом непрерывной дифференцируемости функций
Ф*(/,7) справедливо, что первые производные функций Ф/(/. Yd)). (/. Yd)). к = 1.....п
ограничены одной и той же константой на отрезке [0, / _:
|фк (t, Y(t))| <C2, | Фк (t, ©(t))|< C2, C^ C2(t) = const
при любом выборе t= const, t g [(i. x . для всех к = 1.....n.
В этих условиях коэффициент перед h2 в формуле (13) ограничен на отрезке |0, t _ и имеют место неравенства:
II ©я, II = max | вЕк,\<1 C2h 2, || © e J = max | U1C2h 2. (14)
11 11 1 <k<nf 1 2 II II i<i<„l I 2
С учетом введенных ограничений и предварительных утверждений ниже выводятся условия устойчивости решения системы (1).
Из (10) разность между возмущённым и невозмущённым решением имеет вид точного равенства:
YI+1-YI+1=(E + hA(tI))(YI-YI)+h(F(tI,YI)-F(tI,YI)) + &El, (15)
где ®Ei= ©Ei~ ®Ei . Поскольку || ©e^ || ©Ei || + || ®Ег 11 , то
||©E^<C2 h 2. (16)
Для системы (1) рекуррентное преобразование (15) влечет равенство
i
YM-YM =Y\(E + hA(tl_l))(Y0-Y0)+DEl +SEi, (17)
1=0
где
= if\<P + hA 3 ^ ir-iJr-1 IjF tr_uYr_, J A i ijt X (18)
r=l I-0
i i-r
SZ' = X П( ^^^ ( »®Er-1 + ®Ei , (19)
r = 1 * = 0
где &Ei имеет тот же вид, что в (13), (15).
Переходя к пределу в равенстве (17) при h —> 0 , что равносильно г —» оо , получим
i
Y(t)-Y(t) = \ivaT~\(E + hA( t,_e)) (f0 - 70) + lim DEl + lim SEl, (20)
г—^oo-1- 1 i—>oo г—>oo
¿=0
На основании (16) и того, что норма матрица A(t) ограничена на | /(|. Г] для VT е |(|. со , имеет место
Лемма 1. В рассматриваемых условиях выполняется соотношение
|М| = °( h). (21)
Доказательство. Произведение под знаком суммы в (19) оценивается из неравенства
fj( E + hA (tt_ * ))© Er_ i
*= 0
i—r
<
]~f||(E + hA(tt_*))||-||®ег-i||, r = 1,-, i•
Ы 0
С учётом ограниченности матрицы А по норме и неравенства (16) получим, что
е=о
<c2h2(l+ch У~г+\ г = 1,..., /'.
Отсюда
\sEi\<c2h2 j] < + ch
r=1
или
|SEi|<C2 h2 £ « + Ch Z .
t = 0
Суммирование геометрической прогрессии влечёт
Ст h
sA<^lL{i+chy+l -\ .
Отсюда с учётом (8) следует, что
SEi\\ ^ c2h
ch
-1
где Ст = —.
с
Тем более, верно неравенство
\\SE1 II< c2hisc{ti+l~4)
-1 .
(22)
t-tr
При любом выборе t = /,_, , t = const. / е [(1. х , и при переменном i. h =-- , выполня-
/' +1
ется соотношение c9
V(t-tо)
>
-1 h —> 0, если h —> 0 . Отсюда и из (22) вытекает, что
| \ = О(й). Лемма доказана.
В условиях леммы 1 имеет место равенство
lim SEj = 0 .
/;-> О
С учетом равенства (23) соотношение (20) примет вид
i
Y (t) - Г(0=1ппП(/-.+/7.1( ilW)) (% - Ii,) + lim DEl.
i—Ьсп-*- -L- i —±cr>
t=0
Предположим, что решение задачи (1) устойчиво.
Тогда для V~ > 0 33 > 0, такое, что из неравенства || Y0 - Y01| < 8 следует
||f(0-}'(0||^| V/ е[/0, со). С учетом (24) и в соответствии с (25) с необходимостью выполняется неравенство
г
lim ГГ (E + hA( tt_,))(~o - Yo) + lim DE
i—J—^cf)
Ы0
s < —
для
Vt 6 |0, » . Из (26) следует
lim I),
i
ИтГГ(£ + /7,4(^))(}0-}0)
/'-> м-1 1
f = 0
e < —
(23)
(24)
(25)
(26)
(27)
4+1 to
h
С
Выбирая 3 =
2 св
получим
ИшГТ +
Отсюда и из (27) следует, что
1ш1/)с
г
< —
8
< — ИЛИ 2
Ит Б
Ег
< е, У/е [0
(28)
Условие (28) является необходимым условием устойчивости.
Докажем, что это условие достаточно для устойчивости решений системы (1).
г
Из (24) получаем Ит£>Вг = 7(0-7(0 - ПтГТ(/•; + /?. ((;, , ))(70 -70). Отсюда
7 —^ СП 7 —^ СП -1- А
7 ( г )-7 (г ) <
ИшГТ + (*,._,))
7 -^ СП -А-
¿ = 0
0
7„- 7о +
Ит£>
Ег
V/ е ^ 0, оо
(29)
Выбирая в (2) 3 =-, где сЕ из (13), получим, что величина возмущения в (29) оценива-
2 с„
ется из неравенства
|7(0-7(0||<| + £=у, V/«
2
(30)
Выбирая в соотношениях (25) - (30) е = — г, получим, что
|| 7 (г)-7 (г )||<£, V/ е [¿о, °0). (31)
Таким образом, для произвольного е > 0 указано <5, такое, что из неравенства || 70 -70||< 3, следует неравенство (31). Это по определению означает устойчивость невозмущенного решения системы (1). Достаточность доказана. Тем самым доказана.
Теорема 2. Чтобы в рассматриваемых условиях и при условии устойчивости системы (4) решение задачи (1) было устойчиво, необходимо и достаточно, чтобы для произвольного е >0
нашлось 3 >0 такое, что при ||70 - 701| < 3 выполнялось неравенство (28).
Следствие 1. В условиях теоремы 2 для устойчивости решения системы (1) необходимо и достаточно чтобы выполнялось неравенство
7(0-7(0-ИтТТ (Е + ЬА((,_г))(¥0-¥0)
7 —>00-*- 1
< //
(32)
для V/ е [
Доказательство следствия вытекает из равенства (24) и утверждения теоремы 2.
Полученное условие устойчивости ориентировано на программную реализацию. Выражение, стоящее под знаком нормы в (32), программируется и, таким образом, можно выполнять анализ устойчивости на компьютере. На каждом шаге работы программы будут находиться разностные приближения возмущенного и невозмущенного решений и циклически накапливаться частичное матричное произведение. Через фиксированное (для удобства) количество шагов будет выводиться на экран текущее значение нормы из левой части неравенства (32). По асимптотическому поведению этих значений нормы будет делаться вывод о характере устойчивости решения системы - неограниченный рост нормы будет свидетельствовать о неустойчивости, ограниченность нормы соответствует устойчивости.
2
0
Ниже приводится текст программы на языке программирования Object Pascal в среде Delphi 7: program nelinkriter; {$APPTYPE CONSOLE} uses
SysUtils; const
h=0.00001; n=2; t0=0; TT=1000; eps1=0.0000111; eps2=0.0000111; y01=0; y02=0; yv01=y01+eps1; yv02=y02+eps2; type
matr=array[1..n,1..n] of extended; stolb=array[1..n] of extended; var
A,B,C:matr; delta,yv,ym:stolb; t,s,norma,norma1,norma2:extended; y1,y2,y11,yv1,yv2,yv11:extended; i,j,l:integer; k:longint; procedure matrprav (t:extended; var A:matr); begin
a[1,1]:=sin(t); a[1,2]:=cos(t); a[2,1]:=-cos(t); a[2,2]:=sin(t); end;
function f1 (t,y 1 ,y2 :extended) :extended; begin f1:=sin(t)*y1+cos(t)*y2-t*sin(y1); end; function f2(t,y 1 ,y2 :extended) :extended; begin f2:=-cos(t)*y1+sin(t)*y2-sqr(y1)*y2; end; procedure metod_E (var A:matr); begin
matrprav (t,A); for i:=1 to n do for j:=1 to n do
begin a[i,j]:=h*a[i,j]; if i=j then a[i,j]:=a[i,j]+1; end;
end;
begin
t:=t0; k:=0; metod_E (A); y1:=y01; y2:=y02; yv1:=yv01; yv2:=yv02; delta[ 1] :=yv01-y01; delta[2] :=yv02-y02; repeat
y11:=y1; y1:=y1 +h*f1 (t,y 1,y2); y2:=y2+h*f2(t,y11,y2);
yv11:=yv1; yv1:=yv1+h*f1(t,yv1,yv2); yv2:=yv2+h*f2(t,yv11,yv2);
yv[1]:=yv1-y1; yv[2]:=yv2-y2;
if abs(yv[1]/delta[1])>abs(yv[2]/delta[2]) then
norma1:=abs(yv[1]/delta[1]) else norma1:=abs(yv[2]/delta[2]);
if abs(yv[1])>abs(yv[2]) then norma2:=abs(yv[1]) else
norma2:=abs(yv[2]);
for i:=1 to n do
begin
ym[i]:=0;
for j:=1 to n do ym[i]:=ym[i]+(a[i,j]*delta[j]); end;
for i:=1 to n do yv[i]:=yv [i]-ym[i]; norma:=abs(yv[1]); for i:=1 to n do
if abs(yv[i])>norma then norma:=abs(yv[i]); t:=t+h;
metod_E (B); for i:=1 to n do for j:=1 to n do begin
s:=0;
for l:=1 to n do s:=s+a[i,l] *b [lj];
c[ij]:=s;
end;
for i:=1 to n do for j:=1 to n do a[ij]:=c[ij];
k:=k+1; if k>=5000000 then begin
writeln('t=',t:4:0,' ','norma=',norma:12,' 'norma1=',norma1:12,' ','norma2=',norma2:12); k:=0; end;
until t>=TT;
writeln; writeln;
{Вывод конечной матрицы}
for i:=1 to n do
begin
for j:=1 to n do write(' ',a[ij]:13); writeln; end;
readln; readln; end.
Приближение к
~(t)- Y(t)- lim TT(E + hA(tt_e ))(Y0 " Y0)
i —100 1
1=0
на промежутке моделиро-
вания строится по шагам разностной схемы таким образом, что шаг й не стремится к нулю, а непосредственно фиксируется достаточно малым. Как показал программный и численный эксперимент [2], при этом сохраняется достоверность компьютерного анализа устойчивости.
Полученные условия позволяют определить характер устойчивости, асимптотической устойчивости либо неустойчивости систем ОДУ данного класса без представления решения в аналитической форме непосредственно по значениям разностных приближений. Мультипликативная форма выражений под знаком предела в левой части предоставляет возможность запрограммировать вычисление этих выражений в виде цикла по числу сомножителей. Это влечет предпосылки компьютерного анализа устойчивости в режиме реального времени без обращения к аналитическим методам качественной теории дифференциальных уравнений и к системам компьютерной математики.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Буланов С.Г. Разработка и исследование методов программного моделирования устойчивости систем линейных дифференциальных уравнений на основе матричных мультипликативных преобразований разностных схем: дис. ... канд. техн. наук. Таганрог: ТРТУ. 2006. 232 с.
2. Ромм Я.Е., Буланов С.Г. Метод компьютерного анализа устойчивости систем линейных дифференциальных уравнений / ТГПИ. Таганрог,2009. 119 с. ДЕП в ВИНИТИ 30.04.09, № 268. В 2009.
3. Чезари Л. Асимптотическое поведение и устойчивость решений обыкновенных дифференциальных уравнений. М.: Мир, 1964. 478 с.
