close

Вход

Забыли?

вход по аккаунту

?

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

код для вставкиСкачать
Вестник ТГПИ
Естест венные науки
Раздел IV. Информатика
С.Г. Буланов
СХЕМА КОМПЬЮТЕРНОГО АНАЛИЗА УСТОЙЧИВОСТИ СИСТЕМ
ЛИНЕЙНЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С НЕЛИНЕЙНОЙ ДОБАВКОЙ
Исследования устойчивости по Ляпунову составляют актуальное направление в качественной теории дифференциальных уравнений. Для технических приложений представляется важным
обеспечить возможность компьютерного анализа устойчивости в режиме оперативного контроля
за протеканием моделируемого процесса. Существенно обеспечить мониторинг состояния устойчивости системы в режиме реального времени.
Рассматривается задача Коши для системы дифференциальных уравнений
dY
A ( t )Y
dt
Y ( t 0 ) Y0 ,
F ( t, Y ) ,
(1)
где A (t ) – матрица n n , F ( t , Y ) вектор-функция. Все коэффициенты A (t ) функции ai j (t )
i, j 1,... , n и элементы F ( t , Y ) определены, непрерывны и непрерывно дифференцируемы на
отрезке t0 , T при любом выборе T
const , T
t0 ,
.
n
Всюду ниже используется каноническая норма матрицы
A
1 i n
ная с ней норма вектора Y
a i k и согласован-
max
k 1
max y k .
1 k n
Предполагается, что
0 , такое, что для каждого невозмущенного и возмущенного ре-
шения (1) с начальными данными, удовлетворяющими неравенству
~
Y0 Y0
,
(2)
выполнены все условия существования и единственности решения на t0 ,
.
Требуется исследовать решение задачи (1) на устойчивость в смысле Ляпунова в области
~
~
R : { t0 t
; Y ( t ) , Y ( t ) : Y0 Y0
}.
(3)
Традиционное определение устойчивости заимствовано из [3] с упрощениями [1, 2] на случай рассматриваемой области R .
В данных условиях для T t0 ,
на отрезке t 0 , T выполнено [2]:
1) A ( t )
c , где c
2) пусть g k ( ti , Yi )
c1
const зависит от T ;
ak 1 ( ti ) y1(i ) ... ak n ( ti ) yn (i ) , тогда g k ( t , Y ( t ) )
~
const зависит от T и Y0 , k
~
c1 , g k ( t , Y ( t ) )
c1 , где
1,... , n .
Исследование системы (1) проводится в предположении, что устойчива соответствующая
ей линейная система
dY
A ( t )Y ,
dt
Y ( t 0 ) Y0 .
(4)
Имеет место
Теорема 1 [1, 2]. Чтобы в рассматриваемых условиях решение задачи (4) было устойчиво, необходимо и достаточно выполнение неравенства
118
Раздел IV.
Информатика
i
lim
i
для
t
t0 ,
(E
h A ( ti

))
cE
(5)
const
 0
. Решение асимптотически устойчиво тогда и только тогда, когда выполнено (5) и,
кроме того, имеет место соотношение
i
lim
i
(E
h A ( ti

))
(6)
0
 0
при t
.
Конструируемый метод анализа устойчивости опирается на разностную схему решения
системы (1). В качестве разностной схемы выбран метод Эйлера, который для системы (1) имеет
вид
Yi 1 Yi h A ( ti ) Yi h F ( ti , Yi ) , i 0 ,1,... ,
или
Yi
1
( E h A ( ti ) ) Yi
h F ( t i , Yi ) ,
(7)
где h – шаг разностной схемы, узлы которой здесь и всюду ниже предполагаются равноотстоящими; при любом выборе t
const , t
t0 ,
, h , i и t всегда предполагаются связанными со-
отношениями:
1 t0
, i 0 ,1,... .
i 1
Для возмущѐнного решения системы (1) разностная схема (7) примет вид
~
~
~
Yi 1 ( E h A ( ti ) ) Yi h F ( t i , Yi ) ,
~
~ ~
где Y ( t 0 ) Y0 , Y0 из (2).
t
ti 1 , h
ti
(8)
(9)
Соответствующие приближѐнным решениям (7), (9) точные решения могут быть представлены в форме метода Эйлера с остаточным членом на каждом шаге, остаточный член понимается
как векторная разность между точным решением и его приближением по методу Эйлера.
Иными словами,
~
~
~
~
( E h A ( ti ) ) Yi h F ( ti , Yi )
Yi 1 ( E h A ( ti ) ) Yi h F ( t i , Yi )
(10)
Ei ,
E i , Yi 1
~
где E i , E i – векторы остаточных членов
~
~
~
( E1i ,..., E n i ) , E i ( E1i ,..., E n i ) ,
Ei
компоненты которых описываются ниже.
Для определения и оценки компонент остаточного члена
dyk
k ( t i , Yi ) ,
dt
ak 1 ( ti ) y1 i
k ( t i , Yi )
Ei
запишем систему (1) подробно:
(11)
... ak n ( ti ) y n i
f k ( ti , Yi ) , k
1, ... , n .
Замечание 1. Условия накладываемые на систему (1) в совокупности с ограничениями, накладываемые на элементы матрицы A (t ) , обеспечивают то, что правая часть (11) определена, непрерывна, непрерывно дифференцируема на отрезке
t
t0 ,
Eki
при любом выборе t
const ,
.
Точное решение (11) определяется равенством
yk (i 1) yk i h k ( ti , Yi )
где
t0 , t
Eki
, k 1, ..., n ,
(12)
– остаточный член формулы Тейлора. Для уравнения с номером k остаточный член мож-
но представить в виде
119
Вестник ТГПИ
Естест венные науки
1
2
Eki
k ( ti
ki
h, Y ( ti
ki
h ) ) h2 , 0
1 , k 1, ... , n .
ki
(13)
Аналогично
~
~
~
~
1
2
k i 1 , k 1, ..., n .
k ( ti
k i h, Y ( t i
ki h))h , 0
2
В рассматриваемых ограничениях с учетом непрерывной дифференцируемости функций
~
k ( t , Y ( t ) ) , k 1, ... , n
k ( t , Y ) справедливо, что первые производные функций
k ( t, Y ( t ) ) ,
Eki
ограничены одной и той же константой на отрезке t0 , t :
~
c2 ,
c 2 , c2
k ( t, Y ( t ) )
k ( t, Y ( t ) )
при любом выборе t const , t
c2 ( t ) const
, для всех k 1, ... , n .
t0 ,
В этих условиях коэффициент перед h 2 в формуле (13) ограничен на отрезке
t0 , t
и
имеют место неравенства:
~
~
1
1
(14)
c2 h 2 ,
max E k i
c2 h 2 .
Ei
1 k n
2
2
С учетом введенных ограничений и предварительных утверждений ниже выводятся условия
устойчивости решения системы (1).
Из (10) разность между возмущѐнным и невозмущѐнным решением имеет вид точного равенства:
~
~
~
Yi 1 Yi 1 ( E h A ( ti ) ) ( Yi Yi ) h ( F ( ti , Yi ) F ( ti , Yi ) )
(15)
Ei ,
max
Ei
где
~
Ei
Ei
Ei
Eki
1 k n
~
. Поскольку
Ei
Ei
Ei
, то
c2 h 2 .
Ei
(16)
Для системы (1) рекуррентное преобразование (15) влечет равенство
i
~
Yi
Yi
1
~
( E h A( t i  ) ) ( Y0 Y0 ) DE i
1
SEi ,
(17)
 0
где
i
i r
DE i
E h A ti
~
h F t r 1 , Yr

F t r 1 , Yr
1
~
h F t i ,Y i
1
F t i , Yi
,
(18)
r 1  0
i r
i
SEi
( E h A ( ti  ) )
Er 1
Ei
,
(19)
r 1 0
где
Ei
имеет тот же вид, что в (13), (15).
Переходя к пределу в равенстве (17) при h
~
Y (t ) Y (t )
i
~
( E h A ( t i  ) ) ( Y0 Y0 ) lim DE i
lim
i
0 , что равносильно i
i
 0
, получим
lim S E i ,
(20)
i
На основании (16) и того, что норма матрица A(t ) ограничена на [ t0 , T ] для
T
t0 ,
,
имеет место
Лемма 1. В рассматриваемых условиях выполняется соотношение
SE i
O( h ) .
(21)
Доказательство. Произведение под знаком суммы в (19) оценивается из неравенства
i r
i r
( E h A ( ti  ) )
 0
120
( E h A ( ti  ) )
Er 1
 0
Er 1
, r 1, ... , i .
Раздел IV.
Информатика
С учѐтом ограниченности матрицы A по норме и неравенства (16) получим, что
i r
( E h A ( ti  ) )
c2 h 2 (1 ch ) i
Er 1
r 1
, r 1, ... , i .
 0
Отсюда
i 1
c2 h 2
SEi
i r 1
1 ch
,
r 1
или
i
c2 h 2
SEi

1 ch
.
 0
Суммирование геометрической прогрессии влечѐт
c2 h
(1 c h ) i
c
SEi
1
1 .
Отсюда с учѐтом (8) следует, что
ti
SEi
где c2
1
t0
h
c2 h 1 c h
1 ,
c2
.
c
Тем более, верно неравенство
c 2 h e c ( ti
SEi
При любом выборе t
ется соотношение
c2 e с ( t
ti 1 , t
t0 )
const , t
t0 )
1 .
(22)
, и при переменном i , h
t0 ,
0 , если
1 h
1
t
t0
i 1
, выполня-
h
0 . Отсюда и из (22) вытекает, что

0.
(23)
O ( h ) . Лемма доказана.
SE i
В условиях леммы 1 имеет место равенство
lim S E i
h
0
С учетом равенства (23) соотношение (20) примет вид
~
Y (t ) Y (t )
i
~
( E h A ( t i  ) ) ( Y0 Y0 ) lim DE i .
lim
i
i
 0
Предположим, что решение задачи (1) устойчиво.
~
0 , такое, что из неравенства Y0
Тогда для ~ 0
~
Y (t ) Y (t )
следует
Y0
~
t [ t0 , ) .
2
С учетом (24) и в соответствии с (25) с необходимостью выполняется неравенство
i
lim
i
для
t
t0 ,
~
( E h A( t i  ) ) ( Y0 Y0 )
 0
(24)
(25)
~
lim DE i
i
(26)
2
.
Из (26) следует
i
lim DEi
i
lim
i
 0
~
( E h A ( t i  ) ) ( Y0 Y0 )
~
2
(27)
121
Вестник ТГПИ
Естест венные науки
~
Выбирая
i
2 cE
, получим
i
Отсюда и из (27) следует, что
~
~
( E h A ( t i  ) ) ( Y0 Y0 )
lim
2
 0
lim DE i
i
~
~
2
2
~,
lim DE i
i
t
.
или
.
t0 ,
(28)
Условие (28) является необходимым условием устойчивости.
Докажем, что это условие достаточно для устойчивости решений системы (1).
Из (24) получаем lim DE i
i
~
Y ( t ) Y (t )
i
Выбирая в (2)
i
(E
h A ( ti

~
( E h A( t i  ) ) ( Y0 Y0 ) . Отсюда
 0
i
lim
i
~
Y ( t ) Y ( t ) lim
~
Y0
))
lim D E i ,
Y0
i
 0
~
2c E
t
t0 ,
.
(29)
, где cE из (13), получим, что величина возмущения в (29) оценива-
ется из неравенства
~
~
Y (t ) Y (t )
3~
,
2
~
2
Выбирая в соотношениях (25) – (30) ~
~
Y (t ) Y (t )
t
t0 ,
.
(30)
2
, получим, что
3
,
t [ t0 , ) .
(31)
0 указано
Таким образом, для произвольного
, такое, что из неравенства
~
Y0 Y0
, следует неравенство (31). Это по определению означает устойчивость невозмущен-
ного решения системы (1). Достаточность доказана.
Тем самым доказана.
Теорема 2. Чтобы в рассматриваемых условиях и при условии устойчивости системы (4)
решение задачи (1) было устойчиво, необходимо и достаточно, чтобы для произвольного ~ 0
~
0 такое, что при Y0 Y0
нашлось
выполнялось неравенство (28).
Следствие 1. В условиях теоремы 2 для устойчивости решения системы (1) необходимо и
достаточно чтобы выполнялось неравенство
~
Y ( t ) Y ( t ) lim
i
для
t
t0 ,
i
~
( E h A ( ti  ) ) ( Y0 Y0 )
~,
(32)
 0
.
Доказательство следствия вытекает из равенства (24) и утверждения теоремы 2.
Полученное условие устойчивости ориентировано на программную реализацию. Выражение, стоящее под знаком нормы в (32), программируется и, таким образом, можно выполнять анализ устойчивости на компьютере. На каждом шаге работы программы будут находиться разностные приближения возмущенного и невозмущенного решений и циклически накапливаться частичное матричное произведение. Через фиксированное (для удобства) количество шагов будет выводиться на экран текущее значение нормы из левой части неравенства (32). По асимптотическому
поведению этих значений нормы будет делаться вывод о характере устойчивости решения системы – неограниченный рост нормы будет свидетельствовать о неустойчивости, ограниченность
нормы соответствует устойчивости.
122
Раздел IV.
Информатика
Ниже приводится текст программы на языке программирования 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,y1,y2:extended):extended;
begin f1:=sin(t)*y1+cos(t)*y2-t*sin(y1); end;
function f2(t,y1,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,y1,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;
123
Вестник ТГПИ
Естест венные науки
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[l,j];
c[i,j]:=s;
end;
for i:=1 to n do
for j:=1 to n do
a[i,j]:=c[i,j];
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[i,j]:13);
writeln;
end;
readln; readln;
end.
Приближение к
~
Y (t ) Y (t )
i
lim
i
(E
~
h A ( ti  ) ) ( Y0
Y0 ) на промежутке моделиро-
 0
вания строится по шагам разностной схемы таким образом, что шаг h не стремится к нулю, а непосредственно фиксируется достаточно малым. Как показал программный и численный эксперимент [2], при этом сохраняется достоверность компьютерного анализа устойчивости.
Полученные условия позволяют определить характер устойчивости, асимптотической устойчивости либо неустойчивости систем ОДУ данного класса без представления решения в аналитической форме непосредственно по значениям разностных приближений. Мультипликативная
форма выражений под знаком предела в левой части предоставляет возможность запрограммировать вычисление этих выражений в виде цикла по числу сомножителей. Это влечет предпосылки
компьютерного анализа устойчивости в режиме реального времени без обращения к аналитическим методам качественной теории дифференциальных уравнений и к системам компьютерной
математики.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Буланов С.Г. Разработка и исследование методов программного моделирования устойчивости систем линейных дифференциальных уравнений на основе матричных мультипликативных преобразований разностных схем: дис. … канд. техн. наук. Таганрог: ТРТУ. 2006. 232 с.
2. Ромм Я.Е., Буланов С.Г. Метод компьютерного анализа устойчивости систем линейных дифференциальных уравнений / ТГПИ. Таганрог,2009. 119 с. ДЕП в ВИНИТИ 30.04.09, № 268. В 2009.
3. Чезари Л. Асимптотическое поведение и устойчивость решений обыкновенных дифференциальных
уравнений. М.: Мир, 1964. 478 с.
124
1/--страниц
Пожаловаться на содержимое документа