CC BY
55
УДК 517.929
А.Ф. Зубова, А.В. Зубов, В.И. Зубов, М.В. Стрекопытова
ИССЛЕДОВАНИЕ СИСТЕМЫ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ НА СХОДИМОСТЬ, УСТОЙЧИВОСТЬ И ТОЧНОСТЬ
A.F. Zubova, A.V. Zubov, V.I. Zubov, M.V. Strecopitova
THE INVESTIGATION OF SYSTEMS DIFFERENTIAL EQUATIONS ON INTIMATION, STABILITY AND PRECISION
В статье дан критерий устойчивости линейных систем, а также приведены методы интегрирования матричных дифференциальных уравнений высокой размерности.
МАТРИЦА. УРАВНЕНИЕ. РЕШЕНИЕ. АЛГОРИТМ. ПРОИЗВОДНАЯ. КОРЕНЬ. СИСТЕМА. ПОРЯДОК. ЧИСЛО.
In giving article is give criteria of stability linear systems, and so is bring methods of integration matrix differential equations higher dimension.
MATRIX. EQUATION. SOLUTION. ALGORITHM. DERIVATIVE. ROOT. SYSTEM. ORDER. NUMBER.
При построении вычислительных алгоритмов в случаях асимптотической устойчивости или условной асимптотической устойчивости решения первоначальной вычислительной задачи (решение специально построенной системы дифференциальных уравнений) сходимость и устойчивость алгоритма следуют автоматически, а о точности можно судить по оценке нормы приближения решения специально построенной дифференциальной системы к положению равновесия. Под термином «специально построенная система дифференциальных уравнений» следует понимать такую систему, которая имеет те же корни, что и первоначальная, но порядок которой установлен намного ниже, для того чтобы ее было легче исследовать на сходимость, устойчивость и точность.
Постановка задачи
1. Мы установили, что все решения матричного дифференциального уравнения
X = -А*АХ + А* (1)
сходятся к значению обратной матрицы А-1, если она существует, т. е. если матрица А невы-
рожденная. При этом матричная функция V = = АХ — Е удовлетворяет на решениях системы (1) дифференциальному уравнению
V = -АА>. (2)
Откуда следует экспоненциальная оценка скорости приближения решения Х(^ ,Х0) к значению обратной матрицы А-1. Интегрируя систему (1) методом Эйлера с шагом И , получим итерационный алгоритм
Х к+1 = Х к + И(- Аахк + А*).
Поскольку система (1) линейная, то ее решения легко разложить в ряд и получить эффективный алгоритм Шульца [1] или — представляя этот ряд в виде бесконечного произведения — еще более эффективный алгоритм автора [2]. Использование методов интегрирования более высокого порядка для системы (1) неэффективно.
2. Мы показали также, что все решения нелинейного матричного дифференциального уравнения
Х = Х - ХАХ (3)
также сходятся к значению обратной матрицы А-1, если она существует. Матричная функция
V = X 1 - А удовлетворяет на решениях (3) дифференциальному уравнению
V = —У, (4)
откуда следует оценка величины интервала интегрирования. Интегрируя систему (3) методом Эйлера, получаем итерационный алгоритм
Хк+1 = хk + КХк - х^
который при h = 1 принимает вид
Х k+1 = 2Х k - Х кахк
и называется алгоритмом Шульца — Джонса — Майера ^М). Алгоритм сходится не всегда и требует специальной подготовки начальных данных [3, 4]. Причина расходимости в некоторых случаях алгоритма SJM заключается только в ограниченности машинной арифметики, так как из (4) следует, что он должен сходиться при любых начальных данных. Действительно, представляя решение системы (3) в виде ряда Тейлора при h = 1
X (1) = X (0) + *(0) - Х^)АХ (0) +
X (0) - 3Х (0)АХ (0) + 2Х (0) АХ (0)АХ (0) + 2! +(5)
видим, что присутствующий в силу нелинейности системы (3) множитель АХ(0) может превалировать некоторое время над фактором в знаменателе и, более того, вызвать ошибку арифметического переполнения. Это может происходить только при условии
II АХ(0) ||>> 1.
Отсюда следует и очевидный критерий выбора начального приближения. Метод SJM еще плохо отражен в научной литературе, и впереди — большие численные эксперименты и исследования по его сравнению с другими алгоритмами. Следует отметить также, что в силу нелинейности системы (3) применение неявных или полуявных методов невозможно.
3. Мы установили также, что решение Х(^,Х0) при Х0 = А матричного дифференциального уравнения
Х=Х+Х *-1) 2
(6)
А = Би, Б: {V х еЕл : х Ф 0, х Бх > 0}. Матричная функция V = XX - Е удовлетворяет на решениях системы (6) системе уравнений (4).
4. Мы установили также, что решение Х(^ ,Х0) при Х0 = А*-1 матричного дифференциального уравнения
Х=1(-Х+ХХ*Х) 2
(7)
также сходится к значению ортогональной матрицы и полярного разложения матрицы А. Матричная функция V = (ХХ*)-1 - Е удовлетворяет на решениях системы (7) системе уравнений (4).
5. Нами было показано, что все решения матричного дифференциального уравнения
Х=-Х+Х АА
(8)
сводится к значению ортогональной матрицы
*
и: и и = Е полярного разложения матрицы
сходятся к значению матричной функции ква-
*
дратного корня из матрицы АА при условии, что матрица А — невырожденная. При этом матричная функция V = (1 /2)(Х2 - АА*) удовлетворяет на решениях системы (8) системе уравнений (4) [3]. Из матричных дифференциальных систем (1), (3), (7), (8) только первая линейна, и поэтому для нее годится разложение решения в ряд, который позже представляется в виде бесконечного произведения. Интегрирование остальных систем осуществляется приближенными методами. Причем выбор величины h шага интегрирования, порядка метода и интервала интегрирования осуществляется на основе уравнения (4) и требуемой точности е. Матричная функция V является мерой отклонения решения матричного дифференциального уравнения от решения первоначальной вычислительной задачи. Эта функция удовлетворяет линейному матричному дифференциальному уравнению (4) на решениях построенных дифференциальных матричных систем.
Численный анализ систем большой размерности
Многие системы, с которыми работает исследователь, имеют высокую размерность и большие числа обусловленности. Поэтому не всегда целесообразно приводить систему линейного приближения к нормальной форме. Рассмотрим систему линейных дифференциальных уравнений второго порядка
А0Х + ИА1^ + А2Х = 0:
(9)
где Ао, А1, А2 — «обычные» матрицы, т. е. не имеющие больших чисел обусловленности и не представляющие затруднений в вычислительном аспекте; Н — большой положительный параметр. Система (1) описывает механическую систему, в которой присутствуют быстро вращающиеся массы. Попытаемся отыскать собственные числа уравнений (1). Решение системы (1) мы будем искать в виде
X = ехр(^)С ,
(10)
где С — постоянный вектор. Для существования у системы (1) ненулевого решения вида (2) необходимо и достаточно, чтобы выполнялось соотношение
ёе1(А0 + ИА^ + А2) = 0.
(11)
Будем искать корни последнего уравнения, имеющие порядок величины Н. Для этого в (3) сделаем замену ^ = |И. Получим
ёе1(А0|2И2 + А1цИ2 + А2) = 0.
Полагая а = И
имеем
ёе1;(А0|Г + А1ц + аА2) = 0. (12)
Пусть это уравнение при а = 0 имеет ненулевой корень |0 . Тогда |0 удовлетворяет соотноше-
нию
ёе1(А0|0 + А1) = 0.
(13)
|■ (а) = |0 +ац1 +а2|2 +.
(14)
Получим дД(|0., а)
= |0
д ёе1(А0^0. + А1)
д| ■ д| так как |0. — простой корень уравнения (13).
Из этого следует, что при достаточно малых | а | существует неявная функция |. (а), удовлетворяющая уравнению (12) и условию | . (0) = |0. . Разложимость функции | . (а) в виде ряда (14) следует из теоремы Коши, так как | . (а) удовлетворяет дифференциальному уравнению
■ = ЭД / дД с1 а да д| с голоморфной правой частью в окрестности точки |0 . Следовательно, | . (а) представима
в виде ряда (14), сходящегося при | а |< а0 . Доказательство закончено. Отметим, что корню | . (а) отвечает собственное число а■ =|■ (а)И . Теперь будем искать корни порядка Н—1. Для этого положим в (11) ^ = уИ-1 . Получим уравнение
Ае1(ау 2А0 + уА1 + А2) = 0. (15)
Пусть это уравнение при а = 0 имеет отличный от нуля корень у0 , тогда у0 удовлетворяет соотношению
ёе1(у0А1 + А2) = 0
(16)
В случае, когда матрица Ао невырожденная, число |0 является собственным значением матрицы А-1А1. Обозначим через |01, |02, ..., |0к корни уравнения (13).
Теорема 1. Каждому корню |0 , кратность которого равна единице, отвечает корень | ■ уравнения (12), определяемый по формуле
Обозначим через у01,у02,..., у0п корни этого уравнения.
Теорема 2. Каждому корню у0 ■ кратности единица уравнения (16) отвечает корень уравнения (15) V■ , определяемый по формуле
V.(а) = у0■ + +а у2■ +...
(17)
Ряд (14) сходится при достаточно малых значениях а .
Доказательство. Уравнение (12) определяет неявную функцию |. (а), удовлетворяющую условию | .■ (0) = |0 . Обозначим левую часть (12)
Д(|0■, а).
При этом ряд сходится при достаточно малых значениях а .
Доказательство. Уравнение (15) определяет неявную функцию у(а). Воспользуемся теоремой о существовании неявной функции. Обозначим левую часть (15) через Д^у, а). Найдем частную производную этой функции по V при а = 0:
ЭД1(у0 , ,0) д ёе1:(у0, ,А1 + А2)
Эу
Эу
2
так как у0 — простой корень. Поэтому при малых а существует неявная функция у^ (а). Возможность разложения неявной функции у ^ (а) в ряд (17) можно доказать аналогично теореме 1.
Доказательство закончено.
Если матрица А1 невырожденная, то у0 является собственным числом матрицы -А-1А2 . На основе теорем 1 и 2 могут быть получены конструктивные методы построения собственных значений ц^ у : для этого нужно ряды (14) и (17) подставить соответственно в (12) и (15) и приравнять коэффициенты при одинаковых степенях а . Отметим, что если матрицы А^ А1, А2 невырожденные и матрицы А-^, А-1А2 не имеют кратных собственных чисел, то первая и вторая группы корней содержат по п членов, где п — число уравнений в системе (1). Для нахождения собственных чисел ц0 , у0 потребу-
01 01
ется обращать матрицы А и А1. Здесь можно воспользоваться эффективным методом обращения матрицы. Изложенный подход можно применить к нахождению собственных чисел матрицы вида А + НВ, где Н — большой положительный параметр. Полагая X = Нц , получим, что число ц должно удовлетворять соотношению
ёй(А + НВ - НцЕ) = 0.
Положим в = Н-1, тогда получим соотношение
ёйфА + В-цЕ) = 0. (18)
Пусть ц0 — отличное от нуля собственное число матрицы В. Тогда соотношение (18) определяет неявную функцию ц = ц(в), удовлетворяющую условию ц(0) = ц0 . Если ц0 — простой корень, то выполнено условие
Э ёе1(В -ц0Е) = 0 Эц
и при малых в существует неявная функция ц(в), удовлетворяющая соотношению (18) и определяемая соотношением типа (14). В случае, когда кратность корня ц0 равняется к >1 , введем новую переменную
П = (ц-ц0)к. (19)
Тогда производная левой части (18) по переменной п при в = 0 в точке п = 0 отлична от нуля.
Следовательно, при малых | в | существует неявная функция п = П (в), определяемая рядом типа (14). Функцию ц(в) найдем через функцию п из формулу (19).
Достаточные условия устойчивости линейных систем
Поставим задачу определения и анализа корней уравнения
ёе1:(А -ХЕ) = 0, (20)
где А - (п х п)-матрица. Известно, что вопрос об устойчивости линейной системы обыкновенных дифференциальных уравнений будет решен, если выяснится, что все корни уравнения (20) лежат в левой полуплоскости комплексного переменного X . Следуя [4], отобразим эту левую полуплоскость на единичный круг в плоскости переменного р при помощи преобразования
Х = р=! (21)
р+1
При этом корни уравнения (20) перейдут в корни уравнения
ёе1:(В -рЕ) = 0, (22)
где В = Е + 2(А - Е)-1. Обозначим корни этого уравнения через Р1, р2, ..., рп. Из нашего преобразования следует, что если корни уравнения (20) лежат в левой полуплоскости, то корни уравнения (22) лежат в единичном круге. А тогда,
поскольку матрица Вк имеет собственные чис-
к к к ла Р1, р2, ..., рп , необходимым и достаточным
условием асимптотической устойчивости будет
условие
Вк ^ 0. (23)
к ^^
*
Рассмотрим теперь матрицу В*В . Эта матрица — симметричная, ее собственные числа неотрицательны. Обозначим собственные числа
*
матрицы В В через ц1, ..., цп. Пусть ц1 — наибольшее из этих чисел. Заметим, что из условия IР11< 1 (1 = 1,..., п) не следует | цг1< 1 (г = 1,..., п). При выполнении условия | рг1< 1 (г = 1,..., п) существует целое положительное число к такое, что собственные числа матрицы (Б*Б)к будут меньше единицы. Обратное также верно: если при некотором целом положительном к соб-
ственные числа матрицы (В В) меньше единицы, то собственные числа матрицы В лежат в единичном круге, т. е. имеет место асимптотическая устойчивость. По определению
|| B ||= max^J-^
XeE„ X|
= max-
XeE„
<JX*B*BX
|X|
Введем в рассмотрение величины
Sk = Sp ((B*B)k).
Теорема 3. Если при некотором целом положительном к выполнено неравенство
Sk 1 — <—,
(24)
то собственные числа матрицы В лежат в единичном круге. Следовательно, система с матрицей А асимптотически устойчива.
Доказательство. Вычислим отношение
Бк / ¿1:
^ =^к + |2 +... + |П к-1 (1 + (й-2 /+...)
¿1 11 +|2 +... + 1п 1 (1 + (|2/11) +...) ' Поскольку | — наибольшее собственное число (здесь мы предполагаем, что матрица В В имеет хотя бы одно ненулевое собственное число), то справедливы оценки
..к-1 S
Ml < Sl <.к-1.
n S1
(25)
Отсюда видно, что при выполнении условия (24) будет справедливо неравенство |к-1 < 1. Следовательно, наибольшее собственное значение
*
матрицы В В меньше единицы, а тогда выполнено условие (23), откуда следует, что собственные числа матрицы В лежат в единичном круге. Доказательство закончено.
Теорема является критерием устойчивости
линейных систем. Как отмечено выше, даже при
*
условии | р^ |< 1 матрица В В может иметь собственные числа, превосходящие единицу. Для того чтобы охватить и эти случаи, следует рассматривать вместо матрицы В*В матрицу (В*В)т, где т > 1. Условия теоремы останутся прежними. Для большей эффективности вычислений можно сразу рассматривать матрицу (В*В)т , где т достаточно велико. Следует отметить также, что настоящий критерий является достаточным. Необходимо добавить также, что эффективность предложенных алгоритмов нисколько не снижается в том случае, когда определитель матрицы А близок к нулю.
Работа выполнена при финансовой поддержке РФФИ (пр. № 10-08-000624).
СПИСОК ЛИТЕРАТУРЫ
1. Зубов, А.В. Динамическая безопасность управляемых систем [Текст] / А.В. Зубов, Н.В. Зубов.— СПб.: Изд-во НИИ Химии СПбГУ,— 2009.— 172 с.
2. Зубов, А.В. Математические методы качественного анализа систем управления и устойчивость расчетных движений [Текст] / А.В. Зубов, С.В. Зубов.— СПб.: Изд-во ВВМ, 2011.— 323 с.
3. Зубов, А.В. Управление динамическими системами [Текст]/ А.В. Зубов, О.А. Шабурова.— СПб.: Изд-во НИИ химии СПбГУ, 2005.— 83 с.
4. Зубов, И.В. Анализ управляемых систем и равновесных движений [Текст]/ И.В. Зубов, Н.В. Зубов, М.В. Стрекопытова.— СПб.: Изд-во ВВМ. 2012.— 322 с.
5. Зубова, А.Ф. Математические методы моделирования промышленных процессов и технологий [Текст] / А.Ф. Зубова.— СПб.: Изд-во СПбГУ,— 2004.— 472 с.
6. Зубов, И.В. Анализ систем управления и способы представления программных управлений [Текст] / И.В. Зубов.— СПб.: Изд-во ВВМ, 2014.— 150 с.
7. Стрекопытов, С.А. Теория квазипериодических систем [Текст] / С.А. Стрекопытов.— СПб.: ВВМ, 2014.— 157 с.
8. Стрекопытова, М.В. Анализ равновесных движений [Текст] / М.В. Стрекопытова.— СПб: Изд-во СПбГУ, 2014.— 176 с.
REFERENCES
1. Zubov A.V., Zubov N.V. Dinamicheskaia bezopas-nost' upravliaemykh sistem [Tekst].— SPb.: Izd-vo NII Khimii SPbGU, 2009.— 172 s. (rus.)
2. Zubov A.V., Zubov S.V. Matematicheskie metody kachestvennogo analiza sistem upravleniia i ustoichivost' raschetnykh dvizhenii [Tekst]— SPb.: Izd-vo VVM, 2011.— 323 s. (rus.)
3. Zubov A.V., Shaburova O.A. Upravlenie dinami-cheskimi sistemami [Tekst].— SPb.: Izd-vo NII Khimii SPbGU, 2005.— 83 s. (rus.)
4. Zubov I.V., Zubov N.V., Strekopytova M.V. Analiz upravliaemykh sistem i ravnovesnykh dvizhenii [Tekst].— SPb.: Izd-vo VVM, 2012.— 322 s. (rus.)
5. Zubova A.F. Matematicheskie metody modelirova-niia promyshlennykh protsessov i tekhnologii [Tekst].—
SPb.: Izd-vo SPbGU, 2004.— 472 s. (rus.)
6. Zubov I.V. Analiz sistem upravleniia i sposoby pred-stavleniia programmnykh upravlenii [Tekst].— SPb.: Izd-vo VVM, 2014.— 150 s. (rus.)
7. Strekopytov S.A. Teoriia kvaziperiodicheskikh sistem [Tekst] / SPb.: Izd-vo VVM, 2014.— 157 s. (rus.)
8. Strekopytova M.V. Analiz ravnovesnykh dvizhenii [Tekst].— SPb: Izd-vo SPbGU, 2014.— 176 s. (rus.)
СВЕДЕНИЯ ОБ АВТОРАХ
ЗУБОВА Александра Федоровна — доктор технических наук профессор кафедры теории управления факультета прикладной математики — процессов управления Санкт-Петербургского государственного университета; 191024, ул. Конная, д. 22/5, кв. 4, Санкт-Петербург, Россия; e-mail: ddemidova@ mail.ru
ЗУБОВ Афанасий Владимирович — доктор физико-математических наук профессор кафедры теории управления факультета прикладной математики — процессов управления Санкт-Петербургского государственного университета; 191024, ул. Конная, д. 22/5, кв. 4, Санкт-Петербург, Россия; e-mail: ddemidova@mail.ru
ЗУБОВ Всеволод Иванович — аспирант кафедры теории управления факультет прикладной математики — процессов управления Санкт-Петербургского государственного университета; 191024, ул. Конная, д. 22/5, кв. 4, Санкт-Петербург, Россия; e-mail: ddemidova@mail.ru
СТРЕКОПЫТОВА Мария Владимировна — аспирант кафедры теории управления факультет прикладной математики — процессов управления Санкт-Петербургского государственного университета; 191024, ул. Конная, д. 22/5, кв. 4, Санкт-Петербург, Россия; e-mail: ddemidova@mail.ru
AUTHORS
ZUBOVA Alexandra F. — St. Petersburg State University; 191024, Connya St. 4, St. Petersburg, Russia; e-mail: ddemidova@mail. ru
ZUBOV Afanasyi V. — St. Petersburg State University; 191024, Connya St. 4, St. Petersburg, Russia; e-mail: ddemidova@mail.ru
ZUBOV Vsevolod I. — St. Petersburg State University; 191024, Connya St. 4, St. Petersburg, Russia; e-mail: ddemidova@mail.ru
STRECOPITOVA Maria V. — St. Petersburg State University; 191024, Connya St. 4, St. Petersburg, Russia; e-mail: ddemidova@mail.ru
© Санкт-Петербургский государственный политехнический университет, 2013
